© The Authors 2024

1 Introduction

Massive stars (M_* > 8 M_⊙) shape the morphology of their circumstellar medium (CSM) through ionizing radiation and stellar winds (Freyer et al. 2003; Mackey et al. 2014; Dwarkadas 2022). As massive stars evolve through different stages, mass loss in the form of stellar winds should vary (Maeder & Meynet 1987; Langer 2012; Meynet et al. 2015; Gormaz-Matamala et al. 2022). This variation is crucial for sculpting the environment of massive stars (Dwarkadas 2005, Dwarkadas 2007). The impact of stellar wind properties on the structure of wind bubbles has been explored by many studies. For instance, Garcia-Segura et al. (1996a,b) investigated the structure of wind bubbles around 35 M_⊙ and 60 M_⊙ stars in terms of the interacting stellar winds from different stages of evolution, and the bipolar morphology of nebulae around luminous blue variable (LBV) and Wolf-Rayet (WR) stars as a consequence of an asymmetric ambient medium from prior evolutionary stages has been studied in Dwarkadas & Owocki (2002) and Meyer (2021).

At the end of their lives, massive stars explode as core-collapse supernova remnants (SNRs). Hence, the generated SNR blast waves expand inside the complex wind-blown bubbles. Several SNRs, such as Cas A (van Veelen et al. 2009; Orlando et al. 2016), SN 1987A (Blondin & Lundqvist 1993; Bear & Soker 2018), and G292 + 0.8 (Park et al. 2002; Temim et al. 2022), appear to evolve inside the wind-blown bubble. Evidently, the dynamics of the SNR blast waves are regulated by the CSM parameters (Chevalier & Liang 1989; Ciotti & D Ercole 1989; Dwarkadas 2005, 2007; Meyer et al. 2021, 2022, 2023; Meyer & Meliani 2022; Velázquez et al. 2023), and they should not only differ from that of a uniform interstellar medium (ISM) but can also depend on the type of progenitor star. However, determining the progenitor type of a core-collapse SNR from its morphology is quite difficult, as the observed asymmetries in SNRs can also be developed by the magnetic field and explosion mechanisms, for instance, “ear”-like structures observed in radio and X-ray emission have been interpreted as being indicative of the interplay between the remnants and bipolar CSM created by LBV progenitors (Chiotellis et al. 2021; Ustamujic et al. 2021b), whereas other authors have attributed these structures to the orientation of the external magnetic field (Matt & Balick 2004; Townsend & Owocki 2005; Petruk et al. 2009; 2011) and explosion mechanisms, for example, jet-driven and single-lobed supernova (Khokhlov et al. 1999; Hungerford et al. 2005).

Supernova remnants are assumed to be major sources (Baade & Zwicky 1934; Drury et al. 1994; Büsching et al. 2005; Blasi 2013) of galactic cosmic rays (CRs), where CRs are considered to be accelerated by the diffusive shock acceleration (DSA) process (Fermi 1949; Bell 1978; Drury 1983). In the core-collapse scenario, particle acceleration along with emission (Dwarkadas 2013) should be influenced by the interactions between the SNR and the wind-blown bubble. Particle acceleration during the expansion of the SNR through the wind bubble has been studied assuming Bohm diffusion in Voelk & Biermann (1988); Berezinskii & Ptuskin (1989); Berezhko & Völk (2000). Recently, Sushch et al. (2022), Kobashi et al. (2022) used simplified flow profiles to investigate the influence of an ambient medium shaped by red super giant (RSG) and WR winds on the non-thermal emission from the remnants.

In our previous study (Das et al. 2022), the spectral evolution and emission morphology for SNR-CSM interaction were elaborately investigated using a realistic representation of CSM for the stellar track of a 60 M_⊙ star and for Bohm-like diffusion for energetic particles. In this study, the complete hydrodynamic evolution of the CSM throughout the lifetime of a star has been considered; the treatment of magnetic turbulence was not self-consistent. The impact of the self-generated magnetic field amplification on the maximum attainable particle energies and the resulting softer spectra at higher energies during later times were intricately explored in Brose et al. (2016, 2020) for type-Ia SNR.

Probing the particle acceleration and non-thermal emission in SNRs within wind bubbles by the combined influence of the CSM and the CR-streaming instabilities and their impact on CR diffusion would be desirable. In this paper, we explore the spectral evolution of accelerated particles at forward shocks of SNRs with 20 M_⊙, and 60 M_⊙ progenitors during their expansion inside the wind-blown bubbles while considering the time-dependent evolution of self-generated magnetic turbulence. The 20 M_⊙ star evolves through an RSG phase as the post-main sequence (MS) stage, whereas the 60 M_⊙ star has LBV and WR phases after the MS stage. Our study demonstrates the following aspects: the difference in spectral shapes arising from the morphological dissimilarity of the CSM of both SNRs; in both scenarios, the softening of particle spectra at higher energies during the later stages of evolution, on account of the weak driving of magnetic turbulence and the escape of high-energy particles from the remnants; the temporal evolution of the spectra for different emission channels as well as emission from the escaped particles around the remnants; and the evolution of the morphology for the different energy bands during the lifetime of remnants and its dependence on the structure of the CSM.

2 Numerical methods

In this section, we introduce the numerical methods applied in the presented study. We have modelled the diffusive shockaccel-eration (DSA) at SNR forward shock in test-particle approximation by combining the hydrodynamic evolution of SNR inside CSM, large-scale magnetic field evolution in addition to the time-dependent treatment of the magnetic turbulence, and finally the solution for CR transport equation. We have numerically solved the particle acceleration and hydrodynamics, respectively, with RATPaC (Radiation Acceleration Transport Parallel Code) (Telezhinsky et al. 2012, 2013; Brose et al. 2020; Sushch et al. 2018) and the PLUTO code (Mignone et al. 2007; Vaidya et al. 2018).

2.1 Hydrodynamics

The Euler hydrodynamic equations including an energy source/sink term, considering the magnetic field as dynamically unimportant, can be described as $$\frac{\partial }{\partial t}\left(\begin{array}{c}\rho \\ m\\ E\end{array}\right)+\nabla {\left(\begin{array}{c}\rho u\\ mu+PI\\ \left(E+P\right)u\end{array}\right)}^T=\left(\begin{array}{c}0\\ 0\\ S\end{array}\right)$$(1) $$\frac{\rho u^2}{2}+\frac{P}{\gamma -1}=E;\gamma =\frac{5}{3},$$(2)

where ρ, u, m, P, E, and S are the mass density, flow velocity, momentum density, thermal pressure, total energy density, and source-sink term to include the optically thin cooling and radiative heating for the construction of CSM at the pre-supernova stage, respectively. I is the unit tensor.

2.1.1 Construction of the CSM at the pre-supernova stage

We start with constructing the wind-blown bubble around the progenitor star with 20 M_⊙ at solar metallicity (Z = 0.014) from the zero-age main sequence (ZAMS) to the pre-supernova stage. Das et al. (2022) already simulated the wind bubble for 60 M_⊙ star, and we use the same methodology also for the 20 M_⊙ progenitor. Hence, we will only provide a synopsis here.

We have performed hydrodynamic simulations with PLUTO code in 1-D spherical symmetry with 50000 grid points of uniform spacing out to R_max = 150 pc. To initialise the simulation, the radially symmetric spherical supersonic stellar wind has been injected into a small spherical region of radius 0.06 pc at the origin, using the stellar evolutionary tracks for 20 M_⊙ star (Ekström et al. 2012), and 60 M_⊙ star (Groh et al. 2014). The ISM is assumed to have a constant number density, n_ISM = 1 cm⁻³. The stellar wind density, ρ_w, can be expressed as $${\rho }_{\text{wind\hspace{0.17em}}}=\frac{\dot{M}(t)}{4\pi r^2{V}_{\text{w}}(t)},$$(3)

where r is the radial coordinate, and Ṁ and V_w represent the time-dependent mass-loss rate and the wind velocity, respectively. The wind parameters for the 20 M_⊙ star have been derived from the Geneva library of stellar models (Ekström et al. 2012) and for the 60 M⊙ star, wind properties have been taken from Groh et al. (2014).

Figure 1 depicts the stellar mass in panel a, mass-loss rate in panel b, wind velocity in panel c during stellar evolution at different stages that are used as initial conditions in the hydrodynamic simulation for CSM. The CSM evolution has been simulated from ZAMS to the pre-supernova phase, over 8.64 million years for the 20 M_⊙ star and over 3.95 million years for the 60 M_⊙ star. The 20 M_⊙ star evolves through the MS phase for approximately 6.3 million years followed by the RSG stage and the 60 M_⊙ star evolves through the MS stage for approximately 3 million years followed by LBV phase with the span of approximately 0.2 million years and WR phase. For modelling the evolution of the wind bubble, the second-order Runge-Kutta method has been applied to integrate Eqs. (1) and (2), using the Harten-Lax-Van Leer approximate Riemann Solver (hll) and finite volume methodology. Additionally, optically thin cooling and radiative heating have been included through the source/sink term, S = Φ(Τ, ρ), following the cooling and heating laws, $$\text{Φ}(T,\rho )=n^2(\text{Γ}(T)-\text{Λ}(T)),$$(4)

