© ESO 2021

1 Introduction

Luminous blue variable stars (LBVs) are massive evolved unstable stars that strongly interact with the circumstellar medium (CSM), showing dramatic variations in both their spectra and brightness (Humphreys & Davidson 1994; Humphreys et al. 1999). This class includes stars displaying variability on different timescales and of differing intensity (de Koter et al. 1996). The most typical LBVs are the S Doradus variables, which are characterized by ≈0.5–2 mag quasi-periodic variations on a timescale of years to decades. Giant eruptions instead are much less common than the variability observed in S Doradus variables and determine greater brightness variations (≥1–2 mag) associated with an episode of high mass loss (~10 M_⊙). Examples of these giant eruptions are those occurring in the case of η-Carinae (see Davidson & Humphreys 1997; Smith et al. 2018). Although we still do not fully understand the physical mechanism that drives LBV variability, some progress has been made in recent years. Recently, Grassitelli et al. (2021) developed a model that reproduces the typical observational phenomenology of the S Doradus variability. According to their model, the instability responsible for the observed variability can be triggered when some physical conditions are met, involving inflated envelopes in proximity to the Eddington limit, a temperature range that does not lead to accelerating outflows, and a mass-loss rate that increases with decreasing temperature (see Grassitelli et al. 2021, for more details). Nevertheless, the causes of the sporadic and violent mass-loss events are still poorly understood and a physical driving mechanism has not yet been clearly identified (Smith et al. 2018).

According to the traditional stellar evolution theory, LBVs are presumed to represent a brief transitional phase in the evolution of the most massive stars, between the main sequence O-type stars and the Wolf–Rayet (WR) stars, but in the last decade our theoretical understanding of massive stellar evolution has been radically revised (see Smith 2017b). In particular, in a number of recent works, LBVs have been proposed as the progenitors of some core-collapse supernovae (SNe), particularly Type IIn (e.g., Gal-Yam et al. 2007; Gal-Yam & Leonard 2009; Miller et al. 2010; Smith et al. 2011; Fox et al. 2011; Taddia et al. 2013, 2020; Pastorello et al. 2018) and Type IIb SNe (e.g., Groh et al. 2013; Moriya et al. 2013; Prentice et al. 2020). In many cases, the link with SNe stems from the fact that, in this class of SNe, the blast wave appears to expand in a high-densitymedium, which has been interpreted as resulting from a wind with velocities and mass-loss rates consistent with those expected for massive LBVs. However, discriminating between some stellar transients and the SN events can be a tricky issue in some cases, considering the amount of energy released during some stellar outbursts (see Pastorello & Fraser 2019). Indeed, during giant eruptive episodes, LBVs can reach total luminosities comparable to that of a SN, mimicking in some cases the behavior of real SNe IIn; these events (that have nothing to do with SNe) are known as SN impostors (Van Dyk et al. 2000; Smartt 2009; Van Dyk & Matheson 2012; Tartaglia et al. 2015).

Despite the similarities between the characteristics of some SNe (namely those showing strong shocks resulting from the interaction between their ejecta and the pre-existing dense CSM) and those of LBVs (more specifically, the giant eruption of mass and the highly structured and dense CSM), firm evidence that LBVs are direct core-collapse SNe progenitors is still lacking (Dwarkadas 2011; Smith 2017b; Van Dyk 2017; Nyholm et al. 2020). Only Type IIn SN 2005gl was found to be a valid candidate for an LBV-like progenitor (Gal-Yam et al. 2007; Gal-Yam & Leonard 2009). In fact, a large population of sources closely related to the LBV stage of massive stellar evolution has been identified (Wachter et al. 2010; Gvaramadze et al. 2010, 2015; Smith et al. 2019), but their connection with some Type IIn/IIb SNe is not clear. On the other hand, the information about the environments in which the stars explode may establish some constraints on the evolutionaryphases of the progenitors (Fox et al. 2011; Anderson et al. 2012; Nyholm et al. 2020). However, while observationsshow that some SNe seem to have had LBVs as progenitors, there are only a couple of theoretical models that support LBVs exploding as SNe in the framework of single stellar evolution (Groh et al. 2013; Moriya et al. 2013).

Another interesting clue to be explored is the possible identification of signatures of LBV progenitors in the ejecta distribution and morphology of supernova remnants (SNRs). In fact, the morphology and the distribution of material observed in SNRs is expected to reflect the interaction of the SN blast wave with the ambient environment (e.g., Ustamujic et al. 2021), the physical processes associated with the SN explosion, and the nature of the progenitor star. Being able to disentangle the different effects is a challenging puzzle to explore and solve (e.g., Orlando et al. 2015, 2016, 2019, 2020, 2021; Wongwathanarat et al. 2017; Ferrand et al. 2019, 2021; Tutone et al. 2020; Gabler et al. 2021; Jacovich et al. 2021). Recently, Chiotellis et al. (2021) proposed that the two protrusions that are projected as two “ears” in the morphology of some SNRs could be formed by the interaction of the remnant with the CSM, considering a LBV or a red/yellow supergiant as the most likely progenitors in the case of a core collapse SN event. Other authors instead support the scenario of the jet-driven core-collapse SN mechanisms to explain these elongated features (see Grichener & Soker 2017; Bear & Soker 2018).

In this work, we investigate the physical, chemical, and morphological properties of SNRs of stars exploded as SNe at (or soon after) the LBV evolutionary stage. Generally, LBVs are surrounded by extended circumstellar envelopes that show a wide variety of characteristics reflecting the mass-loss history of the variable stars. For our study, we adopted the LBV candidate Gal 026.47+0.02 as a template (hereafter G26; see Clark et al. 2003, 2005), which exhibits one of the highest observed mass losses from the central object. G26 is located in a very massive nebula that has been extensively studied. Indeed, an accurate description of the structure and density distribution of the nebula has been derived from the analysis of observations (see Paron et al. 2012; Umana et al. 2012). The extreme characteristics of the environment associated with G26 (see Umana et al. 2012) fit quite well with the requirements for the progenitor of Type IIn SN 2010jl according to the model by Andrews et al. (2011). This may indicate that G26 could be the precursor of a very bright Type IIn SN. The central star shows a luminosity log (L∕L_⊙) = 6, a temperature T = 17 000 K, and mass ejection rate Ṁ = 9 × 10⁻⁵ M_⊙ yr⁻¹ (Clark et al. 2003; Wachter et al. 2010). Assuming a distance of 6.5 kpc (Clark et al. 2003), the nebula G26 consists of ≈ 20 M_⊙ of ionized gas distributed in two nested tori (or shells) around a common axis which have been interpreted as the observational evidence of past episodic mass-loss events (Umana et al. 2012). Its temperature and luminosity place G26 in a region of the HR diagram that is populated by relatively typical LBV stars (see Figs. 1 and 3 in Smith 2017b). The high total mass of material detected in G26 could be an indication of a few giant energetic eruptions that occurred in the past or a number of moderate mass-loss events (Umana et al. 2012). Thus, the distribution of material observed in the circumstellar environment reflects the mass-loss history of the star and is therefore particular to G26. According to evolutionary models, G26 evolved from a star with an initial mass of between 60 and 80 M_⊙ (see Fig. 14 in Limongi & Chieffi 2018; see also Smith 2017b).

For our purposes, we selected a grid of massive stars, in agreement with the characteristics of G26, which explode as core-collapse SNe from the stellar models described in Limongi & Chieffi (2018). We then performed three-dimensional (3D) hydrodynamic (HD) simulations that follow the post-explosion evolution of the SN from the breakout of the shock wave at the surface of the LBV progenitor to the interaction of the SNR with the circumstellar environment. We explored a grid of eight models that differ in terms of the progenitor star characteristics and the explosion energy, and considering the same circumstellar environmentas that described in Umana et al. (2012).

The paper is organized as follows. In Sect. 2, we describe the model and the numerical setup; in Sect. 3 we discuss the results; and in Sect. 4 we draw our conclusions.

2 Hydrodynamic model

2.1 Numerical setup

