Rectangular core-collapse supernova remnants: application to Puppis A

ABSTRACT

Core-collapse supernova remnants are the gaseous nebulae of galactic interstellar media (ISM) formed after the explosive death of massive stars. Their morphology and emission properties depend both on the surrounding circumstellar structure shaped by the stellar wind–ISM interaction of the progenitor star and on the local conditions of the ambient medium. In the warm phase of the Galactic plane (⁠|$n\approx 1\, \rm cm^{-3}$|⁠, |$T\approx 8000\, \rm K$|⁠), an organized magnetic field of strength |$7\, \mu \rm G$| has profound consequences on the morphology of the wind bubble of massive stars at rest. In this paper, we show through 2.5D magnetohydrodynamical simulations, in the context of a Wolf–Rayet-evolving |$35\, \rm M_{\odot }$| star, that it affects the development of its supernova remnant. When the supernova remnant reaches its middle age (⁠|$15\!-\!20\, \rm kyr$|⁠), it adopts a tubular shape that results from the interaction between the isotropic supernova ejecta and the anisotropic, magnetized, shocked stellar progenitor bubble into which the supernova blast wave expands. Our calculations for non-thermal emission, i.e. radio synchrotron and inverse-Compton radiation, reveal that such supernova remnants can, due to projection effects, appear as rectangular objects in certain cases. This mechanism for shaping a supernova remnant is similar to the bipolar and elliptical planetary nebula production by wind–wind interaction in the low-mass regime of stellar evolution. If such a rectangular core-collapse supernova remnant is created, the progenitor star must not have been a runaway star. We propose that such a mechanism is at work in the shaping of the asymmetric core-collapse supernova remnant Puppis A.

stars: evolution, stars: massive, ISM: supernova remnants, methods: MHD

1 INTRODUCTION

Core-collapse supernova remnants are the result of the explosive death of massive stars (⁠|$\ge 8\, \rm M_{\odot }$|⁠). These astrophysical events are rare. However, the mass and energy release in supernova explosions triggers a feedback process of substantial importance in the behaviour and evolution of the interstellar medium (ISM). This feedback enriches the ISM with heavy elements and drives turbulence, regulating future star formation processes (Hopkins, Quataert & Murray 2011). Core-collapse supernova remnants result from the propagation of the energetic blastwave produced by supernova explosions into the surroundings of a high-mass star, which prior to the explosion shaped its environment as structured, concentric layers of hot and cold gas, either from the stellar wind or from the ISM, separated by shocks (Woosley & Weaver 1986; Smartt 2009). Hence, the morphology of supernova remnants is a function of both the density distribution of the circumstellar medium (partly reflecting the past evolution of the progenitor star) and the local condition of the ISM. Therefore, studying the morphology of supernova remnants is key to understand stellar evolution, to constrain local properties of the ISM, and to examine the physics of particle acceleration at shocks (Langer 2012). The classic picture for core-collapse supernova remnants can be found in Chevalier (1982), Chevalier & Liang (1989), and Truelove & McKee (1999).

The vast variety of morphologies found in core-collapse supernova remnants is a mystery. Most of them substantially deviate from the spherically symmetric solution (Aschenbach & Leahy 1999; Arias et al. 2019a, b; Domček et al. 2022; Zhou et al. 2022). Many reasons can be invoked to explain this observational trend, like asymmetric stellar winds, or the presence of inhomogeneities in the surrounding medium. Both cases can induce a non-spherical expansion of the supernova shock wave. Indeed, the life span of the progenitor star, the successive stellar evolutionary phases, and the properties of the corresponding stellar winds depend on the star’s initial mass (Maeder & Meynet 2000). The stellar rotation (Szécsi et al. 2022) and the existence of one or more companions (Sana et al. 2012) are additional factors that strongly modify the structure of the stellar wind bubble and, consequently, that of the supernova remnant. The density distribution of the local, primeval ISM, such as the presence of a dense molecular component close to the site of the explosion, is another source of anisotropy prone to affect the propagation of the supernova blastwave (Fesen et al. 2018; Boumis et al. 2022).

Over the past decades, the wind bubble problem has spurred the interest of the numerical astrophysical community (Garcia-Segura, Mac Low & Langer 1996a; Garcia-Segura, Langer & Mac Low 1996b; Freyer, Hensler & Yorke 2003, 2006; Dwarkadas 2005; Eldridge et al. 2006). Moreover, the large morphological variety of supernova remnants provided a natural continuation of that computational effort, both in the high-mass (Dwarkadas 2007; van Veelen et al. 2009) and low-mass (Chiotellis, Schure & Vink 2012; Chiotellis et al. 2013; Broersen et al. 2014; Chiotellis, Boumis & Spetsieri 2020, 2021) regime of the progenitors. Many processes have been identified as sources of asymmetry in the circumstellar medium of massive stars, and consequently in their subsequent supernova remnants. Amongst others, the role of radiation in shaping the surroundings of massive stars is investigated in Toalá & Arthur (2011) and Geen et al. (2015). The runaway motion of a fraction of massive stars (Blaauw 1961; Gies 1987; Moffat et al. 1998) distorts wind bubbles into bow shocks (Gull & Sofia 1979; Wilkin 1996; Baalmann et al. 2020, 2021, 2022; Scherer et al. 2020; Herbst et al. 2022), naturally providing an asymmetric environment to the expanding shock wave and generating asymmetric supernova remnants (Franco et al. 1991; Rozyczka et al. 1993; Brighenti & D’Ercole 1995a, b; Fang et al. 2013; Fang, Yu & Zhang 2017; Zhang, Tian & Wu 2018; Meyer et al. 2020b; Yang, Bao & Liu 2020; Meyer 2021), themselves able to govern the morphology of pulsar wind nebulae (Meyer & Meliani 2022). Recently, Das et al. (2022) noted a softening of the non-thermal particle spectra when the shock wave goes through shocked circumstellar material. Similarly, the role of material expelled throughout pulsating mass-loss events occurring during the last century of the progenitor’s life has been highlighted in the context of the asymmetries in the expanding shock wave of Cassiopea A (Orlando et al. 2021, 2022). More generally, the lengthening effects of an organized magnetic field on the morphology of stellar wind bubble was analysed by van Marle, Meliani & Marcowith (2015).

