# Non-thermal excitation and ionization in supernovae

Chengdong Li, D. John Hillier, Luc Dessart

► To cite this version:

Chengdong Li, D. John Hillier, Luc Dessart. Non-thermal excitation and ionization in supernovae. Monthly Notices of the Royal Astronomical Society, 2012, 426 (2), pp.1671–1686. 10.1111/j.1365-2966.2012.21198.x. hal-01442156

HAL Id: hal-01442156

https://hal.science/hal-01442156

Submitted on 22 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Non-thermal excitation and ionization in supernovae

Chengdong Li,<fnref n="1" />* D. John Hillier<fnref n="1" /> and Luc Dessart<fnref n="2" />

Accepted 2012 April 26. Received 2012 April 25; in original form 2012 March 12

## ABSTRACT

We incorporate non-thermal excitation and ionization processes arising from non-thermal electrons that result from $\gamma$-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. We applied the non-thermal calculations to the blue supergiant explosion model whose early evolution was studied in Dessart & Hillier. Non-thermal processes generally increase excitation and ionization and decrease the temperature of the ejecta. We confirm that non-thermal processes are crucial for modelling 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 $\mu$m provides an excellent opportunity to infer the influence of non-thermal processes. We also discuss the processes controlling the formation of He\,{\i} lines during nebular epochs. Most lines of other species are only slightly affected. We also show that the inclusion of Fe\,{\i} has substantial line-blanketing effects on 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\alpha$. 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 supernova ejecta, whatever the type and the epoch.

**Key words:** radiation mechanisms: non-thermal – stars: atmospheres – stars: evolution – supernovae: general supernovae: individual: SN 1987A.

## 1 INTRODUCTION

Supernova (SN) 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 three broad classes: Type Ib, Type Ic and Type II. Type II SNe exhibit hydrogen lines, while Type I SNe do not.