The model describes the post-explosion evolution of a core-collapse SN from the breakout of the shock wave at the stellar surface (occurring a few minutes after the SN event) to the interaction of the blast wave and ejecta caused by the explosion with the circumstellar environment. We followed the evolution for t ≈ 10 000 yr by numerically solving the full time-dependent HD equations in a 3D Cartesian coordinate system (x, y, z). The HD equations were solved in the conservative form, (1) (2) (3)

where E = ϵ + u²∕2 is the total gas energy (internal energy ϵ, and kinetic energy) per unit mass, t is the time, ρ = μm_Hn is the mass density, μ is the mean atomic mass (μ for the ejecta considers their isotopic composition, whereas μ = 1.29 for the CSM, assuming cosmic abundances), m_H is the mass of the hydrogen atom, n is the total number density, u is the gas velocity, and T is the temperature. We used the ideal gas law, P = (γ − 1)ρϵ, where γ = 5∕3 is the adiabatic index.

The calculations were performed using PLUTO (Mignone et al. 2007), a modular Godunov-type code for astrophysical plasmas. The code provides a multiphysics, multialgorithm modular environment particularly oriented towards the treatment of astrophysical, high-Mach-number flows in multiple spatial dimensions. The code was designed to make efficient use of massive parallel computers using the message-passing interface (MPI) library for interprocessor communications. The HD equations are solved using the HD module available in PLUTO; the integration is performed using the original Piecewise Parabolic Method (PPM) reconstruction by Colella & Woodward (1984, see also Miller & Colella 2002) with a Roe Riemann solver. The adopted scheme is particularly appropriate for describing the shocks formed during the interaction of the remnant with the surrounding inhomogeneous medium, as in our case. A monotonous central difference limiter (the least diffusive limiter available in PLUTO) for the primitive variables is used. The code was extended by additional computational modules to evaluate the deviations from equilibrium of ionization of the most abundant ions (through the computation of the maximum ionization age in each cell of the spatial domain as described in Orlando et al. 2015), and the deviations from temperature equilibration between electrons and ions. For the latter, we included the almost instantaneous heating of electrons at shock fronts up to kT ~ 0.3 keV by lower hybrid waves (see Ghavamian et al. 2007), and the effects of Coulomb collisions for the calculation of ion and electron temperatures in the post-shock plasma (see Orlando et al. 2015, for further details).

2.2 Initial and boundary conditions

We modeled the post-explosion evolution of a core-collapse SN starting immediately after the shock breakout, and followed the transition from the SN to the SNR phase and the interaction of the remnant with the inhomogeneous pre-SN environment. As initial conditions, we adopted the explosive nucleosynthesis models described in Limongi & Chieffi (2018). These authors presented a grid of pre-SN models of massive stars whose mass spans the range between 13 and 120 M_⊙, covering four metallicities (i.e., [Fe/H] = 0, −1, −2, and −3) and three initial rotation velocities (i.e., 0, 150, and 300 km s⁻¹). In particular, we selected the models with solar metallicity and star masses of either 60 or 80 M_⊙, which reproduce evolutionary tracks in agreement with the position of G26 in the Hertzsprung–Russell (HR) diagram (see left upper panel in Fig. 1; see also Fig. 14 in Limongi & Chieffi 2018, and Fig. 3 in Smith 2017b). For the initial rotation velocity, we explored the models with the two extreme values¹, namely either 0 or 300 km s⁻¹. For completeness, for each case, we selected a SN with either low or high explosion energy, considering the total kinetic energy of the ejecta². A summaryof the cases explored is given in Table 1, where we report³: the main sequence mass of the star, M_*; the initial rotation velocity of the star, V_rot; and the energy of the explosion, E_exp.

The initial blast wave is defined from the one-dimensional (1D) profiles of density, pressure, velocity, and abundances of ten species⁴ (⁴He, ¹²C, ¹⁴N, ¹⁶O, ²⁰Ne, ²⁴Mg, ²⁸Si, ⁴⁰Ca, ⁴⁴Ti, ⁵⁶Ni) describing the ejecta after the shock breakout. In Fig. 1, we present the profiles of density for the models outlined in Table 1 (top right panel), and the abundances of the ten species considered for two reference cases, one with low (left bottom panel) and the other with high (right bottom panel) explosion energy, and both for a zero-age main sequence (ZAMS) star of 60 M_⊙. The range of enclosed mass plotted on the x-axis is very different for low- and high-explosion energy models because of the different locations of the mass cut (see Limongi & Chieffi 2018).

In our models, these 1D profiles are mapped in the 3D domain, assuming spherical symmetry, and centered at the origin of the 3D Cartesian coordinate system. We assumed a clumpy initial density structure of the ejecta, as suggested by theoretical and spectropolarimetric studies (e.g., Nagataki 2000; Kifonidis et al. 2006; Wang et al. 2003, 2004; Wang & Wheeler 2008; Gawryszczak et al. 2010; Hole et al. 2010; Wongwathanarat et al. 2015). Thus, after the 1D profiles of ejecta are remapped to the 3D domain, the small-scale structure of the ejecta is modeled as per-cell random density distributions by adopting a power-law probability distribution (see Orlando et al. 2012). In our simulations, the ejecta clumps havethe same initial size (about 2% of the initial remnant radius), and a maximum density perturbation ν_max = 10. In Table 2, we outline the total masses of the chemical elements composing the ejecta for all the models presented in Table 1. In Table 3, we report the total mass of the fallback, M_rem; the mass of the ejecta, M_ej; the initial time of the simulation, t₀; the radius of the sphere containing the ejecta at t₀, R_ej; the extension of the domain, in the first, D_i, and last, D_f, remapping (see later).

The pre-SN environment is parametrized by following the two-shell model proposed for G26 by Umana et al. (2012) based on radio and infrared observations (see Fig. 2; see also Fig. 2 in Umana et al. 2012). According to Umana et al. (2012), the very massive nebula of G26 consists of at least 17 M_⊙ of ionized gas⁵ divided into 7.6 M_⊙ and 9.7 M_⊙ in the innerand outermost nebula, respectively. The estimated total mass of dust is 1.2–3.2 × 10⁻² M_⊙, which is several orders of magnitude lower than the mass of the ionized gas and thus negligible for the purposes of this paper. In light of this,we defined an ambient medium consisting of a spherically symmetric steady wind and two dense nested shells (representing the inner and outermost nebula observed) with a common axis coincident with the z-axis (see bottom panel in Fig. 2). We assumed the two shells (possibly related to different mass-loss episodes occurred in the past) each with a mass of ~ 10 M_⊙, namely slightly higher than the two lower limits (7.6 M_⊙ and 9.7 M_⊙) found by Umana et al. 2012). The spherically symmetric wind is characterized by a gas density proportional to r⁻², defined following a mass-loss rate⁶ Ṁ = 10⁻⁴ M_⊙ yr⁻¹. We fixed a lower threshold of 0.1 cm⁻³ for the pre-SN density of the CSM, which corresponds to assuming a progressive flattening of the wind profile to a uniform density at large radii⁷. We note that the environment outlined in the bottom panel in Fig. 2 is a simplified version of the CSM expected in a LBV exploding star. In particular, we expect to find a more complex and dense CSM close to the star. However, due to the lack of information about the more internal CSM, here we focus on investigating the possible effects of an SN interacting with the toroidal structures identified in G26 by Umana et al. (2012).

The two shells are defined as clumpy structures azimuthally symmetric about the z-axis, centered at the origin of the 3D Cartesian coordinate system (x₀, y₀, z₀) = (0, 0, 0). They follow the shape of an elliptic torus, that is, a surface of revolution produced by rotating an ellipse, in this case, about the z-axis. The dimensions of the shells were chosen following the results of Umana et al. (2012). The radius of the inner and outer shells (i.e., the distances from the center of explosion to the center of the ellipses) are ≈ 0.6 pc and ≈ 1.6 pc, respectively. Assuming a wind velocity of 200 km s⁻¹ (e.g., Smith 2017b), the shells would be the relics of mass-ejection episodes occurring between 3000 and 8000 years before the SN event. The length of the major semi-axis (along the z-direction) is 0.8 and 1.5 pc for the inner and outer shells, respectively; the size of the minor semi-axis (in the xy-plane) is 0.15 pc in both cases.

