© ESO, 2011

1. Introduction

The distance of only  ~50 kpc to SN 1987A makes a number of unique observations possible. In particular, the supernova can after more than 20 years still be observed as resolved ejecta, glowing from radioactive input. About four years after exploding, the supernova entered the ⁴⁴Ti-dominated phase. Since then, its emission has been that of a cool (~10² K) gas powered by radioactive decays, complemented by freeze-out emission from the hydrogen envelope. Owing to the long life-time of ⁴⁴Ti (~85 years), the character of the spectrum has changed little since this transition.

To understand and interpret the late-time emission from the ejecta, detailed spectral synthesis modeling is needed. The spectral formation process in supernovae is complex, especially at late times when all important processes are non-thermal, but radiative transfer effects may still be important. High-energy gamma-rays and positrons from the radioactive decays deposit their energy into free and bound electrons, producing a population of fast electrons which heat, excite and ionize the gas. The ionizations produce secondary electrons that additionally contribute to these processes. The degradation process is quantified by solving the Boltzmann equation for each of the zones present in the supernova (e.g. Xu & McCray 1991; Kozma & Fransson 1992, KF92 hereafter). Once the deposition in the various channels have been determined, the temperature, ionization structure and NLTE level populations can be computed by solving the equations for thermal and statistical equilibrium or their time-dependent variants (e.g. de Kool et al. 1998; Kozma & Fransson 1998a,b, KF98 a,b hereafter). Because the ionization balance affects the solution for the radioactive deposition, the equation systems are coupled and a solution is found by iteration.

The equilibrium also depends on the radiation field, which should simultaneously be solved for. The strong velocity gradients in supernovae imply unique radiative transfer effects. As they travel through the ejecta, photons are continuously red-shifting with respect to the comoving frame, and are therefore exposed to absorption in lines over a broad wavelength range. Because especially iron-group elements have a large number of lines at UV/blue wavelengths, radiation here experiences a complex transfer and a significant fraction of it emerges by fluorescence as a quasi-continuum at longer wavelengths (Li & McCray 1996). The efficient blocking by atomic lines is equivalent to an effective continuum opacity known as line blocking. The emerging UV emission may often be due to “holes” in this line blocking where radiation can escape (Mazzali & Lucy 1993). Although the optical depths decrease over time, the multi-line transfer is important for modeling supernova spectra for many years or even decades after explosion.

The observational determination of the masses of the three main radioactive isotopes ⁵⁶Ni, ⁵⁷Ni, and ⁴⁴Ti constitutes one of the main constraints on explosion models of core collapse supernovae. The production of these isotopes is sensitive to both density and temperature, and thus to the explosion dynamics (e.g. Woosley & Hoffman 1991). During the first  ~500 days the decay of ⁵⁶Co (which is the product of rapid ⁵⁶Ni decay) dominated the bolometric light curve of SN 1987A. From the light curves, the ⁵⁶Ni mass could be determined to be 0.069     M_⊙ (Bouchet et al. 1991). The determination of the ⁵⁷Ni mass was complicated by the fact that when this isotope became important, the bolometric light curve was already affected by the delayed release of ionization energy (freeze-out). A time-dependent modeling that took this into account resulted in an estimated ⁵⁷Ni mass of  ~3.3 × 10^-3   M_⊙ (Fransson & Kozma 1993, FK93 hereafter), which agrees well with the mass derived from observations of the infrared Co II/Fe II lines and the gamma-rays emitted in the decay (Varani et al. 1990; Kurfess et al. 1992).

The determination of the ⁴⁴Ti-mass is complicated by a number of factors, the first being the freeze-out contribution just as for ⁵⁷Ni. Second, at late times the low temperature of the gas and the presence of dust lead to most of the deposited energy being re-radiated in the largely unobservable far-infrared, which prohibits a determination of the bolometric luminosity. The ⁴⁴Ti-mass must therefore be determined from detailed modeling of the fraction of the energy that emerges in the (observable) UV/optical/NIR bands, and must consider the time-dependent freeze-out. From this type of modeling, Fransson & Kozma (2002, FK02 hereafter) determined the ⁴⁴Ti-mass to (0.5−2) × 10^-4 M_⊙, which agrees with the (1−2) × 10^-4 M_⊙ estimated from a nebular analysis of the eight-year spectrum by Chugai et al. (1997, C97 hereafter). Lundqvist et al. (2001) derived an upper limit of 1.1 × 10^-4 M_⊙, based on the non-detection of the [Fe II] 26 μm line by ISO, which rested on the assumption of local positron deposition and the absence of any dust cooling.

Another possibility for determining the ⁴⁴Ti-mass is by direct detection of the gamma-rays produced in the decay. Only in one case, Cas A, has such a detection been made (Iyudin et al. 1994; Vink et al. 2001), with the mass estimated to (Renaud et al. 2006). For Cas A we do not know the masses of the ⁵⁶Ni and ⁵⁷Ni isotopes though. For SN 1987A, the detection limits by INTEGRAL resulted in only a weak upper limit of 1 × 10^-3     M_⊙ for the ⁴⁴Ti-mass (Shtykovskiy et al. 2005).

In this paper we present a detailed model of the UV/optical/NIR spectrum of SN 1987A at an age of eight years, using a 19 M_⊙ explosion model as input. Using a new code with a detailed radiative transfer treatment, our objectives are to see if the spectral model can reproduce the main features in the observed spectrum, to understand the contributions by various zones and elements to the spectral formation process, to quantify the importance of the (non-local) line transfer, and to determine the ⁴⁴Ti-mass as accurately as possible. Previous discussions of the late-time ejecta spectra can be found, e.g., in Wang et al. (1996) and C97. Our paper also serves to describe the code we developed for reference in future modeling.

2. Modeling

In KF98 a,b, a self-consistent model for the nebular phase spectrum of SN 1987A was presented. This included a detailed calculation of the gamma-ray/positron thermalization and a determination of the time-dependent ionization, excitation, and temperature structure of the ejecta. Comparing the resulting spectra from two different explosion models, the evolution of all major emission lines were analyzed. Fransson & Kozma (2002) contains an update of this work, where the modeling of some of the most important lines and the broad band photometry were discussed.

A shortcoming in these models was that (non-local) radiative transfer in atomic lines was not included, which introduced an uncertainty for the internal radiation field and thereby for the ionization balance, and also for the emerging spectrum. Here, we calculate in detail the line scattering and fluorescence throughout the spectrum (Sect. 2.5), which significantly improves the accuracy of the model. In addition, we updated the atomic physics by adding new atomic data and including more and larger model atoms.

While we do a more accurate radiative transfer calculation, this comes at the expense of time-dependent effects as our model assumes steady-state. In FK93 it was shown that in the envelope the recombination time scale became comparable to the radioactive decay time scale at  ~800 days. From that time on, much of the ionization energy was not radiated instantaneously but on a longer time scale, and a time-dependent calculation was therefore necessary. In the ⁴⁴Ti-dominated phase the radioactive time scale becomes much longer,  ~85 years (Ahmad et al. 2006), and is less relevant. Instead, it is the adiabatic expansion time scale, i.e. the age of the supernova, that is the most relevant, and should be compared to the recombination and cooling time scales. In zones where a time-dependent treatment is necessary, we use the solutions obtained by an updated version of the time-dependent code used in KF98 a,b for the ionization balance. As we show below (Sect. 4.2), this is only necessary in the hydrogen envelope.

2.1. Ejecta model

We base our ejecta model on the 19 M_⊙ explosion model computed by Woosley & Heger (2007, WH07 hereafter). From the distribution of elements in this model, we construct seven types of zones; Fe/He, Si/S, O/Si/S, O/Ne/Mg, O/C, He, and H, where the zones are named after their dominant elements. The mass cuts for the zones are taken to be where the most common or second most common element changes. The mass of the Fe/He zone is adjusted to give a total ⁵⁶Ni mass of 0.072 M_⊙ in the ejecta (the value of 0.069 M_⊙ from Bouchet et al. 1991, corrected for the slighly larger distance we use here), and the mass fractions of ⁵⁷Co and ⁴⁴Ti in this zone are then adjusted to give a total ⁵⁷Co mass of of 3.3 × 10^-3 (FK93) and a total ⁴⁴Ti mass of 1.5 × 10^-4   M_⊙, which we later show to be our favored value for this isotope. The H zone mass is adjusted to match the observed estimate of  ~8 M_⊙ (KF98 b). Finally, we replace the He/C zone in the explosion model with a He/N zone, that has not experienced any outward mixing of carbon. This is also based on modeling in KF98 b, where it was found that a significant He/C zone would result in [C I] lines much stronger than the observed ones.

The masses and compositions of the zones in the model are given in Table 1. Abundances are taken as the value in the middle of each zone. While this does not exactly conserve the total element masses, the differences are small compared to using zone-averaged values, and the middle values should give a more consistent description for the typical composition. Because the WH07 models are based on solar metallicity progenitors, we replace the abundances of all elements with mass fractions differing by less than a factor 1.25 from the solar values with the LMC abundances taken from Russell & Dopita (1992). A comparison of models with different metallicities in Woosley & Weaver (1995, WW95 hereafter) shows that most nucleosynthetic yields only weakly depend on metallicity, so the abundances for the processed elements should be very similar for solar and LMC metallicity progenitors. Observations of the circumstellar material indicate that CNO-burning products have been mixed out from the helium core into the hydrogen envelope (Fransson et al. 1989). We use the abundances of He, N and O for the envelope derived by Mattila et al. (2010) and of C by Lundqvist & Fransson (1996, see Table 1).

We model the supernova as a spherically symmetric, homologously expanding nebula, divided into a core and an envelope. The core contains the heavy element zones (Fe/He, Si/S, O/Si/S, O/Ne/Mg, O/C), and fractions of the He and H zones, which are mixed inwards during the explosion. These fractions were chosen to be 25% (see Sect. 2.1.1). The envelope contains the remaining 75% of the He and H zones.

The extent of the core can be estimated from the widths of the lines that are emitted mainly from the heavy element zones, for example the iron lines and [O I] λλ6300, 6364. The oxygen lines have maximum expansion velocities of  ~1700 km s^-1, whereas the iron, cobalt, and nickel lines extend to  ~2500 km s^-1, with some contribution even out to  ~3500 km s^-1 (see McCray 1993, for a review). Our single-core model forces a compromise value, which we choose as V_core = 2000   km s^-1, the same as in KF98 a,b. A higher core velocity means less gamma-ray deposition, but because the positrons dominate the energy input at late times, this choice is not critical for the energy budget here.

2.1.1. The core

We assume that the zones in the core experience complete macroscopic mixing, where they fragment into N_cl clumps that become completely mixed in velocity space (see Fig. 1). Support for this macroscopic mixing comes from theoretical considerations (e.g. Mueller et al. 1991), multi-dimensional simulations (e.g. Herant & Benz 1991; Mueller et al. 1991; Kifonidis et al. 2006; Hammer et al. 2010), the early emergence of X-rays and gamma-rays (Dotani et al. 1987; Pinto & Woosley 1988), and from the similar line profiles of many of the elements (e.g. Spyromilio et al. 1990). We use a value of N_cl = 100 based on the analysis of the iron clumps in Li et al. (1993). In Sect. 4.6 we find that the model spectrum is insensitive to this choice. We further assume that no microscopic mixing occurs (e.g. Fransson & Chevalier 1989; Fryxell et al. 1991), so that the zones retain their compositional integrity.

The clumpy structure of the core calls for a special handling of the computational grid. We use a gridless domain for the core region; when a clump emits photons, the position of the clump within the core is selected by a random draw. When we follow photons through the core (Sect. 2.5), exiting one clump leads to a random draw that determines which type of clump to encounter next. We set the probability of encountering any given type of zone to be proportional to the total surface area of that zone type. A random number also selects the impact parameter on the new clump. The path traveled through zone i is then proportional to A_i·R_i ∝ V_i, which recovers the desired property that a photon should spend its time in a given zone in proportion to the filling factor of the zone, ϵ_i.

