High angular resolution ALMA images of dust and molecules in the SN 1987A ejecta

PHIL CIGAN,¹ MIKAKO MATSUURA,¹ HALEY L. GOMEZ,¹ REMY INDEBETOUW,² FRAN ABELLÁN,³ MICHAEL GABLER,⁴

ANITA RICHARDS,⁵ DENNIS ALP,⁶ TIM DAVIS,¹ HANS-THOMAS JANKA,⁴ JASON SPYROMILIO,⁷ M. J. BARLOW,⁸ DAVID BURROWS,⁹

ELI DWEK,¹⁰ CLAES FRANSSON,¹¹ BRYAN GAENSLER,¹² JOSEFIN LARSSON,⁶ P. BOUCHET,¹³,¹⁴ PETER LUNDQVIST,¹¹

J. M. MARCAIDE,³ C.-Y. NG,¹⁵ SANGWOOK PARK,¹⁶ PAT ROCHE,¹⁷ JACCO TH. VAN LOON,¹⁸ J. C. WHEELER,¹⁹ AND

GIOVANNA ZANARDO²⁰

¹School of Physics and Astronomy, Cardiff University, Cardiff CF24 3AA, UK

²Department of Astronomy, University of Virginia, PO Box 400325, Charlottesville, VA 22904, USA

³Departamento de Astronomía y Astrofísica, Universidad de Valencia, C/Dr. Moliner 50, E-46100 Burjassot, Spain

⁴Max-Planck-Institut für Astrophysik, Karl-Schwarzchild-Straße 1, D-85748 Garching, Germany

⁵JBCA, School of Physics and Astronomy, University of Manchester, Manchester M13 9PL, UK

⁶Department of Physics, KTH Royal Institute of Technology, The Oskar Klein Centre, AlbaNova, SE-106 91 Stockholm, Sweden

⁷European Southern Observatory, Karl-Schwarzschild-Straße 2, 85748 Garching, Germany

⁸Department of Physics and Astronomy, University College London, Gower Street, London, WC1E 6BT, UK

⁹Department of Astronomy and Astrophysics, Pennsylvania State University, University Park, PA 16802, USA

¹⁰Observational Cosmology Lab., Code 665, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

¹¹Department of Astronomy, Stockholm University, The Oskar Klein Centre, AlbaNova, SE-106 91 Stockholm, Sweden

¹²Dunlap Institute for Astronomy and Astrophysics, 50 St. George Street, Toronto, ON M5S 3H4, Canada

¹³DRF/IRFU/DAp, CEA-Saclay, F-91191 Gif-sur-Yvette, France

¹⁴NRS/AIM, Université Paris Diderot, F-9119, Gif-sur-Yvette, France

¹⁵Department of Physics, The University of Hong Kong, Pokfulam Road, Hong Kong

¹⁶Department of Physics, University of Texas at Arlington, Arlington, TX 76019, USA

¹⁷Astrophysics, Department of Physics, University of Oxford, DWB, Keble Road, Oxford OX1 3RH, UK

¹⁸Lennard-Jones Laboratories, Keele University, ST5 5BG, UK

¹⁹Department of Astronomy, The University of Texas at Austin, 1 University Station C1400, Austin, TX 78712-0259, USA

²⁰International Centre for Radio Astronomy Research, M468, University of Western Australia, Crawley, WA 6009, Australia

(Received 12 July, 2019; Revised 16 September, 2019; Accepted 3 October, 2019)

Submitted to The Astrophysical Journal

ABSTRACT

We present high angular resolution (~80 mas) ALMA continuum images of the SN 1987A system, together with CO $J=2 \rightarrow 1$, $J=6 \rightarrow 5$, and SiO $J=5 \rightarrow 4$ to $J=7 \rightarrow 6$ images, which clearly resolve the ejecta (dust continuum and molecules) and ring (synchrotron continuum) components. Dust in the ejecta is asymmetric and clumpy, and overall the dust fills the spatial void seen in Hα images, filling that region with material from heavier elements. The dust clumps generally fill the space where CO $J=6 \rightarrow 5$ is fainter, tentatively indicating that these dust clumps and CO are locationally and chemically linked. In these regions, carbonaceous dust grains might have formed after dissociation of CO. The dust grains would have cooled by radiation, and subsequent collisions of grains with gas would also cool the gas, suppressing the CO $J=6 \rightarrow 5$ intensity. The data show a dust peak spatially coincident with the molecular hole seen in previous ALMA CO $J=2 \rightarrow 1$ and SiO $J=5 \rightarrow 4$ images. That dust peak, combined with CO and SiO line spectra, suggests that the dust and gas could be at higher temperatures than the surrounding material, though higher density cannot be totally excluded. One of the possibilities is that a compact source provides additional heat at that location. Fits to the far-infrared-millimeter spectral energy distribution give ejecta dust temperatures of 18–23K. We revise the ejecta dust mass to $M_{\text{dust}} = 0.2 - 0.4 M_{\odot}$ for carbon or silicate grains, or a maximum of $< 0.7 M_{\odot}$ for a mixture of grain species, using the predicted nucleosynthesis yields as an upper limit.
tion of a large mass of amorphous carbon dust. Though depending on the gas density, CO may be dissociated by electrons produced by radioactive decay (Clayton 2011), and/or (to a lesser extent) X-rays from the ring, however, the models of Sarangi & Cherchneff (2013) and Sarangi & Cherchneff (2015) indicate 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 which 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 two years 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 the Atacama Large Millimeter/Submillimeter Array (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 10352 and 10441 days after the initial explosion. The Band-7 (870µm) and 9 (450µm) integrations utilized between 34 and 36 antennae with baselines spanning 15m to 2.3km. See Table 1 for a summary of the observations.

Each data set was reduced separately with Common Astronomy Software Applications package (CASA<fnref n="1" />, McMullin et al. 2007), version 4.5.1. Once calibrated, the tclean algorithm was used to deconvolve and image the data.

The check source (reference quasar with precisely known position) and phase calibrator coordinates, determined with imfit in CASA, 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 June 2015 observations, with ~1.45 mm precipitable water vapor (PWV), and were therefore repeated in September 2015. 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 Table 2.

### 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⁻<fnref n="1" /> expansion velocity of the ejecta means that the linewidths span a substantial fraction of the observed bands. The continuum bands selected relative to the modelled molecular line emission are shown in Fig. 1, using the ALMA spectra and the emission line model of CO, SiO, SO, and SO₂ from Matsuura et al. (2017). Only windows that were free from molecular line emission (shown by the grey vertical bands) were used to make continuum images, centered at roughly 307, 315, and 679 GHz. The 315 GHz continuum image is shown in Fig. 2. 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 full-width at half-maximum (FWHM) of 0.30.

### 2.3. Molecular Line Data

In this section we present the molecular line data observed in the same blocks as the continuum discussed above: CO J=6→5 with rest frequency 691.47 GHz and SiO J=7→6 at $v_{rest}$ = 303.93 GHz. The 345.80 GHz CO J=3→2 and 347.33 GHz SiO J=8→7 lines were also

<footnote n="1">¹ http://casa.nrao.edu/</footnote>
<table><caption>**Table 1. Observations**</caption><thead><tr><th>Sub-band</th><th>Frequency Range (GHz)</th><th>Baselines (m)</th><th>Angular Scales (arcsec)</th><th>Observation Date</th><th>SN Day</th><th>Bandpass Calibrator</th><th>Phase Calibrator</th><th>Check Source</th><th>Time (min)</th></tr></thead><tbody><tr><td rowspan="6">B7</td><td>A1</td><td>299.88–315.87</td><td>45.4–1574.4</td><td>0.13–4.31</td><td>2015-06-28</td><td>10352</td><td>J0538-4405</td><td>J0635-7516</td><td>J0601-7036</td><td>18.4</td></tr><tr><td>A2</td><td>299.88–315.87</td><td>43.3–2269.9</td><td>0.09–4.52</td><td>2015-09-22</td><td>10438</td><td>J0538-4405</td><td>J0635-7516</td><td>J0601-7036</td><td>18.3</td></tr><tr><td>B</td><td>342.48–358.34</td><td>15.1–1574.4</td><td>0.11–11.46</td><td>2015-07-25</td><td>10379</td><td>J0538-4405</td><td>J0635-7516</td><td>J0601-7036</td><td>20.9</td></tr><tr><td>C</td><td>346.23–362.09</td><td>15.1–1574.4</td><td>0.11–11.34</td><td>2015-07-25</td><td>10379</td><td>J0538-4405</td><td>J0635-7516</td><td>J0601-7036</td><td>19.9</td></tr><tr><td>D1</td><td>303.62–319.48</td><td>45.4–1574.4</td><td>0.13–4.26</td><td>2015-06-28</td><td>10352</td><td>J0538-4405</td><td>J0635-7516</td><td>J0601-7036</td><td>18.8</td></tr><tr><td>D2</td><td>303.62–319.48</td><td>43.3–2269.9</td><td>0.09–4.47</td><td>2015-09-22</td><td>10438</td><td>J0538-4405</td><td>J0635-7516</td><td>J0601-7036</td><td>18.8</td></tr><tr><td rowspan="3">B9</td><td>A</td><td>673.44–681.06</td><td>43.3–2269.9</td><td>0.04–2.10</td><td>2015-09-25</td><td>10441</td><td>J0522-3627</td><td>J0601-7036</td><td>J0700-6610</td><td>12.5</td></tr><tr><td>B</td><td>680.94–688.56</td><td>43.3–2269.9</td><td>0.04–2.07</td><td>2015-09-25</td><td>10441</td><td>J0522-3627</td><td>J0601-7036</td><td>J0700-6610</td><td>12.5</td></tr><tr><td>C</td><td>688.44–696.06</td><td>43.3–2269.9</td><td>0.04–2.05</td><td>2015-09-25</td><td>10441</td><td>J0522-3627</td><td>J0700-6610</td><td>J0450-8101</td><td>12.5</td></tr></tbody></table>

NOTE— Observations for proposal ID 2013.1.00063. Each sub-band is comprised of four 2 GHz blocks of 128 channels (15.625 MHz each). The same flux calibrator, J0519-454, was used for all observation blocks.

<figure><img src="image_1.png" /><figcaption>**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 grey 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 Fig. 2). Data points indicate the flux densities in that band for the ring (green) and the ejecta (yellow).</figcaption></figure>

covered in these observing blocks, but as they are heavily blended we do not consider them in the present work.

The SiO *J*=7→6 and CO *J*=6→5 cubes were created with tclean in CASA, with a spectral resolution of 300 km s⁻¹, which gives a reasonable balance between velocity resolution and sensitivity per channel. CO *J*=6→5 was

imaged with natural weighting to maximize sensitivity per beam. SiO *J*=7→6 was imaged with robust= −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 Fig. 3.
<table><caption>**Table 2.** Imaging Properties</caption><thead><tr><th>Frequency Range<br>(GHz)</th><th>ν<sub>c</sub><br>(GHz)</th><th>Beam FWHM<br>(arcsec<sup>2</sup>)</th><th>Beam PA<br>(deg)</th><th>RMS Noise<br>(mJy bm<sup>-1</sup>)</th></tr></thead><tbody><tr><td><strong>Continuum</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>224.00 – 227.00</td><td>225.50</td><td>0.30 × 0.30</td><td>0.00</td><td>0.12</td></tr><tr><td>238.00 – 243.00</td><td>240.50</td><td>0.30 × 0.30</td><td>0.00</td><td>0.09</td></tr><tr><td>246.00 – 249.00</td><td>247.50</td><td>0.30 × 0.30</td><td>0.00</td><td>0.10</td></tr><tr><td>269.50 – 270.50</td><td>270.00</td><td>0.30 × 0.30</td><td>0.00</td><td>0.18</td></tr><tr><td>278.00 – 280.00</td><td>279.00</td><td>0.30 × 0.30</td><td>0.00</td><td>0.10</td></tr><tr><td>306.06 – 307.47</td><td>306.76</td><td>0.20 × 0.15</td><td>124.32</td><td>0.07</td></tr><tr><td>311.88 – 319.48</td><td>315.68</td><td>0.19 × 0.14</td><td>119.11</td><td>0.04</td></tr><tr><td>673.45 – 685.00</td><td>679.22</td><td>0.08 × 0.06</td><td>74.37</td><td>0.71</td></tr><tr><td><strong>Spectral Lines</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>CO J=2→1</td><td>230.54</td><td>0.06 × 0.04</td><td>27.43</td><td>0.04</td></tr><tr><td>CO J=6→5</td><td>691.47</td><td>0.09 × 0.07</td><td>185.35</td><td>2.82</td></tr><tr><td>SiO J=5→4</td><td>217.10</td><td>0.06 × 0.04</td><td>19.74</td><td>0.05</td></tr><tr><td>SiO J=6→5</td><td>260.52</td><td>0.04 × 0.03</td><td>173.66</td><td>0.06</td></tr><tr><td>SiO J=7→6</td><td>303.93</td><td>0.13 × 0.10</td><td>35.47</td><td>0.47</td></tr></tbody></table>

