© L. Dessart and D. J. Hillier 2020

1 Introduction

Nebular-phase spectroscopy provides critical information on the properties of massive star explosions and type II supernovae (SNe). Originally cloaked by a massive, optically-thick ejecta, the inner metal-rich layers of the SN are revealed after about 100 d as the H-rich material fully recombines and becomes transparent in the continuum. During this phase forbidden line emission, following collisional excitation and non-thermal excitation and ionization, is the dominant cooling process balancing ⁵⁶Co decay heating¹. Such line emission conveys important information on the composition and the yields, the large-scale chemical mixing, the explosion geometry, and, ultimately, the progenitor identity (for a review, see Jerkstrand 2017). In type II SNe, this is also the time when the original ⁵⁶Ni mass can be nearly directly (in the sense that it does not require any modeling of the SN radiation) extracted from the inferred SN luminosity (converting SN brightness to luminosity requires an estimate of the reddening and distance, and of the flux falling outside of the observed spectral range).

Nebular-phase spectroscopic modeling started in earnest with SN 1987A, for which information has been gathered through continuous photometric and spectroscopic monitoring (for reviews, see for example Arnett et al. 1989 and McCray 1993). Numerous studies were focused on the nebular-phase radiation properties of SN 1987A (see, for example, Fransson & Chevalier 1987; Kozma & Fransson 1992; Li et al. 1993; Li & McCray 1992, 1993, 1995), and led to a refined understanding of the physical processes at play at these late times (for an extended discussion, see for example Fransson & Chevalier 1989). The main conclusions from these studies is that under the influence of ⁵⁶Co decay heating, the various metal-rich and H-rich shells of the ejecta cool primarily through forbidden line emission of neutral and once-ionized species. Chemical mixing is inferred, although limited to large-scale macroscopic mixing and absent at the microscopic level. More advanced calculations have been performed since, with a more accurate and richer description of the ejecta composition, the atomic processes involved, and the atomic data employed (Kozma & Fransson 1998; Jerkstrand et al. 2011). A more extensive analysis of SNe II-P and the physics relevant to nebular-phase modeling has been presented in a series of papers by Jerkstrand et al. (2014, 2012, 2015a). This work was also made possible by the acquisition of nebular-phase spectra for nearby SNe, although, even today, the published sample remains limited to a handful of objects (see, for example, the dataset presented by Silverman et al. 2017).

An important conclusion drawn from these works is that the [O I] λλ 6300, 6364 doublet flux, or its ratio with that of [Ca II] λλ 7291, 7323, can be used to constrain the progenitor mass. Applied to the existing dataset, nebular-phase modeling indicates initial progenitor masses below about 17 M_⊙, with a noticeable absence of 20− 25M_⊙ progenitors (e.g., Jerkstrand et al. 2015a; see also Smartt 2009). A notable exception is the type II SN 2015bs, for which a high mass progenitor of 20− 25M_⊙ seems a plausible explanation for its unprecedented large [O I] λλ 6300, 6364 doublet flux (Anderson et al. 2018).

With CMFGEN (Hillier & Dessart 2012), we have conducted numerous simulations for type II SNe (Dessart & Hillier 2011, 2019a; Dessart et al. 2013, 2018; Lisakov et al. 2017) but generally limited these to the photospheric phase. At nebular times, our models match quite closely the observations of standard type II SNe like 1999em (see, for example, Dessart et al. 2013 and Silverman et al. 2017). With the treatment of non-thermal processes, CMFGEN predicts a strong Hα line at nebular times (say 300 d after explosion), while that line is absent in previous CMFGEN simulations in which non-thermalprocesses were ignored (Dessart & Hillier 2011). We shied away from the nebular phase because of our inability, in 1D, to treat chemical segregation satisfactorily. Indeed, in CMFGEN, we simultaneously enforce macroscopic and microscopic mixing, which conflicts with the properties of mixing seen in multidimensional simulations of core-collapse SN explosions (see, for example, Wongwathanarat et al. 2015). This reluctance is, however, questionable. First of all, our models match quite closely most of the observed SNe II at nebular times. Secondly, it is possible to introducetwo different levels of mixing in our simulations. We may for example apply strong mixing of ⁵⁶Ni and daughter isotopes, but impose a very low level of mixing for all other species (in the current context of CMFGEN, this mixing is microscopic and macroscopic). This is not fully satisfactory in that we then microscopically mix ⁵⁶Ni, ⁵⁶Co, and ⁵⁶Fe throughout the ejecta, but of these three species, Fe dominates at nebular times of 300 d². This Fe is already present at the microscopic level with a mass fraction of at least ~ 0.001 throughout most of the ejecta – this floor value is the solar-metallicity value (see Sect. 2.4 for a discussion on the limitations of this statement). Hence, this approach is not so far from what takes place in nature. We believe there is a clear interest in presenting a grid of simulations for SNe II at nebular times and exploring the sensitivity of SN radiation to our various ejecta properties. Such a grid has never been published.

In the next section, we present the numerical setup used for the ejecta and for the radiative transfer modeling. In Sect. 3, we discuss the influence of the adopted Doppler width on the resulting ejecta and radiation properties. For most models in this paper, we adopt a low value of 2 km s⁻¹. In Sect. 4, we describe in detail the results for a SN II from a 20 M_⊙ star on the zero-age main sequence (ZAMS). This case, taken as a reference, is used to discuss the physics controlling nebular-phase spectra, the line formation process at nebular times, and to identify the predicted nebular lines in the optical and near-infrared ranges. Section 5 discusses the critical impact of the O/Ca ratio in the O-rich shell on the [O I] λλ 6300, 6364 doublet strength. Section 6 discusses the impact of the H-rich envelope mass on the nebular-phase properties, which is relevant for comparing the emission properties of SNe II arising from progenitors of lower and higher mass (say between a 12 and a ≥ 25 M_⊙ star). We then discuss the influence of the adopted ⁵⁶Ni mixing and of the ⁵⁶Ni mass on our nebular-phase spectra in Sects. 7 and 8. While all previous simulations were performed at a SN age of 300 d, Sect. 9 discusses the evolution of SN II properties from 150 to 500 d after explosion. We discuss the influence of the progenitor mass (for a fixed ⁵⁶Ni mass and ejecta kinetic energy) on the nebular-phase spectra in Sect. 10. Section 11 presents a succinct comparison of our crafted nebular models with the observations of a few well observed SNe II at about 300 d after explosion. Finally, we present our conclusions in Sect. 12.

2 Context and numerical setup

2.1 Properties of massive star progenitors at core collapse: results from a grid of MESA models

We performed stellar evolution simulations with MESA (Paxton et al. 2011, 2013, 2015, 2018) version 10108 for a 12, 15, 20, 25, 27, and 29 M_⊙ star on the ZAMS. All models were evolved until iron core formation and collapse (i.e., when the maximum infall velocity is 1000 km s⁻¹). We assumed a solar metallicity mixture with Z = 0.014 and no rotation. We used the approx21.net nuclear network. We used the default massive star parameters in MESA with the following exceptions. As in our previous works, we used a mixing length parameter of three (see Dessart et al. 2013), which leads to smaller RSG surface radii. This modification was also used by Paxton et al. (2018) for their simulations of SNe II-P progenitors. We modified the parameters overshoot_f0_above_nonburn_core to be 0.001 instead of 0.0005 (same for overshoot_f0_above_burn_h_core, overshoot_f0_above_burn_he_core and overshoot_f0_ above_burn_z_core). We also modified overshoot_f_ above_burn_z_core and related quantities (counterparts for “h” and “he”) to be 0.004. With these slight enhancements in overshoot, the stability of the code is improved during the advanced burning stages and the nucleosynthesis is somewhat boosted (in both respects, enhanced overshoot has a similar impact as introducing some rotation in the progenitor star on the ZAMS). We used the “Dutch” recipe for mass loss, with a scaling factor “du” that we varied between zero (no mass loss) and one (we use factors of 0, 0.3, 0.6, 0.8, and 1.0). This allows us to gauge the influence of mass loss on the properties of the star (and thus to address the large uncertainty in mass loss rates) and, in particular, the stellar core at death.

In about 50% of the models produced with these parameters, the Si-rich and the O-rich shells merged during Si burning, producing a single Si-rich and O-rich shell with a nearly uniform composition (all elements are microscopically mixed). This has a considerable impact on the appearance of a core-collapse SN at nebular times (Fransson & Chevalier 1989). Indeed, Ca then becomes abundant throughout the merged shell, and the [Ca II] λλ 7291, 7323 doublet becomes the primary coolant, inhibiting the cooling through other lines and in particular the [O I] λλ 6300, 6364 doublet (see Sect. 5). Ca is also produced during the explosion, together with Si and similar elements. However, during the explosion there is no microscopic mixing of the Si-rich material with the O-rich shell because the mixing is exclusively macroscopic.

In this study, we wanted to exclude such configurations to keep this complexity aside. The general wisdom so far has been to assume that these deep convective burning shells do not merge (see for example Fransson & Chevalier 1989, but see discussion in Sects. 5 and 12). Hence, to prevent the merging of the Si-rich shell and the O-rich shell during Si burning with MESA we set min_overshoot_q to 1 and mix_factor to 0 when the central ²⁸Si mass fraction first reaches 0.4. All simulations discussed in this section were produced in this manner and exhibit clearly distinct Si-rich and O-rich shells. The resulting model properties are given in Table 1, including initial and final masses, the main shell masses, and some ratios of mean mass fractions of important species within these shells. Multidimensional simulations of the last burning stages of massive stars prior to collapse are needed to determine the level of mixing, if any, of the Si-rich and O-rich shells (e.g., Meakin & Arnett 2007; Couch et al. 2015; Chatzopoulos et al. 2016; Müller et al. 2017; Yoshida et al. 2019). A recent 3D simulation finds violent merging of the O and Ne shells in a star with an initial mass of 18.88 M_⊙ (Yadav et al. 2020).

With the adopted variation in mass-loss rate, the mass of the H-rich shell (or progenitor envelope) covers a large range from ~ 3 up to ~ 17 M_⊙. Stellar winds in this mass range, metallicity, and adopted mass-loss rate recipes, do not peel the star all the way down to the He core so the He-rich shelland deeper metal-rich shells are not directly affected by mass loss. The mass of the He shell falls in a narrow range between about 1.3 to 3.0 M_⊙. The mass of the O-rich shell follows a much greater variation, growing from around 0.3 in the 12 M_⊙ progenitor and rising to about 6.8 M_⊙ in the 29 M_⊙ progenitor, more than 20 times greater. The Si-rich shell follows a similar trend but increases by only a factor of six between 12 and 29 M_⊙ progenitors (i.e., it goes from about 0.1 to 0.6 M_⊙). Excluding the artificial models that ignore mass loss, these shell masses do not vary much with the different scalings adopted (see third column in Table 1).

