THE EVOLUTION OF DUST MASS IN THE EJECTA OF SN 1987A

First evidence for dust emission from the ejecta of SN 1987A was provided by high-resolution mid-IR images obtained with the Thermal Region Camera and Spectrograph (T-ReCS) on the Gemini South 8 m telescope (Bouchet et al. 2004). Spatially resolved images of SN 1987A at 450 and 870 μm obtained with the Atacama Large Milimmeter/Submillimeter Array (ALMA) between days ∼9200 and ∼9300 after the explosion provided further proof of thermal infrared (IR) emission from the supernova (SN) ejecta (Indebetouw et al. 2014; Zanardo et al. 2014). Combined with almost contemporaneous data obtained on day 9090 after the explosion with the Herschel satellite, the IR observations have been interpreted as evidence for the presence of a large amount of cold dust (Matsuura et al. 2011, 2015). If the emission arises from only silicates in the form of Mg₂SiO₄ then the observations require 2.4 ${M}_{\odot }$ of dust, that is, the precipitation of, respectively, 0.48 and 0.8 ${M}_{\odot }$, of Si and Mg from the gas phase of the ejecta. This mass greatly exceed the SN yield of either Si or Mg (∼0.1 ${M}_{\odot }$, see Table 1 and references therein). The mass of Si required to fit the IR spectrum can be reduced by the inclusion of carbon dust. A combination of silicate and amorphous carbon (AC) reduces the mass of silicates to 0.5 ${M}_{\odot }$. However, it requires an excessive amount of carbon, about (0.3–0.5) ${M}_{\odot }$ or over two to three times the yield of C, to be locked up in dust.

The mass of dust inferred from the late time (t ≳ 8500 day) IR emission also greatly exceeds the ∼10⁻³ ${M}_{\odot }$ of dust derived from early observations of SN 1987A on days 615, 775, and 1144 (Moseley et al. 1989; Dwek et al. 1992; Wooden et al. 1993). This difference has been interpreted as evidence for the continuous growth by accretion of refractory elements onto the seed nuclei that were initially formed in the ejecta (Wesson et al. 2015, and references therein). The ejecta temperature after day ∼1200 has dropped below ∼500 K (Sarangi & Cherchneff 2013). An accretion scenario for the evolution of the dust will therefore require most of the growth to take place by cold accretion, so that the final structure of the dust grains will primarily consist of a mantle, loosely bound with surface energies ≲1 eV, surrounding a small thermally condensed core. Such dust grains will not survive the reverse shock expected to propagate into the ejecta (Nozawa et al. 2007), contrary to the observed ∼0.04 ${M}_{\odot }$ of dust that was subjected to and survived the reverse shock in Cas A, or in the ∼10,000 old Sgr east remnant (Rho et al. 2008; Barlow et al. 2010; Arendt et al. 2014; Lau et al. 2015).

Also, the absence of the 9.7 and 18 μm emission features in the IR spectra on days 615 and 775 has been interpreted as evidence that the newly formed dust consisted of mostly carbon dust (Wooden et al. 1993; Ercolano et al. 2007a, 2007b; Wesson et al. 2015). Using the observed spectra of days 615 and 775 and a fully 3D radiative transfer model, Wesson et al. (2015) concluded that about (4–8) × 10⁻⁴ ${M}_{\odot }$ of carbon condensed in the ejecta of SN 1987A by those epochs. In particular, they derived an upper limit of ∼10⁻⁴ ${M}_{\odot }$ on the mass of silicate dust that could have formed in the ejecta. This low yield is puzzling, considering the preponderance of Si and Mg in the ejecta, and that SNe should be the major source of thermally condensed dust in the local and early universe (Dwek 1998; Todini & Ferrara 2001; Nozawa et al. 2003; Schneider et al. 2004; Dwek et al. 2007, 2014; Zhukovska et al. 2008; Calura et al. 2010, 2014; Cherchneff & Dwek 2010; Dwek & Cherchneff 2011; Gall et al. 2011a, 2011b). Finally, the observed decrease in the [Mg I]0.4571 μm and [Si I]1.65 μm emission lines around day 530 coincided with the first appearance of the IR excess in the SN light curve (Danziger et al. 1991, Lucy et al. 1991). If this decrease is caused by the efficient precipitation of Mg or Si into silicate grains (Lucy et al. 1991, Bouchet et al. 1991), then one would expect a silicate mass of 0.4 ${M}_{\odot }$ for a Mg and Si yield of 0.1 ${M}_{\odot }$ and a silicate composition of MgSiO3.

In this paper, we present an alternative model for the evolution of dust in the ejecta of SN 1987A and show that the currently observed dust mass in the ejecta on days ∼8500–9100 could already have been present at the earlier epochs. The conclusions of Ercolano et al. (2007a) and Wesson et al. (2015) represent just one possible, albeit troublesome, scenario for the evolution of the dust mass in SN 1987A. Large amounts of dust could have been hidden in optically thick regions of the ejecta, a possibility first raised by Lucy et al. (1989, 1991). We emphasize that our scenario for the evolution of dust mass in SN 1987A is not unique, but a viable alternative, which addresses the difficulties in previous models: (1) the large mass of either Si, Mg, or C that needs to be locked up in dust to account for the ALMA observations, (2) the growth of dust mass by cold accretion, creating loosely bound mantles that will not survive the passage of the reverse shock, and (3) the absence of any significant amount of silicate dust in the ejecta.