The core structure is a “virtual” one in the sense that there is no fixed density grid; we have a set of spherical clumps which together make up the total volume of the core, but which obviously cannot be exactly fit together since they are all spheres. There are several advantages to this virtual grid. The extension of the Monte Carlo philosophy to encompass also a randomization of the grid is a powerful and appealing way to make computations in clumpy media such as this. While the radiative transfer can still be done realistically by the use of surface and impact probabilities as described above, a virtual grid naturally offers a way to capture the situation of macroscopic mixing and to parametrize the degree of fragmentation of the core. Finally, the line profiles automatically become smooth rather than jagged, as when using shells.

The filling factors for the oxygen zones are set from the oxygen number density derived from the evolution of the [O I] 6300 Å/6364 Å line ratio, n_O = 6.2 × 10¹⁰   cm^-3 at 100 days (Spyromilio & Pinto 1991; Li & McCray 1992). This allocates  ~10% of the core volume to oxygen clumps.

The iron clumps are believed to expand because of their content of radioactive materials (Woosley 1988; Herant & Benz 1991; Basko 1994), which was verified for SN 1987A by modeling of the cooling lines from the iron clumps (Li et al. 1993). The iron zone was found to have a filling factor of ϵ_Fe/He ≳ 0.3, assuming V_core = 2500 km s^-1. The same density limit using V_core = 2000 km s^-1 corresponds to ϵ_Fe/He ≳ 0.6. We use a value of 0.6 in our modeling here.

In the explosion model, the Si/S zone contains  ~5% of the ⁵⁶Ni, and can therefore be expected to expand as well. For this zone, we therefore use a density in between the oxygen zones and the iron zone, chosen to be ten times lower than the oxygen zone density. This gives a filling factor of  ~3%. In Sect. 4.3.2, we find support from our modeling for a low density in this zone.

Hydrodynamical simulations by Fryxell et al. (1991) find  ~1 M_⊙ of hydrogen mixed inside 2000 km s^-1, while Herant & Benz (1992) find  ~2 M_⊙. From modeling of the Hα line in SN 1987A, KF98 b estimate the presence of 2.2 M_⊙ of H-zone gas inside 2000 km s^-1. Here, we take 2 M_⊙ of H-zone material to be mixed inside the core, which corresponds to  ~25% of the total H zone mass. We assume the same fraction (25%) of the He zone to be mixed into the core as well, giving a He-zone core component of 0.3 M_⊙ and a He-envelope component of 1.0 M_⊙. The total mass of the core, including the mixed-in He and H, is then 5.2 M_⊙.

The filling factors for the H and He zones in the core are set by allocating the remaining volume so that these zones obtain equal number densities. This gives the H component a filling factor of  ~26% and the He component a filling factor of  ~2%. See also KF98 a for further discussions regarding filling factors.

2.1.2. The envelope

Outside the core we attach an envelope consisting of an inner He zone followed by logarithmically spaced (constant dR/R) H-zone shells. As mentioned above, the total He-zone mass is 1.3 M_⊙, of which we put 1.0 M_⊙ in the envelope and 0.3 M_⊙ in the core. The total H-zone mass is taken as 8 M_⊙ from Hα modeling (KF98 b), of which we put 6 M_⊙ in the envelope and 2 M_⊙ in the core (see Sect. 2.1.1).

We use a density distribution of the H-envelope based on hydrodynamical modeling by Shigeyama & Nomoto (1990) (model 14E1). The validity of this model is confirmed by modeling of the Mg II λ2800 and Mg I λ2852 features in C97. The density profile can be described as gradually steepening from an initial −2 power law to an asymptotic  ~− 8 for V > 5000 km s^-1. We terminate the envelope at a velocity of 10 000 km s^-1.

2.2. Energy deposition and degradation

We calculate the deposition of gamma-rays emitted by the clumps containing ⁵⁶Co, ⁵⁷Co and ⁴⁴Ti. Based on the simulations by Colgate et al. (1980); Woosley et al. (1989) and Fransson & Chevalier (1989) we assume absorption with effective opacities of ,  and  cm² g^-1, where and are the average nuclear charges and atomic weights in each zone. We do not include other radioactive isotopes such as ⁶⁰Co and ²²Na, whose abundances are highly uncertain. We discuss the possible effects of these in Sect. 5.

From each zone containing radioactive material, we emit and track ⁵⁶Co, ⁵⁷Co, and ⁴⁴Ti gamma-ray packets as described for the UV/optical/NIR radiation in Sect. 2.5, with the difference that these packets travel in straight paths until they reach the edge of the nebula, depositing a fraction exp(− κρ_il) of their energy for each traveled distance l, where ρ_i is the density of zone i. We use radioactive luminosities taken from KF98, with the ⁴⁴Ti values updated for a life-time of 85 years instead of 78 years.

After  ~2000 days, ⁴⁴Ti completely dominates the deposition. The degradation of the deposited gamma-ray and positron energy into heating, ionizations and excitations is obtained by solving the Spencer-Fano equation using the routine developed by KF92. For this routine, we updated the collisional cross sections for Fe II (Ramsbottom et al. 2005, 2007) and Ca I (Samson & Berrington 2001). Ionizations are generally assumed to leave the ions in their ground states, except for O II, for which we compute specific rates to the first few excited states.

In our standard model, the positrons are assumed to be absorbed on the spot. As noted before, e.g. in C97, this requires a magnetic field. The average positron energy in the ⁴⁴Sc decay is 0.6 MeV and the maximum energy of the positrons is 1.47 MeV (Browne et al. 1986). The energy loss per unit length, dE/dx, is dominated by excitations and ionizations and is given by the Bethe formula (1)where r₀ is the electron radius, β = V/c, I is the effective ionization potential of the atom, f_z is a relativistic correction factor, and all other symbols have their usual meaning. In Fig. 2 we show the stopping range,  (g cm^-2), for the most important elements in the core together with the total range, averaged over the different zones in the core.

These ranges should be compared to the mass column density D of the various zones in the core and envelope. For the total core, this is given by (2)The column density for a core with V_core = 1800 km s^-1 at 7.87 years is shown as the upper dashed line in Fig. 2. For comparison we also show the column density of the Fe/He zone as the lower horizontal line.

That the column density of the Fe/He zone is nearly two orders of magnitude lower than the stopping range for Fe and He means that in the absence of a sufficiently strong magnetic field the positrons will not be trapped in this zone, even if it is not strongly fragmented into many small pieces. The figure also shows, however, that the total core has a sufficiently high column density to trap the positrons even without a magnetic field. There should therefore be no positrons entering the envelope at this or earlier epochs.

The presence of magnetic fields in supernova ejecta has been suggested many times, see e.g. Axelrod (1980). Colgate et al. (1980) argue that any magnetic field present at the time of the explosion becomes radially “combed” and will therefore not be important for the trapping (see also Ruiz-Lapuente & Spruit 1998). While the positrons should be confined to the core at eight years, it is therefore not clear whether they are absorbed by the zones producing them (which is mainly the Fe/He zone, but also the Si/S zone), or if they can propagate into the other core zones before they are absorbed. We will later investigate this question by comparing the spectra produced in the different scenarios (Sect. 4.5).

2.3. Statistical and thermal equilibrium

For each zone we solve the equations for statistical equilibrium with respect to the ionization and excitation structure, as well as the equation of thermal equilibrium. These equations are essentially the ones described in KF98 a, with some modifications and improved rate calculations described below. For each iteration, the algorithm deployed is to first compute the ionization balance, then iterate solutions to excitation structure and thermal equilibrium until the temperature stabilizes, and finally to calculate the radiative transfer (Sect. 2.5). From the radiative transfer, new photoionization rates, photoexcitation rates (although these are set to zero in this work, see Sect. 2.3.2), and radiative heating rates are calculated, and used for the new solutions for ionization, excitation and temperature in the next iteration. Iteration continues until the maximum change in any zone temperature, electron fraction or ion abundance is less than 1%. We then also verify that the change in the emerging radiation field is neglegible.

2.3.1. Ionization balance

Steady-state for the ionization balance is appropriate as long as the recombination time-scale (1/n_eα) is short compared to the evolutionary and radioactive time scales. For the densities at the epoch studied here, this is satisfied (Sect. 4.2) for all core zones as well as for the envelope helium zone. The hydrogen envelope, however, experiences freeze-out after  ~800 days (FK93). For the envelope, we therefore use the solutions for the hydrogen ionization fraction from a time-dependent calculation, based on the model in KF98 a, given in Table 3. We do, however, use steady-state solutions for the other elements in the envelope, because the time-dependent model lacks photoionizations by line radiation, and may therefore significantly underestimate the ionization degree for some of them. We also take the hydrogen envelope temperatures from the time-dependent solution, because temperature and ionization balance are coupled.

For all other zones, the ionization balance is computed from the steady-state variant of Eq. (10) (and its generalized version) in KF98 a. We include the first three ionization stages of all elements listed in Sect. 2.3.2. The main improvement in our treatment here is that we now compute detailed photoionization rates from the Monte Carlo simulation (Sect. 2.5). For zone i and atom k, the photoionization rate is (3)where j is the level, N_i,k,j is the total number of photoionizations for zone i, atom k, level j, in the Monte Carlo simulation, V_i is the zone volume and is the previous solution to the number density of level j.

2.3.2. Excitation structure

Given the ion abundances, steady-state is appropriate for the atomic level populations since the collisional and radiative time-scales are short. We compute NLTE solutions for the neutral and singly ionized elements of H, He, C, N, O, Ne, Na, Mg, Al, Si, S, Ar, Ca, Sc, Ti, V, Cr, Mn, Fe, Co, and Ni, as well as for doubly ionized Fe, using model atoms that include fine-structure levels. With their large number of resonance lines in the UV and optical, iron-group elements are especially important for the radiative transfer, which motivated us to include also those with small abundances. Appendix A gives references for the atomic data that we use.

For the line optical depths, we use the Sobolev approximation (Sobolev 1957) (4)with the corresponding escape probabilities (5)where i and j are the lower and upper levels, respectively. The Sobolev approximation is valid as long as the region of line interaction is much smaller than the region over which physical conditions change, which is generally fulfilled in supernovae. The errors caused by line-overlap should be small in the nebular phase (Li & McCray 1996). For lines with small , we also take into account continuum destruction probabilities (Hummer & Rybicki 1985; Chugai 1987), so that the effective escape probability is (6)where the probabilities are calculated as in KF98 a. The effective radiative deexcitation rates used in the NLTE equations are then .

The statistical equilibrium is obtained by solving Eq. (15) in KF98 a. Deexcitation of the He I(2³S) state also occurs through Penning ionizations of H I. Our treatment here differs mainly in that we now compute detailed photoionization rates P_i,k,j from the Monte Carlo calculation (Sect. 2.5) (7)As in KF98, we include non-thermal ionizations and photoionizations only to the first level of the ionized element. The only exception is for O II, where non-thermal ionizations to the excited states , and ⁴P are also included. Non-thermal ionizations are included only from the ground multiplets in the atoms, whereas photoionizations can occur from excited states where cross sections are available.

In this work, we set all photoexcitation/deexcitation rates to zero. Although these can in principle be computed from the radiative transfer calculation and be included in the population equations, this would force a slow single-stepping Monte Carlo transfer. We also aim to study the scattering/fluorescence process in detail (Sect. 4.3), which is easier if we can track the Monte Carlo packets through all scattering events until they escape. Finally, this approximation also allows us to compute faster (approximate) NLTE solutions including a fewer number of states. At these late times, photoexcitations do not strongly influence the solutions beyond the direct scattering/fluorescence.

