How Do Type Ia Supernova Nebular Spectra Depend on Explosion Properties? Insights from Systematic Non-LTE Modeling

The underlying SN ejecta model is specified by the mass density and elemental abundances in each zone on a 3D Cartesian grid. We assume that the ejecta are in homologous expansion (i.e., velocity proportional to radius), which is appropriate for SNe Ia at a few seconds to days after explosion (e.g., Röpke 2005). The structure of a homologous model at one epoch can be easily scaled to any other time, taking into account compositional changes due to radioactive decay. The models in this paper study SNe Ia out to 400 days and include the following species: C, O, Si, S, Ca, Fe, Co, and Ni (see Appendix A for a description of atomic data sources).

Our calculations assume stationarity—i.e., that the gas temperature and level populations reach equilibrium on a timescale short compared to the ejecta expansion timescale. This assumption is reasonable except at rather late times (≳500 days), when thermal and ionization freeze-out may become important (Fransson & Kozma 1993; Kozma & Fransson 1998a; Sollerman et al. 2004; Fransson & Jerkstrand 2015).

Our transport solver assumes that the ejecta are optically thin to radiation. For many SNe, this is true in the optical and infrared regions at times ≳100 days. However, at blue and ultraviolet wavelengths (≲4000 Å), the ejecta may remain opaque to iron-group lines for hundreds of days (Friesen et al. 2017). This limits the reliability of our nebular models at short wavelengths and may affect the calculated ionization fractions (e.g., Sollerman et al. 2004), issues that will be addressed in the future by incorporating a more general transport solver.

To calculate spectral emission from the SN nebula, one has to determine the level populations of each ion in the gas. Since LTE is a poor approximation at these epochs, we solve for each species the set of NLTE equations expressing statistical equilibrium

$$\begin{eqnarray}&&{\boldsymbol{Mn}}=0,\end{eqnarray} \tag{ 1 }$$

where ${\boldsymbol{n}}$ is a vector of level populations for the species and the matrix ${\boldsymbol{M}}$ encodes the transition rates between the various levels and ionization states (see, e.g., Li & McCray 1993). To the set of equations must be added the constraint of number conservation.

A generalized form for the statistical equilibrium rate equations is

$$\begin{eqnarray}&&{n}_{I,i}\left({R}_{I,i}+\displaystyle \sum _{j\ \ne \ i}{{ \mathcal E }}_{{ij}}\right)=\displaystyle \sum _{k\ \ne \ i}{n}_{k}{{ \mathcal E }}_{{ki}}+{n}_{I+1,{gs}}{n}_{e}{\alpha }_{i},\end{eqnarray} \tag{ 2 }$$

where n_I,i is the level population of the ith level of the ion with Ith ionization state, n_e is the electron density, R_I,i is the total ionization rate from the ith to the I + 1th ionization state (including photoionization, collisional ionization, and nonthermal ionization), and α_i(T) is the total recombination coefficient from ionization state I's ground state to the ith level (including radiative, dielectronic, and three-body recombination). ${ \mathcal E }$ are bound–bound rate coefficients between levels i and j, including spontaneous emission, stimulated absorption/emission, and collisional excitation/de-excitation.

The microphysics relevant to the nebular phase calculation is encoded into the transition matrix of Equation (1). For the SN Ia problem, we treat collisional ("electron-impact") bound–bound and bound–free processes, radiative and dielectronic recombination, and nonthermal deposition of ⁵⁶Co energy into heating, excitation, and ionization channels. Following Nahar & Pradhan (1997), we assume that all collisional ionization rates are ground-state rates and that radiative recombinations come from the ground state. Details of the implementation of these processes are given in Appendix C.

Energy deposition due to the radioactive decay chain ⁵⁶Ni $\to$⁵⁶Co $\to$⁵⁶Fe is the dominant energy source for SNe Ia in the nebular phase (Colgate & McKee 1969; Kuchner et al. 1994). Since ⁵⁶Ni decays on a timescale of ${\tau }_{\mathrm{Ni}}=7.7$ days, its contribution to heating is minor by nebular times, while ⁵⁶Co decay (${\tau }_{\mathrm{Co}}=111.3$ days) is generally the most important source of radioactive energy. Other radioactive isotopes with half-lives in the nebular range (e.g., ²²Na, ³⁵S, ⁴⁵Ca, ⁴⁶Sc, ⁴⁹V, ⁵⁴Mn, ⁵⁵Fe, ⁵⁷Co, ⁵⁸Co, ⁶⁵Zn, and ⁶⁸Ge) are usually not produced in enough abundance to significantly contribute to heating at the epochs we consider here (Iwamoto et al. 1999; Seitenzahl et al. 2013) and are not included in the present calculations.

We determine the radioactive energy deposition rate in each zone using a 3D Monte Carlo radiation transport scheme (Kasen et al. 2006) that samples gamma-ray wavelengths from the radioactive decay lines and treats energy losses due to Compton scattering and photoionization. Positrons from ⁵⁶Co decay are assumed to be trapped locally, in accordance with recent observations of SNe Ia at late times (Leloudas et al. 2009; Kerzendorf et al. 2014; but see Dimitriadis et al. 2017).

The gamma rays from ⁵⁶Co decay produce high-energy electrons (${E}_{0}\sim 1$ MeV) that interact with the ejecta through heating (Coulomb interactions), ionization, and excitation. These interactions affect the temperature and ionization state of the gas. The rates of excitation and ionization, respectively, by nonthermal electrons/positrons are given by