Our alternative model has important implications for the composition of dust in the early universe (z ≳ 7), before carbon-rich AGB stars left the main sequence. Even if only a fraction ≳10% of the silicate dust survives the injection into the ISM, the observations of SN 1987A firmly establish the role of core collapse SNe as the major source of thermally condensed silicate dust in the early universe. Interstellar dust is then silicate-rich, with important observational consequences for galactic extinction corrections and photometric redshift determinations (Dwek et al. 2014).

The development of an IR excess in the light curve of SN 1987A has been observed with ground-based photometry (Suntzeff & Bouchet 1990; Bouchet et al. 1991) and airborne spectrophotometric observations with the Kuiper Airborne Observatory (KAO; Wooden et al. 1993). Of these observations, we selected to perform detailed fits to the KAO data from days 615 and 775, mainly to compare our results with the model of Wesson et al. (2015) for the same epochs

The paper is organized as follows. We first present the basic equations for the escape of IR photons from a homogeneous dusty sphere, the dust model used in our fitting procedure (Section 2). In Section 3, we present the fits to the observations and the evolution of dust mass. We first start with the analysis of the late-time emission at epochs 8515 and 9090, when the ejecta was optically thin, and construct a model in which the IR emission arises from elongated composite grains consisting of a silicate MgSiO₃ matrix with AC inclusions. Such a configuration can arise as a results of the coagulation of silicates and AC grains at late times in the ejecta. The mass of Mg, Si, and C needed to be locked up in the dust is entirely consistent with nucleosynthesis constraints. We then show that the same mass of refractory elements could have already condensed out of the gas at the early epochs (days 615, 775, and 1144). Here the calculations assumed that the ejecta is homogeneous. A clumpy ejecta model is described in Section 4. In Section 5, we calculate the evolution of the IR luminosity and compare it to the energy input by radioactivity. A summary and the astrophysical implications of our results are presented in Section 6.

The general equation for the flux emerging from a homogeneous dusty sphere of radius R at distance D is given by

$$\begin{eqnarray}&&{S}_{\nu }(\lambda )=\displaystyle \frac{4{M}_{{\rm{d}}}\ \kappa (\lambda )}{4\pi {D}^{2}}\ \pi {B}_{\nu }(\lambda ,{T}_{{\rm{d}}})\ {P}_{\mathrm{esc}}[\tau (\lambda )],\end{eqnarray} \tag{ 1 }$$

where M_d is the mass of dust in the sphere, κ(λ) is the mass absorption coefficient of the dust, and B_ν(λ, T_d) is the Planck function at wavelength λ and dust temperature T_d. The function P_esc is the probability for an IR photon to escape from the sphere, given by Cox & Mathews (1969) and Osterbrock & Ferland (2006):

$$\begin{eqnarray}&&{P}_{\mathrm{esc}}(\tau )=\displaystyle \frac{3}{4\tau }\ [1-\displaystyle \frac{1}{2{\tau }^{2}}+(\displaystyle \frac{1}{\tau }+\displaystyle \frac{1}{2{\tau }^{2}}){e}^{-2\tau }]\end{eqnarray} \tag{ 2 }$$

where τ is the radial optical depth of the sphere given by

$$\begin{eqnarray}&&\tau (\lambda )=\displaystyle \frac{3}{4}\ (\displaystyle \frac{{M}_{{\rm{d}}}}{\pi {R}^{2}})\ \kappa (\lambda ).\end{eqnarray} \tag{ 3 }$$

Limiting values for P_esc(τ) are

$$\begin{eqnarray}{P}_{\mathrm{esc}}(\tau ) & \approx & 1-\displaystyle \frac{3\tau }{4}+\displaystyle \frac{2{\tau }^{2}}{5}-\displaystyle \frac{{\tau }^{3}}{6}\;\qquad \mathrm{when}\ \tau \ll 1\\ & \longrightarrow & \displaystyle \frac{3}{4\tau }\qquad \qquad \qquad \qquad \ \ \ \mathrm{when}\ \tau \;\longrightarrow \;\infty .\end{eqnarray} \tag{ 4 }$$

The escape formalism explicitly takes the distribution of optical depth in the sphere into account. As τ becomes very large, P_esc(λ) ∼ 1/κ(λ) (Equations (3) and (4)), with the consequence that the spectral signatures of the dust are cancelled in Equation (1). Furthermore, all information on the dust mass in the sphere becomes indeterminate.

Equations (1) and (3) can be readily generalized to a dust population characterized by different compositions, masses, and temperatures. The flux is then given by the sum over dust components:

