Research Article

Massive-Star Supernovae as Major Dust Factories

See allHide authors and affiliations

Science  14 Jul 2006:
Vol. 313, Issue 5784, pp. 196-200
DOI: 10.1126/science.1128131


We present late-time optical and mid-infrared observations of the Type II supernova 2003gd in the galaxy NGC 628. Mid-infrared excesses consistent with cooling dust in the ejecta are observed 499 to 678 days after outburst and are accompanied by increasing optical extinction and growing asymmetries in the emission-line profiles. Radiative-transfer models show that up to 0.02 solar masses of dust has formed within the ejecta, beginning as early as 250 days after outburst. These observations show that dust formation in supernova ejecta can be efficient and that massive-star supernovae could have been major dust producers throughout the history of the universe.

Millimeter observations of high-redshift (z > 6) quasars have revealed the presence of copious amounts of dust when the universe was as young as 700 million years (1). At the present day, dust in the interstellar medium of the Milky Way and other galaxies is generally thought to be injected mainly by the gentle winds of low-mass stars when they evolve onto the Asymptotic Giant Branch (2). However, stellar-evolution time scales of these low- to intermediate-mass stars are too long for them to have been a major contributor to the dust budget in the early universe (3). Instead, dust in the early universe must reflect the contribution from rapidly evolving (1 to 10 million years) massive stars, which return their nuclear ashes in explosive Type II supernova (SN) events. Theoretical studies have long suggested that dust can condense in the ejecta from core collapse (e.g., Type II) SNe (4), and calculations predict condensation of 0.08 to 1 solar mass (M) of dust within a few years, depending on metallicity and progenitor mass (57). There is also evidence for the origin of some dust in Type II SNe on the basis of the isotopic composition of stardust isolated in meteorites (8).

Direct observational evidence for efficient dust formation in SN ejecta is, however, lacking, largely because SN explosions are rare and far apart. Dust formation was detected in the ejecta of SNe 1987A and 1999em, but only some 10–4 M were inferred for each (911), a factor up to 103 smaller than typical SNe would have had to produce in order to contribute efficiently to the early-universe dust budget (12). Similarly low dust masses have been measured in evolving SN remnants with the recently launched Spitzer Space Telescope (13); however, its mid-infrared (mid-IR) instruments are most sensitive to warm (50 to 500 K) material, whereas dust in these remnants has almost certainly cooled to <30 K. Cold dust has been detected in remnants in the far-IR and sub-mm wavelengths; however, such observations risk strong contamination by cold, unrelated dust clouds along the line of sight (14). As such, the best way to demonstrate dust condensation in SN ejecta is to study them within a few years of their explosion, during the epoch of condensation when the ejecta are much hotter than interstellar dust. The high sensitivity of Spitzer's mid-IR detectors allows us to sample very young core collapse SNe within ∼20 Mpc and opens up the whole nearby universe for such studies. Here we report on such a study of the Type II-P SN 2003gd in the galaxy NGC 628. A rare combination of contemporaneous optical and mid-IR observations of this well-studied SN with a known stellar progenitor mass of Math (15, 16) provides an excellent test case for the efficiency of dust formation in SN ejecta.

Data in support of dust production. As dust condenses in SN ejecta, it increases the internal optical depth of the expanding ejecta, producing three observable phenomena: (i) a mid-IR excess; (ii) asymmetric blue-shifted emission lines, because the dust obscures more emission from receding gas; and (iii) an increase in optical extinction. All of these were observed from 1 to 3 years after outburst in SN 1987A (9, 10). In this section, we present or confirm all three phenomena from SN 2003gd.