In the O-rich shell, O is about 10⁴ times more abundant than Ca, whose mass fraction in that shell is equal to the original metallicity in our MESA simulations (set to solar in this work; see, however, Sect. 2.4). If the Si-rich and O-rich shells were fully mixed, the same models would yield an O/Ca mass ratio in the range 60− 200, so typically 100 times smaller than the O/Ca in the unadulterated O-rich shell.

The Mg to O mass ratio in the O-rich shell is comparable to the Ca to Si mass ratio in the Si-rich shell and is equal to about 0.05−0.1. While the Mg to O mass ratio tends to decrease with main sequence mass the adopted mass-loss rate introduces a small scatter at a given progenitor mass. Finally, the N to He mass ratio in the He-rich shell is around 0.005, with a scatter of about 50% from low to high mass progenitors.

These are representative composition properties for our set of MESA simulations evolved with the network approx21.net, which is routinely used for massive star explosions. We do not attempt in this study to further adjust the composition to reflect additional nuclear processes not accounted for by the network approx21.net. Using state-of-the-art explosion models with a fully consistent composition computed with a huge network is straightforward but is delayed to our next study.

2.2 A simplified description of core-collapse supernova ejecta

Using a physical model of the explosion has the advantage of consistency. The ejecta structure and composition are determined by the laws that govern stellar evolution, stellar structure, and radiation hydrodynamics of stellar explosions. However, important insights can be gained by using an alternative approach in which the ejecta properties are defined analytically by using insights from progenitor and explosion models as a guide. One can then modulate such models without any effort. For example, one can vary the abundance profiles, the abundance ratios, the composition, or the chemical mixing without having to produce a new progenitor with a stellar evolution code and a new ejecta with a radiation hydrodynamics code for each new model.

In this work, we use the shell masses and representative abundance ratios of the main shells (as obtained in the grid of models presented inthe previous section) to craft our ejecta composition. We assume that the ZAMS mass determines the final core properties. We further assume that the effect of mass loss is limited to trimming the star, thus causing no impact on the final core properties. This is largely corroborated by our MESA simulations of single stars – all models in the 10− 25M_⊙ range start losing mass when reaching the RSG phase, and continue to lose mass until core collapse. In practice, mass loss can make a massive starevolve as if it was of a lower mass. However, this latter assumption only implies a shift between the ejecta properties and its corresponding ZAMS mass. Similar shifts can be caused by rotation, in the sense that a rotating star tends to produce core properties similar to those of a more massive but non-rotating star. These slightshifts in shell masses or final masses do not have any significance for the results discussed in this paper.

We generate ejecta models corresponding to 12, 15, 20, and 25 M_⊙ stars. The assigned masses for the He, CO, and Si cores are based on the values obtained in the stellar evolution calculations with MESA and described in the previous section (see also Table 2). We add a fixed H-rich shell of 9 M_⊙ to the He-core (some type II SN models are also done by forcing the total ejecta mass to be 10 M_⊙). Although the Fe-core mass covers the range 1.5−1.97, we set it to 1.5 M_⊙ for all SN models. Table 2 summarizes the basic progenitor properties adopted in our grid of models. We adopt a fixed ejecta kinetic energy of 10⁵¹ erg and a fixed ⁵⁶Ni mass of 0.08 M_⊙ for all models unless otherwise stated (Sect. 8 describes the influence of varying the ⁵⁶Ni mass on the SN radiation properties).

A further simplification is to limit the composition to the main elements witnessed in core-collapse SN spectra at nebular times. The radiative transfer process is the transformation of the decay power absorbed by the ejecta into escaping low-energy radiation. So, the gas properties are controlled by the coolants that balance the instantaneous decay heating (this balance is exact at nebular times). For the present exploration, we limit the composition to those atoms and ions that produce visible features in type II SN spectra at nebular times. In full simulations, left to a forthcoming study, we will include all species that play a role at some depth in the ejecta, even if those species produce no spectral mark.

We thus limit the composition of our ejecta to H, He, N, O, Mg, Si, Ca, Fe, Co, and Ni. We then prescribe the mass ratio for the dominant species of each shell, following the results from our MESA calculations. The H-rich shell is made of 68% H and 28% He (all meant by mass). The He-rich shell is made of 99% He and 0.5% N (some models were inadvertently made with 2% N). The O-rich shell is made of 90% O and 10% Mg. The Si-rich shell is made of 90% Si and 10% Ca. We prescribe a ⁵⁶Ni profile (function of the adopted mixing; see below). At nebular epochs, ⁵⁶Ni has decayed into a mix of ⁵⁶Fe, ⁵⁶Co, and ⁵⁶Ni that have the same distribution in velocity space and a cumulative mass equal to the initial ⁵⁶Ni mass (⁵⁸Ni and ⁵⁹Co are also included as part of the solar metallicity mixture).

In all shells we include the metals at their solar metallicity value whenever their atomic mass is greater than that of the dominant element in that shell (for example, we set Ca at the solar metallicity value in the O-rich shell but the O mass fraction is zero in the Si-rich shell). This property is motivated by our MESA results (see Sect. 2.4 for limitations of this choice). Computing models with a different metallicity is straightforward. Once all these mass fractions are set in each shell, we renormalize to unity at each depth (this introduces a change at the 10% level). The mass ratios and shell masses define the default setup for all models. One can easily adjust a given abundance ratio or a shell mass to test the influence on the resulting SN observables.

Having set the distribution of elements in mass space, we then set the density versus velocity with the constraint that the kinetic energy should be 10⁵¹ erg. We primarily focus on one epoch, 300 d after the explosion. However, we also investigate the evolution with time in Sect. 9. The radius of each ejecta shell is directly set by the velocity for homologous expansion.

The specification of the density versus velocity follows the approach of Chugai et al. (2007). The ejecta density distribution ρ(V) is given by (1)

where ρ₀ and V₀ are constrained by the adopted ejecta kinetic energy E_kin, the ejecta mass M_ej, and the density exponent k (set to eight) through (2)

and where (3)

Finally, our guess for the initial temperature was a uniform value of 5500 K. Under all conditions tested, this allows CMFGEN to converge steadily and robustly to the solution. With this choice, the converged temperature tends to be lower in the inner metal rich regions and higher as we progress in the H-rich shell, though the temperature typically stays within a few 1000 K at most of this initial guess of 5500 K. The final temperature profile (and the ionization, etc.) depends on many ejecta properties including the mass of ⁵⁶Ni and its distribution within the ejecta.

The composition that results from the above prescriptions presents a jump at the edge of each shell. We use a Gaussian smoothing to broaden the composition profiles to mimic mixing. This is a tunable parameter so weak mixing (Gaussian width of 100 km s⁻¹) or strong mixing (Gaussian width of 400 km s⁻¹) were implemented. For greater freedom, the mixing of ⁵⁶Ni and other species is chosen to be independent. The adopted initial profile for ⁵⁶Ni is of the form (4)

and connects continuously to a constant ⁵⁶Ni mass fraction for V < V_Ni. The normalization is set by the specified ⁵⁶Ni mass initially, which is 0.08 M_⊙ by default. By varying V_Ni and ΔV_Ni, we can enforce various levels of mixing of ⁵⁶Ni and its daughterelements, and therefore tune the spatial distribution of the absorbed decay power. Figure 1 illustrates the composition stratification for the type II SN ejecta corresponding to the 20 M_⊙ ZAMS mass.

Fig. 1 Cumulative composition and velocity profiles versus Lagrangian mass in the ejecta of the SN II model based on the 20 M⊙ ZAMS star. The sum for each plotted element i includes the mass fractions for all plotted elements that have a lower atomic mass.

2.3 Radiative transfer modeling during the nebular phase with CMFGEN

Using the ejecta configurations described above, we solve the non-local thermodynamic equilibrium (non-LTE) radiative transfer problem with CMFGEN (Hillier & Miller 1998; Dessart & Hillier 2005; Hillier & Dessart 2012). Unlike in recent years, we here assume steady state so that a large number of simulations can be done independently without any knowledge of the previous evolution. We use the fully relativistic transfer equations which are solved in a similar manner to the transfer equations discussed in Hillier & Dessart (2012) and references therein. During the nebular phase, and at the typical velocities encountered in SNe, these give, for all practical purposes, the same results as found using the time-dependent transfer equation.

We treat non-thermal processes as per normal (Dessart et al. 2012; Li et al. 2012). We limit the radioactive decay to the ⁵⁶Ni chain. For simplicity, we compute the non-local energy deposition by solving the grey radiative transfer equation with a grey absorption-only opacity to γ-rays set to 0.06 Y_e cm² g⁻¹, where Y_e is the electron fraction. The model atoms included are: H I (26,36), He I (40,51), He II (13,30), N I (44,104), N II (23,41), O I (19,51), O II (30,111), Mg I (39,122), Mg II (22,65), Si I (100,187), Si II (31,59), Ca I (76,98), Ca II (21,77), Fe I (44,136), Fe II (275, 827), Fe III (83, 698), Co II (44,162), Co III (33,220), Ni II (27,177), and Ni III (20,107). The numbers in parentheses correspond to the number of super levels and full levels employed (for details on the treatment of super levels, see Hillier & Miller 1998).

Obviously, with the limited composition (and associated limited model atom), some lines will not be predicted in our simulations. This includes, for example, Na I D so that the feature we predict around 5900 Å is primarily due to He I 5875 Å. Line blanketing associated with Ti II is also neglected, although it tends to be much weaker at nebular times. Simulations with a detailed composition and the associated detailed model atom will be used in a forthcoming study.

For the initial conditions needed by CMFGEN, we assume a uniform temperature of 5500 K throughout the ejecta, partial ionization of the gas, and all level populations are at their LTE value. Convergence to the non-LTE solution with a new temperatureand electron density profile takes about 200 iterations, but only 12h of computing time for a Doppler width of 50 km s⁻¹ (see Sect. 3). However, once a given model is converged, variants of that model (for example the same ejecta model but with a different ⁵⁶Ni mass or metallicity, etc.) can be computed quickly by adopting this converged model as an initial guess for the temperature, level populations, etc.

2.4 Caveats

In our simplified approach, the properties of our adopted ejecta are not accurate nor consistent. We use the pre-SN composition of the main shells (which result from hydrostatic burning) to describe the SN ejecta composition (which results in part from explosive burning), with only one change associated with the presence of a ⁵⁶Ni-rich shell. In practice, this simplification amounts to introducing a ⁵⁶Ni-rich shell at the base of the ejected He-core material, leaving the composition of this pre-collapse He core unchanged.