$$\begin{eqnarray}&&{S}_{\nu }(\lambda )=\displaystyle \sum _{j}\ \displaystyle \frac{4{M}_{{\rm{d}},j}\ {\kappa }_{j}(\lambda )}{4\pi {D}^{2}}\ \pi {B}_{\nu }(\lambda ,{T}_{{\rm{d}},j})\ {P}_{\mathrm{esc}}[\tilde{\tau }(\lambda )]\end{eqnarray} \tag{ 5 }$$

where $\tilde{\tau }$ is the optical depth in the presence of multiple dust components:

$$\begin{eqnarray}&&\tilde{\tau }(\lambda )=\displaystyle \frac{3}{4\pi {R}^{2}}\ \displaystyle \sum _{j}{M}_{{\rm{d}},j}\;{\kappa }_{j}(\lambda )\end{eqnarray} \tag{ 6 }$$

where M_d,j, κ_j(λ), and T_d,j are the masses, mass absorption coefficients, and temperatures corresponding to the different grain compositions j.

In the high optical depth limit, S_ν(λ) will still be a blackbody, as long as all dust components are at the same temperature. If the temperatures differ, then the emission features, generated by κ(λ), and the absorption features, generated by P_esc(λ), will not cancel perfectly, and the total spectrum will still retain an imprint of the individual dust spectra. The former case is illustrated by our fit at day 775, and the latter case by fits at day 615 and 1144, as shown later in Figures 3 and 4, and Table 2.

Table 2.  Evolution of Dust Parameters and Ejecta Radius^a

Epoch (day)	615	775	1144	8815	9090
Pure silicates (MgSiO3)
Md(${M}_{\odot }$)	0.40	0.40	0.40	...	...
Td(K)	607	334	145	...	...
Pure amorphous carbon
Md(${M}_{\odot }$)	0.047	0.047	0.047	...	...
Td(K)	454	333	120	...	...
Composite grainsb
Md(${M}_{\odot }$)	...	...	...	0.42 ± 0.07	0.45 ± 0.08
Td(K)	...	...	...	26.3	25.2
Elemental abundances in dust
Mg(${M}_{\odot }$)	0.096	0.096	0.096	0.090	0.095
Si(${M}_{\odot }$)	0.11	0.11	0.11	0.11	0.11
C(${M}_{\odot }$)	0.047	0.047	0.047	0.045	0.048
Rej (1016 cm)c	0.48	0.61	0.90	6.9	7.1
LIR (${L}_{\odot }$)d	3.6 × 105	1.2 × 105	7.0 × 103	310 ± 27	270 ± 19
τ(λ)e	9130	5650	252	0.33	0.33

Notes.

^aEntries in bold were held constant. The spectra from the last two epochs were fit with ellipsoidal composite grains. Grains in the first three epochs consist of a mixture of spherical silicates (MgSiO₃), and amorphous carbon (AC) grains. ^bThe composite grains consist of an MgSiO₃ matrix with AC inclusions occupying 18% of the grain's volume. ^cThe ejecta radius, R_ej, at each epoch was calculated for a constant expansion velocity of 910 km s⁻¹. ^dFor comparison, on day 1144 the U to 20 μm (uvoir) luminosity is 3400 ${L}_{\odot }$ (Suntzeff et al. 1991). ^eThe optical depth, τ(λ) is given by Equation (3), and calculated at around the peak wavelength of the emission for each epoch, at 20, 20, 50, 200, and 200 μm, respectively.

Download table as:  ASCIITypeset image

3.1. Nucleosynthetic Constraints on the Dust Mass

3.2. Abundances Constrained Fits to the Infrared Spectra of SN 1987A

The inner ejecta of the SN was well resolved with the Atacama Large Millimeter/Submillimeter Array (ALMA) at far-IR wavelengths, and found to have an elliptical shape with dimensions corresponding to free expansion velocities of 1350 ± 150 and 750 ± 250 km s⁻¹, along the major and minor axes of the ellipse (Indebetouw et al. 2014). For simplicity, we will adopt a spherical shape for the ejecta with a constant expansion velocity of $v={(1350\times {750}^{2})}^{1/3}=910$ km s⁻¹. The radius of the ejecta at the different epochs is given in Table 2. On day 9000 it is about 7 × 10¹⁶ cm, and its radial 200 μm optical depth is ∼0.3. Since the ejecta is optically thin at these epochs, we start our analysis with the spectra of those epochs.

3.2.1. Days 8515 and 9090

High resolution ALMA images at 450 and 850 μm clearly resolved the emission from the ejecta and the circumstellar ring (Indebetouw et al. 2014; Zanardo et al. 2014). Spectral observations, summarized by Matsuura et al. (2015), showed that the far-IR/submillometer fluxes obtained by Herschel and ALMA observations are dominated by continuum emission from ejecta dust, except for the 70 μm flux, which may be contaminated by the [O I] 63 μm line and by the continuum emission from the hot dust in the circumstellar ring (Dwek et al. 2010). The combined Herschel and ALMA measurements of the far-IR and submillimeter flux densities of SN 1987A used in our analysis were taken from Matsuura et al. (2015, their Table 1). We adopt here their estimate that at least two-thirds of the observed flux in the 70 μm band arises from the ejecta dust.

