© L. Dessart et al. 2021

1 Introduction

The late time radiative properties of core-collapse supernovae (SNe) contain information on the progenitor star at the time of death and the mechanism by which it exploded. This topic has been extensively discussed in the past (for a review, see Jerkstrand 2017). In Dessart & Hillier (2020a), we presented a general study of nebular phase spectra of Type II SNe. Although this study conveyed some insights, it suffered from the adoption of a simplistic-toy setup for the ejecta. The treatment of chemical mixing was inconsistent with current theoretical expectations since the need to account for macroscopic mixing generally meant that we also enforced microscopic mixing, which is unphysical. A trick to avoid this problem was to decouple the mixing of ⁵⁶Ni and that of other species, which was again inconsistent. In Dessart & Hillier (2020b), we introduced a new technique for treating chemical mixing of SN ejecta in grid-based radiative-transfer codes. The method is inherently 1D and consists in shuffling spherical shells of distinct composition in mass space in an ejecta in homologous expansion (e.g., at one year after explosion). At that time, the velocity is essentially proportional to the radius so that this shuffling corresponds to a mixing of material in velocity space while preserving the original chemical mixture for each mass shell. As an introduction to this research topic and a discussion of our modeling technique are given in Dessart & Hillier (2020a,b), they are not repeated here.

In this work, we extend our previous exploration, which was limited to the s15A model of Woosley & Heger (2007, hereafter WH07), to study 23 explosion models taken from WH07 and Sukhbold et al. (2016, hereafter S16). This grid covers the range of progenitor masses most likely to explode (i.e., 9–26 M_⊙) and has two different treatments of the explosion. In WH07, the explosion is generated with a piston whose trajectory is designed to deliver an ejecta kinetic energy of about 10⁵¹ erg at infinity. In contrast, S16 simulates the evolution of the collapsing star through core bounce, shock formation, and the stalling of the shock aswell as its revival through a consistent handling of the neutrino-energy deposition in the infalling mantle. Although the method retains some shortcomings (spherical symmetry, calibration to a core-collapse explosion “engine” meant to reproduce SN 1987A), the approach gains considerably in physical consistency and delivers an unbiased prediction of Type SN II explosion properties. The goal of the present study is to document the nebular-phase spectral properties of these ab initio explosion models fromS16, while the WH07 models are used to test the sensitivity of results to the numerical handling of the explosion.

Our work provides a picture complementary to the results obtained, for the most part, with the code SUMO (Jerkstrand et al. 2011, 2012), which uses a Monte Carlo treatment of the radiative transfer. The probabilistic approach is well suited to handle the stochastic propagation of photons in a complex mixture of clumps whose distinct composition reflects the macroscopically mixed progenitor metal-rich core. Such studies of nearby Type II SNe suggest a lack of high-mass red-supergiant (RSG) progenitors (Jerkstrand et al. 2012, 2014, 2015b; Maguire et al. 2012; Yuan et al. 2016). These studies were all targeted toward the modeling of a specific Type II SN. Using a given explosion model (usually taken from WH07), the ejecta was adjusted so that the expansion rate of the metal-rich layers was compatible with the width of [O I] λλ 6300, 6364, Hα, or [Ca II] λλ 7291, 7323. The adjustment to the initial model was tailored independently for each progenitor mass and guided by the need to match the kinematics inferred from the observations. The present work is therefore complementary to this approach since we take the actual ejecta models of S16 and investigate the range of spectral properties they exhibit at nebular times.

In the next section, we present our numerical setup as well as the explosion models used as initial conditions for the radiative-transfer modeling with CMFGEN (Hillier & Dessart 2012; Dessart & Hillier 2020b). Section 3 presents the qualitative and quantitative results for the simulations, with a particular focus on how the power absorbed in the various ejecta shells is radiated in nebular emission lines. Section 4 discusses the origin of Ca II emission in our TypeII SN simulations and compares our results with those from previous work. The treatment of the ⁵⁶Ni -bubble effect, presented in Sect. 5, reveals only a minor impact on the SN radiative properties. Similarly, tests for material clumping, presented in Sect. 6, reveal a weak sensitivity over the parameter space explored. While our simulations adopt by default a complete macroscopic mixing of the metal-rich core, which is a conservative assumption, we find that the lack of such mixing, as would occur in a quasi-prompt explosion, would significantly impact the nebular spectral properties (Sect. 7). Section 8 discusses the morphology of the strongest optical emission lines in nebular phase spectra of Type II SNe, in particular in connection to line overlap and optical depth effects. A comparison to a selection of observations is presented in Sect. 9. We give our concluding remarks in Sect. 10.

2 Numerical setup

Using the same numerical approach as in Dessart & Hillier (2020b), we study the radiative properties of a diverse set of explosion models derived from 9–29 M_⊙ stars on the main sequence. We consider two sets of explosion models that were taken from WH07 and S16. Both sets were evolved with KEPLER at solar metallicity and without rotation until core collapse. The physics in the KEPLER code has been described by Weaver et al. (1978) and Woosley et al. (2002). In the “s-series” of models used here, convection was included in the Ledoux formalism with a substantial amount of semiconvection (Weaver & Woosley 1993). The rate for ¹² C (α, γ)¹⁶O corresponded to an S factor at 300 keV of 175 keV b, which is moderately high compared with some recent estimates (deBoer et al. 2017) but within experimental error bars. Nucleosynthesis was tracked using an adaptive network with all isotopes appropriate to a given stage (Rauscher et al. 2002) – up to 2000 in the explosion itself.

Upon reaching core collapse, the WH07 models were exploded with a prescribed piston trajectory so that the asymptotic kinetic energy would be about 10⁵¹ erg. In contrast, the S16 models were exploded with the P-HOTB code, which uses a “neutrino engine” whose strength was calibrated to reproduce the elementary properties of the well-studied SN 1987A and SN 1054 (the Crab explosion; Ugliano et al. 2012; Ertl et al. 2016; S16). Between these cases of a higher-mass and a lower-mass progenitor, respectively, the values of the engine parameters were interpolated as a function of stellar core compactness (see S16 for details). Effectively, the engine parameters determine the energy release of the cooling and accreting proto-neutron star (PNS), which, in turn, determines the neutrino luminosity of the high-density PNS core. Neutrino transport in the accretion mantle of the PNS is treated by a gray approximation. The neutrino-driven SN models employed in the present study are spherically symmetric. However, the neutrino-driven mechanism is known to be a multidimensional phenomenon. Although some multi-D effects may be accounted for in an effective way by using the engine calibration, other effects cannot be described by 1D models. For example, long-lasting simultaneous accretion downflows and re-ejected outflows around the PNS in 3D simulations (see Müller et al. 2017a, Bollig et al. 2021) are replaced by an extended post-bounce period of PNS accretion and, after a considerably delayed onset of the explosion, a subsequent neutrino-driven wind, whose strength is overestimated in order to provide a neutrino-heated mass with an energy that is needed to power the SN. Clearly, explosion asymmetries and large-scale mixing processes, which are a natural consequenceof nonradial hydrodynamic instabilities during explosions in 3D, cannot be described in spherical symmetry.

Neutrino-driven explosions modeled in 1D display a highly non-monotonic behavior of the explosion and remnant properties as a function of progenitor mass or compactness (see Ugliano et al. 2012, Ertl et al. 2016, S16). This behavior and most of its consequences were confirmed by independent numerical modeling approaches in 1D (e.g., Pejcha & Thompson 2015; Ebinger et al. 2019) as well as a semi-analytic description that accounts for 3D effects in a parametrized way (Müller et al. 2016). Although some authors express concerns that 3D effects in 1D models are lacking (e.g., Couch et al. 2020), the results of these models are compatible with a number of observational constraints, e.g., galactic chemical abundances, a rough correlation between explosion energy and ⁵⁶Ni mass (Müller et al. 2017b), the spread of observed neutron star masses (within uncertainties), and the observed rarity of Type II SN progenitors above 16–20 M_⊙ (e.g., Smartt 2009; Sukhbold & Adams 2020, and references therein). Since the exact outcome of the 1D modeling depends on the strength of the neutrino engine, which was calibrated in the S16 explosion simulations for different suggestions of the SN 1987A progenitor, further validation will have to come from observations. Ultimately 3D models will need to be employed. However, 3D explosions are computationally expensive, are unsettled, for example, due to uncertainties in nuclear and neutrino physics, and in the properties of the progenitor. Additionally, there is a much larger parameter space making it more difficult to link any one model with an observed SN. Thus 1D models, guided by results from 3D models, still have an important role to play in advancing our understanding of core collapse SNe.

For the first set of models, we selected the WH07 progenitors with an initial mass on the zero-age main sequence (ZAMS) of 12, 15, 20, 25, and 29 M_⊙. In the nomenclature of WH07, these models are s12A, s15A, s20A, s25A, and s29A. The “A” suffix refers to models in which the ejecta kinetic energy at infinity is in the vicinity of 1.2 × 10⁵¹ erg and where the mass cut was placed at the location in the progenitor where the dimensionless entropy S∕N_Ak_B = 4¹. Simulations for other masses were not made because the explosion models of WH07 retain some level of arbitrariness – the models are artificially exploded. For the second set of models, taken from S16, we selected 18 successful explosions (computed with the W18 engine) arising from stars with a ZAMS mass between 9 and 26.5 M_⊙, only excluding a few with which we encountered convergence difficulties in the radiative-transfer calculation. In the nomenclature of S16, these models are s9, s9p5, s10, etc. (where a “p” replaces the dot for noninteger masses). As is apparent from the next section, our sample of S16 models already exhibits some degeneracy in SN properties (i.e., nearly identical explosion and SN radiation properties), so we believe the present grid of models covers the range of nebular properties currently available from state-of-the-art explosion models. For the WH07 models s12A, s15A, s20A, and s25A, there is a model in the S16 sample with the same or a similar ZAMS mass (s12, s15p2, s20p1, s25p2; the models, however, differ in explosion energy, ejecta mass, or yields, including ⁵⁶Ni mass), which allows us to evaluate the impact of using different prescriptions for the explosion and pre-SN evolution (both sets were calculated with KEPLER about ten years apart).

A summary of model properties is given in Table 1, with some properties also shown in Figs 1 and 2. Both WH07 and S16 models show the same trend of increasing O-rich shell mass (defined here as all ejecta depths with an O mass fraction greater than 10%) with progenitor mass. Similar trends would be seen for the He-core or the CO-core mass. The S16 progenitors less massive than 12 M_⊙ on the ZAMS have an O-rich shell mass that drops down to about 0.1 M_⊙ for the lowest-mass massive stars undergoing Fe-core collapse (the minimum mass is not firmly established; Poelarends et al. 2008). The explosions from these lower-mass stars eject mostly H and He, contribute little to the metal enrichment of the Universe, and thus may not exhibit strong metal lines in their associated SN nebular spectra.

The WH07 models have about the same ejecta kinetic energy of 1.2 × 10⁵¹ erg irrespective of ZAMS mass, with a monotonic increase in the ⁵⁶Ni yield with ZAMS mass that reflects the monotonic increase in the progenitor He-core mass or the compactness (this simply reflects the larger amount of dense material above the Fe core at collapse; see S16 for a discussion). This uniform E_kin is not a prediction of the model but is imposed through the prescribed piston trajectory. In practice, this trajectory could have been tuned further so that the E_kin were exactly the same for all models in the “A” series of WH07, irrespective of initial mass. This freedom is not ideal since the explosion energy is thought to depend sensitively on the structure in the inner ~ 2 M_⊙ of the progenitor star at core collapse. This shortcoming is cured in the S16 models by means of a more consistent calculation of the collapse, bounce, and post-bounce evolution until explosion. In these calculations, the progenitor structure is an essential ingredient controlling the asymptotic ejecta kinetic energy (the S16 models still have shortcomings, such as the assumption of spherical symmetry or the calibration of the explosion engine; see the discussion earlier in this section).

In the S16 models, a trend toward higher E_kin and ⁵⁶Ni yield with increasing ZAMS mass can be discerned. There is a modest scatter for the ⁵⁶Ni yield, but there are many outliers for the distribution in E_kin. The non-monotonicity of E_kin with progenitor mass can be understood as follows. On the one hand the neutrino energy input to the explosion grows for more massive neutron stars, which have a tendency to be formed in the bigger He-cores of the most massive of the considered progenitors. On the other hand, the binding energy of the ejected stellar core material grows for stars with bigger He-cores. Because of the saturation of the explosion energy around 2 × 10⁵¹ erg (S16), the asymptotic ejecta kinetic energy attains a maximum at intermediate stellar masses.

