Non-thermal excitation and ionization in supernovae

Non-thermal excitation and ionization in supernovae


We incorporate non-thermal excitation and ionization processes arising from non-thermal electrons that result from -ray energy deposition, into our radiative transfer code cmfgen. The non-thermal electron distribution is obtained by solving the Spencer-Fano equation using the procedure of Kozma & Fransson (1992). We applied the non-thermal calculations to the blue supergiant explosion model whose early evolution was studied in Dessart & Hillier (2010). Non-thermal processes generally increase excitation and ionization and decrease the temperature of the ejecta. We confirm that non-thermal processes are crucial for modeling the nebular spectra. Both optical H i and He i lines are significantly strengthened. While optical He i lines are not easily discerned in observational spectra due to severe blending with other lines, He i 2.058 m provides an excellent opportunity to infer the influence of non-thermal processes. We also discuss the processes controlling the formation of the He i lines during the nebular epoch. Most lines of other species are only slightly affected. We also show that the inclusion of Fe i has substantial line-blanketing effects on the optical spectra. Our model spectra and synthetic light curves are compared to the observations of SN 1987A. The spectral evolution shows broad agreement with the observations, especially H. The uncertainties of the non-thermal solver are studied, and are expected to be small. With this new addition of non-thermal effects in cmfgen, we now treat all known important processes controlling the radiative transfer of a supernova ejecta, whatever the type and the epoch.

radiation transfer – radiation mechanisms: non-thermal – stars: atmospheres – stars: supernovae – stars: evolution

1 Introduction

Supernovae (SNe) explosions are classified as thermonuclear SNe and Core-Collapse SNe (CCSNe), according to their explosive mechanisms. Thermonuclear SNe are also called Type Ia SNe, with light curves that can be standardized to measure cosmological distances. Based on their spectroscopic signatures, CCSNe can be further categorized into 3 broad classes, Type Ib, Type Ic and Type II. Type II SNe exhibit hydrogen lines, while Type I SNe do not.

SNe explosively produce which decays into with a half-life of 6.1 days. further decays into with a half-life of 77 days. Such radioactive decay is the main energy source of Type Ia SNe. While the energy source of Type Ib/c SNe resembles that of Type Ia SNe, Type II SNe are first powered by the shock-deposited energy, but at later times radioactive decay will be the dominant energy source. High-energy -rays are produced during the radioactive decay. The -rays undergo Compton scattering with bound and free electrons, which results in high-energy electrons. These high-energy electrons interact with the medium and deposit energy through heating, excitation, and ionization.

Non-thermal excitation and ionization processes have been suggested to explain some observational features of various SNe of different types. For example, the Bochum event in the H profiles of SN 1987A was explained by non-thermal excitation due to an asymmetric ejection of (Phillips & Heathcote, 1989; Chugai, 1991). Li & McCray (1995) included two-photon processes in models for SN 1987A, and showed that they could fit the evolution of the He i 1.083 m and 2.058 m lines whose upper levels are populated by non-thermal processes induced by -rays.

Several groups (Harkness et al., 1987; Lucy, 1991; Swartz, 1991; Kozma & Fransson, 1992; Swartz et al., 1993; Kozma & Fransson, 1998a, b) studied He i lines in Type Ib SNe and suggested these lines are due to non-thermal excitation. Harkness et al. (1987) carried out LTE calculations for Type Ib SNe and showed that these did not match observations – emission line strengths indicate large departures from LTE for He i. Lucy (1991) treated non-thermal excitation and ionization for He i in their LTE models, and found that it could explain the presence and strength of He i lines. However, these studies did not treat all atoms in non-LTE and/or ignored effects due to radiative transfer. Recently, Dessart et al. (2011) carried out a comprehensive study of SNe which result from quasi-hydrogen-less low-mass cores from 15-25 M main-sequence stars in binary systems using cmfgen. They found that the He i lines could be reproduced at early times if the He mass fraction is greater than 50 per cent, without invoking non-thermal excitation.

Line formation in SNe is inherently complex. The fast expansion velocities exacerbate the influence of lines, by producing broader lines, increasing line overlap, and enhancing line blanketing. Line blanketing strongly influences the optical and UV ranges of SNe spectra, generally producing a flux deficit in those regions. To accurately account for the influence of line blanketing requires that an enormous number of metal lines be included in the model. Because of non-LTE and time dependence, an accurate interpretation of the spectral signatures requires quantitative studies.

For our studies we utilized cmfgen, which is a radiative transfer code developed for modeling Wolf-Rayet stars (Hillier, 1987) that was later enhanced to allow for line blanketing (Hillier & Miller, 1998), to model all early-type stars, and ultimately to model all types of SNe (Hillier & Dessart, 2012). The code includes the following processes: bound-free; bound-bound (including two-photon emission); free-free; dielectronic recombination; electron scattering; Rayleigh scattering by hydrogen; inner-shell ionization by X-rays; collisional ionization/recombination; collisional excitation/de-excitation and charge exchange reactions with H and He. The code does full-non-LTE calculations at great computational cost. Atomic data utilized by cmfgen are listed by Dessart & Hillier (2010).

Time dependent terms in the rate equations are usually neglected in SNe spectral modeling at early times, because the recombination time is much smaller than the age of the SN. However, Utrobin & Chugai (2005) showed that time-dependent ionization of hydrogen is responsible for reproducing the strong H lines at early times in SN 1987A, which was confirmed by Dessart & Hillier (2008). cmfgen takes into account the time-dependent terms in the statistical and radiative equilibrium equations, and in the radiative transfer equations. The code has been extensively tested and used to study various types of SNe (Dessart & Hillier, 2005, 2006, 2008, 2010, 2011; Dessart et al., 2008; Hillier & Dessart, 2012).

At early times, non-thermal effects in Type II-P SNe are small, since the optical depth is large and the ejecta is essentially in LTE in the region where is abundant. However, non-thermal effects will eventually take over other processes, especially during the nebular phase. Further, the study of Jerkstrand et al. (2011) shows that radiative transfer effects are still very important at very late times. In Type Ib SNe, non-thermal processes can become important as early as 10 days (e.g., Lucy, 1991; Dessart et al., 2011, 2012). Thus non-thermal excitation and ionization, as well as radiative transfer effects, need to be taken into account in order to undertake comprehensive studies of SNe from the photospheric through the nebular phase.

With the implementation of non-thermal processes in cmfgen, we are now able to model from just after the explosion, through the photospheric phase, the onset of the nebular phase, and the nebular phase using a single code. The aim of this paper is to describe the implementation of non-thermal excitation and ionization into cmfgen, and to benchmark the non-thermal solver using SN 1987A. The paper is organized as follows: We present our methodology of the non-thermal solver in Section 2, and the hydrodynamical model, utilized in our studies, is summarized in Section 3. Non-thermal effects are discussed in Section 4, while we compare the properties of the non-thermal and thermal models in Section 5. To examine the importance of additional opacity sources we study the influence of Fe i in Section 6. Comparison to the observations is made in Section 7. We explore the uncertainties underlying the non-thermal solver in Section 8, and discuss the effect of mixing and -ray transport in Section 9. Finally, we will summarize our results in Section 10. A parallel paper will discuss the influence of non-thermal processes on Type Ibc spectra (Dessart et al., 2012).

2 The Method

Non-thermal excitation and ionization enter the radiative transfer calculations by changing the level populations. To quantify these excitation and ionization rates, the degradation spectrum, which describes the number of fast electrons as a function of energy, needs to be known. There are several ways to solve for the degradation spectrum (Spencer & Fano, 1954; Shull & van Steenberg, 1985; Fransson & Chevalier, 1989; Xu & McCray, 1991; Kozma & Fransson, 1992). One approach is the Monte Carlo method which, in general, is slow. Our approach is to solve the Spencer-Fano equation following the method of Kozma & Fransson (1992). Once the degradation spectrum is known, all non-thermal excitation and ionization rates are computed and included in the statistical equilibrium equations (SEE).

2.1 The Spencer-Fano Equation

