Radiative-transfer modeling of nebular-phase type II supernovae. Dependencies on progenitor and explosion properties

Nebular phase spectra of core-collapse supernovae (SNe) provide critical and unique information on the progenitor massive star and its explosion. We present a set of 1-D steady-state non-local thermodynamic equilibrium radiative transfer calculations of type II SNe at 300d after explosion. Guided by results for a large set of stellar evolution simulations, we craft ejecta models for type II SNe from the explosion of a 12, 15, 20, and 25Msun star. The ejecta density structure and kinetic energy, the 56Ni mass, and the level of chemical mixing are parametrized. Our model spectra are sensitive to the adopted line Doppler width, a phenomenon we associate with the overlap of FeII and OI lines with Lyalpha and Lybeta. Our spectra show a strong sensitivity to 56Ni mixing since it determines where decay power is absorbed. Even at 300d after explosion, the H-rich layers reprocess the radiation from the inner metal rich layers. In a given progenitor model, variations in 56Ni mass and distribution impact the ejecta ionization, which can modulate the strength of all lines. Such ionization shifts can quench CaII line emission. In our set of models, the OI6300 doublet strength is the most robust signature of progenitor mass. However, we emphasize that convective shell merging in the progenitor massive star interior can pollute the O-rich shell with Ca, which will weaken the OI6300 doublet flux in the resulting nebular SN II spectrum. This process may occur in Nature, with a greater occurrence in higher mass progenitors, and may explain in part the preponderance of progenitor masses below 17Msun inferred from nebular spectra.


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 56 Co decay heating. 1 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 56 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 1 In this study we use the electron energy balance (e.g., Osterbrock 1989;Hillier & Miller 1998) to discuss heating and cooling processes. In this equation (and in addition to the usual terms) the energy absorbed from radioactive decays appears as a heating term, while non-thermal excitation and ionization appear as coolant terms. Emission in Hα does not appear as an explicit coolant. It is primarily produced via recombination but its strength is implicitly linked to non-thermal processes through the coupling between the electron energy balance equation and the rate equations. In practice a significant fraction of the absorbed decay energy is emitted as Hα. 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 andMcCray 1993). Numerous studies were focused on the nebular-phase radiation properties of SN 1987A (see, for example, Fransson & Chevalier 1987;Kozma & Fransson 1992;Li & McCray 1992;Li & McCray 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 56 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 largescale 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 nebularphase modeling has been presented in a series of papers by Jerkstrand et al. (2014Jerkstrand et al. ( , 2012Jerkstrand et al. ( , 2015b. 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 Article number, page 1 of 23 arXiv:2007.02243v1 [astro-ph.SR] 5 Jul 2020 A&A proofs: manuscript no. ms handful of objects (see, for example, the dataset presented by Silverman et al. 2017).
An important conclusion 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 − 25 M progenitors (e.g., Jerkstrand et al. 2015b; see also Smartt 2009). A notable exception is the type II SN 2015bs, for which a high mass progenitor of 20 − 25 M seems a plausible explanation for its unprecedented large [O i] λλ 6300, 6364 doublet flux (Anderson et al. 2018).
With CMFGEN , we have conducted numerous simulations for type II SNe (Dessart & Hillier 2011;Dessart et al. 2013;Lisakov et al. 2017;Dessart et al. 2018;Dessart & Hillier 2019a) 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 andSilverman 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-thermal processes 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 introduce two different levels of mixing in our simulations. We may for example apply strong mixing of 56 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 56 Ni, 56 Co, and 56 Fe throughout the ejecta, but of these three species, Fe dominates at nebular times of 300 d. 2 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 Section 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 section 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 −1 . In Section 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 56 Ni mixing and of the 56 Ni mass on our nebular-phase spectra in sections 7 and 8. While all previous simulations were performed at a SN age of 300 d, section 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 56 Ni mass and ejecta kinetic energy) on the nebular-phase spectra in section 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 section 12.

Context and numerical setup
2.1. Properties of massive star progenitors at core collapse: results from a grid of MESA models We have performed stellar evolution simulations with MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 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 −1 ). 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 section 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 want 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 sections 5 and 12). Hence, to prevent the merging of the Sirich 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 28 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 Sirich 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). Convection is much more vigorous in 3D than 1D, and this leads to greater mixing.
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 shell and deeper metalrich 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 4 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, section 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.

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 in the 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 − 25 M 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 star evolve 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 slight shifts 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 51 erg and a fixed 56 Ni mass of 0.08 M for all models unless otherwise stated (section 8 describes the influence of varying the 56 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 56 Ni profile (function of the adopted mixing; see below). At nebular epochs, 56 Ni has decayed into a mix of 56 Fe, 56 Co, and 56 Ni that have the same distribution in velocity space and a cumulative mass equal to the initial 56 Ni mass ( 58 Ni and 59 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 Table 1. Shell masses and abundance ratios in our grid of massive star models. The quantity {X/Y} corresponds to the total mass of element X over the total mass of element Y in the shell for the corresponding column. Numbers in parenthesis correspond to powers of ten. 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 Section 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 51 erg. We primarily focus on one epoch, 300 d after the explosion. However, we also investigate the evolution with time in section 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 where ρ 0 and V 0 are constrained by the adopted ejecta kinetic energy E kin , the ejecta mass M ej , and the density exponent k (set to eight) through and where C m = π k sin(3π/k) ; C e = π k sin(5π/k) .
(3) Notes: a The type II SN model m12 has a total mass greater than the initial mass by 0.4 M because of the adopted H-rich envelope mass. Its core properties are, however, compatible with the results from our MESA grid. Such a model could arise from a merger (Menon & Heger 2017).
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 56 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 −1 ) or strong mixing (gaussian width of 400 km s −1 ) were implemented. For greater freedom, the mixing of 56 Ni and other species is chosen to be independent. The adopted initial profile for 56 Ni is of the form and connects continuously to a constant 56 Ni mass fraction for V < V Ni . The normalization is set by the specified 56 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 56 Ni and its daughter elements, 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.
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;. 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  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 Li et al. 2012). We limit the radioactive decay to the 56 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 2 g −1 , where Y e is the electron fraction. The model atoms included are: H i 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 temperature and electron density profile takes about 200 iterations, but only 12h of computing time for a Doppler width of 50 km s −1 (see section 3). However, once a given model is converged, variants of that model (for example the same ejecta model but with a different 56 Ni mass or metallicity, etc.) can be computed quickly by adopting this converged model as an initial guess for the temperature, level populations, etc.

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 56 Ni-rich shell. In practice, this simplification amounts to introducing a 56 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 56 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 56 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 preparation).
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, the s-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 10000 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.

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 56 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 56 Ni or sometimes both; the mass of the H-rich envelope in the progenitor; the mass of 56 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.

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 3 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 −1 at the low temperatures of SN ejecta at nebular times. This velocity width scales with 1/ √ A, 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 −1 . This choice is motivated by speed, and the need to have a "reasonable" number of frequency points (still of order 10 5 ) 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 −1 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 −1 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 −1 is used for all species (in cases where the impact on the gas properties was large, this reduction was done in a few steps). 4 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 −1 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 −1 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 nearinfrared Mg i lines) which may have resulted from the increase of the electron density and the H ionization in the recombined Hrich 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 1 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  changes by running test calculations which used a Doppler width of 50 km s −1 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 β.

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 −1 (this mixing is applied to all species), "vni" means that the additional mixing applied to 56 Ni was strong (V Ni = 2500 km s −1 and ∆V Ni = 1000 km s −1 ; see Eq. 4), and "mH9" means that the pre-SN progenitor had a 9 M H-rich outer shell. The ejecta has a kinetic energy of 10 51 erg and contains 0.08 M of 56 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 −1 Å −1 ) 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 56 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 lowenergy, mostly optical, photons -whatever is absorbed is instantaneously radiated and the balance is zero. Because of the strong 56 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 −1 get ∼ 20 %, the O-rich layers between 1000 and 1800 km s −1 get ∼ 40 %, the He-rich layers between 1800 and 2100 km s −1 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: 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,  The bottom panel shows the 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).
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 Herich 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) 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 large peak value at any depth) normalized to the local heating rate. The term "NT" stands for non-thermal processes (in this context non-thermal excitation and ionization) and "COL" stands for collisional processes (i.e., collisional excitation). 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 −1 so the variations of the electron density below 2000 km s −1 reflect the change in ionization and mean atomic weight as we proceed through the He-rich shell, the Orich 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 56 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 56 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 56 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.

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 power emitted is absorbed by the ejecta. Hence, a fundamental aspect of nebular line emission is how the 56 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 Article number, page 10 of 23 Luc Dessart and D. John Hillier: Modeling of nebular phase spectra for type II SNe Let us consider a small gas volume of uniform composition with a gas temperature T (or T 4 when expressed in units of 10 4 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 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). 5 At the critical density, this rate is equal to N u A ul . Thus N e (crit) = 1.16 × 10 7 g u A ul 5 The values listed in the table are those used in the present calculations. Experimental lifetime measurements, and theoretical calculations, suggest that a better estimate for A for the Ca ii transitions is 1.1 (see Meléndez et al. 2007). The value listed in the NIST database is 1.3 (Kramida & NIST ASD Team 2019). It comes from a calculation by Osterbrock (1951) and has an indicated error of greater than 50%. The O i values are slightly lower than those in the NIST database Kramida & NIST ASD Team (2019) (5.6 × 10 −3 and 1.8 × 10 −3 ; E < 7%) however Baluja et al. (1988) suggest even higher values (6.7 × 10 −3 and 2.3 × 10 −3 ). While the adoption of different cross-sections will (slightly) alter line strengths, they will not affect any of the conclusions made in this paper.
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 u A 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 or Similarly, for the Ca ii line (g u /g l = 10/2), we have η(Ca ii) = 5 hν ul A ul N Ca exp −1.4388 λ(µm)T 4 .
Integrating these line emissivities over the O-rich shell volume yields the line luminosities L. Taking the ratio we obtain, Case 2: N e N e (crit). In this case, the line emissivity is since every excitation gives a line photon. Thus, the ratio of line luminosities is now In this case we could add a factor like (T 4 /0.5) 0.5 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.
Article number, page 11 of 23 A&A proofs: manuscript no. ms 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). The shell masses grow with the progenitor mass, with a 2 − 3 times faster increase for the O-rich shell mass.

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 section 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 −5 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).

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 Orich 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 is the 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 section 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 1996or 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 occurs 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.  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 Sirich 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 Sirich 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 56 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 .

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 might be 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 56 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.