While the models (in particular those drawn from the S16 set) have a higher physical consistency than previous ones, one must bear in mind that other explosion models could have been produced by varying the properties of the progenitor on the ZAMS (metallicity, rotation rate), or by adopting different numerical values for parameters that control the progenitor evolution (mixing length, mass loss prescription, overshoot etc.). These variations will alter the progenitor properties, such as the size of the CO and He cores and the extent of the H-rich envelope, at core collapse. Similarly, different explosion engines (see S16 for details) lead to different ejecta properties, most importantly E_kin and the ⁵⁶Ni mass.

Taking the selected WH07 and S16 models at 350 days after explosion, we used the method discussed in Dessart & Hillier (2020b) to introduce macroscopic mixing without any microscopic mixing. The technique amounts to shuffling shells of equal mass but located in different parts of the metal-rich unmixed ejecta at 350 days. By design, this shuffling preserves the kinetic energy, the total mass, and all individual element masses. Only shells within a Lagrangian mass M_sh take part in this shuffling – beyond M_sh the original ejecta composition from WH07 or S16 is left untouched. Because the ejecta is in homologous expansion at this late time, the shell shuffling in mass space is equivalent to a shell shuffling in velocity (or radial) space. Provided there is enough shuffling of the metal-rich shells relative to the ⁵⁶Ni -rich material, the exact shell arrangement in the inner metal-rich ejecta is unimportant (as demonstrated by Dessart & Hillier 2020b). With this provision in mind, what counts is how far out in the ejecta we place the ⁵⁶Ni. To apply a consistent setup for the whole grid of models, we set M_sh ≈ M_He,c + 1.5 M_⊙. This implies that the metal rich core is always “fully” mixed and that about 1.5 M_⊙ of material from the overlying H-rich shell is mixed into the metal-rich core. Because the He-core mass increases with progenitor mass, this criterion implies a greater mixing in velocity (or mass) space (this statement would strictly hold if the final mass at core collapse was independent of ZAMS mass). It also means that in higher mass progenitors, the H-rich material mixed inward represents a smaller fraction of the surrounding metal-rich material. The weakness of the method, which also applies to any similar works, is that we do not know the precise level of mixing in any given SN, nor how it varies with increasing ZAMS mass. The 1D approach also ignores the possibility of large-scale asymmetry, which probably exists at some level in all core-collapse SNe. We show the shuffled-shell composition structure of the s9 and s15p2 models of S16 in Fig. 1.

The CMFGEN simulations assume steady state, and typically use 350 points for the radial grid, uniformly spaced in velocity below M_sh and uniformly spaced on an optical depth scale above that. We treat the ⁵⁶Ni two-step decay chain and ignore the contribution from other unstable isotopes. Nonlocal energy deposition is computed by solving the gray radiative transfer equation (the γ-ray gray absorption-only opacity is set to 0.06 Y_e cm² g⁻¹, where Y_e is the electron fraction). Nonthermal processes are treated as per normal (for details, see Li et al. 2012; Dessart et al. 2012). All simulations presented use a turbulent velocity of 10 km s⁻¹. However, with the new shuffled-shell technique (which does not introduce any artificial microscopic mixing), we find that the choice of turbulent velocity is not as critical for the present simulations of Type II SNe at nebular times (Dessart & Hillier 2020a). Indeed, tests show that the spectra obtained for 10 km s⁻¹ are essentially the same as those obtained for 50 km s⁻¹ for the present set of models.

3 Landscape of nebular phase properties from the S16 grid of explosion models

3.1 Qualitative description

Figure 3 shows the optical spectra obtained with CMFGEN at an epoch of 350 days after explosion. The spectra are normalized so that the maximum flux in the illustrated wavelength range is set to one. The models from the WH07 set are a good place to start the discussion since these models are characterized by a very similar ejecta kinetic energy of 1.1–1.3 × 10⁵¹ erg, and a moderate range in ⁵⁶Ni mass between 0.05 and 0.13 M_⊙. Hence, one expects that the key spectral differences exhibited in the WH07 set primarily reflect the differences in progenitor composition, especially the systematic increase in the O-shell mass with the mass of the progenitor (factor of ~ 10 increase between s12A and s29A).

The WH07 set exhibit a monotonic increase in the strength and width of the [O I] λλ 6300, 6364 doublet line from s12A to s29A; in the same order, the Hα line weakens and broadens, though only most obviously for the highest-mass model s29A; finally, the [Ca II] λλ 7291, 7323 line stays about the same and is never the strongest optical line. Qualitatively, these trends follow expectations and the results from previous calculations (see, for example, Jerkstrand et al. 2012). The [O I] λλ 6300, 6364 line strength in Type II SN spectra increases with progenitor mass because the O-rich shell captures an increasing fraction of the total decay power absorbed by the ejecta and because [O I] λλ 6300, 6364 is the main coolant for the O-rich material. Given the narrow range in ejecta kinetic energy and mass (see Table 1), the variation in line widths reflects the chemical stratification of the progenitor star (Dessart et al. 2010): with increasing progenitor mass, the CO core occupies a growing fraction of the total ejecta mass, which tends to place the O-rich material at larger velocity – the same feature explains the broad Hα in the s29A model. The lack of variation in the [Ca II] λλ 7291, 7323 line strength confirms that it is a poor indicator of progenitor mass (this result is discussed in detail below).

In the larger S16 model set, the spectral evolution shows greater variations and more complicated trends compared to the WH07 model set, with an obvious non-monotonic dependence of the relative line fluxes on progenitor mass. There are several reasons for this. First, the S16 model extends to lower-mass massive stars. S16 predicts a low explosion energy of 1–6 × 10⁵⁰ erg for initial masses 9–12 M_⊙ (with the s12 model having half the kinetic energy of model s12A). Consequently, the nebular spectra of these models exhibit narrower lines than in the rest of the S16 sample or relative to our selection of WH07 models. Second, S16 predict low ⁵⁶Ni masses of 0.004–0.03 M_⊙ for these low-energy explosions. This leads to a fainter SN at nebular times, a reduced heating rate, and potentially lower ejecta ionization (Table 2). Third, the lower expansion rate implies a more compact ejecta at the same SN age, and hence a greater density relative to higher energy explosions. This greater density should favor recombination and hence reduce the ejecta ionization. Finally, the metal yields of lower-mass progenitors is smaller, which leads to a very weak [O I] λλ 6300, 6364 line relative to Hα. As a result of these properties, the S16 set of models s9–s12 stands well apart from the rest of the sample shown in Fig. 3.

For the mass range 12–25 M_⊙, each WH07 model has a counterpart in the S16 set (s12A and s12, s15A and s15p2, s20A and s20p1, s25A and s25p2). We compare the normalized optical spectra of each model pair in Fig. 4 – we omit the s15A and s15p2 pair since they are similar to each other and only slightly different from the 12 M_⊙ pair. Surprisingly, the optical spectra computed for models s12A and s12 are essentially identical, even though the progenitor evolution and the explosion were computed 10 yr apart and in a different fashion (see WH07 and S16 for a discussion). The shuffling is also done independently for each model and this results in a distinct composition structure. Furthermore, s12A has 1.65 times the ejecta kinetic energy of s12 (and roughly the same M_ej) but the only indication of this is the offset in the P Cygni trough associated with Hα and Hβ. Otherwise, emission lines have essentially the same width. The weak sensitivity to the factor of 1.65 difference in E_kin arises because the bulk of the kinetic energy is contained in the outer ejecta layers, which contribute little to the spectrum formation region at 350 d. Instead, the spectrum forms deeper in the ejecta, where the actual offset in expansion rate for the metal-rich layers is more modest². Model s12A also has 60% more ⁵⁶Ni than model s12, and thus a 60% greater luminosity but after normalization, this offset is gone. This offset in ⁵⁶Ni mass is too small to drive a visible change in ionization. Spectral line widths at nebular times may also suffer from a degeneracy between explosion energy and mixing. A larger E_kin may enhance line broadening, but a similar effect may be caused by a greater mixing of ⁵⁶Ni to large velocities. The similar spectral properties at 350 days between models s12 and s12A suggest that slight differences in progenitor, explosion, or mixing (e.g., as treated via our shuffling technique) properties cannot be easily revealed. Such degeneraciesare important limitations of inferences based on nebular phase spectra.

The bottom two panels of Fig. 4 show a comparison of models s20A with s20p1, and of s25A with s25p2, whose ejecta properties (M_ej, E_kin, or yields, including ⁵⁶Ni mass) are very similar (Table 1). While some spectral regions (e.g., below 6000 Å) and some lines (e.g., Hα) appear similar in each pair of models, the strong [O I] λλ 6300, 6364 and [Ca II] λλ 7291, 7323 lines are significantly different. Specifically, the two S16 models exhibit a stronger [Ca II] λλ 7291, 7323 flux and a weaker [O I] λλ 6300, 6364 flux. This difference is due to the much larger Ca to O abundance ratio in the O-rich shell in the s20p1 and s25p2 models, a consequenceof the merging of the Si-rich shell and the O-rich shell during Si burning (see Dessart & Hillier 2020a, and references therein).

Due to the merging, the Ca mass fraction in the O-rich shell rises from ~ 10⁻⁵ to ~ 10⁻³ and the [Ca II] λλ 7291, 7323 doublet becomes a strong coolant of the O-rich material. This occurs even though Ca is still about 1000 times less abundant than O in the O-rich shell. The bulk of the Ca produced in massive stars is located in the Si/S shell (where the Ca mass fraction is about 0.05). Figure B.1 illustrates these composition properties for models s20A and s20p1, in particular the large Ca abundance in the Si/S shells and the contrast in Ca/O abundance ratios between the O-rich shells of these two ejecta models. Hence, despite the marginal change in the total Ca yield, Si-O shell merging can dramatically alter the nebular phase spectral properties and inferences drawn from them. In model s26p5, the Si-rich and O-rich shells merge, but only partly since Si-rich material contaminates only the inner parts of the O/Ne/Mg shell. Consequently, in that model, the [O I] λλ 6300, 6364 is not so visibly reduced and appears both strong and broad (bottom curve in Fig. 3).

Turning to the rest of the S16 models, we obtain very similar properties for progenitors with masses between 12.5 and 18.5. Variations are induced in part by the reduction in E_kin with increasing progenitor mass, which causes a reduction in line width for all strong lines (e.g., compare s12p5 with s18p5). However, in a relative sense, the [O I] λλ 6300, 6364 clearly strengthens relative to Hα in this sequence. Strikingly, these variations, which are related to the progenitor mass, are weaker than those arising from Si–O shell merging. The process of Si–O shell merging is indirectly related to progenitor mass (one expects a greater occurrence of this process in higher mass progenitors) but it is an irony that this effect dominates over that caused by the stark contrast in oxygen yields between lower-mass and higher-mass progenitors (Table 1 and Fig. 2).

Fig. 1 Illustration of the composition profile in the original model and after shuffling individual shells. Left: chemical composition versus Lagrangian mass in the s9 and s15p2 ejecta models of S16 at 350 days (the original mass fraction is shown for 56Ni) – ejecta properties are summarized in Table 1. The name and total yield of each selected species are labeled.An alternating white and gray background is used to distinguish between consecutive shells of distinct composition. Right: corresponding ejecta in which shells have been shuffled in mass space within the Lagrangian mass cut Msh = 2.0 and 5.0 M⊙. Other ejecta models exhibit similar profiles, only modulated by the original composition (in particular the mass of each dominant shell) and the choice of Msh (details on the shuffling procedure are given in the appendix of Dessart & Hillier 2020b).

Fig. 2 From top to bottom: ejecta kinetic energy, the mean expansion rate, the 56Ni mass, the mass of the H-rich shell (i.e., corresponding to the progenitor H-rich envelope), and finally the masses of the O-rich shell (which includes the O/Si, the O/Ne/Mg, and the O/C shells) and of the O/C shell in our selection of WH07 (circles) and S16 (stars) models. The Ekin is prescribed in the WH07 models, but computed ab-initio in the S16 models. This explains the more complicated trends in the latter (e.g., the systematically lower Ekin and M(56Ni) for masses above about 16 M⊙, with some exceptions with the Ekin of s25p2), which includes regions in original mass that lead to no explosion at all (hence these models are not included in our sample; see S16 for a discussion).

3.2 Quantitative description

Figure 5 presents a more quantitative description of the results shown in Fig. 3. For each model, we have measured the bolometric luminosity L_bol, the optical luminosity L_opt (it accountsfor the luminosity falling between 3500 and 9500 Å), the power deposited ė_dep in the various shells (the total power deposited in the ejecta ė_dep(tot) is just L_bol), and the power emitted in the strongest lines (i.e., [O I] λλ 6300, 6364, Hα, and [Ca II] λλ 7291, 7323).