For the treatment of charge transfer, see Sect. 2.4.

2.3.3. Thermal balance

For the core zones, the temperature is computed by balancing heating with cooling. Steady-state is appropriate as long as adiabatic cooling is neglegible. This may not be fulfilled at late times for the oxygen, helium, and hydrogen zones (KF98 a). However, we compared the steady-state temperatures in these zones with those obtained from the time-dependent model (KF98 a,b), and they are found to be very similar. We therefore use the steady-state approximation for all core zones. For the envelope zones, we use the temperatures calculated in the time-dependent model. These are given in Table 3.

Heating has contributions from non-thermal deposition, photoionizations, collisional deexcitations (occuring in the Monte Carlo calculation), charge transfer reactions, free-free absorptions, and Penning ionizations (8)The terms H_pi, H_cd and H_ff are computed from the Monte Carlo simulation (Sect. 2.5) by summing all contributions from Eqs. (23), (24) and (34). The charge transfer heating is (Kingdon & Ferland 1996) (9)where is the reaction rate between atoms k and k′, and ΔE_k,k′ is the energy defect of the reaction, which is negative for endothermic reactions, which act as coolers.

The heating from Penning ionizations is, based on the rate in Bell (1970)(10)Cooling is computed from Eqs. (3)–(6) in KF98 a. It has contributions from line cooling, recombination cooling, and free-free emission.

2.4. Charge transfer

We use charge transfer rates taken from Rutherford et al. (1971), Rutherford et al. (1972), Arnaud & Rothenflug (1985), Pequignot & Aldrovandi (1986), Kimura et al. (1993), Swartz (1994), Kingdon & Ferland (1996), Stancil et al. (1998) and Zhao et al. (2004). For reactions where no measured or calculated rates exist, we use the recipe for estimating rates given by Pequignot & Aldrovandi (1986).

The many unknown or uncertain rates introduce an uncertainty for the level of ionization for some of the elements. For some zones, this prohibits a definite determination of which elements that re-emit the non-thermal and radiative ionization energy. We discuss this and the potential effects on the spectrum in Sects. 4.2 and 4.7.

When the final states are unknown, we assume these to be the ground states of the elements. An exception to this is the important O II + C I  →  O I + C II reaction, where we assume that the oxygen atom is left in the excited 2p(¹D) state, which is the one closest to resonance with an energy defect of 0.4 eV (Sect. 4.2.4). We also assume that the reaction Ca I + O II  →  Ca II + O I occurs to the excited 5p(²P) state in Ca II (Rutherford et al. 1972).

2.5. Radiative transfer

Radiative transfer is treated with a Monte Carlo technique. The application of this method for dealing with differentially expanding astrophysical envelopes was described by Abbott & Lucy (1985), where it was applied to the case of stellar winds. Lucy (1987) and Mazzali & Lucy (1993) describe the method for computing photospheric-phase supernova spectra. The main difficulty working with the equation of radiative transfer is the complexity introduced by the line overlaps that are caused by the differential velocity field. This gives a non-local coupling for which a Monte Carlo treatment is better suited. In addition, Monte Carlo radiative transfer is easily generalized to 2D and 3D, and is well suited for parallelization.

The early Monte Carlo codes mentioned above allowed for electron and resonance line scattering, where deexcitations occured in the same transition as the absorption. Lucy (1999) and Mazzali (2000) improved on this treatment by including fluorescence in all transitions directly connected to the (initial) upper level. Kasen et al. (2006) used the same treatment for the fluorescence in computing photospheric-phase, time-dependent spectra for type Ia supernovae. To allow for a more exact treatment of the fluorescence, but still retain the properties of indivisibility and indestructibility for the photon packets, a Monte Carlo formalism based on the concept of macro-atoms was developed in a series of papers by Lucy (2002, 2003, 2005). Applications and further developments of this method are described in e.g. Sim (2007), Sim et al. (2007) and Kromer & Sim (2009).

Here, we implemented a Monte Carlo algorithm with a similarly high level of physical realism by computing the complete deexcitation cascade following each photon packet absorption. A new photon packet is created for each radiative deexcitation in the cascade, with collisional deexcitations being allowed for as well. The (potential) creation of several packets following the absorption of one packet is handled by a recursive technique. We find that this works well and have therefore not enforced the property of indivisibility otherwise commonly used to simplify the bookkeeping (e.g. Lucy 2002). Also, the property of indestructibility devised in those papers has its main advantage as a Lambda accelerator in the photospheric phase, and is not enforced here.

2.5.1. Packet emission

Below, quantities are denoted by primed symbols in the comoving frame, and by unprimed symbols in the rest frame (observer frame). We ignore terms of order (V/c) or higher in the transfer, except for the critical Doppler shifts. The letter z denotes a random number between 0 and 1, newly sampled for each computation.

Total emissivities from lines and continua (two-photon, free-bound and free-free) are obtained for each zone i and wavelength binj from solutions to the ionization balance, temperature and NLTE level populations (Sect. 2.3). Photon packets are then created with energy (11)where V is the volume of the zone, is the width of the wavelength bin, and N_i,j is the number of packets to split the emission into. The wavelength λ′ is chosen as the maximum wavelength of bin j. This choice ensures that no line self-absorptions, which are already included in the Sobolev escape formalism, occur in the transfer. We currently use a uniform weighting for the number of packets emitted per zone and wavelength bin, N_i,j = N₀. From inspection of the emerging spectrum (which we also smooth, see Sect. 2.5.4), we find N₀ ~ 100 to be sufficient for obtaining neglegible Monte Carlo noise, for binning width Δλ′/λ′ = 10^-3.

The starting radius of the packet is sampled from (12)which corresponds to uniform sampling in a shell bounded by r_i,in = V_i,int and r_i,out = V_i,outt. The shell velocities are listed in Table 2.

Isotropic emission corresponds to sampling a direction cosine μ from (13)If the zone is a core zone, the internal radius in the (spherical) clump is also determined from (14)where R_i,clump is the radius of the clumps of zone type i. An internal direction cosine μ_int is also randomly drawn according to Eq. (13).

2.5.2. Packet propagation

For each packet to be emitted, a random number is drawn to determine at which optical depth the packet should be absorbed (if at all) (15)Next, the distance to the next critical point is computed. The comoving wavelength of the photon packet is continuously red-shifting as the packet travels through the differentially expanding ejecta. From the NLTE solutions, we have a wavelength ordered list of all lines with optical depth τ_ij > 10^-3. The velocity distance to the next line is computed as (16)where is the (atom rest frame) wavelength of the line closest to λ′ in that list.

The distance to the next ionization edge is computed as (17)where is the (atom rest frame) ionization edge closest in wavelength to λ′. Also, the distance to the next shell ΔV_shell is computed, as is the distance to the clump edge ΔV_clump, if the zone is a core zone. The velocity distance to the next critical point is then chosen as the smallest of these four distances (18)The continuum optical depth for this distance is computed as (19)where α_i,dust is the absorption coefficient for the dust (20) is the photoionization absorption coefficient, and α_i,ff(λ′) is the free-free absorption coefficient. Electron scattering is neglegible at the epoch studied here, but can easily be included as well.

The accumulated optical depth τ_tot is increased by τ_cont. If this causes τ_tot to transcend τ_abs, the packet is continuum-absorbed, and propagation terminates. The total number of photoionizations for zone i, atom k, level j, is then increased by (21)where N_packet is the number of photons in the packet (22)The photoelectric heating rate contribution is (23)where χ_k,j is the ionization potential for atom k, level j. The free-free heating rate contribution is (24)If continuum absorption does not happen, the position, flight angle, wavelength and energy of the packet are now updated. The new position is (25)The new flight angle is (26)The new comoving wavelength is (27)and the new comoving energy is (28)If the new point r_new corresponds to the passage of an ionization edge or a shell edge, new velocity distances are computed, starting from Eq. (16), and the procedure is repeated. If r_new corresponds to entering a new core clump, its type is chosen from surface impact probabilities given by (29)which is valid as long as all clumps have the same number of fragments (N_cl). An impact parameter is also drawn as (30)Propagation then continues into the new clump, starting from Eq. (16).

2.5.3. Line interaction

If r_new corresponds to reaching a line with τ_ij ≥ 10^-3, τ_tot is increased by τ_line. If τ_tot then transcends τ_abs, the packet is absorbed in the line.

All downward radiative and collisional transition rates are then computed. If u is the upper level of the line, the radiative deexcitation rates are (i < u) (31)The collisional deexcitation rates are (Osterbrock & Ferland 2006) (32)where Υ_u,i(T) is the effective collision strength. The rates are normalized to obtain the relative transition probabilities, and one of the deexcitation channels is selected by a random draw.

If this is a radiative transition (with a wavelength shorter than a pre-selected cut-off value, chosen as 2.1 μm here) a new photon packet is emitted. The energy of this packet is (33)where ΔE_ui is the energy difference between levels u and i and E_u is the energy of level u. A direction is randomly drawn from Eq. (13) , and the packet is followed, starting from Eq. (15), until it (and its potential offspring packets) are absorbed or escape the ejecta.

If the deexcitation is instead a collisional one, a contribution to the volumetric heating rate (34)is recorded. We ignore the possibility of collisional excitations in the cascade. This is a good approximation at late times when the low temperatures make upward collisions much rarer than downward ones.

New deexcitation probabilities from level i are now computed from Eqs. (31) and (32), and a channel is again randomly selected. The process repeats until the atom has deexcited to the ground state. The energy in each deexcitation from level i to level j is given by Eq. (33) with ΔE_ui replaced by ΔE_ij. This treatment means that (35)While this becomes only an approximate treatment for the energy distribution of the fluorescence when the absorption occurs from an excited state, this is rare at late times when most atoms are in their ground states. This treatment will likely be modified for calculations at earlier epochs when this may not be the case. Note that the outgoing energy (assuming only radiative deexcitations) is generally not conserved in the rest frame, although it is in the comoving frame. This corresponds to adiabatic losses of the radiation field.

2.5.4. Emerging spectrum

When it reaches the edge of the nebula (with expansion velocity V_e), the observed wavelength and energy of the escaping packet are obtained as (36)(37)where μ_e is the direction cosine of the packet upon crossing the outer edge. The escaping packets are binned according to λ_obs, using a wavelength bin size of Δλ/λ = 0.001. The final spectrum is obtained by smoothing this distribution with a box-car filter of width Δλ/λ = 0.002.

An alternative method for computing the emerging spectrum would be to solve for full-size atoms, including photoexcitation rates in the NLTE level population solutions, and compute the formal integral. However, computing large NLTE atoms is time-consuming and also numerically challenging in terms of convergence and accuracy. The line transfer mainly depends on the populations of the ground states and other low-lying states, which are more reliably computed than high-lying states (Lucy 1999). We therefore choose to compute smaller NLTE solutions without photoexcitation/deexcitation rates and treat the scattering/fluorescence process explicitly in the Monte Carlo transfer, which also allows us to investigate its properties from the emerging packet histories.

Because our model assumes steady state, we checked that the typical photon flight times t_p are short compared to the dynamic time scale. Crossing the ejecta takes a time of about  ~V_et/c for the photons, which means that t_p/t ~ V_e/c ≪ 1. However, UV photons may experience multiple line scatterings before emerging and their flight times could therefore become long. Here, the presence of dust in the ejecta limits the total flight times to a maximum of a few crossing times, and steady state is therefore justifiable from the radiative transfer perspective as well.

2.6. Dust