This study is a first step in tackling the problem of the particular appearance of the core-collapse supernova remnant Puppis A, which exhibits a rhomboid shape, inside of which several filaments are arranged orthogonal to each other, providing the overall appearance of nested squares and rectangles. This specific arrangement is observable at many wavebands, spanning from the infrared (Arendt et al. 2010) and the optical/X-rays regimes (Hwang, Flanagan & Petre 2005) to non-thermal radio emission (Hewitt et al. 2012). We employ magnetohydrodynamical (MHD) simulations to explore the single O-type progenitor scenario for Puppis A proposed in Reynoso, Cichowolski & Walsh (2017). In our chosen sequence of events, the supernova explosion takes place within a circumstellar medium that has been self-consistently pre-shaped by the previous release of mass and energy along the successive stellar evolutionary phases (main-sequence, red supergiant, and Wolf–Rayet) into the ISM.

Our work is organized as follows. In Section 2, we review the numerical methods and the initial conditions employed to perform numerical MHD simulations for the circumstellar medium and subsequent supernova remnant of a massive progenitor at rest in a magnetized ISM. The asymmetries in the supernova remnant are directly obtained as a result of the interaction of the supernova shock wave with the pre-shaped magnetized circumstellar medium of the progenitor star (van Marle et al. 2015). The evolution of the remnant is presented in Section 3. We investigate the non-thermal and thermal appearance of our supernova remnant in Section 4, and further discuss our findings in the context of the core-collapse supernova remnant Puppis A. Finally, we summarize our conclusions in Section 5.

2 METHOD

In this section, we provide the reader with a description of the numerical methods utilized to perform simulations of asymmetric core-collapse supernova remnants of rectangular appearance. We also briefly summarize the radiative transfer recipes for non-thermal emission used in the analysis of our simulated models.

2.1 Circumstellar medium

First, the circumstellar medium of a massive star is simulated. The 2.5D models are carried out in cylindrical coordinates (R, z, θ) with a static grid, under the simplifying assumption of z-axis symmetry. The dimensions of the grid is [O, R_max] × [z_min, z_max] × [0, 2π], with a mesh that is uniform along the R- and z-directions. It is constructed of N_R = 2500 and N_z = 5000 grid zones and it extends to |$R_{\rm max}=z_{\rm max}=150\, \rm pc$| and |$z_{\rm min}=-150\, \rm pc$|⁠. The uniform resolution of the mesh is therefore Δ = |z_max − z_min|/N_z = R_max/N_R. Reflective conditions are assigned along the z-axis (i.e. at the boundary R = 0) and outflow conditions at the other boundaries, respectively. The computational domain is initially filled with gas having the properties of the warm phase of the ISM in the Galactic plane of the Milky Way, with number density |$n_{\rm ISM}\approx 0.79\, \rm cm^{-3}$| and |$T_{\rm ISM}\approx 8000\, \rm K$|⁠. The magnetization of the ISM is imposed by the geometry of the domain with magnetic field parallel to the symmetry axis Oz.

Our choice of an organized, parallel magnetic field is justified, since even though the typical description of the ISM is that of a turbulent plasma (Passot, Vazquez-Semadeni & Pouquet 1995; Zank, Nakanotani & Webb 2019), a large-scale coherent magnetic field is present in the spiral disc galaxies (see models and observations of Baryshnikova et al. 1987; Yusef-Zadeh & Morris 1987; Beck 2007; Brown et al. 2007; Jansson & Farrar 2012a). This large-scale field is responsible for, e.g. the alignment of supernova remnants with the enrolled arms of the Milky Way (Gaensler 1998) and affects the development of stellar wind bubbles of massive stars (van Marle et al. 2015). The turbulent component of the magnetic field has a significant strength overall (Jansson & Farrar 2012b), but it includes fluctuations with wavelength up to about 100 pc with an energy density that falls off with wavenumber, k. Radiation modelling of young supernova remnants suggests that the coherence length of the magnetic field can indeed be large compare to the size scale of our simulations (Pais et al. 2020). Then the random magnetic field on small scales acts as an additional isotropic, negligible pressure to the ambient medium in which circumstellar shocks propagate.

The stellar wind is imposed in a sphere of radius r_in = 20 grid zones, following the method of Comerón & Kaper (1998) (see also Mackey et al. 2012; Meyer et al. 2014). The stellar wind density is imposed assuming the profile

$$\begin{eqnarray*} \rho _{w} = \frac{ \dot{M} }{ 4\pi r^{2} v_{\rm w} }, \end{eqnarray*}$$

(1)

where |$r=\sqrt{R^2+z^2}$| is the radial distance from the pre-supernova location, |$\dot{M}$| the mass-loss rate, and v_w the wind terminal velocity. The mass-loss rate |$\dot{M}$| is time-dependently interpolated from the stellar evolutionary track of Ekström et al. (2012); the wind velocity is obtained from the escape velocity using the law of Eldridge et al. (2006). This star successfully evolves through a long hot main-sequence phase, a short cold red supergiant phase and finishes its life via a hot Wolf–Rayet phase. The main properties of the star are presented in great details in Meyer et al. (2021a). The evolution of the wind–ISM interaction is followed throughout the entire star life, from the onset of the main-sequence phase at t = 0 to the pre-supernova time at |$t=5.41\, \rm Myr$|⁠.