This departs from what may occur in realistic explosions, whereby a large fraction of the Si-rich shell collapses into the compact remnant, while the ⁵⁶Ni-rich and Si-rich shells are recreated from explosive burning at the base of the O-rich shell. Strictly speaking, the yields from explosive and hydrostatic burning differ (see, e.g., Arnett 1996), but these differences are small to moderate. For the current study, the essence is simplification so that we can explore a variety of sensitivities of SN radiation to ejecta properties, focusing on the strength of the strongest lines, the influence of ⁵⁶Ni on Ca ionization and other ejecta properties. We have already performed simulations with a state-of-the-art ejecta composition (based on the ejecta models of Woosley & Heger 2007) and these do not yield critical differences that would impact the conclusions drawn in this paper (Dessart & Hillier, in prep.).

Another limitation in our work is that the original shell compositions are based on MESA simulations that are performed with a small nuclear network of only 21 isotopes. Although incomplete, it is standard practice in the community for core-collapse SN simulations since it captures the key physics while keeping the computing time small. It is well known that both resolution and network size alter the results of simulations of massive star evolution (Farmer et al. 2016) and explosion (Paxton et al. 2015), but not to an extent that would alter the results of the present study. For example, the s-process can lower the Ca and Fe abundance within the O-rich shell. The s-process is ignored by our small network and we find instead that the Ca and Fe mass fractions in the O-rich shell are essentially constant in our MESA simulations. In the SN model s15 of Woosley & Heger (2007), based on a 15 M_⊙ star initially, thes-process influences the Ca and Fe abundances so that they show a depression relative to solar by a factor of 5− 10 in the O/C shell, and by a factor of 2−3 in the O/Ne/Mg shell. Ultimately, such features should be properly handled, but we ignore them for the time being. Since the Ca abundance is already 10 000 times lower than the O abundance in the O-rich shell (in the absence of shell merging), a further reduction by a factor of a few does not appear a critical change. When crafting our models, the same assumptions are applied to all ejecta models. Our study focuses on trends, not on accurate quantitative assessments of line fluxes and ejecta yields. Slight offsets in abundances are therefore irrelevant for the conclusions we make.

2.5 Set of simulations

The parameter space that can be explored with our setup is large. One can change parameters in isolation and then try all possible permutations in a systematic way. We choose to be selective and vary parameters that reflect the range of possibilities consistent with our current understanding of stellar evolution and explosion physics. This includes the level of mixing of all species and the separate issue of ⁵⁶Ni mixing, as well as the influence of mass loss in peeling the star’s envelope, yielding a pre-SN star with various amounts of H-rich material.

In the sections below we present results for simulations in which we varied the adopted Doppler width; the level of mixing either of all species or of ⁵⁶Ni or sometimes both; the mass of the H-rich envelope in the progenitor; the mass of ⁵⁶Ni; and the initial star mass. The mass fraction of elements within shells was also varied, in particular to test the influence of the O/Ca mass fraction ratio in the O-rich shell. We varied the metallicity from 0.1 to 2 times solar but found no influence on our results. For all simulations presented here, we set the metallicity to solar. Unless otherwise stated, all simulations adopt a SN age of 300 d.

We do not provide a summary table of all models presented since the models have many properties. Instead, in each section, we emphasize the given parameter that is varied while other parameters are kept identical. Our interest is in the influence of that varying parameter, not the numerous other parameters characterizing each individual model. In each section, the model nomenclature is clear enough to interpret the results and identify each model.

3 Influence of the adopted Doppler width

For nebular phase radiative transfer modeling, most codes use the Sobolev approximation (e.g., Jerkstrand et al. 2011), which is equivalent to assuming an intrinsic line width of zero. In CMFGEN, this approximation is not used and all lines have a finite width.

Two terms control the intrinsic³ broadening of lines. The first one is associated with the thermal velocity of the corresponding atom or ion, and typically amounts to a few km s⁻¹ at the low temperatures of SN ejecta at nebular times. This velocity width scales with , where A is the atomic mass, so the line width of iron group elements is roughly a factor eight narrower than those of hydrogen at the same gas temperature. The other term is associated with turbulence. In our simulations we typically use the same broadening for all lines, and do this by ignoring the thermal contribution and by setting the turbulence to 50 km s⁻¹. This choice is motivated by speed, and the need to have a “reasonable” number of frequency points (still of order 10⁵) to cover the full spectrum from the far ultraviolet to the far infrared.

Tests in this study show that the adopted Doppler width influences the resulting spectra at nebular times The choice of 0 km s⁻¹ implied by the Sobolev approximation is the opposite extreme and is probably not optimal since it prevents any line overlap (and hence interaction) within the Sobolev resonance zone. We adopted a value comparable to that implied by the thermal motions of the gas for IGEs.

In this work, we first converged all models using a fixed Doppler width of 50 km s⁻¹ for all species. This produces a converged model that is pretty close to the true solution (the essential part of the solution being the converged temperature, ionization, and level populations). This converged model is then used as initial conditions for a new CMFGEN model in which a fixed Doppler width of 2 km s⁻¹ is used for all species (in cases where the impact on the gas properties was large, this reduction was done in a few steps)⁴. The frequency grid is set so that all lines are resolved, irrespective of their widths. So, reducing the Doppler width of lines from 50 to 2 km s⁻¹ leads to an increase in the number of frequency points from 49 352 to 841 208 in the present simulations. This model takes much longer per iteration but it requires fewer iterations to converge since it starts with a better guess. With the lower Doppler width, the temperature and ionization are slightly modified and cause a sizable change in the strength of the strongest lines and the Fe II emission forest around 5000 Å. The evolution can be nonlinear, in the sense that reducing the line Doppler width can yield a non-monotonic increase or decrease in the flux of [O I] λλ 6300, 6364, Hα, or [Ca II] λλ 7291, 7323.

Figure 2 shows the influence on the optical and near-infrared spectrum of changing the Doppler width from 50 to 10 and 2 km s⁻¹ in a strongly mixed ejecta model corresponding to a 15 M_⊙ progenitor. Reducing the Doppler width leads to a strengthening of Hα and Pa α (and a weakening of the near-infrared Mg I lines) which may have resulted from the increase of the electron density and the H ionization in the recombined H-rich layers of the ejecta. The same test in a model with weaker mixing would yield a different impact because the modified mixing has a strong influence on the emergent spectrum. So, for example, in a weakly mixed model, the reduction of the Doppler width leads to a reduction of He I 7068 Å and Mg I lines in the near infrared, a strengthening of the Ca II lines in the optical, but has little impact on the H I lines. It is thus important to realize that the adopted Doppler width can influence nebular line ratios.

The influence of line overlap on non-LTE processes is well documented. For example, overlap of an O III line with He II Ly α in planetary nebula leads to anomalous line strengths for several optical O III lines (e.g., Bowen 1934; Osterbrock 1989). In some O stars the chance overlap of Fe IV lines with the He I resonance transition at 584 Å influences the strength of He I singlet transitions in the optical (particularly those involving the 1s 2p ¹P^o level; Najarro et al. 2006). Another example, is the overlap of an O I line with Ly β which can, for example, influence the strength of O I λ 8446 (e.g., Osterbrock 1989). As discussed by Jerkstrand et al. (2012), a significant concern for modeling the nebular phase of Type IIP SNe is the overlap of some Fe II and O I lines with Ly α and Ly β. We verified the cause of the spectral changes by running test calculations which used a Doppler width of 50 km s⁻¹ and modified model atoms for which we artificially reduced the oscillator strengths (effectively setting them to zero) of Fe II and O I lines overlapping with Ly α and Ly β.

Fig. 2 Influence of the adopted Doppler width on the resulting optical and near-infrared spectra for a SN II model arising from a 15 M⊙ progenitor. A weak mixing of all species is used except for 56Ni for which we adopt a strong mixing. The increase in Hα line flux with decreasing VDop is at the expense of a decrease in flux from a forest of Fe II lines in the 4000− 6000Å region. Other lines affected are He I 1.083 μm, O I 1.129 μm, Fe I lines around 1.165 μm, Mg I lines at1.183, 1.502, 1.577, 1.711, H I 1.875 μm. Note that about 70% of the total flux falls in the optical range.

4 Discussion of model results from a 20 M_⊙ type II SN model

We first describe the case of a type II SN ejecta from a 20 M_⊙ ZAMS star. The model is named m20mix100vnimH9, where “mix100” means that the Gaussian smoothing uses a characteristic width of 100 km s⁻¹ (this mixing is applied to all species), “vni” means that the additional mixing applied to ⁵⁶Ni was strong (V_Ni = 2500 km s⁻¹ and ΔV_Ni = 1000 km s⁻¹; see Eq. (4)), and “mH9” means that the pre-SN progenitor had a 9 M_⊙ H-rich outer shell. The ejecta has a kineticenergy of 10⁵¹ erg and contains 0.08 M_⊙ of ⁵⁶Ni initially.

A summary of results for this reference model is shown in Fig. 3. The top two panels show the optical and near-infrared spectral properties. Throughout the paper, we show luminosities L_λ (in erg s⁻¹ Å⁻¹) rather than scaled or normalized fluxes because its integral over wavelength yields the bolometric luminosity and thus the decay power absorbed by the ejecta. Variations in ⁵⁶Ni mass or γ-ray trapping efficiency will yield different brightnesses, which can be assessed from L_λ. The bottom panels show various properties of the ejecta, including the ionization level (electron density, O and Ca ionization), the temperature, the composition (for H, He, O, or Ca, and for the ratio of O/Ca, which is relevant for the [O I] λλ 6300, 6364 doublet strength in core-collapse SNe; see next section), and the profile of the decay power absorbed.

The total radial electron scattering optical depth for this 300 d-old ejecta is 0.4. Under such conditions, the ejecta no longer traps radiation so the energy balance is set by the distribution of decay power absorbed at that time in the ejecta and the reprocessing of this power by the gas in the form of low-energy, mostly optical, photons – whatever is absorbed is instantaneously radiated and the balance is zero. Because of the strong ⁵⁶Ni mixing and some non-local energy distribution, the decay power is absorbed throughout the ejecta. Relative to the total energy absorbed, the Si-rich layers below 1000 km s⁻¹ get ~ 20%, the O-rich layers between 1000 and 1800 km s⁻¹ get ~ 40%, the He-rich layers between 1800 and 2100 km s⁻¹ get ~ 10%, and the outer H-rich layers receive the remaining ~30%. Interestingly, the amount of radiation escaping to infinity and arising from these shells is very different, as shown in Fig. 4. In other words, there is a significant reprocessing of these low-energy photons by the ejecta.

Figure 4 illustrates the spatial origin of the emerging flux and more specifically the quantity δL_λ (R) – it illustrates the location of the last interaction (not necessarily the location where the photon was originally emitted). The fractional contribution to the total flux at wavelength λ from a narrow ejecta shell at radius R is: (5)