There is strong evidence from blue-shifting line peaks, a developing IR excess, and a simultaneous decrease in the optical flux that dust formed in the inner parts of the ejecta, starting around day 530 (Lucy et al. 1989; Wooden et al. 1993), or possibly as early as day 350 (Meikle et al. 1993). By day 700 most of the dust formation was probably over (Lucy et al. 1989). The apparently weak wavelength dependence of the dust extinction and the absence of any dust features in the IR emission suggest that the dust was at least partly formed in optically thick clumps (Lucy et al. 1991).

The last far-IR observations are from day 1731 (Bouchet et al. 1996), when the dust temperature was estimated to be  ~155 K. Fassia et al. (2002) find that essentially all the NIR ejecta lines have blue-shifted peaks even at day 2112, which indicates that the dust was still optically thick at that time. We determined the median value of their measured line-shifts to 500 km s^-1. This agrees well with Wang et al. (1996), who find blue-shifts of  ~400 km s^-1 for several optical lines at 1862 and 2210 days.

Between 2000 and 6000 days there were no instruments capable of detecting the FIR emission by the ejecta. Bouchet et al. (2004) identify a Gemini N filter detection at 6094 days with ejecta dust. Observations at 6526 days were dominated by the heated dust emission from the ejecta – ring collision, and no ejecta dust was detected (Bouchet et al. 2006). But the HST imaging shows a clear “hole” in the central parts of the ejecta even at these late times (Larsson et al. 2011), suggesting that the dust is still optically thick. This may, however, also be a result of external X-ray heating from the ring collisions (Larsson et al. 2011).

We examined the observed line profiles of Mg I] λ4571 and Hα in the eight-year spectrum (see Sect. 3 for a description of the observed spectrum). Both show distinct blue-shifts of their peaks of  ~700 and  ~600 km s^-1, respectively, indicating that the dust is still providing extinction at this epoch. This was also noted from the Hα profile by C97. The constancy of the line shifts compared to earlier epochs suggests that the dust resides in co-expanding clumps with optical depths much higher than unity also at this late epoch.

Assuming emission and absorption in a homogeneous sphere of velocity V_d, the dust radial optical depth τ_d can be found from the blue-shift of the peak ΔV (Lucy et al. 1989) (38)Using this equation for the Mg I] λ4571 and Hα line shifts with V_d = 2000 km s^-1 gives τ_d = 1.0−1.2, similar to its value at earlier epochs (Lucy et al. 1989; Wooden et al. 1993). Using ΔV = 400   km   s^-1, as found from the day 1862 and day 2210 spectra by Wang et al. (1996), gives a value of τ_d = 0.6, and using ΔV = 500   km   s^-1, as determined in Fassia et al. (2002), gives τ_d = 0.8.

Because it is unknown exactly in which zone(s) the dust resides, we model the dust as a constant, gray opacity within the core. From the estimates above, τ_d is likely to be in the 0.6−1.2 range. For our standard model we choose τ_d = 1.

We also note that dust in the ejecta can have strong effects on the cooling (KF98 a,b), with consequences for the spectrum depending on where the dust formed.

3. The observed spectrum

For the UV/optical range, we use the post-COSTAR HST Faint Object Spectrograph observations obtained on January 7 1995 (2875 days or 7.87 years after explosion), which covers the wavelength range from 1650 to 9200 Å. Details on the observations and the data reduction can be found in C97. The flux calibration is stated to be good to  ~5%. The aperture excludes the inner circumstellar ring, but the northern outer loop passes through the same line of sight as the central debris, and therefore contaminates the spectrum with some narrow lines. This is not a major problem because these lines can be easily identified. The strongest ones are C III] λ1909, [O II] λ3727, [O III] λλ4959, 5007, [N II] λλ6548, 6584, and Hαλ6563 (Panagia et al. 1996). We did not remove these narrow lines from the observed spectra displayed in the paper.

For the near-infrared, we use the AAT spectrum from Fassia et al. (2002) from day 2952. The IR spectrum contains contributions from both the ejecta and the inner equatorial ring. Because of the limited spectral resolution of λ/Δλ = 120, narrow lines from the ring can appear almost as broad as the ejecta lines. In addition, the removal of both continuum and line emission from star 3 and possibly star 2 adds a significant uncertainty. The accuracy of the absolute flux levels is stated as  ±40%. The NIR spectrum is therefore useful for line identifications, but only for approximate quantitative comparisons.

To compare the model to the observations, we adopt a distance of 51.4 kpc (Panagia 2005), the extinction curve (for λ < 8000   Å) derived from the neighboring star 2 in Scuderi et al. (1996), which has E_B − V = 0.19. For λ > 8000   Å we use the extinction curve in Cardelli et al. (1989). A different analysis by Fitzpatrick & Walborn (1990) find E_B − V = 0.16, and we take  ±0.03 to be a reasonable estimate for the extinction uncertainty. This dereddening gives a total luminosity in the 1600–8500 Å range of 3.4 × 10³⁵ erg s^-1. The value found by C97 is 3.5 × 10³⁵ erg s^-1, using a distance of 50 kpc and E_B − V = 0.2. We also adjust the model spectrum for interstellar UV absorption lines based on Blades et al. (1988).

The distance uncertainty to SN 1987A is stated as 1.2 kpc in Panagia (2005), which corresponds to a flux uncertainty of  ±4.8%. Dereddening the spectrum with E_B − V = 0.16 or 0.22 instead of 0.19 implies as change of  ~±13% in the total luminosity in the 2000–8000 Å range (which is the range we later use for estimating the ⁴⁴Ti mass). Combined, the observed luminosity then has an estimated uncertainty of  ~(5² + 4.8² + 13²)^1/2 =  ±15%.

4. Results

4.1. Energy deposition and degradation

At eight years, ⁴⁴Ti completely dominates the energy input over ⁵⁶Co and ⁵⁷Co. Using a ⁴⁴Ti mass of 1.5 × 10^-4   M_⊙, which we later show to be the best value for reproducing the observed spectrum, the total energy deposition in the ejecta is 1.9 × 10³⁶ erg s^-1. Of this,  ~85% is positron energy; the other  ~15% is gamma-ray energy, which represents  ~5% of the 5.2 × 10³⁶ erg s^-1 gamma-rays emitted. The other  ~95% of the gamma-rays escape the ejecta.

The dominance of positrons, combined with the assumption that these are locally absorbed, means that most of the energy (~80%) is deposited in the Fe/He clumps. Most of the gamma-rays are absorbed in the H zones, giving them  ~8% of the total energy deposition. The remaining  ~12% are split between the silicon, oxygen, and helium zones. All energy depositions are given in Table 5.

4.2. Temperature, ionization, and emissivities

The solutions for the temperature and ionization in the various zones are given in Table 5. From these steady-state solutions, we need to check the assumption that the recombination time-scale is shorter than the time-scale for density change n/ṅ = t/3 ~ 2.6 years. The last column in Table 2 lists the recombination times (1/[n_eα]) in years for the dominant element in the various zones. We see that a steady-state treatment of all core zones as well as of the helium envelope zone is justified, whereas the hydrogen envelope is experiencing freeze-out as expected.

The temperature is 70–170 K in all core zones (Table 5). The temperatures could be even lower if molecular cooling or dust cooling are important. This would not have any strong influence on the UV/optical/NIR spectra, however, because non-thermal processes dominate the output at those energies. Table 4 lists the dominant coolers in the various zones. For the envelope, we use the temperatures listed in Table 3.

We now discuss the individual zones in more detail to understand the emission lines that are produced. In Sect. 4.5 we discuss how the energy deposition in the various zones changes if there is a non-local positron deposition.

4.2.1. Fe/He zone

The high electron fraction in the Fe/He zone (~0.17) implies that a large part (~59%) of the deposited energy goes into heating, re-emerging mainly in the [Fe II] 26 μm line (Table 5). The reason that this line dominates over [Fe I] 24 μm is that it has a much larger collision strength,  ~6 versus  ~0.02 (Pelan & Berrington 1997; Zhang & Pradhan 1995). Because the Fe/He zone absorbs  ~80% of the deposited energy (assuming on-the-spot positron absorption), and  ~60% of this goes into heating,  ~50% of the total deposited energy or  ~9 × 10³⁵ erg s^-1 (or a factor two or so less considering the dust) should emerge in the [Fe II] 26 μm line at this epoch. We discuss observations of this line in Sect. 5.

About 25% of the energy goes into ionizations, with  ~14% going to He I,  ~7% to Fe I and  ~4% to Fe II. Iron reaches an ionization degree of  ~40%, while helium reaches  ~6%. Fe III is largely neutralized by charge transfer reactions, mainly with Fe I. Nickel, which makes up  ~4% of the zone mass, is almost fully ionized owing to charge transfer with Fe II. He I, Fe I, and Ni I therefore emit strong recombination lines. The recombination line spectrum of Fe I is rather smooth with significant emissivity in the 2000–6000 Å range; this constitutes a major part of the emitted flux in this range (Sect. 4.3). The He I recombination emission is mainly in the form of the two-photon continuum plus He I λ626. This flux is mostly continuum-absorbed by other elements and serves to additionally boost the ionization of iron and other metals (see also KF98 b). In the optical, the strongest recombination line is the 5876 Å line. Helium also emits IR lines at 1.083 μm and 2.058 μm. The 1.083 μm line is mainly excited by recombinations, while the 2.058 μm line arises mainly as a result of non-thermal excitations.

The remaining  ~16% of the deposited energy excites resonance lines from the ground states of Fe I (~9%), Fe II (~3%) and He I (~3%). For the Fe I transitions, we use the approximate Bethe rates, whereas for Fe II and He I we have more accurate cross sections.

Based on modeling of the Fe II lines during the first 800 days (Li et al. 1993), the density of the Fe/He zone may be lower than the one we use here, but not higher. Lowering the density of this zone by a factor of two increases the heating fraction by about five percentage units, at the expense of excitations and ionizations. Assuming that the UV/optical/NIR flux scales in proportion to the excitation plus ionization energy, the UV/optical/NIR emission from the Fe/He zone then falls by  ~12% by such a maneuver, and the total UV/optical/NIR emission by  ~5–10%. Using V_core = 2500 km s^-1 instead of 2000 km s^-1 (which lowers all densities by a factor of  ~2) lowers the total ejecta emission in the 2000–8000 Å range by  ~10%.

4.2.2. Si/S zone

The higher density and the smaller amount of ⁴⁴Ti (~5% of the total) present in this zone compared to the Fe/He zone implies a lower electron fraction (~0.02). A smaller part of the deposited energy therefore goes to heating (31%). Of the rest, 42% goes into ionizations, mainly of the abundant Si I and S I. Still, these elements stay almost fully neutral because of rapid charge transfer with elements such as Fe I, Mg I and Ti I. These recombine by further charge transfer reactions, and it is eventually calcium that is the dominantly ionized species. This means that quite a strong Ca I recombination spectrum is emitted from this zone. In Sect. 4.7, we show that this holds true even without charge transfer, because calcium is also strongly photoionized.

About 27% of the energy goes to non-thermal excitations, mainly for transitions in Si I and S I.

4.2.3. O/Si/S zone

This is the smallest oxygen zone (0.16 M_⊙) but contains some ⁴⁴Ti (~2% of the total). The largest part of the energy (56%) goes into ionizations, which through charge transfer reactions eventually leave iron and aluminium as the dominantly ionized elements. The non-thermal excitations occur in similar amounts for O I, Si I and Mg I.

4.2.4. O/Ne/Mg zone

The O/Ne/Mg-zone is by far the most massive oxygen zone (1.89 M_⊙), but contain no ⁴⁴Ti. The total energy deposition here is 2.8% of the total. About half of the energy goes into ionizations of oxygen. A series of charge transfer reactions neutralize the oxygen ions and transfer the ionization to sodium, which becomes the most abundant ion. This zone therefore emits strong Na I recombination lines, including the doublet at 5890, 5896 Å as well as at 5683, 5688 Å, 6154 Å and 8183, 8194 Å.