We perform two simulation models of similar ISM and stellar wind properties, the difference being the strength of the ISM magnetic field, taken to be |$B_{\rm ISM}=0\, \mu \rm G$| (hydrodynamical model) and |$B_{\rm ISM}=7\, \mu \rm G$| (MHD model), see Table 1.

Table 1.

List of models in this study. All simulations assume a non-rotating static |$35\, \rm M_{\odot }$| star at solar metallicity. The table indicates for each simulation their label, the number of grid points used, and the strength of the ISM magnetic field (in |$\mu \rm G$|⁠).

Model	Grid mesh	BISM (⁠|$\mu \, \rm G$|⁠)	Description
Run-35-HD-0-CSM	2500 × 5000	0	Circumstellar medium in unmagnetized ISM
Run-35-MHD-0-CSM	2500 × 5000	7	Circumstellar medium in magnetized ISM
Run-35-HD-0-SNR	3000 × 6000	0	Supernova remnant in unmagnetized ISM
Run-35-MHD-0-SNR	3000 × 6000	7	Supernova remnant in magnetized ISM

2.2 Supernova remnant

Once we know the structure of the circumstellar medium at the pre-supernova time, we calculate the initial ejecta–wind interaction in spherical symmetry in the radial range [O;r_out]. The ejecta material fills the inner r ≤ r_max < r_out region of the domain. The corresponding density field is defined as a constant density,

$$\begin{eqnarray*} \rho _{\rm core} = \frac{1}{ 4 \pi n } \frac{ \left(10 E_{\rm ej}^{n-5}\right)^{-3/2} }{ \left(3 M_{\rm ej}^{n-3}\right)^{-5/2} } \frac{ 1}{t_{\rm max}^{3} }, \end{eqnarray*}$$

(2)

within r ≤ r_core, and a power-law density profile

$$\begin{eqnarray*} \rho _{\rm max}(r) = \frac{1}{ 4 \pi n } \frac{ \left(10 E_{\rm ej}^{n-5}\right)^{(n-3)/2} }{ \left(3 M_{\rm ej}^{n-3}\right)^{(n-5)/2} } \frac{ 1}{t_{\rm max}^{3} } \bigg (\frac{r}{t_{\rm max}}\bigg)^{-n}, \end{eqnarray*}$$

(3)

in the r_core < r ≤ r_max region, respectively, with n = 11 (see Chevalier 1982; Truelove & McKee 1999).

The ejecta speed follows homologous expansion, v(r) = r/t, and

$$\begin{eqnarray*} v_{\rm core}(r_{\rm core}) = \bigg (\frac{ 10(n-5)E_{\rm ej} }{ 3(n-3)M_{\rm ej} } \bigg)^{1/2}. \end{eqnarray*}$$

(4)

The forward shock of the supernova blastwave at r_max has a flow speed of |$v_{\rm max}=30\, 000\, \rm km\, \rm s^{-1}$|⁠, which constrains the time of the simulation start t_max = r_max/v_max. The value of r_max is calculated using the method described in van Veelen et al. (2009). The value of r_max is arbitrary as long as the mass it embeds is smaller than that of the ejecta, and t_max is determined assuming |$v_{\rm max}=30\, 000\, \rm km\, \rm s^{-1}$|⁠. In equations (2)–(4), |$E_{\rm ej}=5 \times 10^{50}\, \rm erg$| is the total energy released by the supernova explosion, and |$M_{\rm ej}=10.12\, \rm M_{\odot }$| is the ejecta mass, calculated as the pre-supernova stellar mass minus a neutron star mass of |$M_{\mathrm{NS}}=1.4\, \rm M_{\odot }$|⁠. Finally, the 1D ejecta–wind interaction is mapped on to a 2.5D computational domain before the forward shock of the blastwave runs into the termination shock of the stellar wind bubble, see method in Meyer, Petrov & Pohl (2020a) and Meyer et al. (2021a). We assume that supernova ejecta and stellar wind are not magnetized.

2.3 Governing equations

The time evolution of the gas is calculated within the frame of ideal magnetohydrodynamics (Meyer et al. 2020b), described by

$$\begin{eqnarray*} \frac{\partial \rho }{\partial t} + \boldsymbol {\nabla } \cdot \big (\rho \boldsymbol {v}) = 0, \end{eqnarray*}$$

(5)

$$\begin{eqnarray*} \frac{\partial \boldsymbol {m} }{\partial t} + \boldsymbol {\nabla } \cdot \Big (\boldsymbol {m} {\otimes } \boldsymbol {v} {-} \boldsymbol {B} {\otimes } \boldsymbol {B} + \boldsymbol {\hat{I}}p_{\rm t} \Big) \boldsymbol {0}, \end{eqnarray*}$$

(6)

$$\begin{eqnarray*} \frac{\partial E }{\partial t} + \boldsymbol {\nabla } \cdot \Big ((E+p_{\rm t})\boldsymbol {v}-\boldsymbol {B}(\boldsymbol {v}\cdot \boldsymbol {B}) \Big) = \Phi (T,\rho), \end{eqnarray*}$$

(7)

and

$$\begin{eqnarray*} \frac{\partial \boldsymbol {B} }{\partial t} + \boldsymbol {\nabla } \cdot \Big (\boldsymbol {v} {\otimes } \boldsymbol {B} - \boldsymbol {B} {\otimes } \boldsymbol {v} \Big) = \boldsymbol {0}, \end{eqnarray*}$$