NGC 628 was observed by the Spitzer Infrared Nearby Galaxies Survey (SINGS) Legacy program (17) with Spitzer Space Telescope's Infrared Array Camera (IRAC) at 3.6, 4.5, 5.8, and 8.0 μm on 28 July 2004, or day 499 after outburst (18), and with the Multiband Imaging Spectrometer for Spitzer (MIPS) on 23 January 2005 (day 678) at 24 μm; the SN was also observed with IRAC as part of GO-3248 (P. I. W. P. Meikle) on 2005 Jan 15 (day 670). All data were acquired from the Spitzer archive, then spatially enhanced with the SINGS data pipelines to a final resolution of 0.″75 pixel–1. A point source identified in all four IRAC bands from day 499 is consistent to within 0.″17 (0.2 IRAC pixels) with the position of the SN progenitor (Fig. 1), as measured by careful absolute and differential astrometry between the 3.6-μm IRAC image and the archival Hubble Space Telescope (HST) Wide Field and Planetary Camera 2 (WFPC2) data in which the progenitor was identified (16).

Fig. 1.

Hubble and Spitzer Space Telescope images of a 10″by 10″ field centered on the position of SN 2003gd. (A) HST WFPC2 image taken in the F606W filter in which the SN progenitor (arrow) was identified (15, 18), with a resolution of 0.″045 pixel–1. (B and C) False-color SST IRAC images of the SN taken on (B) 28 July 2004 and (C) 15 January 2005, showing the 3.6-μm (blue), 4.5-μm (green), and 8.0-μm images (red). All IRAC images were processed by the SINGS collaboration to a final resolution of 0.″75 pixel–1.

Photometry of the Spitzer data was performed with point-spread-function fitting techniques (19). The resulting flux densities are 20.8 ± 2.6 μJy at 3.6 μm, 73.8 ± 5.6 μJy at 4.5 μm, 64.9 ± 7.3 μJy at 5.8 μm, and 103 ± 22 μJy at 8.0 μm on day 499, and 106 ± 16 μJy at 24 μm on day 678; the SN is nearly undetectable in IRAC on day 670, with 3σ upper limits of 6.0, 10.6, 13.5, and 26.6 μJy at 3.6, 4.5, 5.8, and 8.0 μm, respectively (Fig. 2). An excess at 4.5 μm could be due to CO emitting at 4.6 μm (10). Otherwise, the rising 5.8- to 24-μm flux densities are not expected from the gaseous ejecta, which typically have temperatures of 3000 to 5000 K at these late times (10). It is unlikely that this emission is a thermal light echo, because the time variability of the IRAC fluxes is much faster than that expected for a typical circumstellar dust shell (20). Assuming that the mid-IR emission can be modeled with a single blackbody, the best fit to the 5.8- to 8.0-μmdata from day 499 has a temperature of 480 K, which for an adopted distance to the SN of 9.3 Mpc (18) yields an integrated luminosity of 4.6 × 105 L and an equivalent radius of 6.8 × 1015 cm. The size, temperature, and variability implied by the spectral-energy distributions (SEDs) are thus consistent with lower temperature (T ∼ 500 K) dust that is cooling within the SN ejecta.

Fig. 2.

Spectral-energy distribution of SN 2003gd, showing optical photometry on day 493 (squares) from (18) and extrapolated from day 632 (21) to day 678 (diamonds) based on the evolution of SN 1987A (24); IRAC data from day 499 (circles); upper limits to IRAC from day 670 (arrows); and the MIPS datum from day 678 (triangle). Error bars are computed with a Poisson-noise model that includes detector characteristics, and flat-field and profile uncertainties. Fluxes have been dereddened by E(BV) = 0.14 (18). The curves are MOCASSIN radiative-transfer model fits to the data at day 499 (solid lines) and 678 (dashed lines) based on smoothly distributed (black) and clumpy (gray) dust. See text and Table 2.