Influence of 56 Ni mixing
In this section, we explore the influence of 56 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 56 Ni with V Ni covering from 750 to 2500 km s −1 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 56 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 de-posited below 1000 km s −1 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 line emission. 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 56 Ni mixing is enhanced, it first benefits the Orich shell which produces a stronger [O i] λλ 6300, 6364 doublet. As the 56 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 iii 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 56 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 2+ everywhere beyond 1000 km s −1 for the model with the highest mixing. 6 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 56 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 56 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 56 Ni mixing is a complicated process, with 56 Ni probably having a unique distribution in each SN. This is most likely tied to additional dependencies on progenitor mass and rotation, Hecore 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  that [Ca ii] λλ 7291, 7323 arises from the H-rich material (which they infer from their analysis of SN 1987A) only holds if 56 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 56 Ni and thus absorbs decay power in the absence of strong 56 Ni mixing). A further issue is that Ca overionization can completely inhibit [Ca ii] λλ 7291, 7323 emission from the H-rich layers, no matter how much decay power is absorbed there.

Dependencies on 56 Ni mass
We now discuss the influence of the 56 Ni mass on the resulting optical and near-infrared spectra as well as on the  . 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 56 Ni 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).    ejecta properties. For this exploration, we start off with model m20mix100mH9 (20 M progenitor, weak mixing for all species, no additional mixing of 56 Ni, imposed H-rich envelope of 9 M ), characterized by a 56 Ni mass of 0.08 M and compute additional models in which the 56 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 56 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 56 Ni mass. There are various reasons for this.
For higher 56 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 56 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 56 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 56 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 56 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 56 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 56 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 56 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. 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 3 ) and ionization level (as material cools and recombines) leads to a drop in electron density.

