Leptonic Nonthermal Emission from Supernova Remnants Evolving in the Circumstellar Magnetic Field

Recently, van Marle et al. (2015) performed magneto-hydrodynamical simulations of the wind bubble around the massive star, but ignored the magnetic field in the wind material. Both RSGs and WR stars exhibit strong surface magnetic fields that should be transported with the stellar wind. Additionally, a turbulent dynamo might convert a sizable fraction of the bulk flow energy to magnetic field (Marcowith et al. 2018, and references therein).

In the free-wind zone, the magnetic field is mostly azimuthal due to the stellar rotation. Its amplitude is well described by (Voelk & Forman 1982)

$$\begin{eqnarray}&&{B}_{\mathrm{wind}}(r)={B}_{* }\,\displaystyle \frac{{R}_{* }}{{r}^{2}}\,\sqrt{\displaystyle \frac{{V}_{\mathrm{rot}}^{2}{r}^{2}+{R}_{* }^{2}{V}_{\mathrm{wind}}^{2}}{{V}_{\mathrm{rot}}^{2}+{V}_{\mathrm{wind}}^{2}}},\end{eqnarray} \tag{ 9 }$$

with B_* denoting magnetic field strength at the surface of the star, R_* and V_rot as the radius and the rotation velocity of the star, and V_wind the wind speed. We ignore the very inner region with predominantly radial magnetic field, because that is relevant only during the first days and weeks of the SNR. It is still important to have a rough estimate of the V_wind/ V_rot ratio to correctly scale the magnetic field strength.

The initial rotation velocities of massive stars are quite high, reaching a few hundreds of kilometers per second (Hunter et al. 2008; Maeder & Meynet 2012). At the RSG phase, the large size and high mass-loss rate makes the stars slow rotators with typical velocities, predicted by all stellar models, of a few km s⁻¹ (Maeder & Meynet 2012, and references therein). These estimates lead to

$$\begin{eqnarray}&&{B}_{\mathrm{wind}}(r)=(0.1-0.5){B}_{* }\displaystyle \frac{{R}_{* }}{r},\mathrm{for}\quad r\gtrsim {R}_{* }{V}_{\mathrm{wind}}/{V}_{\mathrm{rot}}\end{eqnarray} \tag{ 10 }$$

for the FW of a RSG.

The rotation velocity of WR stars is very uncertain (see e.g., Crowther 2007). For some stars, the corotation of large-scale structures indicates values ranging from a few tens to a few hundreds of kilometers per second (Chené & St-Louis 2008, 2010; Hénault-Brunet et al. 2011), as do strong deviation from spherical symmetry (Harries et al. 1998). On that basis we find

$$\begin{eqnarray}&&{B}_{\mathrm{wind}}(r)\sim (0.01-0.1){B}_{* }\displaystyle \frac{{R}_{* }}{r}.\end{eqnarray} \tag{ 11 }$$

RSGs exhibit relatively weak surface magnetic field of 1–10 G (Tessore et al. 2017), but a very large size of a few hundred R_⊙. On the contrary, WR stars are argued to be strongly magnetic, up to 1000 G, because the final product of their evolution are highly magnetized pulsars or magnetars (Crowther 2007). There is some evidence that magnetic field in some parts of the wind can reach 100 G implying a surface magnetic field of 1000 G (de la Chevrotière et al. 2014). WR stars are also quite compact, with a radius of only 1–10 R_⊙ (Crowther 2007).

At the transition to the SW zone, we assume a step-like change of the magnetic field. For the WR progenitor, the field strength increases by a factor of 4 at the termination shock and will then linearly increase with radius as the wind speed falls ∝ r⁻². In the RSG scenario, we have a transition to the main-sequence wind with highly uncertain magnetic-field strength. In line with the assumption of B² ∝ ρ we assume a decrease by a factor of 10. Beyond the transition point, we assume that the magnetic field stays constant as there is no flow.

Denoting as B₀ the initial magnetic field close to the star, the circumstellar magnetic field is expressed as