where Δz is the projected shell thickness for a ray with impact parameter p, η is the emissivity along the ray at p and z, and τ is the total ray optical depth at λ (i.e., the integral is performed in the observer’s frame) at the ejecta location (p, z). Figure 4 shows that the bulk of the flux emerges from the H-rich ejecta layers, which is 61% of the total and thus about twice as much as deposited by γ-rays and positrons. The He-rich layers radiate 15% of the total, compared to the 10% of the decay power that they absorb. The O-rich (Si-rich) layers radiate 20% (4%) of the total, compared to the 40% (20%) of the decay power that they absorb. In the absence of optical depth effects, the power from each shell would be equal to the decay power absorbed in each corresponding shell. With a total electron-scattering optical depth of 0.4 at the SN age of 300 d, the ejecta is not thin. It is not optically thick enough to cause radiation storage and diffusion with a sizable delay, but it is optically thick enough to cause strong reprocessing of UV and blue optical photons, which benefits the emission from the outer layers in this model. Such optical depth effects and their wavelength-dependent impact are illustrated in Fig. 5.

Figure 4 also illustrates the complicated line formation process at nebular times. It helps in solving the problem of line overlap since a feature may appear broad with a narrow peak because it forms throughout the ejecta (from small to large velocities) or because it forms in a single narrow shell but overlaps with adjacent lines. Let us consider the formation of the strongest optical lines, which are mostly forbidden. The [O I] λλ 6300, 6364 doublet forms primarily in the O-rich shell (yellow line in Fig. 4), with a small contribution at large velocity from the H-rich shell (blue line) and a small contribution at low velocity from the He-rich shell (red line). This contribution from the He shell arises from our imposed mixing and is confined to the region at the interface between the two shells (Fig. 1).

The [Ca II] λλ 7291, 7323 doublet forms throughout the ejecta – from the outer part of the Si-rich shell and the inner part of the O-rich shell, and from the H-rich shell (very broad component). The [Ca II] λλ 7291, 7323 doublet emission extends up to 7400 Å because of the contribution from the H-rich shell but also because of the influence of electron scattering in the fast moving H-rich layers. The emission further to the red is due to Fe II (strongest component at 7452 Å). The He-rich shell, which is nearly exclusively He with a small amount of N, contributes mostly through the emission of N II 6548 Å, but this emission overlaps with the strong and broader Hα line seen at all times in SNe II. The N II identification is thus nontrivial in SNe II, but it is observed and explained in some SNe IIb (Jerkstrand et al. 2015b). Below 6000 Å, most of the flux emerges from the H-rich layers and stems mostly from Fe I and Fe II line emission (with a mix of permitted and forbidden transitions). The strongest Fe II line is an isolated forbidden line at 7155.2 Å. There is a weak Mg I 4571 Å line from the O-rich shell (as well as the O I 7774 Å further to the red). The nucleosynthetic signatures of the explosion and the pre-SN evolution are thus quite limited, with the main features being [O I] λλ 6300, 6364 and [Ca II] λλ 7291, 7323. A more extensive list of line identifications is provided in the upper panels of Fig. 3.

In the near infrared, the synthetic spectrum shows lines from the Paschen series of H I, He I 1.083 μm, O I 1.129 μm, Mg I 1.488, 1.502 and 1.711 μm, Si I at 1.227 and 1.422 μm, numerous Ca I lines just short of 2 μm, as well as numerous Fe I and Fe II lines spread throughout the near infrared. These identifications are often ambiguous since most features are a composite of different lines.

One aspect that controls the radiative properties of the nebular phase spectrum is the ionization (shown for O and Ca in the bottom panels of Fig. 3). Here, the mass density is constant below about 2000 km s⁻¹ so the variations of the electron density below 2000 km s⁻¹ reflect the change in ionization and mean atomic weight as we proceed through the He-rich shell, the O-rich shell, and the Si-rich shell. Oxygen is primarily neutral in the O-rich shell but Ca is once ionized in the Si-rich shell. The jump in ionization and temperature in the He-rich shell arises because the dominant species in that shell (i.e., He) is a poor coolant.

At all depths in the ejecta, the heating source is radioactive decay. However, because of the stratification in composition and the range of densities between inner and outer ejecta, the sources of cooling vary drastically with depth. In the H-rich layers, we find that the cooling is done through non-thermal excitation and ionization (primarily in association with H I, and by a factor of ten weaker with He I), and collisional excitation of Fe II (and to a lesser extent Mg II and O I). In the He-rich shell, the cooling is done through collisional excitation of Fe II and non-thermal processes tied to He I. In the O-rich shell, cooling arises primarily from collisional excitation of O I and non-thermal processes tied to O I, and then from collisional excitation of Mg I and Fe II. In the Si-rich shell, the cooling is done primarily through collisional excitation of Ca II as well as Ca I, plus non-thermal processes tied to Si I. Figure 6 illustrates these various cooling components at different ejecta depths in this reference model m20mix100vnimH9.

The above results need to be considered in view of our assumptions. First, the ejecta composition is limited to H, He, N, O, Mg, Si, Ca, and the ⁵⁶Ni decay chain elements so the lines present in the models are by design limited to these elements (in their neutral or once ionized state). The second aspect is that we adopt a strong mixing of ⁵⁶Ni with a weak mixing of the other (lighter) elements. This gives some preference to the intermediate and high velocity layers of the ejecta (although we have seen that the outer layers do a lot of reprocessing of the radiation emitted from the inner layers). Finally, because of the weak mixing of species other than ⁵⁶Ni, the metal-rich layers retain an onion-like shell structure. With macroscopic mixing only, material from all shells in the pre-SN star would coexist at a given velocity. The reprocessing of radiation from the inner ejecta by the various lumps of material would be much more complex than in the present configuration of shells of distinct composition stacked on top of each other.

Fig. 3 Illustration of properties for the reference model m20mix100vnimH9 discussed in Sect. 4. In this and similar figures, the top panels show the luminosity Lλ (integrating over wavelength yields the bolometric luminosity) in the optical (upper) and the near infrared (lower). The line identifications are indicative only since in numerous cases multiple lines contribute (we give the primary component to most features). Forbidden lines are shown in red for all ions and atoms apart from Fe (shown in blue). Bottom panels: some ejecta properties computed with CMFGEN.

Fig. 4 Illustration of the spatial regions (here shown in velocity space) contributing to the emergent flux in a 20 M⊙ type II SN model characterized by weak mixing for all non-IGE species, strong mixing for IGE species, and with a 9 M⊙ H-rich envelope mass (model m20mix100vnimH9). Top panel: observer’s frame luminosity contribution δLλ,R (Eq. (5); the map maximum is saturated at 20% of the true maximum to bias against the strong Hα line and better reveal the origin of the weaker emission) versus wavelength and ejecta velocity. The four contributing shells in this type II SN model are clearly seen (see vertical colored stripe at left), although the H-rich layers contribute most of the emergent flux. Bottom panel: emergent luminosity integrated over each shell, together with the total luminosity. In this model, the fraction of the total power emerging from the Si-rich, O-rich, He-rich and H-rich shells is 4, 20, 15, and 61%, respectively. This illustrates the contribution to the emergent total flux from the main ejecta shells (this may also be inferred from the width of lines but complicated when multiple regions contribute or when there is line overlap). Note that with our approach, IGEs are present in all shells and thus Fe I and Fe II flux contribution is present throughout the ejecta. One can also see that the flux from the Si-rich shell is mostly radiated by Ca I and Ca II (together with Fe I and Fe II, rather than by the dominant species Si).

Fig. 5 Illustration of the total (i.e., accounting for all opacity sources) radial optical depth integrated from the outer ejectaboundary to the innermost boundary (i.e., the base of the ejecta; blue), to the He/O shell interface at 1780 km s−1 (red), and to the He/H shell interface at 2150 km s−1 (yellow) for the reference model m20mix100vnimH9.

Fig. 6 Illustration of the main cooling processes balancing the radioactive decay heating at all ejecta depths in the reference model m20mix100vnimH9. 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-thermalexcitation and ionization) and “COL” stands for collisional processes (i.e., collisional excitation).

5 Importance of the Ca mass fraction in the O-rich shell

Converting line strengths into abundances of the associated ion and species is one of the principal goals of nebular phase spectroscopy for any type of SN. There are however many caveats associated with this task. First, the radiation emitted by any ion is very dependent on the atomic physics of that ion. Second, it depends on the efficiency of other species that are radiating from the same region. Third, it depends on where the radioactive decay poweremitted is absorbed by the ejecta. Hence, a fundamental aspect of nebular line emission is how the ⁵⁶Ni is mixed through the ejecta and how far the emitted γ-rays can travel in the ejecta.

Having determined the non-local distribution of this decay power and which fraction of the total is shared between each dominant ejecta shell (H-rich shell, He-rich shell, O-rich shell, and the Si-rich shell), the radiation emitted in each shell will occur through the lines that have the strongest cooling power. The cooling efficiency of a given forbidden line depends of course on the abundance of the corresponding ion (hence a function of mass fraction and ionization level), but also on atomic properties of the ion levels (oscillator strength, critical density, etc.). In practice, type II SN ejecta exhibit a wide disparity in ionization, composition, and density, and the various constituents have widely different atomic properties so that even trace elements can dominate the cooling of the gas (the same holds in H II regions, whose cooling is controlled to a large extent by emission in lines of N II, O II and O III while the dominant species are instead H and He; Osterbrock 1989).

Below, we discuss the case of the doublet forbidden transitions [O I] λλ 6300, 6364 and [Ca II] λλ 7291, 7323 because they are the strongest optical lines in type II SN spectra at nebular times (if we exclude Hα) and also because they are routinely used to set constraints on nucleosynthesis and progenitor properties. We then present simulations in which various amounts of Ca are introduced into the O-rich shell. These simulations show the strong impact Ca can have on O I line emission. We then discuss the implications and comment on previous work.

5.1 The cooling power of [O I] λλ 6300, 6364 and [Ca II] λλ 7291, 7323

Let us consider a small gas volume of uniform composition with a gas temperature T (or T₄ when expressed in units of 10⁴ K) and electron density N_e. Below, we study how the emissivity of (or the cooling rate associated with) the [O I] λλ 6300, 6364 doublet compares with that of the [Ca II] λλ 7291, 7323 doublet. We assume optically-thin emission in this simplified analysis. This approximation may not hold strictly, in particular for large Ca densities.

The collisional de-excitation rate per unit volume per unit time from the upper level “u” to the lower level “l” is given by (6)

where N_u is the upper level population and ϒ_lu is the effective collision strength for the transition (data for collision strengths are taken from Mendoza 1983 for O I and from Meléndez et al. 2007 for Ca II; energy levels are from NIST, while A_ul values for O I are from Osterbrock 1989 and those for Ca II are similar to those in Lambert & Mallia 1969)⁵.

At the critical density, this rate is equal to N_uA_ul. Thus (7)