NOTE— The position angles are counter-clockwise from north. CO J=2→1, SiO J=5→4, and SiO J=6→5 parameters are for the data cubes from Abellán et al. (2017). The $\nu_c$ values listed for the CO and SiO lines are rest frequencies. RMS values for the spectral lines are per velocity channel. For observation dates and epochs, see Table 1.

In addition to the molecular line data described above, we also utilize the CO J=2→1, SiO J=5→4, and SiO J=6→5² 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 signal to noise as they were combined with Cycle 3 data, and due to CO J=6→5 being in a band with poorer atmospheric transmission than CO J=2→1. The additional Cycle 3 data also give their CO J=2→1, SiO J=5→4, and SiO J=6→5 maps finer spatial resolution than the observations presented in the current work, with FWHM 0.04–0.06 (see Table 2). For full details of their data reduction technique, we refer the reader to their Section 2. The channel maps are shown in Figs. A.1 and A.2. 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⁻¹ 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\alpha$ in Fig. 2. 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

Appendix B for details). Our derived ring+ejecta system center used in this work is $\alpha=5^{h}35^{m}27.998$, $\delta=-69^{\circ}16'11.1''107$ (ICRS), $\pm 18$ mas (Fig. B.1). At 315 GHz, the ring is clumpy and the brightness contrast in the east and west components of the ring is different to that observed in the $H\alpha$ 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 § 4.3), and the inner region is thermal dust emission from the SN ejecta (Indebetouw et al. 2014).

Fig. 3 shows an enlarged view of the ejecta images of dust continuum and lines. The majority of the submm 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 lower robust settings in tclean, 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\alpha$ image, like a ‘lock in the keyhole’. This is shown in Fig. 2, where the $3\sigma$ contours highlighting the major continuum features are overlaid onto the continuum and $H\alpha$ images. The alignment accuracy is ~1 pixel in the images, given the astrometric uncertainties discussed in § 2.1 (12 mas for Band-7 continuum) and Appendix B (6 mas for registration of the Hubble 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 Fig. 3). 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 (Panagia 1999), the Band-9 beam FWHM of 0.08 × 0.06 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 Fig. 3 and the dust emission (in red) in Fig. 4 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 Fig. 2 – is consistent with noise.

The molecular images provide a probe of different conditions in the ejecta, where lower transitions probe lower temperature gas (if optically thin, see Section 5). One prominent feature is the central hole seen in the CO J=2→1 and SiO J=5→4 images (middle left and lower left panels of Fig. 3). This was first reported by Abellán et al. (2017) and was seen both in the integrated 2D spatial maps and the 3D data cubes. Although the integrated SiO J=6→5 map (mid-

<footnote n="2">² The images for the SiO J=6→5 transition were described but not shown in Abellán et al. (2017).</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 2.** ALMA 315 GHz (with beam) and 2014 *HST* 625W band image (Fransson et al. 2015), which includes $H\alpha$. The yellow contours display 315 GHz emission at 0.2 mJy/beam. 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) and *HST* image registration based on fitting the ring (6 mas) is of order 1 pixel in these images.</figcaption></figure>

The panel of Fig. 3) does not show the hole clearly in the same manner as SiO $J=5 \to 4$ and CO $J=2 \to 1$, the hole is also visible in the central channels ($v = 0 - 300$ km s$^{-1}$) of the velocity map (Fig. A.2). Because of the additional −600–0 km s$^{-1}$ components located within the same line of sight as the hole (Fig. A.2) in the integrated maps, the hole is not clear in the SiO $J=6 \to 5$ map. The CO and SiO molecular hole is just to the south of the ‘keyhole’ that is seen in $H\alpha$ (Fig. 8 of Fransson et al. (2015); top right panel of our Fig. 3), though the molecular hole appears to be slightly smaller in scale and located on the southern edge of the hole in $H\alpha$ emission. The centers of the holes are offset by ~50 mas, or ~4× the astrometric and alignment errors.

CO $J=2 \to 1$ and SiO $J=5 \to 4$ have similar structures in the integrated images, however the spatial distributions of the higher transitions of each species have some differences. SiO $J=6 \to 5$ is more evenly distributed in a shell pattern while the lower S/N image of CO $J=6 \to 5$ appears clumpy (Fig. 3), though this is likely affected by the noise.

CO $J=6 \to 5$ has emission coincident with the CO $J=2 \to 1$ hole, in that its channel maps (Fig. A.1) show emission around the hole location, albeit at low S/N. However, the integrated spatial distribution appears different from CO $J=2 \to 1$. The brightness peaks are distributed differently, and the hole is not visible in the integrated CO $J=6 \to 5$

map due to some emission at those coordinates in the 600–900 km s$^{-1}$ channels (the far side). The presence of a molecular hole in SiO $J=7 \to 6$ cannot be confirmed in these data, as the systemic line center ($v_{LSRK} \sim 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 SiO $J=7 \to 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 Fig. A.2).

The resolved dust peak (small $5\sigma$ contour in Fig. 3) is co-located 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 SiO $J=5 \to 4$ channels near the systemic velocity. The brightest points of dust emission tend to coincide with relative depressions in the CO $J=6 \to 5$ brightness, giving the appearance of an anti-correlation between the main dust and CO $J=6 \to 5$ features. This is more clearly demonstrated in Fig. 4, where the dust (red) and CO $J=6 \to 5$ (blue) images are overlaid. The individual images were normalized independently to emphasize the main
<figure><img src="image_3.png" /></figure>

**Figure 3.** Multi-wavelength view of the ejecta in SN 1987A. The cyan contours represent the 679 GHz dust emission at $3\sigma$ and $5\sigma$. Beam sizes for individual maps are denoted by the green ellipses. The small cross denotes the system center as defined in Appendix B. The bottom right panel is a 3-color image of CO $J=2 \to 1$ in red, SiO $J=5 \to 4$ in green, and SiO $J=7 \to 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 cophsat, and none matches the distribution of the 679 GHz dust. Some of the CO falls in the H$\alpha$ hole (the lower left), but the majority of the CO peaks on the little H$\alpha$ ‘wing’ to the right of the hole. The small $5\sigma$ cyan contour just northeast of the center of the remnant is the so-called “blob”.

features of each, with the visible colors shown roughly cor-
responding 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 CO J=6→5 peaks
do not generally overlap. Rather, the brightest dust features
are located in areas of relatively faint CO J=6→5 emission
and vice-versa.

To test whether the apparent dust-CO anti-correlation is
an artifact or result of the data reduction or continuum sub-
<figure><img src="image_4.png" /><figcaption>**Figure 4.** The spatial anti-correlation of dust and CO *J*=6→5. In this overlay, CO *J*=6→5 is in blue and 679 GHz is in red, showing their relative spatial distributions. The off-white cross denotes the system center position as defined in Appendix B. The gold line highlights the extent of the major features in the dust, while the teal line demarcates the major CO *J*=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 CO *J*=6→5, respectively. The dust and CO *J*=6→5 emission exhibit a notable anti-correlation.</figcaption></figure>

traction, we performed several checks. First, the Band-
9 dust continuum was reconstructed in different ways, by
imaging (CASA mfs-mode) in a variety of spectral win-
dows and also by making a data cube across the entirety
of Band-9 (including the CO line), and fitting the contin-
uum emission. These techniques gave consistent results. Ini-
tially we used CASA 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 dif-
ferences. This means that the structures seen in our final
CO *J*=6→5 map are robust to variations in how the con-
tinuum is determined and subtracted. The anti-correlation in
CO *J*=6→5 is visible, even before continuum subtraction, in
the CO *J*=6→5 dirty map (i.e., with no cleaning to decon-
volve the interferometer sidelobes). Thus the apparent anti-
correlation seen in the dust and CO distributions is robust
to changes in the data processing. Lastly, we test whether
the anti-correlation is statistically robust by calculating the
weighted version of the normalized cross-correlation func-
tion $\sum \text{cov}_{XY}/\sqrt{\text{cov}_{XX} \cdot \text{cov}_{YY}}$, which returns a standard
correlation measure *r* between the range -1 and +1. The pixel
weights used were the map (S/N)$^2$. The resulting correlation
measure is *r*=-0.30 ± 0.08 – a moderate anti-correlation –
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 in-
vestigating different angular scales using the wavelet analysis
described in Frick et al. (2001) and Arshakian & Ossenkopf
(2016). On scales of 1–2 beam widths (i.e., convolutions
with kernels of those scales), the images start to become in-
creasingly positively correlated, which is expected due to the
peaks being separated by that amount (~5× the astrometry
error). Below these scales, *r* remains negative, so the anti-
correlation is not sensitive to small changes in image resolu-
tion. Using the same wavelet analysis on the other images,
the dust correlates more positively with the other CO and SiO
lines than with CO *J*=6→5. The standard correlation mea-
sures agree, with *r*=+0.04 for CO *J*=2→1 and *r*=+0.36 for
SiO *J*=6→5.