$$\begin{eqnarray}{B}_{\mathrm{CM}}(r)=\left\{\begin{array}{cl}{B}_{0}\displaystyle \frac{{R}_{* }}{r}, & r\leqslant {r}_{\mathrm{TP}},\\ 0.1{B}_{0}\displaystyle \frac{{R}_{* }}{{r}_{\mathrm{TP}}}, & r\gt {r}_{\mathrm{TP}}\end{array}\right.\end{eqnarray} \tag{ 12 }$$

for the RSG and

$$\begin{eqnarray}{B}_{\mathrm{CM}}(r)=\left\{\begin{array}{cl}{B}_{0}\displaystyle \frac{{R}_{* }}{r}, & r\leqslant {r}_{\mathrm{TP}},\\ 4{B}_{0}\displaystyle \frac{{R}_{* }}{{r}_{\mathrm{TP}}}\displaystyle \frac{r}{{r}_{\mathrm{TP}}}, & r\gt {r}_{\mathrm{TP}}\end{array}\right.\end{eqnarray} \tag{ 13 }$$

for the WR star. In the following, we consider three configurations of the adopted values for B₀ and R_* for both scenarios, RSG and WR (see Table 1). HNA and MNA models correspond, respectively, to a high and a moderate magnetic field strength within the range of allowed values discussed above. The AMP model adopts a field strength lying between those of the HNA and MNA models but assumes further amplification of the field upstream of the shock (see Section 3.2).

The standard gas-dynamical equations are solved, meaning we assume that the magnetic field is dynamically unimportant due to its low strength and that the remnant is not in the radiative phase yet (Petruk et al. 2016). The equations are solved in 1D for spherical symmetry.

We initialize the ejecta profile by a constant density, ρ_c, up to the radius r_c, followed by a power-law distribution up to the ejecta radius R_ej with index n = 9 for the core-collapse explosion,

$$\begin{eqnarray}\rho (r)=\left\{\begin{array}{cl}{\rho }_{{\rm{c}}}, & r\lt {r}_{{\rm{c}}},\\ {\rho }_{{\rm{c}}}{\left(\displaystyle \frac{r}{{r}_{{\rm{c}}}}\right)}^{-n}, & {r}_{{\rm{c}}}\leqslant r\leqslant {R}_{\mathrm{ej}},\end{array}\right.\end{eqnarray} \tag{ 14 }$$

beyond which the density profile is that of the circumstellar medium described in Section 2. The velocity of the ejecta is defined as

$$\begin{eqnarray}&&{v}_{\mathrm{ej}}(r)=\displaystyle \frac{r}{{T}_{\mathrm{SN}}},\end{eqnarray} \tag{ 15 }$$

where T_SN = 1 yr is the initial time set for hydrodynamic simulations. Then for the assumed ejecta mass, M_ej = 3 M_⊙, and energy, E_ej = 10⁵¹ erg (for both RSG and WR), and defining the radius of the ejecta as multiple of r_c, R_ej = xr_c, the initial conditions for simulations can be written as

$$\begin{eqnarray}&&{r}_{{\rm{c}}}={\left[\displaystyle \frac{10}{3}\displaystyle \frac{{E}_{\mathrm{ej}}}{{M}_{\mathrm{ej}}}\left(\displaystyle \frac{n-5}{n-3}\right)\left(\displaystyle \frac{1-\tfrac{3}{n}{x}^{3-n}}{1-\tfrac{5}{n}{x}^{5-n}}\right)\right]}^{1/2}{T}_{\mathrm{SN}};\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\rho }_{{\rm{c}}}=\displaystyle \frac{{M}_{\mathrm{ej}}}{4\pi {r}_{{\rm{c}}}^{3}}\displaystyle \frac{3(n-3)}{n}{\left(1-\displaystyle \frac{3}{n}{x}^{3-n}\right)}^{-1};\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{v}_{{\rm{c}}}=\displaystyle \frac{{r}_{{\rm{c}}}}{{T}_{\mathrm{SN}}}.\end{eqnarray} \tag{ 18 }$$

The spatial grid used in the hydrodynamic simulations extends to r = 30 pc with 262,144 linearly spaced grid points.

Upstream of the forward shock of the SNR, the ambient magnetic field is assumed to be amplified by a factor k and to exponentially decline on a length scale Δ l = 0.05R_sh, until it reaches the field strength in the far-upstream medium. Close to the transition point, r_TP, at which the ambient field abruptly changes, the thickness of the precursor region is redefined as Δ l = (r_TP–R_sh). After the transition point it is changed back to 5% of the shock radius.

Assuming a shock compression ratio of 4, the immediate downstream magnetic field is then given by

$$\begin{eqnarray}&&{B}_{d}=\sqrt{11}{{kB}}_{\mathrm{wind}}({R}_{\mathrm{sh}}).\end{eqnarray} \tag{ 19 }$$

The downstream magnetic field is passively transported with the plasma flow, following the induction equation for ideal MHD (Telezhinsky et al. 2013). For HNA and MNA models, we assume that the magnetic field is not amplified upstream of the shock (k = 1), while for the AMP model, the amplification factor is k = 5. The field amplification factor, k, has a particular impact on the rate of synchrotron cooling.

Describing the amplification of the magnetic field only by the parameter k is of course a crude simplification of the microphysics involved. As far as the resonant streaming instability is concerned, the amplification of turbulence should saturate a level of δ B ≈ B₀ (Zweibel 1979). An amplification to $\delta B\approx \sqrt{{v}_{\mathrm{sh}}{U}_{\mathrm{cr}}/c}\geqslant {B}_{0}$ in the precursor is possible, if the nonresonant streaming instability is involved (Lucek & Bell 2000; Bell & Lucek 2001). Cristofari et al. (2021) argue that the saturation level ought to be ∝ ${v}_{\mathrm{sh}}^{3/2}$ for the Bell scenario and by ∝ v_sh for the resonant streaming instability.

The particulars of the nonresonant amplification are clearly beyond the scope of this paper. The nonresonant instability saturates by a back-reaction of the thermal plasma to CR-streaming or a modification of the bulk flow (Riquelme & Spitkovsky 2009; Niemiec et al. 2010; Kobzar et al. 2017). In both cases, the saturation level is independent of the initial level of turbulence or the ambient magnetic field. However, the growth of turbulence is limited by the time that is available for the growth of turbulence, because the turbulence needs to be replenished on timescales of $D(p)/{v}_{\mathrm{sh}}^{2}$ (Brose et al. 2016). For Bell modes at the saturation amplitude $\sqrt{{v}_{\mathrm{sh}}{U}_{\mathrm{cr}}/c}$ throughout the cosmic-ray precursor one would not even get a single exponential growth cycle (Pohl 2021). This was recently shown to be the limiting factor in very young SNRs in dense environments, where the conditions for particle acceleration and turbulence growth were deemed ideal (Inoue et al. 2021). Additionally, enhanced damping might suppress further turbulence growth when the level of δ B/B₀ ≈ 1 is surpassed (Brose et al. 2021). In both cases the amplified field can be approximated by a scaling of the existing, ambient field.

Behind the shock, the amplified magnetic field may be quickly damped (Pohl et al. 2005), and so the field would fall back to the large-scale structure that we describe here.

In Figure 1 we show magnetic-field profiles for the AMP model at different SNR ages for the RSG (top panel) and WR (bottom panel) scenarios. Solid lines show the distribution of the magnetic field at a very early stage of the SNR evolution, dashed lines depict the moment shortly before the transition to the SW and dashed–dotted lines correspond to the evolution in the SW zone. Gray lines illustrate for reference the profile for a constant ambient magnetic field of 5 μG. The magnetic field is generally stronger for a RSG progenitor than for a WR, and we should expect stronger synchrotron cooling in the SNRs of RSG explosions. Even after the transition to the weaker field in the SW, about 1 μG, the downstream magnetic field remains at the level of several hundred microgauss, as is observed in Cas A (Cowsik & Sarkar 1980; Abeysekara et al. 2020).

Zoom In Zoom Out Reset image size

As long as the SNR evolves inside the FW, we see a similar radial variation in the magnetic-field profiles for the RSG and WR scenarios. The circumstellar magnetic field is compressed at the shock and further downstream gradually increases in strength to peak at the CD. Upstream of the reverse shock the magnetic field is very weak.

When the SNR shock moves into the SW, we see differences between the RSG and WR models, because in the case of a RSG the circumstellar magnetic field becomes weaker, whereas it turns stronger for a WR progenitor. Moreover, the posttransition phase of the RSG case model is characterized by an additional hump closely behind the shock (about 0.5 pc inward after 1000 yr, see top panel of Figure 1). This local boost of the magnetic field is caused by compression at a new CD between the accelerated forward shock and a reflected shock. This second CD also coincides with a pile-up of particles, and both effects combined lead to additional ring structure in the synchrotron morphology of the source, which we shall further discuss in Section 4.4.

We simulate the evolution of the differential particle density, N, by solving the transport equation in the form

$$\begin{eqnarray}&&\displaystyle \frac{\partial N}{\partial t}={\rm{\nabla }}(D{\rm{\nabla }}N-{\boldsymbol{v}}N)-\displaystyle \frac{\partial }{\partial p}\left((N\dot{p})-\displaystyle \frac{{\rm{\nabla }}{\boldsymbol{v}}}{3}{Np}\right)+Q,\end{eqnarray} \tag{ 20 }$$

where D is the spatial diffusion coefficient, v is the plasma velocity, $\dot{p}$ represents the energy losses (in our case synchrotron losses), and Q is the source term.

The diffusion is assumed to scale with the Bohm diffusion coefficient

$$\begin{eqnarray}&&D(p)={\eta }_{{\rm{B}}}\frac{{{pc}}^{2}}{3{eB}},\end{eqnarray} \tag{ 21 }$$

where η_B ≥ 1 is the Bohm factor and a measure of acceleration efficiency. A recent systematical study of the acceleration efficiency suggested a possible time evolution of the diffusion scale factor from η_B ≃ 10 at several hundred years to η_B ≈ 1 after several thousand years (Tsuji et al. 2021). The implied extremely high acceleration efficiency obtained for older SNRs is somewhat dubious and fully relies on the assumption that the maximum energy of electrons is limited by synchrotron cooling and the cutoff energy of the X-ray spectrum is independent of the magnetic-field strength (Zirakashvili & Aharonian 2007). Tsuji et al. (2021) motivate this assumption by examining the thin filament-like outer X-ray rims that are supposedly cooling limited. It was, however, shown already on several occasions that the thickness and brightness of X-ray filaments can be equally well explained by the damping of the magnetic field downstream of the shock, i.e., spatial distribution of the magnetic field, rather than its high strength (Pohl et al. 2005; Rettig & Pohl 2012; Sushch et al. 2018). The inaccuracy of the assumption for older SNRs is also supported by estimates of η_B < 1 for two SNRs, namely Vela Jr. and HESS J1731–347, which should be impossible. Therefore, in this study we adopt the value η_B = 10 for the entire evolution of the SNR up to an age of 2000 yr.

We use the thermal leakage injection model (Blasi et al. 2005),

$$\begin{eqnarray}&&Q=\chi {n}_{{\rm{u}}}({V}_{\mathrm{sh}}-{V}_{\mathrm{wind}})\delta (R-{R}_{\mathrm{sh}})\delta (p-{p}_{\mathrm{inj}}),\end{eqnarray} \tag{ 22 }$$

where χ is the injection efficiency parameter, n_u is the plasma number density in the upstream region, V_sh is the shock speed, V_wind is the wind speed upstream of the shock, R_sh is the shock radius, and p_inj = ξ p_th is the injection momentum, defined as a multiple of the thermal momentum in the downstream plasma. The injection efficiency for the compression ratio of 4 is determined as

$$\begin{eqnarray}&&\chi =\displaystyle \frac{4}{\sqrt{\pi }}\displaystyle \frac{{\xi }^{3}}{{e}^{{\xi }^{2}}}.\end{eqnarray} \tag{ 23 }$$

We solve the transport equation for electrons in spherical symmetry using the RATPaC code as described in Telezhinsky et al. (2012a, 2013), Sushch et al. (2018), Brose et al. (2019), taking into account only a forward shock and ignoring a reverse shock. Resulting electron spectra at different moments of time are then used to simulate synchrotron and IC (scattering on CMB photons) emission from the remnant. In this study we do not consider acceleration of protons and also assume that cosmic rays will not dynamically affect the SNR evolution.

First, we simulate the spectra of electrons and the nonthermal emission at the early stage of the SNR evolution, while it still expands in the FW. Figure 2 shows the evolution of the spectrum of accelerated electrons for the four models of the ambient magnetic field for RSG (top panel) and WR (bottom) progenitors. The differential density of electrons is scaled with ${R}_{\mathrm{sh}}^{2}{E}^{2}$ for illustration purposes. For reference, the gray dashed–dotted line shows the evolution of the electron spectrum for the constant ambient magnetic field of 5 μG.

Zoom In Zoom Out Reset image size

For the RSG progenitor, synchrotron cooling can be substantial even if the magnetic field is not amplified upstream of the shock. The blue line depicts the situation for a strong, but unamplified circumstellar magnetic field. The electron spectrum features a cooling break at a few tens of GeV after 100 yr, which eventually moves to higher energies as the field strength decreases. Already for the moderate magnetic field the cooling is inefficient and does not significantly modify the spectrum (green line), but the r⁻¹ distribution of the magnetic field in the free-wind zone can still leave imprints in the time evolution of the maximum electron energy. For a constant ambient magnetic field, the maximum energy is age limited; hence, it evolves as ${E}_{\max }^{\mathrm{age}}\propto {V}_{\mathrm{sh}}^{2}{tB}$, where V_sh is the shock velocity (see e.g., Reynolds 2008). At the free expansion stage the shock radius evolves as R_sh ∝ t^{(n − 3)/(n − m)} = t^6/7 (Truelove & McKee 1999), where the index of the ejecta density profile is n = 9 and the index of the density profile in the wind zone is m = 2. The shock speed declines as V_sh ∝ t^{−1/ 7}, yielding ${E}_{\max }^{\mathrm{age}}\,\propto \,{t}^{5/7}$ for the constant magnetic field, whereas for the realistic scaling, B ∝ r⁻¹, we expect ${E}_{\max }^{\mathrm{age}}\,\propto \,{t}^{-1/7}$.

For the strong circumstellar magnetic field encountered during the early stages of SNR evolution, the maximum energy is limited by synchrotron cooling. The balance of cooling and acceleration yields ${E}_{\max }^{\mathrm{syn}}\propto {V}_{\mathrm{sh}}{B}^{-1/2}$ (see e.g., Reynolds 2008), and hence the maximum energy only slowly grows with time, as ${E}_{\max }^{\mathrm{syn}}\propto {t}^{2/7}$.

For the WR progenitor star (Figure 2, bottom panel), synchrotron cooling is much less efficient even for a very strong magnetic field. Nonamplified field models result in featureless spectra roughly similar to those for a constant, weak ambient field. The amplified field model still shows a cooling break at a few hundred GeV early in the evolution, but already 300 yr after the supernova, the electrons have a featureless uncooled spectrum. The cutoff energy for the circumstellar magnetic field does not significantly increase with time, unlike that for the constant ambient field, indicating that it is still limited by synchrotron cooling.

Figure 3 shows the IC emission generated by CMB radiation scattering off the electrons. The gamma-ray spectra of SNRs resulting from the explosion of RSGs can feature a low-energy break and display spectral softening above the break. The spectral break appears below the GeV scale at early stages and gradually shifts toward higher energies, as the ambient magnetic field decreases. The spectral shape visually resembles the one expected for hadronic gamma-ray emission, introducing new difficulties in discrimination between the leptonic and hadronic scenarios, but it is much smoother then the hard transition in the pion-decay spectra, emphasizing the need for reliable data below 500 GeV.

Zoom In Zoom Out Reset image size

For the WR progenitor (Figure 3, bottom panel), we find only a weak spectral modification, and that only for the amplified field. By and large, there is no significant change of the spectral shape as compared to the constant low ambient magnetic field.

Similarly, the simulated synchrotron spectra (Figure 4) reveal a characteristic cooling break for the strong-field configurations in the RSG scenario, but nothing of that kind in the WR scenario. The flux normalization for circumstellar field models decreases with time, unlike for the constant ambient field, which can be directly attributed to the decrease of the magnetic field itself. This would naturally result in the decrease of the expected radio and X-ray emission, which is further discussed in Section 4.3.

Zoom In Zoom Out Reset image size

The transition to the SW zone is accompanied by an abrupt change of the physical properties of the medium. To investigate the impact of the change of magnetic field on the resulting particle and gamma-ray spectrum, we use AMP models described in Section 2.

The evolution of the electron, gamma-ray, and synchrotron spectra for an RSG progenitor is shown in Figure 5. Before the transition to the SW after nearly 800 yr, the evolution of the electron spectrum is governed by synchrotron cooling and the linear decrease of the magnetic field strength. As the remnant expands, the cooling break shifts to higher energies, while the maximum energy of electrons slowly increases (see previous section). After 800 yr, the dashed–dotted spectrum shows a significant shift of the cutoff to lower energies, which reflects the entrance of the cosmic-ray precursor into the weak-field region beyond the transition point. The diffusion coefficient in the precursor increases, and high-energy particles can effectively escape. As soon as the shock enters the dilute shocked medium, the shock velocity increases considerably, and particles are accelerated to high energies again, despite the weak magnetic field. Already after 900 yr, a second very-high-energy component emerges in the electron spectrum, which becomes more significant with time and eventually starts to dominate (Figure 5 top panel). In the IC spectrum, a hint of this component could be seen already after 800 yr (middle panel). Whereas some of the high-energy particles close to the shock can receive a further boost in energy, low-energy particles get advected downstream, escape the acceleration process, and are subjected to synchrotron cooling. This explains the softening of the spectrum above the cooling break energy at a few hundred GeV shortly after the transition. Later on, the spectrum in this energy range gets harder again due to the rise of the second component. Reaccelerated particles together with freshly accelerated particles constitute this very-high-energy component of the spectrum. Moreover, cooling becomes ineffective close to the shock due to the much lower magnetic field in the shocked medium. The resulting concave spectrum can be described as a superposition of spectra of two electron populations: old cooled electrons with a spectral break at around 100 GeV and new uncooled electrons with no cooling break and higher maximum energy. With time, further cooling of the first population shifts the break energy to lower energies, and the younger particles become more dominant due to their increasing number. This gradually changes the concavity of the spectrum until at some point in time the spectrum becomes featureless again.

Zoom In Zoom Out Reset image size

The concavity of the electron spectrum is reflected in the radiation spectra (Figure 5, middle panel for IC scattering and bottom panel for synchrotron), strongly changing their shape in the GeV band and the keV band. The appearance of the new high-energy component in the electron spectrum causes a strong temporal variation of the spectral shape in the GeV-to-TeV band. Right after the transition to the shocked main-sequence wind the very-high-energy spectrum is soft, but gradually hardens, ending up with a typical featureless IC spectrum when the new high-energy component takes over completely. This evolution takes around 700 yr for our setup, suggesting that the different phases of spectral evolution could potentially be detected. Although the main reason for the concavity of the spectrum is the abrupt change of the magnetic field, the change in density would regulate the timescale on which the gamma-ray spectrum settles to the typical form. The higher the amplitude of the density drop between the RSG and the main-sequence wind, the slower is the evolution of the spectrum, due to the lower injection rate of particles at the shock.

Similar variations are seen in the synchrotron spectrum, the main difference lying in the duration of the process. Since the synchrotron radiation directly depends on the magnetic-field strength, the new high-energy component takes longer to emerge in the synchrotron spectrum than in the IC spectrum. Additionally, the normalization of the spectrum decreases as the magnetic-field strength drops.

The situation with the WR progenitor is quite the opposite and more complicated, but eventually leads to a similar outcome, i.e., a concave spectrum (Figure 6). The magnetic field generated by the WR progenitor is weaker than that generated by the RSG progenitor, implying that synchrotron cooling is less efficient. Indeed, only at very early times a hint of a cooling break can be seen, and it disappears quickly. Nevertheless, the maximum energy of electrons is still cooling-limited, implying a slow rise of the maximum energy (see the previous section).

Zoom In Zoom Out Reset image size

As the forward shock approaches the transition point after nearly 400 yr, as in the RSG case the cosmic-ray precursor extends to the SW and the stronger magnetic field there, leading to a better confinement. The forward shock instantaneously decelerates by a factor of 1.5, when it enters the denser SW, allowing particles of the highest energy to return to the shock and re-enter the acceleration process. A prominent effect is the sharpened cutoff that is seen after 500 yr. As the remnant propagates through the SW with the linearly increasing magnetic field, synchrotron cooling becomes more efficient and imposes a cooling break and a spectral softening at high energies. Particles located sufficiently far downstream still see the lower magnetic field and hence weak cooling (see Figure 1), but they also have a lower maximum energy. This again creates a two-component particle spectrum and concave IC and synchrotron spectra. With time the strong magnetic field encountered in the SW advects to the far downstream and renders synchrotron cooling relevant also for the low energy component. Like for the RSG progenitor, eventually the high-energy component becomes dominant, allowing for considerable evolution of the spectral hardness in the GeV–TeV and eV–keV bands. This spectral modification is slower than for a RSG progenitor, opening an even larger window of opportunity for observing it. As expected, the maximum energy of accelerated electrons slowly decreases with time, as the shock slows down and the magnetic field increases linearly.

The interaction with the termination shock can briefly increase the compression ratio felt by particles and hence induce a spectral bump (Telezhinsky et al. 2013), but the duration of this event is too small to give a significant impact (Das et al. 2021), and therefore it is ignored here. More importantly, such a shock collision would generate a reflected shock propagating in a medium of lower density that might eventually be reflected off the CD of the SNR and collide with the forward shock. This process might repeat several times. The corresponding boost of the shock speed results in an instant increase of the maximum electron energy, modifying the spectrum accordingly. This would result e.g., in brightening and hardening of the X-ray and TeV emission on very short timescales. The recent detection of rapid variations in X-ray filaments of the Tycho's SNR (Okuno et al. 2020) may be attributed to the boosting of the forward shock by a trailing shock. A second possible explanation is dynamic magnetic turbulence (Bykov et al. 2008). After the boosting event, the shock decelerates again, and the maximum energy decreases until the next interaction happens. The reflected shocks bouncing back and forth between the forward shock and the CD would compress and reaccelerate particles, which can further modify the spectrum and the source morphology. For the WR progenitor the reflected shock is responsible for somewhat irregular variations of the electron and radiation spectra that are superposed on the main trend determined by the change of the magnetic field.

In Figure 7 we show the simulated flux in X-rays, high-energy gamma-rays, and very-high-energy gamma rays for our AMP models and RSG and WR progenitors. At the early stages of evolution (in the FW zone) both scenarios, RSG and WR, feature a power-law decrease of the X-ray flux, albeit with different slope. For the RSG progenitor, the dimming of the X-ray emission is determined exclusively by the 1/r scaling of the magnetic field strength. The decreasing upstream density, and likewise the injection rate, is fully compensated by the increasing shock surface, because the synchrotron cutoff does not fall into the energy range in question. For the WR progenitor the expected X-ray flux is sensitive to the change of the cutoff energy of the radiation spectrum, which depends on the maximum electron energy and magnetic field as ${E}_{\mathrm{cut}}^{\mathrm{ph}}\propto {E}_{\max }^{2}B$ (Reynolds 2008; Vink 2012). For the synchrotron-limited particle acceleration in the FW zone the cutoff energy evolves as ${E}_{\mathrm{cut}}^{\mathrm{ph}}\,\propto \,{t}^{-2/7}$. A consequence is faster dimming of the X-ray emission than in the RSG case. When the shock leaves the FW zone, the X-ray flux abruptly decreases in the case of the RSG and increases for the WR progenitor, on account of the jump of the magnetic field that, besides the synchrotron emissivity, also changes the confinement of high energy particles at the shock. In the WR scenario the flux keeps increasing for a few hundred years, which is partly due to reacceleration of old particles to higher energy and partly due to two episodes of shock-shock interaction. These episodes happen at around 560 yr and 690 yr and are clearly reflected in the VHE light curve (Figure 7 bottom panel). Later, the forward shock decelerates, causing a gradual decline of the X-ray flux enriched with sporadic brightenings caused by numerous shock–shock interactions.

Zoom In Zoom Out Reset image size

Unlike synchrotron radiation, the gamma-ray flux generated by IC scattering simply reflects the number of electrons and does not directly depend on the magnetic field strength. As long as the SNR expands into the FW, the HE and the VHE flux grow, because more and more particles get injected into the acceleration process. Like the X-ray flux, the gamma-ray flux abruptly drops at the transition to the main-sequence wind in the RSG scenario. This decline is mainly caused by synchrotron cooling of electrons accelerated in the FW zone that cannot be replaced on account of the abruptly decreased injection. In line with the energy dependence of the cooling time, the TeV-band flux recovers much faster than the GeV-scale flux. Later on the flux increases in both energy bands, reflecting the continuously growing number of accelerated particles.

For the WR progenitor we observe the opposite. At the transition to the SW, the flux in the entire gamma-ray band is boosted by the higher injection of particles and keeps growing until shock deceleration and synchrotron cooling introduce a slow decline, which, naturally, happens earlier in the VHE range than in the HE range. The VHE light curve also features several bumps, which reflect episodes of shock-shock interaction. The boost of the forward shock instantaneously increases the maximum energy of accelerated particles and hence the VHE flux. The counterparts of these bumps in the GeV light curve are much less prominent, if at all visible.

A declining flux of nonthermal X-rays was observed from Cas A (Patnaude et al. 2011) a remnant of the Type-IIb explosion of a RSG (Krause et al. 2008). The 4.2–6 keV flux from the whole remnant decreased at a rate of 1.5% yr⁻¹ over 11 yr, which the authors attributed to a declining maximum energy due to shock deceleration. They argued that the decrease of the magnetic field alone would not be sufficient to explain such a decline (Patnaude et al. 2011). The required deceleration rate of the forward shock considerably exceeds the observed deceleration with the expansion parameter of q≃0.65 for R_sh ∝ t^q (Patnaude & Fesen 2009). Note that this expansion parameter is compatible both with the Sedov–Taylor stage of an SNR in the FW, q = 2/(5 − m) = 2/3 with m = 2 for a ρ ∝ r^−m density profile, and with the free expansion of a core-collapse SNR in a medium of constant density, q = (n − 3)/(n − m) = 2/3 for n = 9 and m = 0. The latter reflects exactly the situation that we expect for the RSG scenario after the transition to the shocked-main-sequence wind. The expansion of the remnant is still ejecta-dominated as it expands in a very dilute medium, and the density of the SW is roughly constant. Interestingly, the average dimming rate of the X-ray emission near the transition time in our simulations, between 770 and 840 yr, is 1.4% yr⁻¹, compatible with that observed for Cas A. In our simulation, the dimming arises from the escape of the high-energy particles upstream of the shock, where the magnetic field abruptly decreases. For about 10 yr during the transition phase, the dimming can reach 10% yr⁻¹, but as soon as the forward shock enters the SW, the dimming normalizes to 1–1.5% yr⁻¹. At this time, the shock expands in a constant-density environment and thus decelerates again. Note that deviations from spherical symmetry can smear out this effect over time. It is plausible that Cas A at its age around 350 yr recently transitioned from the FW to the SW of its progenitor stellar wind bubble, and the decline of the X-ray emission that we observe reflects this transition.

The rate of change in the gamma-ray flux is smaller but can still reach a few % yr⁻¹ near the transition time. Such a flux change is potentially detectable with current instrumentation if we observe the source at the right moment. For WR progenitors, sporadic brightenings arising from shock-shock interactions are also strong enough to be detected. An identification as shock interaction events is likely impossible though, because their duration is longer than the time period available for observations. The general brightening in the TeV band in the RSG case, ∼0.25% yr⁻¹ beyond the transition point, and likewise the dimming for the WR scenario, ∼ −0.15% yr⁻¹, are too slow to be currently detectable.

The two spectral components of electrons, which we saw when the SNR shock left the FW, are not only spectrally distinct but also spatially separated. The pretransition population of electrons, which were accelerated when the shock propagated through the FW zone, are located further inside the remnant, whereas the freshly accelerated, more energetic electrons are found close to the shock. In addition, the speed of the forward shock is modified, changing the separation between the forward shock and the CD, which would likewise be reflected in the morphology. This should be particularly prominent in the radio emission, as low-energy electrons preferentially reside close to the CD, and the magnetic field piles up there (Rosenau & Frankenthal 1976).

In Figures 8 (RSG) and 9 (WR), we show simulated intensity maps for different ages and different photon energies. The top left panels depict the morphology when the SNR shock propagates through the FW zone. Both X-ray and TeV-scale emission are bound to the location of the shock where the high-energy particles reside. To be noted from the figures is that the radio shell is located farther at the CD where the low-energy particles pile up and the magnetic field is high. Time-variable X-ray filaments have been observed in the interior of Cas A (Uchiyama & Aharonian 2008) and in a number of other SNRs. They may be caused by the strong field at the CD, if the site of electron acceleration is very close (Lyutikov & Pohl 2004). Our simulation assumes that most of the acceleration happens at the forward shock, and so we do not see X-ray features associated with the CD.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the WR scenario the rim of 1 GeV emission is roughly co-located with the radio shell, reflecting the high density of low-energy electrons there. In the RSG case, however, the GeV shell is closer to the forward shock, on account of the strong synchrotron cooling of the TeV-scale electrons that produce the GeV-scale IC emission, whereas the radio-emitting GeV-band electrons remain unaffected.

In the RSG case, once the remnant transitions to the shocked medium, the forward shock speeds up. The separation between the shock and the CD increases, and so does that between the inner radio ring and the outer X-ray shell. The acceleration of the forward shock also results in the offset of the 1 GeV ring toward the interior of the remnant, because the production rate of TeV electrons has collapsed. Right after the transition, after 900 yr (upper right map of Figure 8), the highest X-ray and very-high-energy gamma-ray intensity is still seen close to the shock. At later times a second CD is formed, where particles are piled up and the magnetic field is compressed (Section 3.2). After 1500 yr, synchrotron cooling is much less effective, and the GeV–TeV morphology is similar to that seen during the early phase in the WR case, i.e., a shell location gradually shifting from the CD for lower energies to the forward shock for higher energies (bottom left map on Figure 8). The formation of the second CD also results in the appearance of a second radio ring that is prominent after 2000 yr (bottom right map on Figure 8). Initially, high X-ray intensity coincides with this second CD on account of the higher magnetic field, but later, after 2000 yr, a second X-ray shell emerges that is caused by freshly accelerated particles at the shock.

For the WR progenitor, the transition to the shocked medium immediately results in a double-ring structure in the radio and soft X-ray morphology (top left map of Figure 9). These two rings reflect the locations with enhanced magnetic field: the CD and the abruptly increased field at the shock. Later on, as the magnetic field and particles are advected downstream, the synchrotron morphology is mainly determined by the spatial distribution of electrons downstream. Note that, although the radio emission is generated farther in the interior of the remnant, it no longer coincides with the CD, because we do not find the highest electron density and strongest magnetic field there. At later stages of evolution (bottom maps on Figure 9), a second ring of GeV-band emission is formed by the electrons near the CD that gradually cool with time. After 1000 yr (bottom left map on Figure 9) the inner ring appears in the 100 GeV map, and after 2000 yr, it is most prominent at 1 GeV.

X-ray and IC very-high-energy gamma-ray emission are produced by electrons of roughly the same energy. It is therefore tempting to assume that both emissions should be produced in the same regions and with proportional emissivity, hence the gamma-ray morphology would strongly correlate with the X-ray morphology. Fukui et al. (2021) made this assumption to discriminate between the leptonic and hadronic scenario of the gamma-ray emission, likewise assuming that hadronic gamma-ray radiation should correlate with the gas distribution. The latter assumption is dubious, because we measure column densities and not the gas density itself, which determines the hadronic emissivity. Besides, the intensity maps in Figures 8 and 9 show clearly that the X-ray and the leptonic gamma-ray morphology do not necessarily correlate, as the former is also strongly dependent on the local magnetic field. To illustrate this in more detail we calculate projected radial profiles of the gamma-ray-to-X-ray flux ratio (Figures 10 and 11). Following Fukui et al. (2021), we calculate the gamma-ray flux in two energy bands: above 250 GeV (solid lines) and above 1 TeV (dashed lines). The X-ray flux is calculated in the energy range of 1–5 keV. For both types of SNR we do not see a constant flux ratio. Instead, a clear anticorrelation of the gamma-ray and X-ray morphology is visible for certain ages of the SNR at length scales of 5%–15% of the shock radius. For larger SNRs, this could be well resolved even by current generation instruments.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As already emphasized in the previous subsection, two spatially separated populations of electrons are seen after the transition of the forward shock to the shocked medium. We calculated spectral index maps at three different energies, 1, 10, and 100 GeV, to search for corresponding spatial variations of the spectral index, which may be relevant for the discrimination between leptonic and hadronic scenarios. We chose the same ages as for the emission maps in Figure 12 for the RSG case and Figure 13 for the WR progenitor. The index plotted is the photon spectral index, α, defined as

$$\begin{eqnarray}&&{N}_{\gamma }(E)=\displaystyle \frac{{dN}}{{dE}}\,\propto \,{E}^{-\alpha },\end{eqnarray} \tag{ 24 }$$

where N_γ (E) is the differential photon flux. At each location the index for particular energy E is calculated as

$$\begin{eqnarray}&&\alpha (E)=\mathrm{log}\displaystyle \frac{{N}_{\gamma }(E-{dE})}{{N}_{\gamma }(E+{dE})}/\mathrm{log}\displaystyle \frac{E+{dE}}{E-{dE}}.\end{eqnarray} \tag{ 25 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

During the pretransition stage of the RSG scenario, the GeV-band emission spectrum is flat in E² N_γ (E) representation, characteristic for a cooled electron spectrum (left hemisphere of the left map on Figure 12). The spectral index is roughly constant except at 1 GeV, where we see a signature of the cooling break. Indeed, electrons close to the shock have not had the time to cool at the energies relevant for the 1 GeV emission. The formation of two distinct electron populations can be very well seen after 900 yr, in the right hemisphere of the left map of Figure 12, particularly at 1 and 10 GeV. The "old" cooled particles are far behind the accelerated forward shock, and the "new" particles close to it exhibit a much harder spectrum. The difference in spectral index at 10 GeV reaches unity over a distance of 0.2R_sh, which may well be detectable with, e.g., Fermi-LAT, although spatial resolution can be an issue. This would not be a problem for CTA measuring at 100 GeV, but there the spectral index is again spatially constant, suggesting that we already observe "new" particles injected after the transition. Later on, in the right pie of Figure 12, one can see the "new" population of particles gradually taking over. The spatial variations disappear except for a slight (nondetectable) change in the ring related to the second CD.

For the WR case (Figure 13), the GeV spectrum is not immediately affected by the transition. Right before the transition, after 300 yr, and right after it, at 500 yr, the spectral index is roughly constant in the GeV band. A spatial variation of the spectral index emerges when the "old" low-energy component is sufficiently cooled for the cooling signature to appear at photon energies around 1 GeV. After 1000 yr (left hemisphere of the right map on Figure 13), the 10 GeV spectral map exhibits a very similar variation as in the immediate posttransition phase of the RSG model. The 1 GeV map, at the same time, still indicates a featureless uncooled spectrum. The strongest evidence of the formation of two components can be seen at 100 GeV, where the spectral index jumps by Δα ≃ 0.3 near 0.85R_sh. After 2000 yr the situation changes again: the particles have had time to cool, and the characteristic signature shifts to lower energies. The 1 GeV map at 2000 yr resembles that at 10 GeV and 1000 yr, and likewise for the 10 GeV map and the 100 GeV map. After 2000 yr, the 100 GeV map is already dominated by the "new" component.