$$\begin{eqnarray}&&{R}_{\mathrm{ex}}=\displaystyle \frac{{\dot{\epsilon }}_{\mathrm{rad}}{\eta }_{{ij}}}{{\rm{\Delta }}{E}_{{ij}}{n}_{i}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{R}_{\mathrm{ion}}=\displaystyle \frac{{\dot{\epsilon }}_{\mathrm{rad}}{\eta }_{{\rm{k}}}}{{I}_{{\rm{k}}}{n}_{{\rm{k}}}},\end{eqnarray} \tag{ 4 }$$

where the transition from level i to level j (with $j\gt i$) has energy ΔE_ij, n_i is the level population of theith atomic level, ${I}_{{\rm{k}}}$ is the ionization potential of the ion indexed by k having number density ${n}_{{\rm{k}}}$, ${\dot{\epsilon }}_{\mathrm{rad}}$ is the radioactive energy deposition rate per unit volume, and η_ij and ${\eta }_{{\rm{k}}}$ refer to the nonthermal excitation and ionization deposition fractions, respectively, described in Appendix B.

We use an iterative nonlinear solver to determine the free electron density and temperature in each zone. The electron density is constrained by the charge conservation condition that all free electrons come from ionization, i.e.,

$$\begin{eqnarray}&&{n}_{e}=\displaystyle \sum _{{\rm{k}},i}{n}_{{\rm{k}},i}i,\end{eqnarray} \tag{ 5 }$$

where n_k,i represents the population of the ith ionization stage of the kth species (i = 0 is neutral). Given that the nonthermal deposition fractions (η_ij and ${\eta }_{{\rm{k}}}$ in Equations (3)–(4)) are level population dependent, we recalculate them during each iterative step until convergence is reached.

The temperature of each zone is calculated from the balance of heating from radioactive decay and cooling due to line emission. The line emissivity (units erg s⁻¹ cm⁻³) from a transition between the ith and jth levels of an ion (${E}_{j}\gt {E}_{i}$) is

$$\begin{eqnarray}&&{\dot{\epsilon }}_{{ij}}(T)=h{\nu }_{{ij}}{A}_{{ij}}{n}_{j},\end{eqnarray} \tag{ 6 }$$

where ν_ij is the photon frequency of the transition, A_ij is the spontaneous radiative decay rate (Einstein A coefficient), and n_j is the population of the jth level of the ion. The temperature sensitivity of ${\dot{\epsilon }}_{{ij}}$ arises primarily from the temperature-dependent collisional excitation rates, which are important in setting the populations ${n}_{j}$ for low-lying states.

The total emissivity per unit volume is then a sum over all line transitions, which, in equilibrium, is equal to the energy deposition by radioactivity

$$\begin{eqnarray}&&\displaystyle \sum _{i,j}{\dot{\epsilon }}_{{ij}}={\dot{\epsilon }}_{\mathrm{rad}},\end{eqnarray} \tag{ 7 }$$

where {i, j} runs over all transitions. Note that Equation (7) represents a balance of all emission and deposition processes, not just the thermal processes.

We use an iterative Brent–Dekker method to solve the nonlinear Equations (5) and (7) for the temperature and electron density in each zone. For each iteration, we solve the NLTE rate equations to determine the ionization/excitation state and emissivity. The values of T and n_e are then adjusted and the procedure iterated until thermal and ionization equilibrium is reached. A single zone typically converges within 200–400 iterative steps to reach an accuracy of 0.1% in temperature.

Once the line emissivities (Equation (6)) are known for all transitions in a converged model, we can integrate the emission for any arbitrary geometry to determine the nebular spectrum of the SN. Assuming homologous expansion (r = vt) and choosing the direction to the observer to be the z-axis, the nonrelativistic Doppler effect along the line of sight is

$$\begin{eqnarray}&&{\lambda }_{\mathrm{obs}}={\lambda }_{{ij}}\left(1-\displaystyle \frac{z}{{ct}}\right),\end{eqnarray} \tag{ 8 }$$

where ${\lambda }_{\mathrm{obs}}$ is the observed wavelength of a line transition with rest-frame wavelength λ_ij. Equation (8) allows us to associate the observed flux at wavelength ${\lambda }_{\mathrm{obs}}$ with the line emission integrated over a specific plane (perpendicular to the z-axis) sliced through the ejecta. For isolated transitions, this mapping can be used to constrain the geometrical distribution of material along the line of sight (see, e.g., Shivvers et al. 2013). Since the bulk of the emission typically comes from ejecta with velocities ≲10,000 km s⁻¹, using the nonrelativistic Doppler shift leads to errors in z-mapping of only ≲3%.

Assuming that the ejecta are optically thin, the observed specific intensity can be expressed as the integral

$$\begin{eqnarray}&&{L}_{\lambda }(\lambda )d\lambda =\left[\displaystyle \sum _{i,j}\displaystyle \frac{{ct}}{{\lambda }_{{ij}}}\int {\dot{\epsilon }}_{{ij}}(x,y,z){dxdy}\right]d\lambda ,\end{eqnarray} \tag{ 9 }$$

where the sum {i, j} runs over all possible transitions and the integration is over the x–y plane at location $z={ct}({\lambda }_{{ij}}-\lambda )/{\lambda }_{{ij}}$. The spectral flux observed at Earth is simply ${F}_{\lambda }={L}_{\lambda }/(4\pi {D}^{2})$, where D is the distance to the source.

Dust formation might be possible in SNe Ia at late times, but not necessarily in significant amounts (Nozawa et al. 2011). Further, observations of the nearby SN 2011fe show no evidence for dust formation at 930 days past maximum (Kerzendorf et al. 2014). We therefore neglect dust formation and plan on implementing it for future work.

While implementing photoionization and stimulated radiative processes in our code, we neglect them for the work published in this paper under the assumption that there is negligible continuum radiation field.

Charge transfer (CT) is expected to occur in SNe between neutral atoms and ions (Swartz 1994). While CT might affect ionization fractions of SNe Ia at later times, we tested CT between ions of Fe and found that their effect on the nebular spectrum at 200–400 days is negligible. This is because the nebula is primarily ionized gas, and therefore the rate of ionization/recombination due to CT is subdominant. Since CT preferentially ionizes neutral atoms (as Coulomb repulsion suppresses ion–ion interactions), we expect that neutral atoms like O i and Fe i should not contribute much to nebular emission. We plan to add a detailed CT treatment to the code in the near future in order to quantify this effect.

We ran a number of tests to verify the code. We tested the NLTE level population solver in the limit that collisional processes dominate and the limit that the radiation field is the Planck function, both resulting in the expected LTE level populations. We also tested the ionization solution in the collisional ionization equilibrium (CIE) regime and found good agreement with previous results (e.g., Sutherland & Dopita 1993). We further compared total CIE cooling functions for individual ions to published results from the Cloudy code (Gnat & Ferland 2012). Given disparities in the atomic data inputs, exact agreement of the cooling functions is not expected, but we found reasonably similar values and temperature dependences in the regimes of interest. The radiation transport calculation of synthetic spectra was verified by comparing single line profiles to analytic solutions.

We model the ejecta with a broken power-law density profile that is shallow in the core and steep in the outer layers (Chevalier & Soker 1989; Kasen 2010):

$$\begin{eqnarray}\rho (v)=\left\{\begin{array}{ll}{\rho }_{0}{\left(\tfrac{v}{{v}_{{\rm{t}}}}\right)}^{-\delta } & \,v\leqslant {v}_{{\rm{t}}}\\ {\rho }_{0}{\left(\tfrac{v}{{v}_{{\rm{t}}}}\right)}^{-{\rm{n}}} & \,v\gt {v}_{{\rm{t}}}\end{array}\right.,\end{eqnarray} \tag{ 10 }$$

where ${\rho }_{0}$ can be interpreted as the central density of a perfectly flat core profile ($\delta =0$) and ${v}_{{\rm{t}}}$ is the transition velocity marking the interface of the two regions. Integration gives (assuming $\delta \lt 3$ and $n\gt 3$)

$$\begin{eqnarray}&&{\rho }_{0}=\displaystyle \frac{{M}_{\mathrm{ej}}}{4\pi {({v}_{{\rm{t}}}{t}_{\mathrm{ex}})}^{3}}{\left[\displaystyle \frac{1}{3-\delta }+\displaystyle \frac{1}{n-3}\right]}^{-1},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{E}_{{\rm{K}}}=\displaystyle \frac{1}{2}{M}_{\mathrm{ej}}{v}_{{\rm{t}}}^{2}\left[\displaystyle \frac{1}{5-\delta }+\displaystyle \frac{1}{n-5}\right]{\left[\displaystyle \frac{1}{3-\delta }+\displaystyle \frac{1}{n-3}\right]}^{-1},\end{eqnarray} \tag{ 12 }$$

where ${M}_{\mathrm{ej}}$ is the total ejecta mass, ${E}_{{\rm{k}}}$ is the ejecta kinetic energy, and ${t}_{\mathrm{ex}}$ is the time since explosion. The radial density profile is thus completely set by the choice of ${M}_{\mathrm{ej}},\,{E}_{{\rm{k}}}$, and the exponents δ, n. In our calculations, we cut off the model at a radius that encompasses 99% of the total ejecta mass.

We find that the values $\delta =0,\,n=10$ give reasonable fits to the nebular spectra of SNe Ia and so use these values for our fiducial model. In this case, the characteristic velocity and density scales are

$$\begin{eqnarray}&&{v}_{{\rm{t}}}={\rm{10,943}}\,{E}_{51}^{1/2}\,{M}_{1}^{-1/2}\,\mathrm{km}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\rho }_{0}=4.90\times {10}^{-17}\,{E}_{51}^{-3/2}\,{M}_{1}^{5/2}\,{t}_{200}^{-3}\,{\rm{g}}\,{\mathrm{cm}}^{-3},\end{eqnarray} \tag{ 14 }$$

where ${M}_{1}={M}_{\mathrm{ej}}/{M}_{\odot }$, E₅₁ is the kinetic energy in units of 10⁵¹ erg, and ${t}_{200}$ is the time since explosion scaled to 200 days. We explore using different power-law exponents, as well as an exponential density profile, in Section 4.7.

The compositional structure of the ejecta models is assumed to be stratified into three distinct zones (Woosley et al. 2007). The center of the ejecta is assumed to consist of stable IGEs of mass ${M}_{\mathrm{IGE}}$. The stable IGEs are assumed to be composed of a ratio, ${R}_{\mathrm{stb}}$, of ⁵⁴Fe to ⁵⁸Ni. Surrounding the stable IGE region is a zone consisting (initially) mostly of ⁵⁶Ni of mass ${M}_{56\mathrm{Ni}}$. We include a small amount of stable IGEs in this region with mass abundance ${X}_{\mathrm{stb}}$ and the same isotopic ratio ${R}_{\mathrm{stb}}$. Above the ⁵⁶Ni is an outer layer of intermediate-mass elements (IMEs) of mass ${M}_{\mathrm{IME}}$. We compose the IME layer of 70% ²⁸Si, 29% ³²S, and 1% ⁴⁰Ca, roughly consistent with the nucleosynthetic results in the SN Ia explosion models of Plewa (2007) and Seitenzahl et al. (2013). We study the presence of unburned C/O mixed into the nickel zone and IME layer in Section 4.6. The parameters describing the masses of the elements are constrained to add to the total ejecta mass.

The total radioactive energy deposition rate (and hence bolometric luminosity) of a model depends not only on ${M}_{56\mathrm{Ni}}$ but also on the efficiency of the trapping of radioactive decay products. Since the ejecta are largely transparent to gamma rays at nebular phases, the gamma-ray trapping fraction is ${f}_{\gamma ,c}=1-{e}^{-{\tau }_{\gamma }}\approx {\tau }_{\gamma }$, where ${\tau }_{\gamma }$ is the mean optical depth to gamma rays. Taking a typical gamma-ray opacity ${\kappa }_{\gamma }\sim 0.03\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ (Swartz et al. 1995) and integrating the radial optical depth from the center (r = 0) gives an estimate

$$\begin{eqnarray}&&{f}_{\gamma ,c}\approx 0.025\,{E}_{51}^{-1}\,{M}_{1}^{2}\,{t}_{200}^{-2}.\end{eqnarray} \tag{ 15 }$$

Equation (15) presumably overestimates the gamma-ray trapping fraction (since ${\tau }_{\gamma }$ is evaluated at r = 0), but the scaling of ${f}_{\gamma ,c}$ with ${M}_{\mathrm{ej}}$, ${E}_{{\rm{K}}}$, and t will be useful for interpreting how the radioactive deposition rate depends on physical parameters.

To describe the basic features of nebular spectrum formation, we present first a fiducial model with parameters (given in Table 1) typical of standard SN Ia explosion models, i.e., ${M}_{\mathrm{ej}}$ near the Chandrasekhar mass and ${M}_{56\mathrm{Ni}}=0.6\,{M}_{\odot }$. The fiducial model has a transition velocity of ${\rm{10,131}}\,\mathrm{km}\,{{\rm{s}}}^{-1}$, while the interface between the nickel zone and IME layer is near $8800\,\mathrm{km}\,{{\rm{s}}}^{-1}$.

Table 1.  Fiducial Model Parameters

${t}_{\mathrm{ex}}$ (days)	200
${M}_{\mathrm{ej}}$ (${M}_{\odot }$)	1.40
${E}_{{\rm{K}}}$ (1051 erg)	1.2
${M}_{56\mathrm{Ni}}$ (${M}_{\odot }$)	0.6
${M}_{\mathrm{stb}}$ (${M}_{\odot }$)	0.0
${X}_{\mathrm{stb}}$	0.05
${R}_{\mathrm{stb}}$	1.0
${M}_{\mathrm{IME}}$ (${M}_{\odot }$)	0.75
${M}_{\mathrm{CO}}$ (${M}_{\odot }$)	0.0

Note. ${M}_{\mathrm{stb}}$ refers to the mass of the stable IGE region in the core of the explosion; ${X}_{\mathrm{stb}}$ is the mass fraction of stable IGE material mixed into the ⁵⁶Ni region; ${R}_{\mathrm{stb}}$ is the ratio of ⁵⁴Fe/⁵⁸Ni in both core stable IGE and ⁵⁶Ni regions; ${M}_{\mathrm{IME}}$ is the total mass of the IME layer; ${M}_{\mathrm{CO}}$ is the mass of C/O mixed throughout the ejecta.

Download table as:  ASCIITypeset image

Figure 1 shows our calculation of the synthetic nebular spectrum of the fiducial model at 200 days after explosion. We compare to the observed spectrum of the well-studied nearby SN Ia SN 2011fe (Nugent et al. 2011; Shappee et al. 2013; Mazzali et al. 2015). The model spectrum reproduces most of the prominent features.

Zoom In Zoom Out Reset image size

Figure 2 shows a breakdown of the contribution from various ions to the fiducial model spectrum. The strongest features are due to forbidden transitions of Fe ii and Fe iii that are collisionally excited in the nickel zone. Emission due to IME lines is also visible at redder wavelengths. Most of the important individual line transitions are listed in Table 2.

Zoom In Zoom Out Reset image size

Table 2.  Line Identifications for Optical SN Ia Nebular Spectra

Ion	${\lambda }_{{ij}}$ (Å)	${A}_{{ij}}$ (s−1)	Conf. (${}^{2{\rm{S}}+1}{L}_{{\rm{J}}}^{\pi }$)	${E}_{i}$ (eV)	${E}_{j}$ (eV)
S ii	4069	0.2250	${{\rm{S}}}_{3/2}^{{\rm{o}}}{\mbox{--}}^{2}{{\rm{P}}}_{3/2}^{{\rm{o}}}$	0.00	3.05
Fe iii	4658	0.4500	${}^{5}{{\rm{D}}}_{4}^{{\rm{e}}}{\mbox{--}}_{2}^{3}{\rm{F}}{2}_{4}^{{\rm{e}}}$	0.00	2.66
O iii	5007	0.0196	${}^{3}{{\rm{P}}}_{2}^{{\rm{e}}}{\mbox{--}}^{1}{{\rm{D}}}_{2}^{{\rm{e}}}$	0.04	2.51
Fe iii	5011	0.5400	${}^{5}{{\rm{D}}}_{2}^{{\rm{e}}}{\mbox{--}}^{3}{\rm{P}}{2}_{1}^{{\rm{e}}}$	0.09	2.57
Fe ii	5159	0.5805	${{\rm{a}}}^{4}{{\rm{F}}}_{9/2}^{{\rm{e}}}\mbox{--}{{\rm{a}}}^{4}{{\rm{H}}}_{13/2}^{{\rm{e}}}$	0.23	2.63
Fe iii	5270	0.4200	${}^{5}{{\rm{D}}}_{3}^{{\rm{e}}}{\mbox{--}}^{3}{\rm{P}}{2}_{2}^{{\rm{e}}}$	0.05	2.41
Co iii	5888	0.4001	${{\rm{a}}}^{4}{{\rm{F}}}_{9/2}^{{\rm{e}}}\mbox{--}{{\rm{a}}}^{2}{{\rm{G}}}_{9/2}^{{\rm{e}}}$	0.00	2.10
Co iii	6128	0.1100	${{\rm{a}}}^{4}{{\rm{F}}}_{5/2}^{{\rm{e}}}\mbox{--}{{\rm{a}}}^{2}{{\rm{G}}}_{7/2}^{{\rm{e}}}$	0.18	2.20
Co iii	6195	0.1200	${{\rm{a}}}^{4}{{\rm{F}}}_{7/2}^{{\rm{e}}}\mbox{--}{{\rm{a}}}^{2}{{\rm{G}}}_{9/2}^{{\rm{e}}}$	0.10	2.10
Fe ii	7155	0.1495	${{\rm{a}}}^{4}{{\rm{F}}}_{9/2}^{{\rm{e}}}\mbox{--}{{\rm{a}}}^{2}{{\rm{G}}}_{9/2}^{{\rm{e}}}$	0.23	1.96
Ca ii	7291	0.803	${}^{2}{{\rm{S}}}_{1/2}^{{\rm{e}}}{\mbox{--}}^{2}{{\rm{D}}}_{5/2}^{{\rm{e}}}$	0.00	1.70
Ni ii	7378	0.1955	${}^{2}{{\rm{D}}}_{5/2}^{{\rm{e}}}{\mbox{--}}^{2}{{\rm{F}}}_{7/2}^{{\rm{e}}}$	0.00	1.16
Fe ii	7388	0.0435	${}^{4}{{\rm{F}}}_{5/2}^{{\rm{e}}}\mbox{--}{{\rm{a}}}^{2}{{\rm{G}}}_{7/2}^{{\rm{e}}}$	0.35	2.03
Fe ii	7453	0.0485	${{\rm{a}}}^{4}{{\rm{F}}}_{7/2}^{{\rm{e}}}\mbox{--}{{\rm{a}}}^{2}{{\rm{G}}}_{9/2}^{{\rm{e}}}$	0.30	1.96
Fe ii	7638	0.0070	${{\rm{a}}}^{6}{{\rm{D}}}_{7/2}^{{\rm{e}}}\mbox{--}{{\rm{a}}}^{4}{{\rm{P}}}_{5/2}^{{\rm{e}}}$	0.05	1.67
Fe ii	8617	0.0334	${{\rm{a}}}^{4}{{\rm{F}}}_{9/2}^{{\rm{e}}}\mbox{--}{{\rm{a}}}^{4}{{\rm{P}}}_{5/2}^{{\rm{e}}}$	0.23	1.67
Fe iii	8729	0.0495	${}^{3}{\rm{P}}{2}_{2}^{{\rm{e}}}{\mbox{--}}^{3}{{\rm{D}}}_{3}^{{\rm{e}}}$	2.41	3.83
S iii	9069	0.0221	${}^{3}{{\rm{P}}}_{1}^{{\rm{e}}}{\mbox{--}}^{1}{{\rm{D}}}_{2}^{{\rm{e}}}$	0.04	1.40
S iii	9531	0.0576	${}^{3}{{\rm{P}}}_{2}^{{\rm{e}}}{\mbox{--}}^{1}{{\rm{D}}}_{2}^{{\rm{e}}}$	0.10	1.40
S ii	10336	0.1630	${}^{2}{{\rm{D}}}_{3/2}^{{\rm{o}}}{\mbox{--}}^{2}{{\rm{P}}}_{1/2}^{{\rm{o}}}$	1.84	3.04
S ii	10370	0.0779	${}^{2}{{\rm{D}}}_{5/2}^{{\rm{o}}}{\mbox{--}}^{2}{{\rm{P}}}_{1/2}^{{\rm{o}}}$	1.85	3.04

Download table as:  ASCIITypeset image

The spectrum in the nebular phase forms primarily from the collisional excitation of ions by thermal electrons, followed by spontaneous de-excitation via the emission of a line photon. Given the relatively low ejecta temperatures, electrons can only excite low-lying atomic levels, among which radiative transitions are typically electron dipole forbidden (see Table 2). Nevertheless, the rate of collisional de-excitations is generally so small in the low-density nebula that essentially every collisional excitation eventually leads to radiative emission through a forbidden line.

The temperature of the ejecta is determined by the balance of radioactive heating and cooling by line emission. Figure 3 shows the radioactive heating rate and temperature for the fiducial model. The interior ejecta density profile is flat in this model, and so the ejecta temperature is nearly constant at $T\approx 9000$ K in the inner layers. Above the radioactive nickel zone, the temperature drops, due to the declining heating rate, but increases again in the very low density outermost layers, due to the inability of the ejecta to cool efficiently.

Zoom In Zoom Out Reset image size

The important emission lines appearing in the nebular spectra depend on the composition and ionization state of the ejecta. In SNe Ia, ionization is primarily caused by the nonthermal electrons produced by radioactive decay (the collisional ionization rate from thermal electrons is subdominant), which is balanced by radiative recombination. Figure 3 shows the radial dependence of the ionization fractions of iron for the fiducial model. Though Fe iv and Fe v are the most abundant ions, they lack low-lying levels that are able to be excited by the thermal electrons. Thus, the most prominent lines are due to Fe iii and Fe ii.

As a comparison, we also run a model with the same parameters as the fiducial model but with an exponential density profile similar to the commonly used W7 model (Nomoto et al. 1984; Thielemann et al. 1986). Specifically, we use $\rho ={\rho }_{0}{e}^{-v/{v}_{e}}$, where ${\rho }_{0}$ and ${v}_{e}$ are set by total ejecta mass, kinetic energy, and time since explosion. We show the calculated ejecta properties in Figure 4. Higher central density in this model results in higher gamma-ray deposition overall. In general, the degree of ionization increases as the density declines at higher velocities, reflecting the reduced rate of radiative recombination. For the same reason, the overall ionization state of iron is lower in the W7-like model, due to higher overall densities. See Section 4.7 for synthetic spectra of this model and a systematic study of density profiles.

Zoom In Zoom Out Reset image size

As a more comprehensive description of SN Ia nebular spectra, we discuss the features appearing in each key wavelength region seen in Figure 1:

3500–4500 Å: This region contains emission from [S ii] and [Fe ii] transitions. Similar to other studies (e.g., Mazzali et al. 2015; Friesen et al. 2017), our model fails to reproduce all of the observed features, which could be a result of ions missing in the model, uncertainties or incompleteness of the atomic line data, or the neglect of optical depth effects that may produce a pseudo-continuum at the bluest wavelengths. There is some evidence that [Fe ii] emission at 4400 Å contributes to this feature in SN 2011fe (Friesen et al. 2017).

4500–5500 Å: This region is dominated by emission from [Fe iii], which produces prominent features at 4658 and 5270 Å. The latter feature also includes a significant contribution from [Fe ii] λ5159; therefore, the ratio of these two lines is a diagnostic of the ionization state of the gas. The small emission line appearing between the two strong lines is due to [Fe iii] λ5011 and can be washed out if the ejecta velocities are too high.

5500–7000 Å: This region contains a prominent [Co iii] feature near 5888 Å with two smaller [Co iii] features immediately to the red. The line strength depends on the abundance and ionization state of cobalt and typically declines with time as ⁵⁶Co decays to ⁵⁶Fe. The feature near 6500 Å is not well fit by our model (or by previous models; e.g., Mazzali et al. 2015). Because stripped hydrogen from a nondegenerate companion should have velocities $\lesssim 1000\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (e.g., Marietta et al. 2000), this emission is unlikely to be due to Hα. The identification of the feature and the reason for the poor fit are thus unclear.

7000–7700 Å: The feature in this region is a blend of multiple [Fe ii] lines and a broad [Ca ii] line. In addition, [Ni ii] emission can contribute if sufficient stable ⁵⁸Ni is present in the gas, but for our fiducial model [Ni ii] does not dominate this feature. Notably, the fiducial model only contains 0.0075 ${M}_{\odot }$ of calcium, and yet [Ca ii] emission from the IME layer dominates the emission in this region.

8500–11000 Å: This region is made up primarily of [S ii] and [S iii] emitted by the 0.22 ${M}_{\odot }$ of sulfur in the IME layer above the nickel zone. The broad, flat-topped profiles of IMEs are due to their large velocities and absence of IMEs in the core.

Figure 5 shows the nebular spectrum of the fiducial model at times between 200 and 400 days after explosion. As the SN evolves in time, ejecta densities decline and radioactive isotopes decay. The bolometric luminosity drops with time owing to the declining radioactive heating and a decreasing gamma-ray trapping fraction. The relative strength of prominent Fe lines does not evolve significantly, indicating that the Fe iii/Fe ii ionization ratio remains fairly constant over time. However, as gamma-ray trapping becomes inefficient at late times, the relative strength of IME features to IGE features decreases.

Zoom In Zoom Out Reset image size

The biggest change in the spectral features over time is the strength of the [Co iii] emission at 5888 Å, which decreases as ⁵⁶Co decays to ⁵⁶Fe. This behavior has been seen in observations (Childress et al. 2015) and exploited to derive explosion properties such as ⁵⁶Ni mass and total ejecta mass. In addition, the declining gamma-ray deposition in the IME layer results in lowered relative strength of IME features.

The amount of kinetic energy (${E}_{{\rm{K}}}$) imparted to the ejecta by an SN explosion depends on the nucleosynthetic yields and initial binding energy (Woosley et al. 2007). The typical energy of SNe Ia is around a Bethe ($1{\rm{B}}={10}^{51}\,\mathrm{erg}$). Various theoretical models have predicted kinetic energies in the range of 0.87–1.6 B (Gamezo et al. 2005; Golombek & Niemeyer 2005; Plewa 2007; Röpke et al. 2007; Röpke & Niemeyer 2007; Jordan et al. 2008, 2012; Bravo et al. 2009).

Figure 6 shows synthetic spectra of the fiducial model with ${E}_{{\rm{K}}}$ varied between 1 and 2 B. The lower ${E}_{{\rm{K}}}$ models are more efficient at trapping radioactive energy (due to the higher ejecta density, Equation (14)) and so have higher bolometric luminosities.

Zoom In Zoom Out Reset image size

Changing ${E}_{{\rm{K}}}$ has a modest effect on the shape of the spectral features. An increase of ${E}_{{\rm{K}}}$ from 1 to 1.4 B increases the velocity scale by only 18%, which results in a subtle increase in the line widths (as these widths are also set in part by line blending). A larger change of ${E}_{{\rm{K}}}$ by a factor of 2 (from 1 to 2 B) does have noticeable effects, causing the Fe iii/Fe ii complex to be so blended that the small central emission near 5000 Å becomes indistinguishable. Observations of this small feature may therefore be a useful diagnostic of the velocity of the nickel zone in SNe Ia.

Theoretical models of SNe Ia predict an ejected mass in the range of 0.8–2.0 ${M}_{\odot }$, depending on the progenitor scenario. Approximate light-curve modeling studies have suggested that observed SNe Ia could span this entire range (Stritzinger et al. 2006; Scalzo et al. 2010, 2014a, 2014b).

Figure 7 shows synthetic spectra of the fiducial model with ${M}_{\mathrm{ej}}$ varied between 1 and 2 ${M}_{\odot }$. For a fixed kinetic energy, a higher ${M}_{\mathrm{ej}}$ results in higher ejecta densities and lower velocities. As a result, higher-${M}_{\mathrm{ej}}$ models have greater gamma-ray trapping, a brighter bolometric luminosity, and less Doppler-broadened spectral features. This effect of ${M}_{\mathrm{ej}}$ is therefore somewhat degenerate with that of kinetic energy.

Zoom In Zoom Out Reset image size

The relative features in the synthetic spectra of Figure 7 show only subtle variations with ${M}_{\mathrm{ej}}$. For our super-${M}_{\mathrm{Ch}}$ case with ${M}_{\mathrm{ej}}=2.0\,{M}_{\odot }$, the IME features become visibly stronger as a result of the higher total IME mass. The Fe iii/Fe ii line ratio also decreases due to the increased rate of recombination at higher densities.

While we have studied the effect of varying ejecta mass alone in Figure 7, this parameter is likely correlated with kinetic energy and ⁵⁶Ni mass (Woosley et al. 2007); for example, a super-${M}_{\mathrm{Ch}}$ explosion is likely to produce more ⁵⁶Ni and higher kinetic energy. We therefore ran three additional models in which we fixed the ratios ${M}_{56\mathrm{Ni}}/{M}_{\mathrm{ej}}=1/2$ and ${E}_{51}/{M}_{\mathrm{ej}}=1.2/1.4$. Based on the above discussion of gamma-ray trapping (Section 3.1), we expect luminosity to approximately scale as ${L}_{\mathrm{bol}}\sim {M}_{\mathrm{ej}}^{2}$ if these ratios are held constant.

Figure 8 shows synthetic spectra of these scaled models with total ejecta masses of 1.0, 1.4, and $2.0\,{M}_{\odot }$. We find that the spectral features are remarkably unchanged despite a substantial variation in the masses. The only major difference is the bolometric luminosity. This indicates that, in a generic sense, the nebular spectra of SNe Ia are consistent with non-${M}_{\mathrm{Ch}}$ models, provided that ${E}_{{\rm{k}}}$ and ${M}_{56\mathrm{Ni}}$ scale accordingly.

Zoom In Zoom Out Reset image size

A number of observational studies and theoretical models indicate that the ⁵⁶Ni masses of normal SNe Ia range from 0.3 to 1.2 ${M}_{\odot }$(Gamezo et al. 2005; Mazzali et al. 2007b; Plewa 2007; Röpke & Niemeyer 2007; Röpke et al. 2007; Jordan et al. 2008, 2012; Bravo et al. 2009; Raskin et al. 2009; Rosswog et al. 2009; Seitenzahl et al. 2011, 2013), with a typical value near 0.6 ${M}_{\odot }$ (Branch & Khokhlov 1995).

Figure 9 shows synthetic spectra of the fiducial model with ${M}_{56\mathrm{Ni}}$ varied between 0.4 and 0.8 ${M}_{\odot }$. Naturally, the bolometric luminosity increases proportionally with ${M}_{56\mathrm{Ni}}$. The line ratios in these spectra are relatively insensitive to ${M}_{56\mathrm{Ni}}$, with the exception being greater blending around the [Fe iii] λ5011 feature in higher ⁵⁶Ni mass models due to the larger size of the ⁵⁶Ni core. The decrease in IME emission in higher-${M}_{56\mathrm{Ni}}$ models is due to the lower total mass of IMEs in these models by construction (since ${M}_{56\mathrm{Ni}}$ + ${M}_{\mathrm{IME}}$ is held fixed).

Zoom In Zoom Out Reset image size

Stable IGEs are produced in SN Ia ejecta in two distinct ways. The trace presence of neutron-rich isotopes (in particular ²²Ne) due to the metallicity of the progenitor WD and pre-explosion carbon simmering leads to the production of up to ≈25% by mass of neutronized IGEs (in particular ⁵⁴Fe and ⁵⁸Ni) throughout the nickel core (Timmes et al. 2003; Piro & Bildsten 2008; Martínez-Rodríguez et al. 2016). In addition, electron capture occurring in nuclear burning at high central densities can lead to the production of ∼0.05–0.4 ${M}_{\odot }$ of neutronized IGEs (Nomoto et al. 1984; Thielemann et al. 1986; Seitenzahl et al. 2011, 2013). The latter effect only occurs in white dwarfs with ${M}_{\mathrm{ej}}\gtrsim {M}_{\mathrm{ch}}$ (due to their higher central densities) and is also influenced by the timing of a possible deflagration-to-detonation transition (Seitenzahl et al. 2013). 1D ${M}_{\mathrm{Ch}}$ models often predict stable IGEs to be produced at the core (Nomoto et al. 1984; Mazzali et al. 2007b), while multidimensional simulations indicate that buoyancy should mix IGEs throughout the ⁵⁶Ni region (Gamezo et al. 2005; Kasen et al. 2009; Seitenzahl et al. 2011). As a result, nebular spectral indicators of stable IGEs would be valuable for inferring the progenitor scenario.

Figure 10 shows synthetic spectra of the fiducial model varying ${M}_{\mathrm{stb}}$, the mass of a neutron-rich core, between 0.05 and 0.20 ${M}_{\odot }$. The stable IGEs are assumed to be an equal mix of ⁵⁴Fe and ⁵⁸Ni. As expected, the [Ni ii] λ7378 feature becomes apparent with a stable core mass of 0.05 ${M}_{\odot }$, which corresponds to 0.025 ${M}_{\odot }$ of ⁵⁸Ni. [Fe ii] λ7388 emission can also produce a peak near these wavelengths, in some cases dominating the feature, suggesting that [Ni ii] may not even be needed to fit this peak (see Section 4.7).

Zoom In Zoom Out Reset image size

Higher stable core mass also increases the relative fraction of Fe ii compared to Fe iii, which is a result of lowered nonthermal ionization in the nonradioactive core region. Another noteworthy product of a stable core is the flat-topped profile of the [Co iii] λ5888 feature, produced by the absence of cobalt in the core region. The lower flux in IME features is due to the lower mass of IMEs in higher-${M}_{\mathrm{stb}}$ models by construction.

A comparison between two ${M}_{\mathrm{stb}}$ models and the fiducial model with a W7-like exponential density profile is shown in Figure 11. In the ${M}_{\mathrm{stb}}=0.2\,{M}_{\odot }$ model we find an apparent blueshift of the [Ni ii] λ7378 peak by ≈1500 km s⁻¹ due to the relative blending of this line with [Fe ii] λ7155. The feature is also sensitive to [Ca ii] emission, which is stronger in the low-${M}_{\mathrm{stb}}$ model. The [Fe ii] λ7453 feature redward of [Ni ii] λ7378 may also dominate in some models, resulting in an apparent redshift of the composite peak (see also Section 4.7). In the exponential model, both peaks are redshifted owing to line blending, which may lead to misinterpretation as ejecta asymmetry.

Zoom In Zoom Out Reset image size

Maeda et al. (2010a) find shifts in the λ7378 peak of up to 3000 km s⁻¹, which they interpret as indications of ejecta asymmetry. Our models suggest that this geometrical interpretation is complicated by line blending. Isolating the [Ni ii] emission in this feature can been attempted (Maeda et al. 2010b) but is subject to model-dependent uncertainties.

Surprisingly, the mixing of stable nickel through the nickel zone (i.e., varying the parameter ${X}_{\mathrm{stb}}$ between 0.05 and 0.20) does not produce a visible effect on the synthetic spectrum, even for ⁵⁸Ni masses comparable to those in the core stable IGE region of Figure 10. This confirms the findings of Maeda et al. (2010b). This is due to the high level of ionization of Ni in the nickel zone, which suppresses [Ni ii] emission. For the same reason, the ratio of ⁵⁴Fe to ⁵⁸Ni (${R}_{\mathrm{stb}}$) in the nickel zone also had no visible effect on the synthetic spectra in the range (0.5–1.5) that we tested. We predict that up to 0.1 ${M}_{\odot }$ of stable nickel can be hidden in the ejecta of an SN Ia model with fiducial parameters as a result of this ionization effect.

Observational studies have estimated that about 30% of SNe Ia show carbon in their early-time spectra (Maoz et al. 2014, and references therein). The nearby SN 2011fe also showed both carbon and oxygen in very early observations (Nugent et al. 2011; Mazzali et al. 2014), and hydrodynamical simulations have predicted various amounts of unburned C/O material mixed throughout the ejecta (Röpke 2005; Pakmor et al. 2012; Seitenzahl et al. 2013; Moll et al. 2014). 3D nebular modeling by Kozma et al. (2005) showed clear [O i] features for a pure deflagration model of Röpke (2005), which contained 0.6 ${M}_{\odot }$ of unburned C/O material. There have also been detections of possible [O i] emission in subluminous SNe Ia (Kromer et al. 2013; Taubenberger et al. 2013).

Figure 12 shows synthetic spectra of the fiducial model with C/O mass varied between 0.1 and 0.4 ${M}_{\odot }$. We keep a fixed carbon–oxygen ratio of 1:9 and mix C/O into both the nickel zone and IME layer with the same mass fraction, consistent with the expected nucleosynthetic yields of delayed detonation explosions (Seitenzahl et al. 2013).

Zoom In Zoom Out Reset image size

We find a strong contribution of [O iii] at 5007 Å, as well as a weak contribution of [O i] λ6300/λ6364 and a blend of [O ii] features at 7320 Å. The high ionization state of oxygen prevents the [O i] emission seen in Kozma et al. (2005) from contributing significantly, but we expect that [O i] would become stronger if oxygen were more concentrated in the higher-density central regions, or if the mass of ⁵⁶Ni (and hence radioactive deposition) were lower. In addition, our 1D models assume microscopic mixing of the C/O with IGEs; multidimensional models may find that clumps of oxygen isolated from ⁵⁶Co have lower ionization. C/O within the IME layer does not produce a significant nebular feature.

The choice of the ejecta density profile is a critical input into nebular modeling. A broken power law with a relatively flat interior was found to approximate the density structure in some 2D delayed detonation models (Kasen 2010), while a steeper exponential profile more closely fits the structure of the commonly used 1D W7 model (Nomoto et al. 1984; Thielemann et al. 1986). For the case of SN 2003hv, Mazzali et al. (2011) found that a core of lowered density provided a better fit to observed nebular spectra than a W7-like density profile. Here, we attempt to illuminate the effect of steepening our interior power-law density profiles. We also consider an exponential profile of the form $\rho ={\rho }_{0}{e}^{-v/{v}_{e}}$, where ${\rho }_{0}$ and v_e can be determined from ${M}_{\mathrm{ej}}$, ${E}_{{\rm{K}}}$, and ${t}_{\mathrm{ex}}$.

Figure 13 shows synthetic spectra of the fiducial model with varied density profiles. Steeper profiles, which concentrate more mass at low velocities, produce stronger and narrower spectral profiles. In particular, an exponential profile helps resolve the individual peaks in the feature around 7300 Å. Figure 11 shows the nebular emission of the exponential model around 7300 Å decomposed into contributions from individual ions. The bluer peak is dominated by [Fe ii] λ7155, while the redder peak is dominated by [Fe ii] λ7388 and [Fe ii] λ7453, with some contribution from [Ni ii] λ7378. [Ca ii] λ7291 adds an offset to the feature. These are consistent with the findings of Ashall et al. (2016) from their W7-like models for SN 1986G. Therefore, [Fe ii] can account for the emission in the redder peak if stable Ni is absent, and it could potentially redshift that peak if stable Ni is present. In either case, we find that isolation of [Ni ii] λ7378 emission would be difficult.

Zoom In Zoom Out Reset image size

Higher gamma-ray trapping due to higher central densities produces higher bolometric luminosities in steeper density profiles. Furthermore, steeper density profiles lead to lower central temperatures, as cooling becomes more efficient at higher densities. Therefore, the ionization state is lower, which changes the main Fe iii/Fe ii line ratio.

When comparing to the nebular spectrum of SN 2011fe, we find that a broken power law with a flat interior profile ($\delta =0$) best reproduces the shape of the features and the main iron line ratios, whereas an exponential density produces lines that are too centrally peaked and overestimates the Fe ii/Fe iii ratio (indicated by high flux in the [Fe ii] feature at 5159 Å). While there is some degeneracy with other parameters, such as kinetic energy, nebular models may be able constrain the interior ejecta density, which should be useful in testing explosion models.

One important factor for the modeling of nebular spectra is the extensive atomic data inputs. Uncertainties in published atomic data, discrepancies between different sources for the same data, and the need for crude approximations where data are lacking impact the model predictions. Here we attempt to quantify some of the uncertainties by systematically changing the values of three of the most important and uncertain atomic data—the collisional ionization cross sections (${Q}_{{\rm{k}}}$), the thermal collisional excitation rates (${C}_{{ij}}$), and the radiative recombination rates (α).

We show the impact of systematically varying atomic data in Figures 14–15. To explore the effect of the collisional ionization cross sections (Figure 14), we carried out calculations in which (a) all ionization cross sections were increased by a factor of 2, (b) all ionization cross sections were decreased by a factor of 2, (c) the cross section of Fe ii was decreased by a factor of 2 and that of Fe iii was increased by a factor of 2, and (d) the cross section of Fe ii was increased by a factor of 2 and that of Fe iii was decreased by a factor of 2. These variations in Q affect the ionization ratio of Fe iii/Fe ii and can modify the Fe iii/Fe ii line ratio by up to ±25%. There are also indirect effects, since changes to the atomic data of one element can alter the calculated gas temperature and so result in changes in the emission from other species.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To explore the effect of the thermal collisional excitation rates (Figure 15), we carried out calculations in which all ${\sigma }_{{\text{}}\mathrm{ij}}$ were increased or decreased by a factor of 2. These variations in ${\sigma }_{{\text{}}\mathrm{ij}}$ affect the strength of emission features and can modify the Fe iii/Fe ii line ratio by up to ±10%.

We show in Figure 16 the fiducial model calculated using recombination rates from two different databases. Our calculations throughout have used recombination rates from the CHIANTI database, but more recent data for Fe i–v are available from the Nahar OSU Radiative Atomic Database (NORAD; Nahar 1996, 1997; Nahar et al. 1997, 1998). With the latter data set, we have access to state-specific recombination rates, which we neglected in the above treatment. We also neglected CT in the above analyses, which becomes important for a nickel zone with low ionization state (as is the case with NORAD recombination rates). We therefore include CT in Fe i– iv with data from Krstic et al. (1997).

Zoom In Zoom Out Reset image size

While there is reasonable agreement between the synthetic spectra, the strength of emission features does depend on the atomic data set. In particular, the Fe iii/Fe ii line ratio shows a higher Fe ii population when NORAD recombination rates are used, due to higher total recombination rates for Fe i–iii in the NORAD data set, leading to a ∼30% decrease in Fe iii/Fe ii line ratio. Similar variations are seen in most other spectral features.

To determine collisional excitation and de-excitation rates relevant for the thermal pool of electrons, we assume that electrons adhere to a Maxwell–Boltzmann distribution. The subsequent "effective" or "Maxwellian-averaged" collisional strength can be calculated for each transition using

$$\begin{eqnarray}&&{{\rm{\Upsilon }}}_{{ij}}={\int }_{0}^{\infty }{{\rm{\Omega }}}_{{ij}}{e}^{-{\epsilon }_{j}/{kT}}d({\epsilon }_{j}/{kT}),\end{eqnarray} \tag{ 28 }$$

where Ω_ij is the collisional cross section for an electron to excite the i–j transition and T is the electron temperature. These collision strengths can be found in the literature and data tables.

The subsequent rates of collisional excitation and de-excitation, respectively, are

$$\begin{eqnarray}&&{C}_{{ij}}(T)=\displaystyle \frac{8.63\times {10}^{-6}}{{g}_{i}\sqrt{T}}{e}^{-{E}_{{ij}}/{kT}}{{\rm{\Upsilon }}}_{{ij}}{n}_{e}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{C}_{{ji}}(T)=\displaystyle \frac{{g}_{i}}{{g}_{j}}{e}^{{E}_{{ij}}/{k}_{{\rm{B}}}T}{C}_{{ij}},\end{eqnarray} \tag{ 30 }$$

where E_ij is the energy difference between the two levels and $j\gt i$.

In cases where atomic data for ϒ are unavailable, one can use the dipole approximation (van Regemorter 1962; Jefferies 1968), which gives results for neutral atoms and positive ions, respectively,

$$\begin{eqnarray}&&{C}_{{ij}}(T)\\ &&\quad =\,\left\{\begin{array}{l}2.16{f}_{{ij}}{n}_{e}{\left[\displaystyle \frac{{E}_{{ij}}}{{kT}}\right]}^{-1.68}{T}^{-3/2}{e}^{-{E}_{{ij}}/{kT}}\ {{\rm{s}}}^{-1},\,(\mathrm{neutral}\,\mathrm{atoms})\\ 3.9{f}_{{ij}}{n}_{e}{\left[\displaystyle \frac{{E}_{{ij}}}{{kT}}\right]}^{-1}{T}^{-3/2}{e}^{-{E}_{{ij}}/{kT}}\,{{\rm{s}}}^{-1},\,(\mathrm{positive}\,\mathrm{ions})\end{array}\right.,\end{eqnarray} \tag{ 31 }$$

where ${f}_{{ij}}$ is oscillator strength, T is temperature in kelvin, ${E}_{{ij}}$ is transition energy in eV, and ${n}_{e}$ is electron density in cm⁻³.

In practice, atomic data for ${{\rm{\Upsilon }}}_{{ij}}$ are often published in a narrow temperature range; for temperatures outside this range, we use the above approximation to fill in the missing data. We scale Equation (31) so that there is no discontinuity in ${C}_{{ij}}$.

Forbidden transitions for which atomic data are not available should be treated as a special case, since oscillator strengths for these transitions produce artificially low collisional rates. Axelrod (1980) provides the following formulation for forbidden transitions (defined as ${f}_{{ij}}\leqslant 0.001$):

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{ij},\mathrm{forb}}=\left\{\begin{array}{l}0.00375\,{g}_{i}{g}_{j},\,\lambda \leqslant 10\mu \\ 0.0225\,{g}_{i}{g}_{j},\,\lambda \gt 10\mu \end{array}\right..\end{eqnarray} \tag{ 32 }$$

Total recombination rates of relevant ions include a sum over all subshells and incorporate both radiative and dielectronic recombination. We assume that recombinations involve only ground states, a good approximation in regimes where collisional excitation dominates the population of excited states. We use the rates published in the CHIANTI database v.8 (Dere et al. 1997; Del Zanna et al. 2015). Although level-specific recombination rates are available for certain ions (Nahar & Pradhan 1992, 1994, etc.), we opt to use a consistent source for recombination rates for this work. We explore the latter data set in Section 4.8 and plan to implement state-specific recombination rates in any future work where recombination lines might contribute significantly.