In our model, the dominant charge transfer reactant with O II is C I, where a reaction to the excited 2p(¹D) state in O I is close to resonance with an energy defect of only 0.4 eV. With the recipe for estimating rates given in Pequignot & Aldrovandi (1986), we estimate the reaction rate to be 1 × 10^-9 cm³ s^-1, which together with the relatively high carbon abundance in this zone (4.3% by number) makes it the dominant reaction for neutralizing the oxygen ions. Because 2p(¹D) is the upper state for the [O I] λλ6300, 6364 lines, this has the interesting effect of maximizing the doublet emission efficiency by producing a photon for each O I ionization, which is equivalent to  ~7% of the deposited energy. Only  ~10% of radiative recombinations pass through the 2p(¹D) state at these densities (giving  ~0.7% of the deposited energy), and since other competing charge transfer reactions likely do not occur to 2p(¹D), the conversion efficiency is even lower if they would dominate.

The collisional cross sections for the [O I] doublet lines are moderate and only  ~1.5% of the deposited energy goes to direct non-thermal excitation of 2p(¹D). This state also receives some indirect contribution by non-thermal excitations to higher levels, the most important channel being via the 3d(³S) state. Collisions to this state contribute  ~0.4% to 2p(¹D) by cascading. The total energy emitted in [O I] λλ6300, 6364 is then  ~8.9% (7%+1.5%+0.4%) of the deposited energy if the charge transfer reaction with C I dominates, and only  ~2.6% (0.7%+1.5%+0.4%) otherwise.

Magnesium has an abundance (by number) that is  ~25 times lower than that of oxygen. But because the Mg I] λ4571 transition has a large collisional cross-section, it receives about the same number of non-thermal excitations as the [O I] doublet. About 3.7% of the deposited energy goes to ionizing Mg I, but Mg II is neutralized by charge transfer with Ni I, which in turn ionizes Na I.

Sodium is five times less abundant still, but Na I has an even larger cross-section to its first excited state, which also receives  ~1.5% of the energy through non-thermal excitations. This contributes further to the Na I λλ5890, 5896 emission, which is, however, dominated by recombinations.

4.2.5. O/C zone

This is the intermediate-mass oxygen zone (0.58 M_⊙), which receives 0.9% of the total energy deposition. The low abundance of sodium in this zone, compared to the O/Ne/Mg zone, makes magnesium the dominantly ionized element. A strong Mg I] λ4571 recombination line is therefore emitted from this zone, although more magnesium resides in the O/Ne/Mg zone.

4.2.6. He zone

Because of mixing, the helium zone has both a core component and an envelope component (Sect. 2). These components behave similarly. He II is neutralized by charge transfer, mainly with Si I, and ionizations eventually end up in Na II, Mg II, Si II, and Fe II.

Non-thermal excitations go mainly to the He I 2p(¹P) state (~15%), resulting in the 2.058 μm line, with only small contributions from recombinations for this line. As long as the positrons are locally absorbed in the Fe/He zone, all helium emission lines have stronger contributions from the Fe/He zone than from the He zone. As an example, the fraction of the total energy deposited in the ejecta going into non-thermal excitations of 2p(¹P) is 0.16% from the helium zones, and 2.3% from the Fe/He zone, which means a 93% contribution to the 2.058 μm lines from the Fe/He zone. Similar relations hold for all other helium lines as well. In Sect. 5 we explore this point in more detail. Because of the strong helium emission from the Fe/He zone, the He zone is not expected to leave any particular imprint on the spectrum at this epoch.

4.2.7. H zone

As with helium, hydrogen is present both in the envelope and as a mixed-in core component (Sect. 2). The envelope zones are powered mainly by freeze-out recombinations. Here, hydrogen emits a low-temperature case B recombination spectrum. The freeze-out emission from the envelope is generally stronger than the emission from the steady-state core, and it provides the bulk of the hydrogen lines and continua in the model. The core component provides the peaks of the hydrogen lines, but also serves to emit ionizing radiation into the other core zones.

In the core H zone, steady-state still prevails and instantaneous gamma-ray deposition determines the output. Ionization of H I accounts for  ~36% of the deposited energy. Because we are close to case B, most recombinations pass through n = 2. About 28% of the deposited energy goes into non-thermal excitations of the H I 2p state. Both in the core and in the (inner) envelope, 2s and 2p are still populated according to their statistical weights, i.e. as 1:3. In the core component, the Lyα optical depth is  ~4 × 10⁸, which gives a ratio of two-photon to Lyα emissivity of  ~5. The two-photon continuum therefore dominates the decay of the n = 2 state. The two-photon continuum can affect the spectrum in many ways. It can boost the ionization of metals with low ionization treshholds, both in the H zones themselves and in the other zones. It can also provide an important input into optically thick lines and affect the UV/optical/NIR spectrum by fluorescence. In wavelength windows where no or few lines exist, the emission may escape and contribute to the UV continuum. In the envelope, the Lyα optical depth is lower and Lyα dominates over the two-photon continuum.

Hα from the core is powered in similar amounts by non-thermal excitations and recombinations (~4% of the deposited energy each). About 70% of recombinations lead to Hα emission in a Case B scenario at 100 K (Martin 1988).

4.3. The spectrum

The model spectrum, for M(⁴⁴Ti) = 1.5 × 10^-4   M_⊙ and τ_d = 1, is shown together with the observed spectrum in Figs. 3–5. These figures contain five panels. Panel (a) shows the model spectrum (black) overlaid on the observed spectrum (red). The observed spectrum was dereddened and corrected for redshift (Sect. 3). It was also smoothed below 2200 Å and in the 7500–8500 Å range. Panel (b) identifies the contribution by individual elements to the spectrum. These are the elements that emitted the photon if it escaped directly, and the element that produced the last scattering/fluorescence if it did not. Panel (c) shows the part of the flux that was processed by scattering and/or fluorescence in red, and the part that comes from direct emission in blue. Panel (d) shows in which zones the photons were last processed. This is the zone of emission for photons that escaped directly, and the zone of the last scattering/fluorescence for those which did not. Finally, panel (e) shows from which zones the photons were originally emitted (which is the same as in panel (d) for the photons that escaped directly).

We see that the model reproduces most of the features in the observed spectrum, which bolsters our confidence in the fundamental assumptions in the model and to the code. The main exceptions occur in the 3000–3500 Å range, which we discuss below. Apart from this range, all major observed lines are present in the model spectrum, and, equally important, the model does not produce any strong lines that are not observed.

4.3.1. Scattering and fluorescence

Although the line blocking effect in the photospheric phase has been extensively modeled before (e.g. Karp et al. 1977; Lucy 1987; Wagoner et al. 1991; Mazzali & Lucy 1993, an analysis of the actual fluorescence in the nebular phase is rarer. Li & McCray (1996) study the process at 300 and 800 days, using a He I two-photon continuum that scatters in a single-density hydrogen envelope. They conclude that the process transfers a significant fraction of the UV emission into a quasi-continuum in the optical and NIR up to at least 800 days. The importance of the effect was found to be sensitive to the adopted temperature, which determines the population of the lower levels, and thereby the number of optically thick lines.

That especially resonance lines can be optically thick even at eight years can be seen with a simple estimate. The Sobolev optical depth for a line from the ground state can be written as (39)where X is the number fraction of the element and n is the total number density, which is in the range of 10⁴−10⁷   cm^-3 in the core at this epoch. Obviously, even elements with low abundances can still be optically thick in their ground-state resonance lines.

By looking at the panels labeled (b) in Figs. 3–5, we see that scattering and fluorescence are indeed important for the spectral formation at eight years. The fraction of the emerging flux that has been processed by line transfer (red) is  ~60% in the UV (<4000 Å),  ~30% in the optical (4000–7500 Å) and  ~30% in the NIR (7500–21 000 Å). Running the model with the line transfer switched off resulted in a  ~30% increase in the UV flux, a  ~10% decrease in the optical flux, and a  ~30% decrease in the NIR flux.

Few lines are unaffected by the line transfer. Those with the highest fraction of direct emission in the flux are Mg I] λ4571, Hβ, Fe I] λ5012, [Fe I] λ5697, the [O I]/Fe I complex at  ~6300–6400 Å, Hα+Ca I λ6573, He I 1.083 μm, [Si I] 1.60, 1.64 μm, He I 2.058 μm, and a few others in the NIR part of the spectrum.

By inspecting the bottom-most panels in Figs. 3–5, we see that a substantial fraction of the emerging flux, especially in the UV, has its origin in the hydrogen envelope. This is mainly freeze-out energy that is emitted as Lyα and two-photon continuum, and is then processed by scattering and fluorescence. Below  ~2500 Å essentially the whole spectrum originates in the envelope.

An alternative way of illustrating the fluorescence process is to plot the escape wavelengths versus the emitted wavelengths, which we did in Fig. 6 for all photons that have been scattered. The transformation of high-energy photons to low-energy ones can clearly be seen here, and it is clear that most of the fluorescent output originates as UV emission (λ_in < 4000 Å). The low occupation numbers of excited states makes any inverse fluorescence, i.e. conversion to shorter wavelengths, unimportant at this epoch, and the line transfer produces a systematic reddening of the spectrum. The contribution by many input wavelengths to each output wavelength shows that the fluorescence process is insensitive to the input spectrum, as also found by Li & McCray (1996). Some of the more distinct input lines are Lyα and Fe I/Fe II emission lines at  ~2500 Å,  ~2900 Å and  ~3900 Å. Lyα powers the emission in Mg II 2795, 2802 (Sect. 4.3.2).

4.3.2. Line identifications

We now go through the spectrum, identifying and commenting on the most prominent lines. We use the same line peak wavelengths as in C97.

UV.

As panel (e) in Fig. 3 illustrates, most of the emerging far-UV flux (below  ~2400 Å) has its origin as emission from the envelope. Panel (b) shows that most of this flux has been reprocessed by scattering and fluorescence. Longward of  ~2400 Å the core begins to make significant contributions. The feature at 2400 Å is mainly scattering in the Fe II UV 2 multiplet, with similar contributions from the envelope and the iron zone. The feature is blended with an Fe I emission line on the red side. The feature at 2600 Å is also a blend of Fe II scattering (UV 1) and Fe I emission. The strong line at 2513 Å appears to be scattering in Si I UV 1, complemented by Fe I.

The feature at  ~2790 Å is scattering by Mg II λλ2795, 2802 in the envelope. We investigated the formation of this line and found it to be dominantly powered by Lyα, which is absorbed in the Mg II 1240 Å line, which then produces fluorescence at 2795, 2802 Å, 2798 Å, 2928, 2936 Å, 9218 Å, and 10 914 Å. There is some flux missing on the blue side of the Mg II feature, but we find no Fe II contributions that could improve this, as suggested in C97.

Next to the Mg II feature follows emission/scattering in the Mg I λ2852 line. This feature has a broad component that comes from scattering in the envelope (the line is optically thick out to  ~6000 km s^-1 even though magnesium is  ~99% ionized in the envelope), and a narrow component that is dominated by emission from the O/Ne/Mg zone. On the red side there is some contribution from iron and other elements.