To measure the power emitted in lines we undertake a formal solution of the radiative transfer equation ignoring all bound-bound transitions except those of the selected ion. The [O I] λλ 6300, 6364 line flux is then measured from the O I-only spectrum, the Hα line flux from the H I-only spectrum, and the [Ca II] λλ 7291, 7323 line flux from the Ca II-only spectrum. Because there are no neighboring lines that could corrupt the measurement, we can integrate the flux over a broad band (i.e., ±10 000 km s⁻¹ for Hα, ±8000 km s⁻¹ for the other lines) centered around the rest wavelength of the transition (the mean wavelength is used for doublet lines). This technique allows a direct comparison of the line flux with the decay power absorbed in each shell.

To facilitate future discussions we define the H-rich shell as all ejecta locations where the H mass fraction X_H> 0.3, the He-rich shell where X_He > 0.8, the O-rich shell where X_O > 0.6 (for models s9 and s9p5, which have a very low O-shell mass, the criterion is X_O > 0.2 because the peak O-abundance is reduced by the minor mixing we apply to avoid having an O mass fraction profile looking like a Dirac delta function), and the Si-rich shell where X_Si > 0.1.

The total decay power absorbed is L_bol (panel a in Fig. 5) and depends on the ⁵⁶Ni mass and the γ-ray escape fraction. The latter is a function of the ejecta density profile (which depends strongly on E_kin∕M_ej) and the location of the three ⁵⁶Ni-rich regions (there are three such regions because the original unmixed ejecta has only one and our shuffled-shell approach splits each dominant shell into three equal mass subshells). Our models tend to show a decreasing trapping of γ-rays toward lower masses (from s25A to s12A, it drops from 90% to 81%; from s21p5 to s12p5 it drops nearly monotonically from 95% to 67%). In the higher-mass model s29A, the trapping of decay power is only 64% of the total, primarily because of the reduced mass of the H-rich shells (i.e., there is less H-rich material to absorb γ-rays, and the metal-rich core has a lower density, which translates into a greater γ-ray mean free path). In the S16 models below 12 M_⊙ that have a low explosion energy, the trapping is complete. Overall, the trapping efficiency is always greater than 60%, which explains why the distribution of L_bol values (top left of Fig. 5) follows closely the distribution of M(⁵⁶Ni) (top panel of Fig. 2).

The fraction of the model luminosity that falls within the optical range spans the range 66 to 78% (panel b in Fig. 5), with a representative fraction of 71%. Most of the remaining flux is emitted in the IR while less than 5% is emitted in the UV. The optical luminosity at 350 d represents at most 70% of the total decay power emitted, but as low as about 50% in ejecta with the lowest trapping efficiency in our set (e.g., model s29A). The increase above 20 M_⊙ arises from a relative strengthening of the optical emission lines [O I] λλ 6300, 6364 and [Ca II] λλ 7291, 7323 (which fall in the optical range) as the metal-rich core increases in mass while the H-rich ejecta layers become progressively lighter. The 10% increase going from model s20A to model s29A arises from a combination of factors, including the drop by 60% in decay power absorbed, the drop in optical depth by a factor of three, and the drop from 9.7 to 4.2 M_⊙ of the H-rich envelope mass of the progenitor.

Panels c, d and e of Fig. 5 show the fractional decay power absorbed by the H-rich, O-rich, and Si-rich shells. For the H-rich shell, the power is primarily absorbed within the H-rich regions located below or immediately above M_sh, hence in regions that coexist with the O-rich and Si-rich shells and that are not too distant from the ⁵⁶Ni -rich material. In other words, the H-rich material at higher velocities contributes less. Because of optical depth effects, the power absorbed in a given shell does not equal the flux escaping from that shell (Dessart & Hillier 2020a). For example, in model s15p2, about 50% of the total power is absorbed by H-rich material and about 50% of the total escaping flux emerges from the H-rich material (the equality between the two is a coincidence), but a fraction of the power deposited inthe slower-moving H-rich shells is reprocessed by the overlying, faster-moving Fe-rich shells, while the outer H-rich shell, which absorbs about 20% of the total decay power, releases about 30% of the total flux (Fig. 6). Nonetheless, these three panels reveal a clear trend, reflecting the composition structure of the progenitors, and give indicative power limits for spectral lineluminosities.

As the progenitor mass is increased in both the WH07 and S16 model sets, the fractional power absorbed in the H-rich shells decreases (from 80% in s9 to 20% in s29A), On the other hand it increases in the O-rich shell (from nearly zero in s10 up to 50% in model s29A) and Si-rich shell (from nearly zero in s9 up to 40% in s25p2). The highest value reached is for model s25p2, and is partially caused by Si–O shell merging in the progenitor which boosts the mass of the Si-rich shell as defined by our criterion X_Si> 0.1). These trends can be easily understood by considering the variations in the mass of the three shells (Dessart & Hillier 2020a). The low metal yields and the dominance of the H/He shell in lower-mass massive stars naturally implies that the bulk of the decay powerwill be absorbed by H-rich material. As the ZAMS mass is increased, the O-rich and Si-rich shells capture a greater fraction of the decay power. Although the mass of the H-rich shell may be the same for a wide range of progenitor mass (Dessart & Hillier 2019) the ⁵⁶Ni-rich material is closer to the O-rich and Si-rich shells, which benefit from this proximity to absorb a greater fraction of the available decay power. However, for the highest-mass models (e.g., s29A) γ-ray escape from the inner ejecta disfavors the Si-rich shell relative to the massive O-rich shell. The large variations in decay power absorbed are at the origin of the variations in line fluxes, or at least they set the fundamental power limits for the line fluxes.

Panels f–g of Fig. 5 relate line flux to decay power absorbed. We find that the Hα line flux represents about 25% of the decay power absorbed in the H-rich shell. Some models exhibit a higher cooling power from Hα (i.e., s20A and s25A), and we believe this is caused by the somewhat weaker mixing in those models (M_sh ≈ M_He,c + 1 M_⊙) where less H-rich material is mixed with the metal-rich core (the inner H-rich shells absorb a smaller fraction of the total decay power). In that case, the Hα line flux comes from reprocessing of radiation from the inner ejecta. The lower cooling power of Hα in model s9 is caused by the strong Ca II cooling in the H-rich layers of that model (as discussed below, this is caused by the low ⁵⁶Ni mass, which fosters an ionization shift to Ca⁺, compared toCa²⁺ in most other models in our grid; see Table 2). The [O I] λλ 6300, 6364 line flux relative to the decay power absorbed in the O-rich shell covers a much larger range (from 30 to 90% for the bulk of models, with one outlier at 250%!). Model s10 appears as an outlier, and even seems to violate energy conservation. Instead, the apparently large cooling efficiency arises from the formation of the [O I] λλ 6300, 6364 doublet in both the O-rich shell and the H-rich shell (in the latter, it absorbs power that would otherwise have been radiated by [Ca II] λλ 7291, 7323).

The [O I] λλ 6300, 6364 is a strong coolant of the H-rich layers in model s10 in part because O is neutral there, while Ca is twice ionized (so that [Ca II] λλ 7291, 7323 cannot cool these layers). These ionization patterns arise from the larger E_kin of model s10, which causes a lower density for a comparable decay power emitted or absorbed in adjacent models. For example, in model s11, Ca is Ca⁺ and [Ca II] λλ 7291, 7323 cooling dominates that due to [O I] λλ 6300, 6364 in the H-rich layers. This effect is striking in lower-mass models because of the low mass of the O-rich shell (which strongly limits the [O I] λλ 6300, 6364 flux contribution from the O-rich shell). Other effects influencing the [O I] λλ 6300, 6364 emission is the ionization in the O-rich shell, which can affect the cooling power of O and Mg. Furthermore, Si-O shell merging is one way to quench [O I] λλ 6300, 6364 emission to the benefit of [Ca II] λλ 7291, 7323.

The last three panels h–i–j of Fig. 5 give the line fluxes of Hα, [O I] λλ 6300, 6364, and [Ca II] λλ 7291, 7323 as a fraction of the total optical luminosity (or flux). This normalization removes the need for an accurate determination of the ⁵⁶Ni mass, the γ-ray escape, or the distance to the SN – for observations all that is needed to build such ratios is the reddening. For lower-mass progenitors, the fractional Hα line flux is between 14 and 24%, steadily decreasing for initial masses beyond 15 M_⊙ to reach 6% for model s29A. In lower-mass progenitors, the scatter is caused by circumstances in which the radiation leakage in [Ca II] λλ 7291, 7323 is enhanced (models s9 and s11) or inhibited (model s10). For higher masses, the fractional Hα line flux decreases because of the decreasing decay power absorbed by the H-rich material and the decreasing optical depth of the outer H-richejecta.

With the exception of two outliers (those models with Si–O shell merging prior to core collapse), the fractional [O I] λλ 6300, 6364 line flux is an increasing function of the progenitor mass, growing from about 4% in model s9p5 to 25% in model s25A, with a maximum of 37% in model s29A. In the absence of Si–O shell merging, the fractional [O I] λλ 6300, 6364 line flux appears as the best indicator of progenitor mass in both model sets. In contrast to [O I] λλ 6300, 6364, the fractional [Ca II] λλ 7291, 7323 line flux exhibitsa lot of scatter, and shows no correlation with progenitor mass – progenitors with ZAMS masses differing by 20 M_⊙ can exhibit the same flux ratio. The strongest emitters are associated with model s9 (because of the dominance of Ca⁺ in the H-rich zones and because the H-rich and He-rich shells constitute over 97% of the ejecta mass) and models s20p1 and s25p2 (because of Si–O shell merging).

Fig. 3 Montage of spectra for the set of WH07 models (blue label) and S16 models (red label) at 350 days after explosion. The wavelength range is limited to the 5600–7500 Å region to reveal the variations in [O I] λλ 6300, 6364, Hα, and [Ca II] λλ 7291, 7323. We note the non-monotonic evolution of these line fluxes with increasing initial mass in the S16 models, in contrast to the results obtained for the WH07 models.

Fig. 4 Comparison between the S16 (blue) and WH07 (red) models of comparable main sequence mass, with s12 and s12A (top), s20p1 and s20A (middle), and s25p2 and s25A (bottom). In each panel, the flux is divided by the integral of the flux falling between 1000 and 25 000 Å (and then scaled by a factor of 103), which is similar to forcing all models to have the same total decay power absorbed (in reality, the optical flux differs between models of the same mass because of the differences in 56Ni or γ-ray escape; see Table 1). This normalization procedure also allows one to estimate by eye the cooling power of strong lines relative to the total decay power absorbed.

Fig. 5 Line fluxes, powers, and their ratios for our grid of radiative-transfer calculations based on the WH07 (circles) and S16 (stars) ejecta models. Each panel in discussed in turn in Sect. 3.2. A similar figure, including S16 models scaled in kinetic energy, is presented in Fig. A.1.

Fig. 6 Profiles for the composition of dominant species, the cumulative decay power absorbed, and the cumulative frequency-integrated flux versus depth. For the latter two, the sum is performed inward from the outermost ejecta location and normalized for visibility.

4 Origin of Ca II emission

Below we explore in more detail a variety of features discussed in the preceding sections. We start with the origin of the Ca II emission in SN II nebular spectra. As is well known, the [Ca II] λλ 7291, 7323 is a very strong coolant and may dominate the cooling of the gas (Fransson & Chevalier 1989; Li & McCray 1993). Its importance as a coolant depends on the Ca abundance, the Ca ionization, and the competition with cooling from other lines.

In their single zone modeling of SN 1987A, Li & McCray (1993) found that the [Ca II] λλ 7291, 7323 must arise from H-rich material (in which the Ca abundance is given by the LMC metallicity value) rather than from the more abundant and newly synthesized Ca created during the explosion. Jerkstrand et al. (2012) also found that the Ca II emission arises primarily from the H-rich material in their modeling of SN 2004et. In both cases, this H-rich material corresponds to the H-rich material at the base of the progenitor H-rich envelope, or to the H-rich material that was macroscopically mixed down to smaller velocities. An independent confirmation that at least some Ca II emission must arise from H-rich material is given by the observation of a small common bump in emission seen in Hα, the [Ca II] λλ 7291, 7323 and Fe II λ7155 (Spyromilio et al. 1993).

In our simulations for standard-energy Type II SNe, we find that in most models the [Ca II] λλ 7291, 7323 forms primarily in the Si-rich shell, hence where Ca is the most abundant. We found this in previous simulations adopting a simplified ejecta structure (Dessart & Hillier 2020a), in more sophisticated simulations using the shuffled-shell method applied to the s15A model of WH07 (Dessart & Hillier 2020b), and we again find this results in the present extended set of models based on the more physically consistent inputs from S16 (a few exceptions are discussed below). This is most strikingly revealed by Fig. 7, which illustrates the emission sites for the emergent optical radiation. While some [Ca II] λλ 7291, 7323 emission arises from H-rich zones, the bulk of it arises in the Si-rich zone. Type IIb/Ib/Ic SNe also exhibit [Ca II] λλ 7291, 7323 emission (e.g., Shivvers et al. 2019), and in this case it must come from newly synthesized Ca because their progenitor have a low-mass H-rich envelope or no H-rich envelope at all (e.g., Jerkstrand et al. 2015a).