The first indication of such dust formation in SN 2003gd came in a comparison of broad-band photometry and Hα spectra of the SN between days 124 and 493 (18), in which a small decline in the late-time luminosity was accompanied by a slight blueshift in the emission-line peaks. New spectroscopic observations of SN 2003gd were obtained in long-slit mode, covering ∼4500 to 7000 Å with a spectral resolution of ∼7 Å, by the Gemini Multi-Object Spectrograph (GMOS) on Gemini North on 19 August 2004 (day 521). Two spectra of SN 2003gd were obtained, which have been debiased, flattened, wavelength-calibrated, sky-subtracted, extracted, and then combined. The wavelength calibration is accurate to about 2 Å. The Hα and [O I] spectra from days 157 and 493 (18) are compared to the newer data in Fig. 3. Inspection of the lines confirms that the emission peak is indeed blueshifting, whereas the red high-velocity wing seen at the earliest epoch has diminished. Additionally, the most recent spectra show a clear asymmetry in the first few hundred km s–1 redward of each line peak. These blueshifted peaks and asymmetric profiles are consistent with a simple model in which dust with an increasing optical depth is located within an expanding sphere of uniform emission (11); in particular, the day 521 profiles suggest an optical extinction AR < 5.

Fig. 3.

Optical spectra of SN 2003gd showing [O I] λλ6300,6363 Å (top) and Hα (middle and bottom). Spectra have been corrected for the redshift z = 0.00219 of galaxy NGC 628. The gray curves are taken from the data presented in (18), showing the profiles at days 157 (solid line) and 493 (broken line). The solid black curve is from the new data taken on day 521 with GMOS-N on Gemini North. The profiles are normalized to an arbitrary flux scale. A monotonic blueshift in the Hα line peak, first reported by (18), is confirmed. The most recent spectra in all three lines also show a clear profile asymmetry redward of each line peak. This evolution is expected from dust forming homogeneously within the ejecta, which preferentially extinguishes emission from the receding (i.e., redshifted) gas.

The B, V, and R light curves of SN 2003gd through day 493 (18) have been combined with the B- and R-band photometry from day 632 (21) in Fig. 4. Over this 500-day period, the SN is evolving almost identically to SN 1987A, with the exception that SN 2003gd is slightly subluminous in V and R around day 493. The close match between the light curves of these two SNe implies an increase in optical extinction of SN 2003gd as well.

Fig. 4.

Light curves of SN 2003gd, showing the increase of extinction with time. B, V, and R light curves of SN 2003gd compiled from (18) through day 493 after outburst, and (21) for day 632, are plotted and offset as marked. Error bars are all smaller than the point markers. For comparison, the corresponding light curves of SN 1987A (24) are also shown as thick gray lines, and the R-band light curve of SN 1990E is shown as filled stars. Also plotted in black are the light curves expected from Eq. 1 with (dashed) and without (dotted) the effective opacity term, as well as the linear fit to SN 1987A from days 450 to 525 (dashed gray) used by (9).

The observed dust extinction is measured by comparing the change in photometry over time to a standard intrinsic luminosity. The internal extinction in SN 1987A is believed to have increased by AR = 0.8 mag between days 525 and 700, as determined in broad-band photometry by comparing the light curve after day 525 to the best fit line through the data from days 450 to 525 (9). The slope of this line (Fig. 4) is extremely sensitive to the subset of the SN 1987A light curve used, and the resulting dust extinction can vary by >1 mag depending on the days included in the least-squares fit. Thus, we deem this a poor measure of the intrinsic luminosity.

The broad-band evolution of Type II SNe past ∼500 days is poorly documented, because very few SNe have been observed beyond this epoch. At these late times, the light curve is dominated by the energy input of γ-rays from 56Co decay, which decreases with an e-folding time of τ56 = 111.2 days. The R-band photometry of the Type II SN 1990E (22) closely follows this evolution through 540 days (Fig. 4), suggesting that simple Co decay provides a good estimate of the unextinguished R-band light curve for at least that long. However, as the ejecta expand, their opacity to γ-rays is expected to decrease, which results in a modified light curve (23) Math(1) where the term in brackets is the “effective opacity”; κ56 = 0.033 cm2 g–1 is the average opacity to 56Co-decay γ-rays, and ϕ0 = 7 × 104 g cm–2 is the column depth at the fiducial time t0 = 11.6 days, chosen to match the bolometric light curve of SN 1987A (23, 24). Note that Eq. 1 begins fading relative to simple 56Co decay around 500 days.

