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

Introduction

Multiwavelength studies of Supernova 1987A (SN 1987A), located at a distance of $51.4 \pm1.2$ kpc in the Large Magellanic Cloud [102], have provided unprecedented details of how supernova (SN) explosions trigger the dynamical distribution of gas in a supernova remnant (SNR), and how this SN/SNR system evolves over time. The morphology of SN 1987A is well studied (see the recent review in McCray & Fransson [87]), with the system consisting of ejecta, and a bright and distinct equatorial ring (hereafter the ring), together with two fainter outer rings. The ring is composed of circumstellar material that radiates in UV, optical, X-rays, and radio over an extent of $1.\!^{\prime\prime}06$ ($0.3$ pc) (e.g. [16, 121] (Ng et al. 2013)), as well as thermal dust emission due to shock heating of pre-existing dust formed during the red supergiant phase ([12, 30]). The ejecta have a complex morphology. The H$\alpha$ emission, originating from warm gas irradiated by X-rays from the ring ([73, 39]), exhibits an elongated north-south structure and a ‘hole’ in the center. Along with hydrogen lines from the ejecta, near-infrared (NIR) emission from warm ($\sim2000$ K) CO and mid-infrared (MIR) emission from SiO in the SN ejecta were detected early (as early as 112 days) after the explosion (e.g., [122, 110]). After day 9,000, cold expanding CO, SiO and HCO$^{+}$ molecules were detected in the submillimeter (submm) part of the spectrum ([67, 82]), highlighting that a significant part of the ejecta is cold ($13$–$132$ K). Interestingly, the inner ejecta of SN 1987A have not yet mixed with the circumstellar medium (CSM) or interstellar medium (ISM) and the majority has not yet passed through the reverse shock ([37, 38]). Thus this young SNR is an ideal source for studying the footprints of the gas dynamics since the very early days of the SN, as the gas has been assumed to be free-expanding since its explosion ([85]). Indeed, recent high angular resolution emission line images of SiO and CO ([1]) from the Atacama Large Millimeter/submillimeter Array (ALMA) have been used to compare the distribution of the molecular gas ejecta with the predictions from models of the gas dynamics after the SN ([138], Gabler et al., in preparation).

The progenitor of SN 1987A, Sanduleak $-69^{\circ}$ 202, was a blue supergiant ([136, 137, 44, 69]), thought to have had a zero-age main sequence mass of $\sim19\,M_{\odot}$ ([143] (Hashimoto, Nomoto & Shigeyama 1989)), with a mass of $\sim14\,M_{\odot}$ at the time of the explosion ([141, 120, 124]). From its mass, the expectation is that a neutron star should have formed at the time of explosion. Despite prompt neutrino emission observed at the burst ([56]) indicating the formation of a neutron star ([15, 124]), the search for a compact object associated with SN 1987A has been difficult: observational searches have proven unfruitful (e.g. [78, 3, 34, 148]). The possible detection of radio polarization towards the ejecta ([147]) hints at the presence of magnetized shocks, potentially due to a compact object. Alp et al. [3] proposed that a thermally-emitting neutron star could be dust-obscured, and that this may be detectable as a point source in far-infrared (FIR) or submm images of the remnant, though this has not yet been detected.

It is still largely debated whether or not SNe are net dust producers or destroyers in galaxies (e.g., [92, 79, 42, 48, 31, 113, 132, 71, 126] (Michalowski et al. 2015)). Due to its youth and proximity, SN 1987A is an excellent laboratory for studying SN dust. It is also rare, since any dust emission seen in the inner region of the remnant can be attributed unambiguously to dust formed in the supernova ejecta and not from the swept-up CSM/ISM or unrelated foreground/background material (a common issue with Galactic SNRs, e.g., [93, 46, 26, 19]). SN 1987A also provides insight into dust formation at an early stage compared to previously studied Galactic supernova remnants – here we can probe timescales on the order of decades rather than centuries, filling in the gap between very young SNe (e.g., [43]) and historical remnants.

Thermal emission from small quantities ($10^{-4}\,M_{\odot}$) of dust was detected in the early days after the SN explosion (day$\sim$300–600) using MIR observations ([24, 77, 11, 111, 140]). More surprisingly, the Herschel Space Observatory ([105], hereafter Herschel) revealed a large amount of cold dust ($\sim0.5M_{\odot}$) at the location of the remnant ([80, 81]). ALMA resolved the emission from dust in SN 1987A on scales of $0.3^{\prime\prime}$ and confirmed that the $\sim0.5M_{\odot}$ of cold (20 K) dust discovered with Herschel originates from the inner ejecta region ([58, 81]). Dwek & Arendt [32] and Wesson et al. [135] re-visited the dust emission at early times ($<1200$ days). Dwek & Arendt [32] find that a large mass of dust can be present early on ($0.4M_{\odot}$ at $\sim615$ days) with a model of silicates and amorphous carbon. Wesson et al. [135] conclude instead, from comparing radiative transfer models to the optical–IR SED limits, that the dust mass increased more slowly over the first 10 years. This substantial mass of dust observed in the inner debris of SN 1987A demonstrates that a large fraction of the heavy elements ejected in a SN may be locked up in a dust reservoir.

Whether dust grains formed in the ejecta of a supernova are carbon or silicate-rich remains an unanswered question: the models of Cherchneff & Dwek [20, 21] and Sarangi & Cherchneff [115, 116] predict that for abundance ratios C/O $< 1$, carbon atoms will mostly be locked up in CO molecules in the first 1000 days, preventing the formation of a large mass of amorphous carbon dust. Though depending on the gas density, CO may be dissociated by electrons produced by radioactive decay Clayton [23], and/or (to a lesser extent) X-rays from the ring, however, the models of Sarangi & Cherchneff [115] and Sarangi & Cherchneff [116] 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. [81]), 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 [32]. 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 [111].

In this work, we present high angular resolution ALMA (Cycle 2) dust images for SN 1987A, where we resolve dust clumps on scales of $\sim80$ 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.

Data

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 15 m to 2.3 km. See Table 1 for a summary of the observations.

Sub-band		Frequency Range (GHz)	Baselines (m)	Angular Scales (arcsec)	Observation Date	SN Day	Bandpass Calibrator	Phase Calibrator	Check Source	Time (min)
B7	A1	299.88–315.87	45.4–1574.4	0.13–4.31	2015-06-28	10352	J0538-4405	J0635-7516	J0601-7036	18.4
	A2	299.88–315.87	43.3–2269.9	0.09–4.52	2015-09-22	10438	J0538-4405	J0635-7516	J0601-7036	18.3
	B	342.48–358.34	15.1–1574.4	0.11–11.46	2015-07-25	10379	J0538-4405	J0635-7516	J0601-7036	20.9
	C	346.23–362.09	15.1–1574.4	0.11–11.34	2015-07-25	10379	J0538-4405	J0635-7516	J0601-7036	19.9
	D1	303.62–319.48	45.4–1574.4	0.13–4.26	2015-06-28	10352	J0538-4405	J0635-7516	J0601-7036	18.8
	D2	303.62–319.48	43.3–2269.9	0.09–4.47	2015-09-22	10438	J0538-4405	J0635-7516	J0601-7036	18.8