The two dense shells are expected to play a central role in modifying the expansion of the forward shock and in driving a reflected shock through the ejecta. As the geometry and density distribution adopted in the paper are idealized, this may introduce some features in the remnant structure if the shells are assumed to be uniform. In order to get a non-uniform density distribution for the two shells, the material was modeled as a set of randomly distributed spherical clumps of ≈ 0.15 pc in radius mixed with spherical subclumps of ≈0.05 pc in radius (see, for instance, Ustamujic et al. 2021). The density of the plasma in the spheres follows a normal distribution with a mean density of 160 cm⁻³ and 45 cm⁻³ for the innerand outer shells, respectively, in agreement with the values of average electron density estimated by Umana et al. (2012) for the corresponding shells in G26 (assuming a distance of 6.5 kpc, see Clark et al. 2003). The total mass in each modeled shell is ≈ 12 M_⊙, in agreement with the estimates given by Umana et al. (2012). We assume the two dense shells to be in pressure equilibrium with the environment. We note that we neglected any bulk velocity in both components of the pre-SN CSM, namely the wind component and the two dense shells, because it would be much smaller than the velocity of the forward shock anyway. Furthermore, because of a lack of observational constraints and evidence, we neglected mass eruptions that may have occurred in the latest phases of the progenitor evolution before core collapse and that could have generated a dense and inhomogeneous medium in the immediate surrounding of the SN. We have no indications of this material from observations of G26, and we prefer to keep the CSM as simple as possible to investigate the effects of the extended shells observed in G26 on the SNR evolution.

The simulations include passive tracers (C) to follow the evolution of the different plasma components (the ejecta and the two dense shells), to store information on the shocked plasma (time, shock velocity, and shock position when a cell of the mesh is shocked by either the forward or the reverse shock), and to follow the chemical evolution of the ejecta for ten different species (⁴He, ¹²C, ¹⁴N, ¹⁶O, ²⁰Ne, ²⁴Mg, ²⁸Si, ⁴⁰Ca, ⁴⁴Ti, ⁵⁶Ni). The continuity equations of the tracers are solved in addition to our set of HD equations. In the case of tracers associated with the different plasma components (ejecta or the shells), each material is initialized with C_i = 1, while C_i = 0 elsewhere, where the index i refers to the ejecta as a whole, or to the material in the inner or in outer shell. The chemical evolution of the ejecta is followed by adopting a multiple fluids approach (see e.g., Orlando et al. 2016, 2021). The fluids following the evolution of the different species are initialized with the abundances calculated at the shock breakout with the SN models described in Limongi & Chieffi (2018) and calculated from the set of stellar models summarized in Table 1 (see Fig. 1). The different fluids mix together during the evolution and in particular when the ejecta interact with the reverse shock that develops during the expansion of the remnant. The density of a specific element in a fluid cell is calculated as ρ_i = ρ ⋅ C_i, where, in this case, C_i is the mass fraction of each element and the index i refers to theconsidered element. This approach allows us to follow the spatial distribution of the chemical elements both inside and outside the reverse shock during the model evolution. All the other tracers (in particular those that store information on the shocked plasma) are initialized to zero everywhere.

The computational domain is a Cartesian box covered by a uniform grid of 512 × 512 × 512 zones, including the initial remnant defined by the profiles presented in Fig. 1(see Limongi & Chieffi 2018) mapped in 3D. The initial computational domain extends from (−D_i∕2) to (D_i ∕2) in all directions (see Table 3), to cover the sphere containing the ejecta at the beginning of the simulation, leading to a spatial resolution in the range of 4.5 × 10⁹–9.2 × 10¹¹ cm depending on the model. In order to follow the large physical scales spanned during the remnant expansion, we adopted the approach described in Orlando et al. (2019, 2020). During the evolution, the computational domain was gradually extended following the expansion of the remnant through the CSM and the physical quantities were remapped in the new domain. The domain is extended by a factor of 1.2 in all directions when the forward shock reaches one of the boundaries of the Cartesian box. The numberof mesh points is the same at each remapping, and therefore the spatial resolution gradually decreases during evolution. In each remapping, all the physical quantities in the new region – added outside the previous computational domain – are set to the values of the pre-SN CSM. For the models explored here, between 74 and 100 remappings were necessary to follow the interaction of the remnant with the CSM during 10000 yr of evolution. The final domain extends between (−D_f∕2) and (D_f ∕2) in all directions (see Table 3), leading to a spatial resolution in the range 3.0 × 10¹⁷ − 5.5 × 10¹⁷ cm depending on the model. All physical quantities were fixed to the values of the pre-SN CSM at all boundaries.

Fig. 1 Left upper panel: Pre-SN evolution of stars of 60 M⊙ and 80 M⊙ (from the main sequence phase up to the core He depletion stage, see Limongi & Chieffi 2018), with solar (initial) metallicities and initial rotation velocities Vrot = 0 km s−1 (solid lines) and Vrot = 300 km s−1 (dotted lines). The yellow asterisk marks the position of G26 in the HR diagram. Right upper panel: initial radial profiles of density in the models defined in Table 1. Models with initial rotation velocities Vrot = 0 km s−1 and Vrot= 300 km s−1 are represented with solid and dotted lines respectively. Lower panels: initial radial profiles of mass fractions for the species considered for two representative models, one with a low explosion energy (left panel) and the other with a high explosion energy (right panel). In both cases, the model describes a ZAMS star of 60 M⊙. The range of enclosed mass (x-axis) is very different in these models because of the different location of the mass cut (see Limongi & Chieffi 2018).

Fig. 2 Top: map of G26 which combines infrared at 24 μm (blue) and radio (red) images (see Umana et al. 2012). Bottom: schematic view of the CSM around G26, representing the two dense shells (in red and light blue) observed in the radio and infrared bands. The LBV is located in the center of the two nested toroidal shells.

3 Results

3.1 Hydrodynamic evolution

Soon after the shock breakout, the ejecta propagate freely through the spherically symmetric wind, driving a forward shock in the wind and a reverse shock backward through the ejecta. During this phase of evolution, the unshocked ejecta expand almost homologously, thus maintaining their initial structure and chemical stratification. After ≈ 18–40 yr of evolution, depending on the case, the forward shock hits the innermost dense shell of the CSM (see Fig. 2). In this section, we describe the interaction of the modeled remnants (see Table 1) with the two shells and their subsequent evolution for ≈10000 years. In all the figures presented in this section, we assumed the system oriented as the toroidal shells of G26, as deduced from the analysis of observations (see Fig. 2 in Umana et al. 2012). To this end, we rotated the original system about the three axes by the angles i_x = 30°, i_y = 30° and i_z = 25° to fit the orientation of the shells in G26 with respect to the line of sight (LoS).