The Co-decay curves, both with and without the effective-opacity correction, are much more luminous than the aforementioned linear fit to SN 1987A at the time of dust formation (Fig. 4). The slopes of these Co-decay curves also closely resemble that of SN 1987A after day 775, when dust production is believed to have ended (10). Thus, Eq. 1 offers a more realistic standard luminosity for the unextinguished R-band light curve of SN 1987A. Comparison of L56 to the SN 1987A photometry yields extinctions of 1.5 mags around day 700 when effective opacity is included, and 2.5 mag when it is excluded. Because the evolution of SNe 1987A and 2003gd are so similar, Eq. 1 is also used to estimate the extinction of SN 2003gd at each epoch, as listed in Table 1.

Table 1.

R-band extinction of SN 2003gd. Equation 1 is used with and without the effective opacity term to estimate the average and maximal extinction, respectively. The observed values are listed first, and these were used along with the R-band light curve of SN 1987A to extrapolate the extinction of the SN at the epochs of SST observation.

DayAR (mags)
493 0.52 0.73 0.09
632 1.36 1.78 0.12
499 0.53 0.74 0.14
678 1.51 2.13 0.20

Dust-mass analysis and interpretation. The observations presented above overwhelmingly point toward dust forming within the ejecta of SN 2003gd, beginning sometime between 250 and 493 days after outburst. To estimate the mass of dust present, we use the three-dimensional Monte Carlo radiative-transfer code MOCASSIN (25). Briefly, the paths of photon absorption, scattering, and escape are followed from a specified source through a given composition, grain-size distribution, and geometry of dust. The particular choices of these are either constrained a priori or are varied until the model emission and extinction match the observed values.

Because our hypothesis is that dust condenses within the ejecta, the radiative-transfer model is constructed under the initial assumption that the dust and source luminosity are mixed within a spherical, expanding shell with inner radius rin, outer radius Y·rin, and ρ ∝ r–2 density profile, and with the illuminating radiation proportional to the dust density. Initial values for the shell size, source luminosity, and temperature are guided by the blackbody previously fit to the mid-IR data and by models of SN 1987A at similar epochs (10). There are numerous models for dust formation within SN ejecta [for a review, see (26)], most of which predict that grain sizes will remain small. We adopt a standard a–3.5 size distribution (27), with grain radii between 0.005 and 0.05 μm, and the dust composition is taken to be 15% amorphous carbon and 85% silicates (7), with optical constants taken from (28, 29). Finally, because there are very few optical and mid-IR data to constrain a given model, the source luminosity is restricted to evolve according to Eq. 1, whereas its temperature remains constant (10).

Two dust distributions are considered. In the first, “smooth” model, the dust is uniformly distributed throughout the shell according to the adopted density profile. However, as early as a few hours after outburst, postshock ejecta become Rayleigh-Taylor unstable (30, 31), forming an inhomogeneous or “clumpy” distribution, which we model as a two-phase medium, in which spherical clumps with size rc = δ·(Y rin), volume filling factor fc, and density contrast α = ρc/ρ are embedded within an interclump medium of density ρ. This is analogous to the mega-grains approximation (32) with the addition of a radial density profile. Only macroscopic mixing has been found in the clumpy ejecta of SN remnant Cas A (33), which suggests that elemental ejecta layers remain heterogeneous. We therefore assume that the source luminosity is completely separated from the dust clumps. For a given geometry, a clumpy model will always require more mass than a smooth one to fit a given SED, because clumping lowers the overall optical depth for a given mass of dust (32). Rather than explore the extensive parameter space of clumpy models, we study the limiting case where all dust is in clumps, i.e., α→∞, which should provide upper mass limits, whereas the smooth models will provide lower mass limits. Finally, as suggested from hydrodynamic simulations (31), we fix δ = 0.025.