In contrast to the [Ca II] λλ 7291, 7323 transition, the Ca II λλ 8498, 8542, 8662 triplet has a significant contribution from the H-rich shell. In model s12, for example, approximately 50% of the triplet emission originates in the H-rich shell. That the [Ca II] λλ 7291, 7323 and Ca II λλ 8498, 8542, 8662 have different emitting regions is not surprising – they have a different critical density, a different temperature dependence, and the Ca II λλ 8498, 8542, 8662 has a much larger optical depth.

Figure 8 reveals more information on the origin on line emission (bottom panel) by connecting it to both the abundance profile (top panel) and the ionization state (middle panel) for H, He, O, Ca, and Fe. In the Si-rich shell, Ca is once ionized and this is where the bulk of the emission is coming from. In all other shells, it is either twice ionized (Ca²⁺ dominates and Ca II emission is negligible) or underabundant (in the C/O shell). In model s15p2, the main deterrent for Ca II emission in the H-rich shell is therefore overionization (Table 2). This is further seen in Fig. 9where the dominant cooling processes are shown versus velocity. As is apparent, Ca II (together with S I) collisional excitation dominates the cooling of the Si-rich regions, while the main processes cooling the H-rich material are nonthermal excitation of H I and Fe II collisional excitation. This situation holds in essentially all models of the present grid, with only a few exceptions.

In three of our models [Ca II] λλ 7291, 7323 emission also comes from outside the Si-rich shell, giving rise to the outliers in the [Ca II] λλ 7291, 7323 line flux plot shown in the bottom left panel of Fig. 5. The extended emission arises from two effects. First, in the models that undergo Si-O shell merging in the final stages of evolution prior to collapse, the Ca II emission also arises from the O-rich region contaminated by Si-rich shell material. Hence, a fraction of the power deposited in the O-rich shell is then radiated by [Ca II] λλ 7291, 7323, which boosts its strength relative to models without Si–O shell merging. Secondly, a boost to the [Ca II] λλ 7291, 7323 emission can arise if Ca II collisional excitation becomes an important coolant for the H-rich material, as in model s9 (Fig. 10). Because of the low ⁵⁶Ni mass (0.004 M_⊙), the decay heating in s9 is reduced relative to a more standard explosion model like s12p5. Further, because of the low ejecta kinetic energy (1.1 × 10⁵⁰ erg; one tenth of that in model s12p5), the ejecta density is greater at late times. Both effects lead to lower ejecta ionization, and consequently, Ca is roughly singly ionized throughout the ejecta. As a consequence, and because a large fraction of the decay power is absorbed by the H-rich material, the [Ca II] λλ 7291, 7323 forms in both the Si-rich and H-rich layers. An analogous situation affects the O I emission, which comes in part from the H-rich layers in this model. Because of its extended spatial distribution in model s9, the [Ca II] λλ 7291, 7323 is very strong and broader relative to other optical emission lines (see Figs. 3 and 5).

Fig. 7 Illustration of the spatial regions (here shown in velocity space) contributing to the emergent flux in the S16 model s15p2. The grayscale image shows the observer’s frame flux contribution ∂Fλ,V ∕∂V (the map maximum is saturated at 20% of the true maximum to bias against the strong emission lines and better reveal the origin of the weaker emission) versus wavelength and ejecta velocity. The plot at left connects these emission contributions to ejecta shells and the normalized decay-power deposition profile (dashed line) on the same velocity scale. The upper panel shows the corresponding normalized flux Fλ integrated over all ejecta velocities.

Fig. 8 Top: illustration of ejecta and radiation properties for the shuffled-shell model s15p2. From upper to lower panels, we show the profiles versus velocity for the mass fraction of H, He, O, Ca, and Fe, their ionization state (zero for neutral, one for once ionized etc.; see also Table 2), and the formation regions for Na I D, [O I] λ 6364, Hα, He I 7065 Å, and [Ca II] λ 7291. The quantity ∫ ξ d log V is the line equivalent width, which here is normalized to unity for visibility.

5 ⁵⁶Ni bubble effect

Decay heating causes a time integrated energy injection that can approach the local specific kinetic energy of the ⁵⁶Ni -rich regions. In these regions, decay heating can impact the dynamics for days after shock breakout, at times when homologous expansion would otherwise hold (Woosley 1988; Herant et al. 1992; Basko 1994). In practice, these ⁵⁶Ni -rich regions can expand and create low-density “bubbles” surrounded by a dense wall of swept-up material, thereby triggering RT instabilities and additional mixing. Hence, this bubble effect implies a decrease in density of the ⁵⁶Ni -rich regions and a compression of the swept-up material that builds a thin dense wall around these newly created bubbles. This bubble effect depends on the 3D distribution and abundance of ⁵⁶Ni, as well as on the density and expansion rate of the surrounding material. This process is therefore very nonlinear and complex. Three-dimensional simulations over days and weeks after explosion suggest that the 3D ejecta can have density variations of about a factor of ten throughout the volume where the bubbles reside (Gabler et al. 2021). Such inhomogeneities can impact the SN radiation during the photospheric phase (Dessart et al. 2018; Dessart & Audit 2019) as well as the nebular phase properties of core-collapse SNe (Li et al. 1993; Jerkstrand et al. 2011).

In our 1D shuffled-shell method, we can implement this bubble effect by introducing a radial stretching of the zones containing ⁵⁶Ni, and a radial compression of all other zones within M_sh. For simplicity, we assume that this radial stretching and compression is uniform throughout M_sh.

We proceed in the following manner. In the original unmixed ejecta model, we initially compute the volume V_{Ni, 0} occupied by the ⁵⁶Ni -rich shell and the volume V_sh contained within M_sh. We then specify what fraction α_Ni of V_sh should be occupied by the ⁵⁶Ni-rich bubble in the mixed ejecta at late times (i.e., at 350 days in the present simulations). Numerical simulations suggest that α_Ni can be as large as 0.5, that is the bubble takes up half the volume occupied by the progenitor He core.

To conserve mass, the density of the ⁵⁶Ni bubble must be divided by a factor f_Ni = (α_NiV_sh)∕V_{Ni, 0} and the density of the other shells within M_sh must be divided by a factor f_other = (1 − α_Ni)V_sh∕(V_sh − V_{Ni, 0}). Because V_sh ≫ V_{Ni, 0} (in model s15A, the ⁵⁶Ni -rich shell represents only 2.5% of the total volume within M_sh immediately after explosion), this treatment implies a large density reduction of the ⁵⁶Ni -rich material but only a modest compression of the rest of the metal-rich regions (for the s15A model with a bubble having α_Ni = 0.6, the density in the ⁵⁶Ni -rich bubble drops by a factor of 24, while the rest of the material within M_sh is only compressed by a factor of 2.5). While the former is physical, the latter is probably not very realistic since the dense wall around the ⁵⁶Ni-rich bubble is expected to be very dense and the material further away hardly affected by the bubble effect. Hence, in reality, there should be a complicated 3D distribution of compressed and rarefied material within the metal-rich regions of core-collapse SNe.

Having determined the density scalings for the ⁵⁶Ni-rich shell and the other shells within M_sh, we proceed from the innermost zone outward through our shuffled-shell structure and apply a radial stretching or a radial compression to modify the volume of all shells as necessary. When we encounter a ⁵⁶Ni -rich zone of volume δV located between radii r_l and r_u on the original radial grid, we build a new radial grid and set: $$r_{\text{u}}^3=r_l^3+\frac{3}{4\pi }f\delta V,$$

where f is equal to f_Ni in ⁵⁶Ni -rich zones and f_other otherwise. The density of the corresponding zone is scaled by 1∕f so that the mass of the zone is unchanged. The radial grid and the density beyond M_sh are unchanged. Because of homology, the radial stretching or compression is equivalent to a similar distortion in velocity space. With the bubble effect, the ⁵⁶Ni-rich shells cover a larger range of velocities while the rest of the material is confined to narrower velocity regions.

Using the model s15A with M_sh = 5.36M_⊙ as a reference, we computed two other models to test the influence of the bubble effect adopting α_Ni= 0.3 and 0.6 (the ⁵⁶Ni -rich bubble occupies 30 and 60% of the total volume within M_sh). The results for the SN radiation and the ejecta properties are shown in Fig. 11. Interestingly, the impact on the SN radiation is weak, although some slight differences are visible in Ca II lines. Concerning ejecta properties, the bottom panel of Fig. 11 shows that the gas properties are hardly modified by the bubble effect, amounting mostly to a radial or velocity shift of the profiles depicted for the three models, with only an increase in ionization in ⁵⁶Ni -rich regions resulting from the lower density. As said earlier, while the bubble effect leads to a strong reduction of the density in ⁵⁶Ni -rich regions, the adopted uniform compression of the surrounding material is about a factor of two, which is too small to cause an ionization shift, in particular because the ionization is typically low for the dominant species of each zone.

The ⁵⁶Ni bubble effect should have a similarly weak effect in other models. In lower-mass models, the ⁵⁶Ni mass is small (down by a factor of ten or more relative to s15A), but the E_kin is also small (also down by up to a factor of ten relative to s15A) so the bubble effect should be comparable and therefore weak (as for model s15A). In higher mass models, the ⁵⁶Ni mass is within a factor of two of the value for s15A so we expect essentially the same behavior as seen in the tests performed for the model s15A.

Fig. 9 Illustration of the main cooling processes balancing the radioactive decay heating at all ejecta depths in the model s15p2. We show each dominant cooling rate (stepping down from the rate having the largest peak value at any depth) normalized to the local heating rate. The term “NT” stands for nonthermal processes (in this context nonthermal excitation and ionization) and “COL” stands for collisional processes (i.e., collisional excitation).

6 Influence of a uniform clumping

The ⁵⁶Ni-bubble effect described above is expected to be just one component contributing to the complex 3D structure of core-collapse SN ejecta. Besides the general density variations discussed in a simplified manner in the preceding section, the ⁵⁶Ni -bubble effect should lead to additional instabilities causing mixing and clumping of the material at the bubble interface with the surrounding medium. The instabilities driven by the bubble effect, and the likely interaction between distinct bubbles, should also build upon the instabilities associated with the explosion itself and the shock propagation in the progenitor envelope (Gabler et al. 2021). One thus expects a complex distribution of clump densities and composition distributed in a complicated 3D pattern.

To further explore the impact of density variations on the ejecta properties and escaping radiation, we introduce a uniform clumping of the ejecta. Following Dessart et al. (2018), we assume that the ejecta are composed of clumps that occupy a fraction f_vol of the total volume. The clumps are assumed to be small relative to a photon mean free path, and the medium surrounding these clumps is assumed to be void. With these assumptions the initial ejecta model density is simply scaled by a factor of 1∕f_vol, while opacities and emissivities (which are computed with populations and temperature deduced for the clumps) are all scaled by a factor of f_vol. These assumptions leave column densities unchanged but have a direct influence on processes that depend on the density squared (e.g., free-free), and an indirect influence on the radiation field because of the sensitivity of the kinetic equations to density.

Figure 12 shows the results for model s15A in the smooth density case and in the clumped case with a uniform volume filling factor of 0.3 and 0.1. These adopted values correspond to density compressions of 3.3 and 10.0, and may be considered representative of what is seen in 3D simulations of core-collapse SNe (see, e.g., Gabler et al. 2021)³. The impact of this magnitude of clumping on the spectral properties is modest. As expected, clumping causes a reduction in the ionization and the temperature of the gas, but the effect is not strong. It leads to a slightly stronger Ca II NIR triplet (as the Ca²⁺ transitions modestly to Ca⁺). The increased density in the clumped model favors collisional de-excitation, causing a weakening of forbidden lines like [O I] λλ 6300, 6364 and [Ca II] λλ 7291, 7323. The reduction in Mg ionization causes a strengthening of Mg I lines at 4571.1 Å and 1.502 μm. Similarly, the Fe I lines in the red part of the spectrum strengthen with enhanced clumping. For the highest level of clumping, the [O I] λλ 6300, 6364 does show a sizable reduction in flux (by up to 30% in the model with a 10% volume filling factor).

Overall, these changes are modest. A much greater level of clumping would be required to cause a significant ionization shift and a large impact on the spectrum. One limitation for the influence of clumping is that the ionization of the species dominating the cooling in each zone is almost optimal.⁴ For example, O is predominantly neutral in the oxygen emitting zones (see Table 2). By introducing strong clumping, the [O I] λλ 6300, 6364 line strength decreases because of the enhanced recombination of other species like Na or Mg, which take over some of the cooling through Na I D or Mg I 4751 Å. In the H-rich layers, Ca²⁺ still dominates and prevents cooling by [Ca II] λλ 7291, 7323. The H-rich zones might be significantly clumped if mixed inward into the metal-rich regions, but clumping of the H-rich material at large velocity is unlikely to be strong since there is no process to cause it.

