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

TL;DR

This study develops a non-LTE radiative transfer model to analyze how nebular spectra of Type Ia supernovae depend on explosion parameters, revealing spectral insensitivity to many parameters and highlighting diagnostic features.

Contribution

It introduces a systematic non-LTE modeling approach for Type Ia supernova nebular spectra, assessing parameter dependencies and uncertainties, which was not comprehensively done before.

Findings

Spectral features are largely insensitive to most physical parameters.

Degeneracy among parameters prevents unique ejecta property determination.

[SIII] emission near 9500A is a useful nucleosynthesis diagnostic.

Abstract

We present a radiative transfer code to model the nebular phase spectra of supernovae (SNe) in non-local thermodynamic equilibrium (NLTE). We apply it to a systematic study of Type Ia SNe using parameterized 1D models and show how nebular spectral features depend on key physical parameters, such as the time since explosion, total ejecta mass, kinetic energy, radial density profile, and the masses of 56Ni, intermediate mass elements (IMEs), and stable iron-group elements (IGEs). We also quantify the impact of uncertainties in atomic data inputs. Among the results of the study are: (1) The main features of SNe Ia nebular spectra are relatively insensitive to most physical parameters. Degeneracy among parameters precludes a unique determination of the ejecta properties from spectral fitting. In particular, features can be equally well fit with generic Chandrasekhar mass (M_ch), sub-M_Ch, and super-M_Ch models; (2) A sizable (~0.1 Msun) central region of stable IGEs, often claimed as evidence for M_Ch models, is not essential to fit the optical spectra and may produce an unusual flat-top [CoIII] profile; (3) The strength of [SIII] emission near 9500A can provide a useful diagnostic of explosion nucleosynthesis; (4) Substantial amounts (~0.1 Msun) of unburned C/O mixed throughout the ejecta produce [OIII] emission not seen in observations; (5) Shifts in the wavelength of line peaks, sometimes used to infer ejecta geometry, can also arise from line blending effects; (6) The steepness of the ejecta density profile affects the line shapes, with flatter slopes providing better fits to the observations of SN 2011fe, offering a constraint on explosion models; (7) Uncertainties in atomic data affect spectral line ratios by ~30%, a level similar to the effect of varying physical parameters.

Tables

Table 1. Table 1 Fiducial model parameters. M stb subscript 𝑀 stb {M}\rm_{stb}  refers to the mass of the stable IGE region in the core of the explosion; X stb subscript 𝑋 stb {X}\rm_{stb}  is the mass fraction of stable IGE material mixed into the 56 Ni region; R stb subscript 𝑅 stb {R}\rm_{stb}  is the ratio of 54 Fe/ 58 Ni in both core stable IGE and 56 Ni regions; M IME subscript 𝑀 IME {M}\rm_{IME} is the total mass of the IME layer; M CO subscript 𝑀 CO {M}\rm_{CO}  is the mass of C/O mixed throughout the ejecta.