B9	A	673.44–681.06	43.3–2269.9	0.04–2.10	2015-09-25	10441	J0522-3627	J0601-7036	J0700-6610	12.5
	B	680.94–688.56	43.3–2269.9	0.04–2.07	2015-09-25	10441	J0522-3627	J0601-7036	J0700-6610	12.5
	C	688.44–696.06	43.3–2269.9	0.04–2.05	2015-09-25	10441	J0522-3627	J0700-6610	J0450-8101	12.5

Each data set was reduced separately with Common Astronomy Software Applications package (CASA¹, McMullin et al. [88], 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 $\sim$(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 $\sim35\%$. Bands 7A and 7D also suffered from poor weather in the original June 2015 observations, with $\sim1.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.

Frequency Range (GHz)	$\nu_c$ (GHz)	Beam FWHM (arcsec$^2$)	Beam PA (deg)	RMS Noise (mJy bm$^{-1}$)
Continuum
224.00 – 227.00	225.50	0.30 × 0.30	0.00	0.12
238.00 – 243.00	240.50	0.30 × 0.30	0.00	0.09
246.00 – 249.00	247.50	0.30 × 0.30	0.00	0.10
269.50 – 270.50	270.00	0.30 × 0.30	0.00	0.18
278.00 – 280.00	279.00	0.30 × 0.30	0.00	0.10
306.06 – 307.47	306.76	0.20 × 0.15	124.32	0.07
311.88 – 319.48	315.68	0.19 × 0.14	119.11	0.04
673.45 – 685.00	679.22	0.08 × 0.06	74.37	0.71
Spectral Lines
CO $J=2\to1$	230.54	0.06 × 0.04	27.43	0.04
CO $J=6\to5$	691.47	0.09 × 0.07	185.35	2.82
SiO $J=5\to4$	217.10	0.06 × 0.04	19.74	0.05
SiO $J=6\to5$	260.52	0.04 × 0.03	173.66	0.06
SiO $J=7\to6$	303.93	0.13 × 0.10	35.47	0.47

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 $\pm1000$ km s$^{-1}$ expansion velocity of the ejecta means that the linewidth spans a substantial fraction of the observed bands. The continuum bands selected relative to the modelled molecular line emission are shown in Figure 1, using the ALMA spectra and the emission line model of CO, SiO, SO, and SO$_2$ from Matsuura et al. [82]. 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 Figure 2. We also utilize here the Cycle 2 Band–6 imaging data presented by Matsuura et al. [82], 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.

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\to5$ with rest frequency 691.47 GHz and SiO $J=7\to6$ at $\nu_{\mathrm{rest}}$ = 303.93 GHz. The 345.80 GHz CO $J=3\to2$ and 347.33 GHz SiO $J=8\to7$ lines were also 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.

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

The SiO $J=7\to6$ and CO $J=6\to5$ cubes were created with tclean in CASA, with a spectral resolution of 300 km s$^{-1}$, which gives a reasonable balance between velocity resolution and sensitivity per channel. CO $J=6\to5$ was imaged with natural weighting to maximize sensitivity per beam. SiO $J=7\to6$ was imaged with $\mathrm{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 Figure 3.

NOTE— The position angles are counter-clockwise from north. CO $J=2\to1$, SiO $J=5\to4$, and SiO $J=6\to5$ parameters are for the data cubes from Abellán et al. [1]. 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\to1$, SiO $J=5\to4$, and SiO $J=6\to5$² data as described in Abellán et al. [1]. Although both sets were taken in Cycle 2, the molecular line data presented by Abellán et al. [1] have higher signal to noise as they were combined with Cycle 3 data, and due to CO $J=6\to5$ being in a band with poorer atmospheric transmission than CO $J=2\to1$. The additional Cycle 3 data also give their CO $J=2\to1$, SiO $J=5\to4$, and SiO $J=6\to5$ 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$^{-1}$ receding from Earth (Gröningsson et al. [49]).²: The images for the SiO $J=6\to5$ transition were described but not shown in Abellán et al. [1].

Description of Images

The SN 1987A ring and ejecta continuum image at 315 GHz is shown in comparison to H$\alpha$ in Figure 2. These images have been aligned following a technique in Alp et al. [3] where the ring emission is used to derive a reference center, though here we take a simpler approach (see Appendix § 11 for details). Our derived ring+ejecta system center used in this work is \alpha=5^\mathrm{h}35^\mathrm{m}27\fs998, \delta=-69^\circ16^\prime11\farcs107 (ICRS), $\pm18$ 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 ([51, 38]). 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 § 5.3), and the inner region is thermal dust emission from the SN ejecta (Indebetouw et al. [58]). Figure 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 south-west. 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 Figure 2, where the $3\sigma$ contours highlighting the major continuum features are overlaid onto the continuum and H$\alpha$ images. The alignment accuracy is $\sim1$ pixel in the images, given the astrometric uncertainties discussed in § 2.1 (12 mas for Band-7 continuum) and Appendix § 11 (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 Figure 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 $\pm$ 1.2 kpc (Panagia [102]), 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 Figure 3 and the dust emission (red) in Figure 4 have pixel S/N>3, and the surrounding ejecta area has pixel S/N values of $\sim2$ in the 679 GHz image. The area between the ejecta and the ring – outside of the outermost ejecta contour in Figure 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\to1$ and SiO $J=5\to4$ images (middle left and lower left panels of Figure 3). This was first reported by Abellán et al. [1] and was seen both in the integrated 2D spatial maps and the 3D data cubes. Although the integrated SiO $J=6\to5$ map (middle panel of Figure 3 does not show the hole clearly in the same manner as SiO $J=5\to4$ and CO $J=2\to1$, 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\to5$ map. The CO and SiO molecular hole is just to the south of the ‘keyhole’ that is seen in H$\alpha$ (Fransson et al. [40]; top right panel of our Figure 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 $\sim50$ mas, or $\sim4\times$ the astrometric and alignment errors.

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

CO $J=6\to5$ has emission coincident with the CO $J=2\to1$ 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\to1$. The brightness peaks are distributed differently, and the hole is not visible in the integrated CO $J=6\to5$ 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\to6$ cannot be confirmed in these data, as the systemic line center ($v_{\mathrm{LSRK}}\sim300$ 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\to6$ 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 $\sim-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 Figure 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\to4$ channels near the systemic velocity. The brightest points of dust emission tend to coincide with relative depressions in the CO $J=6\to5$ brightness, giving the appearance of an anti-correlation between the main dust and CO $J=6\to5$ features. This is more clearly demonstrated in Figure 4, where the dust (red) and CO $J=6\to5$ (blue) images are overlaid. The individual images were normalized independently to emphasize the main features of each, with the visible colors shown roughly corresponding to areas of S/N$>3$. The gold and teal lines are guides highlighting the highest-S/N features of the dust and CO, respectively, in order to compare peaks in the emission. While there is some overlap in the faint features of the dust and CO, and the southern extent of the dust peak starts to fade to a combined magenta, the dust and CO $J=6\rightarrow5$ peaks do not generally overlap. Rather, the brightest dust features are located in areas of relatively faint CO $J=6\rightarrow5$ emission and vice-versa.

To test whether the apparent dust-CO anti-correlation is an artifact or result of the data reduction or continuum subtraction, we performed several checks. First, the Band–9 dust continuum was reconstructed in different ways, by imaging (CASA mfs-mode) in a variety of spectral windows and also by making a data cube across the entirety of Band–9 (including the CO line), and fitting the continuum emission. These techniques gave consistent results. Initially we 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 differences. This means that the structures seen in our final CO $J=6\to5$ map are robust to variations in how the continuum is determined and subtracted. The anti-correlation in CO $J=6\to5$ is visible, even before continuum subtraction, in the CO $J=6\to5$ dirty map (i.e., with no cleaning to deconvolve 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 function $\sum\operatorname{cov}_{XY}/\sqrt{\operatorname{cov}_{XX}\cdot\operatorname{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 \pm0.08$ – a moderate anti-correlation – using the accuracy estimate from, e.g., Frick et al. [41]. 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 $\pm1$.

We tested the robustness of the correlation measure by investigating different angular scales using the wavelet analysis described in Frick et al. [41] and Arshakian & Ossenkopf [5]. On scales of 1–2 beam widths (i.e., convolutions with kernels of those scales), the images start to become increasingly positively correlated, which is expected due to the peaks being separated by that amount ($\sim5\times$ the astrometry error). Below these scales, $r$ remains negative, so the anti-correlation is not sensitive to small changes in image resolution. Using the same wavelet analysis on the other images, the dust correlates more positively with the other CO and SiO lines than with CO $J=6\to5$. The standard correlation measures agree, with $r=+0.04$ for CO $J=2\to1$ and $r=+0.36$ for SiO $J=6\to5$.

The brightest dust feature is located one beam width northeast of the secondary CO $J=6\to5$ peak, and the brightest CO $J=6\to5$ 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 Figure 1), but there is some fainter (S/N$<5$) CO $J=6\to5$ emission overlapping with some of the dust emission. There is also low-level (S/N$<2$) dust emission that roughly spans the full extent of the ejecta. The other molecular species further complicate the picture, as noted earlier – the peak of SiO $J=7\to6$ is coincident with the dust peak and the SiO $J=5\to4$ hole (as seen in the bottom right panel of Figure 3). At the southern edge of the hole, some faint H$\alpha$ emission appears to be aligned with CO $J=2\to1$ in projection, but their velocity ranges differ. The 3-D view of H$\alpha$ (Larsson et al. [74], their Figure 6) indicates that the H$\alpha$ emission in this region peaks at velocities around –1500 km s$^{-1}$ while the peak CO $J=2\to1$ is between 0 and +200 km s$^{-1}$. While no velocity information is available for the dust continuum emission, it is spatially offset from the nearby H$\alpha$ and CO $J=2\to1$ by $\sim1$ dust resolution element. That is, the CO $J=2\to1$ and H$\alpha$ 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$\times$81 mas. The peaks of the ejecta CO and SiO emission are not cospatial with the peaks of the ejecta dust emission (with anti-correlated CO $J=6\to5$ and dust structures). The small peak/clump in the dust emission revealed in the Band 9 image (the blob) overlaps with holes previously observed in the lower line transitions of SiO and the CO molecular ejecta, and is coincident with some emission observed in the SiO $J=7\to6$ line.

The Spectral Energy Distribution of SN 1987A

The three physical mechanisms primarily responsible for emission in the FIR–radio portion of the continuum are thermal IR greybody emission from dust (from the ejecta region – Matsuura et al. [80]; Indebetouw et al. [58]), non-thermal radio/mm synchrotron emission (from the ringManchester 2007; Potter et al. 2009; Zanardo et al. 2010; Lakicević 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 146 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.

Photometry

The continuum bands are defined as those frequencies that are molecular line free, as demonstrated in Figure 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 § 3 (Appendix § 11), 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_{\mathrm{MAJ}} \times r_{\mathrm{MIN}} = 0.42 \times0.36$ and major axis P.A. of $25^\circ$ (N through E). For the 315 GHz circumstellar ring, an elliptical annulus with inner $r_{\mathrm{MAJ}}$ and $r_{\mathrm{MIN}}$ of $0.45 \times0.42$, outer $r_{\mathrm{MAJ}}$ and $r_{\mathrm{MIN}}$ of $1.45 \times1.35$, and $85^\circ$ 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.

$\nu_c$ (GHz)	$\Delta\nu$ (GHz)	$S_\nu$ (mJy)	RA (deg)	DEC (deg)	$R_{\mathrm{MAJ}}^\dagger$ (arcsec)	$R_{\mathrm{MIN}}^\dagger$ (arcsec)	P.A. (deg)
Ejecta
225.50	3.0	$0.5 \pm0.2^{+0.3}_{-0.0}$	83.866586	-69.269733	0.32	0.29	125
240.50	5.0	$0.7 \pm0.2^{+0.4}_{-0.1}$	83.866600	-69.269720	0.37	0.36	75
247.50	3.0	$1.1 \pm0.4^{+0.6}_{-0.1}$	83.866639	-69.269739	0.41	0.35	125
270.00	1.0	$1.2 \pm0.6^{+0.7}_{-0.1}$	83.866652	-69.269709	0.40	0.31	60
279.00	2.0	$1.3 \pm0.4^{+0.8}_{-0.1}$	83.866662	-69.269747	0.41	0.40	125
306.76	1.4	$1.8 \pm0.6^{+1.1}_{-0.2}$	83.866667	-69.269753	0.42	0.36	125
315.68	7.6	$1.8 \pm0.6^{+1.1}_{-0.2}$	83.866667	-69.269753	0.42	0.36	125
679.22	11.5	$36.2 \pm7.2^{+7.2}_{-7.2}$	83.866728	-69.269740	0.42	0.36	125
Ring
225.50	3.0	$13.0 \pm0.5^{+7.8}_{-1.3}$	83.866585	-69.269731	1.45, 0.35	1.35, 0.33	175
240.50	5.0	$12.5 \pm0.5^{+7.5}_{-1.3}$	83.866600	-69.269722	1.45, 0.42	1.35, 0.39	175
247.50	3.0	$13.4 \pm0.8^{+8.1}_{-1.3}$	83.866630	-69.269736	1.45, 0.43	1.35, 0.40	175
270.00	1.0	$13.7 \pm1.1^{+8.2}_{-1.4}$	83.866592	-69.269720	1.45, 0.45	1.35, 0.42	175
279.00	2.0	$12.9 \pm0.8^{+7.7}_{-1.3}$	83.866661	-69.269747	1.45, 0.44	1.35, 0.41	175
306.76	1.4	$12.0 \pm1.5^{+7.2}_{-1.2}$	83.866666	-69.269753	1.45, 0.45	1.35, 0.42	175
315.68	7.6	$11.3 \pm1.7^{+6.8}_{-1.1}$	83.866667	-69.269753	1.45, 0.45	1.35, 0.42	175
679.22	11.5	$19.4 \pm19.3^{+3.9}_{-3.9}$	83.866720	-69.269737	1.45, 0.45	1.35, 0.42	175
Total System
225.50	3.0	$13.5 \pm0.5^{+8.1}_{-1.3}$	83.866585	-69.269731	1.45	1.35	175
240.50	5.0	$13.2 \pm0.5^{+7.9}_{-1.3}$	83.866600	-69.269722	1.45	1.35	175
247.50	3.0	$14.6 \pm0.7^{+8.7}_{-1.5}$	83.866630	-69.269736	1.45	1.35	175
270.00	1.0	$15.4 \pm0.8^{+9.2}_{-1.5}$	83.866592	-69.269720	1.45	1.35	175
279.00	2.0	$14.2 \pm0.8^{+8.5}_{-1.4}$	83.866661	-69.269747	1.45	1.35	175
306.76	1.4	$13.9 \pm1.7^{+8.3}_{-1.4}$	83.866666	-69.269753	1.45	1.35	175
315.68	7.6	$13.3 \pm1.8^{+8.0}_{-1.3}$	83.866667	-69.269753	1.45	1.35	175
679.22	11.5	$55.5 \pm20.9^{+11.1}_{-11.1}$	83.866720	-69.269737	1.45	1.35	175

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_{\mathrm{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 $\sim15\%$, 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 § 12 for more details.

The ejecta, ring and total system flux densities and uncertainties are shown in Figure 5 as gold circles, green rings, and purple diamonds, respectively, and their values are listed in Table 3. Previous measurements are also shown in Figure 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 § 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 (Lakicević et al. [70]), and 102 GHz (Zanardo et al. [146]) is dominated by synchrotron emission from the ring. The total emission at 213, 345, & 672 GHz (Zanardo et al. [146]), on the other hand, gradually consists of a higher and higher fraction of thermal emission until that is dominant in the submm

and the FIR. As the synchrotron brightness increases in time ([123, 18]), these literature flux densities were scaled to their levels at day 9280 by Zanardo et al. [146] 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.

Modified Blackbody Fits to the Ejecta Dust Emission

Description of the modified blackbody fits

Figure 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. [146] and Matsuura et al. [81].

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. [81] 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 $\mu\mathrm{m}$, as it is possibly contaminated by warm ring dust (Matsuura et al. [83]) and/or [O I] 63 $\mu\mathrm{m}$ emission; and 500 $\mu\mathrm{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 $\sim1300$ 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}\mathrm{Ti}$ decay, which has an estimated lifetime of 85 years ([2, 64], also see later Matsuura et al. [80]). 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 $\sim0.2$K, 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}\mathrm{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. [36]); 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_{\nu}(\lambda)=M_{\mathrm{dust}}\kappa_{\mathrm{abs}}(\lambda)B_{\nu}(\lambda)/d^{2}$, where $S_{\nu}(\lambda)$ is flux density, $\kappa_{\mathrm{abs}}(\lambda)$ is the mass absorption coefficient of dust grains, $B_{\nu}(\lambda)$ is the Planck function and $d$ is the distance to SN 1987A. We assume that the emission is optically thin across this wavelength range. $\kappa_{\mathrm{abs}}(\lambda)$ 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 $\rho$, $\kappa_{\mathrm{abs}}(\lambda)$ is defined as $\kappa=3/(4\rho a)Q$, where $Q$ is the absorption efficiency of the dust, which can be calculated from optical constants with Mie theory. $\kappa_{\mathrm{abs}}$ is often assumed to be a power law defined as $\kappa_{0}(\lambda/\lambda_{0})^{-\beta}$ (where $\kappa_{0}$ is the reference value of $\kappa_{\mathrm{abs}}(\lambda)$ at wavelength $\lambda_{0}$ and $\beta$ is the power law emissivity index of $\kappa_{\mathrm{abs}}(\lambda)$).

Fitting a dust greybody to the FIR-mm SED with the power law approximation to $\kappa_{\mathrm{abs}}$ gives fit parameters of $\beta=2.05^{+0.11}_{-0.10}$, $M_{\mathrm{dust}}=1.53^{+0.13}_{-0.13}\,M_{\odot}$, and $T_{\mathrm{d}}=17.83^{+0.60}_{-0.57}$ K if $\kappa_{0}(850\mu\mathrm{m})=0.07\ \mathrm{m}^{2}\,\mathrm{kg}^{-1}$ (using the empirical measurement of $\kappa_{\mathrm{abs}}$ assuming the fraction of metals locked in dust is constant across the local ISM, James et al. [62]). 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 $\kappa_{\mathrm{abs}}(\lambda)$, following Indebetouw et al. [58] and Matsuura et al. [81].

We perform dust greybody fits to the FIR-mm SED using a wide variety of $\kappa_{\mathrm{abs}}(\lambda)$ profiles directly obtained from the literature or calculated using Mie theory; in total we used 134 $\kappa_{\mathrm{abs}}(\lambda)$ profiles. These include: Weingartner & Draine [133] LMC average (as an approximation to the conditions near SN 1987A in the LMC), Demyk et al. [27] amorphous silicate samples at 30K, Ormel et al. [98] calculations of bare and icy silicate+graphite grains, Ossenkopf & Henning [99] bare and icy mant­le grains (protostellar core coagulation models), and using Mie scattering calculations for amorphous carbon (amC, Zubko et al. [149]; Rouleau & Martin [112]), graphite (Draine & Lee [29]), “cosmic silicates” (amorphous FeMgSiO₄, Jaeger et al. [59]), other silicates including enstatite and forsterite from Henning & Stognienko [53], pure iron (Henning & Stognienko [53]), and 68 profiles from the Jena database (Henning et al. [55]) of minerals including silicates, amorphous carbons, carbides, oxides, and sulfides. Many of the models in the Jena database only extend to 500 $\mu\mathrm{m}$ or 1000 $\mu\mathrm{m}$. In these cases we extrapolate to our longest ALMA wavelength of 1329 $\mu\mathrm{m}$ in log space (as lines, as most models follow power laws in the sub-mm), from the last 300–500 $\mu\mathrm{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 [32]), with a carbon volume filling factor $f_{\mathrm{C}}$ of 18%.

For Mie theory calculations, we adopt a grain radius of $a=0.01\ \mu\mathrm{m}$ for PAHs, and $a=0.1\ \mu\mathrm{m}$ for all other models. For most dust models, Mie-derived FIR $\kappa_{\mathrm{abs}}$ curves are

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 percentile of the distributions of fits from 1000 samplings of the observed flux densities. For $\kappa(\lambda)$ models derived from Mie theory, grains of radius $a=0.1\mu$m were assumed except for the case of PAHs, where $a=0.01\mu$m was used. Grain densities are given for the Mie and CDE cases: amorphous carbon from Zubko et al. [150] (and Rouleau & Martin [112] for their AC1 sample), graphite and SiC from Laor & Draine [72], PAHs from Li & Draine (2001), stoichiometric varieties of olivines and pyroxenes from Henning & Stognienko [53], silicate/carbon composite CDE with carbon volume filling factor $f_{\mathrm{C}}=18\%$ from Dwek & Arendt [32], jadeite from Mutschke et al. [94], and in general from the Jena optical constants database⁴ (Henning et al. [55]). $\kappa_{850}$, the mass absorption constant at 850 $\mu$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 $\chi_{\nu}^{2}<2$.

nearly identical for grain radii of $a < 5\,\mu\mathrm{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. [149] amorphous carbons (amCs), Jäger et al. [61] amorphous silicates, and iron from Henning & Stognienko [53], among others. The fitted masses and temperatures for this sample are listed in Table 4.

Dust	Reference	Grain Density (g cm$^{-3}$)	$M_{\mathrm{dust}}$ ($M_{\mathrm{dust}}$)	T (K)	$\kappa_{850}$ (m$^2$ kg$^{-1}$)	Good Fit
amC (ACH2 sample), Mie 0.1$\mu$m	Zubko et al. [149]	1.81	$1.46^{+0.09}_{-0.08}$	$17.5^{+0.1}_{-0.1}$	0.087	Y
amC (ACAR sample), Mie 0.1$\mu$m	Zubko et al. [149]	1.81	$0.38^{+0.02}_{-0.02}$	$22.0^{+0.2}_{-0.2}$	0.254	N
amC (BE sample), Mie 0.1$\mu$m	Zubko et al. [149]	1.81	$0.77^{+0.05}_{-0.04}$	$20.7^{+0.2}_{-0.2}$	0.141	N
amC (AC1 sample), Mie 0.1$\mu$m	Rouleau & Martin [112]	1.85	$0.43^{+0.03}_{-0.03}$	$21.6^{+0.3}_{-0.3}$	0.203	Y
Cellulose (800K sample), Mie 0.1$\mu$m	Jager et al. [60]	1.81	$0.46^{+0.03}_{-0.03}$	$18.8^{+0.2}_{-0.2}$	0.178	Y
Graphite, Mie 0.1$\mu$m	Draine & Lee [29]	2.26	$1.62^{+0.11}_{-0.10}$	$17.8^{+0.2}_{-0.2}$	0.069	Y
PAH (neutral), Mie 0.01$\mu$m	Laor & Draine [72]	2.24	$1.69^{+0.11}_{-0.10}$	$18.0^{+0.2}_{-0.2}$	0.071	Y
Silicate – Enstatite, Mie 0.1$\mu$m	Jäger et al. [61]	2.71	$4.10^{+0.23}_{-0.24}$	$18.0^{+0.2}_{-0.2}$	0.029	Y
Silicate – Forsterite, Mie 0.1$\mu$m	Jäger et al. [61]	3.2	$4.03^{+0.26}_{-0.25}$	$17.9^{+0.2}_{-0.2}$	0.029	Y
Silicate – “Cosmic”, Mie 0.1$\mu$m	Jaeger et al. [59]	3.2	$3.46^{+0.24}_{-0.22}$	$17.7^{+0.2}_{-0.2}$	0.034	Y
Silicate/Carbon Composite CDE – $f_{\mathrm{C}}=0.18$	Dwek & Arendt [32]	2.95	$0.38^{+0.02}_{-0.02}$	$21.1^{+0.2}_{-0.2}$	0.270	N
Silicate – LMC Average	Weingartner & Draine [133]	$\cdots$	$2.49^{+0.15}_{-0.14}$	$18.0^{+0.2}_{-0.2}$	0.047	Y
Silicate – 30K Average	Demyk et al. [27]	$\cdots$	$0.27^{+0.02}_{-0.02}$	$18.5^{+0.2}_{-0.2}$	0.265	Y
Silicate – Composite Aggregate	Semenov et al. [118]	$\cdots$	$35.06^{+2.18}_{-2.18}$	$21.4^{+0.2}_{-0.2}$	0.003	Y
Silicate – Porous Multilayer Spheres	Semenov et al. [118]	$\cdots$	$6.13^{+0.36}_{-0.37}$	$29.4^{+0.4}_{-0.4}$	0.010	N
Silicate – Bare Grains, 0.03 Myr	Ormel et al. [98]	$\cdots$	$0.50^{+0.03}_{-0.03}$	$20.7^{+0.3}_{-0.2}$	0.180	Y
Silicate – Icy Grains, 0.03 Myr	Ormel et al. [98]	$\cdots$	$0.60^{+0.04}_{-0.04}$	$18.3^{+0.2}_{-0.2}$	0.184	Y
Silicate – Naked Grains, $n_{\mathrm{H}}=10^{5}$	Ossenkopf & Henning [99]	$\cdots$	$0.64^{+0.04}_{-0.04}$	$20.9^{+0.2}_{-0.2}$	0.163	N
Silicate – Thin Ice Mantles, $n_{\mathrm{H}}=10^{5}$	Ossenkopf & Henning [99]	$\cdots$	$0.74^{+0.04}_{-0.04}$	$19.1^{+0.2}_{-0.2}$	0.142	Y
SiC, Mie 0.1$\mu$m	Pegourie [103]	3.22	$1.57^{+0.08}_{-0.09}$	$23.3^{+0.2}_{-0.2}$	0.058	N
FeS, Mie 0.1$\mu$m	Henning & Stognienko [53]	4.83	$0.68^{+0.04}_{-0.04}$	$34.7^{+0.6}_{-0.6}$	0.086	N
FeO, Mie 0.1$\mu$m	Henning et al. [52]	5.7	$0.28^{+0.01}_{-0.02}$	$28.2^{+0.4}_{-0.3}$	0.259	N
SiO$_2$, Mie 0.1$\mu$m	Henning & Mutschke [54]	2.196	$4.41^{+0.39}_{-0.33}$	$17.5^{+0.2}_{-0.2}$	0.022	N
TiO$_2$, Mie 0.1$\mu$m	Posch et al. [106]	3.78	$81.77^{+5.10}_{-4.49}$	$17.7^{+0.2}_{-0.2}$	0.001	Y
Al$_2$O$_3$ “Compact” sample, Mie 0.1$\mu$m	Begemann et al. [8]	3.2	$0.90^{+0.06}_{-0.06}$	$19.0^{+0.2}_{-0.2}$	0.112	Y
NaAlSi$_2$O$_6$, Mie 0.1$\mu$m	Mutschke et al. [94]	2.4	$0.10^{+0.01}_{-0.01}$	$23.1^{+0.3}_{-0.3}$	0.982	Y
MgAl$_2$O$_4$, Mie 0.1$\mu$m	Fabian et al. [35]	3.64	$121.74^{+7.54}_{-7.48}$	$17.7^{+0.2}_{-0.2}$	0.001	Y
Pure Iron, Mie 0.1$\mu$m	Henning & Stognienko [53]	7.87	$3.97^{+0.28}_{-0.23}$	$19.9^{+0.2}_{-0.2}$	0.025	Y

To determine whether the SED fit is ‘good’, we set a quality threshold whereby $\chi_\nu^2 < 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_\nu^2 < 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 § 7.1.

Results of the dust fits

Figure 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.) Figure 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 $\sim18$–$23\,\mathrm{K}$ (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_{\mathrm{abs}}$ values and thus return lower inferred dust masses, $M_{\mathrm{dust}} \sim0.4$–$0.8M_\odot$ (Table 4). Amorphous silicates from Demyk et al. [27], which have more than a factor of 10 larger $\kappa_{850}$ compared to those from Semenov et al. [118], also yield a relatively moderate mass of $\sim0.27M_\odot$. The silicate models of Ormel et al. [98] and Ossenkopf & Henning [99] 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 $\mu\mathrm{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 $\sim0.5$–$0.74\,M_\odot$ for Ormel et al. [98] and Ossenkopf & Henning [99]. 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. [59] or Jäger et al. [61], require even larger masses of $\sim1$–$4\,M_\odot$ to fit the observed SED. The CDE composite model of Dwek & Arendt [32] results in a moderate mass of $0.38M_\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⁵, specifically the value of the mass absorption coefficient ($\kappa_{\mathrm{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 (Figure 7). The inferred masses very clearly follow an inverse linear relation based on the submm $\kappa_{\mathrm{abs}}$ value where $M_{\mathrm{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.

The Ring

The ring flux densities are shown in Figure 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 (Lakicevic et al. [70]). 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. [146]).

The Cycle 2 ring flux densities exhibit more scatter and are lower than the ATCA values by $\sim30\%$. 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. [146]) and ALMA Cycle 0 data (Indebetouw et al. [58]) to be $-0.73 \pm0.02$ at day 9280 (Zanardo et al. [146]). An updated fit, using the ALMA Cycle 2 ring flux densities and the more recent 1–9 GHz measurements from Callingham et al. [17] from day 9686, results in a slightly lower $\alpha= -0.70 \pm0.06$.

The radio ring emission has steadily increased due to the synchrotron-producing electrons being accelerated by expanding shockwaves (Zanardo et al. [144]; Staveley-Smith et al. [123]). Recently, Cendes et al. [18] have fit the 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. [18] model prediction for day 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. [12]; Matsuura et al. [83]).

Analysis and Results of Molecular Lines

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\to4$ and CO $J=2\to1$ images exhibit a hole (Figure 3) where the dust emission peaks, and the SiO line ratio indicates lower SiO $J=5\to4$ brightness with respect to SiO $J=6\to5$ and SiO $J=7\to6$ (seen in the SiO line ratio maps in Figure 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. [130]).

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 [100]) 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 LAMDA database (Schöier et al. [117]), and we use H$_2$–SiO collisional cross sections based on the calculations by Dayou & Balança [25]. In the ISM, H$_2$ is widely assumed to be the dominant collisional partner in molecular clouds; however, that is probably not the case for the ejecta of SNe. As a consequence of a series of nuclear burning processes, the progenitor star’s core will have built up layers of different newly synthesised elements, with hydrogen being depleted. In the two layers containing abundant Si, the major elements are O and S (e.g., Woosley [141]), and the collisional partner of SiO is likely to be O$_2$ and SiS. This would potentially change the collisional cross section by a factor of 1–10, depending on the transitions (Matsuura et al. [82]).

One of the RADEX input parameters is the FWHM of the line width of the Gaussian ($\Delta v$). We adopted a $\Delta v$ of 400 km s$^{-1}$, whose integrated area over the Gaussian profile would be equivalent to that of a box-shaped 300 km s$^{-1}$ line profile. If the line is optically thin, the assumed line profile is not a major issue. However, as we will see later in this analysis, the lines are mildly optically thick at the line center, but not at the side of the line profile, so we therefore make the assumption that the line profile only moderately affects the line ratios and the line intensities in this “mildly” optically thick regime.

The main parameters involved in RADEX calculations are the kinetic temperature ($T_{\rm kin}$), the density of the collisional partner ($n_{\rm coll}$) and the column density ($N_{\rm SiO}$). In the optically thin regime, as found in the calculated parameter range where solutions are found, the SiO line intensities are determined by $N_{\rm SiO}$ together with the area filling factor and the expansion velocity $v_{\rm exp}$, while the SiO line ratios are determined by $T_{\rm kin}$ and $n_{\rm coll}$, independent of $N_{\rm SiO}$. The filling factor is defined as the area of the line emitting fraction within the beam/pixel, following Goldsmith & Langer [45]. Previous analyses of SN 1987A with lower angular resolution suggested the range of 2.5–45% (Kamenetzky et al. [67]; ); Matsuura et al. [82]); thus, we assume that the filling factor of 1–50% is a reasonable range. In this analysis, we adopted a column density grid of $N_{\rm SiO}=10^{12}$–$10^{18}$ cm$^{-2}$ in factor of 10 increments, and searched for the predicted line intensities that can match the measured ones within the assumed range of the filling factor. The adopted ranges of the parameters are $n_{\rm coll}=10^{3}$–$10^{10}$ cm$^{-3}$ and $T_{\rm kin}=10$–200 K. Matsuura et al. [82] suggested the temperature range is below 190 K for SiO, with the CO kinetic temperature between 30–50 K, so we restricted the analysis to below 200 K. Although 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\to4$ and SiO $J=6\to5$ were convolved to match the lowest spatial resolution – that of the SiO $J=7\to6$ beam – with a uniform pixel width of 0\farcs015. The convolved and regridded line ratio maps, made on a channel-by-channel basis, are displayed in Figure 8. Intensities were averaged over $5\times5$ pixels in order to increase the signal to noise ratio in the SiO $J=7\to6$ image. As the minor-axis FWHM of the beam is 0\farcs17 for SiO $J=7\to6$, there is a small loss of spatial information by this averaging.

One caveat: the continuum subtraction for SiO $J=6\to5$ was performed on the final imaged data cube, whereas for SiO $J=5\to4$ and SiO $J=7\to6$ 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\to5$ 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\to5$ 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\to5$ and SiO $J=5\to4$, and 5% for SiO $J=7\to6$ 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\to5$ were higher than the measured value, the line ratio of SiO $J=6\to5$/SiO $J=5\to4$ would increase, but the ratio of SiO $J=7\to6$/SiO $J=6\to5$ 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.

Results of SiO analysis

We selected two representative spatial regions, shown in Figure 9, in order to understand the change in physical parameters ($T_{\rm kin}$, $n_{\rm coll}$) within the ejecta, which in turn affect the SiO emission. One is close to the molecular hole seen in SiO $J=5\to4$ and CO $J=2\to1$ 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\to4$ line intensity, so we chose the nearest possible point for region A. All ratios in this analysis were determined from the 0 km s$^{-1}$ LSRK channel, which has good S/N for all three SiO lines, and the hole is clearly identified.

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

Figure 10 demonstrates the plausible ranges for the kinetic temperature ($T_{\rm kin}$) and collisional partner’s density ($n_{\rm coll}$) at region A, for $N_{\rm SiO}=10^{16}\ {\rm cm}^{-2}$. The collisional partner for the RADEX calculations is H$_2$; see the last paragraph of this section for discussion of this point. In Figure 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. Figure 10 (d) compares the possible ranges of the kinetic temperature ($T_{\rm kin}$) and collision partner density ($n_{\rm coll}$) from the three SiO line ratios. The uncertainty of the SiO $J=7\to6$/SiO $J=5\to4$ ratio has the tightest constraint, so only the uncertainty of this ratio is plotted in Figure 10 (d). Figure 10 (d) shows that the ratios of all three SiO lines can be matched with a very similar set of $T_{\rm kin}$ and $n_{\rm coll}>10^{7}\ {\rm cm}^{-3}$, as indicated by the solid lines of the three different ratios almost overlapping each other. There are multiple scenarios that match the measured SiO line ratios: the first possibility is $T_{\rm kin}=34$ K with LTE conditions at $n_{\rm coll}>10^{7}\ {\rm cm}^{-3}$, the second possibility is a non-LTE condition with a range of ($10^{6}<n_{\rm coll}<10^{7}\ {\rm cm}^{-3}$) and $50<T_{\rm kin}<200$ K, and finally the curves connecting these two conditions. The line center turns optically thick at the column density of $N_{\rm SiO}=10^{16}\ {\rm cm}^{-2}$, but the majority of the lines at off-center velocities are optically thin, so the overall analysis is not strongly affected by optical thickness at this column density. The filling factor at this column density is 9–14%.

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

In the non-LTE range, the required $n_{\mathrm{coll}} > 10^{7}\,\mathrm{cm}^{-3}$ is more or less comparable to that for $N_{\mathrm{SiO}} = 10^{16}\,\mathrm{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_{\mathrm{SiO}} = 10^{16}\,\mathrm{cm}^{-2}$ case.

At region B, the physical parameters are slightly different from those at the hole (region A). Figure 12 shows the possible parameter space which fits the SiO ratios for region B; only $N_{\mathrm{SiO}} = 10^{16}\,\mathrm{cm}^{-2}$ gives a feasible range. The most plausible temperature is $T_{\mathrm{kin}} = 18$ K for LTE conditions ($n_{\mathrm{coll}} > 10^{7}\,\mathrm{cm}^{-3}$), and an alternative possibility is non-LTE with $3 \times10^{5} < n_{\mathrm{coll}} < 10^{7}\,\mathrm{cm}^{-3}$ and $18 < T_{\mathrm{kin}} < 200\,\mathrm{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 \rightarrow4$ and $J=6 \rightarrow5$ 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 & Cherchneff [115]). 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 § 8.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.

CO analysis and results

Although the CO $J=6 \rightarrow5$ 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 Figure 13 shows the location of 20 regions selected across the CO-emitting ejecta, overlaid on top of the CO $J=2 \rightarrow1$ emission. Nine regions, each the size of the CO $J=6 \rightarrow5$ beam, are highlighted as areas of interest, potentially probing different conditions (numbered 1–9 in rows from left to right). The middle panel of Figure 13 compares the summed spectra of CO $J=6 \rightarrow5$ and CO $J=2 \rightarrow1$ (with the latter convolved to the CO $J=6 \rightarrow5$ beam before integration) for these 9 regions showing their location with respect to the CO $J=2 \rightarrow1$ emission. The spectra are in units of mJy, having been spatially integrated over each region.

The bottom panel of Figure 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 of the CO $J=2 \rightarrow1$ and CO $J=6 \rightarrow5$ transitions are similar to each other in the scaled spectra (see for example regions 6, 7, 9 in Figure 13). However, in regions 1 and 8, the CO $J=6 \rightarrow5$ profile is suppressed at negative velocities with respect to the CO $J=2 \rightarrow1$ line. The CO $J=6 \rightarrow5$ emission is also suppressed with respect to CO $J=2 \rightarrow1$ 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 \rightarrow5$ emission at velocities of $-1000$–$+1000\,\mathrm{km\,s}^{-1}$ with respect to CO $J=2 \rightarrow1$ (bottom panel Figure 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 reprepresentative of the general SiO emitting properties across the molecular ejecta based on the SiO line ratios, whereas Figure 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. Figure 14 shows the resulting RADEX calculations of the CO $J=6\to5/2\to1$ 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$^{-2}$ used in the SiO ratio calculations. This is because the CO spectra are compared in velocity space, whereas we use integrated line intensities for the SiO analysis.⁶

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

Discussion

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 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$_2$ (Posch et al. [106]) which requires 82 $M_\odot$ of dust, or MgAl$_2$O$_4$ Fabian et al. [35] which requires 122 $M_\odot$ for example, are clearly untenable as they are larger than the progenitor star mass (18–20 $M_\odot$, Woosley [141]). This rules out a further three dust varieties in Table 4 (all of which satisfy our $\chi_\nu^2$ criteria for a good fit).

Simple upper limits to the dust masses of various dust species can be estimated by calculating the resulting mass for 100% of the elements in the nucleosynthesis models being locked into dust, ignoring all chemistry, mixing, and other physical limitations. Such a highly unrealistic scenario is useful for winnowing out dust models that yield dust masses that are too large. Considering the 18 $M_\odot$ $Z = 0.008Z_\odot$ progenitor model from Nomoto et al. [97], 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. [97] 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 [112], 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. [97] 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. [97] 18 $M_\odot$ nucleosynthesis model is 0.079 $M_\odot$ (${}^{56}\mathrm{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 [53]), which requires 3.97 $M_\odot$ of dust to fit the SED. The Woosley & Weaver [142] 15 $M_\odot$ and 18 $M_\odot$ $Z = 0.1Z_\odot$ progenitor models predict iron masses ranging from 0.14–0.20 $M_\odot$, where roughly half of the iron originates from ${}^{56}\mathrm{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, [53, 52]) or resulting in unrealistically high dust masses for the pure iron model.

Limits from the Nomoto et al. [97] 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. [81]. 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 [141] 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 $\sim$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. [27] model) and amorphous carbons as also concluded by Matsuura et al. [81], limiting the total dust mass to potentially $< 0.7M_\odot$.

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

Our spatially resolved images (Figure 3) show that the dust distribution is clumpy and asymmetric. Comparing the spatial distribution of the dust with the CO $J=6\to5$ reveals a weak anti-correlation with the integrated CO $J=6\to5$ 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\to5$ is suppressed compared with CO $J=2\to1$ 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\to1$ coincides with the dust peak, with relatively strong CO $J=6\to5$ emission observed at $-1000$–$1000\ \mathrm{km\,s^{-1}}$ with respect to CO $J=2\to1$, and this strong CO $J=6\to5$ continuing to its neighboring region (region 5 in Figure 13). We will discuss this region separately in § 8.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., [115]). The outermost region is the H-envelope, followed by the He shell, which can contain more carbon atoms than oxygen atoms ([142, 108]). The He shell can also form CO, as it enriches with C and O (e.g., [115]). 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 ${}^{56}\mathrm{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\to5$ 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 [116] 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. [28] and Clayton [23] 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 [23] calculations involved only a few reaction rates involving C and CO. In contrast, Sarangi & Cherchneff [115] and Sarangi & Cherchneff [116] 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 \times10^{-3}\ M_{\odot}$ of carbonaceous dust out of total of $0.04\ M_{\odot}$ dust mass for a 19 $M_{\odot}$ star). One simple explanation is that Sarangi & Cherchneff [115] and Sarangi & Cherchneff [116] 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, suggesting 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, particularly involving the highly energetic electrons. Together with the recent detection of HCO$^{+}$ (Matsuura et al. [82]), which was largely under-predicted in the chemical model of Sarangi & Cherchneff [115], 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 allow Si and C dust to form in the same regions. Indeed, some 3D hydrodynamical models (e.g., B15, N20; see Utrobin et al. [129]) suggest that $\sim30$–$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 temperature is the consequence of dust formation. Second, the temperature of these regions may have already been lower in the early days when the dust grains were formed, and the temperature reached the dust sublimation temperature while the gas density was reasonably higher, driving more efficient dust condensation. In this case, dust formation is the consequence of cooler gas temperature. Currently, the data do not provide a way to distinguish between these two cases.

We note that we see faint diffuse dust emission that might be more extended beyond the Band-9 dust structure in Figure 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 of 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. [1]). Aside from SN 1987A, we know from the knots and filaments observed across the Galactic SNR Cassiopeia A (e.g., Milisavljevic & Fesen [91]) that the gas ejected in SN explosions can be clumpy. Interestingly, Sarangi & Cherchneff [116] included clumpiness in their chemical models and showed that clumpy gas, compared with smoothly distributed gas, provides density enhancements in the ejecta resulting in larger grain sizes for SN-formed dust.

Interpreting the Bright Point Source (the Blob)

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\to1$ and SiO $J=5\to4$ line-emitting ejecta, and we also see emission located in the hole in the SiO $J=7\to6$ 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 observations (Figure 4 in [1]) show, that in the explosions 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 compresses the material such that the entropy increases significantly and an outwards moving shock forms. This feature is termed the ‘self-reflected’ shock, first discussed in Ertl et al. [33] (with forthcoming details in Gabler et al. in preparation). The self-reflected shock accelerates the innermost ejecta compared to homologous expansion and leaves a region with lower density in the center. When passing more and more of the ejecta, the shock loses energy. It finally dissipates and cannot accelerate the outermost ejecta. This acceleration of the low-velocity ejecta leads to a formation of a higher density shell-like configuration as observed in the emission of the transition lines SiO $J=5\to4$ or CO $J=2\to1$ (Figure 4 in [1]). However, in the hydrodynamic simulations, one model based on the progenitor model from Shigeyama & Nomoto [119] (model N20 in [138, 1]) 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 ([128, 129]) disfavors that model.

The blob: gas and dust heated by the compact source

To explain the higher SiO and CO gas temperature and increased brightness for the observed properties of the blob, we suggest two possibilities: (i) gas heated by a compact object or (ii) a clump heated by ${}^{44}\mathrm{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 radiation from the compact object, with an early development of a pulsar wind nebula.

${}^{44}\mathrm{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 ([64, 74]). SN explosion models show that ${}^{44}\mathrm{Ti}$ is located more or less at a similar radial extent as ${}^{56}\mathrm{Ni}$ and ${}^{28}\mathrm{Si}$ ([142, 138, 139]) with almost identical bulk velocities, though the modeled distribution of ${}^{44}\mathrm{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}\mathrm{Ti}$ ejecta is redshifted, suggesting that the bulk is moving away from the observer ([10]). Since the predicted distribution of ${}^{44}\mathrm{Ti}$ shows qualitatively similar properties as Figure 4 in [1], because it is subject to effectively the same hydrodynamical history (Jerkstrand et al., in preparation), it is unlikely that gas heated by ${}^{44}\mathrm{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 ([3]), 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 ([125, 47, 101]). G54.1+0.3 ([127, 109]), G11.2−0.3, G21.5−0.9 and G29.7−0.3 ([19]). In the Crab Nebula, the radiation from the pulsar wind can partially dissociate gas, contributing to the formation of the molecules OH$^{+}$ and ${}^{36}\mathrm{ArH}^{+}$ ([7]). 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 ([13, 9]), instead of a neutron star, jets from the black hole can also heat up gas and dust locally (e.g., [114]).

Instead of synchrotron, direct thermal radiation from the neutron star can also heat the gas and dust locally ([3]). 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 [87, 4]. 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. [68]) and the redshifted $^{44}\mathrm{Ti}$ spectrum (Boggs et al. [10]), 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. [63]). 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. [3]: $\alpha= 05^\mathrm{h}35^\mathrm{m}27^\mathrm{s}.9875$, $\delta= -69^\circ16'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. [3]. This offset is $\sim$3–5 times the total alignment uncertainty, however, it is still within the range of the neutron star kick. Zanardo et al. [146] 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. [3]. 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. [3]).

Moreover, the location of the compact object and the brightest part of the pulsar wind in the Crab Nebula are also not coincident ([134, 47]). 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 $\sim$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. [149] or to 26 K from 18 K for Jäger et al. [61] dust grains would be required to explain the peak).

We estimate the dust peak flux density as $S_{679\mathrm{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 [14]), 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. [149] ACAR model and a temperature of 33 K, we obtain the light orange shaded region in Figure 15. Alp et al. [3] 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_{\mathrm{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}\mathrm{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}\mathrm{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.

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., [86, 73, 39, 74, 75]). 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_{\mathrm{ext}}$ (related to, but different from the $\kappa_{\mathrm{abs}}$ used in the modBB fits) according to $\tau_\lambda= \int\kappa_{\mathrm{ext},\lambda}\rho\,\mathrm{d}s$. For shell edges corresponding to the features in the images (0''.05–0''.2, or $\sim$0.4–1.5$\times10^{15}$ m), the Zubko et al. [149] ACAR model and dust mass fit give an opacity at 6563 Å of $\tau_{\mathrm{H}\alpha}\sim560$. Silicate dust varieties tend to have lower $\kappa_{\mathrm{ext},\lambda}$ at a given wavelength – the amorphous forsterite model of Jäger et al. [61], for example, gives $\tau_{\mathrm{H}\alpha}\sim400$.

$\kappa_{\mathrm{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$ $\mu\mathrm{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.

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. [67]; Indebetouw et al. [58]; Zanardo et al. [146] and Matsuura et al. [81], and complements the molecular line studies of Matsuura et al. [82] and Abellán et al. [1] 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\sim500$. Dust extinction could be a factor in the appearance of the H$\alpha$ hole.

• We see an anti-correlation between the CO $J=6\to5$ emission and dust, and this anti-correlation is not seen when comparing dust and CO $J=2\to1$. 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.4M$_\odot$ for carbon or Si grains, or a total of $<0.7M_\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/L003059/1). Research at Garching was supported by the ERC through grant ERC-AdG No. 341157-COCO2CASA, and by the Deutsche Forschungsgemeinschaft through Sonderforschungsbereich SFB 1258 “Neutrinos and Dark Matter in Astro- and Particle Physics” (NDM) and the Excellence Cluster “ORIGINS: From the Origin of the Universe to the First Building Blocks of Life” (EXC 2094; http://www.origins-cluster.de). 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 cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

Facility: Atacama Large Millimeter Array

Software: astropy [6], CASA [88], emcee [36], ipython [104], lmfit [95], matplotlib [57], multicolorfits [22], numpy [131], RADEX [130], SAOImage DS9 [66], scipy [65], WWCC [5].

Molecular Line Channel Maps

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. [3] 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. [3]). They derive a ring center position of $\alpha=5^{\mathrm h}35^{\mathrm m}27^{\mathrm s}.9875$, $\delta=-69^\circ16'11''.00107$ (ICRF J2015.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^{\mathrm h}35^{\mathrm m}27^{\mathrm s}.998$, $\delta=-69^\circ16'11''.00107$ (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$\alpha$ emission.

Simulating Potential Flux Loss

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.04^{\prime\prime}$ 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.