The 3000–3500 Å range is not well reproduced by the model, and no definitive identifications can be made. Panels (d) and (e) in Fig. 3 show that most of the flux here originates in the Fe/He zone. The model produces strong lines at 2986 Å (Fe I), 3090 Å (Fe I and Al I), 3160 Å (Ti II scattering), 3270 Å (Fe II), 3310 Å (Ni I scattering), 3400 Å, 3440 Å and 3490 Å (emission and scattering mainly from iron, cobalt, and nickel). The overproduction of several lines in this region could conceivably be due to the omission of trace elements that could provide additional line absorptions. We surveyed the line lists in the 3000–3500 Å range of all omitted elements with Z < 30, however, and there are no lines from these elements that could be important scatterers at this epoch. The discrepancy could also be caused by the lack of accurate collision cross-sections for the non-thermal excitations in Fe I (we use the Bethe approximation). Although several of the lines are produced by scattering/fluorescence, the origin is often an Fe I emission line. As we see in Sect. 4.7, charge transfer can also be responsible by introducing uncertainties in the ionization balance. Finally, it is possible that the dust extinction is wavelength-dependent and should be higher in the UV. This would affect the UV lines coming from the core more than the ones coming from the envelope.

We identify the feature at 3650 Å as  ~50% Balmer continuum and  ~50% scattering by mainly Fe I and Cr I in the iron core. The “contamination” of the Balmer continuum makes the attempts to determine of the hydrogen zone temperature in W96 and C97 unsuitable, since they are based on an assumption of pure Balmer emission. The temperatures derived there are also significantly higher than in the ones calculated in our model (~100 K). Taking the contaminations of both the Balmer continuum and Hα into account (see below), the observed ratio between the Balmer continuum and Hα is 0.25, which is close to the theoretical ratio of 0.32 obtained by using the Case B recombination rates at 80 K in Martin (1988).

The strong observed line at 3730 Å is reproduced by the model and appears to be scattering and emission from Fe I. We find the [O II] λ3727 contribution to be negligible, agreeing with C97 that there is no need for this line, a possibility raised in W96.

Optical.

The 4000–4500 Å range is a blend of Fe I, Ca I, H I, and others. The strongest line in this range is at 4223 Å, which the model identifies as Ca I λ4226 emission from the Si/S zone. As we discussed in Sect. 4.2.2, Ca I is likely to emit most of the ionization energy in this zone, a conclusion which is supported by the presence of the 4226 Å line in the observed spectrum. This line is sensitive to density; Aβ has a value of  ~600 s^-1 for the density we use here (a number density ten times lower than in the O/Ne/Mg zone), which dominates the rate of A_tot = 140 s^-1 for other radiative branching paths. But if the Si/S zone would have the same density as the oxygen zones, Aβ would be only  ~10 s^-1 (calcium is also less ionized at higher density and so the scaling is not proportional), and the other deexcitation paths would then dominate. The presence of this line in the spectrum, at this epoch, therefore requires a number density in the Si/S zone at least  ~10 times lower than the one of the oxygen zones. Because the Si/S zone contains some (~5%) of the ⁵⁶Ni, according to the explosion model we use, this is consistent with the physical picture that radioactive decays cause expansions in the early phase of the supernova, as was found to have happened for the Fe/He clumps by Li et al. (1993).

The rest of the 4200–4500 Å “plateau” is mainly emission and scattering from the Fe/He zone, complemented with Hγ at 4343 Å. The Fe II contribution, which is suggested to dominate the features at 4223 Å, 4339 Å, and 4453 Å in C97, is only at the  ~10% level in our model.

Mg I] λ4571 is one of the lines with the smallest contribution from scattering/fluorescence. The model underproduces the flux by a factor of  ~2–3, possibly as a result of the uncertain ionization balance in the O/Ne/Mg zone (Sect. 4.2.4). Although most of the synthesized magnesium resides in this zone, its contribution to the Mg I] λ4571 line is moderate because Mg II is largely neutralized by charge transfer with nickel and sodium, leaving only non-thermal excitations to power the line. It is possible that we use rates that are too high for some of these charge transfer reactions, and that Mg II instead recombines radiatively, in which case the Mg I] λ4571 line would become stronger. Our model gives emission of similar magnitudes from the O/Ne/Mg zone, the O/C zone, and the He and H zones.

The feature at 4850 Å is  ~2/3 Hβ and  ~1/3 Fe I in the model. Removing the Fe I contamination gives a Balmer decrement of  ~5, which agrees well with a Case B recombination scenario at  ~100 K.

The 5000–5700 Å range is a quasi-continuum dominated by Fe I lines from the iron core, with about 3/4 of the flux being unscattered. The range is well reproduced by the model, apart from an overproduction of the Fe I line at 5007 Å (the narrow observed line does not come from the ejecta, but from the northern loop).

The line at 5900 Å is emission and scattering in Na I 5896 Å. The scattered part is predominantly from emission in Na I 5890 Å, but also from He I 5876 Å emitted in the Fe/He zone. Panel (e) in Fig. 4 shows that most of the energy originates in the O/Ne/Mg zone, which means that the helium contribution is minor. Na I emits strongly from the O/Ne/Mg zone because of a series of charge transfer reactions and a high photoionization rate that make it the dominant ion in this zone (Sect. 4.2.4). Apart from the 5890, 5896 Å doublet, the Na I recombination spectrum also includes lines at 5690 Å, 6154 Å, 8190 Å and 1.14 μm. There are indeed observed lines at  ~5690 and  ~6154 Å, but blending with Fe I and Ca I makes identifications uncertain. At 8190 Å the spectrum is too noisy to clearly detect a line, and in addition there is also contamination from other lines in the model. There is a clear detection at 1.14 μm in the day 2112 spectrum in Fassia et al. (2002) (not shown here), but all in all it is difficult to tell whether a significant Na I recombination spectrum is emitted by the ejecta or not. The model overproduces the Na I λλ5890, 5896 doublet, and it is possible that the ionizations instead go to magnesium in the O/Ne/Mg zone, especially because the Mg I] λ4571 line is underproduced. However, in Sect. 4.7 we show that the Na I doublet is overproduced also if the charge transfer reactions are switched off.

[O I] λλ6300, 6364 is contaminated, which the deviation of the observed ratio from the expected 3:1 value verifies. In the model, most of the contaminating flux is from Fe I. The complex is well reproduced by the model, although the good fit relies on the assumed high rate for the O II + C I  →  O I(2p¹D) + C II charge transfer reaction, as discussed in Sect. 4.2.4. If this reaction is not the dominant channel for O II recombinations, the [O I] λλ6300, 6364 emission would decrease by a factor of  ~4. The uncertainties regarding this reaction combined with the Fe I contamination make the oxygen doublet unsuitable for drawing any strong conclusions at this epoch, including attempts to determine the oxygen mass.

Interestingly, Hα is significantly contaminated by Ca I 6572 Å in the model. This line originates, like the Ca I 4226 Å line, mainly in the Si/S zone, but also in the O/Si/S zone. Being the first transition from the ground state in Ca I, the line has a high effective recombination rate and is about twice as strong as the Hα emission from the core H zone. Note that this component is needed to satisfactorily reproduce the “Hα” profile, which would otherwise be too flat-topped. For the high-velocity part of Hα, the good fit of the model is an important verification that our treatment of the freeze-out in the envelope is accurate.

The [Ca II] λλ7291,7323 doublet is well reproduced and results from pumping in the H and K lines. Figures 4d, e show that the fluorescence is produced by both synthesized and primordial calcium. The H, K lines are optically thick in all core zones, and in the envelope out to  ~7000 km s^-1. Figure 6 shows that most of the photons start as emission close to the H and K lines. The flux in the Ca II λλ8600 IR triplet is similar to the flux in the doublet, as it must be for a H, K pumping scheme.

Infrared.

In the infrared, some of the observed lines are dominated by narrow circumstellar components seen at low resolution. We identify the 1.083 μm line, observed in the day 2112 spectrum (Fassia et al. 2002), with He I emission from the Fe/He zone. The broad component has an observed flux (Table 2 in Fassia et al. 2002) similar to the model flux. Paγλ1.0938 μm also has a strong narrow component. The model shows the ejecta emission here to be dominated by [Si I] 1.099 μm. The model flux is 1.5 times the the observed broad component at 2112 days.

Further, we identify 1.14 μm with Na I, 1.20 μm with Si I, 1.25 μm with [Fe II], and 1.28 μm with Paβ. The strong 1.44 μm line is [Fe I] emission from the core, with a model flux 1.8 times the observed value. Also 1.50 μm and 1.53 μm are Fe I emission. 1.60 and 1.64 μm are dominantly [Si I] emission from the Si/S zone, and 2.058 μm is He I from the iron core. The 2.058 μm line is overproduced in the model. We also refer to Kjær et al. (2010), where we discuss the IR spectrum at an age of 19 years, with most line identification being the same as at this epoch.

Summary of the spectrum.

In general, we find good agreement with the line identifications in W96 and C97. The most important difference is that we find many lines to be from Fe I rather than Fe II. This is what to expect since Fe I emits not only a recombination line spectrum, but also receives more non-thermal excitations because Fe I atoms outnumber Fe II ions by a factor of  ~1.5 (see Sect. 4.2.1). In our model,  ~30% of the flux between 2000 and 8000 Å is from Fe I, and only  ~9% is from Fe II (including both emission and scattering). We do not, however, include ionizations to excited states in Fe II, which causes some underestimate of the Fe II emission.

4.4. The ⁴⁴Ti mass

The character of the spectrum, i.e. the line ratios and the ratio of lines to continuum, is not very sensitive to the absolute level of the radioactive powering. The character is partly determined by the relative amounts of energy that go into excitations and ionizations, which is insensitive to the total deposition. It is also influenced by the ionization balance of the dominant elements in each zone, which is likewise insensitive to the total deposition; the ionization fraction of the dominant element approximately scales as the square root of the deposition. Furthermore, there is a negative feedback present because a higher electron fraction reduces the fraction that goes into ionizations (and excitations) at the expense of the heating fraction.

The ⁴⁴Ti mass therefore mainly determines the absolute flux level of the spectrum, rather than its form. Unfortunately, it is difficult to distinguish the ⁴⁴Ti mass from the uncertain dust properties. We therefore calculated the integrated flux in the 2000–8000 Å range for different combinations of M(⁴⁴Ti) and the dust optical depth τ_d (Sect. 2.6). We did not include the NIR because of the uncertainties in the calibration of these observations (Sect. 3). Figure 7 shows the ⁴⁴Ti mass that reproduces the observed flux in the 2000–8000 Å range as function of τ_d, together with the masses that produce the correct flux plus/minus 10% and 20%. Note that the emerging flux is less than proportional to the ⁴⁴Ti mass. The main reason is that a significant fraction of the flux comes from freeze-out emission rather than from instantaneous ⁴⁴Ti deposition. Another reason is that a higher deposition increases the electron fraction, which increases the fraction of the flux emerging in the FIR.

A ⁴⁴Ti mass of 1.5 × 10^-4   M_⊙ is our best match for τ_d = 1, which is our best estimate for the dust extinction from the line profiles. This result is only valid for the assumption of local positron deposition. The dust optical depth has a  ±30% uncertainty (Sect. 2.6), which according to Fig. 7 gives a  ±20% uncertainty for the best fitting ⁴⁴Ti mass. In Sect. 4.2.1 we noted that the uncertainty in the density of the Fe/He zone adds an uncertainty of  ±10% for the model flux, and in Sect. 3 we saw that the observed flux level has a  ±15% uncertainty from errors in distance, reddening, and flux calibration. Together, these give a (15%² + 10%²)^1/2 = 18% uncertainty for the ratio between modeled and observed flux. According to the contours in Fig. 7, this corresponds to a  ±30% uncertainty in the ⁴⁴Ti mass. Together with the dust error, the total uncertainty for the ⁴⁴Ti-mass becomes  ~(20%² + 30%²)^1/2 = 36%. The ⁴⁴Ti mass is then constrained to (1.0−2.0) × 10^-4 M_⊙.

4.5. Positron trapping