Fig. 10 Same as Fig. 8, but now for the s9 model, which corresponds to the lowest-mass progenitor model in the S16 set.

7 Influence of mixing

In Dessart & Hillier (2020a), we explored the sensitivity of nebular-phase spectra to the adopted ejecta properties, in particular to the level of chemical mixing. Although the mixing procedure was not optimal in that study (since we crafted toy models for simplicity and flexibility), we found that the mixing of ⁵⁶Ni was the primary ingredient influencing spectral properties and in particular emission line fluxes. This arises because the location of ⁵⁶Ni determines the deposition profile of the decay power, tuning the fraction delivered to other ejecta regions. In the absence of mixing in a SN II, the bulk of the decay power would be absorbed in the inner ejecta (rich in metals) and little in the outer ejecta (rich in H and He). As ⁵⁶Ni mixing is enhanced,a greater fraction of the power is deposited further out, biasing against the inner ejecta regions. If only ⁵⁶Ni were mixed whileother shells remained unmixed, this could strongly impact how the ejecta cools and how the spectra appear.

In nature, one expects a significant mixing of the metal-rich core in core-collapse SNe (say below 2000–3000 km s⁻¹ in a standard-energy standard-mass SN II), in combination with the mixing of ⁵⁶Ni. Consequently, the near uniform deposition of decay power in the mixed metal-rich inner ejecta implies that the exact distribution of metals in the inner region is unimportant. We confirmed this expectation using our new mixing technique (i.e., the shuffled-shell approach) in Dessart & Hillier (2020b). These metal-rich regions are bathed in a sea of γ-rays, whose large mean free path allows them to fill the entire volume occupied by the inner ejecta.

The above picture is generally meant in a 1D sense – a spherical average of an explosion should present a high level of chemical mixing. However, numerical simulations of 3D neutrino-driven explosions exhibit a variety of explosion morphologies, all suggestive of strong mixing. For example, in the simulations of Gabler et al. (2021), the model B15 is characterized by a quasi-isotropic distribution of ⁵⁶Ni fingers while the model W15 has a few big ⁵⁶Ni protrusions but shows very weak ⁵⁶Ni mixing along some directions, sometimes covering a wide solid angle.

Our 1D treatment of mixing is probably not accurate since some radial directions may exhibit strong mixing while others may show little. Further, radiative transfer in spherical-averaged ejecta is not the same as performing 2D/3D radiative transfer in realistic 2D/3D geometries. Our approach applies to configurations in which the mixing is quasi-spherical but is inadequate for strongly aspherical composition distributions.

To test the influence of mixing in the present set of simulations, we rerun models s9, s9p5, and s11 using the initial unmixed ejecta conditions from S16. In low-energy explosion models, the shock propagation may be initiated sooner after core bounce, cutting short the development of fluid instabilities that are at the origin of much of the chemical mixing (see, however, Stockinger et al. 2020). Hence, if weak or absent mixing can occur in a general sense (i.e., be weak along all radial directions), it should be in this type of models. But weak mixing may occur along some radial direction in any SN II, so this test is useful to gauge the implications.

The influence of mixing in the model s9p5 is shown in Fig. 13 – models s9 and s11 are not discussed since they show a similar behavior. As discussed above, we see that going from the mixed to the unmixed model, the strength of [O I] λλ 6300, 6364 and [Ca II] λλ 7291, 7323 increases while the strength of Hα decreases. This occurs because a greater fraction of the decay power is absorbed by the metal-rich inner ejecta. Emission is biased toward the inner ejecta so that numerous emission lines appear narrower than in the mixed model s9p5. This is exemplified by the [Ca II] λλ 7291, 7323 line which shows distinct emission for each component of the doublet.

Fig. 11 Illustration of the influence of the bubble effect applied to the reference model s15A. We show the resultsfor the SN radiation and gas properties in model s15A and for its counterparts in which the 56Ni -rich material occupies 30 (model bubble0p3) and 60% (model bubble0p6) of the volume within Msh (see Sect. 5 for a discussion).

Fig. 12 Illustration of the influence of clumping, comparing the reference model s15A and its counterparts s15Afvol0p3 and s15Afvol0p1 in which a uniform volume filling factor of 0.3 and 0.1 is used.

8 Line profile morphology

The ejecta velocity in the spectrum formation region controls the width of the emission lines that are present in the spectra. In our model set, the representative ejecta expansion rate (which we take as $\sqrt{2E_{\text{kin}}/M_{\text{ej}}}$) covers the range from 1230 km s⁻¹ (model s9) to 3790 km s⁻¹ (model s12p5), thusonly a factor of three. Since the spectrum forms in the inner ejecta, one may expect a smaller contrast (say a factor of two) in line width between the weakest and the strongest explosions. In addition to its influence on the intrinsic line width, the ejecta expansion rate tunes the importance of overlap with other lines, either from the same atom or ion, or from other elements. The separation between the components of the O I and Ca II doublets are 3050 and 1330 km s⁻¹, so overlap between the O I doublet components will occur only in the energetic explosion models, while overlap of the Ca II doublet components will occur in all our models.

Figure 14 shows a montage of spectra for the lowest-energy explosion model s9 and the more standard energy explosion model s15p2, and zooming on the [O I] λλ 6300, 6364, the Hα, and the [Ca II] λλ 7291, 7323 spectral regions. In both models Hα is optically thick while the other two lines are marginally thin. In model s9, the total Rosseland-mean optical depth is 0.1 (of which 85% is due to electron scattering), so this model is essentially thin in the continuum, while in model s15p2, this optical depth is 0.4 and will affect spectrum formation. Hence, optical-depth effects cannot be neglected.

With our shuffled-shell mixing technique, the shells of distinct composition are staggered in velocity space so that optically thin emission lines forming in these shells tend to show a boxy profile with steps as one goes from the profile maximum to the blue and red wings. These steps are easily visible in the two doublets, but less obvious in Hα because it forms over a large ejecta volume. The artifact does not affect the total integrated line flux but only its distribution in wavelength.

In all lines shown, the peak emission occurs at or very near the rest wavelength of the line (or its doublet component). In model s9, each component of [O I] λλ 6300, 6364 is well separated but shows two distinct emissions, one occurring at very low velocity (causing a central peak) and an additional emission at larger velocity and producing a ledge on each side of the line. This [O I] λλ 6300, 6364 line flux sits on a quasi-continuum produced by a forest of overlapping Fe I lines. The smaller separation of the components of the [Ca II] λλ 7291, 7323 doublet line produces a profile with three peaks due to overlap in between the two components. In model s15p2, the higher expansion rate of all emitting shells leads to broader features, with overlap even for the [O I] λλ 6300, 6364 line. Other properties are similar to what we find for model s9.

As evidenced in the top two panels of Fig. 14, the flux due to O I exceeds the total flux. This is suggestive of optical depth effects, primarily due to overlap with Fe I lines in model s9, but probably due to both the influence of Fe I and continuum opacity in model s15p2 (the individual components show a red deficit). The large Hα line optical depth in both models is also at the origin of the extended blue absorption as well as the skewness of Hα. The extended Hα trough is caused by the outer, faster expanding H-rich ejecta along the line of sight, absorbing the overlapping Fe I emission. This reciprocal process (i.e., Fe I also attenuates Hα emission nearer line center) likely arises from the interweaving of H-rich and Fe-rich shells, with alternate emission and absorption in selective spectral bands.

In general, our spherically symmetric ejecta models do not exhibit systematic wavelength shifts, relative to the rest wavelength, of the location of maximum line flux (see Fig. 3), in apparent agreement with the observations of most SNe II (see Sect. 9). This is in contrast with SNe Ibc, which reveal systematic shifts (see, for example Taubenberger et al. 2009; Milisavljevic et al. 2010).

Fig. 13 Comparison of model s9p5 with the standard mixing procedure produced by the shuffled-shell technique (see Sect. 2) and the same model in which the original unmixed ejecta has been used without adjustments (model s9p5nomix).

9 Comparison to observations

The SN II models discussed in the preceding section have a specific set of properties, which depend on the progenitor star (primarily its mass), the final properties of the progenitor at core collapse, and its explosion as SN. While observed SNe occupy a wide parameter space in ⁵⁶Ni mass and explosion energy, our set of models define a specific pairing between progenitor mass and SN ejecta. This grid of models nonetheless permits a comparison to a sample of observed SNe. One does not expect a perfect match to any object but if the theoretical framework holds, the distribution of models and observations should compare favorably.

Our selection of observations, their inferred characteristics and associated references are given in Table 3. It includes low-luminosity SNe II-P (SNe 1997D, 2008bk, and 2016aqf), standard-luminosity SNe II-P (SNe 1999em, 2004et, 2012aw, 2013ej, 2015bs, and 2017eaw), and the Type II-peculiar SN 1987A. The inferred ⁵⁶Ni mass for these observations is generally obtained by assuming full trapping and estimating the fraction of the SN luminosity that falls outsideof the observed range. Combined with the uncertainty in SN distance and reddening, these ⁵⁶Ni masses are probably uncertain by a factor of two. In our simulations, our synthetic spectra exhibit mostly a bolometric offset and little variation in relative flux between different spectral regions when the ⁵⁶Ni mass varies by a factor of a few (see also Sect. 8 of Dessart & Hillier 2020a). Hence, an uncertainty or an offset in ⁵⁶Ni mass for both observations and models is of limited impact for the discussion that follows.

Fig. 14 Montage of spectra in the [O I] λλ 6300, 6364 region (top), the Hα region (middle), and the [Ca II] λλ 7291, 7323 region (bottom) for the s9 and the s15p2 models of S16. A thin vertical line indicates the rest wavelength of each component of the doublet linesand Hα. Colored lines indicate the flux arising exclusively from bound-bound transitions associated with a specific ion (other ions and their opacity are then neglected, which explains the flux difference and the excess flux in places). Notice the presence of [Ni II] 7377.8 and [Ni II] 7411.6 Å in the bottom panels.

9.1 Flux comparisons

9.1.1 General comments

The SN luminosity at nebular times is set by the ⁵⁶Ni mass (modulo γ-ray escape). To cancel the offset in ⁵⁶Ni mass between ourmodels and observations, we renormalize the spectra, so that they match at a given optical wavelength. As the ⁵⁶Ni mass can also influence the ionization and temperature of the gas we restrict our comparisons to cases where the offset in ⁵⁶Ni mass between model and observation is moderate. For standard-luminosity SNe II, the WH07 and S16 models have a ⁵⁶Ni mass that is within a factor of two of that inferred for our selection of observations. The offset between the model and observationsis typically larger for low-energy SNe II since our selected observations cover a range from 0.002 to ~ 0.009 M_⊙, while the ⁵⁶Ni mass in our low-energy explosion models is around 0.02 except for the model s9 (with ~ 0.004 M_⊙). In some cases, we opted for the better fit even if that meant using a model with a larger offset in ⁵⁶Ni mass.

The distribution of ⁵⁶Ni within the ejecta (and chemical mixing in general) can also impact line fluxes. The mixing formalism applied in our models is the same irrespective of progenitor mass: the shuffled structure in our models corresponds to a complete mixing of the metal rich core, with some mixing inward of 1–2 M_⊙ of material from the progenitor H-rich envelope. This level of mixing should hold in most SN II and so our models should be in the ballpark of what is required to match observations. What is lacking in our modeling is the consideration of large-scale asymmetries, which probably affect many, if not all, Type II SNe. As large-scale asymmetries require a 2D or 3D treatment, they cannot be addressed in the present work – they are left to a future study.

Using the set of observations logged in Table 3, we measure the flux falling between 4100 and 5500 Å to which we refer as Fe flux (it arises primarily from a forest of Fe II lines, with a smaller contribution from Fe I lines – the latter are more dominant beyond about 6000 Å), the [O I] λλ 6300, 6364 flux, the Hα flux, and the [Ca II] λλ 7291, 7323 flux. In contrast to the modeling Sect. 3.2, we cannot easily estimate the influence of overlapping lines, and hence we simply measure the flux of an emission feature over a fixed passband. We assume that the spectral range around a given line extends slightly beyond the half-width at half-maximum on both sides of that line, or is bounded by the location of minimum flux if there is a neighboring line (for example, as happens in between [O I] λλ 6300, 6364 and Hα). For the three strongest optical emission lines this ensures that we account for the bulk of the flux associated with the feature while minimizing contamination by neighboring emission lines (for example the flux from Ni II or Fe II lines immediately adjacent to [Ca II] λλ 7291, 7323). When treating the models, the measurement is performed in the same manner and is therefore affected by the same bias. These flux measurements, while quantitatively different from those measured with the alternate method of Sect. 3.2, show the same trends. For the discussion, we assume that the contribution to a given line stems primarily from the associated transition (e.g., the [O I] λλ 6300, 6364) and ignore any flux contribution by other species (e.g., Fe I).