The brightest dust feature is located one beam width north-
east of the secondary CO *J*=6→5 peak, and the brightest
CO *J*=6→5 features curve around the main dust peak. The
CO peaks are obvious because the line emission is brighter
than the continuum (c.f. the integrated profile in Fig. 1), but
there is some fainter (S/N<5) CO *J*=6→5 emission overlap-
ping 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 SiO *J*=7→6 is co-
incident with the dust peak and the SiO *J*=5→4 hole (as seen
in the bottom right panel of Fig. 3). At the southern edge of
the hole, some faint Hα emission appears to be aligned with
CO *J*=2→1 in projection, but their velocity ranges differ.
The 3-D view of Hα (Larsson et al. 2016, their Figure 6) in-
dicates that the Hα emission in this region peaks at velocities
around –1500 km s⁻¹ while the peak CO *J*=2→1 is between
0 and +200 km s⁻¹. While no velocity information is avail-
able for the dust continuum emission, it is spatially offset
from the nearby Hα and CO *J*=2→1 by ~1 dust resolution
element. That is, the CO *J*=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 copspatial with the peaks of the
ejecta dust emission (with anti-correlated CO *J*=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 emis-
sion observed in the SiO *J*=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 ther-
mal IR greybody emission from dust (from the ejecta re-
gion - Matsuura et al. 2011; Indebetouw et al. 2014), non-
thermal radio/mm synchrotron emission (from the ring -
Manchester 2007; Potter et al. 2009; Zanardo et al. 2010; Lakićević et al. 2012; Zanardo et al. 2013; Indebetouw et al. 2014) and a lesser contribution from free-free mm/submm 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 Fig. 1 and §2.2, with the chosen frequency windows summarized in Table 3; 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 Section 3 (Appendix B), 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 (Table 3). For the 315 GHz ejecta, this results in an elliptical aperture with semi-major and semi-minor axes of $r_{MAJ} \times r_{MIN} = 0.4'42 \times 0.3'6$ and major axis P.A. of 25° (N through E). For the 315 GHz circumstellar ring, an elliptical annulus with inner $r_{MAJ}$ and $r_{MIN}$ of $0.4'45 \times 0.4'42$, outer $r_{MAJ}$ and $r_{MIN}$ of $1.4'45 \times 1.3'5$, 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 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.

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, $\sigma_{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 $0.1'$, a size which 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 & 7, and 20% in Band-9³. Finally, we considered the possibility that we may be missing diffuse emission from cold dust within the SN structure due to over-resolving 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 Appendix C for more details.

The ejecta, ring and total system flux densities and uncertainties are shown in Fig. 5 as gold circles, green rings, and purple diamonds, respectively, and their values are listed in Table 3. Previous measurements are also shown in Fig. 5 for reference. Preliminary Cycle 6 flux densities (Matsuura et al., 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 mJy and 18.1 mJy, respectively. The uncertainty on each of these flux density measurements is estimated as 0.4 mJy.

We 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 § 2.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 § 2.1, weather affected the phase stability of observations between 346–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 Section 3, 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 grey hexagons, and represent the overall spectral energy distribution (SED) of the system. The total emission at 1.4 GHz, 18 GHz, 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, & 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 submm

<footnote n="3">³ ALMA Cycle 2 Technical Handbook,
https://almascience.eso.org/documents-and-tools/cycle-2/alma-technical-handbook/</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Continuum values from 225–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 10402, and Cycle 6 is day 11522. The higher angular resolution ALMA data confirm the ejecta emission follows a thermal dust profile down to ~200 GHz. The ring emission, on the other hand, shows no evidence of submm dust emission but instead is consistent with a synchrotron emission profile – ATCA observations from Days 9280–9686 (Zanardo et al. 2014; Callingham et al. 2016) combined with the Cycle 2 ALMA ring flux densities give a power-law index of $\alpha = -0.70 \pm 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 (Matsuura et al., in preparation) are shown as the orange crosses and light purple hexagons for the ejecta and total emission, respectively. The positive error bars for the Cycle 2 Band 6&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 empty grey hexagons are measures of the total system flux density in various parts of the spectrum at day 9280: 1.4 GHz, 18 GHz, 44 GHz (Zanardo et al. 2013); 9 GHz (Ng et al. 2013); 94 GHz (Lakićević et al. 2012); 102, 213, 345, & 672 GHz (Zanardo et al. 2014). The empty black hexagons represent the ATCA total system flux densities between 1–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 10402 due to the ALMA band 7 systematic error discussed in Section 4.1. The brown asterisks are the unresolved Herschel 70–500 µm flux densities (Matsuura et al. 2015), and the crimson asterisks are the 2010 HERITAGE flux densities (Meixner et al. 2013; Matsuura et al. 2011). The ring emission model S(ν,t) from Cendes et al. (2018) for day 10402 and $\alpha = -0.70$ is shown by the dashed grey line.</figcaption></figure>
<table><caption>**Table 3. Continuum Photometry**</caption><thead><tr><th>ν<sub>c</sub><br>(GHz)</th><th>Δν<br>(GHz)</th><th>S<sub>ν</sub><br>(mJy)</th><th colspan="2">Aperture Center</th><th>R<sub>MAJ</sub><sup>†</sup><br>(arcsec)</th><th>R<sub>MIN</sub><sup>†</sup><br>(arcsec)</th><th>P.A.<br>(deg)</th></tr><tr><th></th><th></th><th></th><th>RA (deg)</th><th>DEC (deg)</th><th></th><th></th><th></th></tr></thead><tbody><tr><td colspan="8" style="text-align: center; font-weight: bold;">Ejecta</td></tr><tr><td>225.50</td><td>3.0</td><td>0.5 ± 0.2 <sup>+0.3</sup><sub>-0.0</sub></td><td>83.866586</td><td>-69.269733</td><td>0.32</td><td>0.29</td><td>125</td></tr><tr><td>240.50</td><td>5.0</td><td>0.7 ± 0.2 <sup>+0.4</sup><sub>-0.1</sub></td><td>83.866600</td><td>-69.269720</td><td>0.37</td><td>0.36</td><td>75</td></tr><tr><td>247.50</td><td>3.0</td><td>1.1 ± 0.4 <sup>+0.6</sup><sub>-0.1</sub></td><td>83.866639</td><td>-69.269739</td><td>0.41</td><td>0.35</td><td>125</td></tr><tr><td>270.00</td><td>1.0</td><td>1.2 ± 0.6 <sup>+0.7</sup><sub>-0.1</sub></td><td>83.866652</td><td>-69.269709</td><td>0.40</td><td>0.31</td><td>60</td></tr><tr><td>279.00</td><td>2.0</td><td>1.3 ± 0.4 <sup>+0.8</sup><sub>-0.1</sub></td><td>83.866662</td><td>-69.269747</td><td>0.41</td><td>0.40</td><td>125</td></tr><tr><td>306.76</td><td>1.4</td><td>1.8 ± 0.6 <sup>+1.1</sup><sub>-0.2</sub></td><td>83.866667</td><td>-69.269753</td><td>0.42</td><td>0.36</td><td>125</td></tr><tr><td>315.68</td><td>7.6</td><td>1.8 ± 0.6 <sup>+1.1</sup><sub>-0.2</sub></td><td>83.866667</td><td>-69.269753</td><td>0.42</td><td>0.36</td><td>125</td></tr><tr><td>679.22</td><td>11.5</td><td>36.2 ± 7.2 <sup>+7.2</sup><sub>-7.2</sub></td><td>83.866728</td><td>-69.269740</td><td>0.42</td><td>0.36</td><td>125</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Ring</td></tr><tr><td>225.50</td><td>3.0</td><td>13.0 ± 0.5 <sup>+7.8</sup><sub>-1.3</sub></td><td>83.866585</td><td>-69.269731</td><td>1.45, 0.35</td><td>1.35, 0.33</td><td>175</td></tr><tr><td>240.50</td><td>5.0</td><td>12.5 ± 0.5 <sup>+7.5</sup><sub>-1.3</sub></td><td>83.866600</td><td>-69.269722</td><td>1.45, 0.42</td><td>1.35, 0.39</td><td>175</td></tr><tr><td>247.50</td><td>3.0</td><td>13.4 ± 0.8 <sup>+8.1</sup><sub>-1.3</sub></td><td>83.866630</td><td>-69.269736</td><td>1.45, 0.43</td><td>1.35, 0.40</td><td>175</td></tr><tr><td>270.00</td><td>1.0</td><td>13.7 ± 1.1 <sup>+8.2</sup><sub>-1.4</sub></td><td>83.866592</td><td>-69.269720</td><td>1.45, 0.45</td><td>1.35, 0.42</td><td>175</td></tr><tr><td>279.00</td><td>2.0</td><td>12.9 ± 0.8 <sup>+7.7</sup><sub>-1.3</sub></td><td>83.866661</td><td>-69.269747</td><td>1.45, 0.44</td><td>1.35, 0.41</td><td>175</td></tr><tr><td>306.76</td><td>1.4</td><td>12.0 ± 1.5 <sup>+7.2</sup><sub>-1.2</sub></td><td>83.866666</td><td>-69.269753</td><td>1.45, 0.45</td><td>1.35, 0.42</td><td>175</td></tr><tr><td>315.68</td><td>7.6</td><td>11.3 ± 1.7 <sup>+6.8</sup><sub>-1.1</sub></td><td>83.866667</td><td>-69.269753</td><td>1.45, 0.45</td><td>1.35, 0.42</td><td>175</td></tr><tr><td>679.22</td><td>11.5</td><td>19.4 ± 19.3 <sup>+3.9</sup><sub>-3.9</sub></td><td>83.866720</td><td>-69.269737</td><td>1.45, 0.45</td><td>1.35, 0.42</td><td>175</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Total System</td></tr><tr><td>225.50</td><td>3.0</td><td>13.5 ± 0.5 <sup>+8.1</sup><sub>-1.3</sub></td><td>83.866585</td><td>-69.269731</td><td>1.45</td><td>1.35</td><td>175</td></tr><tr><td>240.50</td><td>5.0</td><td>13.2 ± 0.5 <sup>+7.9</sup><sub>-1.3</sub></td><td>83.866600</td><td>-69.269722</td><td>1.45</td><td>1.35</td><td>175</td></tr><tr><td>247.50</td><td>3.0</td><td>14.6 ± 0.7 <sup>+8.7</sup><sub>-1.5</sub></td><td>83.866630</td><td>-69.269736</td><td>1.45</td><td>1.35</td><td>175</td></tr><tr><td>270.00</td><td>1.0</td><td>15.4 ± 0.8 <sup>+9.2</sup><sub>-1.5</sub></td><td>83.866592</td><td>-69.269720</td><td>1.45</td><td>1.35</td><td>175</td></tr><tr><td>279.00</td><td>2.0</td><td>14.2 ± 0.8 <sup>+8.5</sup><sub>-1.4</sub></td><td>83.866661</td><td>-69.269747</td><td>1.45</td><td>1.35</td><td>175</td></tr><tr><td>306.76</td><td>1.4</td><td>13.9 ± 1.7 <sup>+8.3</sup><sub>-1.4</sub></td><td>83.866666</td><td>-69.269753</td><td>1.45</td><td>1.35</td><td>175</td></tr><tr><td>315.68</td><td>7.6</td><td>13.3 ± 1.8 <sup>+8.0</sup><sub>-1.3</sub></td><td>83.866667</td><td>-69.269753</td><td>1.45</td><td>1.35</td><td>175</td></tr><tr><td>679.22</td><td>11.5</td><td>55.5 ± 20.9 <sup>+11.1</sup><sub>-11.1</sub></td><td>83.866720</td><td>-69.269737</td><td>1.45</td><td>1.35</td><td>175</td></tr></tbody></table>

NOTE—Integrated flux densities of the ejecta, ring, and total system for each continuum band with central frequency $\nu_c$ and bandwidth $\Delta\nu_c$. Integrated flux densities are quoted as *value* ± *measurement uncertainty* ± *systematic uncertainty*. The systematic uncertainty includes calibration uncertainties (10% in Bands 6&7 or 20% in Band 9), and an additional 50% on the positive side in Bands 6&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.

† The ring annulus radii are given as ($R_{outer}$, $R_{inner}$).
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

Fig. 5 displays the mm to FIR SED. In this figure, the brown asterisks show the FIR flux densities measured by Herschel for the total SN 1987A (unresolved) system, and the gold circles show the mm 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 mm, 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 the Herschel flux densities from Matsuura et al. (2015) in our model fits since they are measuring the emission from the ejecta dust, albeit unresolved. Two Herschel flux densities are treated as upper limits: 70 µm, as it is possibly contaminated by warm ring dust (Matsuura et al. 2019) and/or [O I] 63 µm emission; and 500 µm, as it was a non-detection. One potential issue with using the Herschel flux densities is that the Herschel 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 from $^{44}$Ti decay, which has an estimated lifetime of 85 years (Ahmad et al. 2006; Jerkstrand et al. 2011, also see later Matsuura et al. (2011)). The predicted decrease in this decay energy between the 2012 Herschel and 2015 ALMA observations is 4.2%. Assuming the FIR luminosity decreased by this amount between the 2012 and 2015 epochs, the reduction in the temperature of a 20K blackbody would be ~0.2K, translating into individual Herschel 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 by $^{44}$Ti decay, the use of the Herschel flux densities with the latest ALMA measurements is valid. An alternative additional heating source will be discussed in Section 6.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 the `emcee` 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&7 flux densities from other cycles (see § 4.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 form $S_ν(λ) = M_{dust}κ_{abs}(λ)B_ν(λ)/d^2$, where $S_ν(λ)$ is flux density, $κ_{abs}(λ)$ is the mass absorption coefficient of dust grains, $B_ν(λ)$ is the Planck function and $d$ 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 radius $a$ and density $ρ$, $κ_{abs}(λ)$ is defined as $κ = 3/4 Qρa$, where $Q$ 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 $κ_0 (λ/λ_0)^{-β}$ (where $κ_0$ is the reference value of $κ_{abs}(λ)$ at wavelength $λ_0$ and $β$ is the power law emissivity index of $κ_{abs}(λ)$).

Fitting a dust greybody to the FIR-mm SED with the power law approximation to $κ_{abs}$ gives fit parameters of $β = 2.05_{-0.10}^{+0.11}$, $M_{dust} = 1.53_{-0.13}^{+0.13} M_⊙$, and $T_d = 17.83_{-0.57}^{+0.60}$ K if $κ_0(850µm) = 0.07 m^2 kg^{-1}$ (using the empirical measurement of $κ_{abs}$ assuming 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 greybody fits to the FIR-mm 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: 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 30K, 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, Zubko et al. 1996; Rouleau & Martin 1991), graphite (Draine & Lee 1984), ‘cosmic silicates’ (amorphous FeMgSiO$_4$, Jaeger et al. 1994), other silicates including enstatite and forsterite from Henning & Stognienko (1996), pure iron (Henning & Stognienko 1996), 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 µm 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 submm–mm), from the last 300–500 µm of each curve to ensure a smooth continuation of the general trend in each. We also consider a “continuous distribution of ellipsoids” (CDE) model for a composite of carbons and silicates (Dwek & Arendt 2015), with a carbon volume filling factor $f_C$ of 18%.

For Mie theory calculations, we adopt a grain radius of $a = 0.01$ µm for PAHs, and $a = 0.1$ µm for all other models. For most dust models, Mie-derived FIR $κ_{abs}$ curves are
<table><caption>**Table 4.** Modified Blackbody Fits</caption><thead><tr><th>Dust</th><th>Reference</th><th>Grain Density<br>(g cm<sup>-3</sup>)</th><th>M<sub>dust</sub><br>(M<sub>dust</sub>)</th><th>T<br>(K)</th><th>κ<sub>850</sub><br>(m<sup>2</sup> kg<sup>-1</sup>)</th><th>Good Fit</th></tr></thead><tbody><tr><td>amC (ACH2 sample), Mie 0.1µm</td><td>Zubko et al. (1996)</td><td>1.81</td><td>1.46<sup>+0.09</sup><sub>-0.08</sub></td><td>17.5<sup>+0.1</sup><sub>-0.1</sub></td><td>0.087</td><td>Y</td></tr><tr><td>amC (ACAR sample), Mie 0.1µm</td><td>Zubko et al. (1996)</td><td>1.81</td><td>0.38<sup>+0.02</sup><sub>-0.02</sub></td><td>22.0<sup>+0.2</sup><sub>-0.2</sub></td><td>0.254</td><td>N</td></tr><tr><td>amC (BE sample), Mie 0.1µm</td><td>Zubko et al. (1996)</td><td>1.81</td><td>0.77<sup>+0.05</sup><sub>-0.04</sub></td><td>20.7<sup>+0.2</sup><sub>-0.2</sub></td><td>0.141</td><td>N</td></tr><tr><td>amC (AC1 sample), Mie 0.1µm</td><td>Rouleau & Martin (1991)</td><td>1.85</td><td>0.43<sup>+0.03</sup><sub>-0.03</sub></td><td>21.6<sup>+0.3</sup><sub>-0.3</sub></td><td>0.203</td><td>Y</td></tr><tr><td>Cellulose (800K sample), Mie 0.1µm</td><td>Jager et al. (1998)</td><td>1.81</td><td>0.46<sup>+0.03</sup><sub>-0.03</sub></td><td>18.8<sup>+0.2</sup><sub>-0.2</sub></td><td>0.178</td><td>Y</td></tr><tr><td>Graphite, Mie 0.1µm</td><td>Draine & Lee (1984)</td><td>2.26</td><td>1.62<sup>+0.11</sup><sub>-0.10</sub></td><td>17.8<sup>+0.2</sup><sub>-0.2</sub></td><td>0.069</td><td>Y</td></tr><tr><td>PAH (neutral), Mie 0.01µm</td><td>Laor & Draine (1993)</td><td>2.24</td><td>1.69<sup>+0.11</sup><sub>-0.10</sub></td><td>18.0<sup>+0.2</sup><sub>-0.2</sub></td><td>0.071</td><td>Y</td></tr><tr><td>Silicate – Enstatite, Mie 0.1µm</td><td>Jäger et al. (2003)</td><td>2.71</td><td>4.10<sup>+0.23</sup><sub>-0.24</sub></td><td>18.0<sup>+0.2</sup><sub>-0.2</sub></td><td>0.029</td><td>Y</td></tr><tr><td>Silicate – Forsterite, Mie 0.1µm</td><td>Jäger et al. (2003)</td><td>3.2</td><td>4.03<sup>+0.26</sup><sub>-0.25</sub></td><td>17.9<sup>+0.2</sup><sub>-0.2</sub></td><td>0.029</td><td>Y</td></tr><tr><td>Silicate – “Cosmic”, Mie 0.1µm</td><td>Jaeger et al. (1994)</td><td>3.2</td><td>3.46<sup>+0.24</sup><sub>-0.22</sub></td><td>17.7<sup>+0.2</sup><sub>-0.2</sub></td><td>0.034</td><td>Y</td></tr><tr><td>Silicate/Carbon Composite CDE – f<sub>C</sub>=0.18</td><td>Dwek & Arendt (2015)</td><td>2.95</td><td>0.38<sup>+0.02</sup><sub>-0.02</sub></td><td>21.1<sup>+0.2</sup><sub>-0.2</sub></td><td>0.270</td><td>N</td></tr><tr><td>Silicate – LMC Average</td><td>Weingartner & Draine (2001)</td><td>...</td><td>2.49<sup>+0.15</sup><sub>-0.14</sub></td><td>18.0<sup>+0.2</sup><sub>-0.2</sub></td><td>0.047</td><td>Y</td></tr><tr><td>Silicate – 30K Average</td><td>Demyk et al. (2017)</td><td>...</td><td>0.27<sup>+0.02</sup><sub>-0.02</sub></td><td>18.5<sup>+0.2</sup><sub>-0.2</sub></td><td>0.265</td><td>Y</td></tr><tr><td>Silicate – Composite Aggregate</td><td>Semenov et al. (2003)</td><td>...</td><td>35.06<sup>+2.18</sup><sub>-2.18</sub></td><td>21.4<sup>+0.2</sup><sub>-0.2</sub></td><td>0.003</td><td>Y</td></tr><tr><td>Silicate – Porous Multilayer Spheres</td><td>Semenov et al. (2003)</td><td>...</td><td>6.13<sup>+0.36</sup><sub>-0.37</sub></td><td>29.4<sup>+0.4</sup><sub>-0.4</sub></td><td>0.010</td><td>N</td></tr><tr><td>Silicate – Bare Grains, 0.03 Myr</td><td>Ormel et al. (2011)</td><td>...</td><td>0.50<sup>+0.03</sup><sub>-0.03</sub></td><td>20.7<sup>+0.3</sup><sub>-0.2</sub></td><td>0.180</td><td>Y</td></tr><tr><td>Silicate – Icy Grains, 0.03 Myr</td><td>Ormel et al. (2011)</td><td>...</td><td>0.60<sup>+0.04</sup><sub>-0.04</sub></td><td>18.3<sup>+0.2</sup><sub>-0.2</sub></td><td>0.184</td><td>Y</td></tr><tr><td>Silicate – Naked Grains, n<sub>H</sub> = 10<sup>5</sup></td><td>Ossenkopf & Henning (1994)</td><td>...</td><td>0.64<sup>+0.04</sup><sub>-0.04</sub></td><td>20.9<sup>+0.2</sup><sub>-0.2</sub></td><td>0.163</td><td>N</td></tr><tr><td>Silicate – Thin Ice Mantles, n<sub>H</sub> = 10<sup>5</sup></td><td>Ossenkopf & Henning (1994)</td><td>...</td><td>0.74<sup>+0.04</sup><sub>-0.04</sub></td><td>19.1<sup>+0.2</sup><sub>-0.2</sub></td><td>0.142</td><td>Y</td></tr><tr><td>SiC, Mie 0.1µm</td><td>Pegourie (1988)</td><td>3.22</td><td>1.57<sup>+0.08</sup><sub>-0.09</sub></td><td>23.3<sup>+0.2</sup><sub>-0.2</sub></td><td>0.058</td><td>N</td></tr><tr><td>FeS, Mie 0.1µm</td><td>Henning & Stognienko (1996)</td><td>4.83</td><td>0.68<sup>+0.04</sup><sub>-0.04</sub></td><td>34.7<sup>+0.6</sup><sub>-0.6</sub></td><td>0.086</td><td>N</td></tr><tr><td>FeO, Mie 0.1µm</td><td>Henning et al. (1995)</td><td>5.7</td><td>0.28<sup>+0.01</sup><sub>-0.02</sub></td><td>28.2<sup>+0.4</sup><sub>-0.3</sub></td><td>0.259</td><td>N</td></tr><tr><td>SiO<sub>2</sub>, Mie 0.1µm</td><td>Henning & Mutschke (1997)</td><td>2.196</td><td>4.41<sup>+0.39</sup><sub>-0.33</sub></td><td>17.5<sup>+0.2</sup><sub>-0.2</sub></td><td>0.022</td><td>N</td></tr><tr><td>TiO<sub>2</sub>, Mie 0.1µm</td><td>Posch et al. (2003)</td><td>3.78</td><td>81.77<sup>+5.10</sup><sub>-4.49</sub></td><td>17.7<sup>+0.2</sup><sub>-0.2</sub></td><td>0.001</td><td>Y</td></tr><tr><td>Al<sub>2</sub>O<sub>3</sub> “Compact” sample, Mie 0.1µm</td><td>Begemann et al. (1997)</td><td>3.2</td><td>0.90<sup>+0.06</sup><sub>-0.06</sub></td><td>19.0<sup>+0.2</sup><sub>-0.2</sub></td><td>0.112</td><td>Y</td></tr><tr><td>NaAlSi<sub>2</sub>O<sub>6</sub>, Mie 0.1µm</td><td>Mutschke et al. (1998)</td><td>2.4</td><td>0.10<sup>+0.01</sup><sub>-0.01</sub></td><td>23.1<sup>+0.3</sup><sub>-0.3</sub></td><td>0.982</td><td>Y</td></tr><tr><td>MgAl<sub>2</sub>O<sub>4</sub>, Mie 0.1µm</td><td>Fabian et al. (2001)</td><td>3.64</td><td>121.74<sup>+7.54</sup><sub>-7.48</sub></td><td>17.7<sup>+0.2</sup><sub>-0.2</sub></td><td>0.001</td><td>Y</td></tr><tr><td>Pure Iron, Mie 0.1µm</td><td>Henning & Stognienko (1996)</td><td>7.87</td><td>3.97<sup>+0.28</sup><sub>-0.23</sub></td><td>19.9<sup>+0.2</sup><sub>-0.2</sub></td><td>0.025</td><td>Y</td></tr></tbody></table>

NOTE— Mass and temperature fits for greybodies 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 radius $a = 0.1µm$ were assumed except for the case of PAHs, where $a = 0.01µm$ was used. Grain densities are given for the Mie and CDE cases: amorphous carbon from Zubko et al. (2004) (and Rouleau & Martin (1991) 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 factor $f_C=18\%$ from Dwek & Arendt (2015), jadeite from Mutschke et al. (1998), and in general from the Jena optical constants databaseª (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 $χ^2 < 2$.

<footnote>ªhttps://www.astro.uni-jena.de/Laboratory/OCDB/index.html</footnote>
nearly identical for grain radii of $a < 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 Table 4.

To determine whether the SED fit is ‘good’, we set a quality threshold whereby $\chi^2_{\nu} < 2$. 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 $\chi^2_{\nu} < 2$ even if it falls below all of the Herschel points due to the nature of the $\chi^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 § 6.1.

### 4.2.2. Results of the dust fits

Fig. 6 shows 28 of the resulting dust fits to the FIR-submm SED (top panel of the figure) using the various dust emissivity curves (bottom panel of the figure) listed in Table 4. These have been grouped into similar dust compositions: amC, graphite, PAHs, silicates, oxides, carbides, sulphides, 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 Table 4.) Fig. 6 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–23K (Table 4); nine of the 28 dust varieties listed here fail our $\chi^2$ criteria for a ‘good’ fit.

Models that represent amorphous carbons and graphites tend to have higher $\kappa_{\text{abs}}$ values and thus return lower inferred dust masses, $M_{\text{dust}} \sim 0.4–0.8 M_{\odot}$ (Table 4). Amorphous silicates from Demyk et al. (2017), which have more than a factor of 10 larger $\kappa_{850}$ compared to those from Semenov et al. (2003), also yield a relatively moderate mass of $\sim 0.27 M_{\odot}$. 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 mm. As a consequence, the inferred dust masses for these grain types are $\sim 0.5–0.74 M_{\odot}$ 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 years, 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 $\sim 1–4 M_{\odot}$ to fit the observed SED. The CDE composite model of Dwek & Arendt (2015) results in

a moderate mass of $0.38 M_{\odot}$, 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, $SiO_2$, organics, and water ice.

The largest source of uncertainty in the observed dust mass is the choice of dust emission profile$^4$, specifically the value of the mass absorption coefficient ($\kappa_{\text{abs}}$) in the submm. We attempted to investigate spatial variations in the fitted dust parameters across the ALMA maps, however the limiting beam size translates to only 2–3 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 (Fig. 7). The inferred masses very clearly follow an inverse linear relation based on the submm $\kappa_{\text{abs}}$ value where $M_{\text{dust}}[M_{\odot}] = 0.117 \times \kappa_{850}^{-1}$. 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.

### 4.3. The Ring

The ring flux densities are shown in Fig. 5, along with previous measurements of the ring from ALMA and the Australia Telescope Compact Array (ATCA): 1.4 GHz, 18 GHz, 44 GHz (Zanardo et al. 2013); 9 GHz (Ng et al. 2013); 94 GHz (Lakićević et al. 2012). As the frequency increases towards the submm 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, & 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 $\sim 30\%$. The integrated ring emission follows a non-thermal power law profile of the form $S_{\nu} \propto \nu^{\alpha}$. The spectral index $\alpha$ was previously estimated from ATCA radio (Zanardo et al. 2014) and ALMA Cycle 0 data (Indebetouw et al. 2014) to be $-0.73 \pm 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 $\alpha = -0.70 \pm 0.06$.

The radio ring emission has steadily increased due to the synchrotron-producing electrons being accelerated by expanding shockwaves (Zanardo et al. 2010; Staveley-Smith et al. 2014). Recently, Cendes et al. (2018) have fit the

<footnote>$^4$ We did not attempt to fit a two-temperature component modBB as the SED shape is narrow, which suggests only one component is necessary (Matsuura et al. 2011, 2015; Mattsson et al. 2015).</footnote>
<figure><img src="image_6.png" /></figure>

**Figure 6.** Top: Modified blackbody fits to the Herschel and ALMA observations of the ejecta continuum in SN 1987A, using the 28 dust varieties and parameters listed in Table 4. The Herschel 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, TiO₂, Al₂O₃, NaAlSi₂O₆, and MgAl₂O₄) in light blue, carbide (SiC) in pink, sulfide (FeS) in green, and pure iron in light purple. Bottom: Dust $\kappa_{\text{abs}}(\lambda)$ curves used for the modified blackbody fits, showing the large variation in emission profiles.

radio emission of the (2D) ring and (3D) torus emission models across many epochs as a power law of the form $S(\nu, t)=K\nu^\alpha(t-t_0)^\beta$, where $\alpha$ is the spectral index of the emission across the spectrum, $\beta$ is the power law slope, and K 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

<figure><img src="image_7.png" /></figure>

**Figure 7.** Fitted dust mass as a function of $\kappa_{850}$ from various dust models. The colors for the different dust varieties are the same as in Fig. 6. The fitted dust mass closely follows a linear inverse trend over 3 orders of magnitude: $M_{\text{dust}} \propto \kappa_{850}^{-1}$.

10402 where $K=1.5\pm0.1$, $\alpha=0.70$ and $\beta=0.59\pm0.02$. We find no evidence of a contribution from dust in the ring to the mm 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 modelled the integrated SED of the dust emission under the assumption of uniform temperature and density within the ejecta. However, in §3, we saw that the SiO $J=5 \to 4$ and CO $J=2 \to 1$ images exhibit a hole (Fig. 3) where the dust emission peaks, and the SiO line ratio indicates lower SiO $J=5 \to 4$ brightness with respect to SiO $J=6 \to 5$ and SiO $J=7 \to 6$ (seen in the SiO line ratio maps in Fig. 8). 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 code RADEX (van der Tak et al. 2007).

### 5.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 the LAMBDA database (Schöier et al. 2005), and we use H₂-SiO collisional cross sections based on the calculations by Dayou & Balança (2006). In the ISM, H₂ is widely assumed to be the dominant collisional partner in molecular clouds; however, that
<figure><img src="image_8.png" /></figure>

**Figure 8.** Molecular line ratios of SiO. The line emission in each channel was converted to flux [W m⁻²] integrated over the 300 km s⁻¹ bin before division. Channel velocity centers are LSRK; for reference, the SN 1987A systematic velocity is 287 km s⁻¹. Left: SiO *J*=7→6 / SiO *J*=5→4. Note that while the 300 km s⁻¹ window covers the line centers, the SiO *J*=7→6 line has poor coverage there due to that frequency falling at the edges of the two observed tunings. Right: SiO *J*=6→5 / SiO *J*=5→4. Intensity maps were convolved to the SiO *J*=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.

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 synthesised elements, with hydrogen being depleted. In the two layers containing abundant Si, the major elements are O and S (e.g., Woosley 1988), and the collisional partner of SiO is likely to be O₂ 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 the RADEX input parameters is the FWHM of
the line width of the Gaussian (Δv). We adopted a Δv of
400 km s⁻¹, whose integrated area over the Gaussian profile
would be equivalent to that of a box-shaped 300 km s⁻¹ 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 anal-
ysis, 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 in RADEX calculations are
the kinetic temperature (T<sub>kin</sub>), the density of the collisional
partner (n<sub>coll</sub>) and the column density (N<sub>SiO</sub>). In the opti-

cally thin regime, as found in the calculated parameter range
where solutions are found, the SiO line intensities are deter-
mined by N<sub>SiO</sub> together with the area filling factor and the
expansion velocity v<sub>exp</sub>, while the SiO line ratios are de-
termined by T<sub>kin</sub> and n<sub>coll</sub>, independent of N<sub>SiO</sub>. The fill-
ing factor is defined as the area of the line emitting frac-
tion 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 fill-
ing factor of 1–50% is a reasonable range. In this analysis,
we adopted a column density grid of N<sub>SiO</sub>=10<sup>12</sup>–10<sup>18</sup> cm<sup>−2</sup>
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 are n<sub>coll</sub>=10<sup>3</sup>–10<sup>10</sup> cm<sup>−3</sup> and T<sub>kin</sub>=10–200 K.
Matsuura et al. (2017) suggested the temperature range is be-
low 190 K for SiO, with the CO kinetic temperature between
30–50 K, so we restricted the analysis to below 200 K. Al-
though we include temperatures up to 200 K in the RADEX
calculations, it is very unlikely that 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 SiO J=5→4 and SiO J=6→5 were convolved to match the lowest spatial resolution – that of the SiO J=7→6 beam – with a uniform pixel width of 0.′015. The convolved and regridded line ratio maps, made on a channel-by-channel basis, are displayed in Fig. 8. Intensities were averaged over 5×5 pixels in order to increase the signal to noise ratio in the SiO J=7→6 image. As the minor-axis FWHM of the beam is 0.′17 for SiO J=7→6, there is a small loss of spatial information by this averaging.

One caveat: the continuum subtraction for SiO J=6→5 was performed on the final imaged data cube, whereas for SiO J=5→4 and SiO J=7→6 the continuum subtraction was done in uv space before imaging (imcontsub vs. uvcontsub in CASA). Continuum subtraction in uv 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 SiO J=6→5 the two methods should give similar results. As discussed in § 3, there was no appreciable difference between the two methods for the CO J=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 SiO J=6→5 and SiO J=5→4, and 5% for SiO J=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 SiO J=6→5 were higher than the measured value, the line ratio of SiO J=6→5/SiO J=5→4 would increase, but the ratio of SiO J=7→6/SiO J=6→5 would decrease. The Cycle 2 flux densities being systematically lower than Cycle 0 as discussed in § 4.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 Fig. 9, in order to understand the change in physical parameters ($T_{\text{kin}}$, $n_{\text{coll}}$) within the ejecta, which in turn affect the SiO emission. One is close to the molecular hole seen in SiO J=5→4 and CO J=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 SiO J=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⁻¹ LSRK channel, which has good S/N for all three SiO lines, and the hole is clearly identified.

<figure><img src="image_9.png" /><figcaption>**Figure 9.** Locations of regions A and B used for the RADEX analysis of the SiO transitions. The background image is the SiO J=5→4 line at 0 km s⁻¹ (LSRK) channel, which clearly shows the molecular hole at region A. Region B is a representative region of the general SiO emitting ejecta.</figcaption></figure>

In region A, which has a hole in SiO J=5→4, we found that a column density of $N_{\text{SiO}} = 10^{16}$ cm⁻² is reasonable. For $N_{\text{SiO}} = 10^{15}$ cm⁻² or below, the predicted line intensities do not reach the measured levels, assuming a filling factor of 50%. Slightly larger values of $N_{\text{SiO}} = 10^{17}$ cm⁻² can match the measurements, but beyond $10^{18}$ cm⁻² the predicted and measured line intensities do not match.

Fig. 10 demonstrates the plausible ranges for the kinetic temperature ($T_{\text{kin}}$) and collisional partner’s density ($n_{\text{coll}}$) at region A, for $N_{\text{SiO}} = 10^{16}$ cm⁻². The collisional partner for the RADEX calculations is H₂; see the last paragraph of this section for discussion of this point. In Fig. 10 (a)–(c), the colored contours show the calculated ratios, the black lines show the measured SiO line ratios, and the $1\sigma$ background uncertainties and $1\sigma$ systematic uncertainties are indicated with dashed and dotted lines, respectively. Fig. 10 (d) compares the possible ranges of the kinetic temperature ($T_{\text{kin}}$) and collision partner density ($n_{\text{coll}}$) from the three SiO line ratios. The uncertainty of the SiO J=7→6/SiO J=5→4 ratio has the tightest constraint, so only the uncertainty of this ratio is plotted in Fig 10 (d). Fig 10 (d) shows that the ratios of all three SiO lines can be matched with a very similar set of $T_{\text{kin}}$ and $n_{\text{coll}} > 10^7$ cm⁻³, 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 is $T_{\text{kin}}=34$ K with LTE conditions at $n_{\text{coll}} > 10^7$ cm⁻³, the second possibility is a non-LTE condition with a range of ($10^6 < n_{\text{coll}} < 10^7$ cm⁻³) and $50 < T_{\text{kin}} < 200$ K, and finally the curves connecting these two conditions. The line center turns optically thick at the column density of $N_{\text{SiO}} = 10^{16}$ cm⁻², but the majority
<figure><img src="image_10.png" /><figcaption>**Figure 10.** (a)–(c) Results of RADEX modeling of the SiO line ratios at region A, representing the molecular “hole”. These calculations assume a column density of $N_{\text{SiO}} = 10^{16} \text{ cm}^{-2}$. (d) The 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 H$_2$ density values to the observed line ratios, while the dashed and dotted curves represent the 1 $\sigma$ background uncertainties and 1 $\sigma$ systematic uncertainties, respectively.</figcaption></figure>

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%.

By increasing the column density from $N_{\text{SiO}} = 10^{16} \text{ cm}^{-2}$ to $N_{\text{SiO}} = 10^{17} \text{ cm}^{-2}$, the SiO ratios change much more gradually as a function of the kinetic temperature ($T_{\text{kin}}$) and collision partner density ($n_{\text{coll}}$) (Fig. 11), expanding the feasible parameter space. This is the effect of higher optical depth. The sets of $T_{\text{kin}}$ and $n_{\text{coll}}$ from the three different calculated SiO ratios do not overlap in Fig. 11, but it is still possible to consider a wide range of $T_{\text{kin}}$ at $n_{\text{coll}} > 10^7 \text{ cm}^{-3}$,

from 19–22 K within the 1 $\sigma$ uncertainty. In the non-LTE range, the required $n_{\text{coll}} > 10^7 \text{ cm}^{-3}$ is more or less comparable to that for $N_{\text{SiO}} = 10^{16} \text{ 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 the $N_{\text{SiO}} = 10^{16} \text{ cm}^{-2}$ case.

At region B, the physical parameters are slightly different from those at the hole (region A). Fig. 12 shows the possible parameter space which fits the SiO ratios for region B; only $N_{\text{SiO}} = 10^{16} \text{ cm}^{-2}$ gives a feasible range. The most plausible temperature is $T_{\text{kin}}=18$ K for LTE conditions ($n_{\text{coll}} > 10^7 \text{ cm}^{-3}$), and an alternative possibility is non-LTE
<figure><img src="image_11.png" /><figcaption>**Figure 11.** The SiO line ratios that match the measured ALMA values at the SiO hole. The same as Fig. 10 (d) but for a column density $N_{\text{SiO}} = 10^{17} \text{ cm}^{-2}$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 12.** The same SiO calculations as in Fig. 10 (d) but for region B with column density $N_{\text{SiO}} = 10^{16} \text{ cm}^{-2}$.</figcaption></figure>

with $3 \times 10^5 < n_{\text{coll}} < 10^7 \text{ cm}^{-3}$ and $18 < T_{\text{kin}} < 200 \text{ 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.

In summary, the difference in the modelled line ratios and intensities near the SiO $J=5 \to 4$ and $J=6 \to 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 H$_2$, but rather other molecules such as O$_2$ or SiS, according to chemical models (Sarangi & Cherneff 2013). Therefore, instead of higher H$_2$ 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 O$_2$, whereas in the hole the dominant partner could be SiS. However, as will be discussed in § 6.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 CO $J=6 \to 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 beam-width) regions to aid our analysis. The top panel of Fig. 13 shows the location of 20 regions selected across the CO-emitting ejecta, overlaid on top of the CO $J=2 \to 1$ emission. Nine regions, each the size of the CO $J=6 \to 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 Fig. 13 compares the summed spectra of CO $J=6 \to 5$ and CO $J=2 \to 1$ (with the latter convolved to the CO $J=6 \to 5$ beam before integration) for these 9 regions showing their location with respect to the CO $J=2 \to 1$ emission. The spectra are in units of mJy, having been spatially integrated over each region.

The bottom panel of Fig. 13 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 CO $J=2 \to 1$ and CO $J=6 \to 5$ transitions are similar to each other in the scaled spectra (see for example regions 6, 7, and 9 in Fig. 13). However, in regions 1 and 8, the CO $J=6 \to 5$ profile is suppressed at negative velocities with respect to the CO $J=2 \to 1$ line. The CO $J=6 \to 5$ emission is also suppressed with respect to CO $J=2 \to 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 CO $J=6 \to 5$ emission at velocities of –1000–+1000 km s$^{-1}$ with respect to CO $J=2 \to 1$ (bottom panel Fig. 13); 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 rep-
<figure><img src="image_13.png" /><figcaption>**Figure 13.** Top: Regions chosen for comparing the CO *J*=2→1 and CO *J*=6→5 fluxes and line profiles; each region spans one CO *J*=6→5 beam. Regions of particular interest are labeled 1–9. The shaded greyscale contours illustrate the CO *J*=2→1 emission (lower left panel of Fig. 3). Middle: A qualitative comparison of the stacked CO *J*=2→1 (blue) and CO *J*=6→5 (red) spectra in each region with a representative scale shown in the bottom left corner. Bottom: A 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⁻¹ (LSRK).</figcaption></figure>

representative of the general SiO emitting properties across the molecular ejecta based on the SiO line ratios, whereas Fig. 13 demonstrates that region 8 has very different CO gas properties to most of the other regions).

Using this information, we carry out an analysis with RADEX similar to that for SiO. Fig. 14 shows the resulting RADEX calculations of the CO *J*=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 to the W m⁻² 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<fnref n="5" />.

As demonstrated originally in Fig. 3 in Kamenetzky et al. (2013), the CO *J*=6→5 line is sensitive to temperature change. We propose here that the CO *J*=6→5 suppression with respect to CO *J*=2→1 indicates the gas is at a lower temperature in regions 1 and 8 (where we also see a peak in the dust emission, § 3, Fig. 4) compared to the surrounding regions. In these regions, the blue wing of the CO *J*=6→5 emission is lower, thus if the dust and CO *J*=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 CO *J*=6→5 line, we cannot specify the exact excitation temperatures in these regions. However, the CO excitation temperature is higher near the CO *J*=2→1 hole (regions 4 and 5) at velocities of −1000–1000 km s⁻¹.

# 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 (Table 4). Inferred dust masses of several tens of solar masses, as in the case of TiO₂ (Posch et al. 2003) which requires 82 *M*☉ of dust, or MgAl₂O₄ Fabian et al. (2001) which requires 122 *M*☉ for example, are clearly untenable as they are larger than the progenitor star mass (18–20 *M*☉, Woosley 1988). This rules out a further three dust varieties in Table 4 (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

<footnote n="5">⁵ The spectra and channel maps in Figs. 13 and A.2 are spectral density units (mJy and mJy per beam). In order to compare the line ratios in the integrated fluxes in W m⁻², the flux densities in mJy per velocity channel units need to be multiplied by a factor of ∂f/∂v = −f₀/c to account for the change from integrating in velocity v to frequency f. For the CO 6→5/2→1 ratio, multiply by f₀CO65/f₀CO21 ~ 3.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 14.** Ratios of CO $J=6\to 5$ and CO $J=2\to 1$ flux densities in Jy units derived using RADEX 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 Fig. 13.</figcaption></figure>

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 M_{\odot}$ $Z = 0.008 Z_{\odot}$ 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 M_{\odot}$, with $1.21 M_{\odot}$ of oxygen. Focusing on individual dust varieties, the limits can be further differentiated. The Nomoto et al. (2013) $18 M_{\odot}$ model predicts $0.149 M_{\odot}$ 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 gives both 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 M_{\odot}$. For silicate dust, good fits to the SED produce masses of $0.3-0.7 M_{\odot}$, 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 M_{\odot}$.

The iron yield from the Nomoto et al. (2013) $18 M_{\odot}$ nucleosynthesis model is $0.079 M_{\odot}$ ($^{56}$Fe only) or $0.085 M_{\odot}$ (including all isotopes); this is orders of magnitude less than the inferred pure iron dust model (Henning & Stognienko 1996), which requires $3.97 M_{\odot}$ of dust to fit the SED. The Woosley & Weaver (1995) $15 M_{\odot}$ and $18 M_{\odot}$ $Z = 0.1 Z_{\odot}$ progeni-

tor models predict iron masses ranging from $0.14-0.20 M_{\odot}$, where roughly half of the iron originates from $^{56}$Ni 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 & Stognienko 1996; Henning et al. 1995) or resulting in unrealistically high dust masses for the pure iron model.

Limits from the Nomoto et al. (2013) $18 M_{\odot}$ explosive synthesis model for some other common dust varieties include $0.337 M_{\odot}$ for forsterite (Mg$_2$SiO$_4$), $0.481 M_{\odot}$ for enstatite (MgSiO$_3$), $0.014 M_{\odot}$ for alumina (Al$_2$O$_3$), and $0.373 M_{\odot}$ for silica (SiO$_2$) – assuming 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 M_{\odot}$ for forsterite, $4.1 M_{\odot}$ for enstatite and $0.9 M_{\odot}$ for alumina (Table 4). 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 M_{\odot}$ given the available metal budget. From this it is possible to immediately rule out a further seven dust varieties listed in Table 4 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+MgSiO$_3$+FeS would give a total metal mass, and an upper limit on the total dust mass, of $0.72 M_{\odot}$. However, taking the yields of the Woosley (1988) $18 M_{\odot}$ $0.1 Z_{\odot}$ progenitor model, for example, would give $0.55 M_{\odot}$ 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 supernovae 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 M_{\odot}$.

## 6.2. The Spatial Distributions of Dust, CO and SiO, and Chemistry Leading to Dust Formation

Our spatially resolved images (Fig. 3) show that the dust distribution is clumpy and asymmetric. Comparing the spatial distribution of the dust with the CO $J=6\to 5$ reveals a
weak anti-correlation with the integrated CO J=6→5 and dust distributions, while there is little spatial correlation between the dust and SiO images. We suggested earlier that this anti-correlation occurs because CO J=6→5 is suppressed compared with CO J=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 SiO RADEX analysis). Here, the hole in SiO and CO J=2→1 coincides with the dust peak, with relatively strong CO J=6→5 emission observed at −1000–1000 km s⁻¹ with respect to CO J=2→1, and this strong CO J=6→5 continuing to its neighboring region (region 5 in Fig. 13). We will discuss this region separately in § 6.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 & Cherchneff 2013).
The outermost region is the H-envelope, followed by the He
shell, which can contain more carbon atoms than oxygen
atoms (Woosley & Weaver 1995; Rauscher et al. 2002). The
He shell can also form CO, as it enriches with C and O (e.g.,
Sarangi & Cherchneff 2013). 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 a ⁵⁶Ni
core, that also contains Si, but very low C or O. Here we
interpret an apparent anti-correlation between the dust and
CO J=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 supernova dust formed in chemical models. Sarangi & Cherchneff (2015) predicts 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) and Sarangi & Cherchneff (2015) included more extensive chemical networks, as well as dissociation of molecules by energetic electrons, and found this resulted in very few carbonaceous dust grains (6 × 10⁻³ M☉ of carbonaceous dust out of total of 0.04 M☉ dust mass for a 19 M☉ star). One simple explanation is that Sarangi & Cherchneff (2013) and Sarangi & Cherchneff (2015) modelled 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, sug-