Evolution during the nebular phase
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 result at 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 section 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 . 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). 7

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 56 Ni and daughter products, which are strongly mixed. As is well known, for the same decay absorbed in the metalrich layers, the strength of the [O i] λλ 6300, 6364 doublet increases with ZAMS mass (see for example Fransson & Chevalier 1989or 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 section 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 56 Ni mixing in this model set, Ca is twice ionized throughout the H-rich layers, which quenches Ca ii line emission from the outer  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 56 Ni mass at a given SN age (see Fig. 11).  ejecta. But because of strong 56 Ni mixing, decay power is 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 56 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).

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, 56 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 15M 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 56 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 4000  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 56 Ni and daughter isotopes. All models in this set have an ejecta kinetic energy of 10 51 erg and a 56 Ni mass of 0.08 M . arises from the insufficient 56 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 Ca ii emission. Weak Ca ii emission is rare but has been observed in SN 2012ec (Jerkstrand et al. 2015b). 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 56 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 Table 4. Characteristics of the observed Type II SNe used in this paper, including the inferred time of explosion, the redshift, the distance, the reddening, and the reference from whence these quantities and observational data were taken. 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. 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 56 Ni mixing in model m15mix100vnimH9). Model m15mix100mH9, used for SN 2012aw, would provide a better match to the observations of SN 1999em.

Conclusions
We have presented 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 − 25 M 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 concerns especially the properties of mixing, the abundance of 56 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 star evolution 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 −1 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 −1 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 −1 .
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 Hrich material before it escapes. Although the conditions are nebular, there are still significant optical depth effects at 300 d after explosion.
The mixing of 56 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 56 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 overionization of Ca and the weakening of the [Ca ii] λλ 7291, 7323 doublet strength (this perhaps takes place in SN 2012ec; Jerkstrand et al. 2015b). This effect may serve as a means to constrain the level of 56 Ni mixing in the ejecta. Clumping may reduce this over-ionization .
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