$t_{\mathrm{ex}}$ (days)	200
$M_{\mathrm{ej}}$ ($M_{\odot }$)	1.40
$E_K$ (${10}^{51}$ erg)	1.2
$M_{56Ni}$ ($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

Table 2. Table 2 Line Identifications for Optical SNe Ia Nebular Spectra

Ion	${\lambda }_{\mathrm{ij}}$ ( Å)	$A_{\mathrm{ij}}$ (s-1)	Conf. (${}_{}{}^{2S+1}L_J^{\pi }$)	$E_i$ (eV)	$E_j$ (eV)
SII	4069	0.2250	$S_{3/2}^o{-}^2{P}_{3/2}^o$	0.00	3.05
FeIII	4658	0.4500	${}_{}{}^{5}D_4^e{-}_2^3{\mathrm{F2}}_4^e$	0.00	2.66
OIII	5007	0.0196	${}_{}{}^{3}P_2^e{-}^1{D}_2^e$	0.04	2.51
FeIII	5011	0.5400	${}_{}{}^{5}D_2^e{-}^3{\mathrm{P2}}_1^e$	0.09	2.57
FeII	5159	0.5805	$a^4{F}_{9/2}^e-a^4{H}_{13/2}^e$	0.23	2.63
FeIII	5270	0.4200	${}_{}{}^{5}D_3^e{-}^3{\mathrm{P2}}_2^e$	0.05	2.41
CoIII	5888	0.4001	$a^4{F}_{9/2}^e-a^2{G}_{9/2}^e$	0.00	2.10
CoIII	6128	0.1100	$a^4{F}_{5/2}^e-a^2{G}_{7/2}^e$	0.18	2.20
CoIII	6195	0.1200	$a^4{F}_{7/2}^e-a^2{G}_{9/2}^e$	0.10	2.10
FeII	7155	0.1495	$a^4{F}_{9/2}^e-a^2{G}_{9/2}^e$	0.23	1.96
CaII	7291	0.803	${}_{}{}^{2}S_{1/2}^e{-}^2{D}_{5/2}^e$	0.00	1.70
NiII	7378	0.1955	${}_{}{}^{2}D_{5/2}^e{-}^2{F}_{7/2}^e$	0.00	1.16
FeII	7388	0.0435	${}_{}{}^{4}F_{5/2}^e-a^2{G}_{7/2}^e$	0.35	2.03
FeII	7453	0.0485	$a^4{F}_{7/2}^e-a^2{G}_{9/2}^e$	0.30	1.96
FeII	7638	0.0070	$a^6{D}_{7/2}^e-a^4{P}_{5/2}^e$	0.05	1.67
FeII	8617	0.0334	$a^4{F}_{9/2}^e-a^4{P}_{5/2}^e$	0.23	1.67
FeIII	8729	0.0495	${}_{}{}^3\mathrm{P2}_2^e{-}^3{D}_3^e$	2.41	3.83
SIII	9069	0.0221	${}_{}{}^{3}P_1^e{-}^1{D}_2^e$	0.04	1.40
SIII	9531	0.0576	${}_{}{}^{3}P_2^e{-}^1{D}_2^e$	0.10	1.40
SII	10336	0.1630	${}_{}{}^{2}D_{3/2}^o{-}^2{P}_{1/2}^o$	1.84	3.04
SII	10370	0.0779	${}_{}{}^{2}D_{5/2}^o{-}^2{P}_{1/2}^o$	1.85	3.04

Table 3. Table 3 List of models included in physical parameter study. The parameter varied from the fiducial value in each model is shown in the model description, and some derived values are shown. v t subscript 𝑣 t {v}\rm_{t}  is the transition velocity of the density profile determined by power-law exponents δ 𝛿 \delta and n. ρ c subscript 𝜌 c \rm\rho_{c} is the central density. The column for M stb , tot subscript 𝑀 stb tot {M}\rm_{stb,tot} , which is the sum of neutron-rich core material and stable isotopes mixed into the nickel zone, reports 54 Fe and 58 Ni masses separately. L bol subscript 𝐿 bol {L}\rm_{bol}  is the total luminosity over all wavelengths. L 4665 / L 5272 subscript L 4665 subscript L 5272 \rm L_{4665}/L_{5272} is the ratio of luminosities at 4665 Å and 5272 Å, the wavelengths at which the prominent [FeIII] and [FeII]/[FeIII] features peak, respectively.

Model	$v_t$	${\rho }_c$	$M_{\mathrm{IME}}$	$M_{\mathrm{stb},\mathrm{tot}}$	$L_{\mathrm{bol}}$	L4665/L5272
	(km/s)	(${10}^{-17}$ g/cm3)	($M_{\odot }$)	Fe, Ni ($M_{\odot }$)	(10${}_{}{}^{40}\mathrm{erg}/s$)
fiducial	10,131	8.65	0.75	0.016, 0.016	7.96	2.30
$t_{\mathrm{ex}}=250\mathrm{days}$	10,131	4.43	0.75	0.016, 0.016	3.93	2.45
$t_{\mathrm{ex}}=200\mathrm{days}$	10,131	2.56	0.75	0.016, 0.016	2.11	2.53
$t_{\mathrm{ex}}=350\mathrm{days}$	10,131	1.61	0.75	0.016, 0.016	1.19	2.54
$t_{\mathrm{ex}}=400\mathrm{days}$	10,131	1.08	0.75	0.016, 0.016	0.70	2.48
$M_{\mathrm{ej}}=1.0M_{\odot }$	12,000	3.73	0.36	0.016, 0.016	5.27	2.75
$M_{\mathrm{ej}}=1.1M_{\odot }$	11,442	4.73	0.46	0.016, 0.016	5.82	2.65
$M_{\mathrm{ej}}=1.2M_{\odot }$	10,954	5.88	0.54	0.016, 0.016	6.58	2.53
$M_{\mathrm{ej}}=1.3M_{\odot }$	10,525	7.19	0.65	0.016, 0.016	7.18	2.41
$M_{\mathrm{ej}}=2.0M_{\odot }$	8,485	21.10	1.33	0.016, 0.016	13.91	1.86
$E_K=1.0B$	9,258	11.37	0.75	0.016, 0.016	8.90	2.15
$E_K=1.1B$	9,710	9.86	0.75	0.016, 0.016	8.33	2.23
$E_K=1.3B$	10,556	7.67	0.75	0.016, 0.016	7.57	2.37
$E_K=1.4B$	10,954	6.86	0.75	0.016, 0.016	7.22	2.44
$E_K=1.5B$	11,339	6.19	0.75	0.016, 0.016	6.96	2.50
$E_K=2.0B$	13,093	4.02	0.75	0.016, 0.016	5.79	2.73
$M_{56Ni}=0.4M_{\odot }$	10,131	8.65	0.96	0.011, 0.011	5.49	2.25
$M_{56Ni}=0.5M_{\odot }$	10,131	8.65	0.85	0.013, 0.013	6.77	2.28
$M_{56Ni}=0.7M_{\odot }$	10,131	8.65	0.64	0.019, 0.019	9.22	2.32
$M_{56Ni}=0.8M_{\odot }$	10,131	8.65	0.55	0.021, 0.021	10.21	2.33
$M_{\mathrm{stb}}=0.05M_{\odot }$	10,131	8.65	0.71	0.040, 0.040	7.75	2.28
$M_{\mathrm{stb}}=0.10M_{\odot }$	10,131	8.65	0.64	0.066, 0.066	7.81	2.25
$M_{\mathrm{stb}}=0.15M_{\odot }$	10,131	8.65	0.59	0.093, 0.093	7.64	2.21
$M_{\mathrm{stb}}=0.20M_{\odot }$	10,131	8.65	0.54	0.116, 0.116	7.55	2.17
$M_{\mathrm{ej}}=1.0\mathrm{scaled}$111In these models, we keep the following ratios fixed: $M_{56Ni}/M_{\mathrm{ej}}=1/2$ and $E_{51}/M_{\mathrm{ej}}=1.2/1.4$. See §4.3 for further details.	10,240	6.15	0.46	0.013,0.013	5.35	2.49
$M_{\mathrm{ej}}=1.4\mathrm{scaled}$	10,240	8.65	0.64	0.019,0.019	9.21	2.32
$M_{\mathrm{ej}}=2.0\mathrm{scaled}$	10,240	12.40	0.91	0.027,0.027	16.88	2.17
$M_{\mathrm{CO}}=0.05M_{\odot }$	10,131	8.65	0.70	0.016, 0.016	7.85	2.25
$M_{\mathrm{CO}}=0.10M_{\odot }$	10,131	8.65	0.64	0.017, 0.017	8.04	2.20
$M_{\mathrm{CO}}=0.15M_{\odot }$	10,131	8.65	0.59	0.018, 0.018	7.90	2.15
$M_{\mathrm{CO}}=0.20M_{\odot }$	10,131	8.65	0.55	0.019, 0.019	7.82	2.09
$M_{\mathrm{CO}}=0.40M_{\odot }$	10,131	8.65	0.33	0.023, 0.023	7.76	1.86
$\delta =0.5$, n = 10222These refer to the power-law exponents in Equation 10. A higher value of $\delta$ produces a steeper density profile.	10,528	66.89	0.75	0.016, 0.016	8.21	2.17
$\delta =1.0$, n = 10	11,098	486.5	0.75	0.016, 0.016	8.95	1.88
exponential333This refers to the exponential density profile explained in §4.7.	N/A	106.6	0.75	0.016, 0.016	12.56	1.63
all $Q_k\times 2$444$Q_k$ is the collisional (electron impact) ionization cross section for the ground state of ion $k$.	10,131	8.65	0.75	0.016, 0.016	7.96	2.49
all $Q_k\times 0.5$	10,131	8.65	0.75	0.016, 0.016	7.96	2.02
$Q_{\mathrm{FeII}}\times 2$, $Q_{\mathrm{FeIII}}\times 0.5$	10,131	8.65	0.75	0.016, 0.016	7.96	2.74
$Q_{\mathrm{FeII}}\times 0.5$, $Q_{\mathrm{FeIII}}\times 2$	10,131	8.65	0.75	0.016, 0.016	7.96	1.72
all ${\sigma }_{\mathrm{ij}}\times 2$555${\sigma }_{\mathrm{ij}}$ is the thermal collisional (electron-impact) excitation rate for a transition $ij$.	10,131	8.65	0.75	0.016, 0.016	7.96	2.10
all ${\sigma }_{\mathrm{ij}}\times 0.5$	10,131	8.65	0.75	0.016, 0.016	7.96	2.52
$\alpha \mathrm{from}\mathrm{NORAD}$	10,131	8.65	0.75	0.016, 0.016	7.96	1.66

Equations

$$\mathbf{M}\vec{n}=0$$

$$n_{\mathrm{I},\mathrm{i}}\left(R_{\mathrm{I},\mathrm{i}}+\sum _{\mathrm{j}\ne \mathrm{i}}{\mathcal{E}}_{\mathrm{ij}}\right)=\sum _{\mathrm{k}\ne \mathrm{i}}{n}_{\mathrm{k}}{\mathcal{E}}_{\mathrm{ki}}+n_{\mathrm{I}+1,\mathrm{gs}}{n}_{\mathrm{e}}{\alpha }_{\mathrm{i}}$$

$$R_{\mathrm{ex}}=\frac{{\dot{ϵ}}_{\mathrm{rad}}{\eta }_{\mathrm{ij}}}{\Delta E_{\mathrm{ij}}{n}_{\mathrm{i}}}$$

$$R_{\mathrm{ion}}=\frac{{\dot{ϵ}}_{\mathrm{rad}}{\eta }_{\mathrm{k}}}{I_{\mathrm{k}}{n}_{\mathrm{k}}}$$

$$n_{\mathrm{e}}=\sum _{\mathrm{k},\mathrm{i}}{n}_{\mathrm{k},\mathrm{i}}i$$

$${\dot{ϵ}}_{\mathrm{ij}}(T)=h{\nu }_{\mathrm{ij}}{A}_{\mathrm{ij}}{n}_{\mathrm{j}}$$

$$\sum _{\mathrm{i},\mathrm{j}}{\dot{ϵ}}_{\mathrm{ij}}={\dot{ϵ}}_{\mathrm{rad}}$$

$${\lambda }_{\mathrm{obs}}={\lambda }_{\mathrm{ij}}\left(1-\frac{z}{ct}\right)$$

$$L_{\lambda }(\lambda )d\lambda =\left[\sum _{\mathrm{i},\mathrm{j}}\frac{ct}{{\lambda }_{\mathrm{ij}}}\int {\dot{ϵ}}_{\mathrm{ij}}(x,y,z)dxdy\right]d\lambda$$

$$\rho(v)=\left\{\begin{array}[]{ll}\rm\rho_{0}\left(\frac{v}{v_{t}}\right)^{-\delta}&\rm~{}v\leq v_{t}\\ \rm\rho_{0}\left(\frac{v}{v_{t}}\right)^{-n}&\rm~{}v>v_{t}\\ \end{array}\right.$$

$${\rho }_0=\frac{{\mathrm{M}}_{\mathrm{ej}}}{4\pi ({\mathrm{v}}_{\mathrm{t}}{\mathrm{t}}_{\mathrm{ex}}{)}^3}{\left[\frac{1}{3-\delta }+\frac{1}{\mathrm{n}-3}\right]}^{-1}$$

$${\mathrm{E}}_{\mathrm{K}}=\frac{1}{2}{\mathrm{M}}_{\mathrm{ej}}{\mathrm{v}}_{\mathrm{t}}^2\left[\frac{1}{5-\delta }+\frac{1}{\mathrm{n}-5}\right]{\left[\frac{1}{3-\delta }+\frac{1}{\mathrm{n}-3}\right]}^{-1}$$

$${\mathrm{v}}_{\mathrm{t}}=10,943{\mathrm{E}}_{51}^{1/2}{\mathrm{M}}_1^{-1/2}\mathrm{km}{\mathrm{s}}^{-1}$$

$${\rho }_0=4.90\times 10^{-17}{\mathrm{E}}_{51}^{-3/2}{\mathrm{M}}_1^{5/2}{\mathrm{t}}_{200}^{-3}\mathrm{g}\mathrm{c}{\mathrm{m}}^{-3}$$

$${\mathrm{f}}_{\gamma ,\mathrm{c}}\approx 0.025{\mathrm{E}}_{51}^{-1}{\mathrm{M}}_1^2{\mathrm{t}}_{200}^{-2}$$

$${\eta }_{\mathrm{h}}=\frac{L_{\mathrm{e}}(E_0)+f{\sum }_{\mathrm{k}}{x}_{\mathrm{k}}{I}_{\mathrm{k}}{Q}_{\mathrm{k}}(E_0)}{A}$$

$${\eta }_{\mathrm{k}}=\frac{x_{\mathrm{k}}{I}_{\mathrm{k}}{Q}_{\mathrm{k}}(E_0)}{A}$$

$${\eta }_{\mathrm{ij}}=\frac{x_{\mathrm{i}}{E}_{\mathrm{ij}}{\sigma }_{\mathrm{ij}}(E_0)}{A}$$

$$A=L_{\mathrm{e}}(E_0)+(1+f)\sum _{\mathrm{k}}{x}_{\mathrm{k}}{I}_{\mathrm{k}}{Q}_{\mathrm{k}}(E_0)+\sum _{\mathrm{ij}}{x}_{\mathrm{ij}}{E}_{\mathrm{ij}}{\sigma }_{\mathrm{ij}}(E_0)$$

$$L_{\mathrm{e}}(E)=-\frac{1}{n}{\left(\frac{dE}{dx}\right)}_{\mathrm{e}}=x_{\mathrm{e}}\frac{2\pi e^4}{E}\ln \left(\frac{4E}{{\zeta }_{\mathrm{e}}}\right)$$

$${\omega }_{\mathrm{p}}(n_{\mathrm{e}})=\sqrt{\frac{4\pi n_{\mathrm{e}}{e}^2}{m_{\mathrm{e}}}}=56414.6\sqrt{\frac{n_{\mathrm{e}}}{\mathrm{c}{\mathrm{m}}^{-3}}}{\mathrm{s}}^{-1}$$

$$Q_{\mathrm{I}}(E)=\frac{1}{u{I}^2}\left[A\left(1-\frac{1}{u}\right)+B{\left(1-\frac{1}{u}\right)}^2+C\ln (u)+D\frac{\ln (u)}{u}\right]\mathrm{c}{\mathrm{m}}^2$$

$${\sigma }_{\mathrm{ij}}=\frac{8\pi }{\sqrt{3}}\frac{1}{k_{\mathrm{i}}^2}\frac{I_{\mathrm{H}}}{\Delta E_{\mathrm{ij}}}{f}_{\mathrm{ij}}\bar{g}\pi a_0^2$$

$$\bar{g}\sim -0.0065{\left(\frac{E}{\Delta E_{\mathrm{ij}}}\right)}^2+0.228\left(\frac{E}{\Delta E_{\mathrm{ij}}}\right)-0.07$$

$$\bar{g}=\frac{\sqrt{3}}{2\pi }\ln \left(\frac{\mathrm{E}}{\mathrm{\Delta }{\mathrm{E}}_{\mathrm{ij}}}\right),E/\Delta {\mathrm{E}}_{\mathrm{ij}}>36$$

$$\bar{g}=0.2,E/\Delta {\mathrm{E}}_{\mathrm{ij}}<\sqrt{2}$$

$${\sigma }_{\mathrm{ij}}=\frac{1}{g_{\mathrm{i}}{k}_{\mathrm{i}}^2}{\Omega }_{\mathrm{ij}}\bar{g}\pi a_0^2$$

$${{\rm Y}}_{\mathrm{ij}}={\int }_0^{\infty }{\Omega }_{\mathrm{ij}}{e}^{-{ϵ}_{\mathrm{j}}/kT}d({ϵ}_{\mathrm{j}}/kT)$$

$$C_{\mathrm{ij}}(T)$$

$$C_{\mathrm{ji}}(T)$$

$$C_{\rm ij}(T)=\begin{cases}\displaystyle 2.16f_{\rm ij}n_{\rm e}\biggl{[}\frac{E_{\rm ij}}{kT}\biggr{]}^{-1.68}T^{-3/2}e^{-E_{\rm ij}/kT}\rm\ s^{-1},~{}(neutral~{}atoms)\rule{0.0pt}{11.19443pt}\\ \displaystyle 3.9f_{\rm ij}n_{\rm e}\biggl{[}\frac{E_{\rm ij}}{kT}\biggr{]}^{-1}T^{-3/2}e^{-E_{\rm ij}/kT}~{}\rm s^{-1},~{}(positive~{}ions)\rule{0.0pt}{11.19443pt}\end{cases}$$

Full text

1. Introduction

Spectra taken of supernovae at late times (in the nebular phase, $\gtrsim$ 100 days after explosion) probe the central regions of the ejecta and thus contain a wealth of information about the explosion, such as nucleosynthetic yields, compositional mixing, and geometry. The quantity and breadth of nebular spectra has grown rapidly in recent years due to international observational efforts. However, further modeling is needed to develop a systematic understanding of how nebular spectra depend on explosion parameters and how atomic data inputs affect spectral modeling.

Nebular spectra can be used to study the uncertain progenitors of Type Ia supernovae (SNe Ia, Sollerman et al., 2004; Maeda et al., 2010a; Mazzali et al., 2011; Mazzali & Hachinger, 2012; Mazzali et al., 2015). While SNe Ia are thought to be the result of the explosion of carbon-oxygen white dwarfs (C/O WDs) in a binary system (Hoyle & Fowler, 1960), their progenitor systems and explosion mechanisms are still unknown. Despite the success of the empirical width-luminosity relation (Phillips, 1993) to calibrate luminosities of SNe Ia, systematic variation due to intrinsic SN Ia diversity is an ongoing challenge for precision cosmology (Howell, 2011).

Several candidates for the progenitors of “normal” SNe Ia have been proposed (see Branch et al., 1993; Wang & Han, 2012; Maoz et al., 2014, and references therein). In the single-degenerate ${M}\rm_{Ch}$ model, the WD gains mass from a non-degenerate binary companion (Whelan & Iben, 1973) and ignites carbon burning when the mass approaches $M_{\rm ch}\approx 1.4M_{\odot}$. In the double-detonation model, a layer of helium gas from a binary companion detonates above the primary C/O WD causing the WD to detonate at a mass ${<M}\rm_{Ch}$ (Livne, 1990; Woosley & Weaver, 1994; Fink et al., 2007; Guillochon et al., 2010; Shen & Bildsten, 2014). In the double-degenerate model, two WDs coalesce or collide and may detonate violently on impact or subsequent to the merger (Iben & Tutukov, 1984; Webbink, 1984; Benz et al., 1989; Raskin et al., 2009; Rosswog et al., 2009; Pakmor et al., 2012; Kushnir et al., 2013; Moll et al., 2014; Kashyap et al., 2015). It may be that SNe Ia come from most, if not all, of these progenitor channels. A main goal of spectral modeling is to help understand the origin of different events.

Nebular spectra reveal emission throughout the entire ejecta and so are a valuable probe of density and compositional structure. If the ejecta mass, kinetic energy, and/or compositional yields can be determined, these can differentiate explosion scenarios. For example, it is thought that ${M}\rm_{Ch}$ models produce more stable IGEs (such as 58Ni and 54Fe) than sub-${M}\rm_{Ch}$ mass models due to burning at higher central densities (Iwamoto et al., 1999). Nebular line profiles are also sensitive to global asymmetries and therefore offer a way to study the imprints of the explosion mechanism on ejecta geometry.

Initial work by Axelrod (1980) enabled a number of codes used to model the nebular spectra of supernovae (Ruiz-Lapuente & Lucy, 1992; Liu et al., 1997; Kozma & Fransson, 1998a, b). Updated codes use improved atomic data and incorporate more sophisticated radiative transport and non-thermal deposition physics (Sollerman et al., 2000; Mazzali et al., 2001; Sollerman et al., 2004; Kozma et al., 2005; Maeda et al., 2006; Mazzali et al., 2007b; Maurer et al., 2010; Dessart & Hillier, 2011; Jerkstrand et al., 2011; Li et al., 2012). Despite the publication of numerous 3D hydrodynamic simulations of SNe Ia (e.g., Hillebrandt et al., 2013), few such explosion models have been analyzed in the nebular phase (but see Kozma et al., 2005). Furthermore, previous nebular modeling work has focused on interpreting events individually (Ruiz-Lapuente & Lucy, 1992; Maeda et al., 2010a, b; Sollerman et al., 2004; Stehle et al., 2005; Maurer et al., 2011; Mazzali et al., 2011; Mazzali & Hachinger, 2012; Taubenberger et al., 2013; Mazzali et al., 2015; Ashall et al., 2016), rather than carrying out parameter studies to identify the general dependencies and degeneracies of nebular spectra.

The aim of this work is to systematically study how variations in explosion properties, density and abundance structure, and atomic data inputs affect the spectra of SNe Ia at late times. To that end, we present a new NLTE code to model the nebular spectra of supernovae. In §2, we present our method of calculating level populations, non-thermal deposition, temperature and ionization balance, and nebular spectra. In §3, we use a fiducial model to describe the physics of nebular spectral formation in SNe Ia. We then vary the parameters of the model to probe the sensitivity of the spectra to ejecta mass, composition, and kinetic energy (§4.1-4.6), as well as density profile (§4.7) and atomic data inputs (§4.8).

2. Methods

We have developed a new 3D radiative transfer tool to model nebular spectra of SNe. Given an initial ejecta model, the code calculates the emissivity of each atomic transition by solving for the temperature, ionization state, and NLTE atomic level populations, including non-thermal effects from radioactive decay products, and generates spectra by integrating the radiative transfer equation in a moving medium.

2.1. Basic Assumptions

The underlying supernova ejecta model is specified by the mass density and elemental abundances in each zone on a 3D Cartesian grid. We assume 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 ($\gtrsim$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 supernovae, this is true in the optical and infrared regions at times $\gtrsim 100$ days. However, at blue and ultraviolet wavelengths ($\lesssim 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, an issue that will be addressed in the future by incorporating a more general transport solver.

2.2. Non-LTE Level Populations

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

[TABLE]

where $\vec{n}$ is a vector of level populations for the species, and the matrix $\mathbf{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

[TABLE]

where ${n}\rm_{I,i}$ is the level population of the i-th level of the ion with I-th ionization state, ${n}\rm_{e}$ is the electron density, ${R}\rm_{I,i}$ is the total ionization rate from the i-th to the I+1-th ionization state (including photoionization, collisional ionization, and non-thermal ionization), and ${\alpha}\rm_{i}$(T) is the total recombination coefficient from ionization state I’s ground state to the i-th level (including radiative, dielectronic, and three-body recombination). ${\mathcal{E}}\rm_{ij}$ 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 eq. 1. For the SN Ia problem, we treat collisional (“electron-impact”) bound-bound and bound-free processes, radiative and dielectronic recombination, and non-thermal deposition of 56Co energy into heating, excitation, and ionization channels. Following Nahar & Pradhan (1997), we assume that all collisional ionizations 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.

2.3. Non-thermal Effects

Energy deposition due to the radioactive decay chain 56Ni$\rightarrow$56Co$\rightarrow$56Fe is the dominant energy source for SNe Ia in the nebular phase (Colgate & McKee, 1969; Kuchner et al., 1994). Since 56Ni decays on a timescale of $\rm\tau_{Ni}=7.7$ days, its contribution to heating is minor by nebular times, while 56Co decay ($\rm\tau_{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., 22Na, 35S, 45Ca, 46Sc, 49V, 54Mn, 55Fe, 57Co, 58Co, 65Zn, and 68Ge) 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 56Co 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 56Co 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 non-thermal electrons/positrons are given by

[TABLE]

[TABLE]

where the transition from level $i$ to level $j$ (with $j>i$) has energy ${\Delta E}\rm_{ij}$, ${n}\rm_{i}$ is the level population of the i-th atomic level, ${I}\rm_{k}$ is the ionization potential of the ion indexed by k having number density ${n}\rm_{k}$, ${\dot{\epsilon}}\rm_{\rm rad}$ is the radioactive energy deposition rate per unit volume, and ${\eta}\rm_{ij}$ and ${\eta}\rm_{k}$ refer to the non-thermal excitation and ionization deposition fractions, respectively, described in appendix B.

2.4. Temperature and Ionization Balance

We use an iterative non-linear 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.

[TABLE]

where ${n}\rm_{k,i}$ represents the population of the i-th ionization stage of the k-th species (i=0 is neutral). Given that the non-thermal deposition fractions (${\eta}\rm_{ij}$ and ${\eta}\rm_{k}$ in eqs. 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 ergs s*-1* cm*-3*) from a transition between the i-th and j-th levels of an ion ($E_{\rm j}>E_{\rm i}$) is

[TABLE]

where ${\nu}\rm_{ij}$ is the photon frequency of the transition, ${A}\rm_{ij}$ is the spontaneous radiative decay rate (Einstein A coefficient), and ${n}\rm_{j}$ is the population of the j-th level of the ion. The temperature sensitivity of $\dot{\epsilon}_{\rm ij}$ arises primarily from the temperature-dependent collisional excitation rates, which are important in setting the populations $n_{\rm 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

[TABLE]

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 non-linear 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 and $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.

2.5. Spectrum Calculation

Once the line emissivities (eq. 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 non-relativistic Doppler effect along the line of sight is

[TABLE]

where $\lambda_{\rm obs}$ is the observed wavelength of a line transition with rest-frame wavelength ${\lambda}\rm_{ij}$. Equation 8 allows us to associate the observed flux at wavelength $\lambda_{\rm 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 $\lesssim$10,000 km/s, using the non-relativistic Doppler shift leads to errors in z-mapping of only $\lesssim 3\%$.

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

[TABLE]

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.

2.6. Physical Processes Neglected

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 SN011fe 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 sub-dominant. Since CT preferentially ionizes neutral atoms (as Coulomb repulsion suppresses ion-ion interactions), we expect that neutral atoms like OI and FeI 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.

2.7. Code Verification

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.

3. Modeling Type Ia nebular spectra

To study the nebular spectra of SNe Ia, we construct spherically symmetric models in which the ejecta properties (e.g., total mass, energy, and abundances) are free parameters. We first describe the general properties of spectrum formation in a “fiducial” model that resembles the ejecta structure expected for normal SNe Ia. In §4 we carry out a parameter survey that demonstrates how the nebular spectra depend on explosion properties.

3.1. Ejecta Modeling

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

[TABLE]

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

[TABLE]

[TABLE]

where ${M}\rm_{ej}$ is the total ejecta mass, $E_{\rm k}$ is the ejecta kinetic energy, and ${t}\rm_{ex}$ is the time since explosion. The radial density profile is thus completely set by the choice of $M_{\rm ej},~{}E_{\rm k}$, and the exponents $\delta$, $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 $\rm\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

[TABLE]

[TABLE]

where $M_{1}=M_{\rm ej}/M_{\odot}$, $E_{51}$ is the kinetic energy in units of $10^{51}$ erg, and ${t}\rm_{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 §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 be consist of stable iron group elements (IGEs) of mass $M_{\rm IGE}$. The stable IGEs are assumed to be composed of a ratio, ${R}\rm_{stb}$, of 54Fe to 58Ni. Surrounding the stable IGE region is a zone consisting (initially) mostly of 56Ni of mass ${M}\rm_{56Ni}$. We include a small amount of stable IGEs in this region with mass abundance ${X}\rm_{stb}$ and the same isotopic ratio ${R}\rm_{stb}$. Above the 56Ni is an outer layer of intermediate mass elements (IMEs) of mass ${M}\rm_{IME}$. We compose the IME layer of 70% 28Si, 29% 32S, and 1% 40Ca, roughly consistent with the nucleosynthetic results in the SNIa 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 §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}\rm_{56Ni}$ 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_{\rm\gamma,c}=\rm 1-e^{-\tau_{\gamma}}\approx\tau_{\gamma}$, where $\tau_{\gamma}$ is the mean optical depth to gamma-rays. Taking a typical gamma-ray opacity $\rm\kappa_{\gamma}\sim 0.03~{}cm^{2}~{}g^{-1}$ (Swartz et al., 1995) and integrating the radial optical depth from the center ($r=0$) gives an estimate

[TABLE]

Eq. 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}\rm_{ej}$, ${E}\rm_{K}$, and $t$ will be useful for interpreting how radioactive deposition rate depends on physical parameters.

3.2. Fiducial Model

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_{\rm ej}$ near the Chandrasekhar mass and $M_{\rm 56Ni}=0.6~{}M_{\odot}$. The fiducial model has a transition velocity of $\rm 10,131~{}km~{}s^{-1}$ while the interface between the nickel zone and IME layer is near $\rm 8,800~{}km~{}s^{-1}$.

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 Type Ia SN2011fe (Nugent et al., 2011; Shappee et al., 2013; Mazzali et al., 2015). The model spectrum reproduces most of the prominent features.

Figure 3 shows a breakdown of the contribution from various ions to the fiducial model spectrum. The strongest features are due to forbidden transitions of FeII and FeIII which 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.

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 2 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.

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 non-thermal electrons produced by radioactive decay (the collisional ionization rate from thermal electrons is sub-dominant), which is balanced by radiative recombination. Figure 2 shows the radial dependence of the ionization fractions of iron for the fiducial model. Though FeIV and FeV 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 FeIII and FeII. In general, the degree of ionization increases as the density declines, reflecting the reduced rate of radiative recombination.

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 [SII] and [FeII] transitions. Similar to other studies (e.g. Mazzali et al., 2015), 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.

4500-5500 Å This region is dominated by emission from [FeIII], which produces prominent features at 4658 Å and 5270 Å. The latter feature also includes a significant contribution from [FeII] 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 [FeIII] 5011 Å, and can be washed out if the ejecta velocities are too high.

5500-7000 Å This region contains a prominent [CoIII] feature near 5888 Å with two smaller [CoIII] features immediately to the red. The line strength depends on the abundance and ionization state of Cobalt, and typically declines with time as 56Co decays to 56Fe. 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 non-degenerate companion should have velocities $\rm\lesssim~{}1000~{}km~{}s^{-1}$ (e.g. Marietta et al., 2000), this emission is unlikely to be due to H$\alpha$. 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 [FeII] lines and a broad [CaII] line. In addition, [NiII] emission can contribute if sufficient stable 58Ni is present in the gas, but for our fiducial model [NiII] does not dominate this feature. Notably, the fiducial model only contains 0.0075$\rm M_{\odot}$ of calcium and yet [CaII] emission from the IME layer dominates the emission in this region.

8500-11000 Å This region is made up primarily of [SII] and [SIII] emitted by the 0.22$\rm 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.

4. Systematic Parameter Study

We present the following systematic study probing the sensitivity of synthetic nebular spectra to model parameters and atomic data inputs. In particular, we calculate how the spectra depend on the total ejecta mass, kinetic energy, 56Ni mass, and time since explosion, as well as density and compositional structure. A full list of models and their derived properties is shown in Table 3.

4.1. Time since explosion

Figure 4 shows the nebular spectrum of the fiducial model at times between 200 and 400 days after explosion. As the supernova evolves in time, ejecta densities decline and radioactive isotopes decay. The bolometric luminosity drops with time due to the declining radioactive heating and a decreasing gamma-ray trapping fraction. The relative strength of features in the spectrum, however, remains fairly constant, as seen in the figure inset. The ratio of the prominent Fe lines does not evolve significantly, indicating that the FeIII/FeII ionization ratio remains fairly constant over time.

The biggest change in the spectral features over time is the strength of the [CoIII] emission at 5888 Å, which decreases as 56Co decays to 56Fe. In addition, the IME emission of [SII] 4000 Å declines at later times as a result of a declining sulfur ionization fraction in the IME layer.

4.2. Explosion kinetic energy

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

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

Changing ${E}\rm_{K}$ has a modest effect on the shape of the spectral features. An increase of ${E}\rm_{K}$ from 1 B 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 B to 2 B) does have noticeable effects, causing the FeIII/FeII 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.

4.3. Total ejecta mass

Theoretical models of SNe Ia predict an ejected mass in the range $0.8-2.0$ $\rm 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 6 shows synthetic spectra of the fiducial model with ${M}\rm_{ej}$ varied between $1-2$ $\rm M_{\odot}$. For a fixed kinetic energy, a higher ${M}\rm_{ej}$ results in higher ejecta densities and lower velocities. As a result, higher ${M}\rm_{ej}$ models have greater gamma-ray trapping, a brighter bolometric luminosity, and less Doppler broadened spectral features. This effect of ${M}\rm_{ej}$ is therefore somewhat degenerate with that of kinetic energy.

The relative features in the synthetic spectra of Figure 6 show only subtle variations with ${M}\rm_{ej}$. For our super-${M}\rm_{Ch}$ case with $M_{\rm ej}=2.0~{}M_{\odot}$, the IME features become visibly stronger due to the higher total IME mass. The FeIII/FeII 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 6, this parameter is likely correlated with kinetic energy and 56Ni mass (Woosley et al., 2007); for example, a super-${M}\rm_{Ch}$ explosion is likely to produce more 56Ni and higher kinetic energy. We therefore ran three additional models in which we fixed the ratios $\rm M_{56Ni}/M_{ej}=1/2$ and $\rm E_{51}/M_{ej}=1.2/1.4$. Based on the above discussion of gamma-ray trapping (§3.1), we expect luminosity to approximately scale as $\rm L_{bol}\sim M_{ej}^{2}$ if these ratios are held constant.

Figure 7 shows synthetic spectra of these scaled model with total ejecta masses 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}\rm_{Ch}$ models, provided ${E}\rm_{k}$ and ${M}\rm_{56Ni}$ scale accordingly.

4.4. Radioactive nickel mass

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

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

4.5. Neutron-rich Iron Group Elements (IGEs)

Stable IGEs are produced in SN Ia ejecta in two distinct ways. The trace presence of neutron-rich isotopes (in particular 22Ne) due to the metallicity of the progenitor WD and pre-explosion carbon simmering leads to the production of up to $\approx 25$% by mass of neutronized IGEs (in particular 54Fe and 58Ni) 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 $\sim 0.05-0.4$ $\rm M_{\odot}$ of neutronized IGEs (Nomoto et al., 1984; Thielemann et al., 1986; Seitenzahl et al., 2011, 2013). This latter effect only occurs in white dwarfs with $\rm M_{ej}\gtrsim M_{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}\rm_{Ch}$ models often predict stable IGEs to be produced at the core (Nomoto et al., 1984; Mazzali et al., 2007a), while multi-D simulations indicate that buoyancy should mix IGEs throughout the 56Ni region (Gamezo et al., 2005; Kasen et al., 2009; Seitenzahl et al., 2011). As a result, nebular spectra indicators of stable IGEs would be valuable for inferring the progenitor scenario.

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

Higher stable core mass also increases the relative fraction of FeII compared to FeIII, which is a result of lowered non-thermal ionization in the non-radioactive core region. Another noteworthy product of a stable core is the flat-topped profile of the [CoIII] 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}\rm_{stb}$ models by construction.

A comparison between two ${M}\rm_{stb}$ models is shown in Figure 10. In the $M_{\rm stb}=0.2M_{\odot}$ model we find an apparent blueshift of the [NiII] 7378 Å peak by $\approx 1500$ km/s due to the relative blending between FeII and NiII components. The feature is also sensitive to [CaII] emission, which is stronger in the low-${M}\rm_{stb}$ model. The [FeII] feature redward of [NiII] 7378 Å may also dominate in some models, resulting in an apparent redshift of the peak.

Maeda et al. (2010b) 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 [NiII] emission in this feature can been attempted (Maeda et al., 2010a) but is subject to model-dependent uncertainties.

Surprisingly, the mixing of stable Nickel through the nickel zone (i.e., varying the parameter ${X}\rm_{stb}$ between $0.05-0.20$) does not produce a visible effect on the synthetic spectrum, even for 58Ni masses comparable to those in the core stable IGE region of figure 9. This is due to the high level of ionization of Ni in the nickel zone, which suppresses [NiII] emission. For the same reason, the ratio of 54Fe to 58Ni (${R}\rm_{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 $\rm M_{\odot}$ of stable nickel can be hidden in the ejecta of a SN Ia model with fiducial parameters due to this ionization effect.

4.6. Mass of Carbon/Oxygen in ejecta

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 SN2011fe 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 [OI] features for a pure deflagration model of Röpke (2005), which contained 0.6$\rm M_{\odot}$ of unburned C/O material. There have also been detections of possible [OI] emission in sub-luminous SNe Ia (Taubenberger et al., 2013; Kromer et al., 2013).

Figure 11 shows synthetic spectra of the fiducial model with C/O mass varied between $0.1-0.4$ $\rm 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).

We find a strong contribution of [OIII] at 5007 Å as well a weak contribution of [OI] 6300/6364 Å and a blend of [OII] features at 7320 Å. The high ionization state of oxygen prevents the [OI] emission seen in Kozma et al. (2005) to contribute significantly, but we expect [OI] to become stronger if oxygen were more concentrated in the higher density central regions, or if the mass of 56Ni (and hence radioactive deposition) were lower. C/O within the IME layer does not produce a significant nebular feature.

4.7. Sensitivity to density profile

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). Here, we attempt to illuminate the effect of steepening our interior power-law density profiles. We also consider an exponential profile of the form $\rm\rho=\rho_{0}e^{-v/v_{e}}$ where $\rho_{0}$ and ${v}\rm_{e}$ can be determined from ${M}\rm_{ej}$, ${E}\rm_{K}$, and ${t}\rm_{ex}$.

Figure 12 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 Å. The two peaks in this model are dominated by [FeII] and [CaII], with little contribution from [NiII].

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 FeIII/FeII line ratio.

When comparing to the nebular spectrum of SN2011fe, 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 FeII to FeIII ratio (indicated by high flux in the [FeII] 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.

4.8. Sensitivity to atomic data uncertainties

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 ($\rm Q_{k}$), the collisional excitation rates ($\rm C_{ij}$), and the radiative recombination rates ($\alpha$).

We show the impact of systematically varying atomic data in Figures 13-14. To explore the effect of the collisional ionization cross sections (Figure 13), 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 FeII was decreased by a factor of 2 and that of FeIII was increased by a factor of 2, and (d) the cross section of FeII was increased by a factor of 2 and that of FeIII was decreased by a factor of 2. These variations in $Q$ affect the ionization ratio of FeIII/FeII and can modify the FeIII/FeII line ratio by up to $\pm$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.

To explore the effect of the collisional excitation rates (Figure 14), we carried out calculations in which all ${\sigma}\rm_{ij}$ were increased or decreased by a factor of 2. These variations in ${\sigma}\rm_{ij}$ affect the strength of emission features and can modify the FeIII/FeII line ratio by up to $\pm$10%.

We show in Figure 15 the fiducial model calculated using recombination rates from two different databases. Our calculations throughout has used recombination rates from the CHIANTI database, but more recent data for FeI-V are available from the Nahar OSU Radiative Atomic Database (NORAD) (Nahar, 1996; Nahar et al., 1997; Nahar, 1997; Nahar et al., 1998). With this latter data set, we have access to state-specific recombination rates, which we neglected in the above treatment. We also neglected charge transfer 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 charge transfer in FeI-IV with data from Krstic et al. (1997).

While there is reasonable agreement between the synthetic spectra, the strength of emission features does depend on the atomic data set. In particular, the FeIII/FeII line ratio shows a higher FeII population when NORAD recombination rates are used, due to higher total recombination rates for FeI-III in the NORAD data set, leading to a $\sim$30% decrease in FeIII/FeII line ratio. Similar variations are seen in most other spectral features.

5. DISCUSSION and conclusion

We have presented a new tool for calculating nebular spectra of SNe and applied it to a systematic parameter study of SNe Ia. We summarize some of our main results as follows:

Robustness and Degeneracy of Spectral Features: On the whole, we found that the features in SN Ia nebular spectra were remarkably insensitive when physical parameters were changed within the range expected for white dwarf explosions. Individually changing the ejected mass, 56Ni mass, or kinetic energy by 50% produced only minor changes in the relative strength of features (although the bolometric luminosity was affected due to changes in gamma-ray trapping). Varying parameters by a larger factor of $\approx 2$ began to show noticeable spectral changes.

In addition, we found degeneracies in the effect of different physical parameters. An increase in the total ejecta mass, for example, could be mostly offset by a corresponding increase in kinetic energy so as to keep the overall density and velocity scale roughly fixed. As a result, quite different sets of physical parameters may be able to fit the same nebular spectrum.

Progenitor Mass: Though some previous studies have favored Chandrasekhar mass models in explaining SN Ia nebular spectra, our model survey demonstrates that observed nebular spectra are equally well fit by generic sub-Chandrasekhar mass, Chandrasekhar mass, and super-Chandrasekhar mass models. The spectral features remained essentially unchanged when the total mass, 56Ni mass, and kinetic energy were scaled up or down in unison (Figure 7). Nebular spectra therefore do not alone constrain the overall mass scale of SNe Ia, and instead are more diagnostic of the relative abundance yields and the density profile.

Stable Iron Group Elements: Previous studies have fit the nebular spectra of SNe Ia with models having a $\sim 0.1-0.2~{}M_{\odot}$ sphere of stable IGEs at the ejecta center. The presence of such an IGE “core” would favor a Chandrasekhar mass model (where electron capture occurs during high-density deflagration burning), while disfavoring double-detonation and violent merger models.

A stable IGE “core” has been claimed necessary to explain the FeIII/FeII line ratio in observed spectra (Maeda et al., 2010a; Mazzali et al., 2015). We show that this line ratio also depends on other explosion parameters, such as the density profile and kinetic energy, and is sensitive to the atomic data inputs. We present models that well fit the SN 2011fe iron line ratio well without invoking any stable IGE “core.”

Our models that did include a large ($\gtrsim 0.1~{}M_{\odot}$) stable IGE “core” produced a flat-topped [CoIII] feature near 5900 Å. Such a flat-top [CoIII] profile is not seen in the nebular spectra of SN 2011fe (Mazzali et al., 2015) indicating that at least some radioactive cobalt exists at the lowest ejecta velocities. Höflich et al. (2004) and Motohara et al. (2006) observe relatively flat-topped profiles of the 1.64 $\mu$m [FeII] feature in SN 2003du and SN 2003hv, respectively, which they take as evidence that SNe Ia have a stable IGE “core” lacking any radioactive heating. The lack of an obvious flat-topped optical [CoIII] feature in SN 2011fe appears to provide evidence to the contrary.

Another diagnostic of stable IGE is the feature at $\sim 7300$ Å usually attributed to [NiII]. Similar to (Mazzali et al., 2015), we find that a small amount ($\sim 0.01~{}M_{\odot}$) of central 58Ni is sufficient to reproduce a peak at 7300 Å. However, we note that in some models the [FeII] 7388 Å line can alone account for this peak, while in general [FeII] 7155 Å and [CaII] emission also shape the overall line profile. The steepness of the inner ejecta density profile also affects how narrow and resolved the separate line peaks of the 7300 Å feature are. We find that additional 58Ni mixed throughout the lower density layers of radioactive 56Ni zone may be too ionized to produce significant [NiII] 7378 Å emission. As a result, we consider it difficult derive a precise 58Ni mass constraint from analysis of the 7300 Å feature.

In sum, our parameter study does not provide strong evidence that SNe Ia possess a substantial stable IGE “core,” nor do we see a robust way of accurately inferring the total stable IGE mass from nebular spectrum analysis. We therefore do not consider the nebular spectra of SNe Ia as providing particularly strong support for the Chandrasekhar mass model over other progenitor scenarios. However, our study did not focus on infrared features (Höflich et al., 2004; Motohara et al., 2006; Gerardy et al., 2007; Telesco et al., 2015) or the effects of positron transport (Penney & Hoeflich, 2014). Further study of stable IGE signatures using specific multi-dimensional explosion models is warranted.

IME indicators: While the nebular spectra of SNe Ia are dominated by iron lines, we emphasize that some features provide constraints on IME abundances. In particular, the feature at 7300 Å is sensitive to [CaII] emission at 7291 Å, while the feature at 9500 Å is dominated by [SIII] emission. We showed that the relative strength of the [SIII] 9500 Å feature to the [FeIII] 4664 Å feature tracked the relative abundance of IGE to IME (see Figure 8). Analysis of this line ratio may thus provide valuable constraints on the nucleosynthetic yields and hence progenitors of SNe Ia. Our successful fit to the calcium and sulfur features in SN 2011fe suggests that the total IME yields may be $\approx 0.75~{}M_{\odot}$.

Line Shifts and Ejecta Geometry: Shifts in the wavelength of nebular line peaks are seen in many SNe Ia (Black et al., 2016) and have been used to deduce asymmetries of, e.g., iron group elements (Maeda et al., 2010a). Though we have only considered spherically symmetric models in this paper, we nevertheless see shifts in the location of emission peaks due to line blending. For example, the feature around 7300 Å is a blend of [NiII], [FeII], and [CaII] lines. Depending on the ejecta compositional structure and time-evolution of line strengths, the location of the composite peak can vary by 1500 km/s in different models. This highlights the difficulty in using the 7300 Å feature shift to derive reliable kinematic measures of asymmetry (though see Maeda et al. (2010a) for attempts separate components in the blend). Our models do not show a strong, progressive redshift evolution in the main IGE features as a function of time. This supports the claim of Black et al. (2016) that these shifts are due to permitted line emission that may help shape the spectrum even at times $\gtrsim 200$ days after explosion.

Carbon/Oxygen mixing: We find that C/O material mixed into a nickel zone is highly ionized. For large enough C/O masses ($\gtrsim 0.1M_{\odot})$ this produces a visible [OIII] feature at 5007 Å, which is not seen in SN 2011fe. This disfavors models leaving significant amounts of oxygen mixed throughout the ejecta. To produce a narrow [OI] feature like that seen by Taubenberger et al. (2013) in the sub-luminous SN 2010lp would presumably require the C/O to to be located in a dense central region with less non-thermal ionization from radioactivity. C/O in the outermost ejecta layers (within or above the IME zone) does not experience significant radioactive heating and has no impact on the nebular spectra.

Density Profile: We have shown that ejecta having a steep interior density profile produce narrower and more central peaked nebular emission features. In fact, the line widths depended more on the density profile than on the kinetic energy, since the former sets the degree of central concentration of mass. Comparing to the spectra of SN 2011fe, we find that a flat density profile (power-law index $\delta=0-1$) best reproduces the observed spectral profiles. This suggests that nebular spectra may be able to test the density structures predicted by detailed SN Ia explosion models.

Atomic Data: Though it is well known that synthetic nebular spectra are affected by uncertainty and incompleteness in the input atomic data, the level of error has not been well quantified. We showed here that factor of 2 errors in some of the key atomic data inputs can result in changes in spectral features at the $\approx 25\%$ level. This level of variation is similar to the level we found when varying the physical ejecta parameters. Given the uncertainties, we advise against over-interpreting model fits to extract quantitative mass estimates from individual observed SNe Ia. The trends seen among a sample of nebular spectra, however, are likely less affected by the atomic data limitations.

Limitations: While the spectral models presented here included a broad range of the most important and complex atomic processes, other physical effects may be relevant and should be addressed in future work. While the ejecta are mostly optically thin at optical wavelengths, there is evidence that some lines are nevertheless optically thick, especially in the ultraviolet (Friesen et al., 2017). Furthermore, dust/molecule absorption and emission may be important at some phases. At very late times, the freezing out of ionization requires that the problem be treated time-dependently (Kozma & Fransson, 1998a; Jerkstrand et al., 2015). The collection of a reliable and complete atomic data sets has been and remains a substantial challenge for all nebular phase modeling efforts.

Finally, while this paper was restricted to 1D parameterized models, the asphericity found in more realistic explosion models, and inferred from polarization observation of SNe, is expected to affect the nebular spectra. Our NLTE nebular spectrum code is 3D, and future work on SNe Ia and other types of SNe will aim to explore the predictions of 3D hydrodynamical models that have complex abundance distributions, clumping on multiple scales, and global ejecta asymmetries.

This material is based upon work supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE 1106400. DK is supported in part by a Department of Energy Office of Nuclear Physics Early Career Award, and by the Director, Office of Energy Research, Office of High Energy and Nuclear Physics, Divisions of Nuclear Physics, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

Appendix A Atomic Data

We use the atomic data compilation of the cmfgen code (Hillier & Miller, 1998; Li et al., 2012) for radiative transition probabilities (Einstein A Coefficients), effective collisional strengths $\Upsilon_{\rm ij}$ (see appendix C.1 for treatment of missing data), and photoionization cross sections. We take ground-state ionization energies from NIST (Kramida et al., 2016). For thermal collisional ionization rates, we use the results of Voronov (1997). For collisional ionization cross sections, we use Mattioli et al. (2007) and references therein. Radiative and di-electronic recombination rates are taken from the CHIANTI database v.8 (Dere et al., 1997; Del Zanna et al., 2015).

Appendix B Nonthermal Deposition Fractions

Following the Continuous Slowing Down Approximation formulation of Xu & McCray (1991), we can write the fraction of deposition into the three channels heating ($\eta_{\rm h}$), ionization ($\eta_{\rm k}$), and excitation ($\eta_{ij}$) as such

[TABLE]

[TABLE]

where $E_{0}$ is the injection energy of electrons, $L_{e}$ is the Coulomb loss function (see below), k labels an ion and its ionization potential $I_{\rm k}$, ionization fraction $x_{\rm k}=n_{\rm k}/n$, and collisional ionization cross section $Q_{\rm k}(E)$, and $f\sim 0.3$ is the mean energy of a secondary electron, and

[TABLE]

where ij labels a transition with energy $E_{\rm ij}=E_{\rm j}-E_{\rm i}>0$, the fractional level population of the lower-energy state is $x_{\rm i}=n_{\rm i}/n$, and the collisional excitation cross section is $\sigma_{\rm ij}$, and

[TABLE]

is the normalization constant. The Coulomb loss function (Schunk & Hays, 1971; Xu & McCray, 1991) for high-energy electrons is

[TABLE]

where $E$ is the electron energy, $e$ is the electron charge, $x_{\rm e}=n_{\rm e}/n$ is the electron fraction, and $\zeta_{\rm e}=2\hbar\omega_{\rm p}$ for plasma frequency given by

[TABLE]

For collisional ionization cross sections, we use the approximate form (Younger, 1981)

[TABLE]

where $u=E/I$, I is the ionization potential, and A, B, C, and D are fitting coefficients. We use the fitting coefficients published by Arnaud & Rothenflug (1985), taking into consideration updates by Arnaud & Raymond (1992) and Mattioli et al. (2007). While more recent calculations are available, Dere (2007) points out that the discrepancies between new and old ionization rates are minor except for NiV-XI, which are not important for this work. Coefficients are unavailable for CoII-IV, so we use FeII-IV cross sections as instead and plan to update this data when it becomes available in the future.

For collisional excitation cross sections of allowed transitions by non-thermal electrons, we use the approximation of van Regemorter (1962),

[TABLE]

where $I_{\rm H}$ is the ionization potential of Hydrogen, $k_{\rm i}^{2}$ is the initial electron energy scaled to 13.60 eV, $f_{\rm ij}$ is the oscillator strength, $a_{0}$ is the Bohr radius, and $\bar{g}$ is given by a quadratic fit to the results of van Regemorter (1962)

[TABLE]

with an imposed high-energy limit of

[TABLE]

for all ions, and a low-energy limit for positive ions

[TABLE]

Forbidden transitions must be treated separately. We use the following form for the collision cross sections:

[TABLE]

where $g_{i}$ is the statistical weight and $\Omega_{ij}$ is the collision strength given by equation C4.

Appendix C Implementation of Atomic Processes

C.1. Collisional excitation and de-excitation

To determine collisional excitation and de-excitation rates, we assume electrons adhere to a Maxwell-Boltzmann distribution. The subsequent “effective” or “Maxwellian-averaged” collisional strength can be calculated for each transition using

[TABLE]

where $\Omega{\rm ij}$ is the collisional cross section for an electron to excite the i-j transition. These collision strengths can be found in the literature and data tables.

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

[TABLE]

where $E_{ij}$ is the energy difference between the two levels and $j>i$.

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

[TABLE]

where $f_{\rm ij}$ is oscillator strength, $T$ is Temperature in Kelvin, $E_{\rm ij}$ is transition energy in eV, and $n_{\rm e}$ is electron density in $\rm cm^{-3}$.

In practice, atomic data for $\Upsilon_{\rm 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 C4 so that there is no discontinuity in $C_{\rm ij}$.

Forbidden transitions for which atomic data is 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_{\rm ij}\leq 0.001$):

[TABLE]

C.2. Recombination

Total recombination rates of relevant ions include a sum over all sub-shells and incorporate both radiative and dielectronic recombination. We assume 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 this latter dataset in §4.8 and plan to implement state-specific recombinations rates in any future work where recombination lines might contribute significantly.