An described in Sect. 2.2, an interesting and important question is whether the positrons are indeed locally absorbed, as assumed so far, or if they leak into the other zones in the core before they are absorbed. We recomputed the spectrum for a model where we assume that the positrons are not absorbed on the spot, but in proportion to the electron content (free and bound) of each zone. This treatment is based on the Bethe stopping formula (Eq. (1)), where we neglect the small differences due to different effective ionization potentials. The number of electrons per zone is proportional to , where M is the zone mass and and are the average nuclear charge and atomic weight, respectively.

For this alternative model, the ⁴⁴Ti mass needed to give the correct flux in the 2000–8000 Å range for τ_d = 1 is only 8 × 10^-5   M_⊙, which we therefore used. The reason for the lower ⁴⁴Ti mass is that a smaller fraction of the positron energy is now emitted by far-IR cooling lines (mainly [Fe II] 26 μm) by the low-density Fe/He zone, and that the oxygen and hydrogen zones have some strong lines in the 2000–8000 Å range. Consequently, less ⁴⁴Ti is needed to get the same UV/optical/NIR luminosity.

In this model, the zones obtain total energy depositions (positrons + gamma-rays) given in Table 6, where we also give the relative change in the deposition fraction with respect to the on-the-spot model. From this, we expect positron leakage to lead to a severe weakening of emission lines from the Fe/He zone, and a boost by up to a factor of  ~10 for emission lines from the O/Ne/Mg, O/C, core He, and core H zones. The ⁴⁴Ti content of the silicon zones (Si/S and O/Si/S) however, is such that the deposition here does not change much between the two scenarios.

Figure 8 shows the resulting model spectrum compared to the standard on-the-spot model (which is the same as in Figs. 3–5). On the whole, the spectra are quite similar, but contain some important differences as well.

The hydrogen emission from the core is boosted by a factor of several, as expected. Hα and the Balmer continuum are now overproduced. However, because the fraction of the hydrogen that we mix into the core (25%) is not tightly constrained, there is some room to improve the H lines by reducing this in-mixing. That would, however, increase the deposition in the other core zones and strengthen the lines from them.

The [O I] λλ6300, 6364 doublet becomes much too strong, but this result is sensitive to the uncertain O II + C I  →  O I (2p¹D) + C II charge transfer rate used, as discussed in Sect. 4.2.4. If this reaction is switched off, the agreement is actually better in this model than in the on-the-spot model. The Mg I] λ4571 line is mainly emitted from the oxygen zones, and becomes about as much overproduced as it was underproduced before. The Mg I λ2852 line, however, now becomes much too strong (it is truncated in the figure). This line is powered by non-thermal excitations, and is therefore independent of the uncertain ionization balance in the O/Ne/Mg zone (magnesium is dominantly neutral in all scenarios). This line should therefore be a more reliable indicator than the recombination-driven 4571 Å line, and all this suggests that too much energy is now being deposited in the oxygen zones. The Na I recombination lines at 5890, 5896 Å, and 6160 Å from the O/Ne/Mg zone become overproduced as well. The 5890, 5896 Å doublet, which was already too strong without leakage, is now a factor  ~4 too strong.

For the rest of the spectrum, the two models produce surprisingly similar spectra. Because fluorescence makes important contributions at almost all wavelengths, turning off the iron emission lines from the Fe/He zone typically lowers the line fluxes by no more than a factor of about two. An example of this situation is the 5000–5500 Å range. The iron lines with the smallest contamination by fluorescence in the UV/optical are the lines at 3270 Å, 3490 Å, 5012 Å, and 5159 Å, and 5697 Å. All of these lines seem somewhat better reproduced in the leakage model.

Because the flux calibration in the NIR is uncertain (Sect. 3), we refrain from comparing the models in this range. All in all, we think the on-the-spot model produces a spectrum that better agrees with the observations, although the uncertainties mentioned above prohibit a definitive ruling out of the leakage scenario. The average χ²-value (per spectral bin) is about twice as high in the leakage model as in the on-the-spot model.

4.6. Fragmentation level

The fragmentation level, N_cl, is a proxy for the hydrodynamical mixing. This parameter determines the length of the paths that the photons travel in their parent zone before they have the chance to enter the other core zones. A high fragmentation level, implying quick exposure to other zones, will lead to a higher degree of radiative energy exchange between the zones. If one were able to constrain this parameter it would be very useful for evaluating multi-dimensional explosion models.

To check the sensitivity of the spectrum to this parameter, we compared spectra for models using N_cl = 10 and N_cl = 10⁵. The results were very similar, and we therefore conclude that the fragmentation parameter has no significant influence on the spectrum at this epoch. One possible explanation for this is that few line optical depths can be expected to be close enough to unity that an abundance change (by switching zones) between a high value in the zones where the element is synthesized to the (low) primordial values where it is not, will make any difference. Another could be that few of the line absorptions occur before the photon has exited its parent zone, even in the low fragmentation limit. Finally, a significant part of the line transfer occurs in the envelope, where there is no mixing in the model.

4.7. Effects of charge transfer

Finally, the uncertainty in many of the charge transfer rates prompted us to examine the influence of these reactions on the spectrum. It is also interesting to see how important these reactions are for supernova spectra in general. We recomputed our standard model (M(⁴⁴Ti) = 1.5 × 10^-4   M_⊙, τ_d = 1, and on-the-spot positron deposition) with all charge transfer reactions switched off. The resulting spectrum is compared to the standard one in Fig. 9.

The total flux in the 2000–8000 Å range changes by only a few percent, so charge transfer should not affect our conclusions regarding the ⁴⁴Ti mass. Several lines differ significantly between the two models though. The [O I] λλ6300, 6364 doublet is weakened for reasons discussed in Sect. 4.2.4, and the Fe I lines now dominate the doublet. At the same time, the O I λ7775 recombination line, which is not observed, emerges in the model. This suggests that the oxygen ions are indeed neutralized by charge transfer reactions, as in our standard model.

The important Ca I lines (4226 Å and 6572 Å) do not decrease by much, which shows that charge transfer is not critical for them. Ca I is significantly ionized by the internal radiation field, and so emits strong recombination lines even without charge transfer.

The Na I λ5896 emission is reduced somewhat, but is still too strong. Just as with Ca I, Na I is significantly ionized by the internal radiation field, and consequently strong recombination emission still occurs. An important conclusion is therefore that the emission lines from elements with low ionization thresholds are is not very sensitive to the charge transfer, because low ionization thresholds also correlate with high radiative ionization rates.

For many lines, it is not possible to directly say why they change. The “picket fence” of absorption lines changes as the ionization balance changes, which opens some wavelength windows and closes others. Consider, for example, the Fe I emission line at  ~2990 Å. This emission line is still there in the model without charge transfer, but now a Ni I absorption line has become optically thick that absorbs it. The UV spectrum is actually somewhat better reproduced in the model without charge transfer, which may indicate that we overestimate some of our adopted rates. Unknown charge transfer rates clearly constitute one of the main obstacles to accurate supernova spectral modeling today, and more calculations and measurements of these would be highly desirable.

5. Discussion

Although most of the optical and infrared emission escapes freely in the nebular phase of supernovae, the UV emission does not. The freeze-out energy in the envelope, as well as non-thermal excitations and ionizations in the core, produce strong UV emission that is significantly absorbed by atomic lines and emerges by fluorescence at longer wavelengths. By this process, the emerging flux in the optical and NIR is enhanced by  ~10% and  ~30% compared to a purely nebular spectrum at the epoch we have looked at here. While the non-local line transfer can have a large impact on individual lines at late times, the fact that it has a moderate impact on the supernova photometry is important for judging the accuracy of models without this effect, which has previously been difficult.

While the fraction of UV photons that are absorbed decreases with time, the fraction of the optical/NIR spectrum that comes from fluorescence may actually increase (Li & McCray 1996). The reason for this is that with time the cooling lines move from the optical/NIR into the FIR, and non-thermal processes do not produce very strong optical/NIR emission. Our finding that fluorescence is important also at eight years illustrates this point well.

We showed that the only FUV flux that escapes comes from far out in the hydrogen envelope. A consequence of this is that the ejecta appear larger in the UV than in the optical (Jakobsen et al. 1994; Wang et al. 1996). At longer UV wavelengths both the envelope and the core contribute to the spectrum.

From a mathematical perspective, it is not a very well-conditioned problem to attempt exact modeling of the UV lines, because the output depends sensitively on the exact shape of the “picket fence” of absorption lines, which in turn depends on many uncertain parameters such as abundances, densities, and charge transfer rates. Considering this complexity, it is probably not realistic or meaningful to attempt any modeling of the UV lines in much more detail than we have done here.

As discussed in Sect. 2.2, the positrons emitted in the ⁴⁴Ti decay will enter the neighboring oxygen, helium, and hydrogen zones at late times unless a magnetic field exists to trap them. This would have significant consequences for the FIR cooling lines. The model flux in the [Fe II] 26 μm line, discussed in Sect. 4.2.1, corresponds to a peak flux of  ~1.8 Jy at day 2875. ISO observations at day 3425 failed to detect this line, giving a upper limit for the peak flux of 0.64 Jy (Lundqvist et al. 2001). Spitzer observations at day 6190 (Bouchet et al. 2006) show a distinct broad feature (FWHM ~ 3800 km s^-1) at 26 μm, but with a peak flux of only  ~0.02 Jy. The slow decay of ⁴⁴Ti should make the situation at 6190 days similar to the one at 2875 days, and it is therefore surprising that the observed flux in this line is about two orders of magnitude weaker than expected. Possible explanations are that the positrons are indeed leaking out into the other core zones, or that dust is cooling the Fe/He zone. In the positron leakage model, the flux in the 26 μm line is  ~0.16 Jy at 2875 days (for a ⁴⁴Ti mass of 8 × 10^-5 M_⊙ and no dust correction). As we found in Sect. 4.5, a positron leakage scenario produces a UV/optical spectrum with several lines that are clearly too strong. Our favored explanation for the weak [Fe II] 26 μm line is therefore that dust cooling occurs in the Fe/He zone. The possibilities of having dust cooling in the ejecta was discussed in KF98 a and in Lundqvist et al. (2001).

In the case of magnetic confinement of the positrons, all helium lines are dominated by helium in the Fe/He zone, which is produced by freeze-out in the explosion. Although we do not identify any distinct observed lines with helium in the eight-year spectrum, the 19-year spectrum in (Kjær et al. 2010) contains a strong and distinct 2.058 μm line, which, as we speculated there, probably shows that α-rich freeze-out has indeed happened in the Fe/He zone; without it the local positron deposition could not produce such a strong line.

Together with the ⁵⁶Ni mass of 0.069     M_⊙ (Bouchet et al. 1991) and the ⁵⁷Ni mass of  ~3.3 × 10^-3   M_⊙ (Fransson & Kozma 1993; Varani et al. 1990; Kurfess et al. 1992), the determination of the ⁴⁴Ti mass should put considerable constraints on supernova explosion models. Our result of   agrees with the result of (1−2) × 10^-4 M_⊙ found in C97, although that result was based on an analysis of [Fe II] lines, whose identifications are different in our model on several instances, and on the [O I] λλ6300, 6364 doublet, whose emission is highly sensitive to uncertain charge transfer reactions, and which is also contaminated by Fe I lines. It is therefore a bit of a coincidence that we get similar results for the ⁴⁴Ti mass. Fransson & Kozma (2002) determine the ⁴⁴Ti mass to (0.5−2) × 10^-4 M_⊙ by comparing the observed B and V band HST photometry between 1800 and 3500 days to the output from a time-dependent model without non-local radiative transfer (as in KF98 a,b). We showed here that this transfer, although it is important for individual lines, only boosts the output in the optical/NIR range by 10–30% at 2875 days (Sect. 4.3.1), which implies that the results in FK02 should have a corresponding accuracy. An inspection of Fig. 4 in the FK02 paper shows that the best-fitting ⁴⁴Ti mass is in the (1.5−2) × 10^-4 M_⊙ range, which agrees well with our results here; ignoring non-local radiative transfer should give a somewhat over-estimated ⁴⁴Ti mass, but not by much. Our result also agrees with the range of (0.8−2.3) × 10^-4 M_⊙ found by Motizuki & Kumagai (2004). Their result is, however, is based on a comparison betwen the deposited energy and the bolometric luminosity at 3600 days, the latter being highly uncertain since most of the flux comes out in the unobservable FIR.