The radiative transfer code cmfgen, which takes into account radioactive decay, previously assumed that all radioactive decay energy is deposited as heat. In reality, fast electrons resulting from Compton scattering between -ray and thermal electrons lose energy through three different channels – heating, excitation and ionization. The fractions of energy for the three channels can be computed if the degradation spectrum is known. This can be done by solving the Spencer-Fano equation (Spencer & Fano, 1954). Following Kozma & Fransson (1992), the Spencer-Fano equation, which balances the number of electrons leaving and entering an energy interval, can be written as:


where , is the maximum energy of non-thermal electrons, and is the ionization potential. is defined as the flux of non-thermal electrons in the energy interval . is the number density of a single species, and the sum over allows for the contribution of different ionization stages and species. is the energy loss of the non-thermal electrons due to Coulomb scattering by the thermal electrons. is the excitation energy, is the excitation cross section of level and is the differential ionization cross section for a fast electron with energy and an energy loss of . Therefore, the ejected electron has energy . The fast electron is called the primary electron and the ejected electron is called the secondary electron. The secondary electron can further degrade through the three same channels as the primary electron. is the source term, and we assume all non-thermal electrons are injected at . The degradation spectrum is normalized to , so that is the degradation spectrum for 1 eV energy input. The first and second term on the left side represent the number of electrons leaving the energy interval by impact excitation and ionization, respectively. The term with is the number of electrons leaving the energy interval by heating the medium. The first term on the right side is the number of electrons entering the energy interval by exciting another ion and losing energy. The second term is the ionization contribution from the primary electrons, and the third term is the ionization contribution from the secondary electrons if applicable.

The Spencer-Fano equation is more conveniently solved numerically using the integral form (Kozma & Fransson, 1992), i.e.,


In Eq. 1, the integral of the differential ionization cross sections needs to be done numerically. However, only the total ionization cross sections appear in the integral form of the Spencer-Fano equation (Eq. 2), and can be computed analytically.

With the solution of the degradation spectrum, the fractional energies entering heating, , excitation of ion to level , , and ionization of ion , , can be computed by




(Kozma & Fransson, 1992). is the mean energy of the initial electrons, and is the lowest excitation or ionization energy.

To save computational effort, we only take into account the dominant ionization stages when computing the degradation spectrum. The other approximation is the neglect of -ray transport. In reality, some -rays may escape the SN, while others may deposit their energy at large distances from where they were created. In these calculations we assume, for simplicity, in-situ deposition. A -ray transport code has been developed for use with cmfgen (Hillier & Dessart, 2012).

We previously mentioned that for the source term we assume that all high energy electrons are injected at . Although not entirely realistic, this kind of source function generally has little influence on the results. High energy electrons “forget” their initial energies after several scatterings, and the source function is quickly washed out. Alternatively we can assume a bell shape source function near . This tends to give a smoother degradation spectrum near .

2.2 Ionization cross sections

Arnaud & Rothenflug (1985) adopted the formula proposed by Younger (1981) for calculating the impact ionization cross sections


where . E is the energy of the impact electron, and is the ionization potential for the element. , , , and are all coefficients. For the meaning of these coefficients, the reader should refer to Arnaud & Rothenflug (1985). In our calculations, we have used the most up-to-date coefficients (Arnaud, private communication).

With the total cross section computed by Eq. 7, the differential cross section can be written as,


where is the distribution of the secondary electrons for a primary electron energy of . Opal et al. (1971) found could be well described by


where is a parameter that gives a best fit to measurements and varies for different species. We adopt , which is shown to be a reasonable approximation (Opal et al., 1971). The secondary distribution function drops quickly so that very few secondary electrons have high energies. Specifically, for H i with  eV, 66.7% of secondary electrons are below the ionization potential if  eV and 69.8% in the case of 10 000 eV. Further, the mean energy of the secondary electrons, , only varies slowly as the energy of the primary electron increases. For hydrogen,  eV for  eV and  eV for 10 000 eV. Secondary electrons with low energies are easily thermalized and deposit energy as heat.

2.3 Excitation cross sections

The electron impact excitation cross sections are only available for a few species, and only for a few levels. To proceed further, we adopt the Bethe approximation as discussed by van Regemorter (1962). With this approximation, is given by


where is the initial electron kinetic energy normalized to hydrogen ionization energy 13.6 eV, is the hydrogen ionization energy, is the energy difference between the two levels, is the oscillator strength, and is the Bohr radius. The Gaunt factor for neutral atoms utilizes Table 1 of van Regemorter (1962), which is shown in Fig. 1. The variable is defined by , where is the kinetic energy of the impact electron and is the transition energy. For positive ions, has a constant value of 0.2 when and is the same as the for neutral atoms when . We also show the second order and the third order polynomial fit to the tabulated data. Although the third order gives a better fit to the data, it drops rapidly at larger values of . Thus, we adopt the second order polynomial fit in our calculations. At high energies, the Bethe approximation shows that the collisional excitation cross-section, like the ionization cross-section, scales as .

Note that the Bethe approximation is only for permitted transitions (electric dipole transitions among states with the same spin). However, impact excitation by non-thermal electrons can occur for all transitions. To make allowance for those forbidden transitions, we use the collision strengths to compute the excitation cross sections, such that


where is the statistical weight of level and is the collision strength. The main difficulty with this approach is that very few collision strengths are available. What we usually have are the thermally averaged effective collision strengths, defined by:


where is the kinetic energy of the outgoing electron. Fortunately, the thermally averaged effective collision strengths for forbidden transitions are almost constant at high temperatures, and hence we adopt such constant effective collision strengths as the collision strengths for impacting non-thermal electrons at all energies. This approximation is unlikely to have a major influence on the results since at high energies collisions for permitted transitions dominate (see Section 5.2), while at low energies there are relatively few non-thermal electrons, and collision rates are dominated by thermal electrons. An important distinction from non-thermal ionization is that non-thermal excitation does not produce secondary electrons.

Figure 1: The tabulated gaunt factor (stars) for neutral atoms as a function of given in van Regemorter (1962). The second order (solid) and the third order (dashed) polynomial fit are also shown.

The accuracies and the uncertainties of the Arnaud Rothenflug cross sections and the Bethe approximation will be discussed in Section 8.1.

3 The hydrodynamical input

The hydrodynamical model ‘lm18a7Ad’ (Woosley, private communication) we use as an initial input is the same as that adopted by Dessart & Hillier (2010). We briefly summarize the main properties of this model below. The hydrodynamical model was produced by the code kepler (Weaver et al., 1978), using a progenitor with main-sequence mass of 18 M. The metallicity of the Large Magellanic Cloud (LMC) was adopted, which is . The star is a BSG when the explosion is initiated.

The hydrodynamical model at 0.27 day was employed, since at this time homology is a good approximation in the outer regions. The rest of the ejecta is enforced to be strictly homologous to accommodate cmfgen, although the innermost part of the ejecta is trimmed since it suffers fallback. Moderate mixing is induced manually to soften the composition gradients. Hydrogen is deficient below 1500 km s, while is primarily found below 2000 km s. The total amount of in the initial hydrodynamical model is 0.084 M. Although a large set of model atoms are utilized, some potentially important species such as Ti i, Fe i, Sc i and Sc ii are missing. The pre-SN steady and explosive nucleosynthesis in the kepler model were undertaken with a 13-isotope network which only captures the composition approximately. It also only computes explicitly the abundance of one unstable isotope () – numerous intermediate-mass elements (IMEs) and iron-group elements (IGEs) are bundled together. Moreover, the hydrodynamical model is enforced to be homologous, which affects the kinematics of the inner ejecta. With the recent availability of a model atom for Fe i we discuss the influence of Fe i on the spectra in Section 6.

The models from day 0.27 to day 20.8 have been presented in Dessart & Hillier (2010), focusing on early evolution at the photospheric phase.

4 The non-thermal model

To study the long-term evolution of our SN model, we evolved the model further in time until day 1000 – well into the nebular phase. After about 40.6 days, the sequence is separated into two branches, a “thermal sequence” with -ray deposition as heat exclusively, and a “non-thermal sequence” with non-thermal excitation and ionization.