gesting that any dissociation process must be more efficient
on longer time scales. Alternatively, the chemical reaction
rate used in the models might have large uncertainties, par-
ticularly involving the highly energetic electrons. Together
with the recent detection of HCO⁺ (Matsuura et al. 2017),
which was largely under-predicted in the chemical model of
Sarangi & Cherchneff (2013), this suggests 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 al-
low 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,
private communication).

An alternative explanation for the anti-correlation 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 anti-correlation 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 due
to radiative cooling via dust emission, i.e., the cooler gas tem-
perature 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 conse-
quence 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 Fig. 3,
for example, in comparison with the 315 GHz image, but the
surface brightness is lower in those locations. The CO-dust
anti-correlation 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 & Fesen 2015) 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 den-
sity 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 CO $J=2 \to 1$ and SiO $J=5 \to 4$ line-emitting ejecta, and we also see emission located in the hole in the SiO $J=7 \to 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
(Gabler et al. in preparation) compared with ALMA obser-
vations (Fig. 4 in Abellán et al. 2017) show, that in the ex-
plosions of three out of the four simulated progenitor models,
the final density distributions of Si and C (each multiplied
with 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 com-
presses the material such that the entropy increases signifi-
cantly and an outwards moving shock forms. This feature is
termed the ‘self-reflected’ shock, first discussed in Ertl et al.
(2016) (with forthcoming details in Gabler et al. in prepa-
ration). The self-reflected shock accelerates the innermost
ejecta compared to homologous expansion and leaves a re-
gion with lower density in the center. When passing more
and more of the ejecta, the shock loses energy. It finally dis-
sipates and cannot accelerate the outermost ejecta. This ac-
celeration of the low-velocity ejecta leads to a formation of
a higher density shell-like configuration as observed in the
emission of the transition lines SiO $J=5 \to 4$ or CO $J=2 \to 1$
(Fig. 4 in Abellán et al. 2017). However, in the hydrody-
namic simulations, one model based on the progenitor model
from Shigeyama & Nomoto (1990) (model N20 in Wong-
wathanarat 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 in-
creased brightness for the observed properties of the blob,
we suggest two possibilities: (i) gas heated by a compact ob-
ject or (ii) a clump heated by $^{44}$Ti 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 radia-
tion from the compact object, with an early development of a
pulsar wind nebula.

$^{44}$Ti 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 that $^{44}$Ti is located more or less at a similar radial extent as $^{56}$Ni and $^{28}$Si (Woosley & Weaver 1995; Wongwathanarat et al. 2015, 2017) with almost identical bulk velocities, though the modeled distribution of $^{44}$Ti could be more uncertain than the other elements due to its greater sensitivity to the explosion physics (Jerkstrand et al., in preparation). Observations show that the $^{44}$Ti ejecta is redshifted, suggesting that the bulk is moving away from the observer (Boggs et al. 2015). Since the predicted distribution of $^{44}$Ti shows qualitatively similar properties as Fig. 4 in Abellán et al. (2017), because it is subject to effectively the same hydrodynamical history (Jerkstrand et al., in preparation), it is unlikely that gas heated by $^{44}$Ti decay would be more centrally distributed and also be co-located 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
& Barlow 2015), 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. 2018). In the Crab Nebula, the radiation from the pulsar
wind can partially dissociate gas, contributing to the forma-
tion of the molecules OH+ and $^{36}$ArH+ (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, & Kushnir 2016), 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 & Fransson 2016; Alp et al. 2018b). The absorbed energy is re-processed 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 redshifted $^{44}$Ti spectrum (Boggs et al. 2015), the direction of the compact object is predicted to be moving toward us and extending along the north-east 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): $\alpha = 05^{h}35^{m}27^{m}.9875$, $\delta = -69^{\circ}16'11''.107$ (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, towards 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).