The flux measurements (normalized by the flux falling between 4000 and 9000 Å; in cases where the spectrum covers a smaller range, we extrapolate using the mean flux near the spectrum edge) are shown for the selected models and observations in Fig. 15. The model and observed distributions of the measurements overlap, although, as discussed in Sect. 3.2, the non-monotonicity of the trends means that a given fractional flux measurement can be compatible with both a low-mass and a high-mass progenitor. In a large part, the region of large scatter in models is associated with low-mass progenitors, characterized by low metal yields but a wide range in explosion energy (or ejecta expansion rate). Furthermore, one needs to be cautious when comparing to observations since one match between a SN model and an observation in one plot (for one diagnostics), might correspond to a mismatch in another. So, we will avoid over-interpreting this figure.

Excluding the lower-mass progenitors, the fractional Fe and Hα luminosity behave in a similar fashion. This suggests the bulk of the Fe emission arises from the H-rich envelope (from Fe at the primordial metallicity) rather than the ⁵⁶Ni-rich layers. This makes sense since the H-rich shells are much more massive than the ⁵⁶Ni -rich layers (by a factor of 100 or more), which implies that it can absorb a large fraction of the decay power. Because of this mass offset, there is about as much Fe in the H-rich layers (where it is primarily in the form of Fe⁺) of SN II as in the zones that used to be ⁵⁶Ni rich (where it ispartially neutral). The scatter in the fractional Fe and Hα fluxes for the models and observations is similar.

Fig. 15 Flux measurements for the Fe II forest between 4100 and 5500 Å, Hα, [O I] λλ 6300, 6364 and [Ca II] λλ 7291, 7323, for our models (shown versus initial mass; symbols) and observations (colored bars at right). In all cases, we normalize to the total flux falling between 4000 and 9000 Å.

9.1.2 Ca II emission

The models show a “flat” distribution of [Ca II] λλ 7291, 7323 line flux as a function of progenitor mass, confirming that it is a poor diagnostics of progenitor properties. The observations fall within a narrow range for the normalized [Ca II] λλ 7291, 7323 line flux, between 0.09 and 0.17 (with SN 2016aqf standing a little off at 0.05). No observed SN in our sample exhibits a [Ca II] λλ 7291, 7323 line as strong as that obtained in models with Si-O shell merging. This process might still operate, but would require a more moderate contamination of Si-rich material in the O-rich shell. Similarly, we see that no low-luminosity SN in our sample exhibits the strong [Ca II] λλ 7291, 7323 line obtained in model s9. In that model, the large strength arises because Ca is once ionized in the H-rich layers, becoming an important coolant for this material. Hence, this may support our finding that the bulk of the [Ca II] λλ 7291, 7323 line emission in SNe II occurs in the Si-rich material, rather than from the H-rich zones.

9.1.3 O I emission

Probably the most interesting and intriguing result of our studies is that the observed fractional [O I] λλ 6300, 6364 line fluxes are compatible with our results for low-mass progenitors (9–12 M_⊙ mass range). The one exception is SN 2015bs, which is known to require a higher-mass RSG progenitor (Anderson et al. 2018). The match of observed [O I] λλ 6300, 6364 line flux ratio by lower-mass RSG progenitors is perplexing. It implies that [O I] λλ 6300, 6364 can be a strong coolant even for low or moderate mass progenitors. In Dessart & Hillier (2020a), we reportedthis result for our toy models, finding that the cooling efficiency of the [O I] λλ 6300, 6364 doublet was systematically higher (although exhibiting the same trend) than obtained in the simulations of Jerkstrand et al. (2015b). At that time, the offset might have arisen because of the simplistic approach, the simplified composition and the adopted model atoms. Here, we use state-of-the-art explosion models as initial conditions, in particular computed with a large nuclear network. We also apply a more suitable chemical mixing, introducing no artificial microscopic mixing.

There are two obvious physical reasons for this offset. First, in Sect. 6, we found that clumping reduces the [O I] λλ 6300, 6364 line flux by as much as 30% for the model with a 10% volume filling factor. While this adopted uniform clumping is probably too large, it suggests that adopting a more realistic ejecta structure (i.e., not smooth) will reduce appreciably the current prediction for the [O I] λλ 6300, 6364 line flux.

Another reason for the offset is our neglect of molecule formation and their potential cooling power. In particular, CO, if it exists, will contribute to the cooling of the O-rich material (specifically the O/C shell, and perhaps the He/C shell if O is abundant there), reducing the amount of power radiated by [O I] λλ 6300, 6364, the primary coolant for these ejecta layers (Fig. 9). In the observations of SN 2017eaw (Rho et al. 2018), the flux associated with CO molecules corresponds to a few percent of the optical flux, and therefore represents a sizeable fraction (a few tens of percent) of the [O I] λλ 6300, 6364 line flux. If we assume that only the O/C shell is dominated by CO cooling (see, e.g., Jerkstrand et al. 2014), the [O I] λλ 6300, 6364 flux in the optical will be reduced by about 25% in model s15p2. Near infrared observations of the fundamental and first overtone CO band are crucial for constraining the amount of CO cooling.

CO fundamental (day 117) and first overtone emission (day 112) was first detected in SN 1987A by Spyromilio et al. (1988). Modeling by these authors suggested a CO mass of about 5 × 10⁻⁵ M_⊙ at day 285, while Liljegren et al. (2020) estimated a much larger mass of 4 × 10⁻³ M_⊙. Because of the reduction in temperature in the C/O region, Liljegren et al. (2020) indicates that the O/C region will not be a significant contribution to the [O I] λλ 6300, 6364 in the optical region. SiO was also present from ~ 160 d, with an estimated mass of about 4 ± 2 × 10⁻⁶ M_⊙ near day 500 (Roche et al. 1991). Observations of the ejecta by the Atacama Large Millimeter/submillimeter Array (ALMA) reveal distinct clumpy and complex torus-like structures for both CO and SiO (Abellán et al. 2017) in SN 1987A. CO and SiO molecular formation is probably ubiquitous in Type II SNe after a few hundred days (Kotak et al. 2005, 2006).

Another potential source of error in our models is explosion asymmetries. An offset in ⁵⁶Ni from the O-rich material (as in a bipolar explosion) will reduce the energy deposition in O-rich material, and potentially weaken the [O I] λλ 6300, 6364 flux fraction.

Irrespective of our models there is a very limited range in the observed [O I] λλ 6300, 6364 flux fractions. The observed ratios differ by less than a factor of three (a factor of two if we exclude SN 2015bs). This can be compared with our models where the change in the flux ratio is a factor of four, with a factor of two change between progenitor masses of 9 and 12 M_⊙.

While determining an accurate progenitor mass is a goal of nebular modeling it is worth mentioning the weak dependence of the [O I] λλ 6300, 6364 flux fraction on the oxygen mass. While the oxygen mass varies by a factor of 50 in our models, and the oxygen mass fraction in the ejecta by a factor of 25, the increase in the [O I] λλ 6300, 6364 flux fraction is only a factor of 4. This insensitivity is not surprising – it arises because the emitted energy is set by the initial mass of ⁵⁶Ni, and because [O I] λλ 6300, 6364 is an important coolant, but only in the oxygen-rich shells. One can use the absolute [O I] λλ 6300, 6364 flux to estimate the oxygen mass, but this requires an accurate temperature estimate, and such an estimate is not easily available from the observations – it typically requires modeling. In principal [O I] λ5577 can be used with the [O I] λλ 6300, 6364 to constrain the temperature, but at the epochs considered in this paper it is very weak and badly blended.

Fig. 16 Left: comparison of model s9p5 (M(56Ni) = 0.017M⊙) at 350 d withthe observations of SN 1997D (M(56Ni) = 0.002M⊙) at 350 d. Right: comparison of model s11 (M(56Ni) = 0.016M⊙) with SN 2008bk (M(56Ni) = 0.0086M⊙), both at 547 d. The fluxes are normalized to unity at 7200 Å in both panels. The observations are corrected for reddening and redshift.

9.2 Specific comparisons to a sample of Type II SNe

9.2.1 Comparison to low-luminosity SNe II 1997D, 2008bk, and 2016aqf

Figures 16 and 17 present a comparison of models s9p5, s11, and s10 with the observed optical spectra of low-luminosity SNe II 1997D, 2008bk (the model was computed for the same epoch as the observations, that is, 547 days after the inferred time of explosion), and 2016aqf. Unless specifically stated, all models correspond to a SN age of 350 days. Rather than give a magnitude offset between models and observations, we just state the model ⁵⁶Ni mass and that inferred for the observed SN. The two spectra shown are normalized at some wavelength (chosen to improve the visibility and facilitate the comparison; this wavelength is often taken at 6800 Å, away from strong lines, in a region dominated by a smooth background flux dominated by Fe I emission).

These low-energy explosions in lower-mass progenitors (E_kin in the range 2–6 × 10⁵⁰ erg) provide a suitable match to the observations, which are characterized by relatively narrow lines, strong Hα, weak metal lines. The model s9p5 yields a satisfactory match to the observations of SN 1997D, although its ⁵⁶Ni mass is ten times greater. Alternatively, we could have used model s9, whose ⁵⁶Ni mass is only twice that inferred for SN 1997D. Model s9 yields a very good fit (not shown) to most of the spectral features of SN 1997D, but strongly overestimates the [Ca II] λλ 7291, 7323 because Ca is mostlyCa⁺ in the H-rich layers (see Table 2). It is unclear what conditions lead to a modest [Ca II] λλ 7291, 7323 in SN 1997D. Ca over-ionization is not expected for a low ⁵⁶Ni mass. A lower metallicity might explain in part this property. This requires further study.

The model s11 at 350 and 547 days does a poor job at matching the [Ca II] λλ 7291, 7323 of SN2016aqf and SN 2008bk, but the rest of the spectra are reasonably well fitted. Model s11 overestimates the strength of [O I] λλ 6300, 6364 and [Ca II] λλ 7291, 7323 compared to SN 2016aqf perhaps because this model predicts a significant contribution from the H-rich material. This contribution may therefore not hold (this extra flux at larger velocity also leads to a mismatch in line width, though this could be cured by using a weaker explosion). The emission lines in the observed spectrum of SN 2016aqf (Müller-Bravo et al. 2020) also seem redshifted relative to the model, which may point to a large-scale asymmetry of the ⁵⁶Ni in that ejecta.

The association of low or moderate mass massive stars exploding as weak explosions with low-luminosity SNe II-P confirms the results from previous studies on SN 1997D (Chugai & Utrobin 2000), SN 2008bk (Maguire et al. 2012; Lisakov et al. 2017), or SN 2016aqf (Müller-Bravo et al. 2020). The S16 model s9 was also studied in detail by Jerkstrand et al. (2018) and found to compare favorably with the observations of low-luminosity SNe 1997D, 2005cs, or 2008bk. This result seems robust, is corroborated by photospheric phase properties, and appears to hold for the whole class of such events (Lisakov et al. 2018).

Fig. 17 Same as Fig. 16, but now for a comparison of model s10 (M(56Ni) = 0.025M⊙) and s11 (M(56Ni) = 0.016M⊙) to the low-luminosity SN 2016aqf (M(56Ni) = 0.008M⊙). The fluxes are normalized to unity at 6800 Å (left panel) and 5400 Å (right panel).

9.2.2 Comparison to standard-luminosity SNe 2012aw and 2013ej

Figure 18 compares the optical spectra of SN 2012aw and model s12. The model matches the Fe II emission forest, the strong O I, H I, and Ca II lines well. The line profiles are also well matched in width, although the model emission has structure not obviously present in the observations. The model Hα lacks a central narrow peak, likely because we do not extend our 1D ejecta models to sufficiently low velocities (the absence of material at low velocity is an artifact of performing the explosion in spherical symmetry). The model [O I] λλ 6300, 6364 doublet line has instead narrow peaks at the rest wavelength of each component, while the observed profile is smoother. This arises because the O-emitting material is more randomly distributed in space than adopted in our simplified shuffled-shell structure. The representative velocity of the O I-emitting volume, as estimated from the full width at half maximum, is, however, well matched. Probably for the analogous reason, the Fe I emission lines in the 8000 Å region are also too resolved in the model, something that could be remedied by a more extended volume of emission from Fe I (this Fe I emission comes primarily from the original ⁵⁶Ni -rich material, which is characterized by a lower Fe ionization level compared to the H-rich material in which Fe is once ionized; see Table 2). This is an obvious, though not critical, shortcoming of using a shuffled-shell structure in place of a more realistic and complex 3D distribution of ⁵⁶Ni.

