The explosion of 9$-$29$M_\odot$ stars as Type II supernovae : results from radiative-transfer modeling at one year after explosion

We present a set of nonlocal thermodynamic equilibrium steady-state calculations of radiative transfer for one-year old type II supernovae (SNe) starting from state-of-the-art explosion models computed with detailed nucleosynthesis. This grid covers single-star progenitors with initial masses between 9 and 29$M_{\odot}$, all evolved with KEPLER at solar metallicity and ignoring rotation. The [OI]$\lambda\lambda$$6300,6364$ line flux generally grows with progenitor mass, and H$\alpha$ exhibits an equally strong and opposite trend. The [CaII]$\lambda\lambda$$7291,\,7323$ strength increases at low $^{56}$Ni mass, low explosion energy, or with clumping. This CaII doublet, which forms primarily in the explosively-produced Si/S zones, depends little on the progenitor mass, but may strengthen if Ca$^+$ dominates in the H-rich emitting zones or if Ca is abundant in the O-rich zones. Indeed, Si-O shell merging prior to core collapse may boost the CaII doublet at the expense of the OI doublet, and may thus mimic the metal line strengths of a lower mass progenitor. We find that the $^{56}$Ni bubble effect has a weak impact, probably because it is too weak to induce much of an ionization shift in the various emitting zones. Our simulations compare favorably to observed SNe II, including SN2008bk (e.g., 9$M_{\odot}$ model), SN2012aw (12$M_{\odot}$ model), SN1987A (15$M_{\odot}$ model), or SN2015bs (25$M_{\odot}$ model with no Si-O shell merging). SNe II with narrow lines and a low $^{56}$Ni mass are well matched by the weak explosion of 9$-$11$M_{\odot}$ progenitors. The nebular-phase spectra of standard SNe II can be explained with progenitors in the mass range 12$-$15$M_{\odot}$, with one notable exception for SN2015bs. In the intermediate mass range, these mass estimates may increase by a few $M_{\odot}$ with allowance for clumping of the O-rich material or CO molecular cooling.


Introduction
The late time radiative properties of core-collapse supernovae (SNe) contain information on the progenitor star at 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 56 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 (say at one year after explosion). At that time, the velocity is essentially proportional to 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, that was limited to the s15A model of Woosley & Heger (2007, WH07), to study 23 explosion models taken from WH07 and Sukhbold et al. (2016, 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 51 erg at infinity. In contrast, S16 simulates the evolution of the collapsing star through core bounce, shock formation, the stalling of the shock, and 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 from S16, while the WH07 models are used to test the sensitivity of results to the numerical handling of the explosion.
Our work provides a complementary picture to the results obtained, for the most part, with the code SUMO (Jerkstrand et al. 2011, 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 macro-Article number, page 1 of 28 arXiv:2105.13029v1 [astro-ph.SR] 27 May 2021 A&A proofs: manuscript no. ms scopically mixed progenitor metal-rich core. Such studies of nearby Type II SNe suggest a lack of high-mass RSG progenitors Maguire et al. 2012;Jerkstrand et al. 2014Jerkstrand et al. , 2015bYuan et al. 2016). These studies were all targeted towards 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 radiativetransfer modeling with CMFGEN 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 Type II SN simulations and compares our results with those from previous work. The treatment of the 56 Ni-bubble effect, presented in Section 5, reveals only a minor impact on the SN radiative properties. Similarly, tests for material clumping, presented in Section 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 impact significantly the nebular spectral properties (Section 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 Section 9. We give our concluding remarks in Section 10.

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 12 C(α, γ) 16 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(Rauscher et al. ) -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 51 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 in-terpolated 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 highdensity PNS core. Neutrino transport in the accretion mantle of the PNS is treated by a gray approximation. The neutrino-driven supernova models employed in the present study are spherically symmetric. However, the neutrino-driven mechanism is known to be a multi-dimensional 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. 2020 are replaced by an extended post-bounce period of PNS accretion and, after a considerably delayed onset of the explosion, a subsequent neutrinodriven wind, whose strength is overestimated in order to provide a neutrino-heated mass with an energy that is needed to power the supernova. Clearly, explosion asymmetries and large-scale mixing processes, which are a natural consequence of 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. 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 56 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 51 erg and where the mass cut was placed at the location in the progenitor where the dimensionless entropy S /N A k B = 4. 1 Simulations for other masses were  Table 1. The name and total yield of each selected species is labeled. An alternating white and grey 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 M sh = 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 M sh (details on the shuffling procedure are given in the appendix of Dessart & Hillier 2020b). 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 non-integer 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 56 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 Fig. 2. Both WH07 and S16 models show the same trend of increasing O-rich shell mass (defined cline at the base of the oxygen shell and modern calculations frequently find a mass cut there (see discussion, for example, in Ertl et al. 2016). 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 lowermass 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 51 erg irrespective of ZAMS mass, with a monotonic increase of the 56 Ni yield with ZAMS mass that reflects the monotonic increase of 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 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 explo- Table 1. Ejecta model properties from WH07 (upper part) and S16 (lower part; the selected models were exploded with the W18 engine), including the ejecta mass, ejecta kinetic energy, the cumulative yields of H, He,O,Mg,Si,Ca,and 56 Ni prior to decay. We specify whether the Si-rich and the O-rich shells merged before core collapse (this did not occur in the WH07 models, and only in three of the S16 models). The last two columns give the choice of mass cut M sh for the shell shuffling, as well as the associated velocity that bounds 99% of the total 56 Ni mass in the corresponding model.

Model
M sion. 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 earlier discussion in this section). In the S16 models, a trend towards higher E kin and 56 Ni yield with increasing ZAMS mass can be discerned. There is a modest scatter for the 56 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 51 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 56 Ni mass.
Taking the selected WH07 and S16 models at 350 d 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 d. 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 56 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 56 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  . From top to bottom, this figures shows the ejecta kinetic energy, the mean expansion rate, the 56 Ni 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 E kin 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 E kin and M( 56 Ni) for masses above about 16 M , with some exceptions with the E kin 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 discussion).
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 Hrich material mixed inwards 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 56 Ni two-step decay chain and ignore the contribution from other unstable isotopes. Non-local energy deposition is computed by solving the grey radiative transfer equation (the γ-ray grey absorption-only opacity is set to 0.06 Y e cm 2 g −1 , where Y e is the electron fraction). Non-thermal 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 −1 . 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 −1 are essentially the same as those obtained for 50 km s −1 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 d 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 51 erg, and a moderate range in 56 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 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 rea-   . 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 25000 Å (and then scaled by a factor of 10 3 ), 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 56 Ni 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. sons for this. First, the S16 model extends to lower-mass massive stars. S16 predicts a low explosion energy of 1 to 6 × 10 50 erg for initial masses 9 to 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 56 Ni masses of 0.004 to 0.03 M for these lowenergy 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 re-A&A proofs: manuscript no. ms . 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 Section 3.2. A similar figure, including S16 models scaled in kinetic energy, is presented in Fig combination 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 to s12 stands well apart from the rest of the sample shown in Fig. 3.
For the mass range 12 to 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 ten years apart and in a different fashion (see WH07 and S16 for 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. 2 Model s12A also has 60% more 56 Ni than model s12, and thus a 60% greater luminosity but after normalization, this offset is gone. This offset in 56 Ni mass is too small to drive a visible change in ionization. Spectral line widths at neb-ular 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 56 Ni to large velocities. The similar spectral properties at 350 d 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 degeneracies are 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 56 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 consequence of the merging of the Sirich 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 −5 to ∼ 10 −3 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). 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 accounts for 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).