In this section, we focus on a typical non-thermal model and discuss the properties of non-thermal processes in the model. We select a non-thermal model at day 127 (hereafter model D127_NT) near the beginning of the nebular phase. It is during the nebular phase that non-thermal effects will be the most important. Moreover, complications arising from dust will not be present since dust formation is not thought to happen until day 300 (Suntzeff & Bouchet, 1990; Gehrz & Ney, 1990).

Figure 2: Illustration of the elemental mass fractions on a logarithmic scale as a function of velocity at day 127. We plot hydrogen (black), helium (red), carbon (green), oxygen (blue) and Co (magenta) – five important elements in the model.

In Fig. 2 we illustrate the elemental mass fractions on a logarithmic scale as a function of velocity at day 127. Hydrogen is an order of magnitude less abundant at velocity 1500 km s and this velocity can be taken to indicate the velocity down to which H has been mixed. Rather than , we show the mass fraction of , since at day 127 only a small amount of is left, and decay is the main energy source. The figure shows that has been mixed out to about 2000 km s in the initial model, although trace amounts also exist at higher velocities.

4.1 The degradation spectrum

A typical degradation spectrum at velocity 1000 km s for the D127_NT model is shown in Fig. 3 (red). The sharp rise at high energy is a result of the adopted source function – we inject all fast electrons at E. The spectrum from E = 1 keV down to 100 eV shows the slowing down of primary electrons, while the rise at low energy indicates a cumulation of secondary electrons. The shape of the degradation spectrum, particularly at low energies, is influenced by the ionization state of the gas.

Uncertainties in the degradation spectrum are introduced by the choice in the number of energy bins, accuracies of the excitation and ionization cross sections, and by the choice of the high energy cut, E. The latter two will be discussed in Section 9. For the number of energy bins, we typically adopt N = 1000 extending from 1 eV to 1 keV. In Fig. 3, we also show another two degradation spectra computed with energy bins N = 100 (black) and N = 10 000 (blue). The flux in the N = 100 degradation spectrum is systematically lower than the other two, except at the high energy edge. The N = 10 000 spectrum is very close to the one with N = 1000, and hence the uncertainty arising from using 1000 energy bins is negligible. We also made spectral comparisons among the above three models. The models with 1000 and 10 000 energy bins showed no distinguishable differences, while the model with 100 energy bins was only slightly different from the other two. Interestingly, we later discovered that a smaller number of bins can be used, provided we scale the degradation spectrum so that the fractional energies entering the three channels sum to unity.

Figure 3: The degradation spectra at a model velocity of 1000 km s computed using different numbers of energy bins: 100 (black), 1000 (red), 10 000 (blue).

4.2 Number density of non-thermal electrons

Knowing the degradation spectrum , the electron spectral density in space can be estimated from , where . The number density of the non-thermal electrons is then given by


Fig. 4 shows the comparison of number density of the non-thermal and thermal electrons at day 127. The number density of the non-thermal electrons is many orders of magnitude smaller than that of the thermal electrons at all depths. The increase of the non-thermal electrons number density at inner regions is consistent with the increase of the Co abundance shown in Fig. 2.

Figure 4: Comparison of the number density of the non-thermal electron (solid) and the thermal electron (dashed).
Figure 5: The fractional energy that goes into heating (red), ionization (blue) and excitation (green) in model D127_NT, as a function of velocity. The ejecta ionization fraction X, which is defined as the ratio of number of ions to number of atoms, is also shown (black).

4.3 Energy fraction of the three channels

With the inclusion of non-thermal excitation and ionization, -rays deposit energy into all three channels. In Fig. 5, we plot the fraction of energy entering the three channels as a function of velocity. We also plot the ejecta ionization fraction X, which shows that the heating fraction is sensitive to X. In the outermost region, the heating fraction reaches unity as X becomes 1 – the highest ionization stage ions are unable to be ionized and excited, quenching the excitation and ionization processes. The high ionization fraction in this region is set by ionization freeze-out due to the time-dependent effects (Dessart & Hillier, 2008). Between 3000 km s and 10 000 km s, non-thermal excitation and ionization become prominent where X becomes very small. As X increases inward from 3000 km s, fractional non-thermal excitation and ionization decline and fractional heating rises and becomes dominant.

4.4 Excitation and ionization

When non-thermal excitation and ionization occur, energy from non-thermal electrons is used to excite or ionize an atom. The Spencer-Fano equation (Eq. 1) implies that the abundances and the cross sections are the key to determine what energy a species “consumes”. In this subsection, we study the fractional energies consumed by a species and compare them to the corresponding species.

Figure 6: The fractional energy (solid lines) consumed by a species in the non-thermal excitation (top) and ionization (bottom) processes as a function of ejecta velocity. Only the largest 5 (summation over all depths) fractions are shown. The fractions of species abundances are plotted in dashed lines for comparison. Importantly, the highest ionization stage of each element is excluded in the calculations of fractional abundances. (see the text for details).

In Fig. 6, we illustrate the fractional energies taken away from the non-thermal electrons by a species in both the excitation and ionization channels. The largest 5 consumers, summing over all depths, are plotted (solid lines). We overplot the fractional abundances of these 5 species (dashed lines). However, the calculations of fractional abundances excludes the highest ionization stages of all elements. These highest ionization stages are the ones in the model, not necessarily the same as the ion of an element with no electron left. For instance, in our models, we utilize detailed model atoms for Fe ii, Fe iii, Fe iv, Fe v, Fe vi and Fe vii, then Fe viii is the highest ionization stage for iron. The reason for excluding the highest ionization stages is that we treat these ionization stages as a single level, and so they don’t have any excitation or ionization routes. In both channels, it is obvious that H i, He i and He ii are the dominant consumers of non-thermal electron energy. He ii is the main consumer above  km s. In the middle region, both H i and He i are prominent. Due to low abundance of H, He i dominates the inner zone, although other species, like C i, C ii, O i and O ii also play a role. Fig. 6 shows that the species abundances, especially H i, He i and He ii, have similar patterns to their fractional energies, in both channels. Such similarity indicates that allocation of the non-thermal energy is generally controlled, as expected, by the species abundances.

5 Non-thermal vs thermal models

In this section, we present a comparison of the temperature and ionization structure between the non-thermal and thermal models. The spectral evolution of the two sequences will also be compared.

5.1 Comparison of the temperature structure

In Fig. 7, we compare the temperature structure between model D127_NT and the thermal model at day 127 (hereafter model D127). As mentioned previously, the non-thermal model generally has a lower temperature in the region between 1000 km s and 5000 km s where the non-thermal effects are crucial. The temperature for the two models above 5000 km s is pretty flat (i.e., almost constant). This is an artifact – a floor temperature of 1500 K is imposed for the models to prevent numerical overflow. The flux mean opacity at 5000 km s is of order , thus the region above 5000 km s has very little influence on spectral formation3.

Figure 7: Comparison of temperature structure between model D127_NT (solid) and model D127 (dashed).

The decrease in the temperature can be understood as follows. The thermal model puts all decay energy into heating, while the non-thermal model only puts a portion of decay energy into heating. Moreover, the coupling between the gas and the radiation at this time is relatively weak, and hence, we expect that only the heating channel will substantially affect the temperature. Energy that has not gone into heat/temperature is stored up in ionization/excitation, driving the material even further away from LTE level populations.

Notice that the temperature of the two models is similar below 1000 km s, despite the existence of a substantial amount of . Even though the region is not in LTE, collisional processes are important, and there is a strong coupling between the radiation field and the gas, and thus energy initially deposited via the ionization and excitation channels can be thermalized. The fractions of energy entering the three channels are shown in Fig. 5. The heating fraction becomes dominant in the inner region, which is consistent with the temperature comparison in Fig. 7.

5.2 Controlling processes for H i and He i lines