Figure 19 compares model s12p5 with SN 2013ej. The particularity of this model is its high V_m, the largest of our model grid, combined with its modest metal yields, typical of a low-mass massive star. This is conducive to a lower-density ejecta and faster expansion in the emitting region, favoring the formation of forbidden transitions and the production of broad emission lines. These characteristics seem to match the properties of the optical spectrum of SN 2013ej well. For example, the Fe II emission below 5500 Å is quite strong relative to the strongest emission lines, which is better explained by a lower-mass progenitor such as s12p5. The [Ca II] λλ 7291, 7323 lines width is underestimated (but other lines seem well matched in width), which might indicate that some Ca II emission arises from H-rich material at large velocities (something that is prevented in most of our models because Ca is twice ionized in such regions).

Previous studies of SNe 2012aw and 2013ej have yielded inconsistent estimates for the ZAMS mass of their progenitors. Nebular phase modeling of 2012aw by Jerkstrand et al. (2014) suggest a ZAMS mass of 15 M_⊙ for the progenitor, while Yuan et al. (2016) suggest 12–15 M_⊙ for 2013ej. Radiation-hydrodynamics modeling of the bolometric light-curve combined with constraints on the photospheric velocity inferred from Doppler-broadened line profiles have produced higher mass estimates for SN 2012aw. Dall’Ora et al. (2014) infer a 20 M_⊙ ejecta, pointing to a massive main sequence star of perhaps 25 M_⊙. Similarly, Morozova et al. (2018) constrain a ZAMS mass of 20 M_⊙ and Nikiforova et al. (2021) a ZAMS mass of 25 M_⊙. Using non-LTE radiative transfer modeling of the multiband light curves and optical spectra of SNe 2012aw and 2013ej, Hillier & Dessart (2019) find that a 15 M_⊙ progenitor (characterized by different pre-SN mass loss rates) with a standard explosion energy can yield a satisfactory match to observations. Given the degeneracy of light curve modeling (Dessart & Hillier 2019; Goldberg et al. 2019), it is not surprising that there isa nonunique solution to the progenitor mass with this approach. The same scatter in inferred progenitor mass from light curve modeling holds for both SN 2012aw and SN 2013ej.

Fig. 18 Same as Fig. 16, but now for a comparison of optical spectra between model s12 (M(56Ni) = 0.032M⊙) and SN 2012aw (M(56Ni) = 0.074M⊙). The fluxes are normalized to unity at 6800 Å.

Fig. 19 Same as Fig. 16, but now for a comparison of model s12p5 (M(56Ni) = 0.05M⊙) to SN 2013ej (M(56Ni) = 0.06M⊙). The flux normalization is done at 5100 Å.

9.2.3 Comparison to the Type II-pec SN 1987A

In Fig. 20, we compare models s12A and s15p2 to the Type II-peculiar SN 1987A. Although both models offer a reasonable match they also have specific discrepancies. In model s12A, the [O I] λλ 6300, 6364, Hα, and [Ca II] λλ 7291, 7323 lines are well matched but the Fe II flux below 5500 Å is overestimated. This probably arises from the greater fraction of decay power absorbed bythe H-rich material relative to the O-rich and Si-rich material (which also stems from the greater V_m). The Ca II triplet is underestimated by at least a factor of three, perhaps somewhat surprising given that our progenitors have solar metallicity. In contrast, in model s15p2, the Fe II flux is weaker and better matched, but then the strength of [O I] λλ 6300, 6364 is overestimated. These subtleties highlight the difficulty of inferring accurately the ejecta properties and constraining the progenitor mass.

This comparison is not fully consistent since our models are at solar, rather than LMC metallicity, and were RSGs (rather than blue supergiants) when they exploded. The actual progenitor of SN 1987A may have had a mass closer to 16–20 M_⊙ (see, for example,Jerkstrand et al. 2011; S16), depending on rotation. We do not think the resulting small difference in core structure matters for the spectrum at 1 year. A more important shortcoming in this comparison is the neglect of clumping and CO molecular cooling.

9.2.4 Comparison to the low-metallicity SN 2015bs

Figure 21 presents the spectral comparison for models s25A and s25p2 at 350 days with the observations of SN 2015bs at 421 days after explosion. SN 2015bs is remarkable in two ways, first because of its location in a very low metallicity environment of about a tenth solar, and second because of its exceptional nebular spectrum exhibiting the strongest [O I] λλ 6300, 6364 line flux of all SNe II for which such data exist (Anderson et al. 2018). The later can only be matched by a higher-than-standard progenitor mass. Here, we show that model s25A reproduces satisfactorily the entire optical spectrum of SN 2015bs (although the spectrum is noisy and prevents a proper comparison of weak lines). As expected, the model s25p2, characterized by Si–O shell merging prior to explosion, exhibits a huge [Ca II] λλ 7291, 7323 flux, in conflict with observations. We thus confirm the previous conclusions reached with more simplistic models (Dessart & Hillier 2020a), as well as the previous work from Anderson et al. (2018).

10 Conclusions

We have presented non-LTE radiative transfer calculations based on state-of-the-art explosion models from WH07 and S16. Rather than trying to adequately adjust the models to match a specific observation, we simply compared the WH07 and S16 models to observations, without any scaling of their yields, ejecta mass, or ejecta kinetic energy. By retaining the global properties of the original models, our work provides the nebular-phase observables characterizing the ab initio explosion models. Although this approach lacks flexibility, it allowed us to construct a comprehensive grid of nebular phase spectra for the explosion of stars with an initial mass between 9 and 29 M_⊙.

Mixing was treated with a shuffled-shell technique (Dessart & Hillier 2020b), with the ansatz that the metal-rich core is fully mixed macroscopically and that 1–2M_⊙ of material from the H-rich envelope is mixed into the metal rich core. While this 1D shuffling implies an adjustment to the initial ejecta models of WH07 and S16, it treats the chemical mixing that must somehow be accounted for to reflect the 3D mixing witnessed inobservations (Abellán et al. 2017) and systematically present in 3D simulations of the neutrino-driven explosion mechanism of massive stars (Gabler et al. 2021).

The general trend of increasing explosion energy, metal yields, and ⁵⁶Ni mass with increasing progenitor mass leads to nebular phase spectra of increasing brightness, line width, and line strength. The low-energy explosion models from 9–11 M_⊙ stars match some of the well-known low-luminosity SNe II-P and confirm the widespread notion that these faint SNe II-P arise from the lowest-mass massive stars undergoing Fe-core collapse. Their low metal yield is an intrinsic characteristic of both the progenitor and the ejected material rather than the result of fall back in a high-mass progenitor. This results in weak and narrow [O I] λλ 6300, 6364, but a very strong narrow Hα testifying for the large fraction of decay power absorbed by H-rich material. Allowing for binarity,the progenitors of low-luminosity SNe II-P are expected to extend below 9 M_⊙ (Zapartas et al. 2017). Stars with an initial mass above about 12 M_⊙ that produce an ejecta kinetic energy of about 10⁵¹ erg are broadly compatible with standard-energy SNe II-P like SN 2012aw or 1987A. The ab initio explosion models of S16 therefore depict a landscape of progenitors and explosions that are broadly consistent with observations of SNe II-P.

However, within this general picture, our results reveal diversity and scatter, even for models that differ little in progenitor mass. This scatter takes its origin from the somewhat chaotic pattern of kinetic energy, ⁵⁶Ni mass, and envelope composition exhibited by the S16 models. In the present grid of models, one ingredient that impacts drastically the spectral properties is the Ca ionization in the H-rich material. In general, Ca²⁺ dominates and this relatively high ionization quenches the [Ca II] λλ 7291, 7323 line emission from the H-rich layers. In that case, [Ca II] λλ 7291, 7323 arises primarily from the Si-rich layers. However, Ca⁺ may dominate in H-rich layers when both the explosion energy and the ⁵⁶Ni mass are low (as in model s9 from S16), causing very strong [Ca II] λλ 7291, 7323 line emission. Observations suggest that this is rarely the case, and hence indicate that [Ca II] λλ 7291, 7323 emission is generally impaired by Ca overionization. Although not explored in this work, a low primordial metallicity could contribute to inhibit metal-line emission from the H-rich envelope.

A second ingredient is Si–O shell merging in the progenitor prior to explosion which may boost the Ca abundance in the O-rich shell by a factor of 100. [Ca II] λλ 7291, 7323 then becomes an important coolant for the O-rich material, weakening [O I] λλ 6300, 6364 emission and strengthening [Ca II] λλ 7291, 7323 emission. The magnitude of the effect, seen in models s20p1 and 25p2, is very strong, and typically incompatible with observations. However, one cannot exclude that it occurs in a milder form, thus allowing a higher mass progenitor to mimic the spectral appearance of lower-mass progenitor stars, which are characterized by a weaker [O I] λλ 6300, 6364.

Most of the erratic behavior concerns [Ca II] λλ 7291, 7323 emission. In contrast,the flux associated with Fe II emission shortward of 5500 Å, the [O I] λλ 6300, 6364 or the Hα line flux exhibit a more consistent evolution with progenitor mass. Excluding the lower-mass progenitors (which are easily identified – see above), models exhibit a strengthening [O I] λλ 6300, 6364 and a weakening Hα with increasing progenitor mass, reflecting the corresponding trend of increasing (decreasing) O (H) mass. While these evolutions are robust, there are not generally strong enough to separate two models unless they differ in progenitor mass by several M_⊙. The most vivid example of this evolution is the extreme case of model s29A, our highest-mass progenitor with a light residual H-rich envelope.

In general, our models produce a [O I] λλ 6300, 6364 line flux that is stronger than previously obtained with SUMO (Jerkstrand et al. 2015b), as already reported in Dessart & Hillier (2020a). Generally we can explain the nebular-phase spectra of observed standard-energy SNe II with a progenitor in the mass range 12 to 15 M_⊙. This offset might indicate that we are missing a coolant from the O-rich zone. One possibility is molecularcooling, which is presently ignored. Clumping is also found to reduce the [O I] λλ 6300, 6364 line flux. Another possibility is the slight contamination of Ca into the O-rich shell, which only occurs (but then strongly) in two of the higher mass progenitor models. A combination (and perhaps all) of these effects probably act at some level in Type II SN ejecta.

We continue to find that SN 2015bs is the current best evidence that higher mass RSG stars of 20–25 M_⊙ can explode. One could wonder whether SiO and CO cooling could affect this inference. This is unlikely because the O/Si and C/O shells, where these molecules may form, represent a smaller fraction of the O-rich shell mass in higher mass progenitors: a growing fraction of the decay power will thus be absorbed by the O/Ne/Mg shell, in which such molecules do not form (see bottom panel of Fig. 2 for the O/C shell mass in our model set).

Metal line emission from the H-rich material suggests that the [O I] λλ 6300, 6364 and [Ca II] λλ 7291, 7323 lines could be stronger (weaker) at higher (lower) metallicity, in particular in lower-mass progenitors because they are characterized by low metal yields. This applies in particular to Fe I, Fe II, and Ca II line emission.

Our model grid suggests that at 350 days, about 70% of the total decay power absorbed is radiated between 3500 and 9500 Å. The fraction of the decay power emitted that escapes the ejecta varies between models, being a few percent at most in weaker explosion and at most 35% in stronger explosions.

To improve the consistency of our ejecta models, we have explored the impact of the ⁵⁶Ni -bubble effect and the influence of clumping on SN radiation. For the present set of Type II SN ejecta models, we find that these processes play a modest role, probably because the ejecta regions concerned by these processes have a modest ionization, or because the effect is to weak. A significant improvement for future work will be to use 3D explosion models as initial conditions for our radiative-transfer calculations, such as those of Gabler et al. (2021).

A major drawback of the current work is the assumption of spherical symmetry for the progenitor evolution, explosion phase, and the radiative transfer. For the first, little can be done. However, 3D explosion models exist and can provide critical insights for radiative transfer. Furthermore, it is possible with CMFGEN to address multi-D effects by performing 1D simulations along multiple radial directions of a 3D explosion model. We can also post-process such simulations with 2D polarized radiation transfer to study the influence of asymmetry on line profile morphology and polarization (see, for example, Dessart et al. 2021).

Fig. 20 Same as Fig. 16, but now for a comparison of model s12A (M(56Ni) = 0.05M⊙) and s15p2 (M(56Ni) = 0.063M⊙) to SN 1987A (M(56Ni) = 0.069M⊙). Both models adopt a solar metallicity, not exactly consistent with the LMC metallicity for SN 1987A. The lower-mass progenitor does a good job at matching the three strongest lines, but severely underestimates the Ca II triplet, and does a poor job with at matching the Fe emission forest shortward of 6000 Å. By contrast, the higher mass model does a better job with the Fe emission, overestimates the O I doublet, and underestimates the Ca II line strengths. The flux normalization is done at 6800 Å.