We estimate the dust peak flux density as $S_{679\text{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 \sim 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. Figure 15 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 & Blandford 2014), 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 Fig. 15. 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 gives $S_{\nu}$ 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 $L_{\text{bol,dust}}=40–90$ L$_{\odot}$ 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 K to 33 K. However, the compact source is most likely additional heating, on top of $^{44}$Ti-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$_{\odot}$.

In 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 of $^{44}$Ti; (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 towards 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 non-thermally 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\alpha$ hole, whether from a physical lack of material (due to the reverse shock or a simple void due to the expansion of ejecta), or 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., McCray 2003; Larsson et al. 2011; Fransson et al. 2013; Larsson et al. 2016, 2019). A simple estimate of the extinction can be made by assuming the dust fills a spherical shell with uniform density, and calculating the optical depth using the mass extinction coefficient $\kappa_{\text{ext}}$ (related to, but different from the $\kappa_{\text{abs}}$ used in the modBB fits) according to $\tau_{\lambda} = \int \kappa_{\text{ext},\lambda} \rho \, ds$. For shell edges corresponding to the features in the images (0.′05–0.′′2, or ~0.4–1.5×10$^{15}$m), the Zubko et al. (1996) ACAR model and its corresponding dust mass fit give an opacity at 6563 Å of $\tau_{H\alpha} \sim 560$. Silicate dust varieties tend to have lower $\kappa_{\text{ext},\lambda}$ at a given wavelength – the amorphous forsterite model of Jäger et al. (2003), for example, gives $\tau_{H\alpha} \sim 400$.
<figure><img src="image_15.png" /></figure>

**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 it corresponds to a temperature of 33K. 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-grey triangles). For comparison, the spectra of the Crab Nebula and the Crab Pulsar (purple and red curves, respectively; Bühler & Blandford 2014, and references therein) are shown scaled to the distance of the LMC.