The evolution of model M60-V0-1Foe is shown in Figs. 3 and 4 at different epochs (increasing from upper left to lower right panels). Figure 3 reports the density distributions for the ejecta in logarithmic scale; Fig. 4 shows the same isosurfaces that appear in Fig. 3, but colored according to the corresponding radial velocity (we note the different color scale in the lower panels). The complete temporal evolution is available as online movies (Movie 1 and Movie 2). The total mass of the ejecta in this model is low, namely ≈ 4 M_⊙ (see Table 3), as most of the stellar mass was lost in the CSM during the pre-SN phases of the star, and in the fallback during the core-collapse event (≈14 M_⊙, see Table 3). The forward shock, whose position is indicated with a semi-transparent surface in both figures, starts to interact with the innermost dense shell (red clumpy semi-transparent structure in the top left panel in Fig. 3) after ~ 40 yr of evolution.The interaction determines a strong slowdown of the forward shock which propagates through the shell (see top left panel in Fig. 4), and, consequently, a strengthening of the reverse shock traveling through the ejecta. Meanwhile, the forward shock continues to expand freely through the wind along the polar directions (roughly along the z-axis). As a result, the initial quasi-spherical distribution of ejecta progressively becomes asymmetric (elongated in the z direction) inthe subsequent evolution (see top central panel in Figs. 3 and 4). In a few years, at t ≈ 115 yr, the forward shock goes beyond the inner shell in the equatorial plane and starts traveling again through the wind. After ≈ 250 yr of evolution, the blast hits the outer dense shell (cyan semi-transparent clumpy structure in the top middle panel in Fig. 3). Similarly to the previously encountered shell, the dense toroidal structure slows down the expansion of the forward shock (see top right panel in Fig. 4), and, again, it produces a strengthening of the reverse shock traveling through the ejecta. This together with the free expansion of the ejecta at the poles further enhances the elongated shape of the remnant (see top right panel in Fig. 3). The reverse shock, powered by the interaction with the two dense shells, heats the internal ejecta and, after ≈ 1000 yr of evolution, refocuses approximately on the z-axis near the center of the explosion (see bottom left panel in Fig. 3). The almost perfect refocusing of the reverse shock is enhanced by the idealized cylindrical symmetry of the system and in particular of the two shells adopted. At this point, the central parts of the ejecta are confined by the material of the shells, while the terminal edges along the poles continue to expand at high velocity (see bottom left panel in Fig. 4). After the interaction with the two shells, the remnant continues to expand through the wind of the progenitor star (see bottom middle panels in Figs. 3 and 4). At the end of the simulation, namely after ≈ 10 000 yr of evolution, the ejecta have slowed down their expansion (see bottom right panel in Fig. 4) and they show an elongated shape due to the interaction with the two dense toroidal shells. At this time, the morphology is characterized by a broad jet-like structure with maximum density along the z-axis, which extends for ≈24 pc from the center of the explosion (see bottom right panel in Fig. 3). It is worth noting that, in our simulations, we assumed that the remnant propagates through an almost uniform ambient environment at large distances from the center of the explosion (namely because of the flattening of the wind density profile to 0.1 cm⁻³; see Sect. 2). However, at these distances, the remnant is expected to propagate through an inhomogeneous ISM that may partially wash out the fingerprints of the previous interaction of the remnant with the two dense shells.

Similarly to the case with low explosion energy reported in Figs. 3 and 4, the evolution of model M60-V0-9Foe is shown in Figs. 5 and 6 at different epochs (increasing from upper left to lower right panels). We note the different scales used in these figures in comparison with the analogous ones for model M60-V0-1Foe. The complete temporal evolution is available as online movies (Movie 3 and Movie 4). The total mass of the ejecta in model M60-V0-9Foe is much higher than in model M60-V0-1Foe, namely ≈13 M_⊙, as only ≈3 M_⊙ were lost in the fallback during the core-collapse event (see Table 3). In this case, the forward shock moves quickly through the CSM powered by the high kinetic energy of the ejecta and starts to interact much earlier with the inner dense shell (red clumpy semi-transparent structure in the top left panels in Figs. 3 and 5), namely after ~ 20 yr of evolution.As in model M60-V0-1Foe, the interaction determines a slowdown of the forward shock which travels through the shell (see top left panel in Fig. 6) and a strengthening of the reverse shock traveling through the ejecta, but in this case the effect is far less pronounced because of the high kinetic energy of the ejecta (see top left and middle panels in Figs. 3–6). The interaction of the remnant with the shell lasts for ~ 38 years. The forwardshock, represented with a semi-transparent surface in Fig. 5, then starts traveling again through the r⁻² wind of the progenitor. After ≈130 yr of evolution, the blast wave hits the outer shell (cyan semi-transparent clumpy structure in the top middle panel in Fig. 5). Also, in this case, the dense shell squeezes the ejecta and pushes them along the z-direction (see top right panels in Figs. 5 and 6). However, because of the high kinetic energy of the ejecta, the deceleration of the forward shock through the dense shell is not as relevant as in model M60-V0-1Foe and the blast moves quickly beyond the shell (see bottom panels in Fig. 6). As a consequence, the relative strengthening of the reverse shock is not as evident as in model M60-V0-1Foe (see top right panel in Fig. 3), and the shock refocuses on the z-axis much later. For instance, at t = 4200 yr the remnant has already reached a maximum expansion of ≈26 pc from the center (see bottom middle panel in Fig. 5) in model M60-V0-9Foe; in other words, the remnant is more extended than that of the analogous model with low explosion energy (M60-V0-1Foe) at the end of the simulation (see bottom right panel in Fig. 3). However, at the same age, the reverse shock in model M60-V0-9Foe has not refocused and the inner ejecta have not yet been heated. At the end of the simulation, namely at ≈ 10 000 yr, the ejecta have an elongated shape forming a broad jet-like structure with maximum density along the z-axis, which extends for ≈38 pc from the center of the explosion (see bottom right panel in Fig. 5).

The evolution and density distribution of the SNRs in models with the same explosion energy but different masses of the progenitor stars (models M60-V0-1Foe and M80-V0-1Foe) are very similar. The reason for this resemblanceis that the increase in mass of the progenitor is not reflected in the mass of the ejecta, as all the additional mass may have been previously expelled into the CSM or lost in the fallback immediately after the SN event (see second and third columns in Table 3). This is also true for models with low explosion energy and a nonzero rotation velocity (i.e., M60-V300-1Foe and M80-V300-1Foe). Thus, the explosion energy plays a crucial role in determining the morphology of the remnant during the interaction with the inhomogeneous CSM.

This is confirmed by the models with a high explosion energy. Specifically, we find that the remnant evolution and its morphology change appreciably from models M60-V0-9Foe and M60-V300-9Foe to models M80-V0-12Foe and M80-V300-12Foe because of ahigher explosion energy in the latter (see Table 3). As a consequence, the ejecta expand faster in models M80-V0-12Foe and M80-V300-12Foe than in the rest of the models and, thus, the reverse shock refocuses at later times. For instance, in model M80-V0-12Foe the internal ejecta have not been shocked by the reverse shock at the end of the simulation (see left panel in Fig. 7), whereas the reverse shock refocused at the center of the explosion at t ≈ 6500 yr in model M60-V0-9Foe. In model M80-V300-12Foe, the ejecta mass is only slightly higher than in model M60-V0-9Foe (see Table 3), and the reverse shock refocuses at t ≈ 9000 yr (instead of t ≈ 6500 yr as in M60-V0-9Foe).

Our simulations show that the remnant morphology is characterized by a large-scale asymmetry due to the interaction with the dense shells, which changes significantly during the evolution. We can describe the degree of asymmetry of the remnant by the parameter R = A_z∕A_xy, where A_z and A_xy are the distances of the forward shock from the center of the explosion along the z-axis and in the xy plane, respectively. In Fig. 8, we plot the parameter R versus the time of evolution (in logarithmic scale) for the different models explored (see Table 1). The figure shows that the models present three well-differentiated phases: (1) the initial free expansion through the stellar wind; (2) the interaction with the two dense shells; and (3) the later expansion through the wind. During the first phase, R ≈ 1, which means that the remnant is almost spherically symmetric: the homologous expansion of the remnant preserves the symmetry of the ejecta in all the models explored. In the subsequent phase, when the ejecta start to interact with the two dense shells (see Fig. 2), the remnant becomes progressively more asymmetric, reaching the maximum asymmetry after the interaction with the outer shell. In fact, when the forward shock starts to interact with the inner dense shell, the z-axis becomes the preferred direction of the remnant expansion and R increases. The growth of R stops shortly after the interaction with the inner shell, and then continues growing while interacting with the outer shell. In Fig. 8, we identify two groups of models: those with low and those with high explosion energy. Models with high explosion energy start to interact earlier with the dense shells and reach a maximum degree of asymmetry (~ 2) at t ≈ 240 yr (see Fig. 8; see also top right panel in Fig. 5); models with low explosion energy start to interact later with the dense shells and reach a maximum degree of asymmetry (~ 2.3) at t ≈ 500 yr (see Fig. 8; see also top right panel in Fig. 3). Models with different initial mass of the progenitor star but with the same explosion energy follow an analogous evolution. More specifically, our simulations show that models with the same explosion energy lead to very similar ejecta masses despite differences in the main sequence mass of the progenitor star (i.e., either 60 or 80 M_⊙; see Table 3). Since the ejecta mass (and not the main sequence mass of the progenitor star) is relevant in the evolution of the remnant for the cases explored here, the explosion energy turns out to be the most relevant factor in the dynamical evolution of the SNR. We also note that models with higher explosion energy evolve on shorter timescales and reach a lower degree of asymmetry compared to those with lower explosion energy; the timescale of evolution reflects the expansion velocity of the remnant, which is higher in models with higher explosion energy. On the other hand, the degree of asymmetry of the remnant depends on both the energy of the explosion and the density contrast of the shells. The changes in the degree of asymmetry of the remnant are less evident than the changes in the timescale because, for the asymmetry, a central role is played by the density contrast of the shells, which is the same in all the models considered here. In any case, models where the other parameters are the same and either with or without rotation follow similar evolution. In the last phase, the ejecta continue to expand through the ambient medium in all directions while again becoming more symmetric. In this phase, the parameter R gradually decreases showing the progressive reduction of the asymmetry. In all cases, the large-scale morphology of the SNR retains an imprint or “memory” of the early interaction of the remnant with the inhomogeneous CSM for almost 10 000 yr (namely the period covered by our simulations).