Quantitative description
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., ± 10000 km s −1 for Hα, ± 8000 km s −1 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 Herich 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 56 Ni mass and the γ-ray escape fraction. The latter is function of the ejecta density profile (which depends strongly on E kin /M ej ) and the location of the three 56 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 sub-shells). Our models tend to show a decreasing trapping of γ-rays towards 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( 56 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 56 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 in the 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 line luminosities.
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 are well 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 power will 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 Hrich shell may be the same for a wide range of progenitor mass ) the 56 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 Orich 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 56 Ni mass, which fosters an ionization shift to Ca + , compared to Ca 2+ 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 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 56 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-rich ejecta.
With the exception of two outliers (those models with Si-O shell merging prior to core collapse), the fractional

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;. 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,  found that the [Ca ii] λλ 7291, 7323 must arise from Hrich 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) 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 non-thermal processes (in this context non-thermal excitation and ionization) and "COL" stands for collisional processes (i.e., collisional excitation). 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 (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 2+ 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. 9 where 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 non-thermal 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 Orich 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 56 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 50 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 Fig. 3 and Fig. 5).

56 Ni Bubble effect
Decay heating causes a time integrated energy injection that can approach the local specific kinetic energy of the 56 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 56 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 56 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 56 Ni, as well as on the density and expansion rate of the surrounding material. This process is therefore very non linear 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 & Audit 2019) as well as the nebular phase properties of core-collapse SNe 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 56 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 Luc Dessart et Fig. 10. Same as Fig. 8, but now for the s9 model, which corresponds to the lowest mass progenitor model in the S16 set. To conserve mass, the density of the 56 Ni bubble must be divided by a factor f Ni = (α Ni V 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 56 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 56 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 56 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 56 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 56 Ni-rich shell and the other shells within M sh , we proceed from the innermost zone outwards 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 56 Ni-rich zone of volume δV located between radii r l and r u on the original Article number, page 13 of 28 A&A proofs: manuscript no. ms radial grid, we build a new radial grid and set : where f is equal to f Ni in 56 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 56 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.36 M as a reference, we computed two other models to test the influence of the bubble effect adopting α Ni = 0.3 and 0.6 (the 56 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 56 Ni-rich regions resulting from the lower density. As said earlier, while the bubble effect leads to a strong reduction of the density in 56 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 56 Ni bubble effect should have a similarly weak effect in other models. In lower-mass models, the 56 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 56 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.

Influence of a uniform clumping
The 56 Ni-bubble effect described above is expected to be just one component contributing to the complex 3D structure of corecollapse supernova ejecta. Besides the general density variations discussed in a simplified manner in the preceding section, the 56 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). 3 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 2+ transitions modestly to Ca + ). The increased density in the clumped model favors colli- 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. 4 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 2+ still dominates and prevents cooling by [Ca ii] λλ 7291, 7323. The Hrich zones might be significantly clumped if mixed inwards 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. 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 56 Ni was the primary ingredient influencing spectral properties and in particular emission line fluxes. This arises because the location of 56 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 56 Ni mixing is enhanced, a greater fraction of the power is deposited further out, biasing  against the inner ejecta regions. If only 56 Ni were mixed while other shells remained unmixed, this could strongly impact how the ejecta cools and how the spectra appear.