where n is the particle number density, Γ is the radiative heating term that represents the effect of stellar photons. In our adopted cooling curve, the cooling term Λ is composed of the contributions from hydrogen (H), helium (He) and metals (Z), whose contributions come from Wiersma et al. (2009) for a plasma at solar helium abundance (Asplund et al. 2009), and it governs the gas cooling at high temperatures. Additionally, for temperatures Τ < 6 × 10⁴ K, a cooling term for H recombination is added following the case B energy-loss coefficient of Hummer (1994). Collisionally excited forbidden lines (mainly from O and C elements) are included as described in Eq. (A9) of Henney et al. (2009), using a relative abundance O/H = 4.89 × 10⁻⁴ in number density (Asplund et al. 2009). The heating rate, Γ, mimics the ionisiation of recombining hydrogenic ions by photons from the photosphere, in which the free electrons receive the energy carried by an ionizing photon once the reionisiation potential of an H atom has been subtracted. In a statistical equilibrium, the rate follows from the recombination coefficient whose values are from the Table 4.4 of Osterbrock (1989). The heating term, Γ, represents the photoelectric heating of dust grains by the Galactic far-UV background. For Τ ≤ 1000 K, we used Eq. (C5) of Wolfire et al. (2003), and the electron-density profile, n_e, follows Eq. (C3) of Wolfire et al. (2003). For Τ > 1000 K we take the value of n_e interpolated from the CIE curve by Wiersma et al. (2009). This cooling strategy has been used in a series of papers (Meyer et al. 2014, 2020, 2023).

Figure 2 illustrates the structure of the pre-supernova CSM at 8.64 million years and 3.95 million years after ZAMS stage for 20 M_⊙ and 60 M_⊙ progenitors, respectively. The wind bubble generated by the 20 M_⊙ progenitor consists of dense RSG material that stretches until around r ≈ 15 pc, and the deceleration of slow RSG wind against the back-flowing MS material causes the formation of the RSG shell. This characteristic structure of CSM arising from RSG-MS interaction has been noted before for a 35 M_⊙ star (Dwarkadas 2007, Sect. 3), and likewise the absence of a wind termination shock in the RSG wind.

The temperature of the CSM around the 60 M_⊙ progenitor is very hot in the entire shocked wind material, extending out to r ≈ 60 pc, whereas the CSM around the lower-mass star shows high temperature only in the region of shocked MS wind. The wind bubble structure around the 60 M_⊙ star has been illustrated before in Das et al. (2022) and is reproduced in Fig. 2 to ease the comparison.

Comparative discussion about wind bubble structure. The structure of wind bubbles shaped during the lifetime of 20 M_⊙ and 60 M_⊙ progenitors reflects the history of the mass-loss rate and stellar-wind velocity as shown in Fig. 1. In literature, a myriad of studies exist concerning the wind bubble structure and the included physics in the numerical modelling. For instance, the properties of the ISM and the prescription of radiative heating and cooling play a significant role in shaping the wind bubbles, along with the stellar-wind parameters as shown for two different stars in Fig. 2. As the cooling timescale of gas depends on ρ⁻¹, for a high ISM density, the cooling timescale is small, which leads to rapid cooling of the gas. For instance, Toalá & Arthur (2011) considered an ISM density that is approximately 100 times higher than that in this paper, and therefore very rapid cooling generates a thin and dense shell, whereas for the bubbles illustrated in Fig. 2 a thin dense shell of shocked ISM does not form. Including only radiative cooling, Dwarkadas (2007) obtained a thin shell of shocked ambient material in the simulation of wind bubble around 35 M_⊙ star. The absence of radiative heating of the gas permits the gas to cool down to a very low temperature, and hence the shell becomes confined by the thermal pressure of the ambient medium. The prescription of radiative heating-cooling as in Eq. (4), with a balance temperature of 10⁴ K in the ISM, resists the cooling of shocked ambient gas in the shell and leads a thick shell of shocked material without any sharp piled-up of shocked ISM material. The shape of the shell made by shocked ISM should also depend on the temperature of the ambient medium as the thermal pressure therein opposes the expansion of the shell. We mention in passing that dense shells at radiative shocks can also be suppressed by magnetic pressure (Petruk et al. 2016).

The structure of this shell also relies on the wind parameters, as with a higher stellar mass, the velocity of stellar wind increases and so does its ram pressure, which supports the formation of prominent shells. The lifetime of stars decreases with stellar mass. Hence, the propagation of the shock front of the shell for more massive stars becomes limited by the available time, whereas for less massive stars, the shock front evolves until the pressure of the shell comes into equilibrium with the ambient thermal pressure, and the shell structure depends on the balance temperature in the simulation.

Therefore, the structure of the shell is quite parameter-dependent, particularly concerning the ambient medium and the prescription for heating-cooling, so care must be exercised when comparing different studies.

In this present study the structure of the CSM beyond the contact discontinuity of the wind bubble has negligible impact on the particle acceleration by the corresponding SNR, as particle acceleration becomes inefficient when the remnant reaches the CSM, and hence we stop our simulation within 1 pc beyond the contact discontinuity.

Fig. 1 Evolution of stellar mass (panel a), mass-loss rate (panel b), and wind velocity (panel c) as a function of stellar age (in Myr) during different evolutionary stages of the 20 M⊙ and 60 M⊙ stars. The evolutionary stages of the 20 M⊙ star are shown with red dashed lines and those of the 60 M⊙ star are shown with blue solid lines. We used these parameters as the initial condition in the hydrodynamic simulation for constructing the corresponding CSM at the pre-supernova stage at 8.64 million years for the 20 M⊙ star and 3.95 million years for the 60 M⊙ star. The stellar parameters for the 20 M⊙ and 60 M⊙ stars are taken from Ekström et al. (2012); Groh et al. (2014), respectively.

Fig. 2 Pre-supernova CSM profiles of the number density (n; panel a), the flow speed (u; in panel b), the thermal pressure (P; panel c), and the temperature (T; panel d) at 8.64 million years after the ZAMS stage for the 20M⊙ progenitor (on the left) and at 3.95 million years for the 60M⊙ progenitor after the ZAMS stage (on the right). The right panels are reproduced from Das et al. (2022). Vertical grey lines mark the boundaries of specific regions. For the 20M⊙ progenitor, these are the free RSG stellar wind (region 1), piled-up RSG wind (region 2), shocked MS wind (region 3), shocked ISM (region 4), and ISM (region 5). The RSG shell indicates the accumulation of decelerated, freely expanding RSG wind; RRSG denotes the transition between RSG and MS wind; and CD represents the contact discontinuity between shocked wind and the ISM. For the 60M⊙ progenitor, we distinguish the free stellar wind (region 1), the shocked LBV and WR wind (region 2), the shocked MS wind (region 3), the ISM (region 4), and the ambient ISM (region 5). The term RWT is the radius of wind termination shock, LBV shell denotes the dense shell created by the interaction between LBV wind and WR wind, and CD represents the contact discontinuity between the shocked wind and shocked ISM.

2.1.2 Supernova ejecta profile