To minimize the mass of dust in the ejecta, Matsuura et al. (2015) considered the possibility that the dust grains are in the shape of elongated ellipsoids. Elongated grains radiate more efficiently at long wavelengths, thereby lowering the mass needed to produce the emission. However, the shape effect in AC grains does not only increase their long wavelength emissivity, it also considerably flattens it with wavelength. Since most of the emission in the Matsuura model is attributed to AC dust, the resulting spectrum became too broad to fit the observations. To circumvent these problems, we adopt a dust model consisting of ellipsoidal composite grains, comprised mostly of a silicate (MgSiO₃) matrix with AC inclusions. The elongated shape increased the dust emissivity, and the addition of AC inclusions further increased the dust emissivity, with only minor flattening of its wavelength dependence.

Figure 1 depicts the various mass absorption coefficients considered in our models. In the Rayleigh limit, where $a\ll \lambda$, the mass absorption coefficient of the dust is independent of grain radius. The left panel shows the wavelength dependence of κ(λ) for pure spherical MgSiO₃ (blue curve) and AC (red curve) grains. Optical constants for the silicate dust were taken from Jäger et al. (2003), and those for the AC dust from (Zubko et al. 1996). The right panel depicts the wavelength dependence of the ellipsoidal composite grains. The blue curve corresponds to pure MgSiO₃ grains. The black curves give the values of κ for the composite grains with different amounts of AC inclusions (see figure caption). The green curve corresponds to the mass absorption coefficient of the composite grains used to fit the 8515 and 9090 epoch spectra.

Zoom In Zoom Out Reset image size

We assumed that the dust radiates at a single temperature and used the IDL routine MPFIT (Markwardt 2009) to simultaneously fit the SN 1987A flux densities for the two epochs, with the dust masses and temperatures as variables. We experimented with different concentrations of AC inclusions. The results presented here are for composite grains with MgSiO₃ and AC volume filling factors of f_s = 0.82 and f_c = 0.18, respectively. The choice of parameters was driven by the constraints on the available mass of Mg and Si in the ejecta. The resulting spectra are shown in Figure 2, and the resulting dust masses and temperatures are summarized in Table 2.

Zoom In Zoom Out Reset image size

The mass fraction of MgSiO₃ and AC in the grain are given by ${f}_{{\rm{s}}}\;({\rho }_{{\rm{s}}}/\bar{\rho })$ and ${f}_{{\rm{c}}}\;({\rho }_{{\rm{c}}}/\bar{\rho })$, respectively, where ρ_s = 3.2 g cm³ and ρ_c = 1.8 g cm³ are, respectively, the mass densities of MgSiO₃ and AC grains, and $\bar{\rho }={f}_{{\rm{s}}}{\rho }_{{\rm{s}}}+{f}_{{\rm{c}}}{\rho }_{{\rm{c}}}$ is the average density of the grain. If M_d is the total dust mass, then the mass of Mg, Si, and C locked up in the dust is given by $0.24\;{f}_{{\rm{s}}}\;({\rho }_{{\rm{s}}}/\bar{\rho }){M}_{{\rm{d}}}$, $0.28\;{f}_{{\rm{s}}}\;({\rho }_{{\rm{s}}}/\bar{\rho }){M}_{{\rm{d}}}$, and ${f}_{{\rm{c}}}\;({\rho }_{{\rm{c}}}/\bar{\rho }){M}_{{\rm{d}}}$, respectively. For an AC filling factor of 0.18, we get $\bar{\rho }=2.95$ g cm³, MgSiO₃ and AC mass fractions of 0.89 and 0.11, respectively, and Mg and Si mass fractions of 0.21 and 0.25, respectively. The elemental masses locked up in the dust are presented in Table 2. About half of the MgSiO₃ mass consists of O which, because of its large abundance in the ejecta (∼1 ${M}_{\odot }$), does not impose any constraints on the mass of dust that can form in the ejecta.

The results show that the late-time spectra require ∼0.09 ${M}_{\odot }$ of Mg, ∼0.1 ${M}_{\odot }$ of Si, and ∼0.046 ${M}_{\odot }$ of C to be locked up in dust. From Table 1, we see that these masses are well within their expected yields, alleviating the abundance problems of previous models.

3.2.2. Days 615 and 775

At these early evolutionary stages we modeled the spectra of SN 1987A using data presented by Wooden et al. (1993). As mentioned before, the lack of spectral features has been construed as evidence for the lack of silicate dust in the ejecta. To test this conclusion, we ran models with pure spherical silicate (MgSiO₃) and AC grains, in which the mass of Mg and Si in the dust was taken to be equal to that derived for epochs 8515 and 9090. For this choice of dust masses, or any other sufficiently large dust mass, the ejecta is optically thick and, consequently, the shape of the dust grains will not affect the blackbody spectrum of the ejecta. Our choice of using distinct silicates and AC dust populations is motivated by the dust growth and coagulation models of Sarangi & Cherchneff (2015). In these models, dust grows primarily by accretion at early epochs, after which the abundance of condensible elements is sufficiently depleted and coagulation becomes the main process for grain growth.