Influence of mixing
In nature, one expects a significant mixing of the metal-rich core in core-collapse SNe (say below 2000 − 3000 km s −1 in a standard-energy standard-mass SN II), in combination with the mixing of 56 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 neutrinodriven 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 56 Ni fingers while the model W15 has a few big 56 Ni protrusions but shows very weak 56 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 Article number, page 15 of 28  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

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 2E kin /M ej ) covers the range from 1230 km s −1 (model s9) to 3790 km s −1 (model s12p5), thus only 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 −1 , 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 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 section 9). This is in contrast with SNe Ibc, which reveal systematic shifts (see, for example Taubenberger et al. 2009;Milisavljevic et al. 2010).

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 56 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 56 Ni mass for these observations is generally obtained by assuming full trapping and estimating the fraction of the SN luminosity that falls outside of the observed range. Combined with the uncertainty in SN distance and reddening, these 56 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 56 Ni mass varies by a factor of a few (see also Section 8 of Dessart & Hillier 2020a). Hence, an uncertainty or an offset in 56 Ni mass for both observations and models is of limited impact for the discussion that follows.

General comments
The SN luminosity at nebular times is set by the 56 Ni mass (modulo γ-ray escape). To cancel the offset in 56 Ni mass between our models and observations, we renormalize the spectra, so that they match at a given optical wavelength. As the 56 Ni mass can also influence the ionization and temperature of the gas we restrict our comparisons to cases where the offset in 56 Ni mass between model and observation is moderate. For standardluminosity SNe II, the WH07 and S16 models have a 56 Ni mass that is within a factor of two from that inferred for our selection of observations. The offset between model and observations is typically larger for low-energy SNe II, since our selected observations cover a range from 0.002 to ∼ 0.009 M , while the 56 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 56 Ni mass. The distribution of 56 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 inwards 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 section 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 section 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 spec-trum 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 section 3.2, the nonmonotonicity 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 56 Ni-rich layers. This makes sense since the H-rich shells are much more massive than the 56 Ni-rich layers (by a factor of one hundred 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 56 Ni rich (where it is partially neutral). The scatter in the fractional Fe and Hα fluxes for the models and observations is similar.

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.

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 much 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 reported this 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 stateof-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 Section 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 10%) 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 −5 M at day 285, while Liljegren et al. (2020) estimated a much larger mass of 4 × 10 −3 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 −6 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(Kotak et al. , 2006. Another potential source of error in our models is explosion asymmetries. An offset in 56 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 56 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.
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 d after the inferred time of explosion), and 2016aqf. Unless specifically stated, all models correspond to a SN age of 350 d. Rather than give a magnitude offset between models and observations, we just state the model 56 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 to 6 × 10 50 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 56 Ni mass is ten times greater. Alternatively, we could have used model s9, whose 56 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 mostly Ca + in the H-rich layers (see Table 2  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 56 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 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).

Comparison to standard-luminosity SNe 2012aw and 2013ej
Figure 18 compares the optical spectra of SN 2012aw and model s12. The model matches well the Fe ii emission forest, the strong O i, H i, and Ca ii lines. 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 low-enough 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 56 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 56 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 lowmass massive star. This is conducive to a lower-density ejecta and faster expansion in the emitting region, favoring the forma- tion of forbidden transitions and the production of broad emission lines. These characteristics seem to match well the properties of the optical spectrum of SN 2013ej. 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 like 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  (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 Goldberg et al. 2019), it is not surprising that there is a non-unique 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.

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 by the 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 red supergiants (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.

Comparison to the low-metallicity SN 2015bs
Figure 21 presents the spectral comparison for models s25A and s25p2 at 350 d with the observations of SN 2015bs at 421 d 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).

Conclusions
We have presented nonLTE 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  . 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. 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 nebularphase 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 − 2 M 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 in observations (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 56 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 51 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, 56 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 2+ 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 56 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 molecular cooling, 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 d, 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 56 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 former, 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. 5, but now for a larger set of models including the 23 models discussed in the main text, as well as the S16 models that were scaled in kinetic energy by a factor of 2 or 0.5 (filled red stars, models labelled S16x).
Article number, page 27 of 28 A&A proofs: manuscript no. ms