SNe explosively produce $^{56}$Ni which decays into $^{56}$Co with a half-life of 6.1 d. $^{56}$Co further decays into $^{56}$Fe with a half-life of 77 d. 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 shock-deposited energy, but at later times radioactive decay will be the dominant energy source. High-energy $\gamma$-rays are produced during radioactive decay. The $\gamma$-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 SNe of different types. For example, the Bochum event in the $H\alpha$ profiles of SN 1987A was explained by non-thermal excitation due to asymmetric ejection of $^{56}$Ni (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 $\mu$m and 2.058 $\mu$m lines whose upper levels are populated by non-thermal processes induced by $\gamma$-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 local thermodynamic equilibrium (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

<footnote n="1">¹PITTSburgh Particle physics, Astrophysics, and Cosmology Center (PITT PACC), Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh, PA 15260, USA</footnote>

<footnote n="2">²Laboratoire d'Astrophysique de Marseille, Université Aix-Marseille & CNRS, UMR7326, 38 rue Frédéric Joliot-Curie, 13388 Marseille, France</footnote>

<footnote>*E-mail: CHL85@pitt.edu</footnote>
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 that result from quasi-hydrogen-less low-mass cores from 15–25 M☉ main-sequence stars in binary systems using CMFGEN. They found that 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 ultraviolet ranges of SN 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, a radiative transfer code developed for modelling Wolf–Rayet stars (Hillier 1987), which 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 rate equations are usually neglected in SN spectral modelling 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 $^{56}$Ni is abundant. However, non-thermal effects will eventually take over other processes, especially during the nebular phase. Further, the study of Jerkstrand, Fransson & Kozma (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 day 10 (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. This 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. A comparison with 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 SN 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.

### 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 γ-rays and thermal electrons lose energy through three different channels – heating, excitation and ionization. The fraction 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

$$ \sum_i n_i y(E) \left[ \sum_j \sigma_{ij}(E) + \int_I^{(I+E)/2} \sigma_c(E, \epsilon) d\epsilon \right] \\ - \frac{d}{dE} [y(E)L_c(E)] = \sum_i n_i \sum_j y(E + E_{ij}) \sigma_{ij}(E + E_{ij}) \\ + \sum_i n_i \int_I^{\lambda} y(E + \epsilon) \sigma_c(E + \epsilon, \epsilon) d\epsilon \\ + \sum_i n_i \int_{2E+I}^{E_{\max}} y(E') \sigma_c(E', E+I) dE' + S(E), \quad (1) $$

where $\lambda = \min(E_{\max} - E, E + I)$, $E_{\max}$ is the maximum energy of non-thermal electrons, and $I$ is the ionization potential. $y(E)dE$ is defined as the flux of non-thermal electrons in the energy interval $(E, E + dE)$. $n_i$ is the number density of a single species, and the sum over $i$ allows for the contribution of different ionization stages and species. $L_c(E)$ is the energy loss of the non-thermal electrons due to Coulomb scattering by the thermal electrons. $E_j$ is the excitation energy, $\sigma_j$ is the excitation cross-section of level $j$ and $\sigma_c(E, \epsilon)$ is the differential ionization cross-section for a fast electron with energy $E$ and an energy loss of $\epsilon$. Therefore, the ejected electron has energy $\epsilon - I$. 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 same three channels as the primary electron. $S(E)$ is the source term, and we assume all non-thermal electrons are injected at $E_{max}$. The degradation spectrum $y(E)$ is normalized to $E_{max}$, so that $y(E)$ is the degradation spectrum for 1 eV energy input. The first and second terms on the left-hand side represent the number of electrons leaving the energy interval $(E, E + dE)$ by impact excitation and ionization, respectively. The term with $L_e(E)$ is the number of electrons leaving the energy interval by heating the medium. The first term on the right-hand side is the number of electrons entering the energy interval $(E, E + dE)$ 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)

$$ 
\begin{aligned}
& \sum_i n_i \sum_j \int_E^{E+E_j} y(E') \sigma_j(E') dE' + y(E)L_e(E) \\
& + \sum_i n_i \int_E^{E_{\max}} y(E') \int_{E'-E}^{(E'+E)/2} \sigma_c(E', \epsilon) d\epsilon \, dE' \\
& = \sum_i n_i \int_{2E+I}^{E_{\max}} y(E') \int_{E+I}^{(E'+I)/2} \sigma_c(E', \epsilon) d\epsilon \, dE' \\
& + \int_E^{E_{\max}} S(E') dE'.
\end{aligned} 
\quad (2) $$

In equation (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 (equation 2), and can be computed analytically.

With the solution of the degradation spectrum, the fractional energies entering heating, $\eta_h$, excitation of ion $i$ to level $j$, $\eta_{ij}$, and ionization of ion $i$, $\eta_{ic}$, can be computed by

$$ \eta_h = \frac{1}{E_{\text{init}}} \int_{E_0}^{E_{\text{max}}} y(E') L_e(E') dE' + \frac{1}{E_{\text{init}}} E_0 y(E_0) L_e(E_0) + \frac{1}{E_{\text{init}}} \int_0^{E_0} N(E') E' dE', \quad (3) $$

$$ \eta_{ij} = \frac{n_{ij} E_{ij}}{E_{\text{init}}} \int_{E_{ij}}^{E_{\max}} y(E') \sigma_{ij}(E') dE', \quad (4) $$

$$ \eta_{ic} = \frac{n_i I_i}{E_{\text{init}}} \int_{I_i}^{E_{\max}} y(E') \sigma_{ic}(E') dE', \quad (5) $$

where

$$ N(E) = \sum_i n_i \left[ \begin{aligned} & \sum_j y(E + E_{ij}) \sigma_{ij}(E + E_{ij}) \\ & + \int_{I_i}^{\lambda_i} y(E + E') \sigma_{ic}(E + E', E') dE' \\ & + \int_{2E + I_i}^{E_{\max}} y(E') \sigma_{ic}(E', E + I_i) dE' \end{aligned} \right] + S(E) \quad (6) $$

(Kozma & Fransson 1992). $E_{init}$ is the mean energy of the initial electrons, and $E_0$ is the lowest excitation or ionization energy.

To save computational effort, we take into account only the dominant ionization stages when computing the degradation spectrum. The other approximation is the neglect of $\gamma$-ray transport. In reality, some $\gamma$-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 $\gamma$-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 $E_{max}$. 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-shaped source function near $E_{max}$. This tends to give a smoother degradation spectrum near $E_{max}$.

## 2.2 Ionization cross-sections

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

$$ Q_i(E) = \frac{1}{u I_i^2} \left[ A_i \left(1 - \frac{1}{u}\right) + B_i \left(1 - \frac{1}{u}\right)^2 + C_i \ln(u) + D_i \frac{\ln(u)}{u} \right], \quad (7) $$

where $u = E/I_i$. $E$ is the energy of the impact electron, and $I_i$ is the ionization potential for the element. $A_i$, $B_i$, $C_i$ and $D_i$ 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 equation (7), the differential cross-section $\sigma_c(E, \epsilon)$ can be written as

$$ \sigma_c(E, \epsilon) = Q(E)P(E, \epsilon - I), \quad (8) $$

where $P(E_p, E_s)$ is the distribution of the secondary electrons for a primary electron energy of $E_p$. Opal, Peterson & Beaty (1971) found $P(E_p, E_s)$ could be well described by

$$ P(E_p, E_s) = \frac{1}{J \arctan\left[\frac{(E_p - I)/2J}{1 + (E_s/J)^2}\right]}, \quad (9) $$

where $J$ is a parameter that gives a best fit to measurements and varies for different species. We adopt $J = 0.6I$, 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 energy. Specifically, for $H_1$ with $I = 13.6$ eV, 66.7 per cent of the secondary electrons are below the ionization potential if $E_p = 1000$ eV and 69.8 per cent in the case of $E_p = 10000$ eV. Further, the mean energy of the secondary electrons, $\bar{v}_{sec}$, only varies slowly as the energy of the primary electron $E_p$ increases. For hydrogen, $\bar{v}_{sec} = 21.53$ eV for $E_p = 1000$ eV and $\bar{v}_{sec} = 33.36$ eV for $E_p = 10000$ eV. Secondary electrons with low energy 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, $Q_{ij}$ is given by

$$ Q_{ij} = \frac{8\pi}{\sqrt{3}} \frac{1}{k_i^2} \frac{I_H}{E_j - E_i} f_{ij} \bar{g} \pi a_0^2, \quad (10) $$

where $k_i^2$ is the initial electron kinetic energy normalized to hydrogen ionization energy of 13.6 eV, $I_H$ is the hydrogen ionization
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The tabulated Gaunt factor (stars) for neutral atoms as a function of $x = \sqrt{mv^2/2\Delta E}$ given in van Regemorter (1962). The second-order (solid curve) and third-order (dashed curve) polynomial fits are also shown.</figcaption></figure>

energy, $E_j - E_i$ is the energy difference between the two levels, $f_{ij}$ is the oscillator strength, and $a_0$ is the Bohr radius. The Gaunt factor $\bar{g}$ for neutral atoms utilizes table 1 of van Regemorter (1962), which is shown in Fig. 1. The variable $x$ is defined by $x = \sqrt{mv^2/2\Delta E}$, where $\frac{1}{2}mv^2$ is the kinetic energy of the impact electron and $\Delta E$ is the transition energy. For positive ions, $\bar{g}$ has a constant value of 0.2 when $x \le 1$ and is the same as the $\bar{g}$ for neutral atoms when $x > 1$. We also show the second- and third-order polynomial fits to the tabulated data. Although the third order gives a better fit to the data, it drops rapidly at larger values of $x$. Thus, we adopt the second-order polynomial fit in our calculations. At high energy, the Bethe approximation shows that the collisional excitation cross-section, like the ionization cross-section, scales as $\log E/E$.

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

$$Q_{ij} = \frac{1}{k_i^2} \frac{1}{w_i} \Omega_{ij} \pi a_0^2 \quad (11)$$

where $w_i$ is the statistical weight of level $i$ and $\Omega_{ij}$ 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

$$\gamma_{ij} = \int_{0}^{\infty} \Omega_{ij} \exp \left( -\frac{E_j}{kT} \right) d \left( \frac{E_j}{kT} \right), \quad (12)$$

where $E_j$ 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.

The accuracies and 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 '1m18a7Ad' (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, Zimmerman & Woosley 1978), using a progenitor with main-sequence mass of $18 M_{\odot}$. The metallicity of the Large Magellanic Cloud was adopted, which is $Z = 0.4 Z_{\odot}$. The star is a blue supergiant (BSG) when the explosion is initiated.

The hydrodynamical model at day 0.27 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 \text{ km s}^{-1}$, while $^{56}\text{Ni}$ is primarily found below $2000 \text{ km s}^{-1}$. The total amount of $^{56}\text{Ni}$ in the initial hydrodynamical model is $0.084 M_{\odot}$. Although a large set of model atoms are utilized, some potentially important species such as TiI, FeI, ScI and ScII are missing. The pre-SN steady and explosive nucleosynthesis in the KEPLER model was undertaken with a 13-isotope network which only captures the composition approximately. It also only computes explicitly the abundance of one unstable isotope ($^{56}\text{Ni}$) – 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 FeI, we discuss the influence of FeI 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 d, the sequence is separated into two branches: a 'thermal sequence' with $\gamma$-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 (Gehrz & Ney 1990; Suntzeff & Bouchet 1990).

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 \text{ km s}^{-1}$ and this velocity can be taken to indicate the velocity down to which H has been mixed. Rather than $^{56}\text{Ni}$, we show the mass fraction of $^{56}\text{Co}$, since at day 127 only a small amount of $^{56}\text{Ni}$ is left, and $^{56}\text{Co}$ decay is the main energy source. The figure shows that $^{56}\text{Co}$ has been mixed out to about $2000 \text{ km s}^{-1}$ 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$^{-1}$ for model D127_NT is shown in Fig. 3 (red curve). The sharp rise at high
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Illustration of the elemental mass fractions on a logarithmic scale as a function of velocity at day 127. We plot five important elements – hydrogen (black curve), helium (red curve), carbon (green curve), oxygen (blue curve) and $^{56}$Co (magenta curve) – in the model.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** The degradation spectra at a model velocity of ~1000 km s⁻¹ computed using different numbers of energy bins: 100 (black curve), 1000 (red curve) and 10 000 (blue curve).</figcaption></figure>

energy is a result of the adopted source function – we inject all fast electrons at $E_{max}$. The spectrum from $E_{max} = 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_{max}$. 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 the number of energy bins $N = 100$ (black curve) and $N = 10 000$ (blue curve). 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.

## 4.2 Number density of non-thermal electrons

Knowing the degradation spectrum $y(E)$, the electron spectral density in space can be estimated from $y(E)/v$, where $v = \sqrt{2E/m_e}$.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Comparison of the number densities of the non-thermal (solid curve) and thermal (dashed curve) electrons.</figcaption></figure>

The number density of the non-thermal electrons is then given by

$$N_e = \int_0^{E_{max}} \frac{y(E)}{v} dE = \int_0^{E_{max}} y(E) \sqrt{\frac{m_e}{2E}} dE. \quad (13)$$

Fig. 4 shows a comparison of the number densities 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 number density of the non-thermal electrons at inner regions is consistent with the increase of $^{56}$Co abundance shown in Fig. 2.

## 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_e$, which shows that the heating fraction is sensitive to $X_e$. In the outermost region, the heating fraction reaches unity as $X_e$ 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 and 10 000 km s⁻¹, non-thermal excitation and ionization become prominent where $X_e$ becomes very small. As $X_e$ increases inwards from 3000 km s⁻¹, fractional non-thermal excitation and ionization decline and fractional heating rises and becomes dominant.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** The fractional energy that goes into heating $\delta E_h$ (red curve), ionization $\delta E_i$ (blue curve) and excitation $\delta E_e$ (green curve) in model D127_NT, as a function of velocity. The ejecta ionization fraction $X_e$, which is defined as the ratio of the number of ions to the number of atoms, is also shown (black curve).</figcaption></figure>
## 5.2 Controlling processes for H I and He I lines

In this section, we compare the hydrogen and helium ionization structures in models D127 and D127_NT 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 model D127 and model D127_NT 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 $v > 5000$ km s⁻¹, the ionization is unchanged, since ⁵⁶Co is deficient there.

While He remains neutral for $v < 10000$ 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).

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 infrared (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, for example, He I 1.083 µm. H I lines in the thermal model (top panel) are generally weak and He I lines are absent. In the non-thermal model (bottom panel), all H I lines are strengthened, with $H\alpha$ 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><img src="image_8.png" /><figcaption>**Figure 8.** Logarithmic species fractions of different ionization stages for hydrogen (top panel) and helium (bottom panel). A comparison is made between models D127_NT (solid curves) and D127 (dashed curves).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 9.** Top panel: model D127 spectrum (black curve) and the model spectra computed by including only bound-bound transitions of He I (red curve) and H I (blue curve) at the optical band. Bottom panel: same as the top panel, but for model D127_NT.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 10.** Illustration of the model spectrum (solid curve) and the model spectrum computed by including only bound-bound transitions of He I (dashed curve) for models D127 (top panel) and D127_NT (bottom panel) between 1.9 and 2.2 µm.</figcaption></figure>

The prominence of He I 1.083 µm and He I 2.058 µm due to non-thermal processes at late times ($t \ge 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 He I 2.058 µm. Similar to the optical band, no He I line is present in model D127 in the wavelength range 1.9–2.2 µm. However, He I 2.058 µm produces a prominent absorption feature in model 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 \le t \le 450$ d.
<figure><img src="image_11.png" /><figcaption>**Figure 11.** Top panel: normalized rates for different processes that populate and depopulate the 1s 2p<sup>1</sup>P<sup>o</sup> state of HeI. Only processes with a fractional rate ≥5 per cent at more than one depth are shown. Bottom panel: same as the top panel, but for the case of the 1s 2p<sup>3</sup>P<sup>o</sup> state of HeI.</figcaption></figure>

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

For the state 1s 2p<sup>3</sup>P<sup>o</sup>, cascades from higher levels (mainly 1s 3s<sup>3</sup>S and 1s 3d<sup>3</sup>D) are the key processes to populate it below 3000 km s<sup>-1</sup>. Recombination is also non-negligible. These are indirectly related to the non-thermal ionization processes. Thus, in contrast to HeI 2.058 µm, non-thermal excitation is largely irrelevant for HeI 1.083 µm. This is due to the difference in the collision strength behaviour of the two transitions. Non-thermal excitation from the ground state (1s 1s<sup>1</sup>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<sup>1</sup>S, 1s 2p<sup>1</sup>P<sup>o</sup>), is similar to that of the triplet transition (forbidden transition), Ω(1s 1s<sup>1</sup>S, 1s 2p<sup>3</sup>P<sup>o</sup>). However, Ω(1s 1s<sup>1</sup>S, 1s 2p<sup>1</sup>P<sup>o</sup>) grows significantly as the impact energy increases, while Ω(1s 1s<sup>1</sup>S, 1s 2p<sup>3</sup>P<sup>o</sup>) 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<sup>3</sup>S to 1s 3p<sup>3</sup>P<sup>o</sup> and subsequent decay to 1s 3s<sup>3</sup>S, an electron cascade from 1s 3s<sup>3</sup>S, rather than 1s 3d<sup>3</sup>D, is the dominant process populating the 1s 2p<sup>3</sup>P<sup>o</sup> state above 3000 km s<sup>-1</sup>.

<figure><img src="image_12.png" /><figcaption>**Figure 12.** Normalized rates for H<sub>I</sub> (top panel) and He<sub>I</sub> (bottom panel), including all processes allowed in the modelling.</figcaption></figure>

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<sub>I</sub> and He<sub>I</sub> ionization balance for model D127_NT. For H<sub>I</sub>, 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<sup>-1</sup>. We check an individual level, n = 1, of H<sub>I</sub> 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 HeI 1s 2p<sup>1</sup>P<sup>o</sup> state, non-thermal ionization and excitation rates in model D127_NT are five orders of magnitude greater than any rate in model D127. Below 3000 km s<sup>-1</sup>, 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 a comparison of the optical and IR spectra between models D127_NT and D127. In the optical range (Fig. 13, top panel), model D127_NT generally produces slightly lower fluxes but with two enhanced prominent features – Hα and HeI 7065 Å. The CaII IR triplet is slightly weakened. In the bottom panel of Fig. 13, most of the lines are strengthened, such as HeI 1.083 µm, Paγ 1.094 µm, OI 1.129 µm, Paβ 1.281 µm and Paα 1.875 µm. However, CaII 1.184 µm and 1.195 µm are weakened.

Dessart & Hillier (2011) found that in their Type II-P SN models the Balmer lines suddenly disappeared at the end of the photospheric
<figure><img src="image_13.png" /><figcaption>**Figure 15.** Left-hand panel: montage of synthetic optical spectra for the thermal sequence. The days since breakout are labelled below each spectrum. Right-hand panel: montage of synthetic spectra at the same epochs as shown in the left-hand 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\alpha$ – its strength persists at all times in the non-thermal sequence but it almost disappears at late times in the thermal sequence.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 16.** Comparison of the optical spectra between models D127 NT and D127 NT FeI.</figcaption></figure>

Both the thermal and non-thermal bolometric light curves predict a higher luminosity than the observed one, partly due to higher $^{56}$Ni mass in our models. Suntzeff & Bouchet (1990) derived a $^{56}$Ni mass of 0.07 $M_{\odot}$, while we have 0.084 $M_{\odot}$ 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 $^{56}$Ni 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 modelling, also introduces uncertainties. Considering that we have no free parameter, the synthetic bolometric light curves agree reasonably well with the observed bolometric light curves. We also plot a light curve resulting from the radioactive decay of 0.084 $M_{\odot}$ of $^{56}$Ni. The agreement between this light curve and the synthetic bolometric

<figure><img src="image_15.png" /><figcaption>**Figure 17.** Comparison of the observed (black stars, Suntzeff & Bouchet 1990) bolometric light curves of SN 1987A to the theoretical bolometric light curves of the thermal (red curve) and non-thermal (blue curve) sequences, with symbols referring to the computed epochs of models. The observational bolometric light curve is constructed by using ultraviolet, optical and IR photometry. A light curve resulting from the conversion of energy from the radioactive decay of 0.084 $M_{\odot}$ (green curve) is also shown for reference.</figcaption></figure>

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 multiband 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
<figure><img src="image_16.png" /><figcaption>**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 the stars, and the synthetic V, R and I magnitudes from the non-thermal sequence are plotted in the squares connected by the dot-dashed curves. 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.</figcaption></figure>

magnitudes are close to the observations ~130 d 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 mag and the I band underestimated by ~0.4 mag. The difference in $^{56}$Ni mass between the models and SN 1987A translates into ~0.2 mag 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 the I band indicate that strong line-blanketing effects are missing in the models. These would redistribute the V-band fluxes to the I-band ones, 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 mag.

## 7.2 Spectral comparison

Fig. 19 shows a comparison of the observed optical spectra against the model spectra of model D127_NT (red curve), model D127 (blue

<figure><img src="image_17.png" /><figcaption>**Figure 20.** Same as Fig. 19, but it now shows a comparison of the Hα profile in the velocity space between observations (solid curve) and model D127_NT (dashed curve).</figcaption></figure>

curve) and the model with Fe I (hereafter model D127_NT_FeI (green curve). 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-mag offset in the V-band magnitudes in the previous section. The inclusion of Fe I improves the fit, mainly in the wavelength range 3000–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-hand emission wing of the profile fits the observation, but the Hα profile in the observation extends at least 3000 km s⁻¹ farther 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 the observational minimum is at 5000 km s⁻¹.

The full width at half-maximum 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 $^{56}$Ni is abundant. The spatial distribution of $^{56}$Ni depends on the degree of mixing. Mixing effects in SN explosions 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 $^{56}$Ni. 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

<figure><img src="image_18.png" /><figcaption>**Figure 19.** Comparison between the 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 photometric (4.37) V magnitudes. The three model spectra are scaled by 0.07/0.084 according to the $^{56}$Ni mass in SN 1987A and in the model.</figcaption></figure>
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 roughly the same time since the explosion as the observations. In Fig. 21, the observations are shown in the left-hand panel and the synthetic spectra are shown in the right-hand 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\alpha$ profile at late times. In the non-thermal sequence, the strength of the $H\alpha$ 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, the Na I D lines show behaviour similar to $H\alpha$, 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

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 which neglects the short range interaction between a perturbing electron and an 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, Burke & Dunseath 1994; Ramsbottom et al. 2007; Ramsbottom 2009).

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.

<figure><img src="image_19.png" /><figcaption>**Figure 21.** Left-hand panel: montage of observed optical spectra of SN 1987A. The days since breakout are shown below each spectrum. Right-hand panel: montage of synthetic spectra of the non-thermal sequence. Each observed spectrum in the left-hand panel 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.</figcaption></figure>
<table><caption>**Table 1.** Combinations of species to scale the cross-sections.</caption><thead><tr><th>Model</th><th>Excitation</th><th>Ionization</th></tr></thead><tbody><tr><td>EH5</td><td>H and He by 5</td><td>Same</td></tr><tr><td>EM5</td><td>Metal species by 5</td><td>Same</td></tr><tr><td>EA5</td><td>All species by 5</td><td>Same</td></tr><tr><td>IH5</td><td>Same</td><td>H and He by 5</td></tr><tr><td>IM5</td><td>Same</td><td>Metal species by 5</td></tr><tr><td>IA5</td><td>Same</td><td>All species by 5</td></tr></tbody></table>

In the case of varying the excitation cross-sections (Fig. 22, left-hand 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_e$ 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-hand 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 and 10 000 km s⁻¹ when the ionization cross-sections increase. This surprising behaviour is related to a low ejecta ionization fraction, $X_e$, shown in Fig. 22 (right-hand panel). At these low ionization fractions, the degradation spectra show that almost all electrons have degraded into very low energies and can no longer cause ionization. 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_e$. As $X_e$ 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 sizeable 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 5 in the uncertainties of the cross-sections will not dramatically influence our results.

## 8.2 $E_{max}$ for high-energy electrons

Radioactive decay from $^{56}$Ni and $^{56}$Co 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_{max}$ = 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 behaviours at high energy, so the uncertainty should be small as long as $E_{max}$ is large enough. We have run two additional models with $E_{max}$ = 2 and 0.5 keV. The difference in the fractional energies entering the three channels between one of the test models and the model with $E_{max}$ = 1 keV is less than 5 per cent, and the difference between the two test models is slightly larger. Importantly, the choice of $E_{max}$ 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_{max}$.

<figure><img src="image_20.png" /><figcaption>**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 FeI. Left-hand panel: the three test cases of scaling the impact excitation cross-sections. The scaling conditions are denoted in the figure. Different colours are fractions of different channels. The solid lines show the fractions of the original model. The dotted, dashed and dot-dashed curves represent the EH5, EM5 and EA5 models, respectively. Right-hand panel: the same as the left-hand panel, but it shows the three cases of scaling the impact ionization cross-sections. The dotted, dashed and dot-dashed curves show the channel fractions of the IH5, IM5 and IA5 models, respectively.</figcaption></figure>
## 8.3 The time-dependent 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 $^{56}$Ni decays is in LTE, and the assumption that all the energy goes into heating is excellent. However, we included non-thermal excitation and ionization at day 40.6, when the time-dependent effect is already functioning at large velocities, in regions of low ionization and underabundant in $^{56}$Ni/$^{56}$Co. We study this uncertainty in two ways. The first method is to compare the thermal and non-thermal models 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\alpha$ profile is only slightly weaker in the model with immediate inclusion of non-thermal processes. The above two tests indicate that 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 $^{56}$Ni is mixed outward to. The latest 3D simulations of CCSNe (Hammer, Janka & Müller 2010) showed that significant amounts of hydrogen are transported into deep layers of the ejecta ($\le 1000$ km s$^{-1}$), and that large amounts of heavy elements are strongly mixed above 2000 km s$^{-1}$. 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) modelled the late-time spectra of SN 1987A and found hydrogen was mixed down to $\le 700$ km s$^{-1}$ (they needed the mixing to fit the observations).

Fassia & Meikle (1999) modelled the He I 1.083 µm line and found some amount of $^{56}$Ni is required to be mixed above 3000 km s$^{-1}$ 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 are due to non-thermal effects, which required a substantial amount of $^{56}$Ni mixing out to 10 000 km s$^{-1}$ 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 $\theta$, $\phi$). 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 indicates macroscopic mixing. Fransson & Chevalier (1989) showed that microscopic mixing can significantly change the observed spectra at late times, and Li, McCray & Sunyaev (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 the explosion with a purely macroscopic mixing model. Furthermore, hydrodynamical models were unable to produce efficient microscopic mixing (Shigeyama & Nomoto 1990; Mueller, Fryxell & Arnett 1991; Herant & Woosley 1994).

Modelling 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 ways to reproduce the Balmer lines. However, the problem is complicated by $\gamma$-ray transport. $\gamma$-rays interact with the medium through three processes: photoelectric ionization, Compton scattering and pair production. The decay of $^{56}$Ni and $^{56}$Co mainly produces $\gamma$-rays with energy $\sim 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, grey opacity $\kappa_{\gamma} = 0.06 Y_e$ cm$^2$ g$^{-1}$ (Swartz, Sutherland & Harkness 1995), where $Y_e$ is the total number of electrons per baryon. The optical depth is calculated by

$$d\tau = -\kappa_{\gamma} \rho \, ds. \qquad (14)$$

Assuming $Y_e = 0.5$ and integrating $\tau$ from velocity 1500 km s$^{-1}$, the place where $^{56}$Co is still abundant, we obtain the mean free path of $\gamma$-rays as $\sim 1.0 \times 10^{14}$ cm, which means the photon can travel from velocity 1500 to 1700 km s$^{-1}$ in model D127_NT. No doubt, very few $\gamma$-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 $\gamma$-ray transport (Hillier & Dessart 2012), which will give a more detailed understanding of the mixing problem. We ran a test with allowance for $\gamma$-ray transport at day 127 and, as expected, there is only a slight enhancement in the line strength and width of $H\alpha$.

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 $\gamma$-rays was made by Matz et al. (1988), shortly after the discovery of hard X-rays. This suggests that $\gamma$-ray transport was beginning to become important at day 127. However, our model at day 154 gives an estimate that $\gamma$-rays can only travel from 1500 to 1800 km s$^{-1}$ using the effective grey opacity $\kappa_{\gamma}$, which indicates an insufficient outward mixing of $^{56}$Ni in the input. In fact, theoretical modelling of bolometric, X-ray and $\gamma$-ray light curves requires substantial outward mixing of $^{56}$Ni into the hydrogen envelope (Kumagai et al. 1988; Arnett & Fu 1989).

## 10 CONCLUSIONS

We developed a non-thermal solver, which takes into account ionization and excitation due to non-thermal electrons created by $\gamma$-rays that arise from the decay of $^{56}$Ni and its daughter products, and incorporated it into our fully non-LTE time-dependent code to model SNe. We use the model ‘Im18a7Ad’ of Woosley at day 0.27, 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_e$. A spectral comparison with the thermal models shows that $H\alpha$ is strongly increased at nebular epochs. At late times, He I lines, which are totally absent in the thermal models, are present in the 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 the 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 ionization balance.

We also compare the non-thermal models with the observations. The rebrightening of the bolometric light curve peaks about 10 d later in our models. An underestimate in the amount of outward mixing of ⁵⁶Ni is the main reason for the delay. Although the synthetic bolometric light curve agrees reasonably with the observations, multiband light curves show discrepancies in various degrees and reveal potential differences between the hydrodynamical model and SN 1987A or possible problems in our modelling. 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\alpha$ maintains a strong profile at the nebular epochs in the non-thermal models, resembling observations of SN 1987A. However, the model $H\alpha$ profile is much narrower than what we observed in SN 1987A, probably due to insufficient outward mixing of ⁵⁶Ni. Although hydrogen is abundant above 2000 km s⁻¹ in the model, the lack of ⁵⁶Ni 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 cut-off $E_{max}$ 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_{max} = 1000$ eV in our model produces almost the same spectra as using $E_{max} = 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 the photospheric to nebular phases continuously. The influence of non-thermal ionization and excitation on Type Ib and Type Ic SN

models is being investigated (Dessart et al. 2012). As more and more nebular spectra are available, 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.

## ACKNOWLEDGMENTS

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 Spectra Database, 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 which these calculations would not have been possible.

## REFERENCES

* Arnaud M., Rothenflug R., 1985, A&AS, 60, 425

* Arnett W. D., Fu A., 1989, ApJ, 340, 396

* Arnett W. D., Bahcall J. N., Kirshner R. P., Woosley S. E., 1989, ARA&A, 27, 629

* Bautista M. A., 1997, A&AS, 122, 167

* Burke P. G., Burke V. M., Dunseath K. M., 1994, J. Phys. B: At. Mol. Phys., 27, 5341

* Cardelli J. A., Clayton G. C., Mathis J. S., 1989, ApJ, 345, 245

* Chugai N. N., 1991, Sov. Astron. Lett., 17, 400

* Colgan S. W. J., Haas M. R., Erickson E. F., Lord S. D., Hollenbach D. J., 1994, ApJ, 427, 874

* Dessart L., Hillier D. J., 2005, A&A, 437, 667

* Dessart L., Hillier D. J., 2006, A&A, 447, 691

* Dessart L., Hillier D. J., 2008, MNRAS, 383, 57

* Dessart L., Hillier D. J., 2010, MNRAS, 405, 2141

* Dessart L., Hillier D. J., 2011, MNRAS, 410, 1739

* Dessart L. et al., 2008, ApJ, 675, 644

* Dessart L., Hillier D. J., Livne E., Yoon S.-C., Woosley S., Waldman R., Langer N., 2011, MNRAS, 414, 2985

* Dessart L., Hillier D. J., Li C., Woosley S., 2012, MNRAS, 424, 2139

* Fassia A., Meikle W. P. S., 1999, MNRAS, 302, 314

* Fransson C., Chevalier R. A., 1989, ApJ, 343, 323

* Gehrz R. D., Ney E. P., 1990, Proc. Natl. Acad. Sci. USA, 87, 4354

* Haas M. R., Erickson E. F., Lord S. D., Hollenbach D. J., Colgan S. W. J., Burton M. G., 1990, ApJ, 360, 257

* Hammer N. J., Janka H.-T., Müller E., 2010, ApJ, 714, 1371

* Harkness R. P. et al., 1987, ApJ, 317, 355

* Herant M., Woosley S. E., 1994, ApJ, 425, 814

* Hillebrandt W., Hoeflich P., Weiss A., Truran J. W., 1987, Nat, 327, 597

* Hillier D. J., 1987, ApJS, 63, 947

* Hillier D. J., Dessart L., 2012, MNRAS, 424, 252

* Hillier D. J., Miller D. L., 1998, ApJ, 496, 407

* Jerkstrand A., Fransson C., Kozma C., 2011, A&A, 530, A45

* Kozma C., Fransson C., 1992, ApJ, 390, 602

* Kozma C., Fransson C., 1998a, ApJ, 496, 946

* Kozma C., Fransson C., 1998b, ApJ, 497, 431

* Kumagai S., Shigeyama T., Nomoto K., Itoh M., Nishimura J., 1988, A&A, 197, L7

* Kurucz R. L., 2009, in Hubeny I., Stone J. M., MacGregor K., Werner K., eds, AIP Conf. Ser. Vol. 1171, Recent Directions in Astrophysical Quantitative Spectroscopy and Radiation Hydrodynamics. Am. Inst. Phys., New York, p. 43

* Landolt A. U., 1992, AJ, 104, 340
Langer N., El Eid M. F., Baraffe I., 1989, A&A, 224, L17

Li H., McCray R., 1995, ApJ, 441, 821

Li H., McCray R., Sunyaev R. A., 1993, ApJ, 419, 824

Lucy L. B., 1991, ApJ, 383, 308

Matz S. M., Share G. H., Leising M. D., Chupp E. L., Vestrand W. T., 1988, Nat, 331, 416

Meikle W. P. S., Spyromilio J., Varani G.-F., Allen D. A., 1989, MNRAS, 238, 193

Mitchell R. C., Baron E., Branch D., Lundqvist P., Blinnikov S., Hauschildt P. H., Pun C. S. J., 2001, ApJ, 556, 979

Mueller E., Fryxell B., Arnett D., 1991, A&A, 251, 505

Opal C. B., Peterson W. K., Beaty E. C., 1971, J. Chem. Phys., 55, 4100

Panagia N., 2005, in Marcaide J.-M., Weiler K. W., eds, IAU Colloq. 192, Cosmic Explosions, On the 10th Anniversary of SN1993J. Springer, Berlin, p. 585

Pelan J., Berrington K. A., 1997, A&AS, 122, 177

Phillips M. M., Heathcote S. R., 1989, PASP, 101, 137

Phillips M. M., Heathcote S. R., Hamuy M., Navarrete M., 1988, AJ, 95, 1087

Podsiadlowski P., Joss P. C., Rappaport S., 1990, A&A, 227, L9

Ramsbottom C. A., 2009, At. Data Nucl. Data Tables, 95, 910

Ramsbottom C. A., Hudson C. E., Norrington P. H., Scott M. P., 2007, A&A, 475, 765

Shigeyama T., Nomoto K., 1990, ApJ, 360, 242

Shull J. M., van Steenberg M. E., 1985, ApJ, 298, 268

Spencer L. V., Fano U., 1954, Phys. Rev., 93, 1172

Suntzeff N. B., Bouchet P., 1990, AJ, 99, 650

Sunyaev R. et al., 1987, Nat, 330, 227

Swartz D. A., 1991, ApJ, 373, 604

Swartz D. A., Filippenko A. V., Nomoto K., Wheeler J. C., 1993, ApJ, 411, 313

Swartz D. A., Sutherland P. G., Harkness R. P., 1995, ApJ, 446, 766

Utrobin V. P., 2004, Astron. Lett., 30, 293

Utrobin V. P., Chugai N. N., 2005, A&A, 441, 271

van Regemorter H., 1962, ApJ, 136, 906

Weaver T. A., Zimmerman G. B., Woosley S. E., 1978, ApJ, 225, 1021

Xu Y., McCray R., 1991, ApJ, 375, 190

Younger S. M., 1981, J. Quant. Spectrosc. Radiat. Transfer, 26, 329

This paper has been typeset from a *TEX/LATEX* file prepared by the author.