$\kappa_{ext,\lambda}$ 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 of $a=5$ µm, the optical depth of ACAR dust at $H\alpha$ is 10. Even this level is optically thick, meaning that dust extinction can play a non-negligible role in the observed $H\alpha$ distribution.

## 7. SUMMARY

We have observed SN 1987A with ALMA in Cycle 2, 10352–10441 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.

* We find that the dust emission in the ejecta is clumpy and asymmetric, fitting within the $H\alpha$ “key hole”, 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 $\tau \sim 500$. Dust extinction could be a factor in the appearance of the $H\alpha$ hole.

* We see an anti-correlation between the CO $J=6 \to 5$ emission and dust, and this anti-correlation is not seen when comparing dust and CO $J=2 \to 1$. Our RADEX analysis suggests 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 due to dissociation of CO, contrary to chemical predictions.

* 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 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.

* We fit the spectral energy distribution of the dust emission from the ejecta with modified blackbody profiles. The derived dust masses and temperatures depend on the submm emissivity of the dust, which is not very well determined. Temperatures from the fits are generally between 18–23K. Amorphous carbon and graphite models have dust with high emissivity and so yield the lowest dust masses, around $0.4–1.6 M_{\odot}$. Silicates return higher dust mass estimates in the range $0.6–4 M_{\odot}$. “Typical ISM” dust varieties from the Milky Way and nearby galaxies give masses between $1–2 M_{\odot}$. Taking the total mass of available metals excluding oxygen predicted to be ejected by a 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 M_{\odot}$ for carbon or Si grains, or a total of $< 0.7 M_{\odot}$ 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.