Thus the only free parameters were the silicate and AC dust temperatures. The resulting spectra are shown in Figure 3. The red and blue curves are, respectively, the carbon and silicate emission components of Equation (5), but only their sum (green curve) is observable. The 20 μm optical depth of the ejecta is ∼10⁴ and 6000 for these two epochs (see also Table 2), giving rise to smooth featureless spectra. They provide a surprisingly good fit to the observations considering the fact that we forced both the silicate and carbon dust to radiate at single (but different) temperatures. The figure shows that the silicate 9.7 and 18 μm emission features have been internally absorbed by the silicate–carbon dust mixture, as evident from the strong absorption features in the carbon dust emission.

Zoom In Zoom Out Reset image size

3.2.3. Day 1144

We selected the 10 and 20 μm photometric data obtained with the CTIO on day 1144 to represent the observations taken during the 935–1352 day interval (Suntzeff et al. 1991). For day 1144, we also have a 50 μm upper limit reported by P. M. Harvey, and presented by Dwek et al. (1992). A low resolution (λ/Δλ = 40) 16–30 μm spectrum of the SN was obtained with the KAO around day 1151 (Dwek et al. 1992). Several emission lines were detected, and the upper limit set on the continuum intensity is consistent with the CTIO 20 μm detection. Figure 4 presents the CTIO data and the 50 μm upper limit. As for days 615 and 775, we assumed that the dust giving rise to the emission consists of a combination of spherical MgSiO₃ and AC grains with masses identical to this in epochs 615 and 775 and allowed the temperature to vary. Figure 4 shows the fit to the observed spectrum. The dust masses and temperatures are presented in Table 2.

Zoom In Zoom Out Reset image size

So far, we assumed that the ejecta can be approximated by a dusty sphere of radius R. However, the early detection of X- and γ-rays from the SN, and the high ionization rate in its envelope, required the macroscopic mixing of clumps containing radioactive Ni, Co, and Fe throughout the ejecta (Graham 1988; Li et al. 1993). Furthermore, the lack of silicate emission features in the early IR spectra (days 615, 775) was attributed to the fact that the dust may have formed in optically thick clumps (Lucy et al. 1989, 1991). In the following, we show that the results of our calculations can be simply extended to a clumpy medium.

We will construct a simple model for a clumpy ejecta. Since the ejecta is optically thin at late epochs (days 8515 and 9090), the model is relevant only for the epochs when the homogeneous ejecta was optically thick. Using the same dust mass and composition as in the homogeneous model, we show that the clumpy model can reproduce the same IR spectrum at each epoch and is also capable of accounting for the UVO region of the SN light curve.

We will assume the presence of N_c identical dusty clumps with a dust mass of m_c = M_d/N_c, where M_d is the total dust mass in the ejecta and radius r_c. We also assume that the clumps maintain pressure equilibrium with their surroundings so that the ratio r_c/R, where R is the outer radius of the ejecta, remains constant with time. We also assume the presence of some dust mass, m_ic, in the interclump medium.

The volume filling factor of the clumps is given by

$$\begin{eqnarray}&&{f}_{V}={N}_{{\rm{c}}}\;{x}^{3}.\end{eqnarray} \tag{ 7 }$$

The number of clumps and their radii are constrained by the requirement that the projected area of the clumps be equal to that of the the homogeneous sphere. This requirement reduces to the constraint that ${N}_{{\rm{c}}}\;{r}_{{\rm{c}}}^{2}=\xi \;{R}^{2}$, where ξ is a factor of the order of unity that compensates for the overlapping of clumps along a given line of sight.

The optical depth of each clump, τ_c, is simply related to τ_h, that of the homogeneous sphere, by

$$\begin{eqnarray}&&{\tau }_{{\rm{c}}}=\displaystyle \frac{3}{4}(\displaystyle \frac{{M}_{{\rm{d}}}/{N}_{{\rm{c}}}}{\pi {r}_{{\rm{c}}}^{2}})\kappa =\displaystyle \frac{{\tau }_{{\rm{h}}}}{\xi }.\end{eqnarray} \tag{ 8 }$$

The clumps are therefore optically thick at UVOIR wavelengths at all early epochs before day 1144, so that their cross section can be taken as equal to the geometrical value π r_c².

The probability that a photon passes a distance z through the ejecta without absorption, and then hits a clump in the z and z + dz interval is given by

$$\begin{eqnarray}&&P(z){dz}=\mathrm{exp}(-z/{\ell })\displaystyle \frac{{dz}}{{\ell }}\end{eqnarray} \tag{ 9 }$$

where ${\ell }={({n}_{{\rm{c}}}\pi {r}_{{\rm{c}}}^{2})}^{-1}$ is the photon mean free path through the ejecta, and n_c is the number density of clumps. A photon incident on the sphere within an annulus 2π b db, where b is the projected distance from the center of the sphere, will strike a clump with a probability given by