For O I we sum the A_ul values for the two transitions, to get a critical density. For Ca II, we treat the upper level as a single level (A = ∑ g_uA_ul∕∑ g_u, and ϒ = ∑ ϒ_lu). In the following, we assume O I and Ca ii are the dominant ionization stages, and that levels other than those listed above can be neglected. We use the atomic properties listed in Table 3.

We may consider two limiting cases for the strength of these forbidden lines, corresponding to electron densities much higher or much lower than the critical density of the transition.

Case 1: N_e ≫ N_e (crit). In this case we can assume the upper level is in LTE relative to the ground state. Thus, the O I line emissivity is given by (8)

or (9)

Similarly, for the Ca II line (g_u∕g_l = 10∕2), we have (10)

Integrating these line emissivities over the O-rich shell volume yields the line luminosities L. Taking the ratio we obtain, (11)

Case 2: N_e ≪ N_e (crit). In this case, the line emissivity is (12)

since every excitation gives a line photon. Thus, the ratio of line luminosities is now (13)

In this case we could add a factor like since the O I collision strength grows faster than that of Ca II. A factor of 2.5 (i.e., 40/16) is introduced if we express Eqs. (11) and (13) with mass fractions.

In the optically-thin limit, the luminosity contrast between these two lines would be huge if Ca was as abundant as O in the O-rich shell (in practice, line optical depth effects would reduce this Ca II emission). Equations (11) and (13) show that the Ca II line would dominate over the O I line if the Ca abundance is at least ~ 1% of the O abundance in the O-rich shell. This situation does not occur if Ca has the solar abundance in the O-rich shell, but the merging of the Si-rich and O-rich shells may change this. The impact would depend on the mass of the Si-rich and O-rich shells in the pre-SN star (and prior to merging). These shells masses grow with the progenitor mass, with a 2− 3 times faster increase for the O-rich shell mass.

5.2 The influence of O/Ca in the O-rich shell: numerical simulations with CMFGEN

Using as a reference the model m20mix100mH9 produced with the prescriptions stipulated in Sect. 2.2, we produce three other models in which the Ca mass fraction in the O-rich shell is raised from its solar value of 6.4 × 10⁻⁵ to 0.0007 (model suffix Cal), 0.007 (suffix Cap), and finally 0.07 (suffix Capp). The CMFGEN results for both the optical and near-infrared radiation and the gas properties are shown in Fig. 7.

As expected, raising the Ca mass fraction in the O-rich shell has a dramatic effect on the [O I] λλ 6300, 6364 doublet strength (note that the distribution of the decay power absorbed is the same in all four models since we switch O for Ca but at constant density and mass; see middle panel). The effect is present even for an O/Ca mass fraction ratio in the O-rich shell of 1000. The influence on Hα is negligible, which is expected since the H-rich layers were not modified. The Ca II lines are strengthened, but for the higher Ca enhancement, their strength decreases to the benefit of Ca I lines (mostly present in the near infrared; bottom spectral panel of Fig. 7). The ionization is still primarily Ca⁺, but there is a small inflection in the ionization, the electron density, and in the temperature in the O-rich shell. The most likely reason for this saturation and even reduction of the Ca II flux is that both doublet components are optically thick (Fig. 8).

5.3 Implications and comparison to previous work

A number of points need to be made at this stage. Rather than being secondary and irrelevant, the Ca mass fraction in the O-rich shell is a critical matter because it affects the O I doublet, which is used for constraining the progenitor mass. In the above experiment, the mass of the O-rich shell is the same in all four ejecta models and yet the O I doublet flux varies by a factor of about five just by tuning the Ca/O ratio in the O-rich shell.

Some of the adopted Ca mass fractions in this experiment are probably too large. However, stellar evolution simulations frequently produce a large Ca mass fraction in the O-rich shell. This occurred in half the MESA simulations we ran for this study. The cause isthe merging of the Si-rich shell and the O-rich shell during Si burning, producing a Ca mass fraction of a few 0.001 in the O-rich shell, so that the O/Ca mass fraction ratio drops by a factor of 100 compared to the case of no shell merging (the total Ca mass is however unchanged by this merging). A similar feature was observed by Fransson & Chevalier (1989) in some of their models (computed with KEPLER and presented in Ensman & Woosley 1988). Collins et al. (2018) report a higher occurrence of merging for the O, Ne, and C burning shells in higher mass progenitors between 16 and 26 M_⊙, yielding a large Si mass fraction in a large part of the O-rich shell – this is the same as what we observe in our MESA simulations. Similarly, Yoshida et al. (2019) obtain distinct composition profiles within the Si- and O-rich shells depending on the adopted convective overshoot strength and progenitor mass.

In the absence of such a merging, the Ca mass fraction in the O-rich shell is typically at the original value on the ZAMS (i.e., at the solar metallicity in our models; see Sect. 2.4 for departures from a solar value), but one can wonder whether this feature would persist in 3D hydrodynamical simulations treating physically the processes of convection and overshoot. Even with strong overshoot from the Si-rich shell, it is not clear that mixing could take place down to the microscopic level given the short time until collapse and explosion (typically a day from the onset of Si burning; see, e.g., Arnett 1996 or Collins et al. 2018) – the mixing might instead be partial and truncated at some spatial scale.

The origin of Ca II emission is a related issue. In the experiment above, if the Ca mass fraction is large in the O-rich shell (say with a ~ 0.01 mass fraction), most of the power absorbed by the O-rich shell is radiated through Ca II lines (this process can be mitigated by line optical depth effects). Otherwise, we obtain Ca II emission from the Si-rich shell, the interface between the Si-rich shell and the O-rich shell (where the Ca mass fraction is around or above 0.01), and from the He-rich and H-rich shells.

Li & McCray (1993) and others argue that the Ca II emission arises primarily from the Ca in the H-rich envelope. The main limitation for Ca emission from the Si-rich shell is that the Si-rich shell is of low mass, and hence it absorbs a small fraction of the total decay power. Even if Ca were the sole coolant in that Si-rich shell, the total power that it would radiate would be a small fraction of the total, typically smaller than the power absorbed in the H-rich layers (though this can depend on the mass of the CO core and the level of mixing). More importantly, the main competitor is the CO core mass, which in massive progenitors may absorb most of the decay power (again a function of mixing, etc.). Whether Ca emission comes from the H-rich envelope is not just a matter of Ca being a strong coolant. It depends primarily on where the decay power is absorbed. For a larger He core mass (a higher mass progenitor), for a larger ejecta mass, or for a weaker level of ⁵⁶Ni mixing, a larger fraction of the decay power may be absorbed in the core rather than the H-rich layers. In this situation, the H-rich ejecta layers would only reprocess the radiation impinging from below.

The origin of the Ca II emission would be altered in the case of Ca pollution into the O-rich shell prior to core collapse. This would enhance the [Ca II] λλ 7291, 7323 doublet flux at the expense of the [O I] λλ 6300, 6364 doublet flux. If the process of shell merging had a higher occurrence rate in higher mass progenitors (Collins et al. 2018), it would provide a natural explanation for the lack of SNe II with strong [O I] λλ 6300, 6364 at nebular times, and inferred progenitor masses below 17 M_⊙.

6 Influence of the H-rich envelope mass

In our toy setup, we either force the ejecta mass to be 10 M_⊙ (and we adjust the H-rich envelope mass) or we force the H-rich envelope mass to be 9 M_⊙ (and we adjust the ejecta mass to be the sum of the H-rich envelope mass plus the He-core mass). The latter mightbe a good description of single stars initially in the mass range between 12 and about 20 M_⊙ and the latter more suitable to higher mass stars but still dying as RSGs (Woosley et al. 2002; Dessart & Hillier 2019b). Since we assume the same ejecta kinetic energy, variations in total mass or envelope structure will impact the distribution of elements in velocity space.

Figure 9 illustrates the impact of the H-rich envelope mass of the progenitor (or the total mass of the ejecta) for a 25 M_⊙ progenitor model (this defines the core properties). The low-envelope mass model is m25mix100vni (total ejecta mass is 9.9 M_⊙) and the high-envelope mass model is m25mix100vnimH9 (total ejecta mass is 16.7 M_⊙) in our nomenclature. With the higher H-rich envelope mass, a greater fraction of the decay power emitted is absorbed (90% compared to 67%). This implies a general flux offset. A sizable fraction of this extra power is radiated in H I lines. Furthermore, for a greater ejecta mass and the same ejecta kinetic energy, the metal-rich regions are located at smaller velocities. This causes a small reduction in the width of lines associated with O I and Mg I (lines forming in the extended H-rich layers are less affected). That effect is exacerbated by the greater ejecta density. This, in turn, enhances the local deposition of the decay power in the slower denser layers that contain more ⁵⁶Ni initially (the γ-ray mean free path is shorter). Hence, the progenitor H-rich envelope not only affects the type II-P SN radiation during the photospheric phase (by modifying its length and brightness), it can also influence the nebular phase properties.

Fig. 7 Similar to Fig. 3, but now showing variants of model m20mix100mH9 in which the Ca mass fraction is progressively increased in the O-rich shell by powers of ten and starting from about 7 × 10−5, the solar metallicity value (the Ca mass fraction profile is shown in the bottom left panel). Notice the dramatic reduction on the [O I]/[Ca II] line ratio as the Ca abundance is increased. Even in the richest Ca model, the Ca mass is still less than 10% of the O mass.

Fig. 8 Illustration of the Sobolev line optical depth for the Ca II 7323.9 Å line for the models shown in Fig. 7.

7 Influence of ⁵⁶Ni mixing

In this section, we explore the influence of ⁵⁶Ni mixing on the SN radiation properties. We use the core properties of a 20 M_⊙ progenitor model and assume a 9 M_⊙ H-rich envelope. Although not essential, we set in all models a 10% H mass fraction throughout the core layers (Si, O, and He-rich shells) to mimic the inward mixing of H-rich material at low velocity (in real explosions, outward mixing of core material occurs simultaneously with inward mixing of envelope material). A weak mixing is applied to all species, and a subsequent additional mixing is applied to ⁵⁶Ni with V_Ni covering from 750 to 2500 km s⁻¹ in four increments (five models). ΔV_Ni is set to V_Ni∕2.5 (see Eq. (4)). Results are shown in Fig. 10.

The ejecta model with weak mixing completely traps the decay power while the ejecta model with the highest mixing of ⁵⁶Ni lets 10% of the decay power escape at 300 d. This effect is difficult to discern from observations since about 30% of the total flux falls outside of the optical range at that time (and is thus rarely measured). Besides facilitating escape, enhanced mixing favors the more extended deposition of decay power. In the weaker mixing model, nearly 70% of the decay power is deposited below 1000 km s⁻¹ and thus within the Si-rich shell (with a small portion in the base layers of the O-rich shell). Consequently, a significant fraction of the SN radiation arises from the Si-rich and O-rich layers, as evidenced by the narrow Si, Ca, and Fe emission lines. Similarly, the faster moving H-rich layers capture little of this decay power so that H I lines are weak and quite narrow, and there is a lack of broad Fe lines.