The ⁴⁴Ti isotope is a result of the α-rich freeze-out (e.g. Woosley et al. 1973), which is a consequence of the slow rate for the triple-alpha reaction, which cannot keep the ⁴He abundance at the NSE value as the temperature and density decreases. The high ⁴He abundance in turn results in a higher abundance of ⁴⁴Ti than is the case in NSE. The ⁴⁴Ti/⁵⁶Ni ratio therefore sensitively depends on the entropy (e.g. Woosley & Hoffman 1991). High temperature and/or low density increases this ratio. The explosion model we use (from WH07) has a mass ratio ⁴⁴Ti/⁵⁶Ni = 3.1 × 10^-4, which gives a ⁴⁴Ti mass of only 2.1 × 10^-5 M_⊙ (for M(⁵⁶Ni) = 0.069 M_⊙), which is clearly ruled out according to our modeling. For a spherically symmetric explosion of a 20 M_⊙ star, WW95 find an ejected mass of (1−2) × 10^-5 M_⊙. On the other hand, Thielemann et al. (1996) find a much larger mass of 1.7 × 10^-4 M_⊙, closer to the value we derive here. The differences are partly caused by the different nuclear reaction rates used, but mainly because WW95 initiate the explosion by injecting momentum (piston-driven), while Thielemann et al. (1996) deposit thermal energy (Timmes et al. 1996). The solar abundance of the decay product ⁴⁴Ca requires a ⁴⁴Ti-production higher than in WW95 by a factor  ~3, averaged over all progenitor masses, suggesting that the yields in those models are indeed low (Timmes et al. 1996).

Nagataki et al. (1998), Nagataki (2000) and Maeda & Nomoto (2003) have investigated the effects of asymmetries in the explosion. These (simplified) models have considerably higher temperatures and entropies in the polar direction. Consequently, the ⁴⁴Ti production occurs mainly in this direction, and the total ⁴⁴Ti mass can be considerably higher than in a spherically symmetric model.

Further, the rates of several of the nuclear reactions involved in the ⁴⁴Ti production are uncertain. Nassar et al. (2006) find an increased rate for the important ⁴⁰Ca(α,γ)⁴⁴Ti reaction, resulting in a factor of two higher ⁴⁴Ti abundance. The theoretical range of ⁴⁴Ti mass is therefore highly uncertain, as discussed by The et al. (2006) for Cas A.

In addition to ⁵⁶Ni, ⁵⁷Ni and ⁴⁴Ti, other radioactive isotopes may be present in the ejecta. In particular, ⁶⁰Co and ²²Na may have sizeable abundances. These isotopes are produced by neutron and proton capture, respectively, during carbon burning, and have uncertain yields that are sensitive to the convective prescription (Timmes et al. 1996).

The ⁶⁰Co-isotope has a decay time of 7.605 years, and may have a similar mass to ⁴⁴Ti (computed as 2.5 × 10^-5 M_⊙ in WW95 and 6.5 × 10^-5 M_⊙ in WH07). The gamma-ray luminosity for these two masses are 56% and 150% of the ⁴⁴Ti gamma-ray luminosity, respectively, and the electron deposition is 12% and 30% of the ⁴⁴Ti positron deposition, respectively.  ⁶⁰Co may therefore make a non-neglegible contribution to the powering at eight years. The β-decays occur as electrons with energies up to 317 keV, with an average energy of 96.5 keV. This is significantly lower than the energy of ⁴⁴Ti-positrons, and the stopping range is therefore much shorter. A major difference is that most of the ⁶⁰Co electron energy should be deposited in the oxygen-rich zones. These account for approximately half of the total column density of the core, and from Fig. 2 it can be seen that they can easily trap most of these electrons. The increased deposition in the oxygen-rich regions means that an appreciable ⁶⁰Co mass may mimic the effects of the positron leakage in the ⁴⁴Ti dominated models (Sect. 4.5). It will in particular increase the O I, Mg I, and Na I line fluxes noticeably, further complicating the question of where the positrons are deposited.

²²Na has a decay time of 3.75 years, emitting a gamma-ray at 1.274 MeV and in 90% of the cases a positron, while 10% of the cases occur by electron capture. The positrons have an average energy of 215.5 keV and a maximum energy of 545.4 keV. The ²²Na mass in the explosion model is 4.5 × 10^-6 M_⊙, in Woosley et al. (1989) it is 2.0 × 10^-6, in WW95 it is 3.0 × 10^-7 M_⊙ and in Thielemann et al. (1996) it is 1.3 × 10^-7 M_⊙. This gives gamma-ray luminosities of 0.5–17% of the ⁴⁴Ti gamma-ray luminosity, and positron depositions of 0.2–5.8% of the ⁴⁴Ti positron deposition. Based on these models it is therefore unlikely that ²²Na is important for the energy budget at the epoch studied here.

It is, unfortunately, difficult to constrain the contribution from ⁶⁰Co and ²²Na to the eight-year-spectrum because of the uncertainties in the modeling of the oxygen zones, and the possibility of non-local positron/electron absorption occurring. We assume both of these isotopes to make neglegible contributions to the powering at eight years. The shorter decay time scales of ⁶⁰Co and ²²Na, compared to ⁴⁴Ti, means that their influence varies with time. The modeling of the spectrum and light curve at other epochs should therefore have the potential to distinguish between the ⁴⁴Ti and the ⁶⁰Co/²²Na input. Unfortunately, only a few spectra exist at epochs later than the one studied here, and they are increasingly contaminated by the flux from the ring collisions. Modeling of the NIR spectrum at 19 years has shown that the combination of M(⁴⁴Ti) = 1 × 10^-4 M_⊙  and τ_d = 1 matches the observed NIR spectrum at this epoch reasonably well (Kjær et al. 2010). However, this spectrum may be affected both by scattered light and X-ray input from the circumstellar ring (Larsson et al. 2011).

Small or moderate masses of ⁶⁰Co and ²²Na are supported by the modeling of the SN 1987A photometry (FK02), which shows that the light curves do not deviate significantly from a pure ⁴⁴Ti input between 1800 and 3500 days.

6. Conclusions

Our main conclusions are:

• A spectral model based on the determination of the positron and gamma-ray degradation, the ionization, excitation and temperature structure of the ejecta, and the multi-line radiative transfer, reproduces most of the observed lines in the eight-year UV/optical/NIR spectrum of SN 1987A to good accuracy. We showed that many regions in the spectrum are produced by Fe I lines. The contribution by Fe II lines is small.

• Even many years after explosion, scattering and fluorescence are important processes for the spectral formation in type-II supernovae. At an age of eight years,  ~30% of the emerging optical/NIR flux and  ~60% of the UV flux in SN 1987A is produced by fluorescence. Non-local line transfer decreases the UV flux by  ~30%, increases the optical flux by  ~10%, and increases the NIR flux by  ~30%, compared to a nebular treatment.

• Most of the flux in the hydrogen lines and continua, as well as significant parts of the UV spectrum, are at eight years produced by freeze-out emission from the envelope rather than by instantaneous ⁴⁴Ti deposition. In the far-UV (≲2500   Å), essentially the whole spectrum is produced in the envelope.

• A scenario of positron leakage from the Fe/He clumps produces too strong emission in several lines from the hydrogen and oxygen zones. We therefore support the conclusion in Chugai et al. (1997) that the positrons remain trapped by a magnetic field.

• The [Fe II] 26 μm line in Spitzer observations at day 6190 is much weaker than expected. Possible explanations are dust cooling or that positron leakage is occuring at that epoch.

• Modeling the dust as a gray absorber with radial optical depth τ_d = 1, which is our best estimate from the line peak blue-shifts, the amount of ⁴⁴Ti synthesized in the explosion that best reproduces the total flux between 2000 and 8000 Å is 1.5 × 10^-4   M_⊙. We estimate the uncertainty in this result to  ±0.5 × 10^-4 M_⊙.

• One of our line identifications is the Ca I λ4226 line. If this identification is correct, the density in the Si/S zone must be at least  ~10 times lower than in the oxygen zones for the line to be as strong as observed at this epoch. Because the Si/S zone contains some of the ⁵⁶Ni, this adds further evidence to the picture that the early radioactive decays cause expansion in the ejecta.

---

Acknowledgments

We thank Alexander Heger for providing us with the explosion model, Peter Meikle for providing the NIR data, Cathy Ramsbottom for providing Fe II collisional cross sections, Sultana Nahar for consultations on radiative oxygen recombination, and Claes-Ingvar Bjornsson for useful discussions. We especially thank Mattias Ergon for comments on the manuscript and many hours of useful discussions. This work is supported by the Swedish Research Council and the Swedish National Space Board.

References

Appendix A: Atomic data

Table A.1 lists a summary of the model atoms used and the sources for the atomic data. The number of levels listed are the number of levels (including fine-structure) whose populations are solved for, followed by the total number of levels included for making the line list. For the population solutions, a size of 50 levels for the neutrals and 10 levels for the ions were typically used, with the larger number for the neutrals because these emit recombination spectra, whereas no ion does in any significance. No lines with optical depths higher than 10^-3 have their lower levels above these limits at this epoch. Elements that have non-thermal excitations or specific recombination rates to higher levels, have more levels included in the population solution. Two examples are Fe I and Fe II, where the populations of all  ~500 levels are calculated.

All atoms are solved for with fine-structure, which is especially appropriate at late times. For allowed transitions that lack published collision strengths, we use the van Regemorter formula (van Regemorter 1962). For forbidden transitions we used the van Regemorter formula with an absorption oscillator strength of 1 × 10^-2   g_u. This treatment is based on the fact that forbidden transitions usually have collision strengths about an order of magnitude smaller than allowed transitions. This treatment also allows a reproduction of the typically very large collision strengths for transitions of the type nl → nl′, since the van Regemorter formula contains a ΔE^-1 dependence. We checked the sensitivity to this treatment by comparing the spectrum with one using a constant Υ_u,l = 0.1g_lg_u for the forbidden lines. Only the Mg I] λ4571 line was sensitive to the treatment. The fine-structure collision strengths in this atom is on the order of 10⁵ (e.g. Sundqvist et al. 2008), and the first treatment is better.

If specific recombination rates are lacking, we allocate 30% of the total recombination rate to the ground multiplet and the remaining 70% in proportion to the statistical weights. If some specific recombination rates are available but not all, we allocate the difference between the total rate and the sum of the known rates to the states with unknown rates in proportion to their statistical weights. If all specific recombination rates are known, but do not match the total rate, we scale them so that they do. We do not include dielectronic recombinations, which are unimportant at the low temperatures here.

All ground-state photo-ionization cross sections are from a routine based on Verner et al. (1996). In addition, cross sections for some selected excited states are included (n = 2 in H I, 2¹S and 2¹P⁰ in He I, ¹D, ¹S and ⁵S in C I, ¹D and ¹S in O I, ³P in Mg I, ²P in Mg II, ¹D in Si I, ¹D in S I, ²D in Ca II, a⁵ F in Fe I and a⁴F in Fe II), taken from TOPBASE¹.

For information on the cross sections used in the Spencer-Fano routine, see KF92 and K98 a,b. An important update is the addition of detailed cross sections for Fe II (Ramsbottom et al. 2005, 2007), and for Ca I (Samson & Berrington 2001).

All Tables

All Figures