*Acknowledgments:* We acknowledge Arka Sarangi, Loretta Dunne, Steve Maddox, and the anonymous reviewer for interesting and useful discussions. PJC and HLG acknowledge support from the European Research Council (ERC) in the form of Consolidator Grant COSMICDUST (ERC-2014-CoG-647939, PI HLGomez). MM 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 Mat-
ter in Astro- and Particle Physics" (NDM) and the Ex-
cellence Cluster "ORIGINS: From the Origin of the Uni-
verse to the First Building Blocks of Life" (EXC 2094;
http://www.origins-cluster.de). MJB acknowledges support
from ERC grant SNDUST (ERC-2015-AdG-694520). JCW
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 co-
operation 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 & Granger 2007), lmfit (Newville et al. 2016), matplotlib (Hunter 2007), multicolorfits (Cigan 2019), numpy (van der Walt et al. 2011), RADEX (van der Tak et al. 2007), SAOImage DS9 (Joye & Mandel 2003), scipy (Jones et al. 2001–), WWCC (Arshakian & Ossenkopf 2016).

REFERENCES

Abellán, F. J., Indebetouw, R., Marcaide, J. M., et al. 2017, ApJL, 842, L24

Ahmad, I., Greene, J. P., Moore, E. F., et al. 2006, PhRvC, 74, 065803

Alp, D., Larsson, J., Fransson, C., et al. 2018, ApJ, 864, 174

Alp, D., Larsson, J., Fransson, C., et al. 2018, ApJ, 864, 175

Arshakian, T. G., & Ossenkopf, V. 2016, A&A, 585, A98

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Barlow, M. J., Swinyard, B. M., Owen, P. J., et al. 2013, Science, 342, 1343

Begemann, B., Dorschner, J., Henning, T., et al. 1997, ApJ, 476, 199

Blum, K., & Kushnir, D. 2016, ApJ, 828, 31

Boggs, S. E., Harrison, F. A., Miyasaka, H., et al. 2015, Science, 348, 670

Bouchet, P., Danziger, I. J., & Lucy, L. B. 1991, AJ, 102, 1135

Bouchet, P., Dwek, E., Danziger, J., et al. 2006, ApJ, 650, 212

Brown, G. E., Bruenn, S. W., & Wheeler, J. C. 1992, Comments on Astrophysics, 16, 153

Bühler, R., & Blandford, R. 2014, Reports on Progress in Physics, 77, 066901

Burrows, A. 1988, ApJ, 334, 891

Burrows, D. N., Michael, E., Hwang, U., et al. 2000, ApJL, 543, L149

Callingham, J. R., Gaensler, B. M., Zanardo, G., et al. 2016, MNRAS, 462, 290

Cendes, Y., Gaensler, B. M., Ng, C.-Y., et al. 2018, ApJ, 867, 65

Chawner, H., Marsh, K., Matsuura, M., et al. 2018, MNRAS, Cherchneff, I., & Dwek, E. 2009, ApJ, 703, 642

Cherchneff, I., & Dwek, E. 2010, ApJ, 713, 1

Cigan, P. 2019, MultiColorFits: Colorize and combine multiple fits images for visually aesthetic scientific plots, ascl:1909.002

Clayton, D. D. 2011, NewAR, 55, 155

Danziger, I. J., Gouiffes, C., Bouchet, P., et al. 1989, IAUC4746, 1

Dayou, F., & Balança, C. 2006, A&A, 459, 297

De Looze, I., Barlow, M. J., Swinyard, B. M., et al. 2017, MNRAS, 465, 3309

Demyk, K., Meny, C., Leroux, H., et al. 2017, A&A, 606, A50

Deneault, E. A.-N., Clayton, D. D., & Meyer, B. S. 2006, ApJ, 638, 234

Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89

Dwek, E., Arendt, R. G., Bouchet, P., et al. 2010, ApJ, 722, 425

Dwek, E., Staguhn, J., Arendt, R. G., et al. 2014, ApJL, 788, L30

Dwek, E., & Arendt, R. G. 2015, ApJ, 810, 75

Ertl, T., Janka, H.-T., Woosley, S. E., et al. 2016, ApJ, 818, 124

Esposito, P., Rea, N., Lazzati, D., et al. 2018, ApJ, 857, 58

Fabian, D., Posch, T., Mutschke, H., et al. 2001, A&A, 373, 1125

Foreman-Mackey, D., Hogg, D. W., Lang, D., et al. 2013, PASP, 125, 306

France, K., McCray, R., Heng, K., et al. 2010, Science, 329, 1624

Frank, K. A., Zhekov, S. A., Park, S., et al. 2016, ApJ, 829, 40

Fransson, C., Larsson, J., Spyromilio, J., et al. 2013, ApJ, 768, 88

Fransson, C., Larsson, J., Migotto, K., et al. 2015, ApJL, 806, L19

Frick, P., Beck, R., Berkhuijsen, E. M., et al. 2001, MNRAS, 327, 1145

Gall, C., Hjorth, J., Andersen, A. C., 2011, ARA&A, 19, 43

Gall, C., Hjorth, J., Watson, D., et al. 2014, Nature, 511, 326

Gilmozzi, R., Cassatella, A., Clavel, J., et al. 1987, Nature, 328, 318

Goldsmith, P. F., & Langer, W. D. 1999, ApJ, 517, 209

Gomez, H. L., Clark, C. J. R., Nozawa, T., et al. 2012, MNRAS, 420, 3557

Gomez, H. L., Krause, O., Barlow, M. J., et al. 2012, ApJ, 760, 96

Gomez, H. 2013, Proceedings of the Life Cycle of Dust in the Universe: Observations, 146
Gröningsson, P., Fransson, C., Leibundgut, B., et al. 2008, A&A,
492, 481

Hashimoto M.-A., Nomoto K., Shigeyama T., 1989, A&A, 210, L5

Helder, E. A., Broos, P. S., Dewey, D., et al. 2013, ApJ, 764, 11

Henning, T., Begemann, B., Mutschke, H., et al. 1995, A&AS,
112, 143

Henning, T., & Stognienko, R. 1996, A&A, 311, 291

Henning, T., & Mutschke, H. 1997, A&A, 327, 743