$$\begin{eqnarray}P(b)\ {db} & = & {\displaystyle \int }_{0}^{2R\sqrt{1-{(b/R)}^{2}}}\ \mathrm{exp}(-z/{\ell })\displaystyle \frac{{dz}}{{\ell }}\\ & = & 1-\mathrm{exp}(-2R\;\sqrt{1-{(b/R)}^{2}}).\end{eqnarray} \tag{ 10 }$$

The fraction f_A of the projected area of the ejecta that is filled by the clumps is then given by the integral:

$$\begin{eqnarray}&&{f}_{{\rm{A}}}={(\pi {R}^{2})}^{-1}{\displaystyle \int }_{0}^{R}\ P(b)2\pi \;b\;{db}.\end{eqnarray} \tag{ 11 }$$

Table 3 summarizes some clump characteristics and the fraction of the projected surface of the ejecta covered by clumps for several values of N_c and ξ. The most desirable combination of these parameters is one that gives a small volume filling factor for the clumps, and a large surface covering factor. The combination of {N_c, ξ} = {1000, 2} fills both requirements. It gives a clump filling factor of 0.09 and a surface covering fraction of 0.82. The value of r_c/R for this model is 0.044, consistent with the scales of clumps predicted by numerical simulations of the instabilities in core collapse supernovae (Wongwathanarat et al. 2015). This means that in the clumpy model, the integrated brightness of the IR emission will be lower by 18% compared to that of the homogeneous sphere. This factor is well within the uncertainties in the fluxes observed on days 615–1144. So all dust masses derived assuming homogeneous ejecta are still valid for the clumpy ejecta model. Figure 5 presents Monte Carlo simulations of a clumpy ejecta for N_c = 1000 clumps and values of ξ = 1.0, 2.0, and 4.0 (top to bottom). The simulations were performed to validate the analytically derived area covering fractions listed in Table 3. From top to bottom, the fraction of the projected area covered by the clumps is 0.59, 0.82, and 0.96, in good agreement with the analytical formula presented in Equation (11).

Zoom In Zoom Out Reset image size

The clumpy model must allow for the escape of UVO photons and any selective extinction in the ejecta. In the following, we use the mega-grain approximation described by Városi & Dwek (1999) to calculate the escape of these photons from the clumpy medium. Equation (59) in their paper can be simplified by assuming that the volume filling factor, f_V, of the clumps is small, that most of the dust resides in the clumps, and that the clumps are solid spheres with a geometrical cross section $\pi {r}_{{\rm{c}}}^{2}$. The radial optical depth, τ_R(λ), of the ejecta is then given by the sum of two terms. The first is wavelength independent and caused by the absorption in the optically thick clumps, and the second is the wavelength dependent absorption by dust residing in the interclump medium.

$$\begin{eqnarray}{\tau }_{{\rm{R}}}(\lambda ) & = & {n}_{{\rm{c}}}\;\pi {r}_{{\rm{c}}}^{2}\;R+{\rho }_{\mathrm{ic}}\;\kappa (\lambda )R\\ & = & (3/4)[\xi +({m}_{\mathrm{ic}}/\pi {R}^{2})\;\kappa (\lambda )]\end{eqnarray} \tag{ 12 }$$

where ρ_ic and m_ic are, respectively, the density and mass of the dust in the interclump region, and V is the volume of the ejecta.

At early times, the IR emission accounts for half the energy deposited in the ejecta by the radioactive decay of ⁵⁶Co. So the unabsorbed UVO luminosity can be at most equal to that of the IR. The observed UVO flux constitutes about 5%–10% of the IR emission from the SN (Lucy et al. 1989, 1991; Wesson et al. 2015, Figures 1 and 2). So at least 5%–10% of the unabsorbed UVO emission must escape the interclump medium. These, and observations of nebular emission lines from the ejecta (Jerkstrand et al. 2011), constrain the mass of dust that can be present in this medium.

For radial optical depths ≳1, the escape probability from the ejecta is approximately given by 3/4τ. For ξ = 2, the first term in Equation (12) contributes 3ξ/4 = 1.5 to τ_R, giving P_esc ≈ 0.5. For P_esc to be larger than 0.05, the optical depth must be smaller than ∼15, setting a limit of τ ≲ 13 on the second term in e.g., (12). The mass absorption coefficient of the dust at 0.55 μm is about 10⁴ cm² g⁻¹ for a mixture of silicate–carbon grains (Figure 1), so the second term becomes

$$\begin{eqnarray}&&(3/4)({m}_{\mathrm{ic}}/\pi {R}^{2})\;\kappa (V)\approx 2\times {10}^{5}\;{m}_{\mathrm{ic}}({M}_{\odot })\lesssim 13\end{eqnarray} \tag{ 13 }$$

where we assumed an ejecta radius of 0.5 × 10¹⁶ cm. So the interclump dust mass is limited to m_ic ≲ 10⁻⁴ ${M}_{\odot }$, allowing UVO emission that is not absorbed in the clumps to escape the ejecta.

In summary, the clumpy model can reproduce most of the early IR emission from the ejecta with the same mass of dust as in the homogeneous model. The clumps are optically thin at late times, so that the mass derived for days 8515 and 9090 is unaltered as well. The mass of dust in the interclump region is small compared to that in the clumps, allowing for the escape of nebular and continuum UVO emission from the ejecta. Because of the small mass of dust in the interclump medium, it makes no significant contribution to the IR emission from the clumps.