Model results are summarized in Fig. 2 and Table 2. A good fit to the day 499 photometric and extinction data was achieved for the smooth model with Y = 7, rin = 5 × 1015 cm, L = 6.6 × 105 L, and T = 5000 K, whereas fitting the day 678 data required changing rin to 6.8 × 1015 cm and L to 9.2 × 104 L. Clumpy models used these same parameters, with fc = 0.02 on day 499, and fc = 0.05 on day 678. A complete exploration of the model parameter space is beyond the scope of this work and will be presented elsewhere. In general, small changes to the model parameters have only modest effects. For example, including maximum grain sizes up to 0.25 μm (typical of dust in the interstellar medium) decreases the dust mass by less than 10%. A 10% change in δ or fc results in a 1 to 5% change in mass for our adopted parameter ranges. Thus, the smooth and clumpy model results shown in Table 2 offer reasonably robust lower and upper mass limits, respectively.

Table 2.

Dust masses and R-band extinction calculated by the radiative-transfer code MOCASSIN.

DayModel A R Mdust (M)
499 Smooth 0.40 2.0 × 10-4
499 Clumpy 0.65 1.7 × 10-3
678 Smooth 1.48 2.7 × 10-3
678 Clumpy 1.22 2.0 × 10-2

These clumpy-model masses, up to 2 10–3 M on day 499 and 2 × 10–2 M on day 678, are much higher than most analytic estimates of the dust mass for SN 2003gd. For example, when the same grain properties as those above are used, 5 × 10–4 and 2 × 10–3 M of dust are required to produce the mid-IR emission at days 499 and 678, if all grains are visible and isothermal (34). When we use the mega-grains approximation for dust uniformly mixed with diffuse emission within a spherical shell (32), the R-band extinction yields masses of only 10–5 and 4 × 10–4 M of smooth dust for days 499 and 678. In contrast, up to 5 × 10–3 M of clumpy dust is deduced from the mega-grains model for day 499, which agrees well with our radiative-transfer model. However, once clumps become optically thick, only geometry (δ, fc) determines their extinction, and thus an arbitrarily large mass of clumpy dust reproduces the extinction from day 678. This behavior of the mega-grains model makes it of limited use in determining dust masses when the observed optical depth reaches unity. We conclude that the most-often used analytic approximations can provide unreliable estimates of dust masses.

Observations similar to those presented here have demonstrated the condensation of dust in the ejecta of SN 1987A and 1999em, but the inferred masses for these SNe were only modest, on the order of 10–4 M (10, 11). Asymmetric Hα line profiles have been detected for Type II SNe 1970G, 1979C, and 1980K (35); however, this phenomenon on its own is also consistent with an expanding ionization front shell of emission catching up with and passing through a preexisting dust shell and does not prove the presence of newly forming dust. For SN 1998S, asymmetric line profiles and near-IR (λ < 4.7 μm) excess emission also point toward dust condensation in the ejecta (36). In all these cases, the absence of contemporaneous mid-IR observations precluded a quantitative estimate of the condensation efficiency. In light of our quantitative clumpy-dust analysis, the amounts of dust believed to have formed in the ejecta of SN 1987A and 1999em—which were based upon analytical estimates—are being carefully revisited.

For a progenitor mass for SN 2003gd between 10 and 12 M (18), roughly 0.16 to 0.42 M of refractory elements are expected to form if the progenitor had solar metallicity (37). Assuming dust formation has finished by day 678, our derived dust mass of 2 × 10–2 M for this progenitor translates into a condensation efficiency, defined here as (mass of refractory elements condensed into dust)/(mass of refractory elements in ejecta), of ≤0.12 for clumpy dust. Having found that analytic analyses of optical and IR observations tend to underestimate the dust mass by an order of magnitude or more, we deem it likely that dust formation in core-collapse SNe is much more efficient than previously believed. In particular, the efficiency implied by SN 2003gd is close to the value of 0.2 needed for SNe to account for the dust content of high-redshift galaxies (12). As noted earlier, too few SNe have been followed sufficiently in time and wave length to establish the frequency with which SNe form dust. We are currently addressing this question through continued, long-term monitoring of a much larger sample of young, Type II SNe. If dust formation is found to be common in core-collapse SNe, then because most dust is expected to survive later passage through high-velocity SN shocks (38), we can conclude that core-collapse SNe played an important role in the production of dust in the early universe.

References and Notes

View Abstract

Stay Connected to Science

Navigate This Article