(8)

which stand for the mass continuity, momentum conservation, energy conservation equations, and for the evolution of the magnetic field, respectively. In the above relations, |$\boldsymbol {m}=\rho \boldsymbol {v}$| is the vector momentum, ρ the gas density, v the gas velocity, p_t the total pressure, and

$$\begin{eqnarray*} E = \frac{p}{(\gamma - 1)} + \frac{ \boldsymbol {m} \cdot \boldsymbol {m} }{2\rho } + \frac{ \boldsymbol {B} \cdot \boldsymbol {B} }{2}, \end{eqnarray*}$$

(9)

the total energy of the gas. The definition of the sound speed |$c_{\rm s} = \sqrt{ \gamma p/\rho }$| closes the system of equations, with γ = 5/3 the adiabatic index. Optically thin radiative cooling and photoheating Φ(T, ρ) are included into the equations via the prescriptions of Meyer et al. (2014), where T is the gas temperature.

The equations are solved using the pluto code¹ (Mignone et al. 2007, 2012; Vaidya et al. 2018). This study uses the Godunov-type solver utilized in Meyer et al. (2021b), made of the shock-capturing HLL Rieman solver (Harten, Lax & van Leer 1983) and the eight-wave finite-volume algorithm (Powell 1997) which ensures that the magnetic field remains divergence-free. The numerical method is made of a third-order Runge–Kutta time integrator, used together with the minmod flux limiter and the WENO3 interpolation scheme. The time-step is controlled by the Courant–Friedrich–Levy, initialized to C_cfl = 0.1 (Meyer 2021).

2.4 Non-thermal emission maps

For further comparison between our MHD models and available observational data, we perform radiative transfer calculations of selected time instances of the supernova remnant evolution, for non-thermal emission such as synchrotron and inverse-Compton emission. Regarding the synchrotron radio emission, we use the method detailed in Meyer et al. (2021a). It first assumes a spectrum of the non-thermal electrons present in the post-shock region of the propagating supernova blastwave,

$$\begin{eqnarray*} N(E) = K E^{-s}, \end{eqnarray*}$$

(10)

where s = 2, K ∝ n with n the gas number density, and with E the electron energy. The emission coefficient at frequency ν then reads

$$\begin{eqnarray*} j_{\rm sync}(\nu) \propto n B_{\perp }^{ (s+1)/2 } \nu ^{ -(s-1)/2 }, \end{eqnarray*}$$

(11)

where B_⊥ stands for the magnetic field component perpendicular to the observer’s line of sight, l. The viewing angle θ_obs is the inclination angle of the supernova remnant with respect to the plane of the sky. The emission maps are calculated using the radiative transfer code radmc-3D² which permits ray-tracing integration of arbitrary emission coefficients, providing non-thermal emission maps of synchrotron intensity,

$$\begin{eqnarray*} I = \int _{\rm SNR} j_{\rm sync}\left(\theta _{\rm obs}\right) \mathrm{ d}l. \end{eqnarray*}$$

(12)

Similarly, we calculate emission maps for inverse-Compton emission by assuming that the electron spectrum (equation 10) extends to a sharp cut-off at a very high energy E_max. We make use of the emission coefficient at a given photon energy ϵ,

$$\begin{eqnarray*} j_{\rm IC} = \int _{0}^{E_{\rm max}} N(E) \Lambda (E) \mathrm{ d}E \propto n \epsilon ^{-(s-1)/2}, \end{eqnarray*}$$

(13)

with Λ(E) an integral that is a function of the target photon field assumed to be the cosmic microwave background, as established in Petruk (2009). This method has been shown to be accurate from the Thomson to the Klein–Nishina regimes. We further assume that the maximum accelerated energy, E_max, is the same everywhere in the supernova remnant. The normalized emission maps are obtained by performing the integral,

$$\begin{eqnarray*} I = \int _{\rm shocks} j_{\rm IC} \mathrm{ d}l, \end{eqnarray*}$$

(14)

for the shocked material in the vicinity of the forward and reverse shocks of the supernova remnant.

These simple recipes permit the calculation of predictive non-thermal emission from supernova remnant simulations (e.g. Castellanos-Ramírez, Velázquez & Cantó 2021). A more careful determination of synchrotron and inverse-Compton fluxes of our supernovae remnants models would require the explicit knowledge of the accelerated particle distribution as done in Brose, Telezhinsky & Pohl (2016), Brose et al. (2019), Sushch et al. (2022), and Das et al. (2022).

2.5 Thermal emission maps

Last, maps of optical H α and soft X-rays emission from the core-collapse supernova remnants are calculated. The procedure is similar as for the non-thermal emission and uses different emission coefficients. The optical H α intensity reads

$$\begin{eqnarray*} I = \int _{\rm SNR} j_{\rm H\alpha } \mathrm{ d}l, \end{eqnarray*}$$

(15)

with

$$\begin{eqnarray*} j_{\rm H\alpha } = \Big (1.21 \times 10^{-22} T^{-0.9} \Big) n_{\rm H}^{2}, \end{eqnarray*}$$

(16)

where n_H is the gas hydrogen number density (Meyer et al. 2015). The soft X-ray emission is calculated as

$$\begin{eqnarray*} I = \int _{\rm SNR} j_{\rm XR} \mathrm{ d}l, \end{eqnarray*}$$

(17)

with

$$\begin{eqnarray*} j_{\rm XR} = \Lambda (T) n_{\rm H}^{2}, \end{eqnarray*}$$

(18)