In this section, we compare the hydrogen and helium ionization structure in the D127 and D127_NT models and study the processes controlling the H i and He i lines in the non-thermal model. The H and He ionization fractions for both the D127 and D127_NT models are illustrated in Fig. 8. In the cobalt region, the ionized H fraction has increased by up to an order of magnitude, while in the region above km s, the ionization is unchanged, since is deficient there.

Figure 8: Logarithmic species fractions of different ionization stages for hydrogen (top) and helium (bottom). Comparison is made between model D127_NT (solid lines) and D127 (dashed lines).

While He remains neutral for  km s, the singly ionized He fraction is now significant. We can also see significant change in the He ionization fraction above the region where cobalt is abundant, which is due to a small amount of cobalt in the model. At large velocities, the ionization of both H and He is frozen, since the expansion time-scale is smaller than the recombination time-scale (Dessart & Hillier, 2010).

Figure 9: Upper panel: The D127 model spectrum (black) and the model spectra computed by including only bound-bound transitions of He i (red) and H i (blue) at optical band. Lower panel: Same as the upper panel, but for the D127_NT model.
Figure 10: Illustration of model spectrum (solid) and model spectrum computed by including only bound-bound transitions of He i (dashed) for model D127 (top panel) and D127_NT (bottom panel) between 1.9 and 2.2 m.

Fig. 9 shows the model spectra and the spectra with only bound-bound transitions of He i and H i at the optical band. Optical, as well as IR, H i and He i lines mainly originate below 2000 km s in the non-thermal model, with a small fraction of the line flux coming from 2000 - 3000 km s for some lines, e.g., He i 1.083 m. H i lines in the thermal model (upper panel) are generally weak and He i lines are absent. In the non-thermal model (lower panel), all H i lines are strengthened, with H particularly boosted by non-thermal processes. He i lines also contribute significantly to the spectrum. While most optical He i lines are contaminated by other lines, He i 7065 Å might be a good diagnostic to infer the influence of non-thermal processes on helium. Notice that He i 4471 Å is excited in the non-thermal model but is absent from the spectra due to line-blanketing.

Figure 11: Top: Normalized rates for different processes that populate and depopulate the 1s 2p P state of He i. Only processes with a fractional rate 5% at more than one depth are shown. Bottom: Same as the top panel, but for the case of 1s 2p P state of He i.

The prominence of He i 1.083 m and He i 2.058 m due to non-thermal processes at late times (t 200 d) was noted by Li & McCray (1995), and observations showed He i 2.058 m is more isolated (Meikle et al., 1989). In Fig. 10, we show the model spectra and the spectra with only bound-bound transitions of He i for models D127_NT and D127, focusing on the He i 2.058 m line. Similar to the optical band, no He i line is present in model D127 in the wavelength range from 1.9 m to 2.2 m. However, He i 2.058 m produces a prominent absorption feature in the D127_NT model. This feature has been seen in observations of SN 1987A (Meikle et al., 1989), providing strong evidence for the influence of non-thermal processes on He i. Li & McCray (1995) showed that He i 2.058 m is mainly due to non-thermal excitation and He i 1.083 m is mainly excited by thermal electrons during 200 d t 450 d.

Figure 12: Normalized rates for H i (top) and He i (bottom), including all processes allowed in the modeling.

In Fig. 11, we illustrate the major fractional rates ( 5% at any depth) that populate and depopulate the states 1s 2p P and 1s 2p P of He i for model D127_NT. Positive normalized rates illustrate the routes that populate the state, while negative ones show the routes that depopulate the state. States 1s 2p P and 1s 2p P are the upper levels for the He i lines 2.058 m and 1.083 m, respectively. For the state 1s 2p P, we confirm that non-thermal excitation is the main excitation mechanism. Cascades from higher levels, especially the 1s 3d D state, also play an important role at this stage.

For the state 1s 2p P, cascades from higher levels (mainly 1s 3s S and 1s 3d D) are the key processes to populate it below 3000 km s. Recombination is also non-negligible. These are indirectly related to the non-thermal ionization processes. Thus, in contrast to the He i 2.058 m line, non-thermal excitation is largely irrelevant for He i 1.083 m. This is due to the difference in the collision strength behavior of the two transitions. Non-thermal excitation from the ground state (1s 1s S) is generally the dominant excitation route, since the ground state is the most populated state. At low energies, the collision strength for the singlet transition (permitted transition), (1s 1s S - 1s 2p P), is similar to that of the triplet transition (forbidden transition), (1s 1s S - 1s 2p P). However, (1s 1s S - 1s 2p P) grows significantly as the impact energy increases, while (1s 1s S - 1s 2p P) is approximately constant at high energies. As a result, non-thermal excitation processes contribute much less to the triplet states than to the singlet states. Due to radiative excitation from 1s 2s S to 1s 3p P and subsequent decay to 1s 3s S, electron cascade from 1s 3s S, rather than 1s 3d D, is the dominant process populating the 1s 2p P state above 3000 km s.

Figure 13: Comparison of optical and IR spectra between models D127_NT and D127.
Figure 14: Comparison of the H profile between models D127_NT (solid) and D127 (dashed).
Figure 15: Left panel: montage of synthetic optical spectra for the thermal sequence. The days since breakout are labeled below each spectra. Right panel: montage of synthetic spectra at the same epochs shown in the left panel for the non-thermal sequence. The synthetic spectra are reddened with E(B-V) = 0.15 and scaled for a distance of 50 kpc. The most prominent difference between the spectra of the two sequences is the evolution of H – its strength persists at all times in the non-thermal sequence but it almost disappears at late times in the thermal sequence.

Despite the importance of non-thermal processes, we find that photoionization and recombination play a crucial, if not dominant, role in populating many other levels. However, in a “thermal” model, the photoionization and recombination rates are too small to produce a prominent effect. The trigger for the large photoionization and recombination rates is non-thermal processes. In Fig. 12 we show the normalized rates affecting the H i and He i ionization balance for model D127_NT. For H i, the fraction of non-thermal energy deposited via the ionization channel is small and photoionization and recombination are the dominant processes and they balance each other below 8000 km s. We check an individual level, , of H i in both models (D127_NT and D127), and find that the photoionization and recombination rates are increased by a factor of 2 when non-thermal processes are included. For the He i 1s 2p P state, non-thermal ionization and excitation rates in model D127_NT are 5 orders of magnitude greater than any rate in model D127. Below 3000 km s, non-thermal ionization becomes significant, and the radiative recombination rate is now balanced by the photoionization and the non-thermal ionization rates. Although the net non-thermal rate is not overwhelmingly dominant, it is crucial for controlling the ionization balance.

5.3 Comparison of optical and IR spectra

Fig. 13 shows the comparison of the optical and IR spectra between models D127_NT and D127. In the optical range (Fig. 13, top panel), D127_NT generally produces slightly lower fluxes but with two enhanced prominent features — H and He i 7065Å. The Ca ii IR triplet is slightly weakened. In the bottom panel of Fig. 13, most of the lines are strengthened, such as He i 1.083 m, Pa 1.094 m, O i 1.129 m, Pa 1.281 m and Pa 1.875 m. However, Ca ii 1.184 m and 1.195 m are weakened.

Dessart & Hillier (2011) found that in their SNe II-P models the Balmer lines suddenly disappeared at end of the photospheric phase – a result of the vanishing of Balmer-continuum photons. However, observations of Type II SNe show a strong H profile during the nebular phase. They attributed this discrepancy to the exclusion of non-thermal excitation and ionization in the models. We present here the comparison of H at day 127 in Fig. 14. While H is almost absent in model D127, it shows a strong P Cygni H profile in model D127_NT. The emission component is stronger and wider and a P Cygni absorption trough is now present. The width of H is related to the outward mixing of . This may be helpful to determine the mixing in Type II SNe, as hydrogen is abundant in these SNe and H is one of the strongest features at nebular epochs.

5.4 Comparison of the spectral evolution