Fig. 3 Density distributions of the ejecta of the model M60-V0-1Foe at different evolution times (increasing from upper left to lower right). The opaque irregular isosurfaces correspond to a value of density at 1% of the peak density with one quadrant cut in order to see the radial distribution. The semi-transparent surface marks the position of the forward shock; the initially toroidal semi-transparent structures in red and cyan colors represent the inner and outer shells in the CSM respectively. The system is oriented as G26 (see Fig. 2), corresponding to the rotation angles ix = 30°, iy = 30°, iz = 25° about the x, y, and z axes, respectively. The evolution time is shown in the upper left corner of each panel. Az indicates the distance of the forward shock from the center of the explosion along the z-axis. The complete temporal evolution is available online (Movie 1).

Fig. 4 Same as in Fig. 3 but the colors give the radial velocity in units of 1000 km s−1 on the isosurface. The complete temporal evolution is available online (Movie 2).

Fig. 5 Same as in Fig. 3 but for model M60-V0-9Foe. The complete temporal evolution is available online (Movie 3).

Fig. 6 Same as in Fig. 5 but the colors give the radial velocity in units of 1000 km s−1 on the isosurface. The complete temporal evolution is available online (Movie 4).

Fig. 7 Density distributions for the ejecta of the models M80-V0-12Foe (left panel) and M80-V300-12Foe (right panel) after ~ 10 000 yr of evolution. The isosurfaces and the orientation of the system are defined as in Fig. 3.

Fig. 8 Degree of asymmetry of the remnant morphology, R, versus time for the models explored (see Table 1): M60-1Foe (black), M80-1Foe (green), M60-9Foe (blue), M80-12Foe (red). In each case, solid lines represent models with Vrot =0, and dotted lines the ones with Vrot = 300 km s−1.

3.2 Mass distribution in velocity space

HD instabilities that develop during the interaction of the ejecta with the reverse shock determine the mixing of the shock-heated ejecta in the region between the forward and the reverse shocks, and thus the distribution of the chemical elements at different evolution times of the SNR. The mixing between layers of different chemical composition during the evolution is reflected in the velocity distributions at different ages of the SNR. Figures 9 and 10 show the mass distributions of selected elements in models M60-V0-1Foe and M60-V0-9Foe, respectively, versus the radial velocity, V_rad (first column), and the velocity along the LoS when the point of view lies either on the y-axis, V_y (second column), or on the z-axis, V_z (third column), at different evolution times (increasing from top to bottom). In all the figures presented, DM_i is the mass of the ith element in the velocity range [v, v + dv], where dv = 100 km s⁻¹ is the velocity binning, and M_i is the total mass of the ith element. In the upper row, we present the distributions after the homologous expansion of the SNR through the stellar wind and immediatelybefore the interaction with the innermost toroidal shell; in the middle row, we show the distributions after the interaction with the two shells when the SNR reaches its maximum degree of asymmetry (see Fig. 8); finally, in the lower row, we plot the distributions at the end of the simulation.

Before the interaction with the two toroidal shells, the mass distributions are similar to those of the initial condition, and the distributions versus V_y and V_z are almost the same (see upper row in Figs. 9 and 10). This is a result of the spherically symmetric initial conditions considered and the homologous expansion of the SNR through the stellar wind. When the remnant starts to interact with the dense shells, the mass distributions of the species versus the LoS velocity is different if the point of view lies on the y or on the z axis. In the first case, the species can reach a maximum velocity which is significantly smaller than in the second case because of the slowdown of the ejecta caused by the interaction with the shells. If the point of view lies on the z axis, the ejecta appear to expand faster, with the z-axis being the preferential direction of expansion. These differences between the mass distributions versus V_y and V_z are maintained also after the interaction with the dense shells (see Figs. 9 and 10), although they reduce following the decrease in the degree of asymmetry of the remnant (see Fig. 8). We note that the mass distributions versus the LoS velocities V_y and V_z are symmetric with respect to the zero velocity. This is due to the idealized density distribution of the CSM (in particular the two nested shells) adopted here, which is characterized by cylindrical symmetry about the z-axis (see Fig. 2 and Sect. 2.2). An asymmetry in the density distribution of the shells would be reflected in asymmetric mass distributions versus the LoS velocities V_y and V_z.

In model M60-V0-1Foe, most chemical elements follow similar shapes and are distributed in a broad maximum below ≈ 10 000 km s⁻¹ after 31 yr of evolution (see upper row in Fig. 9), immediately before the interaction with the innermost dense shell (see left top panel in Fig. 3). The exceptions are the unstable ⁴⁴Ti and ⁵⁶Ni (and their decaying products ⁴⁴Ca, ⁵⁶Co, and ⁵⁶Fe at late evolution times) which follow a narrower distribution with maximum below ≈ 4000 km s⁻¹. However, we note that in all the models with low explosion energy the total masses of ⁴⁴Ti and ⁵⁶Ni (and, therefore, their decaying products) are very low in comparison with lighter elements (see Table 2), as heavy elements go primarily to the fallback during the collapse of the star (Limongi & Chieffi 2018). The distributions of intermediate-mass and light elements at velocities larger than ≈ 7000 km s⁻¹ have similar shapes with a slope much steeper than in the initial condition (soon after the shock breakout). The steepening of their slopes is a sign that these ejecta have already passed through the reverse shock and their similar shape is a sign of efficient mixing in the region between the reverse and forward shocks. The ejecta continue their expansion slowed down as a consequence of their interaction with the dense CSM. After interacting with both dense shells (see top right panel in Fig. 3), the external layers of the ejecta have a maximum radial velocity of ≈ 7000 km s⁻¹ (see left middle panel in Fig. 9). At this time, the remnant reaches its maximum degree of asymmetry (see Fig. 8) with its preferential direction of expansion along the z-axis, being V_y < 4000 km s⁻¹ for all the species (see central and right middle panels in Fig. 9). At t ≈ 473 yr, the reverse shock starts to interact with the internal layers of the ejecta with velocities lower than ≈ 2000 km s⁻¹ and finally focuses on the z axis at t ≈ 1000 yr. Then, the ejecta continue to expand in all directions lowering their degree of asymmetry (see Figs. 3, 4 and 8). At the end of the simulation, V_rad < 2000 km s⁻¹, V_y < 1500 km s⁻¹, and V_z < 2000 km s⁻¹.