As the ⁵⁶Ni mixing is enhanced, it first benefits the O-rich shell which produces a stronger [O I] λλ 6300, 6364 doublet. As the ⁵⁶Ni mixing is enhanced to its maximum, the [O I] λλ 6300, 6364 doublet strength drops again as a greater fraction of the decay power is now absorbed in the H-rich shell. Consequently, H I and Fe I – II lines strengthen and broaden.

Interestingly, the [Ca II] λλ 7291, 7323 doublet weakens and does not broaden, which implies that the greater power absorbed by the H-rich shell is not radiated by [Ca II] λλ 7291, 7323. The enhanced ⁵⁶Ni tends to raise the temperature and ionization (see the electron density) with increasing ejecta velocity. However, the influence on the ionization of a given species is more complicated. In the present set of simulations, the O ionization remains unchanged, but the Ca ionization is enhanced, nearly to Ca²⁺ everywhere beyond 1000 km s⁻¹ for the model with the highest mixing⁶. This increase in Ca ionization quenches the [Ca II] λλ 7291, 7323 doublet emission. This effect was previously seen in a study of Type II-P SN by Chugai (1988) who attribute it to the influence of Ly α photons. The change in electron density is largely driven by the change in H ionization.

Variations in ⁵⁶Ni mixing, and their impact on ejecta ionization, may contribute in part to the diversity of nebular-phase spectra of type II SNe (see, for example, the discussion in Yuan et al. 2016).

By simply varying the ⁵⁶Ni profile, this experiment shows that spectral lines at the nebular epoch can be considerably tuned in strength and width for the same ejecta composition, mass, and kinematics. This is well understood but it is a concern since ⁵⁶Ni mixing is a complicated process, with ⁵⁶Ni probably having a unique distribution in each SN. This is most likely tied to additional dependencies on progenitor mass and rotation, He-core mass, and explosion energy (Wongwathanarat et al. 2015; Sukhbold et al. 2016). It also shows that metal line emission (in particular of Ca and Fe in our limited composition mixture) can arise primarily from the Si-rich shell or from the H-rich shell depending on where the decay power is absorbed. The conclusion of Li & McCray (1993) that [Ca II] λλ 7291, 7323 arises from the H-rich material (which they infer from their analysis of SN 1987A) only holds if ⁵⁶Ni mixing is strong so that most of the decay power gets absorbed by the H-rich material at large velocities (an alternative is that the H-rich material is mixed inwards closer to ⁵⁶Ni and thus absorbs decay power in the absence of strong ⁵⁶Ni mixing). A further issue is that Ca over-ionization can completely inhibit [Ca II] λλ 7291, 7323 emission from the H-rich layers, no matter how much decay power is absorbed there.

8 Dependencies on ⁵⁶Ni mass

We now discuss the influence of the ⁵⁶Ni mass on the resulting optical and near-infrared spectra as well as on the ejecta properties. For this exploration, we start off with model m20mix100mH9 (20 M_⊙ progenitor, weak mixing for all species, no additional mixing of ⁵⁶Ni, imposed H-rich envelope of 9 M_⊙), characterized by a ⁵⁶Ni mass of 0.08 M_⊙ and compute additional models in which the ⁵⁶Ni mass is reset to 0.004, 0.02, 0.16, and 0.32 M_⊙. Results are shown in Fig. 11.

All five ejecta models achieve complete γ-ray trapping so there is a linear scaling between the ⁵⁶Ni mass and the bolometric luminosity. We also find that all five models radiate the same fraction of the total power within the optical range (about 75%). However, the optical spectra vary considerably with ⁵⁶Ni mass. There are various reasons for this.

For higher ⁵⁶Ni mass, the ejecta ionization is greater, yielding amongst other things an increased electron density. Going here from 0.004 to 0.32 M_⊙, which is certainly not a small range, the total radial electron-scattering optical depth to the core increases by more than a factor of ten (it grows from 0.08 to 1.1). One should recall that typical type II SN ejecta turn optically thin within a few months because of recombination. If somehow this recombination can be prevented, as can be done for example by supplying power from a magnetar (Dessart 2018), the ejecta can stay optically thick for two years after explosion.

Here, the effect is not as strong but shows the same trend. Consequently, the model with the lowest ⁵⁶Ni mass is very optically thin (total electron-scattering optical depth of 0.08 and total flux-mean optical depth of 0.8), while the model with the highest ⁵⁶Ni mass is not optically thin (total electron-scattering optical depth of 1.0 and total flux-mean optical depth of 2.2), even at 300 d. In the latter, optical depth effects and collision processes are not negligible. The greater electron density places many forbidden lines above their critical density, which will quench their emission (i.e., upper levels can more frequently de-excite through collisions). A larger ⁵⁶Ni mass also leads to a greater ionization of numerous species such as O and Ca, and to a higher temperature which facilitates populating excited states.

Summarizing, the nebular spectrum in the model with the largest ⁵⁶Ni mass shows a greater flux (built from the combination of many overlapping weak and broad emission lines), and a stronger contribution from the Fe II line forest (many of which are permitted lines). For higher ⁵⁶Ni mass, the forbidden lines of O I and Ca II are relatively weaker. The [O I] λλ 6300, 6364 doublet flux represents 30% (5%) of the optical flux in the model with 0.004 M_⊙ (0.32 M_⊙) of ⁵⁶Ni. Similarly, and in the same order, the [Ca II] λλ 7291, 7323 doublet represents 20% (6%) of the optical flux. In addition to the greater electron density, it is also the shift to a higher element ionization that quenches these lines (see O and Ca ionization in Fig. 11).

This exploration shows that the ⁵⁶Ni mass acts in a very complicated way on the ejecta properties and yields a very non-linear evolution of spectral properties. This result is a warning against using an event like SN 1987A as a reference to which linearly scale various quantities.

Fig. 9 Similar to Fig. 3, but now showing the influence of the progenitor H-rich envelope mass on the synthetic optical and near-infrared spectra at 300 d. The models shown are m25mix100vni and m25mix100vnimH9 (strong mixing of 56Ni but weak mixing of other species) and correspond to a 25 M⊙ progenitor star. There is a greater fraction of the total decay power absorbed in the H-rich layers of the lighter model m25mix100vni, boosting the [Ca II] λλ 7291, 7323 doublet emission, but this has an adverse effect on Hα, which pushes its formation to the outer lower density regions. These differences reflect in part the different distribution of decay-power absorbed since the two ejecta have a different velocity structure – mixing should also differ in nature for ejecta with a similar kinetic energy but a different mass (e.g., as in SNe II and Ibc).

Fig. 10 Similar to Fig. 3, but now showing the influence of the 56Ni mixing on the SN radiation and gas properties at 300 d. The ejecta arises from a 20 M⊙ progenitor, weak mixing for all species, imposed H-rich envelope of 9 M⊙, and a 56Ni mass of 0.08 M⊙. An additional and variable amount of 56Ni mixing is applied increasing VNi from 750 (model 20-vni7p5e7) to 2500 km s−1 (model 20-vni2p5e8). With enhanced 56Ni mixing, more decay power is absorbed in the outer ejecta, causing a strengthening of Hα but a weakening of [Ca II] λλ 7291, 7323 because of the boost to the Ca ionization (the ionization shift also affects Fe I – II and Mg I lines in the optical and the near-infrared). (See Sect. 7 for discussion.)

Fig. 11 Similar to Fig. 3, but now showing the influence of the 56Ni mass on the SN radiation (for better visibility, a logarithmic scale is used for the luminosity). The ejecta arises from a 20 M⊙ progenitor, weak mixing for all species, no additional mixing of 56Ni, imposed H-rich envelope of 9 M⊙, and is characterized by a 56Ni mass of 0.004, 0.02, 0.08, 0.16, and 0.32 M⊙. A greater 56Ni mass boosts the electron density and species ionization, inhibits forbidden line emission (in favor instead, for example, of recombination lines), and strengthens optical depth effects and continuum processes (primarily in association with electron scattering). (See Sect. 8 for discussion.)

Fig. 12 Similar to Fig. 3, but now showing the evolution of the optical and near-infrared properties of model m20mix100vnimH9 from 150 to 500 d after explosion. With increasing time, the conditions become increasingly “nebular”, with a relative strengthening of forbidden line emission, and the progressive reduction of optical depth effects (for example with the disappearance of P-Cygni profiles). This evolution is similar to reducing the 56Ni mass at a given SN age (see Fig. 11).

9 Evolution during the nebular phase

Figure 12 illustrates the evolution of the model luminosity through the optical and near-infrared ranges from 150 to 500 d for the type II SN model m20mix100vnimh9. The continuous drop in luminosity in both spectral domains results from the drop by a factor of about 23 in the decay power emitted and to a lesser extent the increasing γ-ray escape with time: 98% (66%) of the decay power is trapped at 150 d (500 d) after explosion. The simultaneous drop in density (as 1∕t³) and ionization level (as material cools and recombines) leads to a drop in electron density. Since free electrons are the dominant source of continuum opacity in type II SNe ejecta, this leads to a drop in electron-scattering ejecta optical depth, here from a value of 2.8 at 150 d to only 0.073 at 500 d (a value of ~0.3 would resultat constant ionization; the corresponding drop in flux-mean optical depth is from 3.3 at 150 d to 0.37 at 500 d). This change in electron density has numerous implications (see also discussion in Sect. 8). As the electron density drops below their critical densities, forbidden lines strengthen relative to other lines (since the decay power absorbed in the ejecta decreases with time, the strength of all lines tends to decrease). For O I, this is also facilitated by the reduction in ionization which populates the ground state of neutral oxygen.

More subtle properties can be seen from Fig. 13, which focuses on the optical range between the [O I] λλ 6300, 6364 doublet and the Ca II near-infrared triplet. The individual components of the [O I] λλ 6300, 6364 doublet have nearly equal strength at 150 d, which implies that the lines are optically thick, while at 500 d their relative strength is about three (i.e., equal to the ratio of their radiative de-excitation rates A_ul ; see Table 3). The changing flux ratio of the doublet components is discussed in detail by Li & McCray (1992). Time also has a profound impact on the relative strengths of Ca II lines. The [Ca II] λλ 7291, 7323 doublet strengthens considerably with time relative to the Ca II near-infrared triplet, a consequence of the drop in electron density (see also Li & McCray 1993). Although a limitation of assuming spherical geometry, numerous forbidden lines at 500 d have a flat-top profile, indicative of emission from a hollow sphere in velocity space (the sphere is in fact a shell, thus bounded between two large velocities)⁷.