with Λ(T) the |$0.1\!-\!1.0\, \rm keV$| emissivity for diffuse ISM obtained with the xspec software (Arnaud 1996).

3 RESULTS

In this section, we present our (non-)magnetized models for the pre-supernova circumstellar medium of static massive stars and for the evolution of their subsequent supernova remnant, together with predictions by radiative transfer calculation of their thermal and non-thermal appearance.

3.1 Asymmetric circumstellar medium

In Fig. 1, we show the density field (in |$\rm cm^{-3}$|⁠) in the simulations for the circumstellar medium of a static |$35\, \rm M_{\odot }$| star without (panel a, Run-35-HD-0-CSM) and with an ISM magnetic field (panel b, Run-35-MHD-0-CSM), see Table 1 for the models properties. The hydrodynamical simulation (Fig. 1a) exhibits the typical morphology of a stellar wind bubble generated by a massive static star. From the star to the outermost part of the bubble, we have: (1) the freely expanding Wolf–Rayet stellar wind, (2) the unstable shell of shocked Wolf–Rayet wind sweeping the hot diluted cavity of RSG material, (3) the contact discontinuity separating the main-sequence material and the ISM gas, (4) the large layer of cold, dense shocked ISM material, and (5) the forward shock sweeping the unperturbed ambient medium. The overall structure is globally spherically symmetric (Weaver et al. 1977), although the instabilities in the last Wolf–Rayet wind break this symmetry (Garcia-Segura et al. 1996a, b). We refer the reader interested in further details on the physics of stellar wind bubbles around massive stars to the literature devoted to this problem (Freyer et al. 2003, 2006; Dwarkadas 2005).

Figure 1.

Number density field (in |$\rm cm^{-3}$|⁠) in our models for the circumstellar medium of a static |$35\, \rm M_{\odot }$| star, shown at the supernova time, without (a) and with (b) ISM magnetic field. The red contour marks the region of the circumstellar medium where supernova material and energy is injected.

Fig. 1(b) displays the pre-supernova circumstellar medium of the massive star in a magnetized ISM. Only the strength of the magnetic field differs from Run-35-HD-0-CSM, in which it was set to |$B_{\rm ISM}=0\, \mu \rm G$|⁠. The evolution of the contact discontinuity during the main-sequence phase of the stellar evolution is asymmetric and elongated as a consequence of being inhibited along the direction normal to the vertical ISM field lines (van Marle et al. 2015). Note that the density at the forward shock of the wind bubble is smaller in the magnetized case, since the magnetic pressure reduces the compression ratios of shocks (Figs 1a and b). The distribution of the RSG stellar wind is strongly affected by the shape of the contact discontinuity and it adopts the same oblong configuration. Similarly, the expansion of the Wolf–Rayet stellar wind, initially spherical, is modified by the morphology of the distribution of the termination shock of the stellar wind bubble, see also Meyer et al. (2020b).

3.2 Asymmetric supernova remnant

In Fig. 2, we show as a time sequence the evolution of the supernova blastwave into the circumstellar medium of a |$35\, \rm M_{\odot }$| star in a magnetized ISM (model Run-35-MHD-0-CSM). The figures display the number density field of the supernova remnant between |$6.25\, \rm kyr$|⁠, (a), and |$35.0\, \rm kyr$|⁠, (i), after the explosion. The thin black contour marks the location where |$50{{\ \rm per\ cent}}$| of the number density is contributed by ejecta. The supernova shock wave first expands into the freely expanding Wolf–Rayet wind and approaches the unstable shell of shocked Wolf–Rayet gas (Fig. 2a). The shell has begun to interact with the contact discontinuity of the main-sequence wind bubble and loses its spherical shape (van Marle et al. 2015). As the blastwave continues expanding, the shock wave loses sphericity (Fig. 2b) and adopts an elongated morphology mirroring that of the cavity of the stellar wind bubble previously shaped by the progenitor star (Fig. 2c). This is precisely the mechanism that our study aims at highlighting, within the context of the core-collapse supernova remnant Puppis A.

Figure 2.

Number density field (in |$\rm cm^{-3}$|⁠) in our models for the supernova remnant of a static |$35\, \rm M_{\odot }$| star in a magnetized ISM. The black contour marks the region of the remnant made of |$50{{\ \rm per\ cent}}$| of ejecta in number density.

Partial reflection of the blastwave at the shell of Wolf–Rayet material sends gas back towards the centre of the explosion. Shock reflection is at work at the walls of the stellar wind cavity and, to a lesser degree, at the basis of the cylinder of Wolf–Rayet wind, which still expand and interact with the unperturbed wind material (Fig. 2d). A similar mechanism is at play at a stellar wind bow shock (Meyer et al. 2015) or a cold ISM region (Ferreira & de Jager 2008; Castellanos-Ramírez et al. 2021). The tubular morphology of the wind cavity consequently imposes the reflected shocks of supernova ejecta a cylindrical shape, which persists as it propagates towards the centre of the supernova remnant (Fig. 2e). Simultaneously, the forward shock of the supernova is transmitted through the Wolf–Rayet shell and keeps on expanding in the vertical direction of the cavity, as well as in the region of shocked ISM gas (Fig. 2f). After the vertically reflected waves join near the axis of symmetry of the supernova remnant, the process of reflection and transmission of the shock front continues and the supernova ejecta, Wolf–Rayet, red supergiant, and main-sequence materials mix within a thin unstable zone encompassing a dense central region of ejecta (Figs 2g–i).

3.3 Non-thermal emission