As mentioned in the previous section, Dessart & Hillier (2011) were unable to reproduce the strong H at the nebular phase for red supergiant (RSG) models, which was associated with the neglect of non-thermal processes. Moreover, the mixing in their models may have been too weak. The best explanation for the persistence of H at late times is non-thermal excitation and ionization processes combined with mixing of H and . This has been one of the motivations for this study. In Fig. 15, we present the spectral evolution of the non-thermal and the thermal sequences. The disappearance of H is also observed in the thermal model for SNe resulting from a blue supergiant (BSG) progenitor. The strength of H increases till day 96, and then fades quickly. On day 140, no prominent feature is seen. The spectral evolution of the non-thermal sequence shows similarity, except for H, which has a stronger maximum strength at about day 96, decreases gradually after that, and maintains a stronger signature on day 140. This confirms that the persistence of strong H emission during the nebular phase is closely related to non-thermal excitation and ionization. The non-thermal effects on spectra at day 54 are relatively small, confirming that the neglect of the non-thermal processes in earlier models is unlikely to have a significant influence on model spectra.

Figure 16: Comparison of optical spectral between model D127_NT and the model with Fe i (hereafter model D127_NT_FeI).

6 The influence of Fe i

Figure 17: Comparison of the observed (black stars, Suntzeff & Bouchet 1990) bolometric light curve of SN 1987A to the theoretical bolometric light curves of the thermal (red) and non-thermal (blue) sequences, with symbols referring to the computed epochs of models. The observational bolometric light curve is constructed by using ultraviolet, optical and infrared photometry. A light curve resulting from the conversion of energy from the radioactive decay of 0.084 M (green) is also shown for reference.