5.1. The Evolution of the IR Luminosity

5.2. External Energy Sources

5.3. Internal Energy Sources

5.3.1. Compact Central Object

One internal source of energy would be from accretion onto or spin down of a pulsar. Searches for a compact remnant have been unsuccessful (except possibly Middleditch et al. 2000), and have set low luminosity limits on a compact source (≲2 L_⊙; Graves et al. 2005). However, this could be higher if the extinction is far larger than was assumed.

5.3.2. Decaying Radioactive Elements

Observations of the bolometric light curve during the first two years after the explosion, the SN light curve was powered by the energy released in the radioactive decay chain ⁵⁶Ni $\longrightarrow$⁵⁶Co $\longrightarrow$⁵⁶Fe (Suntzeff & Bouchet 1990; Bouchet et al. 1991; Suntzeff et al. 1992). Models show that ∼0.07 ${M}_{\odot }$ of ⁵⁶Ni were required to fit the early UVOIR light curve, before dust has formed in the ejecta (Woosley et al. 1988; Seitenzahl et al. 2014, and references therein). After about 2000 days the SN 1987A light curve is predominantly powered by the decay of ⁴⁴Ti (τ_1/2 = 58.9 year). Efforts to determine the amount of ⁴⁴Ti required to power the late-time light curve have focused on hard X-rays, nebular lines, and UVOIR observations of the SN. A detailed summary is given by Seitenzahl et al. (2014), with highlights given in Table 4 of this paper. Until day ∼1800, the luminosity was determined from the UVOIR observations, which extend out to the 20 μm band. However, most of the dust emission after day ∼1100 is emitted at longer wavelengths (see Figure 4), thus the UVOIR luminosity may be missing a significant fraction of the total energy emitted by the ejecta.

Table 4.  Inferred ⁴⁴Ti Yield from UVOIR and FIR Observations of SN 1987A

M(44Ti)	Observations	Epoch	References
(10−4 ${M}_{\odot }$)		(days)
3.1 ± 0.8	hard X-rays	3600–4000	(1)
1.5	UVOIR lines	∼3000	(2)
0.55 ± 0.17	UVOIR	<1900	(3)
	extrapolated V	<4300
1.6	far-IR dust	8500–9100	this work

References. (1) Grebenev et al. (2012), (2) Jerkstrand et al. (2011), (3) Seitenzahl et al. (2014).

Download table as:  ASCIITypeset image

Figure 6 compares model light curves to the IR luminosities (yellow diamonds). The ⁵⁶Ni, ⁵⁷Ni, ⁵⁵Co, and ⁶⁰Co, and masses were taken to be equal to 0.071, 4.1 × 10⁻³, 9.2 × 10⁻⁶, and 4.5 × 10⁻⁸ ${M}_{\odot }$, respectively (Seitenzahl et al. 2014). The figure shows that at early epochs (days 615 and 775), the IR emission accounts for about half of the radioactive energy deposited by the decay of ⁵⁶Co in the ejecta. At later epochs, the IR luminosity can provide an important measure of the amount of ⁴⁴Ti created in the explosion. The different curves correspond to different masses of ⁴⁴Ti with values of 0.55 × 10⁻⁴ (Seitenzahl et al. 2014, red curve), and 3.1 × 10⁻⁴ (Grebenev et al. 2012, blue curve). Also shown in the figure is the luminosity of 390 ${L}_{\odot }$ on day 6067, derived from the T-ReCS observations of the SN ejecta (Bouchet et al. 2004, open diamond). The black curve provides the best fit to the IR luminosities derived in this paper and requires M(⁴⁴Ti) = 1.6 × 10⁻⁴ ${M}_{\odot }$.

Zoom In Zoom Out Reset image size

This value is in good agreement with the 1.5 × 10⁻⁴ ${M}_{\odot }$ derived by Jerkstrand et al. (2011) from the analysis of UVO emission lines in the ejecta around day 3000. However, in their standard model they adopted a dust optical depth of only 1. The IR optical depth on that day could be around 10 at the peak of the IR emission and significantly higher at UVO wavelengths. The ⁴⁴Ti mass inferred from the energy required to power the nebular emission lines should therefore be increased to account for the effect of dust absorption.

On the other hand, the ⁴⁴Ti mass inferred from the IR emission alone should also be modified: increased to account for any UVO emission, and decreased to account for the contribution of external heating sources to the emission. The mass should definitely be increased by ∼2%, to account for the UVO luminosity of ∼7 ${L}_{\odot }$ observed on day 8328 (Larsson et al. 2013, Figure 5), and at most decreased by ∼20% to account for the possible contribution from external heating sources. So the IR observations provide a robust lower limit of ∼1.3 × 10⁻⁴ ${M}_{\odot }$ on the mass of ⁴⁴Ti synthesized in the explosion.