In model M60-V0-9Foe, the mass distributions versus the velocity are less homogeneous than in models with low explosion energy. Most low- and intermediate-mass elements are distributed in a broad maximum below ≈ 20 000 km s⁻¹ after 16 yr of evolution (see upper row in Fig. 10), immediately before the interaction with the innermost dense shell (see left top panel in Fig. 5). The ¹⁴N, ⁴⁴Ti, and ⁵⁶Ni instead follow a narrower distribution;in particular, for ⁴⁴Ti and ⁵⁶Ni (and their decaying products) which are present mainly in the internal layers of ejecta, the maximum velocity of propagation is ≈ 5000 and ≈ 3000 km s⁻¹, respectively (see upper row in Fig. 10). In this case, the similarities shown by most of the distributions of low- and intermediate-mass elements at velocities larger than ≈ 16 000 km s⁻¹ (see upper row in Fig. 10) indicate an efficient mixing in the region between the reverse and forward shocks. Similarly to the case with low explosion energy, the expansion of the ejecta in the equatorial plane is slowed down during theirinteraction with the dense CSM, and the remnant reaches its maximum degree of asymmetry at t ≈ 240 yr (see Fig. 8 and right top panel in Fig. 5). At this time, V_rad < 15 000 km s⁻¹, V_y < 8000 km s⁻¹, and V_z < 15 000 km s⁻¹ (see middle row in Fig. 10). The reverse shock traveling through the ejecta has not yet reached the internal layers, dominated by ⁴⁴Ti and ⁵⁶Ni and their decaying products (see middle row in Fig. 10). After this interaction phase, the ejecta continue to expand in all directions lowering their degree of asymmetry (see Figs. 5, 6, and 8). At the end of the simulation, V_rad < 3000 km s⁻¹, V_y < 2000 km s⁻¹, and V_z < 3000 km s⁻¹. We find velocities slightly higher than in the case with low explosion energy at t ≈ 10000 yr, but in model M60-V0-9Foe, Az = 38 pc (in M60-V0-1Foe, Az = 24 pc at the same time) and the remnant has already reached a low degree of asymmetry (see Fig. 8 and bottom right panels in Figs. 3 and 5).