Fig. 3 shows the normalized radio synchrotron maps for our simulation Run-35-MHD-0-CSM at the time |$12.5\, \rm kyr$| (left-hand panels), |$16.5\, \rm kyr$| (left middle panels), |$17.5\, \rm kyr$| (right middle panels), and |$25.0\, \rm kyr$| (right-hand panels) after the supernova explosion, respectively. The images are obtained using the emission coefficient described in Section 2 and assuming viewing angles of θ_obs = 0° (top panels) and θ_obs = 45° (bottom panels).

Figure 3.

Normalized non-thermal radio synchrotron intensity maps of our supernova remnant model of a static |$35\, \rm M_{\odot }$| star, seen according viewing angles of θ_obs = 0° (top panels) and θ_obs = 45° (bottom panels). The images are displayed for time instances |$12.5\, \rm kyr$| (left-hand panels) to |$25.0\, \rm kyr$| (right-hand panels) after the supernova explosion. The images are centred on to the location of the supernova explosion.

When the shock wave hits the cavity walls, the radio appearance of the supernova remnant is that of two parallel bright bars aligned with the direction of the ISM magnetic field (Fig. 3a). This loss of sphericity of the expanding blastwave is also clearly visible with a viewing angle of θ_obs = 45° (Fig. 3e). As the blastwave propagates further and is reflected, we would see an outer hexagon produced by the ejecta penetrating the shock ISM material, inside of which two reflected parallel waves demarcate a bright bilateral region (Figs 3b and f). When the supernova blastwave reflects on to the shell of Wolf–Rayet stellar wind, it converges towards to the centre of the explosion. The region of shocked ejecta consequently adopts the shape of a bright rectangle embedded in an hexagon, at least for a viewing angle slightly tilted with respect to the plane of the sky (Figs 3c and g). However, once the reflected waves join in the centre of the supernova remnant, a very bright elongated region of shocked ejecta forms, and the remnant loses its rectangular shape (Figs 3d and h).

Fig. 4 compares synchrotron and inverse-Compton emission for the viewing angles θ_obs = 0° (top) and θ_obs = 45° (bottom) (left columns, respectively) at an age of 17.5 kyr, when the rectangular synchrotron morphology is most pronounced. The ring-like features which appear for a viewing angle θ_obs = 45° are artefacts produced by the 2D nature of the simulations. The synchrotron map with θ_obs = 0° traces the compressed magnetic field in the reflected ejecta and in the ejecta penetrating the shocked ISM bubble (Fig. 4a). The inverse-Compton emission is faint within the central rectangle. In contrast, the brightest inverse-Compton emission originates from the dense material at the unstable ejecta/wind interface, perpendicular to the direction of the ISM magnetic field (Fig. 4b). For both non-thermal emission mechanisms the viewing angle is an important factor for the rectangular morphology (Figs 4e and f).

Figure 4.

Normalized non-thermal radio synchrotron intensity (left) and inverse-Compton (middle left) as well as thermal optical H α (middle right) and thermal X-rays (right) emission maps of our supernova remnant model of a static |$35\, \rm M_{\odot }$| star in a magnetized medium at time |$17.5\, \rm kyr$| after the explosion, seen according to viewing angles of θ_obs = 0° (top) and θ_obs = 45° (bottom). The images are centre on to the location of the supernova explosion.

Fig. 5 displays a series of cross-sections taken vertically (left) and horizontally (right) through the intensity maps in Fig. 4 for a viewing angle of θ_obs = 0°. The top panels of the figures display the cuts for synchrotron (thick dashed red line) and inverse-Compton (thin solid blue line) emission. Fig. 6 is as Fig. 5, but assuming a viewing angle of θ_obs = 45°. The vertical slices show that the synchrotron cavity, extending from |$-20$| to |$20\, \rm pc$|⁠, is centre-filled by bright rings produced by the axisymmetric character of our MHD simulation. This is not seen in the inverse-Compton intensity, regardless of the viewing angle θ_obs (Figs 5a and 6a). The horizontal cross-sections highlight the differential expansion of the supernova blastwave into the tubular magnetized circumstellar medium which induce the rectangular morphology, i.e. the brightest peaks of the horizontal slice sit about |$20\, \rm pc$| from the centre of the explosion, while the vertical extent reaches up to |$40\, \rm pc$|⁠.

Figure 5.

Cuts through the normalized synthetic non-thermal (a,b) and thermal (c,d) emission maps of our rectangular supernova remnant at time |$17.5\, \rm kyr$|⁠, taken through the centre of the explosion, along the vertical direction (left column) and the horizontal direction (right column), under a viewing angle of θ_obs = 0°.

Figure 6.

As Fig. 5, but with θ_obs = 45°.

3.4 Thermal emission

The right-hand series of panels in Fig. 4 display optical H α (panels c, g) and thermal X-ray (panels d, h) emission maps at the time of a prominent rectangular synchrotron morphology. The optical image traces the reverse shock of the blastwave being reverberated towards the centre of the explosion, giving the remnant its peculiar morphology, with a central dense ring that reflects the initial inner ejecta profile. The transmitted shock passing through the magnetized and elongated walls of the stellar wind cavity is not visible despite the high post-shock density, since the high temperature in this region forbids generous optical emission. The soft thermal X-ray emission basically originates from the supernova blastwave which penetrates into the walls of the cavity, where the gas is both dense and hot. The morphology of the projected emission is square, as is the overall shape of Puppis A.

The bottom panels of Figs 5 and 6 plot cross-sections taken vertically (panel c) and horizontally (panel d) through the optical (thick dashed red line) and thermal X-ray (thin solid blue line) emission maps, respectively. Both plots additionally highlight that the H α photons are mostly emitted in the region limited by the reverberated forward shock of the supernova blastwave. Its vertical and horizontal extension away from the centre of the explosion differs greatly so that in projection the region looks like a rectangle. Conversely, the thermal X-ray emission finds its origin in the supernova shock wave penetrating into the shocked material of the stellar wind bubble. Both vertical and horizontal distributions of the emitted material extend to about |$40\, \rm pc$| from the centre of the explosion, and in thermal X-rays the supernova remnant looks like a square surrounding the inner optical rectangle.

