1. Introduction
Multiwavelength studies of SN 1987A, located at a distance of 51.4 ± 1.2 kpc in the Large Magellanic Cloud (Panagia1999), have provided unprecedented details of how supernova (SN) explosions trigger the dynamical distribution of gas in a supernova remnant (SNR) and how this SN/SNR system evolves over time. The morphology of SN 1987A is well studied (see the recent review in McCray & Fransson2016), with the system consisting of ejecta and a bright and distinct equatorial ring (hereafter the ring), together with two fainter outer rings. The ring is composed of circumstellar material that radiates in UV, optical, X-rays, and radio over an extent of 1
6 (0.3 pc) (e.g., Burrows et al.2000; Sonneborn et al.1998; Ng et al.2013), as well as thermal dust emission due to shock heating of preexisting dust formed during the red supergiant phase (Bouchet et al.2006; Dwek et al.2010). The ejecta have a complex morphology. The Hα emission, originating from warm gas irradiated by X-rays from the ring (Larsson et al.2011; Fransson et al.2013), exhibits an elongated north–south structure and a “hole” in the center. Along with hydrogen lines from the ejecta, near-infrared (NIR) emission from warm (∼2000 K) CO and mid-infrared (MIR) emission from SiO in the SN ejecta were detected early (as early as 112 days) after the explosion (e.g., Spyromilio et al.1988; Roche et al.1991). After day 9000, cold expanding CO, SiO, and HCO+ molecules were detected in the submillimeter part of the spectrum (Kamenetzky et al.2013; Matsuura et al.2017), highlighting that a significant part of the ejecta is cold (13–132 K). Interestingly, the inner ejecta of SN 1987A have not yet mixed with the circumstellar medium (CSM) or interstellar medium (ISM), and the majority has not yet passed through the reverse shock (France et al.2010; Frank et al.2016). Thus, this young SNR is an ideal source for studying the footprints of the gas dynamics since the very early days of the SN, as the gas has been assumed to be freely expanding since its explosion (McCray1993). Indeed, recent high angular resolution emission-line images of SiO and CO (Abellán et al.2017) from the Atacama Large Millimeter/submillimeter Array (ALMA) have been used to compare the distribution of the molecular gas ejecta with the predictions from models of the gas dynamics after the SN (Wongwathanarat et al.2015, M. Gabler et al. 2019, in preparation).
The progenitor of SN 1987A, Sanduleak −69° 202, was a blue supergiant (Gilmozzi et al.1987; Kirshner et al.1987; West et al.1987; White & Malin1987), thought to have had a zero-age main-sequence mass of ∼19M⊙ (Hashimoto et al.1989; Woosley et al.1997), with a mass of ∼14M⊙ at the time of the explosion (Woosley1988; Smartt et al.2009; Sukhbold et al.2016). From its mass, the expectation is that a neutron star should have formed at the time of explosion. Despite prompt neutrino emission observed at the burst (Hirata et al.1987) indicating the formation of a neutron star (Burrows1988; Sukhbold et al.2016), the search for a compact object associated with SN 1987A has been difficult: observational searches have proven unfruitful (e.g., Manchester2007; Alp et al.2018a; Esposito et al.2018; Zhang et al.2018). The possible detection of radio polarization toward the ejecta (Zanardo et al.2018) hints at the presence of magnetized shocks, potentially due to a compact object. Alp et al. (2018a) proposed that a thermally emitting neutron star could be dust obscured and that this may be detectable as a point source in far-infrared (FIR) or submillimeter images of the remnant, though this has not yet been detected.
It is still largely debated whether or not SNe are net dust producers or destroyers in galaxies (e.g., Morgan & Edmunds2003; Matsuura et al.2009; Gall et al.2011; Gomez2013; Dwek et al.2014; Rowlands et al.2014; Lakićević et al.2015; Michałowski et al.2015; Temim et al.2015; Watson et al.2015). Due to its youth and proximity, SN 1987A is an excellent laboratory for studying SN dust. It is also rare, since any dust emission seen in the inner region of the remnant can be attributed unambiguously to dust formed in the SN ejecta and not from the swept-up CSM/ISM or unrelated foreground/background material (a common issue with Galactic SNRs, e.g., Morgan et al.2003; Gomez et al.2012a; De Looze et al.2017; Chawner et al.2019). SN 1987A also provides insight into dust formation at an early stage compared to previously studied Galactic SNRs—here we can probe timescales on the order of decades rather than centuries, filling in the gap between very young SNe (e.g., Gall et al.2014) and historical remnants.
Thermal emission from small quantities (10−4M⊙) of dust was detected in the early days after the SN explosion (day ∼ 300–600) using MIR observations (Danziger et al.1989; Lucy et al.1989; Bouchet et al.1991; Roche et al.1993; Wooden et al.1993). More surprisingly, theHerschel Space Observatory (Pilbratt et al.2010; hereafterHerschel) revealed a large amount of cold dust (∼0.5 M⊙) at the location of the remnant (Matsuura et al.2011,2015). ALMA resolved the emission from dust in SN 1987A on scales of 03 and confirmed that the ∼0.5M⊙ of cold (20 K) dust discovered withHerschel originates from the inner ejecta region (Indebetouw et al.2014; Matsuura et al.2015). Dwek & Arendt (2015) and Wesson et al. (2015) revisited the dust emission at early times (<1200 days). Dwek & Arendt (2015) find that a large mass of dust can be present early on (0.4
at ∼615 days) with a model of silicates and amorphous carbon. Wesson et al. (2015) conclude instead, from comparing radiative transfer models to the optical–IR spectral energy distribution (SED) limits, that the dust mass increased more slowly over the first 10 yr. This substantial mass of dust observed in the inner debris of SN 1987A demonstrates that a large fraction of the heavy elements ejected in an SN may be locked up in a dust reservoir.
Whether dust grains formed in the ejecta of an SN are carbon- or silicate-rich remains an unanswered question: the models of Cherchneff & Dwek (2009,2010) and Sarangi & Cherchneff (2013,2015) predict that for abundance ratios C/O < 1, carbon atoms will mostly be locked up in CO molecules in the first 1000 days, preventing the formation of a large mass of amorphous carbon dust. Though CO may be dissociated by electrons produced by radioactive decay (Clayton2011) and/or (to a lesser extent) X-rays from the ring, depending on the gas density, the models of Sarangi & Cherchneff (2013,2015) indicate that the dissociation of CO is insignificant. In contrast, in order to explain the FIR dust emission, a substantial fraction of the dust grains must be composed of amorphous carbon (amC; Matsuura et al.2015), as the emissivity of amC grains is higher than that of silicates in general, thus leaving an unresolved tension between observations and theory. We note that a model that explains the FIR emission and requires only a small amount of amC grains, with the majority of mass in silicates, was proposed by Dwek & Arendt (2015). There they fit the FIR SED with amC–silicate composite grains assuming a “continuous distribution of ellipsoids” (CDE) model and found a reduced dust mass, though the majority of the reduction in mass in this case arises from the inclusion of dust grains with long axial ratios (so-called needles), which allows the CDE model to surpass the FIR emissivity of amorphous carbon. No evidence of the silicate signature was found in the warm dust emission in the first 2 yr after the explosion, suggesting that small silicate grains were not the first condensates (Roche et al.1993).
In this work, we present high angular resolution ALMA (Cycle 2) dust images for SN 1987A, where we resolve dust clumps on scales of ∼80 mas. Here we revisit the dust mass and grain composition using the ALMA photometry. We discuss the implications of our results for the gas-phase chemistry leading to dust formation and find evidence for warmer gas at the center of the inner ejecta, hinting at the possible indirect detection of a compact source.
2. Data
2.1. Observations and Reduction
Our observations of SN 1987A were obtained with ALMA, as part of the Cycle 2 observing program 2013.1.00063.S. The data were taken over several days in the latter half of 2015, between 10,352 and 10,441 days after the initial explosion. The Band 7 (870μm) and 9 (450μm) integrations utilized between 34 and 36 antennae, with baselines spanning 15 m–2.3 km. See Table1 for a summary of the observations.
Table 1. Observations
| Sub- | Frequency | Baselines | Angular Scales | Observation | SN | Bandpass | Phase | Check | Time | |
|---|---|---|---|---|---|---|---|---|---|---|
| band | Range (GHz) | (m) | (arcsec) | Date | Day | Calibrator | Calibrator | Source | (minutes) | |
| B7 | A1 | 299.88–315.87 | 45.4–1574.4 | 0.13–4.31 | 2015 Jun 28 | 10,352 | J0538–4405 | J0635–7516 | J0601–7036 | 18.4 |
| A2 | 299.88–315.87 | 43.3–2269.9 | 0.09–4.52 | 2015 Sep 22 | 10,438 | J0538–4405 | J0635–7516 | J0601–7036 | 18.3 | |
| B | 342.48–358.34 | 15.1–1574.4 | 0.11–11.46 | 2015 Jul 25 | 10,379 | J0538–4405 | J0635–7516 | J0601–7036 | 20.9 | |
| C | 346.23–362.09 | 15.1–1574.4 | 0.11–11.34 | 2015 Jul 25 | 10,379 | J0538–4405 | J0635–7516 | J0601–7036 | 19.9 | |
| D1 | 303.62–319.48 | 45.4–1574.4 | 0.13–4.26 | 2015 Jun 28 | 10,352 | J0538–4405 | J0635–7516 | J0601–7036 | 18.8 | |
| D2 | 303.62–319.48 | 43.3–2269.9 | 0.09–4.47 | 2015 Sep 22 | 10,438 | J0538–4405 | J0635–7516 | J0601–7036 | 18.8 | |
| B9 | A | 673.44–681.06 | 43.3–2269.9 | 0.04–2.10 | 2015 Sep 25 | 10,441 | J0522–3627 | J0601–7036 | J0700–6610 | 12.5 |
| B | 680.94–688.56 | 43.3–2269.9 | 0.04–2.07 | 2015 Sep 25 | 10,441 | J0522–3627 | J0601–7036 | J0700–6610 | 12.5 | |
| C | 688.44–696.06 | 43.3–2269.9 | 0.04–2.05 | 2015 Sep 25 | 10,441 | J0522–3627 | J0700–6610 | J0450–8101 | 12.5 |
Note.
Observations for proposal ID 2013.1.00063. Each sub-band is composed of four 2 GHz blocks of 128 channels (15.625 MHz each). The same flux calibrator, J0519–454, was used for all observation blocks.
Download table as: ASCIITypeset image
Each data set was reduced separately with the Common Astronomy Software Applications package (Casa;21 McMullin et al.2007), version 4.5.1. Once calibrated, thetclean algorithm was used to deconvolve and image the data.
The check source (reference quasar with precisely known position) and phase calibrator coordinates, determined withimfit inCasa, were offset by no more than 0.4 mas from the catalog values. Other measures of astrometric quality for our observing configurations include the ALMA baseline measurement accuracy (2 mas), noise-limited signal error ∼(beam size)/(S/N) (3–4 mas), and the phase transfer error from the measured phase rms (<12 mas), where S/N is the signal-to-noise ratio and rms denotes root mean square. Combining these, the overall astrometric accuracy we assume for the data presented in this work is 10 mas in Band 6, 12 mas in Band 7, and 15 mas in Band 9.
Decorrelation due to factors such as weather was investigated using the flux calibrator, phase calibrator, and check source by phase-averaging over several intervals and integrating the resulting flux densities—a large variation in flux density for different phase-averaging intervals would suggest that decorrelation is pronounced enough to decrease the recovered flux. The variations of the calibrator flux densities in all Band 7 and Band 9 windows were within the systematic uncertainties except for Band-7C (346–362 GHz), which had significantly worse weather than the other segments, with an estimated decorrelation of ∼35%. Bands 7A and 7D also suffered from poor weather in the original 2015 June observations, with ∼1.45 millimeter precipitable water vapor (PWV), and were therefore repeated in 2015 September. These are denoted as 7A2 and 7D2. Despite the poorer quality of the June data, combining them with the September data results in higher-S/N images.
Self-calibration, a common technique for high-S/N data where calibrating the data against itself in successive deconvolution cycles can often result in improved dynamic range, was determined to have a negligible impact on the images. Final images were cleaned with natural weighting applied to the baselines to optimize sensitivity per beam. The imaging parameters, including resolution and sensitivity, are given in Table2.
Table 2. Imaging Properties
| Frequency Range | νc | Beam FWHM | Beam PA | rms Noise |
|---|---|---|---|---|
| (GHz) | (GHz) | (arcsec2) | (deg) | (mJy bm−1) |
| Continuum | ||||
| 224.00–227.00 | 225.50 | 0.30 × 0.30 | 0.00 | 0.12 |
| 238.00–243.00 | 240.50 | 0.30 × 0.30 | 0.00 | 0.09 |
| 246.00–249.00 | 247.50 | 0.30 × 0.30 | 0.00 | 0.10 |
| 269.50–270.50 | 270.00 | 0.30 × 0.30 | 0.00 | 0.18 |
| 278.00–280.00 | 279.00 | 0.30 × 0.30 | 0.00 | 0.10 |
| 306.06–307.47 | 306.76 | 0.20 × 0.15 | 124.32 | 0.07 |
| 311.88–319.48 | 315.68 | 0.19 × 0.14 | 119.11 | 0.04 |
| 673.45–685.00 | 679.22 | 0.08 × 0.06 | 74.37 | 0.71 |
| Spectral Lines | ||||
COJ = 2 1 | 230.54 | 0.06 × 0.04 | 27.43 | 0.04 |
COJ = 6 5 | 691.47 | 0.09 × 0.07 | 185.35 | 2.82 |
SiOJ = 5 4 | 217.10 | 0.06 × 0.04 | 19.74 | 0.05 |
SiOJ = 6 5 | 260.52 | 0.04 × 0.03 | 173.66 | 0.06 |
SiOJ = 7 6 | 303.93 | 0.13 × 0.10 | 35.47 | 0.47 |
Note.
The position angles are counterclockwise from north. COJ = 2
1, SiOJ = 5
4, and SiOJ = 6
5 parameters are for the data cubes from Abellán et al. (2017). Theνc values listed for the CO and SiO lines are rest frequencies. The rms values for the spectral lines are per velocity channel. For observation dates and epochs, see Table1.
Download table as: ASCIITypeset image
2.2. Defining Continuum Wavelength Ranges
The wavelengths covered by these observations include many spectral lines from molecular species—primarily CO and SiO, with contributions from various SO lines and potentially others. The ±∼1000 km s−1 expansion velocity of the ejecta means that the line widths span a substantial fraction of the observed bands. The continuum bands selected relative to the modeled molecular line emission are shown in Figure1, using the ALMA spectra and the emission-line model of CO, SiO, SO, and SO2 from Matsuura et al. (2017). Only windows that were free from molecular line emission (shown by the gray vertical bands) were used to make continuum images, centered at roughly 307, 315, and 679 GHz. The 315 GHz continuum image is shown in Figure2. We also utilize here the Cycle 2 Band 6 imaging data presented by Matsuura et al. (2017), to provide continuum information below 300 GHz. The Band 6 images were restored to a common circular beam with FWHM of 030.
Figure 1. Spectra and integrated continuum flux densities of the ejecta and ring for Bands 6, 7, and 9. The molecular line emission model is taken from Matsuura et al. (2017). Vertical gray bands indicate the portions of the spectrum deemed to be line-free and therefore used in creating the continuum images at 307, 315, and 679 GHz (the 315 GHZ image is shown in Figure2). Data points indicate the flux densities in that band for the ring (green) and the ejecta (yellow).
Download figure:
Standard imageHigh-resolution image
Figure 2. ALMA 315 GHz (with beam) and 2014HST F625W band image (Fransson et al.2015), which includes Hα. The yellow contours display 315 GHz emission at 0.2 mJy beam–1. The 315 GHz continuum in the inner ejecta originates from thermal dust emission, while in the ring it is due to synchrotron emission. The 18 mas uncertainty on the relative alignment due to Band 7 astrometric error (12 mas) andHST image registration based on fitting the ring (6 mas) is of order 1 pixel in these images.
Download figure:
Standard imageHigh-resolution image2.3. Molecular Line Data
In this section we present the molecular line data observed in the same blocks as the continuum discussed above: COJ = 6
5 with rest frequency 691.47 GHz and SiOJ = 7
6 atνrest=303.93 GHz. The 345.80 GHz COJ = 3
2 and 347.33 GHz SiOJ = 8
7 lines were also covered in these observing blocks, but as they are heavily blended, we do not consider them in the present work.
The SiOJ = 7
6 and COJ = 6
5 cubes were created withtclean inCasa, with a spectral resolution of 300 km s−1, which gives a reasonable balance between velocity resolution and sensitivity per channel. COJ = 6
5 was imaged with natural weighting to maximize sensitivity per beam. SiOJ = 7
6 was imaged withrobust= −1 in order to better spatially resolve the central features of interest. A comparison of the integrated (moment-0) maps of the CO and SiO lines is given in Figure3.
Figure 3. Multiwavelength view of the ejecta in SN 1987A. The cyan contours represent the 679 GHz dust emission at 3σ and 5σ. Beam sizes for individual maps are denoted by the green ellipses. The small plus sign denotes the system center as defined in AppendixB. The bottom right panel is a three-color image of COJ = 2
1 in red, SiOJ = 5
4 in green, and SiOJ = 7
6 in blue and highlights how varied the spectrally integrated emission is between the various line transitions. The brightest areas are generally distinct patches of primary color instead of blended, demonstrating that the CO and SiO peaks are not cospatial, and none match the distribution of the 679 GHz dust. Some of the CO falls in the Hα hole (the lower left), but the majority of the CO peaks on the little Hα “wing” to the right of the hole. The small 5σ cyan contour just northeast of the center of the remnant is the so-called “blob.”
Download figure:
Standard imageHigh-resolution imageIn addition to the molecular line data described above, we also utilize the COJ = 2
1, SiOJ = 5
4, and SiOJ = 6
522 data as described in Abellán et al. (2017). Although both sets were taken in Cycle 2, the molecular line data presented by Abellán et al. (2017) have higher S/N, as they were combined with Cycle 3 data, and due to COJ = 6
5 being in a band with poorer atmospheric transmission than COJ = 2
1. The additional Cycle 3 data also give their COJ = 2
1, SiOJ = 5
4, and SiOJ = 6
5 maps finer spatial resolution than the observations presented in the current work, with FWHM 004–006 (see Table2). For full details of their data reduction technique, we refer the reader to their Section 2. The channel maps are shown in AppendixA. The given velocities are the observed values, not shifted to the reference frame of SN 1987A. The systemic velocity (kinematic local standard of rest frame; LSRK) of SN 1987A is 287 km s−1 receding from Earth (Gröningsson et al.2008).
3. Description of Images
The SN 1987A ring and ejecta continuum image at 315 GHz is shown in comparison to Hα in Figure2. These images have been aligned following a technique in Alp et al. (2018a) where the ring emission is used to derive a reference center, though here we take a simpler approach (see AppendixB for details). Our derived ring+ejecta system center used in this work isα = 5h35m27
998,δ = −69°16′11107 (International Celestial Reference System; ICRS), ±18 mas (Figure18). At 315 GHz, the ring is clumpy and the brightness contrast in the east and west components of the ring is different from that observed in the Hα ring emission. The brighter emission observed in the NE and SW regions of the ring in the radio is similar to that seen in hard X-rays (Helder et al.2013; Frank et al.2016). The ejecta are located at the center of the image inside of the ring structure. The ring emission at 315 GHz is attributed to synchrotron (see Section4.3), and the inner region is thermal dust emission from the SN ejecta (Indebetouw et al.2014).
Figure3 shows an enlarged view of the ejecta images of dust continuum and lines. The majority of the submillimeter ejecta continuum is distributed in a roughly symmetrical ellipsoid, with fainter asymmetrical emission protruding west and southwest. At 315 GHz, the ejecta are moderately resolved and show a conspicuously separate clump of emission south of the main body of the ejecta. This clump persists in images produced with lowerrobust settings intclean, where sensitivity is lower and spatial resolution is higher. Both the primary ejecta material and the smaller clump as observed in the 315 GHz image appear to fill in the gaps seen in the Hα image, like a “lock in the keyhole.” This is shown in Figure2, where the 3σ contours highlighting the major continuum features are overlaid onto the continuum and Hα images. The alignment accuracy is ∼1 pixel in the images, given the astrometric uncertainties discussed in Section2.1 (12 mas for Band 7 continuum) and AppendixB (6 mas for registration of theHubble Space Telescope [HST] image to the 315 GHz ALMA image).
The 679 GHz image provides the highest-resolution view of the continuum (top left panel of Figure3). This figure shows that the dust is asymmetrically distributed and is composed of several clumps, with the brightest feature (hereafter the “blob”) just northeast of the center of the remnant. The beam resolution provides a limit on the clump size—assuming a distance of 51.4 ± 1.2 kpc (Panagia1999), the Band 9 beam FWHM of 008 × 006 corresponds to a physical scale of 0.020 × 0.016 pc, or 4125 × 3230 au. Nevertheless, the resolved 679 GHz image indicates that dust is not smoothly distributed across the ejecta, and the locations of dust clumps are not identical to clumps in the CO or SiO. The S/N in the 679 GHz image is moderate—the outer cyan contours in Figure3 and the dust emission (in red) in Figure4 have pixel S/N > 3, and the surrounding ejecta area has pixel S/N values of ∼2 in the 679 GHz image. The area between the ejecta and the ring—outside of the outermost ejecta contour in Figure2—is consistent with noise.
Figure 4. Spatial anticorrelation of dust and COJ = 6
5. In this overlay, COJ = 6
5 is in blue and 679 GHz is in red, showing their relative spatial distributions. The off-white plus sign denotes the system center position as defined in AppendixB. The gold line highlights the extent of the major features in the dust, while the teal line demarcates the major COJ = 6
5 features. The contrast of each component image was set independently to emphasize its major features—visible blue and red features roughly correspond to S/N > 3 for the dust and S/N > 5 for COJ = 6
5, respectively. The dust and COJ = 6
5 emission exhibit a notable anticorrelation.
Download figure:
Standard imageHigh-resolution imageThe molecular images provide a probe of different conditions in the ejecta, where lower transitions probe lower-temperature gas (if optically thin, see Section5). One prominent feature is the central hole seen in the COJ = 2
1 and SiOJ = 5
4 images (middle left and bottom left panels of Figure3). This was first reported by Abellán et al. (2017) and was seen in both the integrated 2D spatial maps and the 3D data cubes. Although the integrated SiOJ = 6
5 map (middle panel of Figure3) does not show the hole clearly in the same manner as SiOJ = 5
4 and COJ = 2
1, the hole is also visible in the central channels (v = 0–300 km s−1) of the velocity map (Figure17). Because of the additional −600–0 km s−1 components located within the same line of sight as the hole (Figure17) in the integrated maps, the hole is not clear in the SiOJ = 6
5 map. The CO and SiO molecular hole is just to the south of the “keyhole” that is seen in Hα (Figure 8 of Fransson et al.2015; top right panel of our Figure3), though the molecular hole appears to be slightly smaller in scale and located on the southern edge of the hole in Hα emission. The centers of the holes are offset by ∼50 mas, or ∼4× the astrometric and alignment errors.
COJ = 2
1 and SiOJ = 5
4 have similar structures in the integrated images; however, the spatial distributions of the higher transitions of each species have some differences. SiOJ = 6
5 is more evenly distributed in a shell pattern, while the lower-S/N image of COJ = 6
5 appears clumpy (Figure3), though this is likely affected by the noise.
COJ = 6
5 has emission coincident with the COJ = 2
1 hole, in that its channel maps (Figure16) show emission around the hole location, albeit at low S/N. However, the integrated spatial distribution appears different from COJ = 2
1. The brightness peaks are distributed differently, and the hole is not visible in the integrated COJ = 6
5 map owing to some emission at those coordinates in the 600–900 km s−1 channels (the far side). The presence of a molecular hole in SiOJ = 7
6 cannot be confirmed in these data, as the systemic line center (vLSRK ∼ 300 km s−1) falls at the edges of two sidebands observed separately, which were concatenated during reduction, and suffers from roll-off at the edge of the spectral window; the resulting S/N is poor in that channel. The other molecular lines do not share this limitation, as they fell well within the sideband spectral windows. We do note a peak of SiOJ = 7
6 emission, however—the brightest source of emission in the entire cube—overlapping with the spatial location of the hole and the dust blob but offset from the systemic velocity by ∼−400 km s−1 (this corresponds to the 0 km s−1 channel of Figure17).
The resolved dust peak (small 5σ contour in Figure3) is colocated with the molecular hole in the low transitions of CO and SiO and slightly extends to the north and east into the relative depression visible in the SiOJ = 5
4 channels near the systemic velocity. The brightest points of dust emission tend to coincide with relative depressions in the COJ = 6
5 brightness, giving the appearance of an anticorrelation between the main dust and COJ = 6
5 features. This is more clearly demonstrated in Figure4, where the dust (red) and COJ = 6
5 (blue) images are overlaid. The individual images were normalized independently to emphasize the main features of each, with the visible colors shown roughly corresponding to areas of S/N > 3. The gold and teal lines are guides highlighting the highest-S/N features of the dust and CO, respectively, in order to compare peaks in the emission. While there is some overlap in the faint features of the dust and CO, and the southern extent of the dust peak starts to fade to a combined magenta, the dust and COJ = 6
5 peaks do not generally overlap. Rather, the brightest dust features are located in areas of relatively faint COJ = 6
5 emission and vice versa.
To test whether the apparent dust–CO anticorrelation is an artifact or result of the data reduction or continuum subtraction, we performed several checks. First, the Band 9 dust continuum was reconstructed in different ways, by imaging (Casamfs-mode) in a variety of spectral windows and also by making a data cube across the entirety of Band 9 (including the CO line) and fitting the continuum emission. These techniques gave consistent results. Initially we usedCasa to subtract a (zeroth-order) continuum in the visibility plane. We compared this with an order-0 subtraction in the image plane and found no significant differences. This means that the structures seen in our final COJ = 6
5 map are robust to variations in how the continuum is determined and subtracted. The anticorrelation in COJ = 6
5 is visible, even before continuum subtraction, in the COJ = 6
5 dirty map (i.e., with no cleaning to deconvolve the interferometer sidelobes). Thus, the apparent anticorrelation seen in the dust and CO distributions is robust to changes in the data processing. Lastly, we test whether the anticorrelation is statistically robust by calculating the weighted version of the normalized cross-correlation function
, which returns a standard correlation measurer between the range of −1 and +1. The pixel weights used were the map (S/N)2. The resulting correlation measure isr = −0.30 ± 0.08—a moderate anticorrelation—using the accuracy estimate from, e.g., Frick et al. (2001). Due to the relatively low S/N in these images and therefore scatter in pixel-by-pixel comparisons,r will always be pulled closer to 0 and will not approach ±1.
We tested the robustness of the correlation measure by investigating different angular scales using the wavelet analysis described in Frick et al. (2001) and Arshakian & Ossenkopf (2016). On scales of one to two beamwidths (i.e., convolutions with kernels of those scales), the images start to become increasingly positively correlated, which is expected owing to the peaks being separated by that amount (∼5× the astrometry error). Below these scales,r remains negative, so the anticorrelation is not sensitive to small changes in image resolution. Using the same wavelet analysis on the other images, the dust correlates more positively with the other CO and SiO lines than with COJ = 6
5. The standard correlation measures agree, withr = +0.04 for COJ = 2
1 andr = +0.36 for SiOJ = 6
5.
The brightest dust feature is located one beamwidth northeast of the secondary COJ = 6
5 peak, and the brightest COJ = 6
5 features curve around the main dust peak. The CO peaks are obvious because the line emission is brighter than the continuum (see the integrated profile in Figure1), but there is some fainter (S/N < 5) COJ = 6
5 emission overlapping with some of the dust emission. There is also low-level (S/N <2) dust emission that roughly spans the full extent of the ejecta. The other molecular species further complicate the picture, as noted earlier—the peak of SiOJ = 7
6 is coincident with the dust peak and the SiOJ = 5
4 hole (as seen in the bottom right panel of Figure3). At the southern edge of the hole, some faint Hα emission appears to be aligned with COJ = 2
1 in projection, but their velocity ranges differ. The 3D view of Hα (Larsson et al.2016, their Figure 6) indicates that the Hα emission in this region peaks at velocities around −1500 km s−1, while the peak COJ = 2
1 is between 0 and +200 km s−1. While no velocity information is available for the dust continuum emission, it is spatially offset from the nearby Hα and COJ = 2
1 by approximately one dust resolution element. That is, the COJ = 2
1 and Hα in this region are offset in velocity, and the dust peak is spatially offset from both.
To summarize, we find that the dust emission is clumpy. The Band 9 image enables us to resolve the dust in the ejecta to angular scales of 62 × 81 mas. The peaks of the ejecta CO and SiO emission are not cospatial with the peaks of the ejecta dust emission (with anticorrelated COJ = 6
5 and dust structures). The small peak/clump in the dust emission revealed in the Band 9 image (the blob) overlaps with holes previously observed in the lower line transitions of SiO and the CO molecular ejecta and is coincident with some emission observed in the SiOJ = 7
6 line.
4. The Spectral Energy Distribution of SN 1987A
The three physical mechanisms primarily responsible for emission in the FIR–radio portion of the continuum are thermal IR graybody emission from dust (from the ejecta region; Matsuura et al.2011; Indebetouw et al.2014), nonthermal radio/millimeter synchrotron emission (from the ring; Manchester2007; Potter et al.2009; Zanardo et al.2010,2013; Lakićević et al.2012; Indebetouw et al.2014) and a lesser contribution from free–free millimeter/submillimeter bremsstrahlung emission from hot ionized material (see Zanardo et al.2014 for a full review of the different components). In this section we measure the photometry, analyze the emission from dust in the ejecta using the ALMA data, investigate the properties derived using a variety of dust models from the literature, and investigate the synchrotron emission in the ring.
4.1. Photometry
The continuum bands are defined as those frequencies that are molecular line free, as demonstrated in Figure1 and Section2.2, with the chosen frequency windows summarized in Table3; these bands are different from the default ALMA wide-band continuum. The centers of the apertures used for deriving photometry are the same as described in Section3 (AppendixB), with elliptical apertures selected to encompass the ejecta and the ring annulus with varying sizes in each ALMA band in order to include only the relevant signal (Table3). For the 315 GHz ejecta, this results in an elliptical aperture with semimajor and semiminor axes ofrMAJ ×rMIN = 042 × 036 and major-axis P.A. of 25° (N through E). For the 315 GHz circumstellar ring, an elliptical annulus with innerrMAJ andrMIN of 045 × 042, outerrMAJ andrMIN of 145 × 135, and 85° P.A. from N was chosen. The extents of the regions were selected independently across the different bands to best match the features in each image, but they only vary slightly. For comparison with previous lower spatial resolution observations, the total system emission is also calculated by summing emission within an ellipse defined by the outer ring extent above. This total system integration includes the contribution from the gap between the ring and the ejecta.
Table 3. Continuum Photometry
| νc | Δν | Sν | Aperture Center | RMAJa | RMINa | P.A. | |
|---|---|---|---|---|---|---|---|
| (GHz) | (GHz) | (mJy) | R.A. (deg) | Decl. (deg) | (arcsec) | (arcsec) | (deg) |
| Ejecta | |||||||
| 225.50 | 3.0 | 0.5 ± 0.2 | 83.866586 | −69.269733 | 0.32 | 0.29 | 125 |
| 240.50 | 5.0 | 0.7 ± 0.2 | 83.866600 | −69.269720 | 0.37 | 0.36 | 75 |
| 247.50 | 3.0 | 1.1 ± 0.4 | 83.866639 | −69.269739 | 0.41 | 0.35 | 125 |
| 270.00 | 1.0 | 1.2 ± 0.6 | 83.866652 | −69.269709 | 0.40 | 0.31 | 60 |
| 279.00 | 2.0 | 1.3 ± 0.4 | 83.866662 | −69.269747 | 0.41 | 0.40 | 125 |
| 306.76 | 1.4 | 1.8 ± 0.6 | 83.866667 | −69.269753 | 0.42 | 0.36 | 125 |
| 315.68 | 7.6 | 1.8 ± 0.6 | 83.866667 | −69.269753 | 0.42 | 0.36 | 125 |
| 679.22 | 11.5 | 36.2 ± 7.2 | 83.866728 | −69.269740 | 0.42 | 0.36 | 125 |
| Ring | |||||||
| 225.50 | 3.0 | 13.0 ± 0.5 | 83.866585 | −69.269731 | 1.45, 0.35 | 1.35, 0.33 | 175 |
| 240.50 | 5.0 | 12.5 ± 0.5 | 83.866600 | −69.269722 | 1.45, 0.42 | 1.35, 0.39 | 175 |
| 247.50 | 3.0 | 13.4 ± 0.8 | 83.866630 | −69.269736 | 1.45, 0.43 | 1.35, 0.40 | 175 |
| 270.00 | 1.0 | 13.7 ± 1.1 | 83.866592 | −69.269720 | 1.45, 0.45 | 1.35, 0.42 | 175 |
| 279.00 | 2.0 | 12.9 ± 0.8 | 83.866661 | −69.269747 | 1.45, 0.44 | 1.35, 0.41 | 175 |
| 306.76 | 1.4 | 12.0 ± 1.5 | 83.866666 | −69.269753 | 1.45, 0.45 | 1.35, 0.42 | 175 |
| 315.68 | 7.6 | 11.3 ± 1.7 | 83.866667 | −69.269753 | 1.45, 0.45 | 1.35, 0.42 | 175 |
| 679.22 | 11.5 | 19.4 ± 19.3 | 83.866720 | −69.269737 | 1.45, 0.45 | 1.35, 0.42 | 175 |
| Total System | |||||||
| 225.50 | 3.0 | 13.5 ± 0.5 | 83.866585 | −69.269731 | 1.45 | 1.35 | 175 |
| 240.50 | 5.0 | 13.2 ± 0.5 | 83.866600 | −69.269722 | 1.45 | 1.35 | 175 |
| 247.50 | 3.0 | 14.6 ± 0.7 | 83.866630 | −69.269736 | 1.45 | 1.35 | 175 |
| 270.00 | 1.0 | 15.4 ± 0.8 | 83.866592 | −69.269720 | 1.45 | 1.35 | 175 |
| 279.00 | 2.0 | 14.2 ± 0.8 | 83.866661 | −69.269747 | 1.45 | 1.35 | 175 |
| 306.76 | 1.4 | 13.9 ± 1.7 | 83.866666 | −69.269753 | 1.45 | 1.35 | 175 |
| 315.68 | 7.6 | 13.3 ± 1.8 | 83.866667 | −69.269753 | 1.45 | 1.35 | 175 |
| 679.22 | 11.5 | 55.5 ± 20.9 | 83.866720 | −69.269737 | 1.45 | 1.35 | 175 |
Notes. Integrated flux densities of the ejecta, ring, and total system for each continuum band with central frequencyνc and bandwidth Δνc. Integrated flux densities are quoted as value ± measurement uncertainty ± systematic uncertainty. The systematic uncertainty includes calibration uncertainties (10% in Bands 6 and 7 or 20% in Band 9), and an additional 50% on the positive side in Bands 6 and 7 due to the systematic offset from Cycle 0 and Cycle 6 levels. Apertures for the ejecta and total system are ellipses; apertures for the ring are annuli. The center coordinates are in ICRS.
aThe ring annulus radii are given as (Router,Rinner).Download table as: ASCIITypeset image
The Cycle 2 images, aside from Band 9, exhibit a slight decrease in integrated flux density for radii just beyond the outer edge of the ring, which could be due to undersampled flux or errors in calibration or deconvolution. We take the rms of flux densities in background pixels from a large annulus beyond the ring to estimate the level of these effects, and the resulting uncertainty contribution is typically of order a few percent of the flux density.
As the reconstruction of the images may propagate systematic as well as random noise in the background, we have used our images to make a series of measurements using the same aperture as for the source. The distribution of these measurements is roughly Gaussian, and thus we adopt the rms of this distribution as our aperture error,σAP.
An additional empirical uncertainty component was included to account for the potential smearing of ring emission into the ejecta aperture. This was estimated by taking the average deviation in the flux density after expanding and shrinking the semimajor and semiminor axes by 01, a size that covers reasonable large differences in aperture choice yet avoids significant overlap between the ejecta and ring.
All of these uncertainties were added in quadrature to estimate the overall uncertainty in a given band (typically ∼15%, dominated by the random-position aperture uncertainty). We include an additional 10% uncertainty for systematic (calibration) error in Bands 6 and 7, and 20% in Band 9.23 Finally, we considered the possibility that we may be missing diffuse emission from cold dust within the SN structure as a result of overresolving an extended source. To address this issue, we simulated observations for multiple synthetic sources resembling SN 1987A but with varying extended ellipse components and found that this effect is at a level below the ALMA systematic uncertainties. The reader is referred to AppendixC for more details.
The ejecta, ring, and total system flux densities and uncertainties are shown in Figure5 as gold circles, green rings, and purple diamonds, respectively, and their values are listed in Table3. Previous measurements are also shown in Figure5 for reference. Preliminary Cycle 6 flux densities (M. Matsuura et al. 2019, in preparation) from 11,522 days after the explosion are included here. The ejecta flux density is 1.6 mJy at 252.4 GHz and 1.7 mJy at 254.3 GHz. The total system flux densities at these frequencies are 17.9 and 18.1 mJy, respectively. The uncertainty on each of these flux density measurements is estimated as 0.4 mJy.
Figure 5. Continuum values from 225 to 679 GHz for the integrated ejecta (yellow), ring (green), and total system (purple) from this work, for observations taken an average of 10,402 days after the SN explosion. For reference, ALMA Cycle 0 corresponds to day 9280, Cycle 2 is day 10,402, and Cycle 6 is day 11,522. The higher angular resolution ALMA data confirm that the ejecta emission follows a thermal dust profile down to ∼200 GHz. The ring emission, on the other hand, shows no evidence of submillimeter dust emission but instead is consistent with a synchrotron emission profile—ATCA observations from day 9280 to day 9686 (Zanardo et al.2014; Callingham et al.2016) combined with the Cycle 2 ALMA ring flux densities give a power-law index ofα = −0.70 ± 0.06. The dark-blue and cyan lines show SED fits for amorphous carbon (ACAR sample; Zubko et al.1996) and silicate (forsterite; Jäger et al.2003) dust emission models, respectively, demonstrating that disparate models give reasonable fits to the data. The previous ALMA ejecta measurements from Zanardo et al. (2014) are shown as red stars, and previous ALMA ring measurements from Zanardo et al. (2014) are denoted as blue stars. Preliminary Cycle 6 ALMA flux densities (M. Matsuura et al. 2019, in preparation) are shown as the orange plus signs and light-purple hexagons for the ejecta and total emission, respectively. The positive error bars for the Cycle 2 Band 6 and 7 data include an additional 50% of the flux density values to reflect the observed systematic offset from the Cycle 0 and Cycle 6 levels. The open gray hexagons are measures of the total system flux density in various parts of the spectrum at day 9280: 1.4, 18, and 44 GHz (Zanardo et al.2013); 9 GHz (Ng et al.2013); 94 GHz (Lakićević et al.2012); and 102, 213, 345, and 672 GHz (Zanardo et al.2014). The open black hexagons represent the ATCA total system flux densities between 1 and 9 GHz at day 9686 (Callingham et al.2016). The ring flux density at day 9280 (blue stars) is larger than ALMA total emission (purple diamonds) at day 10,402 owing to the ALMA Band 7 systematic error discussed in Section4.1. The brown asterisks are the unresolvedHerschel 70–500μm flux densities (Matsuura et al.2015), and the crimson asterisks are the 2010 HERITAGE flux densities (Matsuura et al.2011; Meixner et al.2013). The ring emission modelS(ν,t) from Cendes et al. (2018) for day 10,402 andα = −0.70 is shown by the dashed gray line.
Download figure:
Standard imageHigh-resolution imageWe note that our Cycle 2 total flux densities in Bands 6 and 7 are systematically lower than the ALMA Cycle 0 and 6 flux densities, though they agree within the error bars. They are typically around 50% lower than the equivalent levels from Cycles 0 and 6; therefore, we include an additional 50% to their positive systematic uncertainties. Potential causes, as noted in Section2.1, include decorrelation from poor weather or a mis-scaling of the flux calibrator. The integrated 350 and 360 GHz ejecta flux densities are particularly low, either due to inherently low flux at this epoch or due to data quality. As mentioned in Section2.1, weather affected the phase stability of observations between 346 and 362 GHz, which can result in decorrelation and therefore reduced flux recovery. As these measurements are less reliable, they have been omitted from the remainder of this study. We note that the systematic offset will not have affected the analysis of the resolved dust distribution discussed in Section3, since the offset is not seen in the Band 9 data, where the dust peak (the blob) was identified.
The literature values for the total SN 1987A flux densities at various wavelengths are shown as gray hexagons and represent the overall SED of the system. The total emission at 1.4, 18, and 44 GHz (Zanardo et al.2013); 9 GHz (Ng et al.2013); 94 GHz (Lakićević et al.2012); and 102 GHz (Zanardo et al.2014) is dominated by synchrotron emission from the ring. The total emission at 213, 345, and 672 GHz (Zanardo et al.2014), on the other hand, gradually consists of a higher and higher fraction of thermal emission until that is dominant in the submillimeter and the FIR. As the synchrotron brightness increases in time (Staveley-Smith et al.2014; Cendes et al.2018), these literature flux densities were scaled to their levels at day 9280 by Zanardo et al. (2014) to match the average epoch of the ALMA cycle 0 observations. The details of the ejecta and ring portions of the SED will be discussed in turn in the following two sections.
4.2. Modified Blackbody Fits to the Ejecta Dust Emission
4.2.1. Description of the Modified Blackbody Fits
Figure5 displays the millimeter to FIR SED. In this figure, the brown asterisks show the FIR flux densities measured byHerschel for the total SN 1987A (unresolved) system, and the gold circles show the millimeter flux densities from this work measured with ALMA for the resolved ejecta. The shape of the SED shows that the ejecta emission arises from thermal (dust) radiation, all the way into the millimeter, confirming the results of Zanardo et al. (2014) and Matsuura et al. (2015).
The next step is to fit the thermal dust emission using dust models. In order to cover the peak of the thermal emission, we use theHerschel flux densities from Matsuura et al. (2015) in our model fits since they are measuring the emission from the ejecta dust, albeit unresolved. TwoHerschel flux densities are treated as upper limits: 70μm, as it is possibly contaminated by warm ring dust (Matsuura et al.2019) and/or [Oi] 63μm emission; and 500μm, as it was a nondetection. One potential issue with using theHerschel flux densities is that theHerschel measurements were obtained at an average of ∼1300 days before the ALMA cycle 2 data, and the FIR emission could potentially vary over time. The heating source of the ejecta was suggested to be primarily from44Ti decay, which has an estimated lifetime of 85 yr (Ahmad et al.2006; Jerkstrand et al.2011; also see later Matsuura et al.2011). The predicted decrease in this decay energy between the 2012Herschel and 2015 ALMA observations is 4.2%. Assuming that the FIR luminosity decreased by this amount between the 2012 and 2015 epochs, the reduction in the temperature of a 20 K blackbody would be ∼0.2 K, translating into individualHerschel flux density decrements of 2%–5%. This is several times smaller than the uncertainties on the PACS and SPIRE flux densities. Therefore, if the ejecta heating is dominated by44Ti decay, the use of theHerschel flux densities with the latest ALMA measurements is valid. An alternative additional heating source will be discussed in Section6.3.
We tested the robustness of the SED results using three common parameter estimation techniques: using a maximum likelihood estimation (MLE) with uncertainties determined by bootstrap resampling; Bayesian estimation with Markov chain Monte Carlo (MCMC) posterior distributions, using theemcee package (Foreman-Mackey et al.2013); and finally by checking with ordinary least-squares (OLS) regression. The OLS, MLE, and MCMC routines all yield consistent fits for a given dust emission profile. In order to take into account the systematic offset in our Band 6 and 7 flux densities from other cycles (see Section4.1), we determine our best fits and uncertainties from resampling of the flux densities within their error bars. The best fit and uncertainties are taken to be the 50th, 16th, and 84th percentiles of the distributions from 1000 samplings.
The modified blackbody (modBB) function we use follows the formSν(λ) = Mdustκabs(λ)Bν(λ)/d2, whereSν(λ) is flux density,κabs(λ) is the mass absorption coefficient of dust grains,Bν(λ) is the Planck function, andd is the distance to SN 1987A. We assume that the emission is optically thin across this wavelength range.κabs(λ) can be directly obtained from the literature for some cases, but for the majority of cases, assuming spherical dust grains of radiusa and densityρ,κabs(λ) is defined asκ = (3/4)Qρa, whereQ is the absorption efficiency of the dust, which can be calculated from optical constants with Mie theory.κabs is often assumed to be a power law defined as
(whereκ0 is the reference value ofκabs(λ) at wavelengthλ0 andβ is the power-law emissivity index ofκabs(λ)).
Fitting a dust graybody to the FIR–millimeter SED with the power-law approximation toκabs gives fit parameters of
,Mdust
, andTd
K ifκ0 (850μm) = 0.07 m2 kg−1 (using the empirical measurement ofκabs assuming that the fraction of metals locked in dust is constant across the local ISM; James et al.2002). As this is a significant amount of dust, here we also fit the SED using a wide variety of compositions and the full characterization ofκabs(λ), following Indebetouw et al. (2014) and Matsuura et al. (2015).
We perform dust graybody fits to the FIR–millimeter SED using a wide variety ofκabs(λ) profiles directly obtained from the literature or calculated using Mie theory; in total we used 134κabs(λ) profiles. These include the following: Weingartner & Draine (2001) LMC average (as an approximation to the conditions near SN 1987A in the LMC), Demyk et al. (2017) amorphous silicate samples at 30 K, Ormel et al. (2011) calculations of bare and icy silicate+graphite grains, Ossenkopf & Henning (1994) bare and icy mantle grains (protostellar core coagulation models), and using Mie scattering calculations for amorphous carbon (amC, Rouleau & Martin1991; Zubko et al.1996), graphite (Draine & Lee1984), “cosmic silicates” (amorphous FeMgSiO4; Jaeger et al.1994), other silicates including enstatite and forsterite from Henning & Stognienko (1996), pure iron (Henning & Stognienko1996), and 68 profiles from the Jena database (Henning et al.1999) of minerals including silicates, amorphous carbons, carbides, oxides, and sulfides. Many of the models in the Jena database only extend to 500 or 1000μm. In these cases we extrapolate to our longest ALMA wavelength of 1329μm in log space (as lines, as most models follow power laws in the submillimeter–millimeter), from the last 300–500μm of each curve to ensure a smooth continuation of the general trend in each. We also consider a “CDE” model for a composite of carbons and silicates (Dwek & Arendt2015), with a carbon volume filling factorfC of 18%.
For Mie theory calculations, we adopt a grain radius ofa = 0.01μm for PAHs anda = 0.1μm for all other models. For most dust models, Mie-derived FIRκabs curves are nearly identical for grain radii ofa < 5μm. A representative sample of 28 dust models was selected from this list to span a wide variety of dust types. Selected grain types include Zubko et al. (1996) amorphous carbons (amCs), Jäger et al. (2003) amorphous silicates, and iron from Henning & Stognienko (1996), among others. The fitted masses and temperatures for this sample are listed in Table4.
Table 4. Modified Blackbody Fits
| Dust | Reference | Grain Density | Mdust | T | κ850 | Good Fit |
|---|---|---|---|---|---|---|
| (g cm−3) | (Mdust ) | (K) | (m2 kg−1) | |||
| amC (ACH2 sample), Mie 0.1μm | Zubko et al. (1996) | 1.81 | 1.46 | 17.5 | 0.087 | Y |
| amC (ACAR sample), Mie 0.1μm | Zubko et al. (1996) | 1.81 | 0.38 | 22.0 | 0.254 | N |
| amC (BE sample), Mie 0.1μm | Zubko et al. (1996) | 1.81 | 0.77 | 20.7 | 0.141 | N |
| amC (AC1 sample), Mie 0.1μm | Rouleau & Martin (1991) | 1.85 | 0.43 | 21.6 | 0.203 | Y |
| Cellulose (800 K sample), Mie 0.1μm | Jager et al. (1998) | 1.81 | 0.46 | 18.8 | 0.178 | Y |
| Graphite, Mie 0.1μm | Draine & Lee (1984) | 2.26 | 1.62 | 17.8 | 0.069 | Y |
| PAH (neutral), Mie 0.01μm | Laor & Draine (1993) | 2.24 | 1.69 | 18.0 | 0.071 | Y |
| Silicate—enstatite, Mie 0.1μm | Jäger et al. (2003) | 2.71 | 4.10 | 18.0 | 0.029 | Y |
| Silicate—forsterite, Mie 0.1μm | Jäger et al. (2003) | 3.2 | 4.03 | 17.9 | 0.029 | Y |
| Silicate—“cosmic,” Mie 0.1μm | Jaeger et al. (1994) | 3.2 | 3.46 | 17.7 | 0.034 | Y |
| Silicate/carbon composite CDE—fC = 0.18 | Dwek & Arendt (2015) | 2.95 | 0.38 | 21.1 | 0.270 | N |
| Silicate—LMC average | Weingartner & Draine (2001) | ⋯ | 2.49 | 18.0 | 0.047 | Y |
| Silicate—30 K average | Demyk et al. (2017) | ⋯ | 0.27 | 18.5 | 0.265 | Y |
| Silicate—composite aggregate | Semenov et al. (2003) | ⋯ | 35.06 | 21.4 | 0.003 | Y |
| Silicate—porous multilayer spheres | Semenov et al. (2003) | ⋯ | 6.13 | 29.4 | 0.010 | N |
| Silicate—bare grains, 0.03 Myr | Ormel et al. (2011) | ⋯ | 0.50 | 20.7 | 0.180 | Y |
| Silicate—icy grains, 0.03 Myr | Ormel et al. (2011) | ⋯ | 0.60 | 18.3 | 0.184 | Y |
Silicate—naked grains, | Ossenkopf & Henning (1994) | ⋯ | 0.64 | 20.9 | 0.163 | N |
Silicate—thin ice mantles, | Ossenkopf & Henning (1994) | ⋯ | 0.74 | 19.1 | 0.142 | Y |
| SiC, Mie 0.1μm | Pegourie (1988) | 3.22 | 1.57 | 23.3 | 0.058 | N |
| FeS, Mie 0.1μm | Henning & Stognienko (1996) | 4.83 | 0.68 | 34.7 | 0.086 | N |
| FeO, Mie 0.1μm | Henning et al. (1995) | 5.7 | 0.28 | 28.2 | 0.259 | N |
| SiO2, Mie 0.1μm | Henning & Mutschke (1997) | 2.196 | 4.41 | 17.5 | 0.022 | N |
| TiO2, Mie 0.1μm | Posch et al. (2003) | 3.78 | 81.77 | 17.7 | 0.001 | Y |
| Al2O3 “compact” sample, Mie 0.1μm | Begemann et al. (1997) | 3.2 | 0.90 | 19.0 | 0.112 | Y |
| NaAlSi2O6, Mie 0.1μm | Mutschke et al. (1998) | 2.4 | 0.10 | 23.1 | 0.982 | Y |
| MgAl2O4, Mie 0.1μm | Fabian et al. (2001) | 3.64 | 121.74 | 17.7 | 0.001 | Y |
| Pure iron, Mie 0.1μm | Henning & Stognienko (1996) | 7.87 | 3.97 | 19.9 | 0.025 | Y |
Note. Mass and temperature fits for graybodies for a selection of 28 of the dust models discussed in the text. amC: amorphous carbon. Fit values and uncertainties are from bootstrap resampling of the data within their error bars and are determined from the 50th, 84th, and 16th percentiles of the distributions of fits from 1000 samplings of the observed flux densities. Forκ(λ) models derived from Mie theory, grains of radiusa = 0.1μm were assumed except for the case of PAHs, wherea = 0.01 μm was used. Grain densities are given for the Mie and CDE cases: amorphous carbon from Zubko et al. (2004) (and Rouleau & Martin1991 for their AC1 sample), graphite and SiC from Laor & Draine (1993), PAHs from Li & Draine (2001), stoichiometric varieties of olivines and pyroxenes from Henning & Stognienko (1996), silicate/carbon composite CDE with carbon volume filling factorfC = 18% from Dwek & Arendt (2015), jadeite from Mutschke et al. (1998), and in general from the Jena optical constants database (https://www.astro.uni-jena.de/Laboratory/OCDB/index.html) (Henning et al.1999).κ850, the mass absorption constant at 850μm, is listed for each model, extrapolated as power laws for models where wavelength coverage falls short. The quality of fit is denoted in the rightmost column, where a good fit is defined as
.
Download table as: ASCIITypeset image
To determine whether the SED fit is “good,” we set a quality threshold whereby
. This on its own can be an insufficient indicator of fit quality, however—for example, a fit that closely matches the majority of the ALMA flux densities can still satisfy
even if it falls below all of theHerschel points owing to the nature of theχ2 metric used in our three fitting methods. We note that these fit criteria are purely formal and do not consider availability of mass from nucleosynthesis yields; those physical limits are discussed in Section6.1.
4.2.2. Results of the Dust Fits
Figure6 shows 28 of the resulting dust fits to the FIR–submillimeter SED (top panel of the figure) using the various dust emissivity curves (bottom panel of the figure) listed in Table4. These have been grouped into similar dust compositions—amC, graphite, PAHs, silicates, oxides, carbides, sulfides, and pure iron—to show the qualitative differences for fits using the various composition types. (Individual models are not labeled; for quantitative measures of specific models, we refer the reader to Table4.) Figure6 illustrates that vastly different dust properties can still translate to relatively similar fits to the observed flux densities. All models that result in good fits to the SED have temperatures between ∼18 and 23 K (Table4); 9 of the 28 dust varieties listed here fail ourχ2 criteria for a “good” fit.
Figure 6. Top: modified blackbody fits to theHerschel and ALMA observations of the ejecta continuum in SN 1987A, using the 28 dust varieties and parameters listed in Table4. TheHerschel 70μm and 500 μm flux densities are used as upper limits. Most dust models give reasonable fits to the observed flux densities. The colors denote the mineralogy of the dust models: amorphous carbon models are shown in red, graphite in orange, PAH in yellow, silicates and silicate composites in blue, oxides (including FeO, TiO2, Al2O3, NaAlSi2O6, and MgAl2O4) in light blue, carbide (SiC) in pink, sulfide (FeS) in green, and pure iron in light purple. Bottom: dustκabs(λ) curves used for the modified blackbody fits, showing the large variation in emission profiles.
Download figure:
Standard imageHigh-resolution imageModels that represent amorphous carbons and graphites tend to have higherκabs values and thus return lower inferred dust masses,Mdust ∼ 0.4–0.8
(Table4). Amorphous silicates from Demyk et al. (2017), which have more than a factor of 10 largerκ850 compared to those from Semenov et al. (2003), also yield a relatively moderate mass of ∼0.27
. The silicate models of Ormel et al. (2011) and Ossenkopf & Henning (1994) also give moderate dust masses and represent grains that originate from coagulation processes as proposed in dense molecular clouds of the ISM. These aggregated grains, some of which are coated with ice, can be as large as 100μm in radius, and the mass absorption coefficients increase in the FIR and millimeter. As a consequence, the inferred dust masses for these grain types are ∼0.5–0.74
for Ormel et al. (2011) and Ossenkopf & Henning (1994). Because the dust grains required for ISM coagulation are calculated to form over timescales of tens of thousands to millions of years, it is unclear whether such icy grains can be formed in the SN ejecta in such a short timescale as <30 yr; however, we present the results of their fits here for comparison. Models of larger composite silicate grains, such as those of Jaeger et al. (1994) or Jäger et al. (2003), require even larger masses of ∼1–4
to fit the observed SED. The CDE composite model of Dwek & Arendt (2015) results in a moderate mass of 0.38
, due to the increased emissivity of the carbon inclusion.
Dust varieties that generally satisfy our criteria for a good fit include several amorphous carbon and silicate models (amorphous pyroxene and olivine varieties), graphite, PAHs, and alumina. Varieties that tend to fit the data poorly in the optically thin limit include FeO, FeS, SiC, SiO2, organics, and water ice.
The largest source of uncertainty in the observed dust mass is the choice of dust emission profile,24 specifically the value of the mass absorption coefficient (κabs) in the submillimeter. We attempted to investigate spatial variations in the fitted dust parameters across the ALMA maps; however, the limiting beam size translates to only two to three independent elements across the ejecta, and we found that the differences in these flux ratios are smaller than their uncertainties.
We compare the inferred dust masses from different dust absorption coefficients (Figure7). The inferred masses very clearly follow an inverse linear relation based on the submillimeterκabs value where![${M}_{\mathrm{dust}}[{M}_{\odot }]=0.117\times {\kappa }_{850}^{-1}$](/image.pl?url=https%3a%2f%2fcontent.cld.iop.org%2fjournals%2f0004-637X%2f886%2f1%2f51%2frevision1%2fapjab4b46ieqn180.gif&f=jpg&w=240)
. The significance of this is that for the current SN 1987A SED, the total ejecta dust mass can be estimated based on a single representative value of the desired dust mass absorption coefficients.
Figure 7. Fitted dust mass as a function ofκ850 from various dust models. The colors for the different dust varieties are the same as in Figure6. The fitted dust mass closely follows a linear inverse trend over 3 orders of magnitude:
.
Download figure:
Standard imageHigh-resolution image4.3. The Ring
The ring flux densities are shown in Figure5, along with previous measurements of the ring from ALMA and the Australia Telescope Compact Array (ATCA): 1.4, 18, and 44 GHz (Zanardo et al.2013); 9 GHz (Ng et al.2013); and 94 GHz (Lakićević et al.2012). As the frequency increases toward the submillimeter and FIR, the contribution of the thermal ejecta emission to the total emission becomes more significant, so the ring emission at day 9280 was estimated from the total flux densities at 102, 213, 345, and 672 GHz by scaling and subtracting an ejecta model component in Fourier space based on the Band 9 ejecta flux density at each frequency (Zanardo et al.2014).
The Cycle 2 ring flux densities exhibit more scatter and are lower than the ATCA values by ∼30%. The integrated ring emission follows a nonthermal power-law profile of the formSν ∝να. The spectral indexα was previously estimated from ATCA radio (Zanardo et al.2014) and ALMA Cycle 0 data (Indebetouw et al.2014) to be −0.73 ± 0.02 at day 9280 (Zanardo et al.2014). An updated fit, using the ALMA Cycle 2 ring flux densities and the more recent 1–9 GHz measurements from Callingham et al. (2016) from day 9686, results in a slightly lowerα = −0.70 ± 0.06.
The radio ring emission has steadily increased owing to the synchrotron-producing electrons being accelerated by expanding shock waves (Zanardo et al.2010; Staveley-Smith et al.2014). Recently, Cendes et al. (2018) have fit the radio emission of the (2D) ring and (3D) torus emission models across many epochs as a power law of the formS(ν,t) =
, whereα is the spectral index of the emission across the spectrum,β is the power-law slope, andK is the offset constant for the given model. The ALMA ring flux densities are generally in excellent agreement (within 5%) with the Cendes et al. (2018) model prediction for day 10,402 whereK = 1.5 ± 0.1,α = 0.70, andβ = 0.59 ± 0.02. We find no evidence of a contribution from dust in the ring to the millimeter wavelength flux densities (see, e.g., Bouchet et al.2006; Matsuura et al.2019).
5. Analysis and Results of Molecular Lines
5.1. Analysis of Molecular Lines Using RADEX
In the previous section, we modeled the integrated SED of the dust emission under the assumption of uniform temperature and density within the ejecta. However, in Section3, we saw that the SiOJ = 5
4 and COJ = 2
1 images exhibit a hole (Figure3) where the dust emission peaks, and the SiO line ratio indicates lower SiOJ = 5
4 brightness with respect to SiOJ = 6
5 and SiOJ = 7
6 (seen in the SiO line ratio maps in Figure8). In order to understand the spatial distribution of the line ratios and intensities qualitatively, we use the non-LTE (local thermal equilibrium) line radiative transfer coderadex (van der Tak et al.2007).
Figure 8. Molecular line ratios of SiO. The line emission in each channel was converted to flux (W m−2) integrated over the 300 km s−1 bin before division. Channel velocity centers are LSRK; for reference, the SN 1987A systematic velocity is 287 km s−1. Left: SiOJ = 7
6/SiOJ = 5
4. Note that while the 300 km s−1 window covers the line centers, the SiOJ = 7
6 line has poor coverage there owing to that frequency falling at the edges of the two observed tunings. Right: SiOJ = 6
5/SiOJ = 5
4. Intensity maps were convolved to the SiOJ = 7
6 beam size before dividing. Cyan contours are 679 GHz dust 3σ and 5σ levels. The apparently high ratios at the edges of the map features are primarily due to reduced S/N in those outer regions.
Download figure:
Standard imageHigh-resolution image5.1.1. Description of the Procedure Adopted for SiO
This subsection describes the analysis of the SiO lines and intensities in detail. The analysis for CO was carried out in a similar manner; thus, it is described only briefly in that subsection.
radex calculates the molecular line intensities, using the escape probabilities from Osterbrock (1989), and the uniform sphere method for the gas distribution was chosen for this calculation. The code involves calculations of level populations, using the Einstein coefficients and collisional cross sections of molecular lines assembled by theLAMBDA database (Schöier et al.2005), and we use H2–SiO collisional cross sections based on the calculations by Dayou & Balança (2006). In the ISM, H2 is widely assumed to be the dominant collisional partner in molecular clouds; however, that is probably not the case for the ejecta of SNe. As a consequence of a series of nuclear burning processes, the progenitor star’s core will have built up layers of different newly synthesized elements, with hydrogen being depleted. In the two layers containing abundant Si, the major elements are O and S (e.g., Woosley1988), and the collisional partner of SiO is likely to be O2 and SiS. This would potentially change the collisional cross section by a factor of 1–10, depending on the transitions (Matsuura et al.2017).
One of theradex input parameters is the FWHM of the line width of the Gaussian (Δv). We adopted a Δv of 400 km s−1, whose integrated area over the Gaussian profile would be equivalent to that of a box-shaped 300 km s−1 line profile. If the line is optically thin, the assumed line profile is not a major issue. However, as we will see later in this analysis, the lines are mildly optically thick at the line center, but not at the side of the line profile, so we therefore make the assumption that the line profile only moderately affects the line ratios and the line intensities in this “mildly” optically thick regime.
The main parameters involved inradex calculations are the kinetic temperature (Tkin), the density of the collisional partner (ncoll), and the column density (NSiO). In the optically thin regime, as found in the calculated parameter range where solutions are found, the SiO line intensities are determined byNSiO together with the area filling factor and the expansion velocityvexp, while the SiO line ratios are determined byTkin andncoll, independent ofNSiO. The filling factor is defined as the area of the line-emitting fraction within the beam/pixel, following Goldsmith & Langer (1999). Previous analyses of SN 1987A with lower angular resolution suggested the range of 2.5%–45% (Kamenetzky et al.2013; Matsuura et al.2017); thus, we assume that the filling factor of 1%–50% is a reasonable range. In this analysis, we adopted a column density grid ofNSiO = 1012–1018 cm−2 in factor of 10 increments and searched for the predicted line intensities that can match the measured ones within the assumed range of the filling factor. The adopted ranges of the parameters arencoll = 103–1010 cm−3 andTkin = 10–200 K. Matsuura et al. (2017) suggested that the temperature range is below 190 K for SiO, with the CO kinetic temperature between 30 and 50 K, so we restricted the analysis to below 200 K. Although we include temperatures up to 200 K in theradex calculations, it is very unlikely that the majority of SiO gas has such a high temperature, and most likely, the overall SiO gas should have a temperature close to the CO temperature. We searched for matching SiO line ratios within these parameter ranges.
The ALMA data of the three SiO transitions have different beam sizes and orientations, so SiOJ = 5
4 and SiOJ = 6
5 were convolved to match the lowest spatial resolution—that of the SiOJ = 7
6 beam—with a uniform pixel width of 0015. The convolved and regridded line ratio maps, made on a channel-by-channel basis, are displayed in Figure8. Intensities were averaged over 5 × 5 pixels in order to increase the S/N in the SiOJ = 7
6 image. As the minor-axis FWHM of the beam is 017 for SiOJ = 7
6, there is a small loss of spatial information by this averaging.
One caveat: the continuum subtraction for SiOJ = 6
5 was performed on the final imaged data cube, whereas for SiOJ = 5
4 and SiOJ = 7
6 the continuum subtraction was done inuv space before imaging (imcontsub vs.uvcontsub inCasa). Continuum subtraction inuv space is generally considered preferable if the continuum dominates the line emission, and the difference could slightly affect the ratios, but at high S/N, as in the case of SiOJ = 6
5, the two methods should give similar results. As discussed in Section3, there was no appreciable difference between the two methods for the COJ = 6
5 data.
The uncertainties in the measured SiO lines and line intensities are dominated by calibration uncertainties, and we adopted 10% of the line intensities for this analysis. The intensity uncertainties that were measured as a fluctuation of the “blank” sky level are 3% for SiOJ = 6
5 and SiOJ = 5
4 and 5% for SiOJ = 7
6 after 5-pixel averaging, so that the uncertainties based on the observations are smaller than the systematic uncertainties from calibration errors. Because the dominant uncertainties are systematic, these propagate in an asymmetric way, i.e., if the actual line intensity of SiOJ = 6
5 were higher than the measured value, the line ratio of SiOJ = 6
5/SiOJ = 5
4 would increase, but the ratio of SiOJ = 7
6/SiOJ = 6
5 would decrease. The Cycle 2 flux densities being systematically lower than Cycle 0 as discussed in Section4.1 is not an issue for this analysis, as the important result here is that the relative change of temperature and density within the ejecta can explain the difference in molecular line ratios.
5.1.2. Results of SiO Analysis
We selected two representative spatial regions, shown in Figure9, in order to understand the change in physical parameters (Tkin,ncoll) within the ejecta, which in turn affect the SiO emission. One is close to the molecular hole seen in SiOJ = 5
4 and COJ = 2
1, which we call region A, and the other is representative of the neighboring “typical” SiO ejecta, named region B. The actual hole itself has almost negligible SiOJ = 5
4 line intensity, so we chose the nearest possible point for region A. All ratios in this analysis were determined from the 0 km s−1 LSRK channel, which has good S/N for all three SiO lines, and the hole is clearly identified.
Figure 9. Locations of regions A and B used for theradex analysis of the SiO transitions. The background image is the SiOJ = 5
4 line at the 0 km s−1 (LSRK) channel, which clearly shows the molecular hole at region A. Region B is a representative region of the general SiO-emitting ejecta.
Download figure:
Standard imageHigh-resolution imageIn region A, which has a hole in SiOJ = 5
4, we found that a column density ofNSiO = 1016 cm−2 is reasonable. ForNSiO = 1015 cm−2 or below, the predicted line intensities do not reach the measured levels, assuming a filling factor of 50%. Slightly larger values ofNSiO = 1017 cm−2 can match the measurements, but beyond 1018 cm−2 the predicted and measured line intensities do not match.
Figure10 demonstrates the plausible ranges for the kinetic temperature (Tkin) and collisional partner’s density (ncoll) at region A, forNSiO = 1016 cm−2. The collisional partner for theradex calculations is H2; see the last paragraph of this section for discussion of this point. In Figures10(a)–(c), the colored contours show the calculated ratios, the black lines show the measured SiO line ratios, and the 1σ background uncertainties and 1σ systematic uncertainties are indicated with dashed and dotted lines, respectively. Figure10(d) compares the possible ranges of the kinetic temperature (Tkin) and collision partner density (ncoll) from the three SiO line ratios. The uncertainty of the SiOJ = 7
6/SiOJ = 5
4 ratio has the tightest constraint, so only the uncertainty of this ratio is plotted in Figure10(d). Figure10(d) shows that the ratios of all three SiO lines can be matched with a very similar set ofTkin andncoll > 107 cm−3, as indicated by the solid lines of the three different ratios almost overlapping each other. There are multiple scenarios that match the measured SiO line ratios: the first possibility isTkin = 34 K with LTE conditions atncoll > 107 cm−3, the second possibility is a non-LTE condition with a range of 106 cm−3 < ncoll < 107 cm−3 and 50 < Tkin < 200 K, and finally the curves connecting these two conditions. The line center turns optically thick at the column density ofNSiO = 1016 cm−2, but the majority of the lines at off-center velocities are optically thin, so the overall analysis is not strongly affected by optical thickness at this column density. The filling factor at this column density is 9%–14%.
Figure 10. (a–c) Results ofradex modeling of the SiO line ratios at region A, representing the molecular “hole.” These calculations assume a column density ofNSiO = 1016 cm−2. (d) Combination of all three best-fit curves. For each line, the solid curve denotes the best-fit locus of the grid of kinetic temperature and H2 density values to the observed line ratios, while the dashed and dotted curves represent the 1σ background uncertainties and 1σ systematic uncertainties, respectively.
Download figure:
Standard imageHigh-resolution imageBy increasing the column density fromNSiO = 1016 cm−2 toNSiO = 1017 cm−2, the SiO ratios change much more gradually as a function of the kinetic temperature (Tkin) and collision partner density (ncoll) (Figure11), expanding the feasible parameter space. This is the effect of higher optical depth. The sets ofTkin andncoll from the three different calculated SiO ratios do not overlap in Figure11, but it is still possible to consider a wide range ofTkin atncoll > 107 cm−3, from 19 to 22 K within the 1σ uncertainty. In the non-LTE range, the requiredncoll > 107 cm−3 is more or less comparable to that forNSiO = 1016 cm−2. In order to reproduce the line intensities with this higher column density, the filling factor of the beam area is 0.8%–2%, much lower than for theNSiO = 1016 cm−2 case.
Figure 11. SiO line ratios that match the measured ALMA values at the SiO hole. Same as Figure10(d), but for a column densityNSiO = 1017 cm−2.
Download figure:
Standard imageHigh-resolution imageAt region B, the physical parameters are slightly different from those at the hole (region A). Figure12 shows the possible parameter space that fits the SiO ratios for region B; onlyNSiO = 1016 cm−2 gives a feasible range. The most plausible temperature isTkin = 18 K for LTE conditions (ncoll > 107 cm−3), and an alternative possibility is non-LTE with 3 × 105 cm−3 < ncoll < 107 cm−3 and 18 < Tkin < 200 K. In the optically thin regime, the filling factor and the column density are inversely correlated; thus, the accuracy of the filling factor is limited by our column density grid.
Figure 12. Same SiO calculations as in Figure10(d), but for region B with column densityNSiO = 1016 cm−2.
Download figure:
Standard imageHigh-resolution imageIn summary, the difference in the modeled line ratios and intensities near the SiOJ = 5
4 andJ = 6 
5 hole (region A) with respect to the “representative” ejecta region B can be explained in the following three ways. The first possibility is that the hole region has a higher temperature (35 K) than the surrounding locations (19–22 K), with both having LTE conditions—i.e., high density of the collisional partner. The second possibility also requires LTE conditions, but instead of high temperature, the hole region has a higher column density in a much smaller area, represented by a small beam filling factor. The third possibility is that the entire area is non-LTE, i.e., a lower density of the collisional partner, but with the hole having a factor of a few to a few tens higher density of the collisional partner than the surrounding region. These three explanations are not exclusive to each other; a mixture of these three scenarios is possible.
The uncertainties involved with this analysis arise from uncertainties in the collisional cross section. The collisional partner is most likely not H2, but rather other molecules such as O2 or SiS, according to chemical models (Sarangi & Cherchneff2013). Therefore, instead of higher H2 density at the hole, the collisional partner may be different in the hole and in region B—for example, region B may be dominated by collisions with O2, whereas in the hole the dominant partner could be SiS. However, as will be discussed in Section6.3.1, hydrodynamic simulations predict that such spatial distributions for the Si and O atoms are unlikely; therefore, our conclusion that there is a higher temperature, column density, or density at region A is still valid even with this uncertainty.
5.1.3. CO Analysis and Results
Although the COJ = 6
5 line does not have enough S/N for a quantitative analysis on a pixel-by-pixel basis, we can sum the spectra in independent (single beamwidth) regions to aid our analysis. The top panel of Figure13 shows the location of 20 regions selected across the CO-emitting ejecta, overlaid on top of the COJ = 2
1 emission. Nine regions, each the size of the COJ = 6
5 beam, are highlighted as areas of interest, potentially probing different conditions (numbered 1–9 in rows from left to right). The middle panel of Figure13 compares the summed spectra of COJ = 6
5 and COJ = 2
1 (with the latter convolved to the COJ = 6
5 beam before integration) for these nine regions showing their location with respect to the COJ = 2
1 emission. The spectra are in units of mJy, having been spatially integrated over each region.
Figure 13. Top: regions chosen for comparing the COJ = 2
1 and COJ = 6
5 fluxes and line profiles; each region spans one COJ = 6
5 beam. Regions of particular interest are labeled 1–9. The shaded grayscale contours illustrate the COJ = 2
1 emission (bottom left panel of Figure3). Middle: qualitative comparison of the stacked COJ = 2
1 (blue) and COJ = 6
5 (red) spectra in each region, with a representative scale shown in the lower left corner. Bottom: zoomed-in view of the spectra for regions 1–9 shown in the middle panel. The vertical dashed line is the systemic velocity of SN 1987A, 287 km s−1 (LSRK).
Download figure:
Standard imageHigh-resolution imageThe bottom panel of Figure13 provides a zoomed-in view of these spectra for a more detailed comparison. Across the CO ejecta, we see that the majority of the line profiles in the COJ = 2
1 and COJ = 6
5 transitions are similar to each other in the scaled spectra (see, e.g., regions 6, 7, and 9 in Figure13). However, in regions 1 and 8, the COJ = 6
5 profile is suppressed at negative velocities with respect to the COJ = 2
1 line. The COJ = 6
5 emission is also suppressed with respect to COJ = 2
1 in region 3, though this is across the entire velocity profile. At the location of the CO molecular hole (regions 4 and 5), we see strong COJ = 6
5 emission at velocities of −1000 to +1000 km s−1 with respect to COJ = 2
1 (bottom panel of Figure13); this continues to neighboring region 2. We note that regions A and B from the SiO analysis fall within the CO map regions 4 and 8, respectively, but they are not interchangeable. They were selected independently and are centered at different locations and serve different purposes (region B is representative of the general SiO-emitting properties across the molecular ejecta based on the SiO line ratios, whereas Figure13 demonstrates that region 8 has very different CO gas properties from most of the other regions).
Using this information, we carry out an analysis withradex similar to that for SiO. Figure14 shows the resultingradex calculations of the COJ = 6
5/2 
1 ratios. The three black and white lines show the curves for flux ratio values of 38, 20, and 3, respectively, corresponding to the higher, intermediate, and lower ends of the line ratios observed across the 20 regions. We note that the units of the flux densities used to derive the line ratios in this (CO) section are in Jy and are therefore different from the W m−2 used in the SiO ratio calculations. This is because the CO spectra are compared in velocity space, whereas we use integrated line intensities for the SiO analysis.25
Figure 14. Ratios of COJ = 6
5 and COJ = 2
1 flux densities in Jy units derived usingradex for the optically thin case. The black and two white lines show ratio values of 38, 20, and 3, respectively, indicating high, intermediate, and low CO 6–5/2–1 end ratios observed in the 20 regions in Figure13.
Download figure:
Standard imageHigh-resolution imageAs demonstrated originally in Figure 3 in Kamenetzky et al. (2013), the COJ = 6
5 line is sensitive to temperature change. We propose here that the COJ = 6
5 suppression with respect to COJ = 2
1 indicates that the gas is at a lower temperature in regions 1 and 8 (where we also see a peak in the dust emission; Section3, Figure4) compared to the surrounding regions. In these regions, the blue wing of the COJ = 6
5 emission is lower; thus, if the dust and COJ = 6
5 are spatially coincident in regions 1 and 8, this could imply that the CO and dust originate from a discrete region on the near side of the ejecta, though this is speculative, as we have no velocity information on the dust. Due to the low S/N of the COJ = 6
5 line, we cannot specify the exact excitation temperatures in these regions. However, the CO excitation temperature is higher near the COJ = 2
1 hole (regions 4 and 5) at velocities of −1000 to 1000 km s−1.
6. Discussion
6.1. Dust Properties from the Integrated SED
Any dust model that produces a mass higher than the total abundance of metals formed in the SN ejecta is clearly unphysical. Here we discuss whether the observed, or predicted, metal yields formed in the ejecta of SN 1987A can be used to rule out some of the dust varieties and compositions that produce good fits to the SED (Table4). Inferred dust masses of several tens of solar masses, as in the case of TiO2 (Posch et al.2003), which requires 82
of dust, or MgAl2O4 (Fabian et al.2001), which requires 122
, for example, are clearly untenable, as they are larger than the progenitor star mass (18–20
; Woosley1988). This rules out a further three dust varieties in Table4 (all of which satisfy our
criteria for a good fit).
Simple upper limits to the dust masses of various dust species can be estimated by calculating the resulting mass for 100% of the elements in the nucleosynthesis models being locked into dust, ignoring all chemistry, mixing, and other physical limitations. Such a highly unrealistic scenario is useful for winnowing out dust models that yield dust masses that are too large. Considering the 18
Z = 0.008 Z⊙ progenitor model from Nomoto et al. (2013), the total mass of the key limiting metals C, N, Ne, Mg, Si, S, and Fe is 0.77
, with 1.21
of oxygen. Focusing on individual dust varieties, the limits can be further differentiated. The Nomoto et al. (2013) 18
model predicts 0.149
of carbon, putting a limit on the mass of graphite and amorphous carbon grains, as well as a rough limit for PAH varieties since their masses are dominated by carbon. The carbonaceous grain model that both gives a good fit to the SED and produces the nearest fitted mass to the predicted carbon yield is the amC (AC1 sample) from Rouleau & Martin (1991), with a dust mass of 0.43
. For silicate dust, good fits to the SED produce masses of 0.3–0.7
, though the yields from Nomoto et al. (2013) would result in a maximum combined silicate metal mass (and therefore dust mass, limited by the available Mg) of <0.4
.
The iron yield from the Nomoto et al. (2013) 18
nucleosynthesis model is 0.079
(56Fe only) or 0.085
(including all isotopes); this is orders of magnitude less than the inferred pure iron dust model (Henning & Stognienko1996), which requires 3.97
of dust to fit the SED. The Woosley & Weaver (1995) 15 and 18
Z = 0.1Z⊙ progenitor models predict iron masses ranging from 0.14 to 0.20
, where roughly half of the iron originates from56Ni decaying to Fe. This is an order of magnitude less mass than the fit requires. We can therefore rule out a scenario where iron-rich grains in SN 1987A are producing the bulk of the thermal emission, due to either not fitting the SED (in the case of FeO and FeS, Henning et al.1995; Henning & Stognienko1996) or resulting in unrealistically high dust masses for the pure iron model.
Limits from the Nomoto et al. (2013) 18
explosive synthesis model for some other common dust varieties include 0.337
for forsterite (Mg2SiO4), 0.481
for enstatite (MgSiO3), 0.014
for alumina (Al2O3), and 0.373
for silica (SiO2)—assuming that 100% of all isotopes are locked in dust grains. These are all notably lower than the dust masses resulting from the modBB fits—e.g., 4.0
for forsterite, 4.1
for enstatite, and 0.9
for alumina (Table4). Given the discussion above, one can place a rough upper limit on the total mass of dust that could form in the ejecta of SN 1987A of <1.5
given the available metal budget. From this it is possible to immediately rule out a further seven dust varieties listed in Table4 as producing unphysical dust masses.
Many of the fits to the SED of SN 1987A require more mass in dust grains than the mass of available metals for the corresponding species. This can be explained if the SED is made up of a mixture of several species, each contributing to the overall dust budget, as originally proposed in Matsuura et al. (2015). For example, locking all available mass into a combination of C+MgSiO3+FeS would give a total metal mass, and an upper limit on the total dust mass, of 0.72
. However, taking the yields of the Woosley (1988) 18
0.1Z⊙ progenitor model, for example, would give 0.55
of total metals available for dust formation for this combination of species, or ∼30% less mass than the simple sum of the total species indicates.
Using the predicted metal yields as an upper limit to rule out dust varieties and determine the mass of the ejecta dust also has its own challenges in that model abundances for core-collapse SNe vary with different models, uncertainties in the nuclear process assumed, rotation, and the implementation of the artificially induced shock explosion model. We therefore caution that the model abundance yields can only provide loose upper limits. Considering the various limitations and caveats for the different dust models considered in this study, a likely overall dust composition, based on the measured SED and nucleosynthesis limits, is a combination of amorphous silicates (especially those of reasonably high emissivity, such as the Demyk et al.2017 model) and amorphous carbons as also concluded by Matsuura et al. (2015), limiting the total dust mass to potentially <0.7
.
6.2. The Spatial Distributions of Dust, CO, and SiO, and Chemistry Leading to Dust Formation
Our spatially resolved images (Figure3) show that the dust distribution is clumpy and asymmetric. Comparing the spatial distribution of the dust with the COJ = 6
5 reveals a weak anticorrelation with the integrated COJ = 6
5 and dust distributions, while there is little spatial correlation between the dust and SiO images. We suggested earlier that this anticorrelation occurs because COJ = 6
5 is suppressed compared with COJ = 2
1 in the dust-bright regions, indicating that the excitation temperature of CO is lower than in other neighboring regions. This provides new information about the chemistry and physics involved in the formation of dust.
Note that we see an exception to this in one region (region 4 in our CO analysis and roughly corresponding to region A in our SiOradex analysis). Here, the hole in SiO and COJ = 2
1 coincides with the dust peak, with relatively strong COJ = 6
5 emission observed at −1000 to 1000 km s−1 with respect to COJ = 2
1, and this strong COJ = 6
5 continuing to its neighboring region (region 5 in Figure13). We will discuss this region separately in Section6.3.1.
The SN chemistry after the explosion inherits a series of nuclear synthesis processes at the stellar core prior to and during the SN explosion (e.g., Sarangi & Cherchneff2013). The outermost region is the H envelope, followed by the He shell, which can contain more carbon atoms than oxygen atoms (Woosley & Weaver1995; Rauscher et al.2002). The He shell can also form CO, as it enriches with C and O (e.g., Sarangi & Cherchneff2013). The inner region is roughly represented by an O+Ne zone, which can also contain C, followed by an O+Mg+S+Si zone, and finally with a56Ni core, which also contains Si, but very low C or O. Here we interpret an apparent anticorrelation between the dust and COJ = 6
5 spatial distributions as the result of both CO and dust components having originated from the He envelope or O+Ne nuclear burning zone containing both C and O prior to the explosions. Thus, we propose that the dust grains formed in this region could be carbonaceous.
Our suggestion of carbonaceous dust contradicts predictions of SN dust formed in chemical models. Sarangi & Cherchneff (2015) predict that SiO molecules formed in the O+Mg+S+Si zone and SiO molecules directly condense into silicate dust. The majority of nuclear zones have more O atoms than C atoms, and the formation of CO blocks the formation of graphite or amorphous carbon. Deneault et al. (2006) and Clayton (2011) previously discussed the formation of carbonaceous grains by dissociating CO via highly energetic electrons, making unbound carbon available for carbonaceous dust formation. However, the Clayton (2011) calculations involved only a few reaction rates involving C and CO. In contrast, Sarangi & Cherchneff (2013,2015) included more extensive chemical networks, as well as dissociation of molecules by energetic electrons, and found that this resulted in very few carbonaceous dust grains (6 × 10−3M⊙ of carbonaceous dust out of a total of 0.04M⊙ dust mass for a 19M⊙ star). One simple explanation is that Sarangi & Cherchneff (2013,2015) modeled the SN chemistry only up to day 1500, and the dust composition might have changed since then. However, it still requires CO to be dissociated in order to remove the blockage of carbonaceous dust formation, suggesting that any dissociation process must be more efficient on longer timescales. Alternatively, the chemical reaction rate used in the models might have large uncertainties, particularly involving the highly energetic electrons. Together with the recent detection of HCO+ (Matsuura et al.2017), which was largely underpredicted in the chemical model of Sarangi & Cherchneff (2013), this suggests that some tensions exist in molecular chemistry models of SNe. However, this tension between our proposed C-rich grain formation and the lack of C-rich grains “grown” in the SN dust models could be alleviated if macroscopic mixing is efficient enough to allow Si and C dust to form in the same regions. Indeed, some 3D hydrodynamical models (e.g., B15, N20; see Utrobin et al.2019) suggest that ∼30%–70% of the Si can be mixed out of the central regions into the C shell (A. Wongwathanarat 2019, private communication).
An alternative explanation for the anticorrelation between dust and CO suggested in this work is that carbonaceous dust could originate from the He layer while the CO gas is restricted to the C + O core, possibly resulting in a projected anticorrelation between dust and CO. However, this does not explain why the gas temperature is lower at the dust-emitting region (though the presence of dust can lower the gas temperature; see below). In this scenario, the CO might have originated from two different layers of nuclear burning zones.
Our ALMA observations have also shown that regions of bright dust emission have lower CO excitation temperature than other regions. There are two possibilities to explain this. First, more dust may have led to cooler gas temperatures owing to radiative cooling via dust emission, i.e., the cooler gas temperature is the consequence of dust formation. Second, the temperature of these regions may have already been lower in the early days when the dust grains were formed, and the temperature reached the dust sublimation temperature while the gas density was reasonably higher, driving more efficient dust condensation. In this case, dust formation is the consequence of cooler gas temperature. Currently, the data do not provide a way to distinguish between these two cases.
We note that we see faint diffuse dust emission that might be more extended beyond the Band 9 dust structure in Figure3, for example, in comparison with the 315 GHz image, but the surface brightness is lower in those locations. The CO–dust anticorrelation is not as obvious in these fainter, extended regions; therefore, we do not claim that the entire dust content of the SN 1987A ejecta is carbonaceous in nature; some of the dust components within the system could be associated with silicates as well.
The high-resolution ALMA observations in this work show that the dust distribution in the ejecta of SN 1987A is clumpy even at small scales (as also seen in the CO and SiO ejecta shown originally by Abellán et al.2017). Aside from SN 1987A, we know from the knots and filaments observed across the Galactic SNR Cassiopeia A (e.g., Milisavljevic & Fesen2015) that the gas ejected in SN explosions can be clumpy. Interestingly, Sarangi & Cherchneff (2015) included clumpiness in their chemical models and showed that clumpy gas, compared with smoothly distributed gas, provides density enhancements in the ejecta, resulting in larger grain sizes for SN-formed dust.
6.3. Interpreting the Bright Point Source (the Blob)
6.3.1. The Blob: Predictions from Hydrodynamical Models
The bright dust peak (the blob) is observed at the location of the hole in the COJ = 2
1 and SiOJ = 5
4 line-emitting ejecta, and we also see emission located in the hole in the SiOJ = 7
6 transition. To explain this additional SiO emission, our analysis of the SiO line ratios gives three possibilities: a higher SiO temperature, a high column density with high area filling factor, and high density of the collisional partner. The higher intensity of the emission of the dust blob compared to its surroundings could be explained by higher dust (“column”) density or by higher dust temperatures inside the blob. For both dust and SiO molecules, the two key physical parameters are temperature and density, where the latter is also associated with the column density.
We first discuss the possibility of enhanced density (of dust, as well as SiO and CO) in the blob based on model predictions. To form the SiO and dust, a significant amount of Si must be present in the blob region. The dust could also be carbon based, such that instead of Si, C could also be enhanced. However, recent hydrodynamical simulations (M. Gabler et al. 2019, in preparation) compared with ALMA observations (Figure 4 in Abellán et al.2017) show that in the explosions of three out of the four simulated progenitor models, the final density distributions of Si and C (each multiplied by the oxygen density) have a void in the center. This void is caused by a strong reverse shock from the He/H interface within the ejecta immediately after the explosion. This shock first slows down the expansion speed of the inner ejecta (also containing Si and C) when passing them, but then it compresses the material such that the entropy increases significantly and an outward-moving shock forms. This feature is termed the “self-reflected” shock, first discussed in Ertl et al. (2016) (with forthcoming details in M. Gabler et al. 2019, in preparation). The self-reflected shock accelerates the innermost ejecta compared to homologous expansion and leaves a region with lower density in the center. When passing more and more of the ejecta, the shock loses energy. It finally dissipates and cannot accelerate the outermost ejecta. This acceleration of the low-velocity ejecta leads to the formation of a higher-density shell-like configuration as observed in the emission of the transition lines SiOJ = 5
4 or COJ = 2
1 (Figure 4 in Abellán et al.2017). However, in the hydrodynamic simulations, one model based on the progenitor model from Shigeyama & Nomoto (1990) (model N20 in Wongwathanarat et al.2015; Abellán et al.2017) leads to a weaker reverse shock and, hence, a weaker self-reflected shock. The latter then is not able to significantly accelerate the central ejecta. Therefore, the central region of this explosion model still has similar densities of Si and C compared to that of the fastest-moving Si- and C-rich ejecta. While the N20 model may offer the possibility of a high central density, the light curve (Utrobin et al.2015,2019) disfavors that model.
6.3.2. The Blob: Gas and Dust Heated by the Compact Source
To explain the higher SiO and CO gas temperature and increased brightness for the observed properties of the blob, we suggest two possibilities: (i) gas heated by a compact object or (ii) a clump heated by44Ti decay. Here we argue that the most probable explanation for the detected dust blob is that the innermost part of dust and gas is heated by radiation from the compact object, with an early development of a pulsar wind nebula.
44Ti was synthesized at the time of the SN explosion, and its decay energy is the main source of the heating of the inner ejecta (Jerkstrand et al.2011; Larsson et al.2016). SN explosion models show that44Ti is located more or less at a similar radial extent to56Ni and28Si (Woosley & Weaver1995; Wongwathanarat et al.2015,2017) with almost identical bulk velocities, though the modeled distribution of44Ti could be more uncertain than the other elements owing to its greater sensitivity to the explosion physics (A. Jerkstrand et al. 2019, in preparation). Observations show that the44Ti ejecta is redshifted, suggesting that the bulk is moving away from the observer (Boggs et al.2015). Since the predicted distribution of44Ti shows qualitatively similar properties to those of Figure 4 in Abellán et al. (2017), because it is subject to effectively the same hydrodynamical history (A. Jerkstrand et al. 2019, in preparation), it is unlikely that gas heated by44Ti decay would be more centrally distributed and also be colocated with a small blob of gas and dust at one small region in the very inner ejecta.
Therefore, we suggest that the central blob is due to warmer ejecta, and the most likely source of the heating energy is from the compact object (i.e., possibility (i)). The dust and gas in the blob could be directly heated by X-rays from the surface of the compact object (Alp et al.2018a), or the dust could be heated by synchrotron radiation generated by the compact object.
In the latter case, previously it was shown that synchrotron radiation from neutron stars can heat up the dust grains in pulsar wind nebulae, as seen in Galactic SNRs, including the Crab Nebula (Temim et al.2006; Gomez et al.2012b; Owen & Barlow2015), G54.1+0.3 (Temim et al.2017; Rho et al.2018), G11.2−0.3, G21.5−0.9, and G29.7−0.3 (Chawner et al.2019). In the Crab Nebula, the radiation from the pulsar wind can partially dissociate gas, contributing to the formation of the molecules OH+ and36ArH+ (Barlow et al.2013). SN 1987A may currently be undergoing this phase, i.e., just beginning to develop a pulsar wind nebula in the innermost region. If the compact source is a black hole (Brown et al.1992; Blum & Kushnir2016), instead of a neutron star, jets from the black hole can also heat up gas and dust locally (e.g., Russell et al.2006).
Instead of synchrotron, direct thermal radiation from the neutron star can also heat the gas and dust locally (Alp et al.2018a). The SN 1987A compact object is still surrounded by a dense metal-rich gas. Metals absorb and scatter (soft) X-ray, UV, and optical light efficiently, so that it would be challenging to detect light from the compact source directly (McCray & Fransson2016; Alp et al.2018b). The absorbed energy is reprocessed into longer wavelengths and eventually ends up heating dust and molecules. However, this picture does not consider any radio emission produced directly by the compact object.
From the neutron star kick inferred by the distribution of intermediate-mass elements (e.g., Katsuda et al.2018) and the redshifted44Ti spectrum (Boggs et al.2015), the direction of the compact object is predicted to be moving toward us and extending along the northeast direction in the sky projection (e.g., Janka et al.2017). The location of this blob is consistent with that prediction. There is an offset between the dust blob and the estimated location of the progenitor star as derived by Alp et al. (2018a):α = 05h35m27
9875,
(ICRS). The brightest pixel in the 679 GHz emission is offset 72 mas to the east and 44 mas to the north of the position of the progenitor star from Alp et al. (2018a). This offset is ∼3–5 times the total alignment uncertainty; however, it is still within the range of the neutron star kick. Zanardo et al. (2014) proposed that the compact object may have traveled 20–80 mas from the site of the SN, toward the west, in comparison to the circular radius of 100 mas used in Alp et al. (2018a). The peak of the dust blob from our work falls within this range (though it is at the upper end, at approximately 700 km s−1 in our data, assuming the SN 1987A position of Alp et al.2018a).
Moreover, the location of the compact object and the brightest part of the pulsar wind in the Crab Nebula are also not coincident (Weisskopf et al.2000; Gomez et al.2012b). Although the Crab’s pulsar is powerful enough to affect its environment dynamically, whereas the compact object in SN 1987A would be at an earlier evolutionary stage, nevertheless, a misalignment between the estimated location of the progenitor star and the location of the dust blob is not unprecedented.
The approximate temperature of the dust of the blob can be estimated from the 679 GHz flux density. The peak is a factor of ∼2 brighter than the surrounding ejecta continuum. A factor of two increase in flux density at the dust peak can be explained if the dust is at a higher temperature compared to the global dust ejecta (an increase to 33 K from 22 K for Zubko et al.1996 or to 26 K from 18 K for Jäger et al.2003 dust grains would be required to explain the peak), though this temperature is much lower than predicted for dust heated by a pulsar in SN 1987A (Omand et al.2019).
We estimate the dust peak flux density asS679 GHz = 3–5 mJy, depending on how a 2D Gaussian is placed on the 679 GHz map; the peak pixel flux densities are at S/N ∼ 7 (above the rms level), and the uncertainty is caused by the source being blended with other nearby features. This flux density contains contamination from the underlying continuum emission, which can contribute about 2 mJy; thus, the estimated flux density of the compact source is of order 1–2 mJy. Figure15 shows our estimated 679 GHz flux density range of 1–2 mJy (orange bar), along with spectra for the Crab Nebula and its central pulsar (Bühler & Blandford2014), scaled to the distance of SN 1987A, for comparison. The estimated blob flux density falls between the total spectrum of the pulsar wind nebula and the sole pulsar spectrum, using the Crab as a template. Using our 1–2 mJy range and assuming a dust model, one can estimate the millimeter spectrum and the luminosity of the compact source. For the Zubko et al. (1996) ACAR model and a temperature of 33 K, we obtain the light-orange shaded region in Figure15. Alp et al. (2018a) derived upper limits to the flux densities of the compact object, shown here as blue triangles. At 213 GHz, our ACAR modBB curve givesSν of order 0.1 mJy, which is consistent with their flux density limits around that frequency. An alternative limit on the compact object emission can be estimated from the bolometric luminosities of different components of the ejecta. Integrating the blob modBB curve assuming ACAR dust results in a localized bolometric luminosity
= 40–90 L⊙ for 679 GHz flux densities of 1–2 mJy. This is an order-of-magnitude estimate, incorporating uncertainties in the flux density measurements and temperature estimate. This luminosity range is an upper limit assuming that the compact source heats dust from 0 to 33 K. However, the compact source is most likely additional heating, on top of44Ti-heated 22 K dust, which might be partly subtracted as the underlying 2 mJy continuum, but some of this contribution might not be subtracted, so that the power coming from the compact source could be lower than 40–90 L⊙.
Figure 15. Limits on the luminosity of the compact object in SN 1987A. The flux densities for the entire ejecta are shown by the yellow curve. The estimated 679 GHz compact object flux density range (1–2 mJy) is denoted by the vertical orange bar. The light-orange shaded region shows the Zubko et al. (1996) ACAR model for this flux density range and assumes that it corresponds to a temperature of 33 K. The flux densities assuming this ACAR model and 1–2 mJy at 679 GHz are consistent with the limits from Alp et al. (2018a) (blue-gray triangles). For comparison, the spectra of the Crab Nebula and the Crab pulsar (purple and red curves, respectively; Bühler & Blandford2014, and references therein) are shown scaled to the distance of the LMC.
Download figure:
Standard imageHigh-resolution imageIn summary, we suggest that the dust blob seen in the ALMA Cycle 2 Band 9 images could be due to dust heated by the compact object and potentially an emergence of a pulsar wind nebula based on the following arguments: (i) we expect a compact source to be present; (ii) we see one and only one blob, which is difficult to reconcile with the expected geometry of44Ti; (iii) this scenario would produce a temperature increase at the location of the blob as proposed in this work; and (iv) the position of the dust blob is within the predicted SN kick, though toward the high end. However, we caution that we only have one frequency band, and as such, the nature of the dust peak is not clear: this argument is only valid if the blob is thermally or nonthermally heated emission from dust. Alternative possibilities, such as the direct detection of the compact object spectrum, cannot be ruled out in this work.
6.4. Dust as a Source of Extinction
The precise origin of the Hα hole, whether from a physical lack of material (due to the reverse shock or a simple void due to the expansion of ejecta), or from illumination of the ring X-rays brightening the outer rim of the ejecta, or from dust extinction, has been discussed for many years (e.g., McCray2003; Larsson et al.2011,2016,2019; Fransson et al.2013). A simple estimate of the extinction can be made by assuming that the dust fills a spherical shell with uniform density and calculating the optical depth using the mass extinction coefficientκext (related to, but different from, theκabs used in the modBB fits) according to
. For shell edges corresponding to the features in the images (005–02, or ∼(0.4–1.5) × 1015 m), the Zubko et al. (1996) ACAR model and its corresponding dust mass fit give an opacity at 6563 Å ofτHα ∼ 560. Silicate dust varieties tend to have lower
at a given wavelength—the amorphous forsterite model of Jäger et al. (2003), for example, givesτHα ∼ 400.
at optical wavelengths can depend strongly on the assumed grain size for some models, typically decreasing as grain size increases. Yet even for very large grains ofa = 5μm, the optical depth of ACAR dust at Hα is 10. Even this level is optically thick, meaning that dust extinction can play a nonnegligible role in the observed Hα distribution.
7. Summary
We have observed SN 1987A with ALMA in Cycle 2, 10,352–10,441 days after the explosion, in Bands 7 and 9. This follows on from the first ALMA results from Cycle 0 in Kamenetzky et al. (2013), Indebetouw et al. (2014), Zanardo et al. (2014), and Matsuura et al. (2015) and complements the molecular line studies of Matsuura et al. (2017) and Abellán et al. (2017) from Cycles 2 and 3. In this paper we describe the observations, data reduction, calibration, and photometry in the ALMA bands at the highest angular resolution to date for the continuum of SN 1987A.
- 1.We find that the dust emission in the ejecta is clumpy and asymmetric, fitting within the Hα “keyhole,” with a peak of emission that we name “the blob.” The dust ejecta region is smaller in scale than the cool, clumpy CO and SiO ejecta regions. Dust, in the amounts we fit here in a simple uniform spherical shell geometry, is optically thick at optical wavelengths, withτ ∼ 500. Dust extinction could be a factor in the appearance of the Hα hole.
- 2.We see an anticorrelation between the COJ = 6
5 emission and dust, and this anticorrelation is not seen when comparing dust and COJ = 2 1. Ourradex analysis suggests that this is the result of a lower CO gas temperature where the dust emission is stronger, compared to the surrounding ejecta, hinting that the dust may be C-rich and may have formed as a result of dissociation of CO, contrary to chemical predictions. - 3.We observe a dust peak (the blob) at the location of the molecular hole detected at lower CO and SiO transitions, and the higher CO and SiO line transitions are stronger in that location than the emission in the lower transitions. We suggest that this is the result of warm gas and dust at the location of the blob, and we discuss the possibility that this could be due to slow-moving reverse shock material originating from the explosion (the self-reflected shock), heating from a high concentration of radioactive decay, or the compact source. The most likely scenario is an indirect detection of the compact source.
- 4.We fit the SED of the dust emission from the ejecta with modified blackbody profiles. The derived dust masses and temperatures depend on the submillimeter emissivity of the dust, which is not very well determined. Temperatures from the fits are generally between 18 and 23 K. Amorphous carbon and graphite models have dust with high emissivity and so yield the lowest dust masses, around 0.4–1.6 M⊙. Silicates return higher dust mass estimates in the range of 0.6–4 M⊙. “Typical ISM” dust varieties from the Milky Way and nearby galaxies give masses between 1 and 2 M⊙. Taking the total mass of available metals excluding oxygen predicted to be ejected by an SN with progenitor mass appropriate for SN 1987A, we rule out several grain models and compositions for the ejecta dust. We revise the possible range of dust masses in the ejecta to 0.2–0.4 for carbon or Si grains, or a total of <0.7 for a mixture of grain species. A mixture of dust species, including silicates and carbonaceous grains, seems necessary to reconcile the continuum SED, nucleosynthesis model yields, and the molecular line analysis.
We acknowledge Arka Sarangi, Loretta Dunne, Steve Maddox, and the anonymous reviewer for interesting and useful discussions. P.J.C. and H.L.G. acknowledge support from the European Research Council (ERC) in the form of Consolidator GrantCosmicDust (ERC-2014-CoG-647939, PI H L Gomez). M.M. acknowledges support from the UK STFC Ernest Rutherford fellowship (ST/L003597/1). Research at Garching was supported by the ERC through grant ERC-AdG No. 341157-COCO2CASA, and by the Deutsche Forschungsgemeinschaft through Sonderforschungsbereich SFB 1258 “Neutrinos and Dark Matter in Astro- and Particle Physics” (NDM) and the Excellence Cluster “ORIGINS: From the Origin of the Universe to the First Building Blocks of Life” (EXC 2094;http://www.origins-cluster.de). M.J.B. acknowledges support from ERC grant SNDUST (ERC-2015-AdG-694520). J.C.W. is supported in part by NSF AST-1813825. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00063.S, ADS/JAO.ALMA#2013.1.00280.S, ADS/JAO.ALMA#2015.1.00631.S, ADS/JAO.ALMA#2016.1.00077.S, and ADS/JAO.ALMA#2017.1.00221.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.
Facility: Atacama Large Millimeter Array. -
Software:astropy (Astropy Collaboration et al.2018),Casa (McMullin et al.2007),emcee (Foreman-Mackey et al.2013),ipython (Pérez & Granger2007),lmfit (Newville et al.2016),matplotlib (Hunter2007),multicolorfits (Cigan2019),numpy (Van der Walt et al.2011),radex (van der Tak et al.2007),SAOImage DS9 (Joye & Mandel2003),scipy (Jones et al.2001),WWCC (Arshakian & Ossenkopf2016).
Appendix A: Molecular Line Channel Maps
The channel maps for molecular transitions discussed in this work are given here. The CO maps are shown in Figure16, and the SiO maps are shown in Figure17.
Figure 16. CO channel maps in 300 km s−1 steps, LSRK; for reference, the SN 1987A systematic velocity is 287 km s−1. Left: COJ = 2
1 as presented by Abellán et al. (2017), with the notable central hole in the molecular emission that persists all the way through the line of sight. Right: COJ = 6
5 from the present work. Even with the poorer S/N due to atmospheric transmission in Band 9, the emission profile is noticeably different from that of COJ = 2
1. Cyan contours are Band 9 dust at 3σ and 5σ levels, and the beam size is given by the green ellipse in the lower right panel.
Download figure:
Standard imageHigh-resolution image
Figure 17. SiO channel maps in 300 km s−1 steps, LSRK; for reference, the SN 1987A systematic velocity is 287 km s−1. Top left: continuum + SiOJ = 5
4 as presented by Abellán et al. (2017), also exhibiting the central molecular hole seen in COJ = 2
1. Top right: continuum + SiOJ = 6
5 from Abellán et al. (2017). Bottom: SiOJ = 7
6 from the present work. The emission profiles of the three lines tend to be similar for a given channel, with the conspicuous exception of an excess in SiOJ = 7
6 at the spatial position of the molecular hole but slightly in front (at 0 km s−1). Cyan contours are Band 9 dust at 3σ and 5σ levels, and the beam size is given by the green ellipse in the lower right panel.
Download figure:
Standard imageHigh-resolution imageAppendix B: Defining the Center for Photometric Measurements
The inferred center of the SN 1987A system appears to vary slightly in different parts of the spectrum, and it is important to determine it carefully for analysis of the expanding ejecta material. Alp et al. (2018a) used the hot spots in the ring to determine the center of SN 1987A. This was done by fitting 2D Gaussians to the hot spots in theHST/ACSR-band image from 2006 and then fitting 1D Gaussians along the same directions in all the other 32R andB images from 2003 to 2016 (Table 8 of Alp et al. (2018a)). They derive a ring center position ofα = 5h35m279875,δ = –69°16′11107 (ICRF J2015.0), with uncertainty from bootstrapping the hot spot locations of (11, 4) mas. Here we explore similar methods where we fit the ring emission to determine the central position of the SN 1987A system, but we use the Cycle 2 ALMA data at 315 GHz because it has the highest spatial resolution and S/N ring emission among the new images presented in this work. First, the ejecta and all emission clearly exterior to the ring were masked. Then, a ridge of emission peaks around the ring was determined. This was carried out by starting from the rough center of the map, and at 50 different angles the pixel with the brightest flux along each ray was located. Finally, an ellipse was fit to the ridge using three methods (Figure18): (i) fitting using a quadratic curve method (red), (ii) standard least-squares minimization of the parametric ellipse equation (blue), and (iii) weighted least squares, using the inverse of the squared pixel intensities as the weights (green) so as to favor bright regions over tenuous emission. The last method weights the pixels higher on lines of sight with brighter values than those lines of sight with very faint emission. All three methods, though the resulting ellipses have slightly different shapes, returned center positions within about 10 mas of each other—less than 1 pixel (Figure18). The final weighted fit gives a ring emission center ofα = 5h35m27998,δ = –69°16′11107 (ICRS). We estimate an uncertainty in these fit coordinates of 5.9 mas by varying the number of search angles and by averaging the results from the three fit methods. Combined with the 12 mas astrometric uncertainty for the 315 GHz image (Section2.1), the total uncertainty on the inferred center position is 18 mas, slightly larger than 1 pixel width. The main (systematic) uncertainty in the position comes from assuming that the progenitor/SN should coincide with the center of the ring: the optical ring probes high-density gas, whereas ALMA probes lower-density gas at higher latitudes. The revised center from the ALMA data used here falls on the southern edge of the central hole of the Hα emission.
Figure 18. Illustrating the three methods used for determining the center of the SN using the ring emission as seen in the ALMA 315 GHz image using the quadratic curve method (red), parametric ellipse fitting (blue), and least-squares fitting (green). The centers derived by fitting the ATCA radio ring continuum in Potter et al. (2009) and by fitting toHST image ring hot spots in Alp et al. (2018a) are also shown.
Download figure:
Standard imageHigh-resolution imageAppendix C: Simulating Potential Flux Loss
Simulations exploring the possibility of overresolving extended dust emission were performed withsimobserve andsimanalyze inCasa. The same antenna locations and integration times from the observations were used, to ensure consistency. The most extreme case we considered is a uniform ellipse spanning to the edge of the ring (24 across). This is similar to the largest angular scales these observations were sensitive to (see Table1). Even in this worst-case scenario, the differences in the integrated ring and ejecta flux densities between the simulated model and true observations were roughly 10% different when run through the same photometry prescription. A more realistic input model is shown in Figure19: a uniform broad annulus for the ring plus a fainter diffuse ellipse spanning the entire ejecta, plus more compact clumps; their integrated flux densities match within a few percent.
Figure 19. Casa simulations of ALMA observations at 315 GHz to test whether diffuse flux is lost in the observing process. Modeled emission that resembles the true image is recovered qualitatively and quantitatively and suggests that we are not missing a significant component of extended diffuse emission.
Download figure:
Standard imageHigh-resolution imageNotably, the uniform broad ellipse model does not translate to a simulated model that resembles the true observations; we therefore discount the possibility of the real submillimeter source being more diffuse and uniform as highly unlikely. By matching the input model more closely to the emission seen in the ALMA observations, we obtain reasonable-looking simulated maps of SN 1987A with a flux density difference of several percent. We conclude that there may be some missed extended flux, on the level of a few percent, if the true source has a slightly broader distribution. However, this is below the uncertainty level in the ALMA photometry.
An overly conservative approach to this issue would be to include an additional 10% systematic uncertainty in quadrature to photometry values for this effect. However, the uncertainty estimation presented in the main text already accounts somewhat for differences in faint features because it includes the standard deviation of flux densities from same-sized apertures placed at many random locations in each map.