Fig. 13 Asfor Fig. 12, but using a linear scale for the region covering the [O I] λλ 6300, 6364 doublet and the Ca II near-infrared triplet.

10 Signatures of main-sequence mass for type II SNe and uncertainties

Figure 14 shows the model luminosity through the optical and near-infrared ranges for the type II SN models arising from ZAMS masses of 12, 15, 20, and 25 M_⊙. The model names are, in this order, m12mix100vnimH9, m15mix100vnimH9, m20mix100vnimH9, and m25mix100vnimH9, and all are characterized by weak mixing for all species except for the ⁵⁶Ni and daughter products, which are strongly mixed.

As is well known, for the same decay absorbed in the metal-rich layers, the strength of the [O I] λλ 6300, 6364 doublet increases with ZAMS mass (see for example Fransson & Chevalier 1989 or Jerkstrand et al. 2012). We obtain this result because we did not force strong mixing of elements (preventing the pollution of Ca from the Si-rich layers with the O-rich material; see Sect. 5). The systematic increase of the [O I] λλ 6300, 6364 doublet luminosity with ZAMS mass arises from the increase in decay power absorbed by the O-rich layers relative to other layers – the O-rich layers absorb a greater share of the total power available. This occurs because of the increase in CO core mass with ZAMS mass. This increase is exacerbated by the reduced O ionization with increasing progenitor mass, so that O is essentially neutral in the O-shell (in fact throughout the ejecta) in the 25 M_⊙ model. Interestingly, probably because of the strong ⁵⁶Ni mixing in this model set, Ca is twice ionized throughout the H-rich layers, which quenches Ca II line emission from the outer ejecta. But because of strong ⁵⁶Ni mixing, decay poweris preferentially absorbed outside of the Si-rich layers, which quenches Ca II line emission in the inner ejecta. Consequently, Ca II lines are weak in this model set.

One lesson to draw from this and previous sections is that the [O I] λλ 6300, 6364 doublet emission is a much more robust indicator than the [Ca II] λλ 7291, 7323 doublet because [O I]λλ 6300, 6364 predominantly arises from the O-rich shell and because the O yield strongly scales with ZAMS mass. Ca is most abundant in the low-mass Si-rich shell (whether in the progenitor star or after explosive nucleosynthesis), which absorbs little power, and the associated Ca II emission is often spread over multiple shells/regions in the ejecta. The spreading depends on the ⁵⁶Ni mixing, which is not well constrained from observations. Finally, the mass of the explosively produced Si-rich shell is sensitive to the physics and the properties of the explosion (see, for example, Woosley & Heger 2007).

11 Comparison to observations

Having produced a grid of models, we can now compare our synthetic spectra to a few well observed type II SNe in the nebular phase. We have selected a few events with standard photometric and spectroscopic properties (brightness during the photospheric phase, ejecta mass and kinetic energy, ⁵⁶Ni mass) so that they likely reside in a similar parameter space to our models. Our selection is detailed in Table 4 and includes SNe 1987A, 1999em, 2012aw, 2004et, 2013ej, and 2015bs. Figure 15 shows the comparison between observations (spanning a SN age of 300 to 420 d) to the model spectra (all at 300 d after explosion).

All type II SNe in our sample can be fitted with one or several of our 15 M_⊙ models except for SN 2015bs, which requires both a more massive progenitor and peculiar properties. The models with weak mixing (mix100) yield satisfactory matches to observations, for example for SN 2012aw. The width of the [O I] λλ 6300, 6364 doublet is well matched which suggests that the O-rich material was not extensively mixed. In the corresponding model m15mix100mH9, there is also weak ⁵⁶Ni mixing so that the bulk of the decay power is absorbed in the metal-rich core. The presence of a massive H-rich envelope in the progenitor implies that the H-rich material extends down to low velocities (composition stratification), irrespective of the weak mixing. This model also fits reasonably well SN 1987A, although some lines are too narrow or lack extended wings (the emission from the H-rich envelope is underestimated, which arises from the insufficient ⁵⁶Ni mixing). Because there is little mixing of H-rich material down to low velocities, the Hα profile generally lacks emission at low velocity (the line appears less triangular than observed).

In model m15mix100vnimH9 (used for the comparison to SN 1999em), the enhanced mixing produces broader lines overall, but the smaller power absorbed in the inner metal rich regions leads to weaker Ca II emission from the inner ejecta (the narrow part of the line) while the greater power absorbed in the H-rich layers boosts the Ca ionization and quenches CaII emission. Weak Ca II emission is rare but has been observed in SN 2012ec (Jerkstrand et al. 2015a). Because of the similar spectral appearance of SN 1999em and SN 2012aw, model m15mix100mH9 would yield a satisfactory match to SN 1999em (the goal here is not to reproduce observations but to identify trends and processes).

The broad O I and Ca II lines in SN 2013ej are hard to match with our grid of models. We find that the models with strong mixing of all species fare best. Since such a mixing may not occur in nature, this match may be artificial. It is remarkable to see that despite the very crude setup for our ejecta, most models appear rather close to observations. One intriguing example is SN 2015bs for which model m20mix100vni reproduces quite closely the large widths of [O I] λλ 6300, 6364, [Ca II] λλ 7291, 7323, and Hα. The reason behind this is the strong⁵⁶Ni mixing and the low H-rich envelope mass, allowing the abundant metal-rich material to absorb a large fraction of the decay power and expand much faster than the other SNe in the sample. In addition, the low level of mixing ensures that hydrogen is located at large velocities. Irrespective of these mixing properties, only a high mass progenitor (here a 20 M_⊙ model is used) can match these properties, as proposed by Anderson et al. (2018).

As mentioned earlier, the idea was not to provide a match to all existing type II SN spectra in the nebular phase since we craft our models. The dependencies we find are however indicative of the ejecta properties and processes that control the appearance of Type II SN spectra. Our future work will use physical models for both the progenitor and the explosion physics (as previously done, for example, in Dessart et al. 2013), as well as a more suitable technique for treating the process of mixing with greater physical consistency.

Fig. 14 Similar to Fig. 3, but now showing the optical and near-infrared spectra for type II SN models associated with different ZAMS masses from 12 to 25 M⊙. These models are all characterized by the same H-rich envelope mass of 9 M⊙ in the progenitor, by weak mixing of non-IGE species (100 km s−1), and by strong mixing of 56Ni and daughter isotopes. All models in this set have an ejecta kinetic energy of 1051 erg and a 56Ni mass of 0.08 M⊙.

12 Conclusions

In this paper, we present a grid of non-LTE steady-state calculations for type II SNe in order to examine the influence of various parameters on nebular spectra. Results from stellar evolution calculations were used to craft ejecta models. Critical properties at core collapse of stars in the mass range 12− 25M_⊙ include the monotonically increasing He-core mass with ZAMS mass or the typical abundance ratios within the main shells (i.e., the Si-rich, O-rich, He-rich, and H-rich shells). Variations were introduced in ejecta properties that may show an erratic behavior from SN to SN. This especially concerns the properties of mixing, the abundance of ⁵⁶Ni, or the mass of the H-rich envelope at the time of explosion. With this flexible approach, we examined the dependencies of nebular-phase spectra of type II SNe on variations in ejecta properties but keeping control of the robust differences seen in massive starevolution models.

We first identified a sensitivity of our results to line overlap which is enhanced when we adopt a Doppler width of 50 km s⁻¹ for all species (which is the standard procedure in our SN CMFGEN calculations). With this Doppler width some Fe II and O I lines overlap with Ly α and Ly β. Reducing the Doppler width to 2 km s⁻¹ reduces the overlap, and leads to significant changes in the resulting spectrum, and in particular to an increase in the Hα line strength. Other lines are also affected but their behavior is model-dependent. Most simulations presented in this study were run with a fixed Doppler width of 2 km s⁻¹.

The spectral properties at nebular times are complex and sensitive to numerous effects. Most of the low-energy radiation (falling primarily in the optical and near-infrared ranges) in our simulations emerges from the H-rich shell, even when the majority of the decay power is absorbed by the metal-rich layers, located in the deeper layers of the ejecta. Radiation below about 6000 Å and emitted deep in the ejecta is reprocessed by the H-rich material before it escapes. Although the conditions are nebular, there are still significant optical depth effects at 300 d after explosion.

The mixing of ⁵⁶Ni governs the spatial distribution of the decay power absorbed. Our simulations are particularly sensitive to this distribution because we retain the original progenitor shell stratification, even when mixing is applied. Thus, enhanced ⁵⁶Ni mixing favors deposition in the H-rich shell at the expense of the Si-rich material. This effect is a function of progenitor mass since the metal-rich material is less abundant (He-core mass is smaller) in a lower mass massive star. Adopting a complete mixing of the metal-rich layers would probably reduce this effect since this stratification would disappear (this aspect will be examined in a forthcoming study). Furthermore, enhanced deposition in the H-rich envelope systematically leads to an over-ionization of Ca and the weakening of the [Ca II] λλ 7291, 7323 doublet strength (this perhaps takes place in SN 2012ec; Jerkstrand et al. 2015a). This effect may serve as a means to constrain the level of ⁵⁶Ni mixing in the ejecta. Clumping may reduce this over-ionization (Dessart et al. 2018).

Varying the mass of the H-rich shell (as arises from variations in progenitor mass loss or progenitor mass) while keeping the same ejecta kinetic energy changes the chemical stratification in velocity space, the density in the metal-rich core, the absorbed decay power in the H-rich shell, or the trapping efficiency of γ-rays. This impacts primarily the H I lines in our calculations.

Varying the ⁵⁶Ni mass has far reaching consequences since it changes the luminosity, the electron density (and thus the optical depth), and the ionization. The spectral appearance is thus strongly altered. For a high ⁵⁶Ni mass, the electron density may become so large that it inhibits the formation of forbidden lines. It is therefore questionable to use SN 1987A as a template for estimating the yields or the progenitor mass of other SNe II – there is no linear scaling with ⁵⁶Ni nor with the [O I] λλ 6300, 6364 doublet flux. In the nebular phase, the impact of a varying ⁵⁶Ni mass is similar to the impact of time passing for a given SN model.

Although previously reported, for example, by Fransson & Chevalier (1989), we reemphasize the influence of the O/Ca abundance ratio in the O-rich shell. Stellar evolution calculations with MESA frequently exhibit the merging of the Si-rich and the O-rich shell during Si burning (see also Collins et al. 2018), causing the Ca abundance to rise by a factor of about 100 in the O-rich shell. Because [Ca II] λλ 7291, 7323 is a much stronger coolant than [O I] λλ 6300, 6364, this inhibits the production of [O I] λλ 6300, 6364 at nebular times, destroying any robust relationship between the [O I] λλ 6300, 6364 doublet fluxand the O-rich shell mass (and thus progenitor mass). In nature, the complete merging of the Si-rich and O-rich shells is probably unlikely, but some contamination of the O-rich shell by Ca-rich material from the Si-rich shell is a possibility (caused by convection and overshoot) and could introduce some variations in the [O I] λλ 6300, 6364 doublet strength for a given progenitor composition. If the process of shell merging had a higher occurrence rate in higher mass progenitors (owing to the more violent convection that takes place in their interiors), it would provide a natural explanation for the lack of SNe II with strong [O I] λλ 6300, 6364 emission at nebular times, and inferred progenitor masses below 17 M_⊙.