4 DISCUSSION

The caveats and limitations of our study are presented in this section and the results are discussed in the context of the rectangular supernova remnant Puppis A.

4.1 Limitation of the model

Two main caveats affect our simulation method. First of all, regarding our choice in terms of stellar wind boundaries, we neglect the magnetization of the progenitor stellar wind, as in the precedent papers of this series (Meyer et al. 2020a, 2021a), and as discussed in great details in van Marle et al. (2015). Nevertheless, this is a pre-requisite in the careful estimation of cosmic-rays propagation through core-collapse supernova remnants (Das et al. 2022; Sushch et al. 2022). The adopted scenario for the evolution of the progenitor is arbitrarily taken to be that of a non-rotating |$35\, \rm M_{\odot }$| star with solar metallicity (Ekström et al. 2012). Hence, this first qualitative study for the formation of rectangular middle-age core-collapse supernova remnants is not fine-tuned to a particular object (see Section 4.2).

Secondly, we neglect any turbulence in the medium into which the star blows its wind. In reality, a surrounding H ii region, potentially trapped into the dense layer of the growing stellar wind bubble (Weaver et al. 1977; van Marle, Langer & García-Segura 2004), develops velocity fluctuations when expanding into an inhomogeneous molecular cloud (Medina et al. 2014). Furthermore, the dense environment hosts star formation processes and typically is highly magnetized (Hennebelle & Falgarone 2012), influencing the development of circumstellar structures (Pardi et al. 2017) and supernova remnants (Korpi et al. 1999a, b), which in turn generate and maintain the level of turbulence of the ISM (Seifried et al. 2018). A more realistic depiction of the supernova–wind interaction in a magnetized ISM would require full 3D MHD simulations, better modelling of the circumstellar medium around the massive progenitor star (Meyer et al. 2021b), as well as the propagation of the supernova blastwave in it (Orlando et al. 2021, 2022).

4.2 Interpretation and comparison with observations

4.2.1 A new model for rectangular core-collapse remnants

Our study reveals how core-collapse supernova remnants can reveal a rectangular appearance, by a succession of shock reflections on the walls of the stellar wind cavity, provided their progenitor star blows its pre-supernova material in a medium that is uniformly magnetized (van Marle et al. 2015). This mechanism is the adapted version of the interacting-winds model employed for explaining the formation of bipolar and elliptical planetary nebulae, in a sense that the final outcome is attributed to the interaction of a rather spherical outflow with a pre-existing, non-spherical circumstellar structure.

In that mass regime, other models explain the formation of bipolar and multipolar planetary nebulae, characterized by showing point-symmetric features. These models consider a binary system inside (Bond, Liller & Mannery 1978; Livio, Salzman & Shaviv 1979; Soker & Livio 1994). One of the binary components launches a pair of precessing jets (Soker & Rappaport 2000) which interact with the AGB wind of the other star. This interaction produces point-symmetric planetary nebulae as observed in the Red Rectangle (Cohen et al. 2004; Velázquez et al. 2011), Hen 3-1475 (Riera et al. 2014), and IC 4634 nebulae (Guerrero et al. 2008). The interacting wind model considers the collision of an isotropic, fast, and low-density stellar wind (coming from a low-mass star) with a slow and high-density AGB wind, with a high-density equatorial region (shaped by stellar rotation). The wind–wind interaction induces an asymmetric flow and opens the stellar magnetic field lines (García-Segura & López 2000). Analogically, in our model involving a massive progenitor star, the spherical outflow is represented by the supernova ejecta, while the non-spherical circumstellar medium is sculptured by both the red supergiant and Wolf–Rayet stellar winds under the action of ISM magnetic fields.

We interpret the formation of such rectangular core-collapse supernova remnants as indication of slow or absent motion of the progenitor star. If the massive star were a fast-moving object, the cavity of stellar wind would be strongly asymmetric, and so would be the propagation of the shock wave (Meyer et al. 2015).

4.2.2 Comparison with Puppis A

Puppis A is the archetype of core-collapse supernova remnant with a rectangular morphology. It displays an unusual angular shape shown in a composite X-ray rendering (see Fig. 7) and a rectangular morphology in radio continuum (see fig. 1 in Reynoso et al. 2017), with an intensity enhancement towards the Eastern region. The massive nature of Puppis A’s progenitor has been constrained due to (i) its O and Fe enriched interior (Charles, Culhane & Zarnecki 1978; Katsuda et al. 2013) and (ii) the runaway pulsar it hosts (Holland-Ashford et al. 2017; Vogt et al. 2018; Mayer et al. 2020). This supernova remnant has mainly been observed in the radio (Reynoso et al. 2017), infrared (Arendt et al. 1991, 2010; Arendt, Dwek & Petre 1991), optical (Goudis & Meaburn 1978), and X-Ray (Dwek et al. 1987; Dubner et al. 2013) band, providing a multiwavelength picture of the rectangular shape. Infrared emission from scattered light on dust trapped into the supernova remnant reveals that the supernova blastwave is not expanding into a stratified medium, except in the eastern region (Arendt et al. 1990), which is in accordance with our model.

Figure 7.