Henning, T., Il'In, V. B., Krivova, N. A., et al. 1999, A&AS, 136,
405

Hirata, K., Kajita, T., Koshiba, M., et al. 1987, PhRvL, 58, 1490

Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90

Indebetouw, R., Matsuura, M., Dwek, E., et al. 2014, ApJL, 782, L2

Jaeger, C., Mutschke, H., Begemann, B., et al. 1994, A&A, 292, 641

Jager, C., Mutschke, H., & Henning, T. 1998, A&A, 332, 291

Jäger, C., Dorschner, J., Mutschke, H., et al. 2003, A&A, 408, 193

James, A., Dunne, L., Eales, S., et al. 2002, MNRAS, 335, 753

Janka, H-Th., Gabler, M., Wongwathanarat, A., 2017, IAUS 331, 148

Jerkstrand, A., Fransson, C., & Kozma, C. 2011, A&A, 530, A45

Jones, E., Oliphant, T., Peterson, P., et al. 2001–, SciPy: Open
source scientific tools for Python,
http://www.scipy.org/

Joye, W. A., & Mandel, E. 2003, Astronomical Data Analysis
Software and Systems XII, 295, 489

Kamenetzky, J., McCray, R., Indebetouw, R., et al. 2013, ApJL,
773, L34

Katsuda, S., Morii, M., Janka, H.Th. et al., 2018, ApJ 856, 18

Kirshner, R. P., Sonneborn, G., Crenshaw, D. M., et al. 1987, ApJ,
320, 602

Lakićević, M., Zanardo, G., van Loon, J. T., et al. 2012, A&A, 541, L2

Lakićević, M., van Loon, J. T., Meixner, M., et al. 2015, ApJ, 799, 50

Laor, A., & Draine, B. T. 1993, ApJ, 402, 441

Larsson, J., Fransson, C., Östlin, G., et al. 2011, Nature, 474, 484

Larsson, J., Fransson, C., Spyromilio, J., et al. 2016, ApJ, 833, 147

Larsson, J., Spyromilio, J., Fransson, C., et al. 2019, ApJ, 873, 15

Li, A., & Draine, B. T. 2001, ApJ, 554, 778

Lucy, L. B., Danziger, I. J., Gouiffes, C., et al. 1989, IAU Colloq.
120: Structure and Dynamics of the Interstellar Medium, 164

Manchester, R. N. 2007, Supernova 1987A: 20 Years After:
Supernovae and Gamma-ray Bursters, 134

Matsuura, M., Barlow, M. J., Zijlstra, A. A., et al. 2009, MNRAS,
396, 918

Matsuura, M., Dwek, E., Meixner, M., et al. 2011, Science, 333, 1258

Matsuura, M., Dwek, E., Barlow, M. J., et al. 2015, ApJ, 800, 50

Matsuura, M., Indebetouw, R., Woosley, S., et al. 2017, MNRAS,
469, 3347

Matsuura, M., De Buizer, J. M., Arendt, R. G., et al. 2019,
MNRAS, 482, 1715

Mattsson, L., Gomez, H. L., Andersen, A. C., et al. 2015, MNRAS,
449, 4079

McCray, R. 1993, ARA&A, 31, 175

McCray, R. 2003, Supernovae and Gamma-ray Bursters, 219

McCray, R., & Fransson, C. 2016, ARA&A, 54, 19

McMullin, J. P., Waters, B., Schiebel, D., et al. 2007, Astronomical
Data Analysis Software and Systems XVI, 127

Meixner, M., Panuzzo, P., Roman-Duval, J., et al. 2013, AJ, 146, 62

Michałowski, M. J., Murphy, E. J., Hjorth, J., et al. 2015, A&A,
522, A15

Milisavljevic, D., & Fesen, R. A. 2015, Science, 347, 526

Morgan, H. L., Edmunds, M.G., 2003, MNRAS, 343, 427

Morgan, H. L., Dunne, L., Eales, S. A., et al. 2003, ApJ, 597, L33

Mutschke, H., Begemann, B., Dorschner, J., et al. 1998, A&A, 333, 188

Newville, M., Stensitzki, T., Allen, D. B., et al. 2016, Lmfit:
Non-Linear Least-Square Minimization and Curve-Fitting for
Python, ascl:1606.014

Ng, C.-Y., Zanardo, G., Potter, T. M., et al. 2013, ApJ, 777, 131

Nomoto, K., Kobayashi, C., & Tominaga, N. 2013, ARA&A, 51, 457

Ormel, C. W., Min, M., Tielens, A. G. G. M., et al. 2011, A&A,
532, A43

Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943

Osterbrock, D. E. 1989, Astrophysics of Gaseous Nebulae and
Active Galactic Nuclei

Owen, P. J., & Barlow, M. J. 2015, ApJ, 801, 141

Panagia, N. 1999, New Views of the Magellanic Clouds, 190, 549

Pegourie, B. 1988, A&A, 194, 335

Pérez, F. & Granger, B. E. 2007, Computing in Science &
Engineering, 9, 21

Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A,
518, L1

Posch, T., Kerschbaum, F., Fabian, D., et al. 2003, ApJS, 149, 437

Potter, T. M., Staveley-Smith, L., Ng, C.-Y., et al. 2009, ApJ, 705, 261

Rauscher, T., Heger, A., Hoffman, R. D., et al. 2002, ApJ, 576, 323

Rho, J., Gomez, H. L., Boogert, A., et al. 2018, MNRAS, 479, 5101

Roche, P. F., Aitken, D. K., & Smith, C. H. 1991, MNRAS, 252, 39P

Roche, P. F., Aitken, D. K., & Smith, C. H. 1993, MNRAS, 261, 522

Rouleau, F., & Martin, P. G. 1991, ApJ, 377, 526

Rowlands, K., Gomez, H. L., Dunne, L., et al. 2014, MNRAS, 441, 1040
Russell, D. M., Fender, R. P., Hynes, R. I., et al. 2006, MNRAS, 371, 1334

Sarangi, A., & Cherchneff, I. 2013, ApJ, 776, 107

Sarangi, A., & Cherchneff, I. 2015, A&A, 575, A95

Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., et al. 2005, A&A, 432, 369

Semenov, D., Henning, T., Helling, C., et al. 2003, A&A, 410, 611

Shigeyama, T., & Nomoto, K. 1990, ApJ, 360, 242

Smartt, S. J., Eldridge, J. J., Crockett, R. M., et al. 2009, MNRAS, 395, 1409

Sonneborn, G., Pun, C. S. J., Kimble, R. A. et al., 1998, ApJ 492, L139

Spyromilio, J., Meikle, W. P. S., Learner, R. C. M., et al. 1988, Nature, 334, 327

Staveley-Smith, L., Potter, T. M., Zanardo, G., et al. 2014, Supernova Environmental Impacts, 15

Sukhbold, T., Ertl, T., Woosley, S. E., et al. 2016, ApJ, 821, 38

Temim, T., Gehrz, R. D., Woodward, C. E., et al. 2006, AJ, 132, 1610

Temim, T., Dwek, E., Tchernyshyov, K., et al. 2015, ApJ, 799, 158

Temim, T., Dwek, E., Arendt, R. G., et al. 2017, ApJ, 836, 129

Utrobin, V. P., Wongwathanarat, A., Janka, H.-T., et al. 2015, A&A, 581, A40

Utrobin, V. P., Wongwathanarat, A., Janka, H.-T., et al. 2019, A&A, 624, A116

van der Tak, F. F. S., Black, J. H., Schöier, F. L., et al. 2007, A&A, 468, 627

Van der Walt, S., Colbert, S., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22

Watson, D., Christensen, L., Knudsen, K. K., et al., 2015, Nature, 519, 327

Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296

Weisskopf, M. C., Hester, J. J., Tennant, A. F., et al. 2000, ApJL, 536, L81

Wesson, R., Barlow, M. J., Matsuura, M., et al. 2015, MNRAS, 446, 2089

West, R. M., Lauberts, A., Jorgensen, H. E., et al. 1987, A&A, 177, L1

White, G. L., & Malin, D. F. 1987, Nature, 327, 36

Wongwathanarat, A., Müller, E., & Janka, H.-T. 2015, A&A, 577, A48

Wongwathanarat, A., Janka, H.-T., Müller, E., et al. 2017, ApJ, 842, 13

Wooden, D. H., Rank, D. M., Bregman, J. D., et al. 1993, ApJS, 88, 477

Woosley, S. E. 1988, ApJ, 330, 218

Woosley, S. E., & Weaver, T. A. 1995, ApJS, 101, 181

Woosley, S. E., Heger, A., Weaver, T. A., et al. 1997, arXiv e-prints, astro-ph/9705146

Zanardo, G., Staveley-Smith, L., Ball, L., et al. 2010, ApJ, 710, 1515

Zanardo, G., Staveley-Smith, L., Ng, C.-Y., et al. 2013, ApJ, 767, 98

Zanardo, G., Staveley-Smith, L., Indebetouw, R., et al. 2014, ApJ, 796, 82

Zanardo, G., Staveley-Smith, L., Gaensler, B. M., et al. 2018, ApJL, 861, L9

Zhang, S.-B., Dai, S., Hobbs, G., et al. 2018, MNRAS, 479, 1836

Zubko, V. G., Mennella, V., Colangeli, L., et al. 1996, MNRAS, 282, 1321

Zubko, V., Dwek, E., & Arendt, R. G. 2004, ApJS, 152, 211
# APPENDIX

## A. MOLECULAR LINE CHANNEL MAPS

<figure><img src="image_16.png" /></figure>

**Figure A.1.** CO channel maps in 300 km s⁻¹ steps, LSRK; for reference, the SN 1987A systematic velocity is 287 km s⁻¹. Left: CO J=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: CO J=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 CO J=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.
<figure><img src="image_17.png" /></figure>

**Figure A.2.** SiO channel maps in 300 km s⁻¹ steps, LSRK; for reference, the SN 1987A systematic velocity is 287 km s⁻¹. Top Left: Continuum + SiO *J*=5→4 as presented by Abellán et al. (2017), also exhibiting the central molecular hole seen in CO *J*=2→1. Top Right: Continuum + SiO *J*=6→5 from Abellán et al. (2017). Bottom: SiO *J*=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 SiO *J*=7→6 at the spatial position of the molecular hole but slightly in front (at 0 km s⁻¹). 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.
## 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 hotspots in the ring to determine the center of SN 1987A. This was done by fitting 2D Gaussians to the hotspots in the HST/ACS R-band image from 2006 and then fitting 1D Gaussians along the same directions in all the other 32 R and B images from 2003 to 2016 (Table 8 of Alp et al. (2018a)). They derive a ring center position of $\alpha=5^{h}35^{m}27^{s}9875$, $\delta=-69^{\circ}16'11''107$ (ICRFJ2015.0), with uncertainty from bootstrapping the hotspot 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 (Fig. B.1) (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 latter 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 one pixel (Figure B.1). The final weighted fit gives a ring emission center of $\alpha=5^{h}35^{m}27^{s}998$, $\delta=-69^{\circ}16'11''107$ (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 (§ 2.1), the total uncertainty on the inferred center position is 18 mas, slightly larger than one 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><img src="image_18.png" /></figure>

**Figure B.1.** 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 to HST image ring hotspots in Alp et al. (2018a) are also shown.
C. SIMULATING POTENTIAL FLUX LOSS

<figure><img src="image_19.png" /><figcaption>**Figure C.1.** CASA simulations of ALMA observations at 315 GHz to test if diffuse flux is lost in the observing process. Modelled emission that resembles the true image is recovered qualitatively and quantitatively and suggests we are not missing a significant component of extended diffuse emission.</figcaption></figure>

Simulations exploring the possibility of over-resolving extended dust emission were performed with simobserve and simanalyze in CASA. 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 (2′′4 across). This is similar to the largest angular scales these observations were sensitive to (see Table 1). 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 Fig. C.1: 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.

Notably, 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.