Because we simplified the composition and the associated model atoms, we will have missed some coolants and thus may overestimate the cooling power of some of the lines that we do include. For example, our simulations show a systematic increase in the [O I] λλ 6300, 6364 doublet luminosity with ZAMS mass, largely irrespective of the adopted mixing, from 2% (model m12), to 5% (model m15), 12% (model m20), and ~ 16% (model m25) of the bolometric luminosity (or ⁵⁶Co decay power absorbed). This is typically a factor of two greater than obtained by Jerkstrand et al. (2015a).

Overall, our simulations suggest that the [O I] λλ 6300, 6364 doublet strength is the most robust indicator of progenitor mass. The O-rich shell is the most massive metal-rich shell in 12-25 M_⊙ progenitors. Its mass grows considerably with ZAMS mass and thus its associated material captures a large fraction of the decay power. Furthermore, oxygen is generally neutral under a wide range of ejecta conditions (mixing, ⁵⁶Ni mass, progenitor mass), while we obtain a very complicated behavior for Ca (strength dependent on ionization, ⁵⁶Ni mixing, line optical depth, etc.). We must however tone down this conclusion since the potential Ca pollution from the Si-rich shell into the O-rich shell can mitigate the [O I] λλ 6300, 6364 doublet flux at nebular times, making a type II SN appear as if it arose from a lower mass star without such Ca pollution. The implication of this process for progenitor mass estimates is profound and currently understated.

Fig. 15 Comparison of optical spectra for a sample of observations and a selection of models presented in this work – we make no attempt at obtaining best fits. Apart from SN 2015bs, which stands apart, most of these well-observed SNe II show similar spectral properties at 300− 400 d after explosion. A variety of models is used to show the offsets caused by certain assumptions (for example, the weak [Ca II] λλ 7291, 7323 that follows from Ca over-ionization caused by strong 56Ni mixing in model m15mix100vnimH9). Model m15mix100mH9, used for SN 2012aw, would provide a better match to the observationsof SN 1999em.

Acknowledgements

We thank Bill Paxton and Frank Timmes for their help with MESA, specifically in setting up the proper controls that prevent shell merging during Si burning. Hillier thanks NASA for partial support through the astrophysical theory grant 80NSSC20K0524. We thank the referee, Anders Jerkstrand, for their thoughtful comments on an earlier version of this paper. This work was granted access to the HPC resources of CINES under the allocation 2018 – A0050410554 and 2019 – A0070410554 made by GENCI, France. This research has made use of NASA’s Astrophysics Data System Bibliographic Services.

References

All Tables

All Figures

	Fig. 1 Cumulative composition and velocity profiles versus Lagrangian mass in the ejecta of the SN II model based on the 20 M⊙ ZAMS star. The sum for each plotted element i includes the mass fractions for all plotted elements that have a lower atomic mass.
In the text

Fig. 2 Influence of the adopted Doppler width on the resulting optical and near-infrared spectra for a SN II model arising from a 15 M⊙ progenitor. A weak mixing of all species is used except for 56Ni for which we adopt a strong mixing. The increase in Hα line flux with decreasing VDop is at the expense of a decrease in flux from a forest of Fe II lines in the 4000− 6000Å region. Other lines affected are He I 1.083 μm, O I 1.129 μm, Fe I lines around 1.165 μm, Mg I lines at1.183, 1.502, 1.577, 1.711, H I 1.875 μm. Note that about 70% of the total flux falls in the optical range.

In the text

Fig. 3 Illustration of properties for the reference model m20mix100vnimH9 discussed in Sect. 4. In this and similar figures, the top panels show the luminosity Lλ (integrating over wavelength yields the bolometric luminosity) in the optical (upper) and the near infrared (lower). The line identifications are indicative only since in numerous cases multiple lines contribute (we give the primary component to most features). Forbidden lines are shown in red for all ions and atoms apart from Fe (shown in blue). Bottom panels: some ejecta properties computed with CMFGEN.

In the text

Fig. 4 Illustration of the spatial regions (here shown in velocity space) contributing to the emergent flux in a 20 M⊙ type II SN model characterized by weak mixing for all non-IGE species, strong mixing for IGE species, and with a 9 M⊙ H-rich envelope mass (model m20mix100vnimH9). Top panel: observer’s frame luminosity contribution δLλ,R (Eq. (5); the map maximum is saturated at 20% of the true maximum to bias against the strong Hα line and better reveal the origin of the weaker emission) versus wavelength and ejecta velocity. The four contributing shells in this type II SN model are clearly seen (see vertical colored stripe at left), although the H-rich layers contribute most of the emergent flux. Bottom panel: emergent luminosity integrated over each shell, together with the total luminosity. In this model, the fraction of the total power emerging from the Si-rich, O-rich, He-rich and H-rich shells is 4, 20, 15, and 61%, respectively. This illustrates the contribution to the emergent total flux from the main ejecta shells (this may also be inferred from the width of lines but complicated when multiple regions contribute or when there is line overlap). Note that with our approach, IGEs are present in all shells and thus Fe I and Fe II flux contribution is present throughout the ejecta. One can also see that the flux from the Si-rich shell is mostly radiated by Ca I and Ca II (together with Fe I and Fe II, rather than by the dominant species Si).

In the text

	Fig. 5 Illustration of the total (i.e., accounting for all opacity sources) radial optical depth integrated from the outer ejectaboundary to the innermost boundary (i.e., the base of the ejecta; blue), to the He/O shell interface at 1780 km s−1 (red), and to the He/H shell interface at 2150 km s−1 (yellow) for the reference model m20mix100vnimH9.
In the text

Fig. 6 Illustration of the main cooling processes balancing the radioactive decay heating at all ejecta depths in the reference model m20mix100vnimH9. 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-thermalexcitation and ionization) and “COL” stands for collisional processes (i.e., collisional excitation).

In the text

Fig. 7 Similar to Fig. 3, but now showing variants of model m20mix100mH9 in which the Ca mass fraction is progressively increased in the O-rich shell by powers of ten and starting from about 7 × 10−5, the solar metallicity value (the Ca mass fraction profile is shown in the bottom left panel). Notice the dramatic reduction on the [O I]/[Ca II] line ratio as the Ca abundance is increased. Even in the richest Ca model, the Ca mass is still less than 10% of the O mass.

In the text

	Fig. 8 Illustration of the Sobolev line optical depth for the Ca II 7323.9 Å line for the models shown in Fig. 7.
In the text

Fig. 9 Similar to Fig. 3, but now showing the influence of the progenitor H-rich envelope mass on the synthetic optical and near-infrared spectra at 300 d. The models shown are m25mix100vni and m25mix100vnimH9 (strong mixing of 56Ni but weak mixing of other species) and correspond to a 25 M⊙ progenitor star. There is a greater fraction of the total decay power absorbed in the H-rich layers of the lighter model m25mix100vni, boosting the [Ca II] λλ 7291, 7323 doublet emission, but this has an adverse effect on Hα, which pushes its formation to the outer lower density regions. These differences reflect in part the different distribution of decay-power absorbed since the two ejecta have a different velocity structure – mixing should also differ in nature for ejecta with a similar kinetic energy but a different mass (e.g., as in SNe II and Ibc).

In the text

Fig. 10 Similar to Fig. 3, but now showing the influence of the 56Ni mixing on the SN radiation and gas properties at 300 d. The ejecta arises from a 20 M⊙ progenitor, weak mixing for all species, imposed H-rich envelope of 9 M⊙, and a 56Ni mass of 0.08 M⊙. An additional and variable amount of 56Ni mixing is applied increasing VNi from 750 (model 20-vni7p5e7) to 2500 km s−1 (model 20-vni2p5e8). With enhanced 56Ni mixing, more decay power is absorbed in the outer ejecta, causing a strengthening of Hα but a weakening of [Ca II] λλ 7291, 7323 because of the boost to the Ca ionization (the ionization shift also affects Fe I – II and Mg I lines in the optical and the near-infrared). (See Sect. 7 for discussion.)

In the text

Fig. 11 Similar to Fig. 3, but now showing the influence of the 56Ni mass on the SN radiation (for better visibility, a logarithmic scale is used for the luminosity). The ejecta arises from a 20 M⊙ progenitor, weak mixing for all species, no additional mixing of 56Ni, imposed H-rich envelope of 9 M⊙, and is characterized by a 56Ni mass of 0.004, 0.02, 0.08, 0.16, and 0.32 M⊙. A greater 56Ni mass boosts the electron density and species ionization, inhibits forbidden line emission (in favor instead, for example, of recombination lines), and strengthens optical depth effects and continuum processes (primarily in association with electron scattering). (See Sect. 8 for discussion.)

In the text

Fig. 12 Similar to Fig. 3, but now showing the evolution of the optical and near-infrared properties of model m20mix100vnimH9 from 150 to 500 d after explosion. With increasing time, the conditions become increasingly “nebular”, with a relative strengthening of forbidden line emission, and the progressive reduction of optical depth effects (for example with the disappearance of P-Cygni profiles). This evolution is similar to reducing the 56Ni mass at a given SN age (see Fig. 11).

In the text

	Fig. 13 Asfor Fig. 12, but using a linear scale for the region covering the [O I] λλ 6300, 6364 doublet and the Ca II near-infrared triplet.
In the text

Fig. 14 Similar to Fig. 3, but now showing the optical and near-infrared spectra for type II SN models associated with different ZAMS masses from 12 to 25 M⊙. These models are all characterized by the same H-rich envelope mass of 9 M⊙ in the progenitor, by weak mixing of non-IGE species (100 km s−1), and by strong mixing of 56Ni and daughter isotopes. All models in this set have an ejecta kinetic energy of 1051 erg and a 56Ni mass of 0.08 M⊙.

In the text

Fig. 15 Comparison of optical spectra for a sample of observations and a selection of models presented in this work – we make no attempt at obtaining best fits. Apart from SN 2015bs, which stands apart, most of these well-observed SNe II show similar spectral properties at 300− 400 d after explosion. A variety of models is used to show the offsets caused by certain assumptions (for example, the weak [Ca II] λλ 7291, 7323 that follows from Ca over-ionization caused by strong 56Ni mixing in model m15mix100vnimH9). Model m15mix100mH9, used for SN 2012aw, would provide a better match to the observationsof SN 1999em.

In the text