The density of supernova ejecta, ρ_ej(r), is modelled as (Truelove & McKee 1999) $${\rho }_{\text{ej}}\left(r\right)=\{\begin{array}{ll}{\rho }_{\text{c}},\hfill & r\le r_c\hfill \\ {\rho }_{\text{c}}{\left(\frac{r}{r_c}\right)}^{-n}\hfill & r_{\text{c}}<r\le R_{\text{ej}},\hfill \end{array}$$(5)

where n = 9, as conventionally used for core-collapse explosions. The velocity profile for ejecta u_ej(r) follows a homologous expansion, $$u_{\text{ej}}=\frac{r}{T_{\text{SN}}},$$(6)

at the starting time of the hydrodynamic simulation, T_SN = 3 yr. The initial ejecta temperature is set to 10⁴ K.

The expressions for r_c and ρ_c as a function of the ejecta mass, M_ej, and explosion energy, E_ej will be $$r_{\text{c}}={\left(\frac{10E_{\text{ej}}}{3M_{\text{ej}}}\frac{n-5}{n-3}\frac{n-3x^{3-n}}{n-5x^{5-n}}\right)}^{1/2}{T}_{\text{SN}},$$(7) $${\rho }_{\text{c}}=\frac{M_{\text{ej}}}{4\pi r_{\text{c}}^3}3(n-3){\left(n-3x^{3-n}\right)}^{-1},$$(8)

where R_ej = xr_c and x = 2.5¹. Here, x is a free parameter that determines the ratio between the flat and steep parts of the ejecta density profile and its choice is motivated by the comparability of the ejecta and ambient densities at R_ej and numerical stability of simulations. The choice of x has a negligible effect on the dynamics of the shock; this was thoroughly tested through comparison of numeric simulations with analytic solutions (Chevalier 1982; Truelove & McKee 1999). Further, in comparison to Chevalier (1982) where the ejecta distribution extends to infinity, the truncation at R_ej provides a limit to the ejecta speed and feedback from pushing away the CSM that originally resided very close to the star. The explosion energy is fully transferred to ejecta, E_ej = 10⁵¹ erg. The ejecta mass, M_ej, for the 20 M_⊙ star is 3.25 M_⊙, and M_ej = 11.75 M_⊙ for the 60 M_⊙ star, assuming a left-over compact object on the neutron-star mass scale following $$M_{\text{ej}}=M_{\star }-{\displaystyle {\int }_{t_{{\text{t}}_{\text{ZAMS}}}}^{t_{\text{preSN}}}\dot{M}}(t)\text{d}t-M_{\text{CompactObject\hspace{0.17em}}}\left(1.4M_{\odot }\right).$$(9)

To launch the supernova explosions, the supernova ejecta profiles have been inserted in the pre-calculated pre-supernova CSM profiles, covering the inner 56 to 108 grid points, depending on the progenitor type and also on the ejecta mass. Then Eqs. (1) and (2) have been solved in the absence of heating or cooling (S = 0) using a Harten-Lax-Van Leer approximate Riemann Solver that employs the middle contact discontinuity (hllc), finite-volume methodology, and a second-order Runge-Kutta method. The numerical simulations of the supernova remnant with the PLUTO code have been carried out in 1-D spherical symmetry with a spatial resolution of about 0.0004 pc, about a factor 7.5 finer than that used for the wind bubble prior to the supernova explosion, using linear interpolation of the pre-supernova CSM grid. The proper solution of the cosmic-ray transport equation near the shock requires an exquisite spatial resolution, otherwise the cosmic-ray precursor is not resolved and the velocity jump at the shock is poorly reconstructed. The higher resolution of the grid describing the gas flow also helps in defining a sharp shock in the inhomogeneous spatial grid that we use for the cosmic-ray transport equation (cf. Sect. 2.3). To ensure a sharp transition in hydrodynamic parameters - density, velocity, pressure and temperature from upstream to downstream values at the shock, the hydro data from PLUTO code needs resharpening. Resharpening is done by blinding hydro data for density, velocity, pressure and temperature at 4 bins upstream and downstream of the shock and replacing these values with the linear extrapolation using further downstream and upstream values (Brose 2020, Chapter 4). These resharpened profiles are used for further necessary calculations in the simulation.

2.2 Magnetic field

The total magnetic field strength in the entire system is given by $$B_{\text{tot\hspace{0.17em}}}=\sqrt{B_0^2+B_{\text{turb\hspace{0.17em}}}^2},$$(10)

where B₀ and B_turb are the large-scale and turbulent magnetic field, respectively.

2.2.1 Large-scale field profile

Simulating the CSM magnetic field for the entire lifetime of progenitors by solving magneto-hydrodynamic (MHD) equations is out of the scope of this paper. Therefore, we have parametrised the magnetic field in the CSM around 20 M_⊙ star, similarly as we constructed the CSM magnetic field for the 60 M_⊙ progenitor (Das et al. 2022). The magnetic field in the stellar wind of a rotating star will be toroidal except for very close to the stellar surface. In the equatorial plane of rotation, it can be expressed as (Ignace et al. 1998; García-Segura et al. 1999; Chevalier & Luo 1994) $$B_{\varphi }=B_{\star }\frac{u_{\text{rot}}{R}_{\star }}{u_{\text{wind\hspace{0.17em}}}r}=\frac{B_{\star ,0}{R}_{\star }}{r}r>>R_{\star },$$(11)

where Β_* and R_* are the stellar surface magnetic field and radius, respectively, u_rot and M_wind represent the surface rotational velocity in the equatorial plane and the radial wind speed, respectively. A 20 M_⊙ star evolves through MS and RSG phases. At the RSG phase, the star reaches a very large size up to a few hundred times R_⊙, but becomes a slower rotator (Maeder & Meynet 2012), and the wind speed is 20–50 km s⁻¹. Measurements by Aurière et al. (2010); Tessore et al. (2017) suggest that a RSG has a weak surface field with 1–10G, and therefore we chose B_*.0(R_*/R_⊙) = 750 G. Equation (11) then defines the field strength for both free RSG wind and piled-up RSG wind, region 1 and 2 as marked in Fig. 2. At R_RSG = 16.3 pc ~ 7 × 10⁸ R_⊙, there is the transition to MS wind, which has a very uncertain magnetic field. We assume a decrease in magnetic field strength by a factor of 3 following $B\propto \sqrt{\rho }$ and the density jump at R_RSG. Throughout the MS wind the magnetic field is supposed to have a constant strength. In the ISM, the field strength has been chosen as 4 μG to provide a super-Alfvenic flow to the shell, marked as region 4 in Fig. 2. The field in the shell has been calculated assuming flux conservation (van Marle et al. 2015). In total, the magnetic field strength ($B_{0,20M_{\odot }}$) in the different regions mentioned in Fig. 2 of the wind bubble created by 20 M0 star can be expressed as $$B_{0,20M_{\odot }}=\{\begin{array}{ll}(1.07\text{\hspace{0.17em}μG})\frac{R_{\text{RSG}}}{r}\hfill & \text{\hspace{0.17em}regions\hspace{0.17em}}1\text{\hspace{0.17em}and\hspace{0.17em}}2\hfill \\ 0.35\text{\hspace{0.17em}μG}\hfill & \text{\hspace{0.17em}region\hspace{0.17em}}3\hfill \\ 4.68\text{\hspace{0.17em}μG}\hfill & \text{\hspace{0.17em}region\hspace{0.17em}}4\hfill \\ 4.0\text{\hspace{0.17em}μG}\hfill & \text{\hspace{0.17em}region\hspace{0.17em}}5.\hfill \end{array}$$(12)

For the 60 M_⊙ star, we used the same methodology as in Das et al. (2022) but changed the field strength slightly in order to maintain the magnetic field as dynamically unimportant. Hence, the magnetic field in the different regions of the CSM formed by a 60 M_⊙ progenitor can be written as $$B_{0,60M_{\odot }}=\{\begin{array}{ll}(0.63\text{\hspace{0.17em}\hspace{0.17em}μG})\frac{R_{\text{WT}}}{r}\hfill & \text{\hspace{0.17em}region\hspace{0.17em}}1\hfill \\ 2.52\text{\hspace{0.17em}μG}\hfill & \text{\hspace{0.17em}regions\hspace{0.17em}}2\text{\hspace{0.17em}and\hspace{0.17em}}3\hfill \\ 14.8\text{\hspace{0.17em}\hspace{0.17em}μG}\hfill & \text{\hspace{0.17em}region\hspace{0.17em}}4\hfill \\ 4.3\text{\hspace{0.17em}\hspace{0.17em}μG}\hfill & \text{\hspace{0.17em}region\hspace{0.17em}}5.\hfill \end{array}$$(13)

We determined the magnetic field initially in the supernova ejecta following (Telezhinsky et al. 2013, Sect. 3). Finally, the time-dependent computation of the frozen-in large-scale magnetic field transported passively with the hydrodynamic flow is given by the induction equation in 1D spherical symmetry (Telezhinsky et al. 2013): $$\frac{\partial B_0}{\partial t}=\nabla \times \left(u\times B_0\right).$$(14)

This method mimics MHD for negligible magnetic pressure.

2.2.2 Magnetic turbulence

The temporal and spatial evolution of the magnetic turbulence spectrum can be described by a continuity equation for magnetic spectral energy density per logarithmic bandwidth, E_w(r, k, t) (Brose et al. 2016): $$\frac{\partial E_w}{\partial t}=-\nabla \cdot \left(u{E}_w\right)-k\frac{\partial }{\partial k}\left(k^2{D}_k\frac{\partial }{\partial k}\frac{E_w}{k^3}\right)+2\left({\Gamma }_g-{\Gamma }_d\right)E_w,$$(15)

where k denotes the wavenumber, D_k represents the diffusion coefficient in wavenumber space describing cascading, and Γ_g, Γ_d are the growth and damping rates, respectively. This transport equation for the magnetic turbulence spectrum has been solved in 1D spherical symmetry, considering Alfvén waves as scattering centres for CRs.

Amato & Blasi (2009) suggested that although non-resonant CR streaming instabilities are likely the dominant way of magnetic field amplification during the free expansion phases and early Sedov-Taylor phase of SNR, when the SNR shock is relatively fast, at later times resonant modes provide efficient amplification. It is to be noted that non-resonant amplification is a very complex issue, and some of the commonly quoted concepts may be misperceptions. The magnetic field in the free wind, that the shock passes through early in the evolution, is mostly perpendicular to the shock normal, whereas the usual calculations of the growth rate for resonant (Skilling 1975; Bell 1978) and non-resonant (Bell 2004) modes is performed for streaming parallel to the magnetic field. Here we may have oblique streaming, but that implies that even fewer growth times are available before the plasma passes through the shock. Already for parallel shocks, one expects only very few growth times, that is, a few exponential growth cycles (Niemiec et al. 2008). This is one of the reasons why the peak amplitude of the non-resonant mode is generally less than naively estimated (Pohl 2021). The entire process, and hence its spatial profile, is highly nonlinear. It would be important to understand the cascading properties of the mode and the dependence on the wave spectrum of the CR diffusion coefficient, but that is beyond the scope of this paper. In the literature, cascading is often neglected and Bohm diffusion is simply assumed (e.g. Zirakashvili & Ptuskin 2008; Cristofari et al. 2021). Here, we enhanced by a factor A = 10 the growth term for the resonant streaming instability: $${\text{Γ}}_g=A\frac{v_A{p}^{2}v}{3E_w}\left|\frac{\partial N}{\partial r}\right|,$$(16)

where ν and p are the particle velocity and momentum respectively, ν_A is the Alfvén speed, N is the differential number density of CRs. The value A = 10 agrees with the observations of historical SNRs (Brose et al. 2021) and estimates of the growth rates operating at the early stages of CR-acceleration (Marcowith et al. 2018).

The spectral energy transfer process from a small wavenum-ber scale to a large wavenumber scale through turbulence cascading can be empirically described as a diffusion process in wavenumber space with coefficient (Zhou & Matthaeus 1990; Schlickeiser 2002): $$D_k=k^3{v}_A\sqrt{\frac{E_w}{2B_0^2}}.$$(17)

This form of the diffusion coefficient conforms with Kolmogorov scaling in the case of pure cascading from large scales, Ε_w ∝ k^−2/3. Furthermore, since ν_A ∝ B_tot and $B_{\text{turb\hspace{0.17em}}}^2=4\pi {\displaystyle \int \text{d}}\ln k{E}_w(k)$, D_k is more sensitive to Ε_w, when the turbulent field is amplified beyond the amplitude of the large-scale field, $$D_k=\{\begin{array}{lll}\sqrt{E_w}\hfill & \text{\hspace{0.17em}for\hspace{0.17em}}\hfill & E_w\ll B_0^2/8\pi \hfill \\ E_w\hfill & \text{\hspace{0.17em}for\hspace{0.17em}}\hfill & E_w\gg B_0^2/8\pi .\hfill \end{array}$$(18)

2.2.3 Diffusion coefficient

The diffusion coefficient governs the rate of particle acceleration, the maximum attainable energy of the particles (Lagage & Cesarsky 1983; Schure et al. 2010), and the spatial distribution of the accelerated particles in the remnant. Since we have considered Alfvén waves as the scattering centres for CRs, the resonance condition is $$k_{\text{res}}=\frac{\left|q\right|B_0}{pc},$$(19)

where k_res represents the resonant wavenumber and q is the particle charge. Then, the diffusion coefficient for CRs coupled to E_w can be described as (Bell 1978; Blandford & Eichler 1987), $$D_{\text{r}}=\frac{4v}{3\pi }{r}_g\frac{U_B}{E_w},$$(20)

where U_B is the energy density of the large-scale magnetic field and r_g represents the gyro-radius of particles in the total magnetic field (B_tot). The growth of Alfvén waves, damping, spectral energy transfer through cascading, and the spatial transport of waves suggest that the magnetic turbulence spectra should exhibit a complicated shape rather than a featureless and flat spectrum described by the Bohm-like diffusion coefficient. In addition, the time-dependent calculation of E_w, following Eq. (20), provides the self-consistent diffusion coefficient, hence resulting in a more realistic picture of CR acceleration in SNR, in comparison to the simplified Bohm-like diffusion elucidated in our previous study with core-collapse SNR (Das et al. 2022).

As an initial condition, we assumed a magnetic turbulence spectrum that yields a diffusion coefficient a factor of 100 lower than that for Galactic propagation of cosmic rays (Trotta et al. 2011): $$D_0=\left({10}^{27}\hspace{0.17em}{\text{cm}}^2\hspace{0.17em}{\text{s}}^{-1}\right){\left(\frac{pc}{10\text{\hspace{0.17em}GeV}}\right)}^{\frac{1}{3}}{\left(\frac{B}{3\text{\hspace{0.17em}μG}}\right)}^{-\frac{1}{3}}.$$(21)

2.3 Particle acceleration

The time-dependent transport equation for the differential number density of CRs, N(p), can be written as $$\frac{\partial N}{\partial t}=\nabla \left(D_{\text{r}}\nabla N-uN\right)-\frac{\partial }{\partial p}\left(\dot{p}N-\frac{\nabla \cdot u}{3}Np\right)+Q,$$(22)

where $\dot{p}$ corresponds to the energy loss rate for synchrotron and inverse Compton losses of electrons (Brose et al. 2021), and Q denotes the source term (Skilling 1975).

This equation has been solved in test-particle approximation for spherical symmetry in RATPaC, applying implicit finite-difference algorithms implemented in the FiPy package (Guyer et al. 2009). In this paper, the non-linear shock modification by CR pressure is negligible and ignored, because the cosmic-ray pressure has always constrained below 10% of the shock ram pressure (Kang & Ryu 2010). We solved Eq. (22) on an inhomogeneous spatial grid, r′ = r/R_sh(t), and this co-ordinate transformation provided excellent spatial resolution near the shock (Brose et al. 2020; Das et al. 2022): $$r^{\prime }-1={\left(x_{*}-1\right)}^3,$$(23)

where R_sh is the shock radius.

Fig. 3 Temporal evolution of the forward shock parameters. The parameters include radius, Rsh; speed, Vsh; and the sub-shock compression ratio, Cr. In the upper panel, we provide the angular scale, θsh, for an SNR located 1000 pc away. We also mark the interactions of the forward shocks with discontinuities in the wind bubbles of the two types of progenitor stars.

2.4 Injection of particles

The source term, Q, in the transport equation, can be expressed as $$Q={\eta }_{\text{inj}}{n}_{\text{u}}\left(V_{\text{sh}}-u_{\text{u}}\right)\delta \left(R-R_{\text{sh}}\delta \right)p-p_{\text{inj}},$$(24)

where η_inj is the injection efficiency, n_u and u_u are the upstream plasma number density and velocity in observer’s frame, respectively, V_sh is the shock velocity in the observer’s frame, and $p_{\text{inj\hspace{0.17em}}}=\xi p_{\text{th}}=\xi \sqrt{2m{k}_B{T}_d}$, represents the momentum of injected particles. The injection efficiency is defined following the thermal leakage model (Blasi et al. 2005), $${\eta }_{\text{inj}}=\frac{4}{3{\pi }^{1/2}}\left(R_{\text{sub}}-1\right){\xi }^3\exp \left(-{\xi }^2\right),$$(25)

where R_sub represents the sub-shock compression ratio, the ratio of upstream and downstream flow speed in the shock rest frame where V_sh is calculated by mass flux through the shock in the observer’s frame and u_u and downstream plasma velocity in observer’s frame (u_d) are calculated at the immediate upstream and downstream at the shock using resharpened flow velocity profile. We have used ξ = 4.2 in our simulations, which is consistent with the observed radiation flux from SN1006 (Brose et al. 2021, Appendix A) and conforms with the test-particle limit. Furthermore, as we have injected electrons and protons at the same ξ, the electron-to-proton ratio at higher energy can be determined by their mass ratio, $K_{ep}~\sqrt{m_e/m_p}$ (Pohl 1993).

3 Results

We follow the evolution of the SNRs of the 20 M_⊙, and 60 M⊙ progenitor stars for 30000 years and 110000 years, respectively. We discuss the behaviour of parameters for the SNR forward shock along with the spectra for accelerated particles. Furthermore, we compare the spectra, specifically at later times, obtained from the self-consistent calculations with explicit turbulence transport with those derived for Bohm-like diffusion, illustrated in Das et al. (2022). Additionally, we present the spectra and morphology of non-thermal emissions from the SNRs.

3.1 Shock parameters

Figure 3 illustrates the forward shock parameters for both SNRs.

3.1.1 SNR with 20 M_⊙ progenitorTEXT

During the expansion through free RSG stellar wind region, the speed of the forward shock gradually decreases from 6300km s⁻¹ to 5300 km s⁻¹, and the sub-shock compression ratio is approximately 4. When the forward shock starts interacting with the dense RSG shell, the shock speed falls to about 2000 km s⁻¹, and the sub-shock compression ratio slightly diverges from the value 4 and reaches 3.7. For brief moments around 1500 years and 3600 years the compression ratio becomes approximately 4.2, when the forward shock propagates along a steeply declining density, once after reaching the peak of RSG shell and then at the transition from the piled-up RSG wind to the shocked MS wind. During this period, the upstream density decreases more rapidly in comparison to the downstream density, which consequently results in a slightly higher compression ratio at the forward shock. Thus, the fluctuations in compression ratio reflect interactions with discontinuities in the CSM. The hot MS wind material in region 3 causes reduction in the sonic Mach number, Ms, and consequently, the compression ratio becomes approximately 3.5.

Fig. 4 Proton (PR) and electron (EL) spectra volume-averaged downstream of the forward shock at different regions of the corresponding wind bubble for 20 M⊙ progenitor (cf. Fig. 2).

3.1.2 SNR with 60 M_⊙ progenitorTEXT

The evolution of the forward shock parameters has been described in Das et al. (2022, Sect. 3.1). The most prominent feature is that the sub-shock compression ratio reaches approximately 1.5 while the SNR expands through the hot shocked wind.

The SNR shock interactions with any structures in CSM give rise to reflected shocks that propagate through the downstream of SNR forward shock and eventually hit other structures present there and reflected back in the outward direction. These shocks collide with SNR forward shock and produce instantaneous spikes in the shock velocity, as shown in Fig. 3. This feature has already been described elaborately in Das et al. (2022).

3.2 Particle acceleration and escape

To understand the spectral evolution of accelerated particles, the volume-averaged downstream spectra for protons and electrons are illustrated in Figs. 4 and 7 for the two types of SNRs, and the corresponding spectral indices are shown in Figs. 5 and 8, respectively. The selected ages correspond to the forward shock passing through the different regions of wind bubbles. The earlier study of type-Ia SNR with self-consistent amplification of Alfvénic turbulence (Brose et al. 2020) indicated that the spectral shape is controlled by the transport of turbulence and Alfvénic diffusion. The continuous deceleration of the shock decreases the CR flux and hence reduces the CR gradient, leading to a weaker driving of turbulence. For the core-collapse SNRs discussed here, the resulting particle spectra should get affected both by the complicated hydrodynamics of the wind bubbles and by the dynamics of the self-consistent turbulence. Further, the profiles resulting from the passive transport of CSM magnetic field (cf. Sect. 2.2.1) can also influence the particle acceleration.

Fig. 5 Variation of the spectral index for downstream protons at different ages with momentum.

3.2.1 SNR with 20 M_⊙ progenitorTEXT

In this section, we describe the particle acceleration in the different regions of the wind bubble. For the free RSG wind, during the propagation of the SNR forward shock inside this region, the compression ratio became four, so the particle spectral index was around two. At the end of the evolution inside this region, the maximum attainable energy for the protons reached 20 TeV.

For the piled-up RSG wind, during the interaction between the forward shock and dense RSG shell, the maximum attainable energy of protons slightly decreases to 6 TeV, on account of the decrease in shock speed. During this time, illustrated at 1100 years and 1500 years, spectral softening is observed, and the spectral index reaches 2.2 at energies above 100 GeV for protons, but the spectral index returns to 2 when the forward shock has emerged from the RSG shell. This brief spectral softening reflects the variation of the sub-shock compression ratio depicted in Fig. 3. After climbing the RSG shell, the shock speed again increases and so does the maximum energy of particles. During this time the compression ratio becomes 4.2, which results in spectral hardening to an index of 1.9 at higher energies, shown at 2500 years. Further, the collision between the piled-up RSG wind and the forward shock gives rise to an inward-moving reflected shock that merges with the SNR reverse shock and propagates towards the interior with a speed of ~3000 km s⁻¹. We do not inject particles at this shock, but this inward-moving shock can re-accelerate energetic particles. As we display the total volume-averaged particle spectra downstream of the forward shock, at higher energies the hardest contribution should dominate (Brecher & Burbidge 1972; Büsching et al. 2001).

For the shocked MS wind, in this region, the forward shock passed through hot MS wind material, and consequently the forward shock had a reduced compression ratio and produced relatively soft particle spectra, but these particles are so few that when volume-averaged over the downstream region, we did not observe a spectral softening (e.g. at 5000 years in Fig. 4). The high temperature in the shocked MS wind also provided a high injection momentum, and the volume-averaged spectrum was dominated by previously accelerated particles at all energies. The freshly accelerated particles in the MS wind could not penetrate deep into the interior on account of the strong magnetic field in the piled-up RSG wind.

Regarding the shocked ISM, after passing through the contact discontinuity between shocked MS wind and shock ISM, the forward shock became inefficient at accelerating particles up to a very high energy. In the shocked ISM region, the maximum attainable energy for protons falls from 20 TeV to 50 GeV at the end of our simulation. The forward shock encounters dense and cool material in this region that lowers the injection momentum and increases the injection rate into DSA. But the forward shock is too weak to accelerate these particles to very high energy. Hence, a prominent spectral break near 1 GeV is noticeable in Fig. 4 at 11 000 years and later. Later, the multiple merging of fast shocks with the weak forward shock increases the acceleration efficiency, and the step-like spectral feature emerges at slightly higher energies. The electron spectra roughly follow those of protons, except for an additional softening arising from the radiative cooling.

In our time-dependent treatment of the transport of CR and Alfvénic turbulence, the driving of turbulence diminishes at later stages of SNR evolution on account of the decrease in CR pressure gradient (Brose et al. 2020). Therefore, the diffusion coefficient as well as the acceleration timescale (Drury 1991) increase. Hence, particles with a higher energy in the shock precursor are slow to return and participate in further shock acceleration. A break in the downstream particle spectra should appear at the currently achievable maximum energy, above which particles escape to the far-upstream region. After 30 000 years, this softening is clearly visible in Fig. 5, where the spectral index reaches approximately 2.6 at high energy, starting from about 2.2 above 10 GeV.

Total production spectra of protons, including particles outside of the SNR, are depicted in Fig. 6 at different times. As we solve the CR transport equation out to 65 Rsh, all particles reside inside the simulation domain. Therefore, integrating the particle spectra over the whole simulation domain yields the total CR production of the SNR, and that is spectrally harder than the component inside the remnant. At 5000 years, protons are accelerated up to 20 TeV energy while at 11 000 years protons above 200 GeV are preferentially found in the upstream region, rendering the downstream spectra softer with the spectral index of ~2.2, shown in Fig. 6. At a later stage, the transition energy shifts to 30 GeV and the spectral index reaches ~2.6 above this escape energy. At any time, the break appears at the energy up to where the particles can currently be accelerated. It reflects the interplay between the reduction in maximum energy and the escape of particles from the accelerator (Ohira et al. 2010; Celli et al. 2019; Brose et al. 2020).

Fig. 6 Proton number-spectra at later times. Solid lines represent the total spectra, and dashed and dotted lines indicate the forward shock upstream and the downstream spectra, respectively. The grey vertical lines at 200 GeV and 30 GeV denote the escape energies at 11 000 years and 30 000 years, respectively. The arrows point out the energy bands with spectral indices of 2.2 and 2.6.

Fig. 7 Proton (PR) and electron (EL) spectra volume-averaged downstream of the forward shock at different times and regions of the corresponding wind bubble of the 60 M0 progenitor (cf. Fig. 2).

Fig. 8 Variation of the spectral index for protons at different ages with momentum.

3.2.2 SNR with 60 M_⊙ progenitorTEXT

Particle acceleration by the SNR forward shock at different regions of the wind bubble was explored in our previous work (Das et al. 2022) for Bohm-scaling of the diffusion coefficient. Here we investigate the spectral modifications deriving from the explicit treatment of turbulence transport, considering the CR streaming instability and the growth, damping, and cascading of the waves. The propagation of the SNR forward shock through the hot shocked regions, region 2 and region 3 in Fig. 2, leads to softened spectra below roughly 10 GeV, but not as severe as for Bohm scaling, in particular at higher energies. Later, after 28 000 years, the spectral index reaches 2.75 at lower energies, agreeing with the result for Bohm-like diffusion. Additionally, the introduction of turbulence adds non-linearity to the calculation of the diffusion coefficient. While the SNR shock expands through the shocked ISM, with the self-consistent turbulence model the spectral softening at higher energies arises from the evanescence of Alfvén waves. After 110 000 years, the spectral index reaches approximately 2.2 above ~ 10 GeV, illustrated in Fig. 9, on account of the escape of particles from the remnant. Further, the proton production spectra shown in Fig. 9 suggest that the particles above 70 GeV start leaving the remnant already after 28 000 years when the forward shock resides inside the shocked MS wind. Then, this typical energy of escaping CRs is reduced to about 1.5 GeV after 110 000 years as a consequence of the inefficient confinement of high-energy particles in the Alfvénic diffusion scenario.

At SNR with 20 M_⊙ progenitor the total proton production spectra have an index of ~2.1–2.2 above approximately 100 GeV energies at 11 000 years and 30 000 years. For the 60 M_⊙ progenitor, the index of the total proton production spectra is ~1.9 near the cut-off at later times. Both spectral shapes are quite different from the production index for type-Ia SNRs, about 2.4 above 10 GeV (Brose et al. 2020). There, the downstream spectra are shaped by the temporal decline of the maximum energy, E_max, to which the shock can accelerate (see also Celli et al. 2019). In contrast, the downstream spectra at core-collapse SNRs also reflect the effect of the CSM hydrodynamics, besides the time evolution of turbulence spectra. For example, the SNR with 60 M_⊙ provides downstream spectra with index ~1.9 at higher energies, before the escape starts. Although at later times the downstream spectra become softer as the consequence of particle escape, the proton production spectra show a slight hardness near cut-off, shown in Fig. 9.

Our results suggest that the cascading and decay of turbulence are crucial in the formation of soft particle spectra with spectral breaks for the older remnant. In contrast, with simple Bohm-like diffusion particles near the cut-off energy may escape the shock environment only at later stages, yielding insignificant spectral modifications, as shown in Fig. 10 for SNR with 20 M_⊙ progenitor as well as for 60 M_⊙ progenitor described by Das et al. (2022). We emphasise that the spectral shape with Bohm-like diffusion only reflects the influence of the CSM flow profiles, while in the Alfvénic scenario the spectra are a function of the hydrodynamic parameters and the properties of the turbulence. Our investigation on particle acceleration in SNRs with 20 M_⊙ and 60 M_⊙ progenitors illustrate the differences in spectral features arising from the environment of the SNRs. This is to be noted that the flow structure of wind bubbles becomes eventually very complex, on account of multiple reflected and transmitted shocks. We have solved the CR transport equation on a non-uniform grid centered on the forward shock to obtain the required fine resolution near the shock, described in Sect. 2.3. Simultaneously refining the grid also for the other shocks located downstream is desirable but out of the scope of this paper. Further, acceleration of particles from very low energies is possible only for the forward shock, so we injected the particles into the DSA only there. Re-acceleration of pre-energised particles at the other shocks, in particular the reverse shock, is included though, and we have seen negligible effects of that in the particle spectra. The higher fraction of reaccelerated particles in Das et al. (2022) is entirely caused by their choice of a Bohm-like diffusion coefficient. We obtain softer spectra at high energies for the SNR with 20 M_⊙ progenitor, during the propagation of the SNR forward shock in the shocked ISM, while for the very massive progenitor the signature of softness is observed already when the SNR shock is in the shocked wind, on account of particle escape. Here, the difference between a Bohm-like diffusion coefficient (Das et al. 2022) and the self-consistent treatment in this work becomes particularly evident, as the softness at low energies is present in both cases, due to the hot shocked wind, but at higher energies considerably softer spectra only develop due to particle escape. The enhanced driving of turbulence in the scenario with the 20 M_⊙ progenitor arises from the higher normalisation of the CR spectrum (Brose et al. 2020), that is caused by roughly four orders of magnitude higher density in the RSG wind. Additionally, the shock radius for the 20 M_⊙ progenitor is smaller than for the 60 M_⊙ progenitor, which enhances the gradient in the cosmic-ray distribution and hence the driving rate of turbulence. Figure 7 suggests that the maximum cut-off energy for protons is E_max ≈ 0.5 TeV for the 60 M_⊙ progenitor, which is significantly less than that for the 20 M_⊙ progenitor, E_max ≈ 50 TeV. A contributing factor could be our estimation of ejecta mass, M_ej = 11.75 M_⊙ for SNR with 60 M_⊙ progenitor, assuming the compact object as a neutron star. But, in reality, a black hole might be formed in the explosion of this very massive star, which would result in a lower ejecta mass. The ejecta speed and the shock speed would then be higher, leading to a maximum energy of particles that is higher by a factor of a few, but not more, because the acceleration rate scales with $v_{\text{sh}}^2$, and hence, it scales inversely with the ejecta mass for a given explosion energy. In addition, the result from SNR with 20 M_⊙ progenitor is consistent with the estimation of maximum energy for type II progenitor discussed in Cristofari et al. (2020). Further, the spectral index at high energies obtained for the SNR with 20 M_⊙ progenitor at later times is roughly comparable to those observationally deduced for IC443 and W44 (Ackermann et al. 2013).

Fig. 9 Proton number-spectra at later times. The solid lines represent the total spectra. The dashed and dotted lines indicate the spectrum upstream and downstream of the forward shock, respectively. The grey vertical lines at 70 GeV, 5 GeV, and 1.5 GeV denote the escape energies at 28 000 years, 44 000 years, and 110 000 years, respectively. The arrow points out the energy bands with a spectral index of 2.2.

Fig. 10 Proton number-spectra for Bohm-like (left panel) and self-consistent (right panel) diffusion.

Fig. 11 Synchrotron emission from the SNR at different ages. The upper boundaries of the shaded regions indicate the total emission from the remnant, while the lower boundaries denote emission downstream of the SNR forward shock. The brown band indicates the 50 MHz–10 GHz range and the blue band denotes 0.1 keV−40 keV.

3.3 Non-thermal emission

The processes of non-thermal radiation include synchrotron emission, inverse Compton scattering, and the decay of neutral pions. In our results, the synchrotron cut-off energy reflects the evolution of the electron acceleration efficiency that depends on the magnetic field profile and hydrodynamics. In this paper we only consider the cosmic microwave background to derive the inverse-Compton emission from the remnants. Porter et al. (2006) estimated that the contribution of infrared and optical photon fields to inverse-Compton emission is not significant to that of CMB photons except for SNRs residing near the Galactic centre. Hence, we are aware that these additional photon fields might be needed to consider during the modelling of specific SNRs by taking into account the respective background photon fields. The emission flux from different processes is calculated for remnants at 1 kpc distance. Figs. 11 and 16 illustrate the synchrotron spectra from the SNRs with 20 M_⊙ and 60 M_⊙ progenitors respectively, and Figures 12 and 17 demonstrate the gamma-ray spectra from both the remnants. These illustrations depict the time evolution of the total emission spectra, as well as the part originating downstream of the SNR shock. Figure 13 demonstrates as a function of time the radio flux at 5 GHz, the non-thermal X-ray flux in the 0.1 ke V−10 keV range, as well as the gamma-ray emission at high energy (0.1 GeV−100 GeV) and very high energy (>1 TeV) for the SNR with 20 M_⊙ progenitor.We also calculate intensity profiles at characteristic times for synchrotron (Figs. 14 and 18) and gamma-ray emission (Figs. 15 and 19), following the method described in Telezhinsky et al. (2013). Although total X-ray emission may be largely thermal during the interaction of the remnant with dense molecular clouds (Ustamujic et al. 2021a), efficient CR acceleration may suppress the post-shock temperature and hence the thermal X-ray flux (Drury et al. 2009). Calculating the thermal X-ray emission from SNRs is out of the scope of this study and only the non-thermal emission is presented.

Regarding the SNR with a 20 M_⊙ progenitor, below the age of about 5000 years, the entire non-thermal emission is essentially produced in the interior of the remnant. For the free RSG wind, at the early stages of evolution, the interaction of the SNR with the strong magnetic field in the free wind region, B₀ ∝ 1/r, combined with strong amplification, yields the considerable X-ray flux at this stage, shown at 500 years in Fig. 11. The simulated non-thermal X-ray light curve shown in Fig. 13 indicates declining X-ray emission during the SNR evolution inside free RSG wind. Initially, the remnant expands through dense material, ρ ∝ r⁻² in the free RSG wind, and consequently pion-decay emission dominates over inverse-Compton scattering for the first 200 years.

During this stage, the strongest magnetic field is found near the contact discontinuity between SNR forward and reverse shock (cf. Lyutikov & Pohl 2004). The turbulent magnetic field generated in the upstream and at the shock is eventually damped in the downstream region as a consequence of the weak driving of turbulence (Pohl et al. 2005). Hence, the contact discontinuity between forward and reverse shock of the SNR produces the highest intensity in the radio and X-ray band, shown in Fig. 14. Further, Fig. 15 indicates that at this stage, the peak intensity of pion-decay and inverse-Compton emission also coincides with the contact discontinuity between the forward and reverse shock. The interior of the remnant also appears inverse-Compton bright at 1 TeV, because the high-energy electrons can penetrate the deep interior.

Concerning the piled-up RSG wind, during the collision of the SNR with the RSG shell, at an age between 750 years and 1600 years, the radio spectra start to soften to a spectral index α ≈ 0.54 (S_v ∝ v^−α), well in agreement with the low-energy spectral index of particles, such as those seen in Fig. 5. But after the crossing of the RSG shell, the radio spectral index resumes to being α = 0.5.

In addition, the pion-decay flux enhances on account of the interaction with dense RSG shell. Further, during this interaction, the spectral index of pion-decay emission reaches 2.2. It is important to note here, that the softness of the pion-decay emission is extending all the way to low gamma-ray energies, unlike the all-downstream particle spectrum, as the bulk of the hadronic emission originates from the shocked RSG-shell immediately downstream of the shock. As a result, the low-energy gamma-ray spectra appear even softer than the radio emission that is produced by electrons of comparable energy. This stage of SNR evolution may be comparable to Cas A, which may be expanding inside the dense RSG wind (Chevalier & Oishi 2003). Although the ambient CSM of Cas A may differ from that of a 20 M_⊙ progenitor, and the reverse shock of Cas A is an efficient accelerator (Uchiyama & Aharonian 2008), the spectral index for accelerated protons obtained in this study is comparable with that estimated by Saha et al. (2014). The synchrotron flux in the radio and X-ray band as well as the hadronic emission flux reaches their maximum at this time when the forward shock crosses the dense RSG shell. For the next 200 years after that, the X-ray flux and high-energy pion-decay flux decrease at rates of ~0.3% yr⁻¹ and ~ 0.7% yr⁻¹, respectively, as a consequence of the declining density and the deceleration of the remnant.

For Cas A, a reduction rate of 1 .5% yr⁻¹ was observed for non-thermal X-rays in the band 4.2 keV−6keV (Patnaude et al. 2011). So, it is possible that the forward shock of Cas A is propagating through piled-up RSG wind, and the difference in density, age, and shock radius may reflect a much lower progenitor mass for Cas A than the 20 M_⊙ assumed here.

Further, the synchrotron morphology is centre-filled in this region. After the encounter with the RSG shell, the velocity of forward shock starts to increase, and the reverse shock and the contact discontinuity move towards the interior. After 2500 years, two radio shells are visible, the inner shell at the contact discontinuity and the outer one at the shocked RSG shell. Additionally, the brightest X-ray band is created near the forward shock. In the gamma-ray band, the brightest gamma-ray emission emanates from the region near the contact discontinuity between forward and reverse shock, and reverse shock. In reality, the intensity peak is likely smeared out, on account of the Rayleigh-Taylor instability of the contact discontinuity. Furthermore, at 2500 years the entire remnant looks bright in inverse-Compton emission, specifically at 1 TeV, as the energetic particles reach the reverse shock and may be re-energised there.

At 5000 years, when the SNR forward shock is inside this region and about to collide with the contact discontinuity between the shocked MS wind and shocked ISM, Fig. 11 indicates that the total X-ray synchrotron flux is slightly higher than that coming from the downstream region only. The intensity of X-ray emission at this time, shown in Fig. 14, indicates that the upstream emission is generated near or at the wind bubble contact discontinuity on account of the very strong magnetic field there, which is almost 15 times stronger than that in the shocked MS wind. At the same time, the highest radio intensity comes from the piled-up RSG wind behind the forward shock and also from the region inside of the contact discontinuity between the SNR forward and reverse shock. In this stage the pion-decay flux decreases on account of the lower density of the medium, and the very high energy inverse-Compton emission flux dominates over the high-energy flux, as the maximum achievable energy for electrons increases on account of the increase in shock velocity by 2000 km s⁻¹ in this region. The normalised intensity profile at 5000 years suggests that both the leptonic and hadronic gamma-ray emission emerges from deep downstream.

During the forward-shock passage through the shocked ISM, particles of higher energy can escape the remnant, as discussed in Sect. 3.2, and a significant fraction of the synchrotron flux is produced in the upstream region, depicted at 30 000 years in Fig. 11. Additionally, the pion-decay emission flux enhances on account of the shock propagating in dense material. At this age, a fraction of the gamma-ray emission is produced around the remnant, on account of particle escape. The spectral index for pion-decay emission reflects the soft proton spectra, and it reaches 2.4–2.6 above 10 GeV. This kind of soft gamma-ray spectra has been observed from SNRs (e.g. IC443, W44, and G 39.2−0.3) expanding in or near dense molecular clouds (Fang et al. 2013; Cardillo et al. 2014; de Oña Wilhelmi et al. 2020). The old remnant appears shell-like in pion-decay emission, whereas the inverse-Compton emission is centre-filled, contrary to the synchrotron intensity profile. This morphology of old core-collapse remnant is consistent with that of old type-Ia SNRs (Brose et al. 2021).

Regarding the SNR with a 60 M_⊙ progenitor, for self-consistent turbulence, the flux evolution of non-thermal emission in the different energy bands is similar to that for Bohm-like diffusion, which has been described in our previous study Das et al. (2022). Therefore, in this work, we only summarise the difference in emission spectra and morphology on account of self-consistenceturbulence andthe time-dependentevaluationof diffusion coefficient.

Considerable synchrotron and gamma-ray emission emerge from the upstream region of the SNR forward shock already during the propagation of SNR in the shocked wind, on account of the escaped particles arising from the weak driving of turbulence, as indicated by Fig. 9. The contact discontinuity between SNR forward and reverse shock appears bright in the radio and the X-ray band before the SNR approaches the shocked ISM, clearly visible in the X-ray morphology at 28 000 years, illustrated in Fig. 18. Around 28 000 years, following the collision of the SNR forward shock with the LBV shell, the wind-bubble contact discontinuity appears X-ray bright on account of the very strong magnetic field. The brightest radio emission still emerges from the region near the contact discontinuity between the SNR forward and reverse shock. The synchrotron morphology, illustrated in Fig. 18 modelled with self-consistent turbulence is similar to that for Bohm-like diffusion (cf. Das et al. 2022), except the contribution from escaped particles around the SNR. This is to be noted that the radio spectra are softer with spectral index, α ≈ 0.8, where energy flux, S_v = v^−α, during the shock passage through the shocked wind and the shocked ISM, whichis quite consistentwiththe observed indices for many old Galactic SNRs (Baars et al. 1977; Green 2009; Urošević 2014).

During the early stages of evolution, the strong magnetic field in the free wind permits efficient proton acceleration to very high energy, and the gamma-ray emission is predominantly hadronic, independent of the diffusion model. At later stages, when the SNR expands in the shocked wind, inverse-Compton emission dominates and pion-decay emission is diminished until particle escape becomes significant. Around 28 000 years, high-energy protons residing in the shock precursor can reach the high-density material behind the contact discontinuity of the wind bubble, which causes bright hadronic TeV-scale emission from the periphery of the SNR. Therefore, we obtained significant pion-decay emission flux from the upstream of the remnant, demonstrated in Fig. 17 at 28 000 years and also visible in the morphology at 1 TeV, shown in Fig. 19 where the pion-decay emission emerges from the shell around the wind-bubble contact discontinuity ahead of the SNR shock. This situation is comparable with the scenario of gamma-ray emission from the interaction of escaped protons with ambient molecular clouds, that has been suggested for IC 443, W44, G 39.2 − 0.3 and G 106.3 + 2.7 (Fang et al. 2013; Yang et al. 2022; de Oña Wilhelmi et al. 2020). Pion-decay dominates the gamma-ray emission from old remnants whose forward shock propagates in the shocked ISM, and the spectral index (2.2–2.4) for pion-decay emission above 10 GeV reflects the softness of the proton spectra.

The gamma-ray emission morphology is shell-like at early stages, dominated by the emission from the contact discontinuity between forward and reverse shock. After the collision with the wind bubble contact discontinuity, the velocity of SNR forward shock decreases, and it can no longer accelerate particles at very high energies. At that time, confined particles of very high energy can reach, and be re-accelerated by, the reverse shock that itself is energised by multiple reflected shocks resulting from the SNR-CSM interaction. It is already evident that after 28 000 years, the inverse-Compton emission mainly emerges from the region around the reverse shock, and the morphology eventually becomes centre-filled. This late-time morphology is similar to that obtained for the SNR with 20 M_⊙ progenitor.

Fig. 12 Gamma-ray emission by pion-decay (PD) and inverse Compton (IC) scattering at different ages. The boundaries of the shaded regions indicate the total emission and that from the interior, as in Fig. 11.

Fig. 13 Evolution of the energy flux (Φ) during the lifetime of the SNR for synchrotron emission and gamma-ray emission at specific energy ranges inside the different regions shown in Fig. 2 of the wind-bubble formed by a 20 M⊙ progenitor. We specify that the shown X-ray flux only refers to the non-thermal X-ray.

Fig. 14 Normalised intensity profiles for synchrotron emission. For each specified profile, the intensity has been normalised to its peak value.

Fig. 15 Normalised intensity profiles for gamma-ray emission. Emissions for pion-decay (PD) and inverse-Compton emissions (IC) are shown. For each specified profile, the intensity has been normalised to its peak value.

Fig. 16 Synchrotron spectra from the SNR at different ages. The boundaries of the shaded regions indicate the total emission and that from the interior of the SNR, as in Fig. 11. The brown band indicates the 50 MHz−10 GHz range, and the blue band denotes 0.1 keV−40 keV.

Fig. 17 Gamma-ray emission by pion-decay (PD) and inverse Compton (IC) scattering at different ages. The boundaries of the shaded regions indicate the internal and total emission, as in Fig. 11.

Fig. 18 Normalised intensity profiles for synchrotron emission. For each specified profile, the intensity has been normalised to its peak value.

Fig. 19 Normalised intensity profiles for gamma-ray emission. Emissions for pion-decay (PD) and inverse Compton (IC) emissions are shown. For each specified profile, the intensity has been normalised to its peak value.

4 Conclusions

We have probed the structures of the ambient medium at the pre-supernova stage for 20 M_⊙, and 60 M_⊙ progenitors by solving the hydrodynamic evolution of the CSM from the ZAMS stage to pre-supernova stage. We have explored the interaction of the SNRs with the so-modelled CSM as well as the acceleration and transport of energetic particles. For that purpose, we solved the time-dependent transport equations of CRs in the test-particle limit and of magnetic turbulence as well as the induction equation for the large-scale magnetic field, all in parallel with the hydrodynamic equations for the SNR and the CSM. The self-consistent turbulence module provides a time-and momentum-dependent spatial CR diffusion coefficient that is far more realistic than the oversimplified Bohm-scaling of the diffusion coefficient.

We have demonstrated that inefficient confinement of high-energy particles eventually causes spectral breaks at Giga-electron volt energies, above which the spectral index is 2.2–2.6. The simulations of particle acceleration at the SNR forward shock indicate that the spectra of the particles and their emission products are significantly affected by the structure of the wind bubble. Our morphological analysis of two SNRs with 20 M_⊙ and 60 M_⊙ progenitors has revealed dissimilarities in various frequency bands that reflect the differences in the wind bubbles.

In the scenario with the 20 M_⊙ star, we observed softer transient spectra with a spectral index of 2.2 specifically at higher energy during the collision between the SNR-RSG shell. Beyond this, the flow profiles of the ambient medium did not induce any spectral softness. It is the weak driving of turbulence that for old remnants renders proton spectral indices of around 2.6 at high energies. Similar soft proton spectra have been deduced from observations of the SNRs evolving inside dense molecular clouds. The synchrotron flux depends on the total magnetic field intensity in the different regions of the wind bubble as well as the maximum achievable electron energy. At the later stage, although the magnetic field strength is very high in the shocked ISM, the paucity of high-energy electrons in the shock environment causes the synchrotron cut-off energy to decline from 10 keV to 0.1 keV. The flux of pion-decay emission varies throughout the evolution, whereas the inverse Compton flux shows a steady trend and slightly increases until the SNR enters the shocked ISM. The SNRs with a 20 M_⊙ progenitor have a shell-like morphology in the X-ray band and pion-decay emission, except at middle age, while in the radio band, the SNR appears more centre-filled, except at a very old age. The inverse-Compton morphology is that of a thick shell, and it transitions to a more centre-filled configuration at later times.

For the SNR with a 60 M_⊙ progenitor star, we obtained soft spectra also at low energies, on account of the low sonic Mach number inside the hot shocked wind regions. At later times, the spectra become soft at higher energies as well, reflecting inefficient driving of turbulence and the associated rapid decline in the currently achievable maximum particle energy, which we already described for the SNR with a 20 M_⊙ progenitor. The synchrotron cut-off frequency increases with time on account of the high magnetic field intensity in the shocked wind until the high-energy electrons start to escape from the remnant. Inverse-Compton emission dominates the gamma-ray output until the SNR reaches the vicinity of the contact discontinuity of the wind bubble, and the pion-decay emission is prominent for the old remnant. For both types of SNRs, the gas density of the ambient medium determines the dominant contribution in the gamma-ray band (cf. Yuan et al. 2012). The X-ray morphology resembles a thick shell, whereas the radio emission evolves with time from having a shell-like to a centre-filled configuration. In the gamma-ray band, the pion-decay intensity profile is shell-like, and that of inverse-Compton emission eventually transitions from shelllike to a centre-filled structure. For both SNRs, the morphology looks similar for old remnants, and in this stage, it also resembles that of type-Ia SNRs (Brose et al. 2021).

In conclusion, it is evident that the spectra of accelerated particles are shaped by both the hydrodynamics of the ambient medium and the time-dependent diffusion coefficients. Our results suggest that non-thermal emission and its morphology can provide information about the progenitor stars and the current state of evolution of the remnant, at least until it reaches the shocked ISM. An SNR with a lower-mass progenitor star (20 M_⊙) is more likely to be detected with the current-generation observations on account of the high density of the RSG wind.

Acknowledgements

The authors acknowledge the North-German Supercomputing Alliance (HLRN) for providing HPC resources to perform the hydrodynamic simulation from ZAMS to pre-supernova stage of the star. SD acknowledges support from the Collaborative Research Center SFB1491 Cosmic Interacting Matters – From Source to Signal, funded by the Deutsche Forschungsgemeinschaft DFG under the grant number 445052434.

References

All Figures

Fig. 1 Evolution of stellar mass (panel a), mass-loss rate (panel b), and wind velocity (panel c) as a function of stellar age (in Myr) during different evolutionary stages of the 20 M⊙ and 60 M⊙ stars. The evolutionary stages of the 20 M⊙ star are shown with red dashed lines and those of the 60 M⊙ star are shown with blue solid lines. We used these parameters as the initial condition in the hydrodynamic simulation for constructing the corresponding CSM at the pre-supernova stage at 8.64 million years for the 20 M⊙ star and 3.95 million years for the 60 M⊙ star. The stellar parameters for the 20 M⊙ and 60 M⊙ stars are taken from Ekström et al. (2012); Groh et al. (2014), respectively.

In the text

Fig. 2 Pre-supernova CSM profiles of the number density (n; panel a), the flow speed (u; in panel b), the thermal pressure (P; panel c), and the temperature (T; panel d) at 8.64 million years after the ZAMS stage for the 20M⊙ progenitor (on the left) and at 3.95 million years for the 60M⊙ progenitor after the ZAMS stage (on the right). The right panels are reproduced from Das et al. (2022). Vertical grey lines mark the boundaries of specific regions. For the 20M⊙ progenitor, these are the free RSG stellar wind (region 1), piled-up RSG wind (region 2), shocked MS wind (region 3), shocked ISM (region 4), and ISM (region 5). The RSG shell indicates the accumulation of decelerated, freely expanding RSG wind; RRSG denotes the transition between RSG and MS wind; and CD represents the contact discontinuity between shocked wind and the ISM. For the 60M⊙ progenitor, we distinguish the free stellar wind (region 1), the shocked LBV and WR wind (region 2), the shocked MS wind (region 3), the ISM (region 4), and the ambient ISM (region 5). The term RWT is the radius of wind termination shock, LBV shell denotes the dense shell created by the interaction between LBV wind and WR wind, and CD represents the contact discontinuity between the shocked wind and shocked ISM.

In the text

	Fig. 3 Temporal evolution of the forward shock parameters. The parameters include radius, Rsh; speed, Vsh; and the sub-shock compression ratio, Cr. In the upper panel, we provide the angular scale, θsh, for an SNR located 1000 pc away. We also mark the interactions of the forward shocks with discontinuities in the wind bubbles of the two types of progenitor stars.
In the text

	Fig. 4 Proton (PR) and electron (EL) spectra volume-averaged downstream of the forward shock at different regions of the corresponding wind bubble for 20 M⊙ progenitor (cf. Fig. 2).
In the text

	Fig. 5 Variation of the spectral index for downstream protons at different ages with momentum.
In the text

	Fig. 6 Proton number-spectra at later times. Solid lines represent the total spectra, and dashed and dotted lines indicate the forward shock upstream and the downstream spectra, respectively. The grey vertical lines at 200 GeV and 30 GeV denote the escape energies at 11 000 years and 30 000 years, respectively. The arrows point out the energy bands with spectral indices of 2.2 and 2.6.
In the text

	Fig. 7 Proton (PR) and electron (EL) spectra volume-averaged downstream of the forward shock at different times and regions of the corresponding wind bubble of the 60 M0 progenitor (cf. Fig. 2).
In the text

	Fig. 8 Variation of the spectral index for protons at different ages with momentum.
In the text

Fig. 9 Proton number-spectra at later times. The solid lines represent the total spectra. The dashed and dotted lines indicate the spectrum upstream and downstream of the forward shock, respectively. The grey vertical lines at 70 GeV, 5 GeV, and 1.5 GeV denote the escape energies at 28 000 years, 44 000 years, and 110 000 years, respectively. The arrow points out the energy bands with a spectral index of 2.2.

In the text

	Fig. 10 Proton number-spectra for Bohm-like (left panel) and self-consistent (right panel) diffusion.
In the text

	Fig. 11 Synchrotron emission from the SNR at different ages. The upper boundaries of the shaded regions indicate the total emission from the remnant, while the lower boundaries denote emission downstream of the SNR forward shock. The brown band indicates the 50 MHz–10 GHz range and the blue band denotes 0.1 keV−40 keV.
In the text

	Fig. 12 Gamma-ray emission by pion-decay (PD) and inverse Compton (IC) scattering at different ages. The boundaries of the shaded regions indicate the total emission and that from the interior, as in Fig. 11.
In the text

	Fig. 13 Evolution of the energy flux (Φ) during the lifetime of the SNR for synchrotron emission and gamma-ray emission at specific energy ranges inside the different regions shown in Fig. 2 of the wind-bubble formed by a 20 M⊙ progenitor. We specify that the shown X-ray flux only refers to the non-thermal X-ray.
In the text

	Fig. 14 Normalised intensity profiles for synchrotron emission. For each specified profile, the intensity has been normalised to its peak value.
In the text

	Fig. 15 Normalised intensity profiles for gamma-ray emission. Emissions for pion-decay (PD) and inverse-Compton emissions (IC) are shown. For each specified profile, the intensity has been normalised to its peak value.
In the text

	Fig. 16 Synchrotron spectra from the SNR at different ages. The boundaries of the shaded regions indicate the total emission and that from the interior of the SNR, as in Fig. 11. The brown band indicates the 50 MHz−10 GHz range, and the blue band denotes 0.1 keV−40 keV.
In the text

	Fig. 17 Gamma-ray emission by pion-decay (PD) and inverse Compton (IC) scattering at different ages. The boundaries of the shaded regions indicate the internal and total emission, as in Fig. 11.
In the text

	Fig. 18 Normalised intensity profiles for synchrotron emission. For each specified profile, the intensity has been normalised to its peak value.
In the text

	Fig. 19 Normalised intensity profiles for gamma-ray emission. Emissions for pion-decay (PD) and inverse Compton (IC) emissions are shown. For each specified profile, the intensity has been normalised to its peak value.
In the text