The mass distributions for the rest of models explored with low and high explosion energy are very similar to those presented for models M60-V0-1Foe and M60-V0-9Foe, respectively. In Fig. 11, we plot the mass distributions versus V_rad for models M80 after the interaction with the two toroidal shells, at the time they reach their maximum degree of asymmetry (see Fig. 8). Models M80-V0-1Foe and M80-V0-12Foe presented in the left panels of Fig. 11 show very similar distributions to those presented for models M60-V0-1Foe and M60-V0-9Foe, respectively (see left middle panels in Figs. 9 and 10. In the analogous models but with an initial rotation velocity of the star (see right panels in Fig. 11), the only difference we observe is the higher abundance of ¹⁴N present in the He core due to the effect of the mechanical instabilities induced by rotation (Limongi & Chieffi 2018).

Fig. 9 Mass distributions of 4He, 12C, 14N, 16O, 20Ne, 24Mg, 28Si, 40Ca, 44Ti, and 56Ni versus the radial velocity, Vrad (first column), the velocity along the y-axis, Vy (second column), or the z-axis, Vz (third column), for model M60-V0-1Foe at different evolution times (increasing from top to bottom). Top row: immediately before the interaction with the innermost shell. Middle row: at the maximum degree of asymmetry after the interaction with the two shells. Bottom row: near the end of the evolution.

3.3 Spatial distribution and chemical composition of the ejecta

We explored the chemical evolution of the ejecta for the ten different species considered in this work (⁴He, ¹²C, ¹⁴N, ¹⁶O, ²⁰Ne, ²⁴Mg, ²⁸Si, ⁴⁰Ca, ⁴⁴Ti, ⁵⁶Ni). The most abundant species in all the models analyzed are ¹⁶O and ¹²C (see Table 2). The total mass of the ejecta is much lower than the mass of the fallback in all the models with low explosion energy (see Table 3). Thus, in all these cases, the masses of the heaviest species, namely the ⁴⁴Ti and the ⁵⁶Ni (and, therefore, their decaying products) are very low (see Table 2). In models with low explosion energy, we therefore find quite homogeneous distributions of low- and intermediate-mass elements in the ejecta, with a very residual presence of heavier elements.

In models with high explosion energy instead, we find a less homogeneous distribution of chemical elements in the ejecta. In these cases, the mass of the ejecta is much larger than the mass of the fallback (see Table 3) and a significant amount of ⁴⁴Ti and ⁵⁶Ni (and their decaying products) is present in the internal layers (see Table 2). In Fig. 12, we show the density distribution for the ejecta rich in Fe⁸, Si, and O for the four models with high explosion energy explored at different evolution times (increasing from top to bottom). After the interaction with the two toroidal shells, at t ≈ 200 yr (see upper panels in Fig. 12) the four distributions show similar characteristics: the O, distributed in the external layers of the ejecta and already heated by the reverse shock, forms an elongated structure along the z axis; the Si (the spherical green surface) is starting to interact with the reverse shock (the yellow transparent surface visible in the upper panels) and the Fe (the orange spherical most internal surface) is still unshocked. In model M80-V300-12Foe (see right top panel in Fig. 12), the total mass of Fe is higher (see Table 2) and the Si has already been partially shocked by the reverse shock. The forward and reverse shocks continue to propagate in opposite directions and after ~ 2000 years of evolution the ejecta start to form a Si-rich jet-like structure along the z-axis in most of the models (see middle panels in Fig. 12). In model M80-V0-12Foe, this occurs later (~ 4000 yr) due to the higher mass of the ejecta and the slower propagation of the reverse shock (see Sect. 3.1). At the end of the simulation, all the models with high explosion energy show Fe-rich internal ejecta distributions surrounded by an elongated Si-rich structure, itself surrounded by a more diffuse O-rich ejecta (see lower panels in Fig. 12). The Fe-rich internal structure is more extended and elongated in models M80.

The distribution of the material that we were able to observe in Fig. 12 can be very different depending on the LoS considered. For instance, when the LoS lies on the equatorial plane (the xy plane), the ejecta projected in the plane of the sky form a jet-like structure rich in Fe and Si surrounded by a more diffuse O-rich material in the external layers. When integrating along the z-axis instead, the projected ejecta form a spherical central structure with mixed composition surrounded by a more diffuse O-rich ring. Finally, in the case of the LoS of G26 (see Fig. 12), we expect to observe a structure of projected ejecta intermediate between the two extreme cases described above.

Fig. 10 Same as in Fig. 9 but for model M60-V0-9Foe.

Fig. 11 Mass distributions of 4He, 12C, 14N, 16O, 20Ne, 24Mg, 28Si, 40Ca, 44Ti, and 56Ni versus radialvelocity, Vrad, for models M80 at the time they reach their maximum degree of asymmetry. The figure shows models either with low (upper panels) or high (lower panels) explosion energy, and models either without (on the left) or with (on the right) initial rotation velocity of the star.

4 Summary and discussion

We investigated what the SNR would be like if the LBV candidate G26 exploded as a SN. To this end, we modeled the evolution of a LBV exploding as a core-collapse SN, starting immediately after the shock breakout. We then followed the transition from the SN to the SNR phase and the interaction of the remnant with the inhomogeneous pre-SN environment for ≈ 10 000 yr and for eight different parent SNe from LBV-like stars compatible with G26 (see Sect. 2.2). The SNR models differ from each other in terms of explosion energy and the characteristics of the progenitor star (see Table 1), which have been selected from the pre-SN models investigated by Limongi & Chieffi (2018) and prescribed as initial conditions in this work. The pre-SN environment is the same for all the models and is based on infrared and radio observations of G26 (see Fig. 2).

Our models show three well-differentiated phases in the evolution of the remnant: (1) its expansion through the innermost CSM, modeled here as a r⁻² stellar wind;(2) its interaction with the dense CSM characterized by two dense toroidal shells; and (3) the expansion of the remnant through an almost uniform ambient environment. After the homologous expansion through the stellar wind, at t ≈ 18–40 yr the forward shock starts to interact with the innermost dense toroidal shell characterizing the CSM (see Fig. 2). The interaction determines a strong slowdown of the forward shock in the equatorial plane, and a strengthening of the reverse shock traveling through the ejecta (see Figs. 4 and 6). During this phase the z-axis becomes the preferred direction of expansion of the blast wave and the remnant progressively becomes more asymmetric, reaching its maximum degree of asymmetry (R) after the interaction with the outermost dense shell (see Fig. 8) at t ≈ 200–500 yr. Finally, the remnant continues to expand through the wind while R decreases. After ≈ 10 000 yr of evolution the ejecta have slowed down their expansion and they show an elongated shape forming a broad jet-like structure with maximum density along the z-axis, which extends for ≈ 24–38 pc from the center of the explosion (see Figs. 3 and 5). These jet-like structures are the relic of the early interaction of the remnant with the inhomogeneous CSM that was sculpted by violent mass-loss events occurring in the latest phases of evolution of the progenitor star.

Among the cases explored (see Table 1), there are two classes of models: those with low explosion energy (10⁵¹ erg) and those with high explosion energy (9–12 × 10⁵¹ erg; see Fig. 8). The initial mass of the progenitor and the explosion energy determine the mass, the expansion velocity, and the chemical composition of the ejecta. Models with higher explosion energy evolve on shorter timescales and reach a lower degree of asymmetry compared to those with low explosion energy. More specifically, in the latter models, the remnant starts to interact later with the dense shells but is more affected by the CSM (which limits the ejecta motion in the equatorial plane and triggers a strengthening of the reverse shock traveling through the ejecta) and reaches a higher degree of asymmetry (see Fig. 8). We note that the explosion energy influences the timescale of evolution more than the remnant asymmetry because the timescale reflects the expansion velocity of the remnant (which is higher for higher explosion energy), whilst the asymmetry depends on both the explosion energy and the density contrast of the shells, with the latter being the same in all the models. Models with low explosion energy all have a relatively low total mass of ejecta (< 4 M_⊙) as most of the material was previously expelled into the CSM or lost in the fallback during the SN event (see Table 3). As a consequence, the mass of the heavier species (namely ⁴⁴Ti and ⁵⁶Ni and their decaying products) is very low in all the models with low explosion energy (see Table 2). The spatial distribution of the low- and intermediate-mass elements is quite uniform in all the models with low explosion energy, with ¹²C and ¹⁶O being the most abundant from the species considered. All the models show an efficient mixing in the region between the reverse and forward shocks (see Fig. 9). The only difference we observe between analogous models either with or without rotation is the higher abundance of ¹⁴N present in the He core due to the effect of mechanical instabilities induced by rotation in the models with V_rot = 300 km s⁻¹ (see Fig. 11; see also Limongi & Chieffi 2018).

In the models with high explosion energy, the expansion of the ejecta in the xy-plane is also slowed down due to the interaction with the dense shells but this effect is less pronounced than in models with low explosion energy. As aconsequence, the degree of asymmetry reached in these models is lower (see Fig. 8) and the reverse shock refocuses at later times (see Fig. 5). The models with high explosion energy show a higher mass of ejecta and aless homogeneous chemical distribution of ejecta than the models with low explosion energy. After ≈ 2000–4000 yr of evolution,the reverse shock starts to interact with the innermost layers of ejecta, forming a Si-rich jet-like structure along the z-axis. At the end of the simulation (t ≈ 10 000 yr), all the models with high explosion energy show Fe-rich internal ejecta distributions surrounded by an elongated Si-rich structure, itself surrounded by more diffuse O-rich ejecta (see lower panels in Fig. 12). The Fe-rich internal structure is more extended and elongated in models M80, namely those with 80 M_⊙ on the ZAMS. We stress here that these elongated features originate from the interaction of the remnant with the inhomogeneous ambient medium and do not reflect large-scale asymmetries left from the earliest phases of the SN explosion and in some cases developed from stochastic processes (convective overturn and the standing accretion shock instability) during the first seconds of the SN blast (e.g., Wongwathanarat et al. 2015, 2017). In all the cases examined, we therefore found that the remnant morphology retains an imprint or “memory” of the early interaction of the remnant with the inhomogeneous CSM (the shells), even thousands of years after the SN. In other words, the effects of mass-loss events occurring in the latest phases of the progenitor star evolution could still be encoded in the asymmetries of remnants up to 10 000 yr old (the time covered by our simulations).

In the present study, we examined the characteristics of remnants of core-collapse SNe from LBV stars, which show a strong interaction with their CSM. Given the complex structure and high density of the CSM around LBVs, the characteristics of the CSM where a LBV goes SN play a fundamental role in determining the properties of its remnant. Fixing the structure of the CSM consistently with that inferred from observations of G26, we find a common morphology of the SNR for all the progenitors and explosion energies explored: elongated ejecta with an internal jet-like structure, which is the result of the interaction with the highly inhomogeneous CSM.

It is worth emphasizing that we adopted the LBV candidate G26 as a template because the structure and density distribution of the massive nebula in which the star is located are well characterized by observations (see Umana et al. 2012). This allowed us to define an idealized geometry for the dense CSM of G26 (see Fig. 2), consisting of two axially symmetric shells, each with ~10 M_⊙. As expected, the two dense shells play a central role in modifying the expansion of the forward shock and in driving a reflected shock through the ejecta. However, the CSM structure adopted here is specific to G26, but a wide range of masses and different geometries may characterize LBVs. With the information gathered from our simulations, we now know what to expect if the mass of the shells or the geometry of the CSM are different from those of G26. If we consider the same geometry adopted here but more massive shells, we expect the forward shock to be more slowed down by the interaction with the denser shells, a stronger reflected shock driven through the internal ejecta, and a remnant with a morphology that is similar to that found here but with a higher degree of asymmetry. Conversely, if the remnant interacts with less massive shells, we expect the forward shock to expand quickly through the low-density shells, a fainter reflected shock driven backward through the ejecta, and a remnant morphology with a lower degree of asymmetry. As for the geometry of the CSM, we find that the two axially symmetric shells characterizing the CSM of G26 lead to a broad jet-like structure in the morphology of the remnant. Thus, we expect that large-scale inhomogeneities in the pre-SN CSM, possibly different from that adopted here, may be reflected in the final morphology and large-scale asymmetries of the remnant.

We expect SNRs from LBV progenitors to be relatively uncommon, because LBVs are considered to be brief transitional phase in the evolution of the quite rare most massive stars. Nevertheless, there are a few examples of remnants in the literature showing morphologies analogous to those modeled here, namely W50 (Dubner et al. 1998), SNR G309.2-00.6 (Gaensler et al. 1998), SNR W44 (Shelton et al. 2004), and SNR S 147 (Drew et al. 2005; Gvaramadze 2006). In all these cases, the remnants appear elongated along a preferential direction of expansion, and in some cases they present some hints of a jet-like structure. The morphology characterized by “ears” of these SNRs has been interpreted by some authors as the result of jet-driven core-collapse SN mechanisms (see Grichener & Soker 2017). If that were the case, the remnant could be inflated by jets that are launched during the explosion. Here, we show that similar morphologies can also be reproduced by the interaction of the remnant with a highly inhomogeneous and dense CSM such as that surrounding LBVs. This could therefore be the case for some of the remnants cited above. The role of the CSM in shaping the remnant morphology was also recently investigated by Chiotellis et al. (2021).

Luminous blue variable stars as explosive transients constitute a puzzling and still poorly understood category. According to our models, a distinctive property of SNe from LBV progenitors could be the significant fallback of matter soon after the core-collapse, especially for those with a low-energy explosion. This implies the formation of a black hole instead of a neutron star. In this case, we do not expect the SNRs from LBV progenitors to normally host a detectable compact object (a neutron star) in their interior. LBVs are also extremely interesting SN progenitors because of the strong interaction between their ejecta and the pre-existing slower and dense CSM. In this work, the pre-existing CSM is described following the two-shell environment identified in G26 by Umana et al. (2012), which in turn is a simplified version of the CSM expected in a LBV exploding star. Considering a more realistic description of the star wind and the CSM close to the progenitor in future models could account for the high luminosities observed during the early phases of Type IIn SNe (see Smith 2017a).

The interaction of the blast wave with the dense shells can also have important consequences for the acceleration of cosmic rays (CRs). During the interaction, the remnant can become a strong γ-ray source and may provide evidence for hadronic acceleration of CRs (Bykov et al. 2018). Furthermore, these remnants can be factories of very energetic particles possibly up to PeV energies. These very energetic particles are observed in the galactic CR spectrum but are not inferred from observations of SNRs. As the blast wave from the SNe of LBVs is very energetic during the first decades of evolution of the remnant, it could be possible that PeV particles are produced in the interaction of the remnant with the dense shells (e.g., Zirakashvili & Ptuskin 2016). In this case, we again expect γ-rays that could be detected with current and future instruments (e.g., the Cherenkov Telescope Array). Self-consistent models suchas those analyzed here could allow us to disentangle the effects of interaction of the remnant with an inhomogeneous CSM from those of the structure of the progenitor star. This could shed some light on the last phases of evolution of massive stars and on the link between LBVs and Type IIn/IIb SNe, which remains weak.

Fig. 12 Density distributions for the ejecta rich in Fe (orange), Si (green), and O (blue) for the models with high explosion energy at different evolution times (increasing from top to bottom). The opaque irregular isosurfaces correspond to a value of density at 1% of the peak density in each case. The semi-transparent surface marks the position of the forward shock; the yellow semi-transparent surface barely visible in the top panels represents the position of the reverse shock. The orientation of the system is the same as G26 (see Fig. 2).

Movies

Movie 1: Density distribution for the ejecta of the model M60-V0-1Foe. The opaque irregular isosurface corresponds to a value of density which is at 1% of the peak density with one quadrant cut in order to see the radial distribution. The semi-transparent surface marks the position of the forward shock; the initially toroidal semi-transparent structures in red and cyan colors represent the inner and outer shells in the CSM, respectively. The system is oriented as G26, corresponding to the rotation angles i_x = 30°, i_y = 30°, i_z = 25° about the x, y, and z axes, respectively. Access here

Movie 2: Same as in Movie 1 but for the colors giving the radial velocity in units of 1000 km s⁻¹ on the isosurface. Access here

Movie 3: Same as in Movie 1 but for model M60-V0-9Foe. Access here

Movie 4: Same as in Movie 3 but for the colors giving the radial velocity in units of 1000 km s⁻¹ on the isosurface. Access here

Acknowledgements

We thank the referee for useful comments and suggestions that allowed us to improve the manuscript. We acknowledgethe computing centre of Cineca and INAF, under the coordination of the “Accordo Quadro MoU per lo svolgimento di attivit‘a congiunta di ricerca Nuove frontiere in Astrofisica: HPC e Data Exploration di nuova generazione”, and the HPC facility (SCAN) of the INAF – Osservatorio Astronomico di Palermo for the availability of computing resources and support. The PLUTO code, used in this work, was developed at the Turin Astronomical Observatory in collaboration with the Department of General Physics of Turin University and the SCAI Department of CINECA. We acknowledge financial contribution by the INAF PRIN 2019 grant “From massive stars to supernovae and supernova remnants: driving mass, energy and cosmic rays in our Galaxy” and the INAF mainstream program “Understanding Particle Acceleration in Galactic Sources in the CTA era”.

References

All Tables

All Figures

Fig. 1 Left upper panel: Pre-SN evolution of stars of 60 M⊙ and 80 M⊙ (from the main sequence phase up to the core He depletion stage, see Limongi & Chieffi 2018), with solar (initial) metallicities and initial rotation velocities Vrot = 0 km s−1 (solid lines) and Vrot = 300 km s−1 (dotted lines). The yellow asterisk marks the position of G26 in the HR diagram. Right upper panel: initial radial profiles of density in the models defined in Table 1. Models with initial rotation velocities Vrot = 0 km s−1 and Vrot= 300 km s−1 are represented with solid and dotted lines respectively. Lower panels: initial radial profiles of mass fractions for the species considered for two representative models, one with a low explosion energy (left panel) and the other with a high explosion energy (right panel). In both cases, the model describes a ZAMS star of 60 M⊙. The range of enclosed mass (x-axis) is very different in these models because of the different location of the mass cut (see Limongi & Chieffi 2018).

In the text

	Fig. 2 Top: map of G26 which combines infrared at 24 μm (blue) and radio (red) images (see Umana et al. 2012). Bottom: schematic view of the CSM around G26, representing the two dense shells (in red and light blue) observed in the radio and infrared bands. The LBV is located in the center of the two nested toroidal shells.
In the text

Fig. 3 Density distributions of the ejecta of the model M60-V0-1Foe at different evolution times (increasing from upper left to lower right). The opaque irregular isosurfaces correspond to a value of density at 1% of the peak density with one quadrant cut in order to see the radial distribution. The semi-transparent surface marks the position of the forward shock; the initially toroidal semi-transparent structures in red and cyan colors represent the inner and outer shells in the CSM respectively. The system is oriented as G26 (see Fig. 2), corresponding to the rotation angles ix = 30°, iy = 30°, iz = 25° about the x, y, and z axes, respectively. The evolution time is shown in the upper left corner of each panel. Az indicates the distance of the forward shock from the center of the explosion along the z-axis. The complete temporal evolution is available online (Movie 1).

In the text

	Fig. 4 Same as in Fig. 3 but the colors give the radial velocity in units of 1000 km s−1 on the isosurface. The complete temporal evolution is available online (Movie 2).
In the text

	Fig. 5 Same as in Fig. 3 but for model M60-V0-9Foe. The complete temporal evolution is available online (Movie 3).
In the text

	Fig. 6 Same as in Fig. 5 but the colors give the radial velocity in units of 1000 km s−1 on the isosurface. The complete temporal evolution is available online (Movie 4).
In the text

	Fig. 7 Density distributions for the ejecta of the models M80-V0-12Foe (left panel) and M80-V300-12Foe (right panel) after ~ 10 000 yr of evolution. The isosurfaces and the orientation of the system are defined as in Fig. 3.
In the text

	Fig. 8 Degree of asymmetry of the remnant morphology, R, versus time for the models explored (see Table 1): M60-1Foe (black), M80-1Foe (green), M60-9Foe (blue), M80-12Foe (red). In each case, solid lines represent models with Vrot =0, and dotted lines the ones with Vrot = 300 km s−1.
In the text

Fig. 9 Mass distributions of 4He, 12C, 14N, 16O, 20Ne, 24Mg, 28Si, 40Ca, 44Ti, and 56Ni versus the radial velocity, Vrad (first column), the velocity along the y-axis, Vy (second column), or the z-axis, Vz (third column), for model M60-V0-1Foe at different evolution times (increasing from top to bottom). Top row: immediately before the interaction with the innermost shell. Middle row: at the maximum degree of asymmetry after the interaction with the two shells. Bottom row: near the end of the evolution.

In the text

	Fig. 10 Same as in Fig. 9 but for model M60-V0-9Foe.
In the text

	Fig. 11 Mass distributions of 4He, 12C, 14N, 16O, 20Ne, 24Mg, 28Si, 40Ca, 44Ti, and 56Ni versus radialvelocity, Vrad, for models M80 at the time they reach their maximum degree of asymmetry. The figure shows models either with low (upper panels) or high (lower panels) explosion energy, and models either without (on the left) or with (on the right) initial rotation velocity of the star.
In the text

Fig. 12 Density distributions for the ejecta rich in Fe (orange), Si (green), and O (blue) for the models with high explosion energy at different evolution times (increasing from top to bottom). The opaque irregular isosurfaces correspond to a value of density at 1% of the peak density in each case. The semi-transparent surface marks the position of the forward shock; the yellow semi-transparent surface barely visible in the top panels represents the position of the reverse shock. The orientation of the system is the same as G26 (see Fig. 2).

In the text