In this paper, we presented a new scenario for the evolution of dust in SN 1987A. In this scenario, the observed IR spectrum on days 8515 and 9090 arises from ∼0.45 ${M}_{\odot }$ of composite ellipsoidal grains, consisting of a mixture of MgSiO₃ and AC with volume filling factors of 82% and 18%, respectively. In contrast with previous models (Matsuura et al. 2015), this model requires only ∼0.1 ${M}_{\odot }$ of Mg and a similar mass of Si to be locked up in dust, consistent with the post-explosive mass of these elements in the SN ejecta. The mass of carbon locked up in dust is only 0.04 ${M}_{\odot }$, in contrast to the 0.5 ${M}_{\odot }$ of C in the Matsuura et al. (2015) models.

We also show that the same mass of Mg, Si, and C could have already been locked up in the form of spherical MgSiO₃ and AC grains since day 615 in IR-optically thick clumps. The low mass of dust and the upper limit on the silicate dust in previous models were inferred from the absence of the 9.7 and 18 μm silicate emission features, assuming optically thin clumps. However, we show that the absence of the silicate absorption features at the early epochs is the result of the large optical depth and self-absorption in the clumps. Figures 2–4 and Table 2 show, respectively, the fit to the observed spectra, and the dust parameters corresponding to the fits.

The model presented in this paper thus solves three problems that plague previous models (Matsuura et al. 2015; Wesson et al. 2015, and references therein):

• 1.  
  the excessive mass of silicate or AC dust required to fit the late-time spectrum,
• 2.  
  the need for cold accretion to grow the grains from ∼10⁻³ to ∼0.5 ${M}_{\odot }$, and
• 3.  
  the absence of SN-condensed silicate dust.

The results presented in this paper are significantly different from the conclusion reached by Ercolano et al. (2007a) and Wesson et al. (2015), who firmly ruled out the possible presence of more than ∼3 × 10⁻³ ${M}_{\odot }$ of dust at the early epochs before ∼ day 1200. The main difference between our respective models is that their models were constructed to simultaneously fit the observed UVOIR as well as the IR dust emission. Without this constraint, our models can fit the mid- to far-IR emission at early times with dust having higher temperatures and densities, being optically thick through far-IR wavelengths. At the high densities prevailing at early epochs, the dust is thermally coupled to the gas, and more likely to be collisionally heated by the ambient hot gas, rather than radiatively heated by UVO photons. In either case, the ultimate heating sources are the decaying radioactivities in the ejecta. This scenario seems to be confirmed by the fact that the late-time IR spectra, when the ejecta is optically thin, were fitted with a very narrow range in dust temperatures. Such temperature distributions can be more readily explained if the dust is collisionally heated in the ejecta. In our scenario, and UVOIR emission solely arises from the sources in the interclump region, where the optical depth is very low. Only ∼10% of the energy is emitted at these wavelengths after the dust formation epoch.

The IR observations provide a robust lower limit of ∼1.3 × 10⁻⁴ ${M}_{\odot }$ on the mass of radioactive ⁴⁴Ti that has formed in the explosion, with a more likely value of 1.6 × 10⁻⁴ ${M}_{\odot }$, depending on the contribution of external sources in heating the dust.

The proposed scenario requires the very rapid evolution of dust in the SN ejecta. Theoretical calculations show that only ∼3 × 10⁻³ of dust can form by day 700 assuming a smooth ejecta, a value that increases to ∼0.08 ${M}_{\odot }$ for clumpy ejecta (e.g., Sarangi & Cherchneff 2015). Further increase in the dust mass can be achieved by increasing the density in the clumps, and by the inclusion of molecular and dust-continuum cooling, processes that were ignored during the dust formation epoch.

Finally, our model casts doubts on interpretations for the existence of an evolutionary trend in the dust mass from young SNe to older remnants. Infrared observations of young SNe need to carefully determine the origin of the emission, excluding any possible contributions from IR echoes and shocked interstellar dust. An extensive review of the observations of the IR emission from SNe, and their role as dust sources in the early universe was presented by Gall et al. (2011b). In particular, Figure 2 of Gall et al. (2011b) and Figure 4 in Gall et al. (2014) show that the inferred dust mass in SNe that are less than ∼3 year old, is between ∼(0.5–50) × 10⁻⁴ ${M}_{\odot }$, whereas the mass of SN condensed dust inferred from observations of young unmixed remnants, such as Cas A and the Crab Nebula is about 0.1 ${M}_{\odot }$ (Barlow et al. 2010; Gomez et al. 2012; Temim & Dwek 2013; Arendt et al. 2014). This "trend" could be interpreted as an evolutionary one, suggesting the growth of dust grains in the ejecta with time. However, as noted by Gall et al. (2011b), the low dust mass in young SNe could also be an observational selection effect, since it is inferred from only near-IR observations, leaving a large part of the dust spectrum unexplored. In this paper, we showed that even when the mid-IR region of the spectrum is observed, optical depth effects can result in a large underestimate of the dust mass in the ejecta.

We thank Arkaprabha Sarangi for helpful discussions and the referee for useful criticism that led to the improvement of the manuscript. This work was supported by NASA Astrophysical Data Analysis Program 13-ADAP13-0094.