The source of opacity in SNe is complicated, since lots of IMEs and IGEs are synthesized during the explosion. Although a huge number of levels are included in the model, we are still missing some model atoms and the amount of levels of current model atoms may also be insufficient. Many of the missing model ions are the lowest ionization stages. These do not influence the spectra at early times when the ejecta is hot, but may significantly influence the opacity at late times when the ejecta becomes cold. Here, we explore the influence of Fe i. Oscillator strengths are from Kurucz (2009)4, photoionization cross sections are from Bautista (1997) (accesed through Nahar-OSU-Radiative-Atomic-Data web site), while collision rates among the lowest 10 levels are from Pelan & Berrington (1997) (accessed from TIPBASE at

We included Fe i at day 5 and ran a sequence of models with this new atom and with non-thermal excitation and ionization. Fig. 16 illustrates the comparison of the optical spectra at day 127. With the inclusion of Fe i, the emergent flux decreases appreciably between 3000 Å and 5500 Å, with the energy redistributed to longer wavelengths above 6000 Å. Many other atoms, such as Co i, Ti i and Sc ii, may cause effects similar to Fe i. We will include more complete model atoms, as they become available, to better address line blanketing.

7 Comparison with the observations

Figure 18: Comparison of the V-, R-, and I-band light curves of SN 1987A to the corresponding synthetic light curves. The observed V, R, and I magnitudes are shown in stars, and the synthetic V, R, and I magnitudes from the non-thermal sequence are plotted in squares connected by dashed-dotted lines. A scaling is applied to the R and I band magnitudes, both observed and synthetic, to optimize visualization. The synthetic photometry is computed by convolving the synthetic spectra with the Landolt filter bandpasses (Landolt, 1992). We applied a reddening of E(B-V)=0.15 and adopted a distance of 50 kpc.

SN 1987A is still one of the best observed and studied SNe, and it is about 50 kpc away from the earth. Although its progenitor was a BSG, there is still debate about the evolutionary channel that led to a SN explosion as a BSG (Hillebrandt et al., 1987; Langer et al., 1989; Podsiadlowski et al., 1990). In this section, we present photometric and spectroscopic comparison with the observations of SN 1987A. Our observational data is taken from CTIO (Phillips et al., 1988). For the synthetic spectra, we adopt the extinction curve of Cardelli et al. (1989) and use a reddening of E(B-V). Panagia (2005) derived the distance to SN 1987A to be 51.4 kpc with an uncertainty of 1.2 kpc. We adopt  kpc for consistency with Dessart & Hillier (2010).

Figure 19: Comparison between observed spectrum of SN 1987A and the model spectra at day 127. A scaling is applied to the observed spectrum to correct for the difference between the spectroscopic (4.30) and the photometric (4.37) V-magnitude. The three model spectra are scaled by 0.07/0.084 according to the mass in SN 1987A and in the model.

7.1 The synthetic and observed light curves

The comparison of the bolometric light curves is illustrated in Fig. 17, while Fig. 18 illustrates the comparison of the V-, R-, and I-band light curves. The theoretical bolometric light curve rebrightens at a later time than observation, which is probably due to insufficient outward mixing of (Utrobin, 2004; Dessart et al., 2012). The delay of the predicted peak luminosity supports this idea. Both the thermal and non-thermal bolometric light curves predict a higher luminosity than the observed one, partly due to higher mass in our models. Suntzeff & Bouchet (1990) derived a mass of 0.07 M, while we have 0.084 M in the model, which is the amount that comes out of the explosion model. As the purpose of this work is to explore the influence of non-thermal processes, and not to fit the spectra of SN 1987A, we did not adjust the progenitor model, nor did we explore the influence of mixing, explosion energy, and mass cut in the progenitor model. The difference in the amount of only affects our results quantitatively, not qualitatively. The hydrodynamical model, such as the size of the helium core and the enforcement of homology at the beginning of the modeling, also introduces uncertainties. Considering that we have no free parameter, the synthetic bolometric light curves agrees reasonably well with the observed bolometric light curves. We also plot a light curve resulting from the radioactive decay of 0.084 M of . The agreement between this light curve and the synthetic bolometric light curve during the nebular epoch demonstrates that the energy source of the ejecta is totally dominated by radioactive decay.

Although the synthetic bolometric light curve is reasonably well reproduced, the synthetic multi-band light curves show varying degrees of systematic offsets. We are mainly interested in the nebular epochs when the impact of mixing is less important. The R-band magnitudes are close to the observations 130 days after explosion. However, the V- and I-band magnitudes show significant discrepancies (Fig. 18) on the tail of the light curve, with the V-band overestimated by 1 magnitude and the I-band underestimated by 0.4 magnitude. The difference in the mass between the models and SN 1987A translates into 0.2 magnitude in all bands. This does not remove the discrepancy in the V and I bands. However, the overestimation in the V-band and the underestimation in I-band indicate that strong line-blanketing effects are missing in the models. These would redistribute the V-band fluxes to the I-band, possibly alleviating the discrepancy. Future studies with more complete model atoms may quantify such redistribution. In fact, the non-thermal sequence with Fe i reduces the V magnitudes by 0.2 magnitude.

Figure 20: Same as Fig. 19, but it now shows the comparison of the H profile in the velocity space between observation (solid) and model D127_NT (dashed).

7.2 Spectral comparison

Fig. 19 shows a comparison of the observed optical spectra against the model spectra of D127_NT (red), D127 (blue) and D127_NT_FeI (green). The H profile is too weak in model D127, while both non-thermal models, D127_NT and D127_NT_FeI, show considerably stronger H. Models D127_NT and D127 show much stronger fluxes between 3000 Å and 6000 Å, which is reflected by the 1-magnitude offset in the V-band magnitudes in the previous section. The inclusion of Fe i improves the fit, mainly in the wavelength range from 3000 Å to 4500 Å.

In Fig. 20, we show the theoretical and observed H profiles in velocity space. While the peak of H in the non-thermal model is comparable to that of the observation the profile is narrower. The left emission wing of the profile fits the observation, but the H profile in the observation extends at least 3000 km s further to the red than does the theoretical profile. Further, the P Cygni absorption component does not match — the minimum of the model absorption sits at 3000 km s, while observational minimum is at 5000 km s.

Figure 21: Left panel: montage of observed optical spectra of SN 1987A. The days since breakout are shown below each spectra. Right panel: montage of synthetic spectra of the non-thermal sequence. Each observed spectrum on the left is compared to a synthetic spectrum at roughly the same time since explosion. The synthetic spectra are reddened with E(B-V) = 0.15 and scaled for a distance of 50 kpc.

The full width at half maximum (FWHM) of H in the synthetic spectra is about 4000 km s, indicating that the main region that contributes to H is the core zone below 2000 km s. This is coincident with the region in which is abundant. The spatial distribution of depends on the degree of mixing. Mixing effects in SNe explosion are widely seen in 2D and 3D models, however in pure 1D hydrodynamic models mixing must be induced artificially. According to the observed broad H profile, our 1D hydrodynamical model has insufficient mixing of . To constrain the mixing is not the goal of this paper, but the influence of mixing will be discussed in Section 9. Another possible explanation is -ray transport. We made the assumption that all -rays deposit energy locally, but -ray transport makes it possible to deposit more energy at larger velocities. A discussion of the possible influence of -ray transport will be carried out in Section 9.

7.3 Comparison of the spectral evolution

To better compare the spectral features between the observations and the synthetic spectra, we selected the synthetic spectra which are about 10 days later than the observations. In Fig. 21, the observations are shown on the left panel and the synthetic spectra are shown on the right panel. Recall that Balmer lines disappear suddenly at the end of the photospheric phase in the thermal sequence, thus the biggest improvement of the non-thermal sequence is the persistence of the strong H profile at late times. In the non-thermal sequence, the strength of the H profile first increases, and then decreases slowly, which is consistent with the observations. The Na i D lines are contaminated by He i 5875 Å. However, Na i D lines show behavior similar to H, being narrower and weaker than the observations. This is also related to the non-thermal effects, since higher H ionization produces more electrons, which leads to an increase in neutral Na (Dessart & Hillier, 2008). The synthetic spectra are rich with weak lines around 4000 Å, while the observations show few weak lines there. Fe ii is the main contributor to these lines. With the inclusion of Fe i, this discrepancy becomes much smaller. Moreover, the model does not contain scandium and barium, which may explain a few missing P Cygni profiles in the synthetic spectra. The absence of Ba ii 6142 Å is suggestive. The absence of Sc ii in our model likely leads to an underestimate of line blanketing around 4000 Å (Dessart & Hillier, 2011).

8 Uncertainties

8.1 Impact excitation and ionization cross sections

(a) Scale excitation cross sections
(b) Scale ionization cross sections
Figure 22: Comparison of the energy fractions of the three -ray decay channels in models in which we have artificially scaled the excitation or ionization cross sections. The model used for testing is the one at day 127 without including Fe i. Left: The 3 test cases of scaling the impact excitation cross sections. The scaling conditions are denoted in the figure. Different colors are fractions of different channels. The solid lines shows the fractions of the original model. The dotted lines, the dashed lines and the dashed-dotted lines represent the EH5 model, the EM5 model and the EA5 model, respectively. Right: The same as the left, but it shows the 3 cases of scaling the impact ionization cross sections. The dotted lines, the dashed lines and the dashed-dotted lines show the channel fractions of the IH5, IM5 and IA5 model, respectively.

The uncertainty of the non-thermal solver mainly comes from the uncertainties of the electron impact excitation and ionization cross sections. As mentioned previously, the impact excitation cross sections are computed with the Bethe approximation, and the impact ionization cross sections utilize the analytical formula of Arnaud & Rothenflug (1985) with updated coefficients. The Bethe approximation is a high-energy approximation derived from the Born approximation that neglects the short range interaction between the perturbing electron and the atomic electron. van Regemorter (1962) corrected the Gaunt factor according to existing experimental data and accurate calculations to allow the formula to be used at low energies. The Bethe approximation generally works better at high energies. For some less understood species, the Bethe approximation may break down even at high energies. The Arnaud & Rothenflug formula may have troubles for elements heavier than helium. Both experimental and theoretical cross sections are required to evaluate their uncertainties, but both of them are very difficult to estimate (Burke et al., 1994; Ramsbottom et al., 2007; Ramsbottom, 2009).

Model excitation ionization
EH5 H and He by 5 Same
EM5 Metal species by 5 Same
EA5 All species by 5 Same
IH5 Same H and He by 5
IM5 Same Metal species by 5
IA5 Same All species by 5
Table 1: Combinations of species to scale the cross sections

To investigate the influence of the uncertainties, we ran six test models by artificially scaling the impact excitation or ionization cross sections for different combinations of species. These test models were all based on the same input model computed using standard cross-sections, and were run until the populations and the radiation field were fully converged. Table 1 lists the combinations and scalings applied to the cross sections. In Fig. 22, we compare the fractional energies of the three channels in these six tests. In the case of varying the excitation cross sections (Fig. 22, left panel), increasing the excitation cross sections of H and He results in an increase of fractional energy in the excitation channel and a decrease of fractional energy in the ionization channel. The change in the heating fraction is relatively small. An increase in the excitation cross sections of metal species has no effect above 2000 km s, which is obviously due to their abundance deficiency. Below 2000 km s, this also tends to increase the fraction of energy in the excitation channel and to decrease the fraction of energy in the ionization channel. The ejecta ionization fraction X in all cases is roughly the same, implying that the ionization structure of the ejecta is hardly affected. Varying the ionization cross sections (Fig. 22, right panel) has the reverse effect of varying the excitation cross sections – the increase of ionization cross sections results in an increase in the fractional energy of the ionization channel and a decrease in the fractional energy of the excitation energy. This is what we expect from the Spencer-Fano equation, given that the source term S(E) is fixed.

The fractional heating remains almost unchanged between 2500 km s and 10 000 km s when the ionization cross sections increase. This surprising behavior is related to a low ejecta ionization fraction, X, seen in Fig. b. At these low-ionization fractions, the degradation spectra show that almost all electrons have degraded into very low energies and can no longer cause ionizations. Since an enhancement in the ionization cross-section does not dramatically change the number of secondary electrons produced, the ionization channel gains energy at the expense of the excitation channel. Conversely, a change in the excitation cross-sections does change the amount of energy entering the thermal channel. Crudely, excitation causes higher energy losses before the creation of a secondary electron via ionization. The competition between various processes is extremely sensitive to X. As X becomes appreciable (e.g., above 10 000 km s.), many degraded electrons possess high energy. Increasing either the ionization or the excitation cross sections strengthens non-thermal processes at the expense of thermal energy deposition.

There are two types of uncertainties introduced by the approximations in computing the excitation and ionization cross sections. One is the uncertainty in the electron degradation spectrum. The second is the uncertainty in the rate for an individual process. For hydrogen these are explicitly coupled, but this is not necessarily so for other species since such species have a negligible influence on the degradation spectrum. In the six test cases, it turns out that the model spectra show no sizable difference. The H profile shows a slightly stronger dependence on the impact cross sections of hydrogen. However, the two approximations adopted to compute the excitation and ionization cross sections produce fairly good calculations for hydrogen, which reduces our concern about such uncertainties. Generally, a factor of five in the uncertainties of the cross sections will not dramatically influence our results.

8.2 E for high energy electrons

Radioactive decay from and generally produces fast electrons with energy 1 MeV. However, in the non-thermal solver, we assume a source function that injects all electrons with an energy E = 1 keV, which is very small compared to reality. While this high energy cut also introduces another source of uncertainty, the Bethe approximation, the Arnaud Rothenflug formula, and the electron thermal term, have similar asymptotic behaviors at high energy, so the uncertainty should be small as long as E is large enough. We have run two additional models with E = 2 keV and E = 0.5 keV. The difference in the fractional energies entering the three channels between one of the test models and the model with E = 1 keV is less than 5%, and the difference between the two test models is slightly larger. Importantly, the choice of E makes no difference to the spectra, despite the effects on the deposition fractions. Previous work by Xu & McCray (1991) showed that the deposition fractions only have a weak sensitivity to E.

8.3 The time-dependence effects on non-thermal processes

Another uncertainty comes from the late inclusion of non-thermal excitation and ionization. At very early times the region where decays is in LTE, and the assumption that all the energy goes into heating is excellent. However, we included non-thermal excitation and ionization on day 40.6, when the time-dependent effect is already functioning at large velocities, in regions of low ionization and underabundant in Ni/Co. We study this uncertainty in two ways. The first method is to compare the thermal and the non-thermal model at day 40.6. At this epoch, non-thermal excitation and ionization are still unimportant. The model spectrum is hardly affected by the inclusion of non-thermal processes. This is also the main reason why we start the non-thermal sequence at this epoch. The second method is to compare two non-thermal models at a later time – one is evolved from day 40.6 and the other with immediate inclusion of non-thermal processes. At day 127, we found very subtle differences between the spectra of the two models. The H profile is only slightly weaker in the model with immediate inclusion of non-thermal processes. The above two tests indicate that the late inclusion of non-thermal excitation and ionization only introduces very small effects on the nebular spectra. Future studies will start with non-thermal processes, thereby removing such uncertainty.

9 Discussion

The amount of mixing in SN 1987A has been widely debated (See Arnett et al. 1989 for a review). The two fundamental questions are what velocity hydrogen is mixed down to, and what velocity is mixed outward to. The latest 3D simulations of CCSNe (Hammer et al., 2010) showed that significant amounts of hydrogen are transported into deep layers of the ejecta ( 1000 km s), and that a large amount of heavy elements are strongly mixed above 2000 km s. These results are thought to be very similar to the case of SN 1987A. However, it is difficult to draw a conclusion, since the simulations have a different progenitor from SN 1987A. Kozma & Fransson (1998b) modeled the late time spectra of SN 1987A and found hydrogen was mixed down to 700 km s(they needed the mixing to fit the observations).

Fassia & Meikle (1999) modeled the He i 1.083 m line and found some amount of is required to be mixed above 3000 km s to obtain a best fit. Mitchell et al. (2001) synthesized the spectrum of SN 1987A and argued that the strong Balmer lines at early times is due to non-thermal effects, which required a substantial amount of mixing out to 10 000 km s in the hydrogen envelope – velocities much larger than predicted in hydrodynamic simulations. However, Utrobin & Chugai (2005) and Dessart & Hillier (2008) showed time-dependent ionization is the main reason for the strong Balmer lines. The inward mixing of hydrogen and outward mixing of heavy elements in SN 1987A are still open questions.

The scale of mixing also plays an important role. Microscopic mixing results in homogeneous compositions in the ejecta, while macroscopic mixing can lead to large-scale compositional inhomogeneities. With microscopic mixing the mixed species can interact directly (e.g., charge exchange of H with Fe). With macroscopic mixing, this will not occur even if H and Fe are at the same radius (but different , ). Microscopic mixing can be treated in 1D, while an accurate treatment of macroscopic mixing may require 3D models. In SN 1987A, the high velocity feature of iron-group emission lines (e.g., [Fe ii] 17.94 m and 25.99 m (Haas et al., 1990); [Ni ii] 6.6 m (Colgan et al., 1994)), at late times strongly indicate macroscopic mixing. Fransson & Chevalier (1989) showed that microscopic mixing can significantly change the observed spectra at late times, and Li et al. (1993) emphasized the importance of a clumpy structure on the Fe/Co/Ni lines. Jerkstrand et al. (2011) were able to reproduce reasonably the SN 1987A spectra 8 years after explosion with a purely macroscopic mixing model. Furthermore, hydrodynamical models were unable to produce efficient microscopic mixing (Shigeyama & Nomoto, 1990; Mueller et al., 1991; Herant & Woosley, 1994).

Modeling the late time spectra of SN 1987A provides insight into the amount of mixing. The disappearance of Balmer continuum photons at the nebular phase makes non-thermal excitation and ionization the only way to reproduce the Balmer lines. However, the problem is complicated by the -ray transport. -rays interact with the medium through three processes: photoelectric ionization, Compton scattering and pair production. The decay of and mainly produces -ray with energy of 1 MeV. At this energy range, Compton scattering is the dominant process. To give a rough estimate of the mean free path of a high energy electron before it is Compton scattered by atoms, we adopt an effective, purely absorptive, gray opacity  cm g (Swartz et al., 1995), where is the total number of electrons per baryon. The optical depth is calculated by,


Assuming = 0.5 and integrating from velocity 1500 km s, the place where is still abundant, we obtain the mean free path of -ray is  cm, which means the photon can travel from velocity 1500 to 1700 km s in model D127_NT. No doubt, very few -rays travel a long way before they interact with the medium and deposit their energies. A Monte Carlo code has been developed for computing the -ray transport (Hillier & Dessart, 2012), which will give a more detailed understanding of the mixing problem. We ran a test with allowance for -ray transport at day 127 and, as expected, there is only a slight enhancement in the line strength and width of H.

The first discovery of hard X-rays in SN 1987A was made by Sunyaev et al. (1987), approximately 5 months after the explosion, and the first detection of -ray was made by Matz et al. (1988), shortly after the discovery of the hard X-rays. This suggests that -ray transport was beginning to become important at day 127. However, our model at day 154 gives an estimate that -rays can only travel from 1500 to 1800 km s using the effective gray opacity , which indicates an insufficient outward mixing of in the input. In fact, theoretical modeling of the bolometric, X-ray and -ray light curves require substantial outward mixing of into the hydrogen envelope (Kumagai et al., 1988; Arnett & Fu, 1989).

10 Conclusion

We developed a non-thermal solver, which takes into account ionization and excitation due to non-thermal electrons created by -rays that arise from the decay of and its daughter products, and incorporated it into our fully non-LTE time-dependent code to model SNe. We use the model ‘lm18a7Ad’ of Woosley on 0.27 day, enforced into homology as an initial input, and benchmark the non-thermal solver by comparing model results with the observations of SN 1987A.

Non-thermal models have lower temperature and more excited/ionized material in the region where the non-thermal processes are crucial. Fractional energy of non-thermal excitation/ionization is prominent at places with very small ejecta ionization fraction X. The spectral comparison with the thermal models shows that H is strongly increased at nebular epochs. At late times, He i lines, which are totally absent in the thermal models, are present in non-thermal models. While most He i lines are significantly contaminated by other lines, He i 2.058 m provides an excellent opportunity to infer the influence of non-thermal processes on helium. He i 7065Å is a possible optical line that can be used to infer the influence of non-thermal effects in the nebular epoch. There are many O i lines locating near 1.129 m, and they are also strengthened by non-thermal processes. Most of other lines are only weakly affected at the epochs considered here.

Optical and IR He i lines mainly originate below 2000 km s. We confirm that non-thermal excitation is the most important process for He i 2.058 m. However, He i 1.083 m is due to cascades from higher levels, which indirectly relate to non-thermal ionization. Although photoionization and recombination are prominent processes for populating many levels, non-thermal excitation and ionization are the processes controlling the ionization balance.

We also compare the non-thermal models with the observations. The re-brightening of the bolometric light curve peaks about 10 days later in our models. An underestimate in the amount of outward mixing of is the main reason for the delay. Although the synthetic bolometric light curve agrees reasonably with the observations, multi-band light curves show discrepancies in various degrees and reveal potential differences between the hydrodynamical model and SN 1987A or possible problems in our modeling. The opacity issue, related to limited levels in the atomic data or some missing species, is one possibility, considering that the V-band is overestimated and the I-band is underestimated. H maintains a strong profile at the nebular epochs in the non-thermal models, resembling observations of SN 1987A. However, the model H profile is much narrower than what we observed in SN 1987A, probably due to insufficient outward mixing of . Although hydrogen is abundant above 2000 km s in the model, the lack of in this region and the assumption of local deposition of -rays make the non-thermal effects negligible there. We also compared the optical spectral evolution to the observations and found a satisfactory agreement.

The uncertainties introduced by various approximations and assumptions in the non-thermal solver have been investigated. The primary uncertainties are due to the not-well-known electron impact excitation and ionization cross sections of metal species. However, our results for a BSG model are largely dependent on accuracies of the cross sections for hydrogen and helium, since they are the most abundant two elements. The Bethe approximation and the Arnaud & Rothenflug formula produce very good estimates of the excitation and ionization cross sections for these two elements. Therefore, our conclusions are not affected by these uncertainties. Another source of uncertainty comes from the computation of the degradation spectrum, which arises from the choice of the number of energy bins and the energy cutoff E at the high end of the degradation spectrum. We demonstrate that our choice for the number of energy bins, N = 1000, is sufficient to give a reliable degradation spectrum and using E = 1000 eV in our model produces almost the same spectra as using E = 2000 eV. We also explore the time dependence of non-thermal processes and find there is only a very subtle effect on predicted spectra.

With the non-thermal solver, we are able to simulate SNe from photospheric to nebular phases continuously. The influence of non-thermal ionization and excitation on Type Ib and Type Ic models is being investigated (Dessart et al., 2012). As more and more nebular spectra are available, the comparison of observations with models will allow us to place constraints on the hydrodynamic models and nucleosynthesis. The non-thermal solver also provides opportunities to constrain mixing effects in SNe.


CL acknowledges support from NASA theory grant NNX10AC80G. DJH acknowledges support from STScI theory grant HST-AR-11756.01.A and NASA theory grant NNX10AC80G. LD acknowledges financial support from the European Community through an International Re-integration Grant, under grant number PIRG04-GA-2008-239184. We would also like to thank Bob Kurucz, the people involved with the NIST Atomic Data Base, and the participants of the Opacity and Iron Projects (with special thanks to Keith Butler, Sultana Nahar, Anil Pradhan and Peter Storey) for computing and supplying atomic data to the astrophysical community without their work these calculations would not be possible.


  1. pagerange: Non-thermal excitation and ionization in supernovaeReferences
  2. pubyear: 2012
  3. This is confirmed by computing two models with different floor temperatures. The synthetic spectra show no distinguishable difference between models with floor temperatures of 1500 K and 4000 K.
  4. Data from R. L. Kurucz is available online at


  1. Arnaud M., Rothenflug R., 1985, A&APS, 60, 425
  2. Arnett W. D., Bahcall J. N., Kirshner R. P., Woosley S. E., 1989, ARA&A, 27, 629
  3. Arnett W. D., Fu A., 1989, ApJ, 340, 396
  4. Bautista M. A., 1997, A&APS, 122, 167
  5. Burke P. G., Burke V. M., Dunseath K. M., 1994, Journal of Physics B Atomic Molecular Physics, 27, 5341
  6. Cardelli J. A., Clayton G. C., Mathis J. S., 1989, ApJ, 345, 245
  7. Chugai N. N., 1991, Soviet Astronomy Letters, 17, 400
  8. Colgan S. W. J., Haas M. R., Erickson E. F., Lord S. D., Hollenbach D. J., 1994, ApJ, 427, 874
  9. Dessart L. et al., 2008, ApJ, 675, 644
  10. Dessart L., Hillier D. J., 2005, A&A, 437, 667
  11. Dessart L., Hillier D. J., 2006, A&A, 447, 691
  12. Dessart L., Hillier D. J., 2008, MNRAS, 383, 57
  13. Dessart L., Hillier D. J., 2010, MNRAS, 405, 2141
  14. Dessart L., Hillier D. J., 2011, MNRAS, 410, 1739
  15. Dessart L., Hillier D. J., Li C., Woosley S., 2012, MNRAS, (submitted)
  16. Dessart L., Hillier D. J., Livne E., Yoon S.-C., Woosley S., Waldman R., Langer N., 2011, MNRAS, 414, 2985
  17. Fassia A., Meikle W. P. S., 1999, MNRAS, 302, 314
  18. Fransson C., Chevalier R. A., 1989, ApJ, 343, 323
  19. Gehrz R. D., Ney E. P., 1990, Proceedings of the National Academy of Science, 87, 4354
  20. Haas M. R., Erickson E. F., Lord S. D., Hollenbach D. J., Colgan S. W. J., Burton M. G., 1990, ApJ, 360, 257
  21. Hammer N. J., Janka H.-T., Müller E., 2010, ApJ, 714, 1371
  22. Harkness R. P. et al., 1987, ApJ, 317, 355
  23. Herant M., Woosley S. E., 1994, ApJ, 425, 814
  24. Hillebrandt W., Hoeflich P., Weiss A., Truran J. W., 1987, Nature, 327, 597
  25. Hillier D. J., 1987, ApJS, 63, 947
  26. Hillier D. J., Dessart L., 2012, MNRAS, (in press)
  27. Hillier D. J., Miller D. L., 1998, ApJ, 496, 407
  28. Jerkstrand A., Fransson C., Kozma C., 2011, A&A, 530, A45
  29. Kozma C., Fransson C., 1992, ApJ, 390, 602
  30. Kozma C., Fransson C., 1998a, ApJ, 496, 946
  31. Kozma C., Fransson C., 1998b, ApJ, 497, 431
  32. Kumagai S., Shigeyama T., Nomoto K., Itoh M., Nishimura J., 1988, A&A, 197, L7
  33. Kurucz R. L., 2009, in American Institute of Physics Conference Series, Vol. 1171, Recent Directions in Astrophysical Quantitative Spectroscopy and Radiation Hydrodynamics, Hubeny I., Stone J. M., MacGregor K., Werner K., eds., pp. 43–51
  34. Landolt A. U., 1992, AJ, 104, 340
  35. Langer N., El Eid M. F., Baraffe I., 1989, A&A, 224, L17
  36. Li H., McCray R., 1995, ApJ, 441, 821
  37. Li H., McCray R., Sunyaev R. A., 1993, ApJ, 419, 824
  38. Lucy L. B., 1991, ApJ, 383, 308
  39. Matz S. M., Share G. H., Leising M. D., Chupp E. L., Vestrand W. T., 1988, Nature, 331, 416
  40. Meikle W. P. S., Spyromilio J., Varani G.-F., Allen D. A., 1989, MNRAS, 238, 193
  41. Mitchell R. C., Baron E., Branch D., Lundqvist P., Blinnikov S., Hauschildt P. H., Pun C. S. J., 2001, ApJ, 556, 979
  42. Mueller E., Fryxell B., Arnett D., 1991, A&A, 251, 505
  43. Opal C. B., Peterson W. K., Beaty E. C., 1971, JChPh, 55, 4100
  44. Panagia N., 2005, in IAU Colloq. 192: Cosmic Explosions, On the 10th Anniversary of SN1993J, J.-M. Marcaide & K. W. Weiler, ed., pp. 585–+
  45. Pelan J., Berrington K. A., 1997, A&APS, 122, 177
  46. Phillips M. M., Heathcote S. R., 1989, PASP, 101, 137
  47. Phillips M. M., Heathcote S. R., Hamuy M., Navarrete M., 1988, AJ, 95, 1087
  48. Podsiadlowski P., Joss P. C., Rappaport S., 1990, A&A, 227, L9
  49. Ramsbottom C. A., 2009, Atomic Data and Nuclear Data Tables, 95, 910
  50. Ramsbottom C. A., Hudson C. E., Norrington P. H., Scott M. P., 2007, A&A, 475, 765
  51. Shigeyama T., Nomoto K., 1990, ApJ, 360, 242
  52. Shull J. M., van Steenberg M. E., 1985, ApJ, 298, 268
  53. Spencer L. V., Fano U., 1954, Physical Review, 93, 1172
  54. Suntzeff N. B., Bouchet P., 1990, AJ, 99, 650
  55. Sunyaev R. et al., 1987, Nature, 330, 227
  56. Swartz D. A., 1991, ApJ, 373, 604
  57. Swartz D. A., Filippenko A. V., Nomoto K., Wheeler J. C., 1993, ApJ, 411, 313
  58. Swartz D. A., Sutherland P. G., Harkness R. P., 1995, ApJ, 446, 766
  59. Utrobin V. P., 2004, Astronomy Letters, 30, 293
  60. Utrobin V. P., Chugai N. N., 2005, A&A, 441, 271
  61. van Regemorter H., 1962, ApJ, 136, 906
  62. Weaver T. A., Zimmerman G. B., Woosley S. E., 1978, ApJ, 225, 1021
  63. Xu Y., McCray R., 1991, ApJ, 375, 190
  64. Younger S. M., 1981, JQSRT, 26, 329
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description