Fig. 21 Same as Fig. 16, but now for the comparison between model s25A (M(56Ni) = 0.157M⊙; left) and s25p2 (M(56Ni) = 0.112M⊙; right) and observations of SN 2015bs (M(56Ni) = 0.048M⊙). Both models adopt a solar metallicity, not exactly consistent with the one tenth solar metallicity inferred for SN 2015bs (Anderson et al. 2018). The flux normalization is done at the Hα rest wavelength in both panels. Some gaussian smoothing has been applied to the observed spectrum to reduce the level of noise.

Acknowledgements

This work was supported by the “Programme National de Physique Stellaire” of CNRS/INSU co-funded by CEA and CNES. Hillier thanks NASA for partial support through the astrophysical theory grant 80NSSC20K0524. This work was granted access to the HPC resources of CINES under the allocation 2019 – A0070410554 and 2020 – A0090410554 made by GENCI, France. This research has made use of NASA’s Astrophysics Data System Bibliographic Services. H.T.J. acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through Sonderforschungsbereich (Collaborative Research Centre) SFB-1258 “Neutrinos and Dark Matter in Astro- and Particle Physics (NDM)” and under Germany’s Excellence Strategy through Cluster of Excellence ORIGINS (EXC-2094)-390783311.

Appendix A Results for S16 models with a scaled kinetic energy

Figure A.1 presents results for the measurements of various powers and line luminosities for a larger set of models. Besides the selected 23 models from WH07 and S16, we include 18 model counterparts from the S16 set in which the ejecta kinetic energy was scaled up or down by a factor of two. In practice, models with a lower (higher) kinetic energy relative to the models of similar initial mass were scaled in velocity by a factor $\sqrt{2}$ ($1/\sqrt{2}$). The radius was also adjusted to retain the same SN age, and the density was scaled so that R³ ρ remains unchanged. Such scaled versions are artificial. For example, we retain the same chemical composition from explosive nucleosynthesis although the explosion energy was modified. Nonetheless, they offer a means to gauge the sensitivity of spectral signatures to controlled variations in ejecta properties.

Fig. A.1 Same as Fig. 5, but now for a larger set of models including the 23 models discussed in the main text, as wellas the S16 models that were scaled in kinetic energy by a factor of 2 or 0.5 (filled red stars, models labeled S16x).

The factor of two variation in kinetic energy that was artificially introduced leads to small variations in the results, as was anticipated from the earlier discussion in Sect. 3.1 (see, for example, the spectral differences between model s12 and s12A in Fig. 4).

Appendix B Comparison of models s20A and s20p1

Figure B.1 compares the chemical composition of the inner ejecta in the models s20A from WH07 and s20p1 from S16. The latter underwent a Si–O shell merging in the last stages of evolution prior to core collapse, while the former was not affected by this phenomenon. The corresponding spectral properties are discussed in the main text, and shown in the middle panel of Fig. 4. With Si–O shell merging, the [Ca II] λλ 7291, 7323 doublet becomes a strong coolant for the O-rich material, which quenches the [O I] λλ 6300, 6364 line flux.

Fig. B.1 Chemical composition versus Lagrangian mass in the s20A model (left) of WH07 and the s20p1 model (right) of S16 at 350 days (the original mass fraction of 56Ni is shown together with that of Fe). The origin of the x-axis is at the neutron-star center in the left panel, or at the neutron star surface (at a Lagrangian mass of 1.82 M⊙) in the right panel. Si–O shell merging is present in the s20p1 model, but absent in the s20A model.

References

All Tables

All Figures

Fig. 1 Illustration of the composition profile in the original model and after shuffling individual shells. Left: chemical composition versus Lagrangian mass in the s9 and s15p2 ejecta models of S16 at 350 days (the original mass fraction is shown for 56Ni) – ejecta properties are summarized in Table 1. The name and total yield of each selected species are labeled.An alternating white and gray background is used to distinguish between consecutive shells of distinct composition. Right: corresponding ejecta in which shells have been shuffled in mass space within the Lagrangian mass cut Msh = 2.0 and 5.0 M⊙. Other ejecta models exhibit similar profiles, only modulated by the original composition (in particular the mass of each dominant shell) and the choice of Msh (details on the shuffling procedure are given in the appendix of Dessart & Hillier 2020b).

In the text

Fig. 2 From top to bottom: ejecta kinetic energy, the mean expansion rate, the 56Ni mass, the mass of the H-rich shell (i.e., corresponding to the progenitor H-rich envelope), and finally the masses of the O-rich shell (which includes the O/Si, the O/Ne/Mg, and the O/C shells) and of the O/C shell in our selection of WH07 (circles) and S16 (stars) models. The Ekin is prescribed in the WH07 models, but computed ab-initio in the S16 models. This explains the more complicated trends in the latter (e.g., the systematically lower Ekin and M(56Ni) for masses above about 16 M⊙, with some exceptions with the Ekin of s25p2), which includes regions in original mass that lead to no explosion at all (hence these models are not included in our sample; see S16 for a discussion).

In the text

Fig. 3 Montage of spectra for the set of WH07 models (blue label) and S16 models (red label) at 350 days after explosion. The wavelength range is limited to the 5600–7500 Å region to reveal the variations in [O I] λλ 6300, 6364, Hα, and [Ca II] λλ 7291, 7323. We note the non-monotonic evolution of these line fluxes with increasing initial mass in the S16 models, in contrast to the results obtained for the WH07 models.

In the text

Fig. 4 Comparison between the S16 (blue) and WH07 (red) models of comparable main sequence mass, with s12 and s12A (top), s20p1 and s20A (middle), and s25p2 and s25A (bottom). In each panel, the flux is divided by the integral of the flux falling between 1000 and 25 000 Å (and then scaled by a factor of 103), which is similar to forcing all models to have the same total decay power absorbed (in reality, the optical flux differs between models of the same mass because of the differences in 56Ni or γ-ray escape; see Table 1). This normalization procedure also allows one to estimate by eye the cooling power of strong lines relative to the total decay power absorbed.

In the text

	Fig. 5 Line fluxes, powers, and their ratios for our grid of radiative-transfer calculations based on the WH07 (circles) and S16 (stars) ejecta models. Each panel in discussed in turn in Sect. 3.2. A similar figure, including S16 models scaled in kinetic energy, is presented in Fig. A.1.
In the text

	Fig. 6 Profiles for the composition of dominant species, the cumulative decay power absorbed, and the cumulative frequency-integrated flux versus depth. For the latter two, the sum is performed inward from the outermost ejecta location and normalized for visibility.
In the text

Fig. 7 Illustration of the spatial regions (here shown in velocity space) contributing to the emergent flux in the S16 model s15p2. The grayscale image shows the observer’s frame flux contribution ∂Fλ,V ∕∂V (the map maximum is saturated at 20% of the true maximum to bias against the strong emission lines and better reveal the origin of the weaker emission) versus wavelength and ejecta velocity. The plot at left connects these emission contributions to ejecta shells and the normalized decay-power deposition profile (dashed line) on the same velocity scale. The upper panel shows the corresponding normalized flux Fλ integrated over all ejecta velocities.

In the text

Fig. 8 Top: illustration of ejecta and radiation properties for the shuffled-shell model s15p2. From upper to lower panels, we show the profiles versus velocity for the mass fraction of H, He, O, Ca, and Fe, their ionization state (zero for neutral, one for once ionized etc.; see also Table 2), and the formation regions for Na I D, [O I] λ 6364, Hα, He I 7065 Å, and [Ca II] λ 7291. The quantity ∫ ξ d log V is the line equivalent width, which here is normalized to unity for visibility.

In the text

Fig. 9 Illustration of the main cooling processes balancing the radioactive decay heating at all ejecta depths in the model s15p2. We show each dominant cooling rate (stepping down from the rate having the largest peak value at any depth) normalized to the local heating rate. The term “NT” stands for nonthermal processes (in this context nonthermal excitation and ionization) and “COL” stands for collisional processes (i.e., collisional excitation).

In the text

	Fig. 10 Same as Fig. 8, but now for the s9 model, which corresponds to the lowest-mass progenitor model in the S16 set.
In the text

	Fig. 11 Illustration of the influence of the bubble effect applied to the reference model s15A. We show the resultsfor the SN radiation and gas properties in model s15A and for its counterparts in which the 56Ni -rich material occupies 30 (model bubble0p3) and 60% (model bubble0p6) of the volume within Msh (see Sect. 5 for a discussion).
In the text

	Fig. 12 Illustration of the influence of clumping, comparing the reference model s15A and its counterparts s15Afvol0p3 and s15Afvol0p1 in which a uniform volume filling factor of 0.3 and 0.1 is used.
In the text

	Fig. 13 Comparison of model s9p5 with the standard mixing procedure produced by the shuffled-shell technique (see Sect. 2) and the same model in which the original unmixed ejecta has been used without adjustments (model s9p5nomix).
In the text

Fig. 14 Montage of spectra in the [O I] λλ 6300, 6364 region (top), the Hα region (middle), and the [Ca II] λλ 7291, 7323 region (bottom) for the s9 and the s15p2 models of S16. A thin vertical line indicates the rest wavelength of each component of the doublet linesand Hα. Colored lines indicate the flux arising exclusively from bound-bound transitions associated with a specific ion (other ions and their opacity are then neglected, which explains the flux difference and the excess flux in places). Notice the presence of [Ni II] 7377.8 and [Ni II] 7411.6 Å in the bottom panels.

In the text

	Fig. 15 Flux measurements for the Fe II forest between 4100 and 5500 Å, Hα, [O I] λλ 6300, 6364 and [Ca II] λλ 7291, 7323, for our models (shown versus initial mass; symbols) and observations (colored bars at right). In all cases, we normalize to the total flux falling between 4000 and 9000 Å.
In the text

	Fig. 16 Left: comparison of model s9p5 (M(56Ni) = 0.017M⊙) at 350 d withthe observations of SN 1997D (M(56Ni) = 0.002M⊙) at 350 d. Right: comparison of model s11 (M(56Ni) = 0.016M⊙) with SN 2008bk (M(56Ni) = 0.0086M⊙), both at 547 d. The fluxes are normalized to unity at 7200 Å in both panels. The observations are corrected for reddening and redshift.
In the text

	Fig. 17 Same as Fig. 16, but now for a comparison of model s10 (M(56Ni) = 0.025M⊙) and s11 (M(56Ni) = 0.016M⊙) to the low-luminosity SN 2016aqf (M(56Ni) = 0.008M⊙). The fluxes are normalized to unity at 6800 Å (left panel) and 5400 Å (right panel).
In the text

	Fig. 18 Same as Fig. 16, but now for a comparison of optical spectra between model s12 (M(56Ni) = 0.032M⊙) and SN 2012aw (M(56Ni) = 0.074M⊙). The fluxes are normalized to unity at 6800 Å.
In the text

	Fig. 19 Same as Fig. 16, but now for a comparison of model s12p5 (M(56Ni) = 0.05M⊙) to SN 2013ej (M(56Ni) = 0.06M⊙). The flux normalization is done at 5100 Å.
In the text

Fig. 20 Same as Fig. 16, but now for a comparison of model s12A (M(56Ni) = 0.05M⊙) and s15p2 (M(56Ni) = 0.063M⊙) to SN 1987A (M(56Ni) = 0.069M⊙). Both models adopt a solar metallicity, not exactly consistent with the LMC metallicity for SN 1987A. The lower-mass progenitor does a good job at matching the three strongest lines, but severely underestimates the Ca II triplet, and does a poor job with at matching the Fe emission forest shortward of 6000 Å. By contrast, the higher mass model does a better job with the Fe emission, overestimates the O I doublet, and underestimates the Ca II line strengths. The flux normalization is done at 6800 Å.

In the text

Fig. 21 Same as Fig. 16, but now for the comparison between model s25A (M(56Ni) = 0.157M⊙; left) and s25p2 (M(56Ni) = 0.112M⊙; right) and observations of SN 2015bs (M(56Ni) = 0.048M⊙). Both models adopt a solar metallicity, not exactly consistent with the one tenth solar metallicity inferred for SN 2015bs (Anderson et al. 2018). The flux normalization is done at the Hα rest wavelength in both panels. Some gaussian smoothing has been applied to the observed spectrum to reduce the level of noise.

In the text

	Fig. A.1 Same as Fig. 5, but now for a larger set of models including the 23 models discussed in the main text, as wellas the S16 models that were scaled in kinetic energy by a factor of 2 or 0.5 (filled red stars, models labeled S16x).
In the text

Fig. B.1 Chemical composition versus Lagrangian mass in the s20A model (left) of WH07 and the s20p1 model (right) of S16 at 350 days (the original mass fraction of 56Ni is shown together with that of Fe). The origin of the x-axis is at the neutron-star center in the left panel, or at the neutron star surface (at a Lagrangian mass of 1.82 M⊙) in the right panel. Si–O shell merging is present in the s20p1 model, but absent in the s20A model.

In the text