X-ray view of the supernova remnant Puppis A as a mosaic of XMM–Newton and Chandra data, taken from Dubner et al. (2013), where the red (⁠|$0.3\!-\!0.7\, \rm keV$|⁠), green (⁠|$0.7\!-\!1.0\, \rm keV$|⁠), and blue (⁠|$1.0\!-\!8.0\, \rm keV$|⁠) colours indicate each emission energy band. The supernova remnant exhibits a morphology strongly deviating from sphericity, appearing as an overall rectangular shape inside of which multiple square structures are imbricated.

High-energy data provided by the Fermi-LAT and H.E.S.S. facilities indicate that Puppis A is the site of particle acceleration up to the GeV band only, which supports the scenario of the supernova blastwave not yet massively interacting with the surrounding dense molecular clouds (Hewitt et al. 2012; H. E. S. S. Collaboration 2015; Xin et al. 2017). This fact favours the rectangular morphology being principally produced by the blastwave colliding with circumstellar material, as proposed our model. In particular, a blastwave passing through hot, shocked circumstellar material tends to produce soft cosmic ray spectra (Das et al. 2022), leading to much weaker emission in the TeV band than at GeV energies, as is observed. The sole signs of isolated cloud collision in the bright eastern knot of Puppis A (Paron et al. 2008), not included into our simulation, are therefore not responsible for its overall shaping.

Our numerical model is consistent with the observations of Puppis A, as it displays an overall squared X-ray shape (Figs 4d and h) with smaller imbricated optical rectangles (Figs 4c and g). This finding permits us to propose that the progenitor of Puppis A was a massive star which evolved in a background medium of organized magnetization in agreement with the scenario presented by Reynoso, Velázquez & Cichowolski (2018).

Throughout stellar evolution, the many stellar evolutionary phases released first main-sequence material that first shape an elongated cavity in which an evolved red supergiant wind is blown (Reynoso et al. 2017, 2018), and, depending on the mass of the progenitor, it might have further evolved up to the Wolf–Rayet phase. An extra numerical effort is necessary to fine-tune simulations to the case of Puppis A and its shaping.

5 CONCLUSION

This study investigates the possibility that core-collapse supernova remnants adopt a rectangular emission morphology. We perform 2.5D MHD simulations of the circumstellar medium of a |$35\, \rm M_{\odot }$| single massive star at rest (Ekström et al. 2012) evolving and dying in the warm phase of the Galactic plane of the Milky Way. The numerical models are performed with the pluto code (Mignone et al. 2007, 2012; Vaidya et al. 2018), which has been previously used to study the surroundings of massive stars and their associated supernova remnants (Meyer et al. 2014).

It is known that magnetic field of order of few |$\mu$|G is not effective to modify the shape of supernova remnant because the ratio of the magnetic pressure to the thermal pressure is small and magnetic field is not dynamically important in evolution of the remnant. However, it is very effective to determine the structure of the ambient density where this ratio is considerably larger than unity (van Marle et al. 2015). We demonstrate that such circumstellar medium of the progenitor star, profoundly pressure-supported by an organized background magnetic field, subsequently constrains the propagation of the supernova shock wave, inducing multiple shock reflections along the directions normal to the ISM magnetic fields. Eventually, in the middle age of its evolution (⁠|$15\!-\!20\, \rm kyr$|⁠), the remnant adopts the unusual rectangular morphology due to the isotropically expanding supernova ejecta interacting with the tubular structure of the progenitor surroundings.

We calculated intensity maps of radio synchrotron and inverse-Compton emission using radiative transfer calculations. At some viewing angles, these maps show a rectangular morphology for the age window between the reflection of the supernova shock wave off the circumstellar medium and the moment it reaches the centre of the remnant. This mechanism is the adapted version of the shaping of bipolar and elliptical planetary nebula by the interaction of two stellar winds (see García-Segura & López 2000). Our scenario for rectangular core-collapse supernova remnants requires that its progenitor was not a runaway star, for which the remnant will shape following the scenario of Meyer et al. (2021a). We suggest that the asymmetric supernova remnant Puppis A, depicting a rectangular morphology, is at least partly shaped by the mechanism described here.

Further investigations are necessary to validate our model for rectangular core-collapse supernova remnants, both on the observational and on the numerical side. The former to update our knowledge on the abundance of rectangular remnants from massive stars, the latter to understand the detailed morphology of Puppis A (Reynoso et al. 2017).

ACKNOWLEDGEMENTS

The authors thank the anonymous referee for comments which improved the quality of the paper. DMAM thanks L. Oskinova for discussion on X-ray observational data. The authors acknowledge the North-German Supercomputing Alliance (HLRN) for providing HPC resources that have contributed to the research results reported in this paper. MP acknowledges the Max Planck Computing and Data Facility (MPCDF) for providing data storage resources and HPC resources which contributed to test and optimize the pluto code. PFV, JCT-R, and AC-R acknowledge the financial support for PAPIIT-UNAM grants IA103121 and IG100422. AC-R acknowledges support from CONACYT postdoctoral fellowship. AC-F acknowledges financial support by the Spanish Ministry of Science and Innovation through the research grant PID2019-107427GB-C31. EMR is a member of the Carrera del Investigador Científico of CONICET, Argentina, and is partially funded by CONICET grant PIP 112-201701-00604CO.

DATA AVAILABILITY

This research made use of the pluto code developed at the University of Torino by A. Mignone (http://plutocode.ph.unito.it/) and of the radmc-3d code developed at the University of Heidelberg by C. Dullemond (https://www.ita.uni-heidelberg.de/~dullemond/software/radmc-3d/). The figures have been produced using the matplotlib plotting library for the Python programming language (https://matplotlib.org/). The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

1

http://plutocode.ph.unito.it/

2

https://www.ita.uni-heidelberg.de/~dullemond/software/radmc-3d/