The timing and location of dust formation in the remnant of SN 1987A
The discovery with the Herschel Space Observatory of bright far infrared and submm emission from the ejecta of the core collapse supernova SN 1987A has been interpreted as indicating the presence of some 0.4–0.7 M of dust. We have constructed radiative transfer models of the ejecta to fit optical to far-infrared observations from the literature at epochs between 615 days and 24 years after the explosion, to determine when and where this unexpectedly large amount of dust formed.
We find that the observations by day 1153 are consistent with the presence of 310M of dust. Although this is a larger amount than has previously been considered possible at this epoch, it is still very small compared to the amount present in the remnant after 24 years, and significantly higher dust masses at the earlier epochs are firmly ruled out by the observations, indicating that the majority of the dust must have formed at very late times. By 8515-9200 days after the explosion, 0.6–0.8 M of dust is present, and dust grains with radii greater than 2 m are required to obtain a fit to the observed SED. This suggests that the dust mass increase at late times was caused by accretion onto and coagulation of the dust grains formed at earlier epochs.
These findings provide further confirmation that core collapse supernovae can create large quantities of dust, and indicate that the reason for small dust masses being estimated in many cases is that the vast majority of the dust forms long after most supernovae have been detectable at mid-infrared wavelengths.
keywords:supernovae: individual: SN 1987A – ISM: supernova remnants – Magellanic Clouds – Supernovae: general
The discovery that quasars at redshifts 5 can already contain up to 10M of dust (e.g. Bertoldi et al. 2003; Dwek et al. 2007) triggered vigorous research aimed at identifying the sources of the dust. Massive stars have long been considered the most likely source due to their rapid progress to the late stages of stellar evolution where they may become significant dust producers (Cernuschi et al. 1967). In particular, dust is predicted to be formed in core collapse supernovae (CCSN), in quantities sufficient to account for the dust observed in the distant universe: a number of authors have modeled dust formation processes inside CCSN ejecta, concluding that at least 0.1 M of dust can condense (e.g. Todini & Ferrara 2001; Cherchneff & Dwek 2010; Sarangi & Cherchneff 2013), while other investigators (e.g. Morgan & Edmunds 2003; Dwek et al. 2007) have estimated that CCSN dust yields of between 0.1 and 0.5 M per CCSN are necessary to account for the observed quantities of dust in high-z galaxies. However, many mid-infrared studies of warm dust formation in core collapse supernova ejecta in the local universe have found that the amount of dust forming in the first few years appeared to be much smaller than required (Ercolano et al. 2007a; Andrews et al. 2010, 2011; Gallagher et al. 2012; Pozzo et al. 2004; Smith et al. 2008).
Dust formation in the ejecta of CCSNe may be detected by a) a steepening of the decline in the optical brightness, as the extinction increases due to new dust formation; b) the development of an infrared excess due to optical photons being absorbed and re-emitted at longer wavelengths; and c) a blue-shifting of emission lines from the ejecta, due to the red-shifted emission from the far side of the ejecta being absorbed by dust along the line of sight. SN 1987A was one of the first CCSNe in which dust formation was definitively detected (Lucy et al. 1989; Wooden et al. 1993; Meikle et al. 1993; Bouchet & Danziger 1993). Early studies assuming a smooth dust distribution estimated the presence of about 10M of dust. However, evidence for clumping of ejecta was discussed by Lucy et al. (1991), and Ercolano et al. (2007a) carried out an analysis using 3D radiative transfer models of clumped ejecta, finding somewhat more dust than analyses assuming a smooth distribution. However, the amount was still far lower than theoretical predictions and, if typical of CCSNe, would mean that they could not contribute significantly to the dust budget of distant quasars.
SN 1987A was unexpectedly detected in Herschel Space Observatory observations of the Large Magellanic Cloud (LMC) taken in 2010 (Matsuura et al. 2011). The far infrared and submm emission from the position of the supernova was consistent with continuum emission from between 0.4 and 0.7 M of dust, vastly more than was thought to have formed within a few years of the explosion. More recent observations with the Atacama Large Millimeter Array (ALMA) confirmed that the emission was located within the expanding supernova ejecta and not in the pre-existing circumstellar material (Indebetouw et al. 2014), while follow-up spectroscopy using the Fourier-Transform Spectrograph of the SPIRE instrument on board Herschel showed that the SPIRE broad-band images were dominated by continuum emission, with an upper limit of 8.4% on the contribution of emission lines to the measured broadband fluxes (Kamenetzky et al. 2013). Independent evidence for dust in the ejecta at these late times is provided by the asymmetry in Balmer line profiles observed by Fransson et al. (2013) on days 7982–7998, indicating that internal dust was preferentially absorbing optical emission from the far side of the ejecta.
The recent unambiguous detections of large amounts of dust in this and other supernova remnants (e.g. Gomez et al. 2012, Barlow et al. 2010) show that, in at least some cases, supernovae produce quantities of dust easily sufficient to explain the quantities of dust found in distant quasars. However, the discovery raises the question of when and how this large amount of dust formed in SN 1987A, and whether similar amounts of dust could have formed in the ejecta of other supernovae, despite the generally small quantities measured within 5 years of the explosions for CCSNe previously studied at mid-IR wavelengths. Are the dust masses of 10–10M determined at these early times unreliable for any reason, with large amounts of dust forming unseen within the first few years after the outburst? Or does the majority of the dust form at later times?
In this study we present new radiative transfer models of the ejecta of SN 1987A at seven epochs between 3 and 24 years after its outburst. The self-consistent models at each epoch demonstrate the physical plausibility of the large dust mass estimated 24 years after the explosion, and place strong constraints on when this dust could have formed.
2 Observational data
We construct models to fit the observations at seven epochs. We use 0.3–30 m spectrophotometry from Wooden et al. (1993), obtained on days 615 and 775 after the explosion, and we supplement these data with far infrared observations from Harvey et al. (1989), who observed the supernova at 100 m on day 626, and at both 50 m and 100 m on day 793. We dereddened the Wooden et al. (1993) data in the same way as described in Ercolano et al. (2007a).
In addition to these two epochs which were modelled by Ercolano et al. (2007a), we also created models for five subsequent epochs, 1153, 1300, 1650, 8515 and 9200 days after the explosion. Models for the ejecta 1153 days after the explosion were constructed to match the data described in Dwek et al. (1992), which consists of Kuiper Airborne Observatory spectrophotometry from 16–30 m, combined with photometry in filters from U to Q taken at the Cerro Tololo Interamerican Observatory (U to I; Suntzeff & Bouchet (1990)) and the European Southern Observatory (J to Q; Bouchet & Danziger (1993)). We measured the fluxes plotted in Figure 1 of Dwek et al. (1992), excluding the points longward of 24.5 m which are affected by noise and line emission, to derive a flux of 0.530.10 Jy at a mean wavelength of 21.23.7 m. Finally at this epoch we construct our model to respect an upper limit at 50 m of 1.65 Jy at days 1179-83 (Harvey P., priv. comm.)
SN 1987A was detected at a wavelength of 1.3mm with the Swedish-ESO Submillimetre Telescope (SEST) on days 1285 and 1645 (Biermann et al. 1992); we supplemented these data with J-Q infrared photometry from Bouchet & Danziger (1993) at day 1317, and radio observations at 1.4, 2.4, 4.8, and 8.6 GHz on days 1270-1306 and 1637-1662 from Zanardo et al. (2010). At these radio wavelengths the emission from the remnant is dominated by synchrotron radiation. The synchrotron flux is observed to increase with time, while dust emission fades as the remnant expands and cools, so that synchrotron radiation may contribute significantly to longer wavelength fluxes at late times. We therefore used the observed radio fluxes to extrapolate a power law spectrum which we then added to the model predictions for dust emission at epochs after day 1150, to obtain a total predicted flux.
Finally, we construct models to match the observations presented in Matsuura et al. (2011), which were obtained with the Herschel Space Observatory on 30 April and 5 August 2010, days 8467 and 8564 after the explosion was observed; and newly obtained Herschel photometry from days 9090 and 9122 after the explosion (Matsuura et al., submitted) together with ALMA 440 and 870 m photometry obtained between days 9273 and 9387 (Indebetouw et al. 2014).
The observational data are summarised in Table 1.
|615, 632, 638||0.3-30||Wooden et al. (1993), Moseley et al. (1989)||KAO spectrophotometry|
|623-631||100||Harvey et al. (1989)||0.6 Jy|
|775||0.3-30||Wooden et al. (1993)||KAO spectrophotometry|
|791-795||50||Harvey et al. (1989)||2.5 0.85 Jy|
|100||Harvey et al. (1989)||1.35 Jy|
|1110||J-Q||Bouchet & Danziger (1993)||ESO|
|1144||U-I||Suntzeff & Bouchet (1990)||ESO/CTIO|
|1153||16-29||Dwek et al. (1992)||0.530.10 Jy ; KAO spectrophotometry|
|1179-1183||50||P. Harvey (priv. comm.)||1.65 Jy; KAO|
|1270-1306||1.4–4.8GHz||Zanardo et al. (2010)||ATCA|
|1283-1290||1300||Biermann et al. (1992)||7.59 2.45 mJy|
|1296||V-R||Walker & Suntzeff (1991)||CTIO|
|1317||J-L||Bouchet & Danziger (1993)||ESO|
|10.4 (N)||Bouchet & Danziger (1993)||0.062 Jy|
|18.5 (Q)||Bouchet & Danziger (1993)||0.82 Jy|
|1637-1662||1.4–8.6GHz||Zanardo et al. (2010)||ATCA|
|1644-1646||1300||Biermann et al. (1992)||18.8 4.0 mJy|
|8014||1.4–4.8GHz||Zanardo et al. (2010)||ATCA; flux densities at subsequent epochs obtained by extrapolation.|
|8467-8564||100–500||Matsuura et al. (2011)||Herschel photometry|
|9090-9122||70–500||Matsuura et al. (submitted)||Herschel photometry|
|9273-9387||450, 850||Indebetouw et al. (2014)||ALMA band 7 & 9 photometry|
We analyse the formation of dust in the ejecta from SN 1987A by constructing self consistent three dimensional models, illuminated by physically plausible heating sources, using the radiative transfer code mocassin (Ercolano et al. 2003, Ercolano et al. 2005, Ercolano et al. 2008). We model the ejecta as being highly clumped, with clump sizes of between 1-3% of the ejecta radius. We adopt this geometry based on the expectation that Rayleigh-Taylor instabilities with characteristic size of this order will rapidly develop at the onset of the explosion (Arnett et al. 1989). The hydrodynamical instabilities may be driven by the passage of the shockwave through density jumps in the first few thousand seconds after the explosion, or after a few days by the bursting of bubbles of gas heated by the decay of Ni (Basko 1994). Observational evidence for clumping in the ejecta of SN 1987A has been discussed by Lucy et al. (1991) and Biermann et al. (1990).
Ercolano et al. (2007a) developed two varieties of clumpy models, the first (referred to as clumpy I) representing the case in which the dust and heating source are present in both the clumps and the interclump medium, and the second (clumpy II) in which dust is located only in the clumps, and the heating source is located only in the interclump medium. In this work we consider only the clumpy II scenario; the dust mass that can be contained in optically thick clumps is maximised in this scenario, and thus we can place an upper limit on the amount of dust formed at each epoch prior to the Herschel observations, and thus determine the minimum amount of additional dust which must have formed during the 20 years between our third and fourth epochs. In our models, clumps are distributed stochastically with a number density proportional to r, and the clumps themselves have a density profile, modelled using a 5x5x5 subgrid inside the mother grid cell in which the clump is located. The number of clumps in the model is specified via a volume filling factor, , which represents the volume of the expanding shell occupied by the clumps. We investigated volume filling factors between =0.025 and 0.2. We also considered clump number densities proportional to r and r where r is the radius of the ejecta; the predicted optical and infrared peaks of the SED vary by 20%, with a steeper density profile resulting in a higher optical flux and lower infrared flux. The predicted SED was found to be insensitive to the density distribution in the clumps themselves.
We investigated various mixtures of amorphous carbon and silicate dust, using optical constants from Draine & Lee (1984) for silicates and Hanner (1988) for amorphous carbon. For all our models we adopted a distance to the LMC of 49.97 kiloparsecs (Pietrzyński et al. 2013). A list of all the model parameters is given in Table 2
|Model parameter||Parameter space|
|Inner radius||510cm (day 615) – 7.510cm (day 9200)|
|Outer radius/inner radius||5|
|Clump radius||R/30, R/45, R/60|
|Clump number density||n(r) r|
|Clump density profile||n r|
|Clump volume filling factor||0.025–0.2|
|Luminosity of heating source||200–5.710 L|
|Temperature of heating source||3000–7000 K|
|Spectrum of heating source||Blackbody|
|Dust mass||10–2.0 M|
|Fraction in interclump medium||0%|
|Fraction in clumps||100%|
|Dust species||Amorphous carbon, astronomical silicates, iron|
|Distribution of grain radii||n(a)(a-x), with (x+0.005)a(x+0.25) m; x=0,1,2,3,4,5 m|
We began by developing viable models for day 615. We then created models for subsequent epochs by evolving the day 615 models assuming that all parts of the ejecta expanded with constant velocity between epochs. Thus, the inner edge of the ejecta is expanding at 941 kms. We then varied the luminosity and dust parameters to achieve a fit to the observations of the subsequent epoch. The primary dust variable is the mass; if varying the dust mass could not give a good fit, we also varied the dust composition and the grain size distribution.
The dust heating source at the first three epochs is the radioactive decay of Co, with contributions from Ti and other radionuclides after day 900. The luminosity used in the models for these epochs was taken from the bolometric light curve of Suntzeff et al. (1991) which covers days 254–1352.
By day 8515, the decay of Ti is the dominant heating source, with possible additional contributions from X-rays produced in shock interactions between the ejecta and pre-existing circumstellar material, and the diffuse interstellar field. We discuss the heating luminosity at this epoch in Section 3.8.1.
3.1 Reproducing previous results
Using mocassin v2.02.70, we first attempted to reproduce the results of Ercolano et al. (2007a). In an erratum, Ercolano et al. (2007b) noted that mocassin version 2.02.21, which they had used to calculate all their models, was afflicted by a bug which affected the calculation of the IR SED. Probability density functions used for calculating dust emission were incorrect, while those used for absorption calculations were unaffected. This meant that while the shape of the calculated IR excess was incorrect, the optical SED was unaffected.
The bug was fixed with the release of v2.02.41 of mocassin; we recalculated the models presented in Ercolano et al. (2007a) using mocassin v2.02.70. Figure 1 shows the results of model calculations for amorphous carbon, created using the input parameters specified in Table 1 of Ercolano et al. (2007a), compared to the output obtained with mocassin v2.02.21. The optical portion of the reproduced SED is very similar to the original, while the shape of the IR SED differs, as expected from the nature of the bug. Figure 1 shows that the reproduced models give higher mid-infrared fluxes than the original models, and the shape of the SED significantly diverges from the original at wavelengths longer than about 10 m. The lower fluxes at longer wavelengths allow the possibility that a higher dust mass than previously reported could be necessary to fit the data. We therefore began by searching for ejecta parameters which better fit the data at day 615.
3.2 Day 615
Using the models originally presented in Ercolano et al. (2007a) as a basis, we varied the model parameters to regain a good fit to the IR SED. We used the same R and Y(=R/R) as adopted by Ercolano et al. (2007a) (5.0e15 cm and 5 respectively), as well as the same grain size distribution, with the number density of dust particles with a radius being proportional to between limits of 0.005 and 0.25 m. We then ran models for volume filling factors, , of 0.025, 0.05, 0.1 and 0.2, and for dust masses between 10 and 1.0 M. At this epoch, we investigated different clumping scales, creating models with clump radii of R/30, R/45 and R/60. Varying the clump size did not have a major effect on the SED. The number of grid cells required for a clumping scale of R/x is proportional to x and so the computing time required becomes prohibitively expensive as the clump size is decreased; at subsequent epochs we therefore modelled ejecta only with a clump size of R/30. This size is in agreement with theoretical predictions of the characteristic scale of the Rayleigh-Taylor instabilities predicted to form in supernova ejecta (eg Herant & Woosley 1994).
We found that an improved fit to the observed day 615 SED could be obtained by increasing the dust mass in the Ercolano et al. (2007a) models. The volume filling factor f is well constrained - as described in Ercolano et al. (2007a), for optically thick clumps in the clumpy II scenario, the emergent optical flux is determined solely by f, and as noted in Ercolano et al. (2007b), the optical flux in the models produced with mocassin 2.02.21 was not affected by the bug in that version. Therefore it remains the case that f=0.1 provides the best fit to the optical SED. Fits to the IR SED are then obtained by varying the dust mass only.
3.2.1 The carbon and silicate fractions
We investigated dust compositions at this epoch of pure carbon, pure silicates, and mixtures of silicates:carbon of 15:85 and 50:50. Plots of the predicted SED at day 615 for a dust mass of 0.001 M with each composition are shown in Figure 2. As can be seen in the figure, pure silicate models are ruled out by the data; at dust masses of 0.001 M, a reasonable fit to the 10–100 m region of the SED is obtainable, but the emission between 3–10 m is strongly underpredicted. Increasing the dust mass to 0.01 M results in an improved fit to the shorter wavelength region, but the 10-100 m emission is then strongly overpredicted.
Models with a 50:50 mixture of silicates and amorphous carbon provide a somewhat better fit to the data, but cannot fit all parts of the IR excess; lower dust masses result in an underpredicted 3–10 m region while higher dust masses overpredict the 10–100 m region. Any dust mass higher than 0.01 M overpredicts the 100 m observation.
Pure carbon models can reproduce all parts of the SED well. Models with a 15:85 ratio of silicates to carbon provide a slightly improved fit to the observations at 2-4 m. Figure 3 shows our best fitting model at day 615. The model has a volume filling factor of 0.1 and a clump size of R/30. The dust mass is 1.010M, approximately a factor of 5 higher than that found by Ercolano et al. (2007a); the best fitting amorphous carbon model at day 615 in Ercolano et al. (2007a) had a dust mass of 2.210M. From our models with slightly higher and lower dust masses compared to the best fitting model, we estimate an uncertainty on our derived dust mass of 2.010M.
Regardless of dust composition, the day 626 upper limit at 100 m provides a strong constraint on the dust mass, requiring that it be lower than 0.01 M in all cases. For the 15:85 composition which provides the best overall fit to the SED, the 100 m flux upper limit of 0.6 Jy (Harvey et al. 1989) is incompatible with any dust mass higher than 0.005 M (Figure 4).
3.3 Day 775
To fit the observations at day 775, we radially expanded the grid for the best fitting day 615 model by a factor of (775/615) and then varied the dust mass to achieve a fit. We find that this alone enabled a reasonable fit to the day 775 data, with no need to vary the volume filling factor. The dust mass required to match the data is 2.010M, consistent with further dust formation having occurred between days 615 and day 775, in the same region in which it had formed by day 615. As with our estimated day 615 dust mass, this mass is approximately a factor of 5 higher than found by Ercolano et al. (2007a) at this epoch with the same dust composition.
A slightly better fit to the optical SED is achieved with a volume filling factor of f=0.2 (see Figure 5). This could suggest that between days 615 and 775, dust formed in clumps which had previously been dust free. However, models with the same value of f as in day 615 are not ruled out and thus it is also possible that the extra dust mass came from further dust forming in the same clumps in which it had previously been present.
3.4 Dust mass limits at early epochs
With the discovery of the presence of large amounts of dust in the remnant of SN 1987A 23 years after the explosion, it is of great importance to determine when this large amount of dust formed. Studies which have quantified the amount of dust present in young (3 years old) supernova remnants including SN 1987A have typically found dust masses of 0.001 M or less. This raises the question of whether those low dust masses might have been underestimates.
In this study we find that in our best fitting model of the ejecta at day 615, the infared emission at 10 m is optically thick and insensitive to an increase in the total mass of dust for dust masses greater than about 0.001 M. However, the far infrared data points at days 615 and 775 provide strong upper limits on the dust mass that could be hidden inside dense clumps at these epochs, as shown in Figure 4. To obtain a flux below the 3 upper limit of 0.6 Jy at 100 m on day 615 with a dust mass as high as 1 would require that the volume occupied by the clumps be very small; the emission at 100 m is still optically thin with f=0.025 and M=0.01 M, and already exceeds the observed flux by a factor of 3.5. But models with volume filling factors below 0.025 significantly overpredict the optical fluxes. Models in which the volume filling factor gives optical fluxes consistent with the observations place a strong upper limit of 0.002 M on the dust mass at day 615. Figure 4 shows the model outputs for a volume filling factor of 0.1 and dust masses of 0.005, 0.01, 0.1 and 1.0 M at day 615, showing that the emission at 100 m becomes optically thick only for dust masses greater than 0.01 M, where the emission at 100 m already significantly exceeds the observed value.
3.5 Day 1153
To model the day 1153 observations we proceeded as before, expanding the day 615 grid by a factor of (1153/615) and varying the dust mass. The uncertainties on the L and M band fluxes reported by Bouchet & Danziger (1993) imply that the detections were significant only at the level of 2.5 and 1.6 respectively, and we therefore use 3 upper limits as constraints on our models. In their analysis of the day 1153 observations of SN 1987A, Dwek et al. (1992) took N and Q band fluxes from Bouchet & Danziger (1993), obtained at ESO, and rescaled them to match photometry obtained at the CTIO. Without this rescaling, we find that the measured fluxes cannot be matched by any model, with no combination of dust mass and luminosity predicting a flux as high as those observed. We therefore use the rescaled fluxes. The systematic uncertainties are thus likely to be considerably larger than the uncertainties reported by Bouchet & Danziger (1993).
We find that a uniformly expanded grid with a dust mass of 0.003 M gives the best overall match to the SED (Figure 6), although it underpredicts the N-band flux by approximately a factor of three. A lower dust mass results in a better fit to this point, with the predicted SED for M=0.001 M matching the observation well; however, the optical fluxes are then overpredicted significantly, and the dust mass is lower than found at the previous epoch. Increasing the dust mass results in cooler dust, and the N-band flux is further underpredicted.
We investigated whether varying the dust composition could improve the overall fit to the SED. While at the two previous epochs, the shape of the SED requires that the dust be predominantly carbon, the decrease of silicon fine structure emission line strengths as well as the disappearance of iron fine structure emission have been attributed to the formation of silicate and iron dust at later epochs (Dwek et al. 1992). The decrease in the fluxes of forbidden lines may also be due simply to the declining temperature of the gas in the expanding remnant; to investigate whether the observed SEDs could be consistent with dust of different species forming, we created models in which part of the dust was in the form of iron or silicates. A moderate silicate fraction is possible at this epoch, but the lower opacity of silicate dust in the optical compared to carbon dust means that when the silicate fraction is higher than 25%, the dust mass required to match the infrared SED results in the optical SED being overpredicted. Increasing the dust mass to provide sufficient opacity to match the optical SED results in an overprediction of the infrared fluxes. The silicate fraction at day 1153 must therefore still be much lower than the carbon fraction.
The mid-infrared fluxes at this epoch place strong constraints on the amount of dust that could be present. Figure 7 shows the predicted SED for dust masses between 0.002 and 0.008 M, in the clumpy shell consistent with expansion from previous epochs. As with the upper limit to the dust mass at earlier epochs, the dependence of the predicted optical flux on the covering factor of the dusty clumps constrains that parameter to a value inconsistent with optically thick clumps at long wavelengths. The dust mass is thus strongly limited to less than 0.006 M at day 1153, and we consider the most likely value to be 0.0030.001 M.
3.6 Day 1300
The remnant of SN 1987A was detected at a wavelength of 1.3mm, 1300 days after the explosion (Biermann et al. 1992, with a measured flux was 7.592.45 mJy. Models constructed by expanding the previous grids at constant velocity and increasing the dust mass cannot satisfactorily fit this point, if it is assumed that all of the 1.3mm flux arises from dust continuum emission; even for a dust mass of 0.5 M, the flux at 1.3mm is still underpredicted by a factor of 7. Besides dust continuum emission, synchrotron radiation and free-free emission may also contribute to the flux at 1.3mm. Using radio observations of the remnant from Zanardo et al. (2010) at 1.4–5.6GHz at which frequencies the emission is dominated by synchrotron emission, we estimate a synchrotron flux in Janskys at a wavelength in microns given by F()=0.00138 Jy. From this we find that the predicted synchrotron flux at 1.3mm is 0.03 mJy.
To estimate the free-free emission at day 1300, we use the following expression adapted from Wooden (1989):
where 1- is the ionised fraction of the gas, T is the gas temperatured, is the wavelength in microns and t is the time in days since the supernova explosion. Wooden (1989) estimated the free-free flux as coming from a partially ionised 10 M hydrogen envelope. The 3.6–1300 m SED can be well fitted by assuming a gas temperature of 1800 K and an ionised fraction of (1-)=0.037; in particular this reproduces the observation that the L band flux is considerably higher than the JHK fluxes. The 1.3mm flux is then dominated by free-free emission, and the amount of dust is poorly constrained. To fit the near-infrared observations, a dust mass of at least 0.02 M is required; lower dust masses do not give enough extinction and the JHK fluxes are overpredicted. The 1.3mm flux constrains the dust mass to less than 0.5 M, for the case of amorphous carbon. We therefore find that at day 1300, 0.02M0.5 M (Figure 8). Our estimated lower limit is already a factor of five higher than the 0.004 M estimated on day 1153.
3.7 Day 1650
SN 1987A was detected again at 1.3mm in observations obtained on days 1644-1646 by Biermann et al. (1992), with a measured flux of 18.84 mJy. As at day 1300, we used radio data from Zanardo et al. (2010) to estimate the synchrotron emission at shorter wavelengths; we find that F()=0.00741Jy, giving a predicted synchrotron flux at 1.3mm of 0.11 mJy. We also consider the likely free-free flux at this epoch: the free-free spectrum is flat between 1.3mm and the lower frequencies observed by Zanardo et al. (2010), and so the free-free flux is limited to less than 2.5 mJy, the flux reported at 8.6GHz. The flux that could be due to dust continuum emission is then 16.24 mJy. This is considerably higher than the fluxes predicted by any of our models even under extreme assumptions; for a dust mass of 1.0 M consisting of dust grains as large as 3.25 m in radius, the predicted dust continuum emission at 1.3mm is only 4.5 mJy (Figure 9). We conclude that the day 1644..1646 1300 m detection reported by Biermann et al. (1992) could not have been due to dust or free-free emission. A similar conclusion may then apply to their 1300 m detection on day 1300 (Section 3.6),
3.8 Day 8515
The discovery by Matsuura et al. (2011) of far-IR to submm emission at the position of SN 1987A on days 8467 and 8564 was unexpected. The emission was interpreted as being dust continuum emission, and analytical fits to the observed SED gave an estimated dust mass of 0.4–0.7 M, depending on the composition of the dust. The 440- and 870-m ALMA observations of Indebetouw et al. (2014) between days 9273 and 9387 confirmed the presence of at least 0.2 M of dust in the inner ejecta. Given the nucleosynthesis predictions for the estimated progenitor mass of 19 M (Thielemann et al. 1990, Woosley 1988), this implies that refractory elements in the ejecta condense into dust with an efficiency approaching 100%.
Subsequent observations with the SPIRE spectrometer on board Herschel have confirmed that the emission is indeed continuum emission. CO 6-5 and 7-6 line emission was detected with SPIRE, and the observations place an upper limit to the line contribution to the broadband fluxes of 8.4%. At the same time, line and continuum observations using ALMA at high spatial resolution confirm that this emission arises from within the ejecta, definitively ruling out any thermal echos from pre-existing structures (Kamenetzky et al. 2013; Indebetouw et al. 2014). The presence within the ejecta of a large amount of dust is thus confirmed. By considering our models at earlier epochs and expanding them appropriately, we can now determine whether the high dust mass estimated by Matsuura et al. (2011) is consistent with our own modelling of the evolution of the shell, and furthermore place constraints on when and where this dust formed.
3.8.1 The heating source
In previous epochs, the dust heating source was assumed to be from the radioactive decay of short-lived radioisotopes including Co. By 24 years after the explosion, much longer-lived isotopes, primarily Ti, dominate the radioactive heating, and other energy sources may be significant; additional heating is expected from the X-rays emitted by the interaction between supernova ejecta and circumstellar material, and from the ambient interstellar UV field.
The luminosity expected from the decay of Ti is approximately 400 L (Jerkstrand et al. 2011). The X-ray luminosity from ejecta-CSM interaction is approximately 500 L, although little of this X-ray emission is expected to go into heating the newly formed dust (Sturm et al. 2010). We constructed models for day 8515 using luminosities between 225 and 900 L. The infrared SED is insensitive to the assumed temperature of the heating source – models with heating source temperatures of 7000 K (used at previous epochs) and 3000 K predicted different optical SEDs but very similar IR SEDs.
3.8.2 Dust grain size distribution
To model the emission at day 8515, we expanded the day 615 grid as before, this time by a factor of (8515/615), and then varied other parameters to achieve a fit. Zanardo et al. (2010) presented radio observations up to day 8014; from the exponential fitting parameters in their Table 2, we extrapolated to day 8515, and find that the synchrotron flux at this epoch in Jy at a wavelength is given by F()=0.315 Jy. This predicts a synchrotron flux at 500 m of 11 mJy, compared to the upper limit from the Herschel observations of 29 mJy, and 9 mJy at 350 m compared to the observed flux of 49 mJy. The synchrotron flux is thus small but not negligible and we subtracted it from the fluxes tabulated in Matsuura et al. (2011) before proceeding.
The volume filling factor is not uniquely constrained by the observations at this epoch alone; when the SED is fitted with a given combination of and , a similar fit to the observations can be obtained by increasing the volume filling factor while reducing the luminosity or vice versa. Given that at each of the previous epochs, the volume filling factor was constrained by the observations to be 0.1, we set it to this value at day 8515.
Nucleosynthesis models appropriate for SN 1987A predict that if dust condensation from refractory elements were 100% efficient, then the maximum amount of dust that could be present in the remnant would be around 0.7 M (Thielemann et al. 1990; Woosley et al. 1988). When all parameters are kept the same as at previous epochs except for the luminosity and the dust mass, we are unable to fit the data well with any physically plausible combination of the two; the predicted dust SED was always warmer than observed unless the dust mass is greater than 1.0 M. We therefore considered the possibility that the large increase in dust mass by day 8515 was due to previously existing grains growing rather than new grains condensing.
It can be shown that the incremental growth of a grain by gas accretion, from a radius a to a larger radius a, is independent of the intial radius a, so that as accretion progresses, the grain size distribution will narrow if no fragmentation is also occurring. We therefore considered models which assume only gas accretion, linearly increasing by the same amount the size of all grains in the distribution used at earlier epochs, as well as models in which smaller dust grains remain present, either due to the formation of new dust grains, or the fragmentation of existing large grains.
The effect on the predicted SED of increasing the dust size is shown in Figure 10. We find that increasing the grain size can result in cooler dust as required by the observations, and a match to the observed SED can then be obtained by varying the dust mass and luminosity. As the grain size is increased, the mass of dust required decreases, while the luminosity required increases. For dust masses less than 1.0 M, a grain size distribution in which grains larger than 2 m in radius are present is required. Lower dust masses can also fit the observed SED but with even larger grain sizes required; to fit the data with a dust mass of 0.5 M requires grain radii of 3.005-3.25 m, while models with 0.4 M require grains larger than 5 m in radius. 0.4 M is the lowest possible dust mass which can still give a fit to the data; for lower dust masses, the luminosity required to match the Herschel observations is implausibly high, and also results in dust too hot, even for extremely large grain size distributions. Dust masses higher than 0.7M can give good fits to the SED with smaller grains; for a mass of 0.8M, micron-sized grains are required to fit the observed SED. Figure 11 shows fits obtained with models containing 0.4-1.0 M of dust.
The best fitting day 615 model contained 4.1510 dust grains, while the best fitting models at day 8515 contain a much lower number: 0.4M of dust with grain sizes between 5.005 and 5.25 m consists of 6.410 grains, while 1.0M of dust grains with 2.005a2.25 m consists of 2.2610 grains. This suggests that the increase in grain size must be due not only to accretion, but also to coagulation of grains.
The large grains in these models have a narrow range of temperatures. The distribution of grain temperatures throughout the model grids typically has a standard deviation of 1.5K. In a day 9200 model (see Section 3.9) containing grains around 5 m in radius, consisting of 85% amorphous carbon and 15% silicates, with a mass of 0.8M, the dust temperatures were 15.71.4 K and 19.21.5 K for the carbon and silicates respectively. The SED can thus be well fitted by a modified black body at a single temperature. In comparison, in the best fitting model at day 615 (Section 3.2), the temperature distributions were 25229 K for carbon, and 31631 K for silicates.
Having found that large dust grains are required for the dust in the models to be as cool as that observed at day 8515, we then investigated whether large grains could have been present at earlier epochs. Rerunning our best fitting models at earlier times with the enlarged grain sizes required to fit later epochs, we find that the models no longer fit the data. Figure 12 shows that at day 615, using the enlarged grain size distribution required by late-time data results in less optical extinction and cooler dust, resulting in over-predicted optical fluxes and under-predicted near-IR fluxes. The situation could be remedied by including more clumps in the model; however, one would then have to assume that the number of dusty clumps in the ejecta was not constant with time, and that more dusty clumps were present at earlier times than at later times. It seems unlikely that dust could form in some clumps but then be destroyed, while in other clumps dust formation could be prolific.
The largest mass of dust that could still match the day 1153 observations was 0.003 M, while the minimum amount of dust that must be present by day 8515 is 0.4M. This means that more than 99% of the dust in the ejecta formed after day 1153. This observational result stands in potentially sharp contrast to recent theoretical work by Sarangi & Cherchneff (2013), who predict that the bulk of dust formation should be complete within five years of the explosion. For the two results to be reconciled, dust formation in the ejecta would have had to have occurred almost entirely in the narrow time window between 3.2 and 5 years after the explosion.
3.8.3 Dust composition
Although there is observational evidence for iron and silicaceous dust formation in the remnant of SN 1987A, our models show that any non-carbon dust which has formed by day 1153 must be less than 25% of the total mass. For the Herschel wavelengths observed on day 8515, including iron and/or silicate dust causes the dust to become slightly warmer, so that a higher dust mass is required to restore the fit to the SED. A mass greater than 0.7 M may be difficult to reconcile with nucleosynthesis predictions for the explosion (Thielemann et al. 1990; Woosley et al. 1988). The dust is therefore likely to still be predominantly carbon at the latest epochs. The models shown in Figure 11 used dust consisting of 85% amorphous carbon and 15% silicates, implying a total carbon mass of 0.42 M and a silicate mass of 0.075 M for the 0.5 M model.
We considered the effect on the derived mass of different optical constants for amorphous carbon; all masses quoted thus far were calculated from models in which the optical constants of Hanner (1988) were used. To see if our results were sensitive to the choice of optical constants we re-ran all models using the optical constants from Zubko et al. (1996), specifically the ACH2 sample produced by arc discharge between amorphous carbon electrodes in an H atmosphere. We find that the difference at early epochs is significant, with a fit to the observed SED being possible with a lower dust mass. The Zubko optical constants have lower optical extinction such that the dust is cooler than found with the Hanner optical constants, and a smaller R and R are required to fit the SED. At day 615, a fit to the observations can be obtained with M=510 M and R and R both 0.7 the value previously derived values. At late times, however, the fit to the SED is less sensitive to the choice of optical constants. For the parameters obtained using Hanner optical constants, re-running the model with Zubko optical constants does not change the predicted SED significantly. For the models with smaller R and R found necessary to fit the day 615 observations using Zubko optical constants, when expanded to day 8515 we find that these models still require dust masses of the order of 0.60.2 M to fit the observed SED, and it remains the case that to obtain low enough dust temperatures while respecting the constraints on geometry obtained from the day 615 models, grain growth is required. However, unless the luminosity is considerably lower than found using Hanner optical constants, the dust is too warm to fit the observations. Our best fitting Zubko model had a luminosity of 315L, a dust mass of 0.65 M, and grain sizes between 3.005 and 3.25 m.
Thus, while the amount of dust present at each epoch is fairly sensitive to the choice of optical constants, our findings that grain growth has occurred within the ejecta, and that the vast majority of the dust formed more than 1000 days after the explosion, are not affected by the choice of optical constants. Figure 13 shows the predicted SEDs for both sets of optical constants, at days 615 and 8515 after the explosion.
|Epoch (days)||R (10 cm)||R (10 cm)||Grain size distribution (a (m), a, (m), power)||Estimated dust mass (M)|
|615||5.0||25.0||0.005-0.25, -3.5||0.0010 0.0002|
|775||6.3||31.5||0.005-0.25, -3.5||0.0020 0.0002|
|1153||9.4||46.9||0.005-0.25, -3.5||0.003 0.001|
|8515||69.2||346.1||3.005-3.25, -3.5||0.6 0.2|
|9200||74.8||374.0||3.005-3.25, -3.5||0.8 0.2|
At the final two epochs, the dust grain size distribution was created by increasing all grain radii by a constant amount, which means that the distribution of grain radii is proportional to (a-3.0) rather than .
3.9 Day 9200
We constructed a model at day 9200 after the explosion, to match Herschel photometry obtained by Matsuura et al. (submitted) on days 9090-9122 and ALMA 440 and 870 m photometry from days 9273-9387 (Indebetouw et al. 2014). We estimated the synchrotron flux at submm wavelengths as before, by extrapolating the observations of Zanardo et al. (2010) to 9200 days, and then deriving a power law from those fluxes. We find that the synchrotron flux in Jy can be approximated by the expression F(=0.411 Jy; we add the predicted synchrotron flux from this expression to our model fluxes for dust emission to obtain the total flux to be compared to the observations. Once again, we restricted the volume filling factor of the models to the value of 0.1 required at early epochs, and we obtained a fit to observed fluxes by varying the dust mass, luminosity, and grain size distribution.
We find that a slightly larger dust mass is required to fit the day 9200 observations than was required at day 8515. Models with dust masses between 0.6 and 1.0M can fit the observed data when the luminosity and grain size distribution are varied, and as before, the lower the dust mass, the higher the luminosity and the larger the dust grains required to fit the SED. A good fit to the observations can be obtained with grains 3 m in radius, a luminosity of 360 L and a mass of 0.8 M; fits can be obtained with dust masses up to 0.2 M above and below this value by varying the luminosity and the grain size distribution. However, for dust masses as low as 0.6 M, the luminosity required is 525L, which is higher than the luminosities required for any model on day 8200 and is thus probably unrealistic. Figure 14 shows the predicted SEDs for three models.
3.9.1 Clump size
For all the models discussed so far, all parts of the ejecta expanded uniformly between epochs. We also considered models in which the clump radius is held constant, while the clumps themselves move outwards with the rest of the ejecta. This may be a more realistic scenario, as in the case of uniform expansion the outer edge of the clumps would be moving at some 140 km s relative to the centre of the clump, and after 23 years of expansion the density might be too low for continuing dust formation. For the models with non-expanding clumps, a match to the observed SED requires a higher dust mass and a higher luminosity compared to the models in which the ejecta expands uniformly. The greater density of the clumps results in cooler dust, and while an increase in grain radii is still necessary to fit the SED, the increase required is less than the case of uniformly expanding ejecta. Figure 15 shows the predicted SEDs for two models at day 9200: for expanding clumps, the model contains 0.8 M of dust, with grain sizes from 3.005–3.25m, and the luminosity is 360 L, while in the case of non-expanding clumps, a similar fit to the SED is obtained with 1.0 M of dust, with grain sizes from 1.005–1.25m, and a luminosity of 900 L. Such a high luminosity is difficult to account for: it would require that in addition to the 400 L from the decay of Ti, all of the 500 L X-ray luminosity from ejecta-CSM interaction goes into heating the dust. This suggests that the clumps must be expanding but at a lower velocity than in the case of uniform expansion. The dust masses derived assuming uniform expansion are thus likely to be somewhat underestimated.
By constructing realistic three-dimensional models of the dusty ejecta of SN 1987A, we have provided new constraints on the location and timing of dust formation in the supernova ejecta. We find that the very large dust mass deduced by Matsuura et al. (2011) from the far-infrared and submm emission detected by Herschel is physically plausible, and that the dust geometry determined from day 615 observations, expanded at constant velocity and with its total dust mass increased, then also gives good fits to the subsequent epochs. At days 615, 775 and 1153, the observational data is sufficient to place strong upper limits on the mass of dust formed, and we thus show that the overwhelming majority of the dust mass present by day 8515 must have formed after day 1153. The dust masses we derive at each epoch are summarised in Table 3.
Very few supernovae have been observed at the mid-infrared and longer wavelengths necessary to quantify dust production at epochs later than 1000 days after explosion. In addition, at such late epochs, optical and infrared fluxes are commonly elevated by scattering and thermal echoes of the supernova flash from pre-existing circumstellar or interstellar material, which may dominate the observed fluxes at late times. Otsuka et al. (2012) presented optical to mid-IR photometry of 6 supernovae at epochs 500 days after explosion, and found evidence of light echoes in 5 of the sample supernovae. Given that we have shown that in the case of SN 1987A, at least 99% of the dust must have formed after day 1153, the much lower masses of dust found at mid-IR wavelengths for most core collapse supernovae observed within the first 1000 days clearly cannot rule out the later formation of a much greater mass of dust.
Sarangi & Cherchneff (2013) predicted that most dust formation in the ejecta of core collapse supernovae must occur at relatively early times (5 years after the explosion), while for SN 1987A we have shown that the vast majority of the dust must have formed more than 3.2 years after the explosion. When the current large quantity of dust formed, and over how long, is not well constrained beyond the fact that it must have been after day 1153. However, we find that the evolution of the dust mass with time can be approximated by a sigmoid function of the form
From a least squares minimisation, we obtained coefficients of a=1.0 M (representing the maximum dust mass), b=-8.53 and c=-0.000366 to fit the dust masses estimated from our modelling. The dust masses derived at each epoch, together with the sigmoid function, are shown in Figure 16. Such a function may be physically reasonable, with a slow increase in dust mass at first due to marginal fulfilment of the necessary physical conditions for dust formation, and a slow increase in the later stages of dust formation due to depletion of refractory elements, with the highest rate of dust formation occurring at an intermediate epoch. If the dust formation can realistically be approximated by a function such as this, then even by 3000 days after the explosion, the dust mass was only 5% of its final value, and 90% of the dust would form between days 3000 and 14000.
We estimated the dust emission expected 7 and 10 years after the explosion by expanding the grid as for the other epochs, and using the dust mass estimated from the sigmoid function. 2500 days after the explosion, the expected dust mass is 0.03 M, while 3650 days after the explosion, the function predicts a dust mass of 0.11 M. We ran models for luminosities of 750 and 1500 L, and for three grain size distributions, in which grain radii were increased by 0, 0.05 and 0.5 m from the originally assumed grain size distribution. We find that at day 2500, the SED peaks at around 100 m with fluxes ranging from 0.05–0.2 Jy, while at day 3650, the predicted peak fluxes are 0.2–0.45 Jy.
We have found that to obtain dust cool enough to match the observed SED at day 8515, large grains are required, with radii of up to several microns. Considering the standard grain growth equation for the time taken for a grain to grow by accretion to a radius of from an initial radius :
then for a time of 21.5 years between days 615 and 8515, a grain bulk density of 3.3 gcm, a sticking coefficient of 0.5, particles of atomic weight 12, and a gas temperature of 20K, the mean velocity of accreting particles is 2.310cm s, and the relation between the number density of accreting particles and the growth in grain radius in microns is given by:
For grains to grow in radius by 1 m in 21 years, the carbon atom number density required would be 3.610 cm. This would go up or down by a factor of ten for a final grain size an order of magnitude larger or smaller respectively. At day 615, the clumps in our model have a radius of 8.310cm and a volume of 2.410cm; in the models with the best-fitting volume filling factor of 0.1 there are 12368 such clumps. Allowing for other lighter species in the metal-rich regions of the ejecta, the total gas density would be (1-3)10 cm. An order of magnitude estimate for the total gas mass required in the clumps would then be 22.4 M. This exceeds the progenitor mass of 19 M, again suggesting that the large grains present in the ejecta at very late times cannot have formed solely by accretion onto existing grains but must also have formed by coagulation.
Larger grains are more likely to survive the passage of the reverse shock through the ejecta; Nozawa et al. (2007) find that grains larger than 0.2 m in size would survive into the ISM essentially unchanged. In our best fitting day 8515 and day 9200 models, all the grains are larger than this, and the majority of the dust could thus potentially survive the reverse shock. A significant fraction of the dust can survive even in the models with smaller maximum radii, which fit the SED less well; for a model with grain radii in the range 0.05 a 0.5 m and n(a)a, for example, 2.8% of the grains are larger than 0.2 m, but they contain 9.6% of the mass. This suggests, therefore, that a significant fraction of the dust formed in the ejecta of SN 1987A will survive the passage of the reverse shock.
Recently, Gall et al. (2014) have found evidence for large grains (0.7 m) being formed in SN2010jl. They find evidence for large grains within a short time of the supernova explosion, in contrast to our finding for SN1987A that the larger grains present at day 8500 could not have been present before day 1000. Nevertheless, in both cases, the presence of large dust grains that will probably survive the passage of the reverse shock through the ejecta supports the hypothesis that core-collapse supernova can be a major source of dust in galaxies.
One of the long standing questions regarding SN 1987A is the nature of the compact remnant. Despite extensive searches, no evidence for a pulsar has been found at any wavelength, despite strong theoretical predictions that one should have formed in the explosion (McCray 1993; Arnett et al. 1989). Graves et al. (2005) found an upper limit to the integrated flux between 0.29–0.965 m of F1.610 erg scm, and used this limit to rule out scenarios including accretion onto a neutron star or black hole, based on an assumption of relatively low dust absorption in the remnant, with A=0.595. The large amount of dust now known to be present means that the optical extinction could be much higher than previously assumed. If one or more dense clumps of dust were to lie along the line of sight, the upper limit on the optical flux of the compact remnant could be much higher than the value found by Graves et al. (2005). In a day 8515 model with 0.6M, each clump has a mass of 4.8510M, and an edge-to-edge optical depth of 12.2, equivalent to A=13.3; assuming that the dust clumps expand with , the optical depth will decline with 1/. If a single clump were to lie directly along the line of sight, then with no dust formation or destruction, or alteration of the dust geometry apart from expansion, the lower extinction assumed by Graves et al. (2005) would not be reached for approximately another 85 years. In fact, several clumps are likely to lie along the line of sight; in the geometry adopted, the clump covering factor is almost unity. From the point of view of the central source, the average number of clumps intercepted by an arbitrary line of sight is 3.9, and the chance of a line of sight not intercepting any clumps is 0.1%. We find that over all lines of sight, the median value of A is 4.8, and it exceeds 15 in approximately 5% of directions. The vast dust reservoir of SN 1987A may thus keep any compact remnant hidden from view at optical wavelengths for a very long time.
We thank Paul Harvey and Allan Meyer for providing further information about the day 1179 and 1183 KAO 50-m observations of SN 1987A. This work was co-funded under the Marie Curie Actions of the European Commission (FP7-COFUND).
- Andrews et al. (2010) Andrews J. E., Gallagher J. S., Clayton G. C., Sugerman B. E. K., Chatelain J. P., Clem J., Welch D. L., Barlow M. J., Ercolano B., Fabbri J., Wesson R., Meixner M., 2010, ApJ, 715, 541
- Andrews et al. (2011) Andrews J. E., Sugerman B. E. K., Clayton G. C., Gallagher J. S., Barlow M. J., Clem J., Ercolano B., Fabbri J., Meixner M., Otsuka M., Welch D. L., Wesson R., 2011, ApJ, 731, 47
- Arnett et al. (1989) Arnett D., Fryxell B., Mueller E., 1989, ApJ, 341, L63
- Arnett et al. (1989) Arnett W. D., Bahcall J. N., Kirshner R. P., Woosley S. E., 1989, ARA&A, 27, 629
- Barlow et al. (2010) Barlow M. J., Krause O., Swinyard B. M., Sibthorpe B., Besel M.-A., Wesson R., Ivison R. J., Dunne L., Gear W. K., Gomez H. L., Hargrave P. C., Henning T., Leeks S. J., Lim T. L., Olofsson G., Polehampton E. T., 2010, A&A, 518, L138
- Basko (1994) Basko M., 1994, ApJ, 425, 264
- Bertoldi et al. (2003) Bertoldi F., Cox P., Neri R., Carilli C. L., Walter F., Omont A., Beelen A., Henkel C., Fan X., Strauss M. A., Menten K. M., 2003, A&A, 409, L47
- Biermann et al. (1992) Biermann P. L., Chini R., Haslam C. G. T., Kreysa E., Lemke R., Sievers A., 1992, A&A, 255, L5
- Biermann et al. (1990) Biermann P. L., Chini R., Haslam G., Kreysa E., Lemke R., 1990, A&A, 236, L17
- Bouchet & Danziger (1993) Bouchet P., Danziger I. J., 1993, A&A, 273, 451
- Cernuschi et al. (1967) Cernuschi F., Marsicano F., Codina S., 1967, Annales d’Astrophysique, 30, 1039
- Cherchneff & Dwek (2010) Cherchneff I., Dwek E., 2010, ApJ, 713, 1
- Draine & Lee (1984) Draine B. T., Lee H. M., 1984, ApJ, 285, 89
- Dwek et al. (2007) Dwek E., Galliano F., Jones A. P., 2007, ApJ, 662, 927
- Dwek et al. (1992) Dwek E., Moseley S. H., Glaccum W., Graham J. R., Loewenstein R. F., Silverberg R. F., Smith R. K., 1992, ApJ, 389, L21
- Ercolano et al. (2005) Ercolano B., Barlow M. J., Storey P. J., 2005, MNRAS, 362, 1038
- Ercolano et al. (2003) Ercolano B., Barlow M. J., Storey P. J., Liu X.-W., 2003, MNRAS, 340, 1136
- Ercolano et al. (2007a) Ercolano B., Barlow M. J., Sugerman B. E. K., 2007a, MNRAS, 375, 753
- Ercolano et al. (2007b) Ercolano B., Barlow M. J., Sugerman B. E. K., 2007b, MNRAS, 379, 1248
- Ercolano et al. (2008) Ercolano B., Young P. R., Drake J. J., Raymond J. C., 2008, ApJS, 175, 534
- Fransson et al. (2013) Fransson C., Larsson J., Spyromilio J., Chevalier R., Gröningsson P., Jerkstrand A., Leibundgut B., McCray R., Challis P., Kirshner R. P., Kjaer K., Lundqvist P., Sollerman J., 2013, ApJ, 768, 88
- Gallagher et al. (2012) Gallagher J. S., Sugerman B. E. K., Clayton G. C., Andrews J. E., Clem J., Barlow M. J., Ercolano B., Fabbri J., Otsuka M., Wesson R., Meixner M., 2012, ApJ, 753, 109
- Gomez et al. (2012) Gomez H. L., Krause O., Barlow M. J., Swinyard B. M., Owen P. J., Clark C. J. R., Matsuura M., Gomez E. L., Rho J., Besel M.-A., Bouwman J., Gear W. K., Henning T., Ivison R. J., Polehampton E. T., Sibthorpe B., 2012, ApJ, 760, 96
- Graves et al. (2005) Graves G. J. M., Challis P. M., Chevalier R. A., Crotts A., Filippenko A. V., Fransson C., Garnavich P., Kirshner R. P., 2005, ApJ, 629, 944
- Hanner (1988) Hanner M., 1988, Technical report, Grain optical properties
- Harvey et al. (1989) Harvey P., Lester D., Dinerstein H., Smith B., Colome C., 1989, in Bulletin of the American Astronomical Society Vol. 21 of Bulletin of the American Astronomical Society, Far-IR Monitoring of SN1987A. p. 1215
- Herant & Woosley (1994) Herant M., Woosley S. E., 1994, ApJ, 425, 814
- Indebetouw et al. (2014) Indebetouw R., Matsuura M., Dwek E., Zanardo G., Barlow M. J., Baes M., Bouchet P., Burrows D. N., et al. 2014, ApJ, 782, L2
- Jerkstrand et al. (2011) Jerkstrand A., Fransson C., Kozma C., 2011, A&A, 530, A45
- Kamenetzky et al. (2013) Kamenetzky J., McCray R., Indebetouw R., Barlow M. J., Matsuura M., Baes M., Blommaert J. A. D. L., Bolatto A., et al. 2013, ApJ, 773, L34
- Lucy et al. (1989) Lucy L. B., Danziger I. J., Gouiffes C., Bouchet P., 1989, in Tenorio-Tagle G., Moles M., Melnick J., eds, IAU Colloq. 120: Structure and Dynamics of the Interstellar Medium Vol. 350 of Lecture Notes in Physics, Berlin Springer Verlag, Dust Condensation in the Ejecta of SN 1987 A. p. 164
- Lucy et al. (1991) Lucy L. B., Danziger I. J., Gouiffes C., Bouchet P., 1991, in Woosley S. E., ed., Supernovae Dust Condensation in the Ejecta of Supernova 1987A - Part Two. p. 82
- Matsuura et al. (2011) Matsuura M., Dwek E., Meixner M., Otsuka M., Babler B., Barlow M. J., et al. 2011, Science, 333, 1258
- McCray (1993) McCray R., 1993, ARA&A, 31, 175
- Meikle et al. (1993) Meikle W. P. S., Spyromilio J., Allen D. A., Varani G.-F., Cumming R. J., 1993, MNRAS, 261, 535
- Morgan & Edmunds (2003) Morgan H. L., Edmunds M. G., 2003, MNRAS, 343, 427
- Moseley et al. (1989) Moseley S. H., Dwek E., Glaccum W., Graham J. R., Loewenstein R. F., 1989, Nature, 340, 697
- Nozawa et al. (2007) Nozawa T., Kozasa T., Habe A., Dwek E., Umeda H., Tominaga N., Maeda K., Nomoto K., 2007, ApJ, 666, 955
- Pietrzyński et al. (2013) Pietrzyński G., Graczyk D., Gieren W., Thompson I. B., Pilecki B., Udalski A., Soszyński I., Kozłowski S., et al. 2013, Nature, 495, 76
- Pozzo et al. (2004) Pozzo M., Meikle W. P. S., Fassia A., Geballe T., Lundqvist P., Chugai N. N., Sollerman J., 2004, MNRAS, 352, 457
- Sarangi & Cherchneff (2013) Sarangi A., Cherchneff I., 2013, ApJ, 776, 107
- Smith et al. (2008) Smith N., Foley R. J., Filippenko A. V., 2008, ApJ, 680, 568
- Sturm et al. (2010) Sturm R., Haberl F., Aschenbach B., Hasinger G., 2010, A&A, 515, A5
- Suntzeff (1997) Suntzeff N. B., 1997, ArXiv Astrophysics e-prints
- Suntzeff & Bouchet (1990) Suntzeff N. B., Bouchet P., 1990, AJ, 99, 650
- Suntzeff et al. (1991) Suntzeff N. B., Phillips M. M., Depoy D. L., Elias J. H., Walker A. R., 1991, AJ, 102, 1118
- Thielemann et al. (1990) Thielemann F.-K., Hashimoto M.-A., Nomoto K., 1990, ApJ, 349, 222
- Todini & Ferrara (2001) Todini P., Ferrara A., 2001, MNRAS, 325, 726
- Walker & Suntzeff (1991) Walker A. R., Suntzeff N. B., 1991, PASP, 103, 958
- Wooden (1989) Wooden D. H., 1989, PhD thesis, California Univ., Santa Cruz.
- Wooden et al. (1993) Wooden D. H., Rank D. M., Bregman J. D., Witteborn F. C., Tielens A. G. G. M., Cohen M., Pinto P. A., Axelrod T. S., 1993, ApJS, 88, 477
- Woosley (1988) Woosley S. E., 1988, ApJ, 330, 218
- Woosley et al. (1988) Woosley S. E., Pinto P. A., Weaver T. A., 1988, Proceedings of the Astronomical Society of Australia, 7, 355
- Zanardo et al. (2010) Zanardo G., Staveley-Smith L., Ball L., Gaensler B. M., Kesteven M. J., Manchester R. N., Ng C.-Y., Tzioumis A. K., Potter T. M., 2010, ApJ, 710, 1515
- Zubko et al. (1996) Zubko V. G., Mennella V., Colangeli L., Bussoletti E., 1996, MNRAS, 282, 1321