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

Phil Cigan School of Physics and Astronomy, Cardiff University, Cardiff CF24
3AA, UK<EMAIL_ADDRESS>Mikako Matsuura School of Physics and Astronomy,
Cardiff University, Cardiff CF24 3AA, UK  Haley L. Gomez School of Physics and
Astronomy, Cardiff University, Cardiff CF24 3AA, UK  Remy Indebetouw
Department of Astronomy, University of Virginia, PO Box 400325,
Charlottesville, VA 22904, USA  Fran Abellán Departamento de Astronomía y
Astrofísica, Universidad de Valencia, C/Dr. Moliner 50, E-46100 Burjassot,
Spain  Michael Gabler  Max-Planck-Institut für Astrophysik, Karl-Schwarzchild-
Straße 1, D-85748 Garching, Germany  Anita Richards JBCA, School of Physics
and Astronomy, University of Manchester, Manchester M13 9PL, UK  Dennis Alp
Department of Physics, KTH Royal Institute of Technology, The Oskar Klein
Centre, AlbaNova, SE-106 91 Stockholm, Sweden  Tim Davis School of Physics and
Astronomy, Cardiff University, Cardiff CF24 3AA, UK  Hans-Thomas Janka Max-
Planck-Institut für Astrophysik, Karl-Schwarzchild-Straße 1, D-85748 Garching,
Germany  Jason Spyromilio European Southern Observatory, Karl-Schwarzschild-
Straße 2, 85748 Garching, Germany  M. J. Barlow Department of Physics and
Astronomy, University College London, Gower Street, London, WC1E 6BT, UK
David Burrows Department of Astronomy and Astrophysics, Pennsylvania State
University, University Park, PA 16802, USA  Eli Dwek Observational Cosmology
Lab., Code 665, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
Claes Fransson Department of Astronomy, Stockholm University, The Oskar Klein
Centre, AlbaNova, SE-106 91 Stockholm, Sweden  Bryan Gaensler Dunlap Institute
for Astronomy and Astrophysics, 50 St. George Street, Toronto, ON M5S 3H4,
Canada  Josefin Larsson Department of Physics, KTH Royal Institute of
Technology, The Oskar Klein Centre, AlbaNova, SE-106 91 Stockholm, Sweden  P.
Bouchet  DRF/IRFU/DAp, CEA-Saclay, F-91191 Gif-sur-Yvette, France  NRS/AIM,
Université Paris Diderot, F-9119, Gif-sur-Yvette, France  Peter Lundqvist
Department of Astronomy, Stockholm University, The Oskar Klein Centre,
AlbaNova, SE-106 91 Stockholm, Sweden  J. M. Marcaide  Departamento de
Astronomía y Astrofísica, Universidad de Valencia, C/Dr. Moliner 50, E-46100
Burjassot, Spain  C.-Y. Ng Department of Physics, The University of Hong Kong,
Pokfulam Road, Hong Kong  Sangwook Park Department of Physics, University of
Texas at Arlington, Arlington, TX 76019, USA  Pat Roche  Astrophysics,
Department of Physics, University of Oxford, DWB, Keble Road, Oxford OX1 3RH,
UK  Jacco Th. van Loon Lennard-Jones Laboratories, Keele University, ST5 5BG,
UK  J. C. Wheeler Department of Astronomy, The University of Texas at Austin,
1 University Station C1400, Austin, TX 78712-0259, USA  Giovanna Zanardo
International Centre for Radio Astronomy Research, M468, University of Western
Australia, Crawley, WA 6009, Australia

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

###### Abstract

We present high angular resolution ($\sim$80 mas) ALMA continuum images of the
SN 1987A system, together with CO $J$=2 $\\!\rightarrow\\!$ 1, $J$=6
$\\!\rightarrow\\!$ 5, and SiO $J$=5 $\\!\rightarrow\\!$ 4 to $J$=7
$\\!\rightarrow\\!$ 6 images, which clearly resolve the ejecta (dust continuum
and molecules) and ring (synchrotron continuum) components. Dust in the ejecta
is asymmetric and clumpy, and overall the dust fills the spatial void seen in
H$\alpha$ images, filling that region with material from heavier elements. The
dust clumps generally fill the space where CO $J$=6 $\\!\rightarrow\\!$ 5 is
fainter, tentatively indicating that these dust clumps and CO are locationally
and chemically linked. In these regions, carbonaceous dust grains might have
formed after dissociation of CO. The dust grains would have cooled by
radiation, and subsequent collisions of grains with gas would also cool the
gas, suppressing the CO $J$=6 $\\!\rightarrow\\!$ 5 intensity. The data show a
dust peak spatially coincident with the molecular hole seen in previous ALMA
CO $J$=2 $\\!\rightarrow\\!$ 1 and SiO $J$=5 $\\!\rightarrow\\!$ 4 images.
That dust peak, combined with CO and SiO line spectra, suggests that the dust
and gas could be at higher temperatures than the surrounding material, though
higher density cannot be totally excluded. One of the possibilities is that a
compact source provides additional heat at that location. Fits to the far-
infrared–millimeter spectral energy distribution give ejecta dust temperatures
of 18–23K. We revise the ejecta dust mass to
$\mathrm{M_{dust}}=0.2-0.4$$M_{\odot}$ for carbon or silicate grains, or a
maximum of $<0.7$$M_{\odot}$ for a mixture of grain species, using the
predicted nucleosynthesis yields as an upper limit.

††journal: The Astrophysical Journal††software:  astropy (Astropy
Collaboration et al., 2018), Casa (McMullin et al., 2007), emcee (Foreman-
Mackey et al., 2013), ipython (Pérez & Granger, 2007), lmfit (Newville et al.,
2016), matplotlib (Hunter, 2007), multicolorfits (Cigan, 2019), numpy (van der
Walt et al., 2011), radex (van der Tak et al., 2007), SAOImage DS9 (Joye &
Mandel, 2003), scipy (Jones et al., 2001–), WWCC (Arshakian & Ossenkopf,
2016).

##  1 Introduction

Multiwavelength studies of Supernova 1987A (SN 1987A), located at a distance
of $51.4\pm 1.2$ kpc in the Large Magellanic Cloud (Panagia, 1999), 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, 2016), 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$\farcs$6 (0.3 pc) (e.g. Burrows et al., 2000; Sonneborn et al., 1998; Ng et
al., 2013), as well as thermal dust emission due to shock heating of pre-
existing dust formed during the red supergiant phase (Bouchet et al., 2006;
Dwek et al., 2010). The ejecta have a complex morphology. The H$\alpha$
emission, originating from warm gas irradiated by X-rays from the ring
(Larsson et al., 2011; Fransson et al., 2013), exhibits an elongated north-
south structure and a ‘hole’ in the center. Along with hydrogen lines from the
ejecta, near-infrared (NIR) emission from warm ($\sim$2000 K) CO and mid-
infrared (MIR) emission from SiO in the SN ejecta were detected early (as
early as 112 days) after the explosion (e.g., Spyromilio et al., 1988; Roche
et al., 1991). After day 9,000, cold expanding CO, SiO and HCO+ molecules were
detected in the submillimeter (submm) part of the spectrum (Kamenetzky et al.,
2013; Matsuura et al., 2017), highlighting that a significant part of the
ejecta is cold (13–132 K). Interestingly, the inner ejecta of SN 1987A have
not yet mixed with the circumstellar medium (CSM) or interstellar medium (ISM)
and the majority has not yet passed through the reverse shock (France et al.,
2010; Frank et al., 2016). Thus this young SNR is an ideal source for studying
the footprints of the gas dynamics since the very early days of the SN, as the
gas has been assumed to be free-expanding since its explosion (McCray, 1993).
Indeed, recent high angular resolution emission line images of SiO and CO
(Abellán et al., 2017) from the Atacama Large Millimeter/submillimeter Array
(ALMA) have been used to compare the distribution of the molecular gas ejecta
with the predictions from models of the gas dynamics after the SN
(Wongwathanarat et al., 2015, Gabler et al., in preparation).

The progenitor of SN 1987A, Sanduleak –69∘ 202, was a blue supergiant (West et
al., 1987; White & Malin, 1987; Gilmozzi et al., 1987; Kirshner et al., 1987),
thought to have had a zero-age main sequence mass of $\sim 19\,\rm M_{\odot}$
(Woosley et al., 1997; Hashimoto, Nomoto & Shigeyama, 1989), with a mass of
$\sim 14\,\rm M_{\odot}$ at the time of the explosion (Woosley, 1988; Smartt
et al., 2009; Sukhbold et al., 2016). From its mass, the expectation is that a
neutron star should have formed at the time of explosion. Despite prompt
neutrino emission observed at the burst (Hirata et al., 1987) indicating the
formation of a neutron star (Burrows, 1988; Sukhbold et al., 2016), the search
for a compact object associated with SN 1987A has been difficult:
observational searches have proven unfruitful (e.g. Manchester, 2007; Alp et
al., 2018a; Esposito et al., 2018; Zhang et al., 2018). The possible detection
of radio polarization towards the ejecta (Zanardo et al., 2018) hints at the
presence of magnetized shocks, potentially due to a compact object. Alp et al.
(2018a) proposed that a thermally-emitting neutron star could be dust-
obscured, and that this may be detectable as a point source in far-infrared
(FIR) or 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., Morgan & Edmunds, 2003; Matsuura et al., 2009;
Gall et al., 2011; Gomez, 2013; Dwek et al., 2014; Rowlands et al., 2014b;
Michałowski et al., 2015; Watson et al., 2015; Lakićević et al., 2015; Temim
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.,
Morgan et al., 2003; Gomez et al., 2012a; De Looze et al., 2017; Chawner et
al., 2018). 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., Gall et al., 2014) and historical remnants.

Thermal emission from small quantities ($10^{-4}\,\rm M_{\odot}$) of dust was
detected in the early days after the SN explosion (day$\sim$300–600) using MIR
observations (Danziger et al., 1989; Lucy et al., 1989; Bouchet et al., 1991;
Roche et al., 1993; Wooden et al., 1993). More surprisingly, the Herschel
Space Observatory (Pilbratt et al., 2010, hereafter Herschel) revealed a large
amount of cold dust ($\sim 0.5\rm M_{\odot}$) at the location of the remnant
(Matsuura et al., 2011, 2015). ALMA resolved the emission from dust in SN
1987A on scales of $0.3^{\prime\prime}$ and confirmed that the $\sim\rm
0.5\,M_{\odot}$ of cold (20 K) dust discovered with Herschel originates from
the inner ejecta region (Indebetouw et al., 2014; Matsuura et al., 2015). Dwek
& Arendt (2015) and Wesson et al. (2015) re-visited the dust emission at early
times ($<$1200 days). Dwek & Arendt (2015) find that a large mass of dust can
be present early on (0.4$M_{\odot}$ at $\sim$615 days) with a model of
silicates and amorphous carbon. Wesson et al. (2015) conclude instead, from
comparing radiative transfer models to the optical–IR 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
(2009, 2010) and Sarangi & Cherchneff (2013, 2015) predict that for abundance
ratios $\rm 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, 2011), and/or (to a lesser
extent) X-rays from the ring, however, the models of Sarangi & Cherchneff
(2013) and Sarangi & Cherchneff (2015) indicate the dissociation of CO is
insignificant. In contrast, in order to explain the FIR dust emission, a
substantial fraction of the dust grains must be composed of amorphous carbon
(amC, Matsuura et al. 2015), as the emissivity of amC grains is higher than
that of silicates in general, thus leaving an unresolved tension between
observations and theory. We note that a model which explains the FIR emission
and requires only a small amount of amC grains, with the majority of mass in
silicates, was proposed by Dwek & Arendt (2015). There they fit the FIR SED
with amC-silicate composite grains assuming a “continuous distribution of
ellipsoids” (CDE) model and found a reduced dust mass, though the majority of
the reduction in mass in this case arises from the inclusion of dust grains
with long axial ratios (so-called needles), which allows the CDE model to
surpass the FIR emissivity of amorphous carbon. No evidence of the silicate
signature was found in the warm dust emission in the first two years after the
explosion, suggesting that small silicate grains were not the first
condensates (Roche et al., 1993).

In this work, we present high angular resolution ALMA (Cycle 2) dust images
for SN 1987A, where we resolve dust clumps on scales of $\sim$80 mas. Here, we
revisit the dust mass and grain composition using the ALMA photometry. We
discuss the implications of our results for the gas phase chemistry leading to
dust formation, and find evidence for warmer gas at the center of the inner
ejecta hinting at the possible indirect detection of a compact source.

##  2 Data

###  2.1 Observations and Reduction

Table 1:  Observations

| Sub- | Frequency | Baselines | Angular Scales | Observation | SN | Bandpass | Phase | Check | Time
---|---|---|---|---|---|---|---|---|---|---
| band | Range (GHz) | (m) | (arcsec) | Date | Day | Calibrator | Calibrator | Source | (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

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

Figure 1:  Spectra and integrated continuum flux densities of the ejecta and
ring for Bands 6, 7, and 9. The molecular line emission model is taken from
Matsuura et al. (2017). Vertical grey bands indicate the portions of the
spectrum deemed to be line-free and therefore used in creating the continuum
images at 307, 315, and 679 GHz (the 315 GHZ image is shown in Fig. 2). Data
points indicate the flux densities in that band for the ring (green) and the
ejecta (yellow).  Table 2:  Imaging Properties

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

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

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$\mu$m) and 9 (450$\mu$m) integrations utilized between 34 and 36
antennae with baselines spanning 15m to 2.3km. See Table 1 for a summary of
the observations.

Each data set was reduced separately with Common Astronomy Software
Applications package (Casa111http://casa.nrao.edu/, McMullin et al., 2007),
version 4.5.1. Once calibrated, the tclean algorithm was used to deconvolve
and image the data.

The check source (reference quasar with precisely known position) and phase
calibrator coordinates, determined with imfit in Casa, were offset by no more
than 0.4 mas from the catalog values. Other measures of astrometric quality
for our observing configurations include the ALMA baseline measurement
accuracy (2 mas), noise-limited signal error $\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 $\sim$35%. Bands 7A and 7D also suffered from poor
weather in the original June 2015 observations, with $\sim$1.45 mm
precipitable water vapor (PWV), and were therefore repeated in September 2015.
These are denoted as 7A2 and 7D2. Despite the poorer quality of the June data,
combining them with the September data results in higher S/N images.

Self-calibration, a common technique for high S/N data where calibrating the
data against itself in successive deconvolution cycles can often result in
improved dynamic range, was determined to have a negligible impact on the
images. Final images were cleaned with natural weighting applied to the
baselines to optimize sensitivity per beam. The imaging parameters, including
resolution and sensitivity, are given in Table 2.

###  2.2 Defining Continuum Wavelength Ranges

The wavelengths covered by these observations include many spectral lines from
molecular species – primarily CO and SiO, with contributions from various SO
lines and potentially others. The $\pm\sim$1000 km s-1 expansion velocity of
the ejecta means that the linewidths span a substantial fraction of the
observed bands. The continuum bands selected relative to the modelled
molecular line emission are shown in Fig. 1, using the ALMA spectra and the
emission line model of CO, SiO, SO, and SO2 from Matsuura et al. (2017). Only
windows that were free from molecular line emission (shown by the grey
vertical bands) were used to make continuum images, centered at roughly 307,
315, and 679 GHz. The 315 GHz continuum image is shown in Fig. 2. We also
utilize here the Cycle 2 Band–6 imaging data presented by Matsuura et al.
(2017), to provide continuum information below 300 GHz. The Band–6 images were
restored to a common circular beam with full-width at half-maximum (FWHM) of
0$\farcs$30.

###  2.3 Molecular Line Data

In this section we present the molecular line data observed in the same blocks
as the continuum discussed above: CO $J$=6 $\\!\rightarrow\\!$ 5 with rest
frequency 691.47 GHz and SiO $J$=7 $\\!\rightarrow\\!$ 6 at
$\nu_{\mathrm{rest}}$ = 303.93 GHz. The 345.80 GHz CO $J$=3
$\\!\rightarrow\\!$ 2 and 347.33 GHz SiO $J$=8 $\\!\rightarrow\\!$ 7 lines
were also covered in these observing blocks, but as they are heavily blended
we do not consider them in the present work.

The SiO $J$=7 $\\!\rightarrow\\!$ 6 and CO $J$=6 $\\!\rightarrow\\!$ 5 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 $\\!\rightarrow\\!$ 5 was imaged with natural weighting
to maximize sensitivity per beam. SiO $J$=7 $\\!\rightarrow\\!$ 6 was imaged
with robust$=-1$ in order to better spatially resolve the central features of
interest. A comparison of the integrated (Moment-0) maps of the CO and SiO
lines is given in Fig. 3.

In addition to the molecular line data described above, we also utilize the CO
$J$=2 $\\!\rightarrow\\!$ 1, SiO $J$=5 $\\!\rightarrow\\!$ 4, and SiO $J$=6
$\\!\rightarrow\\!$ 5222The images for the SiO $J$=6 $\\!\rightarrow\\!$ 5
transition were described but not shown in Abellán et al. (2017). data as
described in Abellán et al. (2017). Although both sets were taken in Cycle 2,
the molecular line data presented by Abellán et al. (2017) have higher signal
to noise as they were combined with Cycle 3 data, and due to CO $J$=6
$\\!\rightarrow\\!$ 5 being in a band with poorer atmospheric transmission
than CO $J$=2 $\\!\rightarrow\\!$ 1. The additional Cycle 3 data also give
their CO $J$=2 $\\!\rightarrow\\!$ 1, SiO $J$=5 $\\!\rightarrow\\!$ 4, and SiO
$J$=6 $\\!\rightarrow\\!$ 5 maps finer spatial resolution than the
observations presented in the current work, with FWHM 0$\farcs$04–0$\farcs$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., 2008).

##  3 Description of Images

Figure 2:  ALMA 315 GHz (with beam) and 2014 HST 625W band image (Fransson et
al., 2015), which includes H$\alpha$. The yellow contours display 315 GHz
emission at 0.2 mJy/beam. The 315 GHz continuum in the inner ejecta originates
from thermal dust emission, while in the ring it is due to synchrotron
emission. The 18 mas uncertainty on the relative alignment due to Band–7
astrometric error (12 mas) and HST image registration based on fitting the
ring (6 mas) is of order 1 pixel in these images.  Figure 3:  Multi-wavelength
view of the ejecta in SN 1987A. The cyan contours represent the 679 GHz dust
emission at 3$\sigma$ and 5$\sigma$. Beam sizes for individual maps are
denoted by the green ellipses. The small cross denotes the system center as
defined in Appendix B. The bottom right panel is a 3-color image of CO $J$=2
$\\!\rightarrow\\!$ 1 in red, SiO $J$=5 $\\!\rightarrow\\!$ 4 in green, and
SiO $J$=7 $\\!\rightarrow\\!$ 6 in blue, and highlights how varied the
spectrally-integrated emission is between the various line transitions. The
brightest areas are generally distinct patches of primary color instead of
blended, demonstrating that the CO and SiO peaks are not cospatial, and none
matches the distribution of the 679 GHz dust. Some of the CO falls in the
H$\alpha$ hole (the lower left), but the majority of the CO peaks on the
little H$\alpha$ ‘wing’ to the right of the hole. The small $5\sigma$ cyan
contour just northeast of the center of the remnant is the so-called “blob”.
Figure 4:  The spatial anti-correlation of dust and CO $J$=6
$\\!\rightarrow\\!$ 5. In this overlay, CO $J$=6 $\\!\rightarrow\\!$ 5 is in
blue and 679 GHz is in red, showing their relative spatial distributions. The
off-white cross denotes the system center position as defined in Appendix B.
The gold line highlights the extent of the major features in the dust, while
the teal line demarcates the major CO $J$=6 $\\!\rightarrow\\!$ 5 features.
The contrast of each component image was set independently to emphasize its
major features – visible blue and red features roughly correspond to S/N$>$3
for the dust and S/N$>$5 for CO $J$=6 $\\!\rightarrow\\!$ 5, respectively. The
dust and CO $J$=6 $\\!\rightarrow\\!$ 5 emission exhibit a notable anti-
correlation.

The SN 1987A ring and ejecta continuum image at 315 GHz is shown in comparison
to H$\alpha$ in Fig. 2. These images have been aligned following a technique
in Alp et al. (2018a) where the ring emission is used to derive a reference
center, though here we take a simpler approach (see Appendix B for details).
Our derived ring+ejecta system center used in this work is
$\alpha$=5h35m27$\fs$998, $\delta$=–69∘16′11$\farcs$107 (ICRS), $\pm$18 mas
(Fig. B.1). At 315 GHz, the ring is clumpy and the brightness contrast in the
east and west components of the ring is different to that observed in the
H$\alpha$ ring emission. The brighter emission observed in the NE and SW
regions of the ring in the radio is similar to that seen in hard X-rays
(Helder et al., 2013; Frank et al., 2016). The ejecta are located at the
center of the image inside of the ring structure. The ring emission at 315 GHz
is attributed to synchrotron (see § 4.3), and the inner region is thermal dust
emission from the SN ejecta (Indebetouw et al., 2014).

Fig. 3 shows an enlarged view of the ejecta images of dust continuum and
lines. The majority of the submm ejecta continuum is distributed in a roughly
symmetrical ellipsoid, with fainter asymmetrical emission protruding west and
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 Fig. 2, where the 3$\sigma$ contours highlighting the major
continuum features are overlaid onto the continuum and H$\alpha$ images. The
alignment accuracy is $\sim$1 pixel in the images, given the astrometric
uncertainties discussed in § 2.1 (12 mas for Band–7 continuum) and Appendix B
(6 mas for registration of the Hubble Space Telescope (HST) image to the 315
GHz ALMA image).

The 679 GHz image provides the highest resolution view of the continuum (top
left panel of Fig. 3). This figure shows that the dust is asymmetrically
distributed and is composed of several clumps, with the brightest feature
(hereafter the “blob”) just northeast of the center of the remnant. The beam
resolution provides a limit on the clump size – assuming a distance of 51.4
$\pm$ 1.2 kpc (Panagia, 1999), the Band–9 beam FWHM of
0$\farcs$08$\times$0$\farcs$06 corresponds to a physical scale of
0.020$\times$0.016 pc, or 4125$\times$3230 au. Nevertheless, the resolved 679
GHz image indicates that dust is not smoothly distributed across the ejecta,
and the locations of dust clumps are not identical to clumps in the CO or SiO.
The S/N in the 679 GHz image is moderate – the outer cyan contours in Fig. 3
and the dust emission (in red) in Fig. 4 have pixel S/N$>$3, and the
surrounding ejecta area has pixel S/N values of $\sim$2 in the 679 GHz image.
The area between the ejecta and the ring – outside of the outermost ejecta
contour in Fig. 2 – is consistent with noise.

The molecular images provide a probe of different conditions in the ejecta,
where lower transitions probe lower temperature gas (if optically thin, see
Section 5). One prominent feature is the central hole seen in the CO $J$=2
$\\!\rightarrow\\!$ 1 and SiO $J$=5 $\\!\rightarrow\\!$ 4 images (middle left
and lower left panels of Fig. 3). This was first reported by Abellán et al.
(2017) and was seen both in the integrated 2D spatial maps and the 3D data
cubes. Although the integrated SiO $J$=6 $\\!\rightarrow\\!$ 5 map (middle
panel of Fig. 3) does not show the hole clearly in the same manner as SiO
$J$=5 $\\!\rightarrow\\!$ 4 and CO $J$=2 $\\!\rightarrow\\!$ 1, the hole is
also visible in the central channels ($v=0-300$ km s-1) of the velocity map
(Fig. A.2). Because of the additional $-$600–0 km s-1 components located
within the same line of sight as the hole (Fig. A.2) in the integrated maps,
the hole is not clear in the SiO $J$=6 $\\!\rightarrow\\!$ 5 map. The CO and
SiO molecular hole is just to the south of the ‘keyhole’ that is seen in
H$\alpha$ (Fig. 8 of Fransson et al. (2015); top right panel of our Fig. 3),
though the molecular hole appears to be slightly smaller in scale and located
on the southern edge of the hole in H$\alpha$ emission. The centers of the
holes are offset by $\sim$50 mas, or $\sim$4$\times$ the astrometric and
alignment errors.

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

CO $J$=6 $\\!\rightarrow\\!$ 5 has emission coincident with the CO $J$=2
$\\!\rightarrow\\!$ 1 hole, in that its channel maps (Fig. A.1) show emission
around the hole location, albeit at low S/N. However, the integrated spatial
distribution appears different from CO $J$=2 $\\!\rightarrow\\!$ 1. The
brightness peaks are distributed differently, and the hole is not visible in
the integrated CO $J$=6 $\\!\rightarrow\\!$ 5 map due to some emission at
those coordinates in the 600–900 km s-1 channels (the far side). The presence
of a molecular hole in SiO $J$=7 $\\!\rightarrow\\!$ 6 cannot be confirmed in
these data, as the systemic line center ($v_{\mathrm{LSRK}}\sim$300 km s-1)
falls at the edges of two sidebands observed separately, which were
concatenated during reduction, and suffers from roll-off at the edge of the
spectral window; the resulting S/N is poor in that channel. The other
molecular lines do not share this limitation as they fell well within the
sideband spectral windows. We do note a peak of SiO $J$=7 $\\!\rightarrow\\!$
6 emission, however – the brightest source of emission in the entire cube –
overlapping with the spatial location of the hole and the dust blob but offset
from the systemic velocity by $\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 Fig. 3) is co-located with
the molecular hole in the low transitions of CO and SiO, and slightly extends
to the north and east into the relative depression visible in the SiO $J$=5
$\\!\rightarrow\\!$ 4 channels near the systemic velocity. The brightest
points of dust emission tend to coincide with relative depressions in the CO
$J$=6 $\\!\rightarrow\\!$ 5 brightness, giving the appearance of an anti-
correlation between the main dust and CO $J$=6 $\\!\rightarrow\\!$ 5 features.
This is more clearly demonstrated in Fig. 4, where the dust (red) and CO $J$=6
$\\!\rightarrow\\!$ 5 (blue) images are overlaid. The individual images were
normalized independently to emphasize the main features of each, with the
visible colors shown roughly corresponding to areas of S/N$>$3\. The gold and
teal lines are guides highlighting the highest-S/N features of the dust and
CO, respectively, in order to compare peaks in the emission. While there is
some overlap in the faint features of the dust and CO, and the southern extent
of the dust peak starts to fade to a combined magenta, the dust and CO $J$=6
$\\!\rightarrow\\!$ 5 peaks do not generally overlap. Rather, the brightest
dust features are located in areas of relatively faint CO $J$=6
$\\!\rightarrow\\!$ 5 emission and vice-versa.

To test whether the apparent dust-CO anti-correlation is an artifact or result
of the data reduction or continuum 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 $\\!\rightarrow\\!$ 5 map are robust to variations in how the continuum
is determined and subtracted. The anti-correlation in CO $J$=6
$\\!\rightarrow\\!$ 5 is visible, even before continuum subtraction, in the CO
$J$=6 $\\!\rightarrow\\!$ 5 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\mathrm{cov}_{XY}/\sqrt{\mathrm{cov}_{XX}\cdot\mathrm{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 $\pm$ 0.08 – a moderate anti-correlation – using the accuracy
estimate from, e.g., Frick et al. (2001). Due to the relatively low S/N in
these images and therefore scatter in pixel-by-pixel comparisons, $r$ will
always be pulled closer to 0 and will not approach $\pm$1.

We tested the robustness of the correlation measure by investigating different
angular scales using the wavelet analysis described in Frick et al. (2001) and
Arshakian & Ossenkopf (2016). On scales of 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 ($\sim 5\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 $\\!\rightarrow\\!$ 5. The standard correlation measures agree, with
$r$=+0.04 for CO $J$=2 $\\!\rightarrow\\!$ 1 and $r$=+0.36 for SiO $J$=6
$\\!\rightarrow\\!$ 5.

The brightest dust feature is located one beam width northeast of the
secondary CO $J$=6 $\\!\rightarrow\\!$ 5 peak, and the brightest CO $J$=6
$\\!\rightarrow\\!$ 5 features curve around the main dust peak. The CO peaks
are obvious because the line emission is brighter than the continuum (c.f. the
integrated profile in Fig. 1), but there is some fainter (S/N$<$5) CO $J$=6
$\\!\rightarrow\\!$ 5 emission overlapping with some of the dust emission.
There is also low-level (S/N$<2$) dust emission that roughly spans the full
extent of the ejecta. The other molecular species further complicate the
picture, as noted earlier – the peak of SiO $J$=7 $\\!\rightarrow\\!$ 6 is
coincident with the dust peak and the SiO $J$=5 $\\!\rightarrow\\!$ 4 hole (as
seen in the bottom right panel of Fig. 3). At the southern edge of the hole,
some faint H$\alpha$ emission appears to be aligned with CO $J$=2
$\\!\rightarrow\\!$ 1 in projection, but their velocity ranges differ. The 3-D
view of H$\alpha$ (Larsson et al., 2016, 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 $\\!\rightarrow\\!$ 1 is between 0 and +200 km s-1.
While no velocity information is available for the dust continuum emission, it
is spatially offset from the nearby H$\alpha$ and CO $J$=2 $\\!\rightarrow\\!$
1 by $\sim$1 dust resolution element. That is, the CO $J$=2
$\\!\rightarrow\\!$ 1 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
$\\!\rightarrow\\!$ 5 and dust structures). The small peak/clump in the dust
emission revealed in the Band 9 image (the blob) overlaps with holes
previously observed in the lower line transitions of SiO and the CO molecular
ejecta, and is coincident with some emission observed in the SiO $J$=7
$\\!\rightarrow\\!$ 6 line.

##  4 The Spectral Energy Distribution of SN 1987A

The three physical mechanisms primarily responsible for emission in the
FIR–radio portion of the continuum are thermal IR greybody emission from dust
(from the ejecta region - Matsuura et al., 2011; Indebetouw et al., 2014),
non-thermal radio/mm synchrotron emission (from the ring - Manchester, 2007;
Potter et al., 2009; Zanardo et al., 2010; Lakićević et al., 2012; Zanardo et
al., 2013; Indebetouw et al., 2014) and a lesser contribution from free-free
mm/sub-mm bremsstrahlung emission from hot ionized material (see Zanardo et
al. 2014 for a full review of the different components). In this Section we
measure the photometry, analyze the emission from dust in the ejecta using the
ALMA data, investigate the properties derived using a variety of dust models
from the literature, and investigate the synchrotron emission in the ring.

###  4.1 Photometry

Table 3:  Continuum Photometry

$\nu_{c}$ | $\Delta\nu$ | Sν | Aperture Center | RMAJ${\dagger}$${\dagger}$The ring annulus radii are given as (Router, Rinner). | RMIN${\dagger}$${\dagger}$The ring annulus radii are given as (Router, Rinner). | P.A.
---|---|---|---|---|---|---
(GHz) | (GHz) | (mJy) | RA (deg) | DEC (deg) | (arcsec) | (arcsec) | (deg)
Ejecta
225.50 | 3.0 | 0.5 $\pm$ 0.2 ${}^{+0.3}_{-0.0}$ | 83.866586 | -69.269733 | 0.32 | 0.29 | 125
240.50 | 5.0 | 0.7 $\pm$ 0.2 ${}^{+0.4}_{-0.1}$ | 83.866600 | -69.269720 | 0.37 | 0.36 | 75
247.50 | 3.0 | 1.1 $\pm$ 0.4 ${}^{+0.6}_{-0.1}$ | 83.866639 | -69.269739 | 0.41 | 0.35 | 125
270.00 | 1.0 | 1.2 $\pm$ 0.6 ${}^{+0.7}_{-0.1}$ | 83.866652 | -69.269709 | 0.40 | 0.31 | 60
279.00 | 2.0 | 1.3 $\pm$ 0.4 ${}^{+0.8}_{-0.1}$ | 83.866662 | -69.269747 | 0.41 | 0.40 | 125
306.76 | 1.4 | 1.8 $\pm$ 0.6 ${}^{+1.1}_{-0.2}$ | 83.866667 | -69.269753 | 0.42 | 0.36 | 125
315.68 | 7.6 | 1.8 $\pm$ 0.6 ${}^{+1.1}_{-0.2}$ | 83.866667 | -69.269753 | 0.42 | 0.36 | 125
679.22 | 11.5 | 36.2 $\pm$ 7.2 ${}^{+7.2}_{-7.2}$ | 83.866728 | -69.269740 | 0.42 | 0.36 | 125
Ring
225.50 | 3.0 | 13.0 $\pm$ 0.5 ${}^{+7.8}_{-1.3}$ | 83.866585 | -69.269731 | 1.45, 0.35 | 1.35, 0.33 | 175
240.50 | 5.0 | 12.5 $\pm$ 0.5 ${}^{+7.5}_{-1.3}$ | 83.866600 | -69.269722 | 1.45, 0.42 | 1.35, 0.39 | 175
247.50 | 3.0 | 13.4 $\pm$ 0.8 ${}^{+8.1}_{-1.3}$ | 83.866630 | -69.269736 | 1.45, 0.43 | 1.35, 0.40 | 175
270.00 | 1.0 | 13.7 $\pm$ 1.1 ${}^{+8.2}_{-1.4}$ | 83.866592 | -69.269720 | 1.45, 0.45 | 1.35, 0.42 | 175
279.00 | 2.0 | 12.9 $\pm$ 0.8 ${}^{+7.7}_{-1.3}$ | 83.866661 | -69.269747 | 1.45, 0.44 | 1.35, 0.41 | 175
306.76 | 1.4 | 12.0 $\pm$ 1.5 ${}^{+7.2}_{-1.2}$ | 83.866666 | -69.269753 | 1.45, 0.45 | 1.35, 0.42 | 175
315.68 | 7.6 | 11.3 $\pm$ 1.7 ${}^{+6.8}_{-1.1}$ | 83.866667 | -69.269753 | 1.45, 0.45 | 1.35, 0.42 | 175
679.22 | 11.5 | 19.4 $\pm$ 19.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 $\pm$ 0.5 ${}^{+8.1}_{-1.3}$ | 83.866585 | -69.269731 | 1.45 | 1.35 | 175
240.50 | 5.0 | 13.2 $\pm$ 0.5 ${}^{+7.9}_{-1.3}$ | 83.866600 | -69.269722 | 1.45 | 1.35 | 175
247.50 | 3.0 | 14.6 $\pm$ 0.7 ${}^{+8.7}_{-1.5}$ | 83.866630 | -69.269736 | 1.45 | 1.35 | 175
270.00 | 1.0 | 15.4 $\pm$ 0.8 ${}^{+9.2}_{-1.5}$ | 83.866592 | -69.269720 | 1.45 | 1.35 | 175
279.00 | 2.0 | 14.2 $\pm$ 0.8 ${}^{+8.5}_{-1.4}$ | 83.866661 | -69.269747 | 1.45 | 1.35 | 175
306.76 | 1.4 | 13.9 $\pm$ 1.7 ${}^{+8.3}_{-1.4}$ | 83.866666 | -69.269753 | 1.45 | 1.35 | 175
315.68 | 7.6 | 13.3 $\pm$ 1.8 ${}^{+8.0}_{-1.3}$ | 83.866667 | -69.269753 | 1.45 | 1.35 | 175
679.22 | 11.5 | 55.5 $\pm$ 20.9 ${}^{+11.1}_{-11.1}$ | 83.866720 | -69.269737 | 1.45 | 1.35 | 175

Note. — Integrated flux densities of the ejecta, ring, and total system for
each continuum band with central frequency $\nu_{c}$ and bandwidth
$\Delta\nu_{c}$. Integrated flux densities are quoted as value $\pm$
measurement uncertainty $\pm$ systematic uncertainty. The systematic
uncertainty includes calibration uncertainties (10% in Bands 6&7 or 20% in
Band 9), and an additional 50% on the positive side in Bands 6&7 due to the
systematic offset from Cycle 0 and Cycle 6 levels. Apertures for the ejecta
and total system are ellipses, apertures for the ring are annuli. The center
coordinates are in ICRS.

Figure 5:  Continuum values from 225–679 GHz for the integrated ejecta
(yellow), ring (green), and total system (purple) from this work, for
observations taken an average of 10,402 days after the SN explosion. For
reference, ALMA Cycle 0 corresponds to day 9280, Cycle 2 is day 10402, and
Cycle 6 is day 11522. The higher angular resolution ALMA data confirm the
ejecta emission follows a thermal dust profile down to $\sim$200 GHz. The ring
emission, on the other hand, shows no evidence of submm dust emission but
instead is consistent with a synchrotron emission profile – ATCA observations
from Days 9280–9686 (Zanardo et al., 2014; Callingham et al., 2016) combined
with the Cycle 2 ALMA ring flux densities give a power-law index of
$\alpha=-0.70\pm 0.06$. The dark blue and cyan lines show SED fits for
amorphous carbon (ACAR sample, Zubko et al., 1996) and silicate (forsterite,
Jäger et al., 2003) dust emission models, respectively, demonstrating that
disparate models give reasonable fits to the data. The previous ALMA ejecta
measurements from Zanardo et al. (2014) are shown as red stars, and previous
ALMA ring measurements from Zanardo et al. (2014) are denoted as blue stars.
Preliminary Cycle 6 ALMA flux densities (Matsuura et al., in preparation) are
shown as the orange crosses and light purple hexagons for the ejecta and total
emission, respectively. The positive error bars for the Cycle 2 Band 6&7 data
include an additional 50% of the flux density values to reflect the observed
systematic offset from the Cycle 0 and Cycle 6 levels. The empty grey hexagons
are measures of the total system flux density in various parts of the spectrum
at day 9280: 1.4 GHz, 18 GHz, 44 GHz (Zanardo et al., 2013); 9 GHz (Ng et al.,
2013); 94 GHz (Lakićević et al., 2012); 102, 213, 345, & 672 GHz (Zanardo et
al., 2014). The empty black hexagons represent the ATCA total system flux
densities between 1–9 GHz at day 9686 (Callingham et al., 2016). The ring flux
density at day 9280 (blue stars) is larger than ALMA total emission (purple
diamonds) at day 10402 due to the ALMA band 7 systematic error discussed in
Section 4.1. The brown asterisks are the unresolved Herschel 70–500 $\mu$m
flux densities (Matsuura et al., 2015), and the crimson asterisks are the 2010
HERITAGE flux densities (Meixner et al., 2013; Matsuura et al., 2011). The
ring emission model S($\nu$,t) from Cendes et al. (2018) for day 10402 and
$\alpha=-0.70$ is shown by the dashed grey line.

The continuum bands are defined as those frequencies that are molecular line
free, as demonstrated in Fig. 1 and §2.2, with the chosen frequency windows
summarized in Table 3; these bands are different from the default ALMA wide
band continuum. The centers of the apertures used for deriving photometry are
the same as described in Section 3 (Appendix B), with elliptical apertures
selected to encompass the ejecta and the ring annulus with varying sizes in
each ALMA band in order to include only the relevant signal (Table 3). For the
315 GHz ejecta, this results in an elliptical aperture with semi-major and
semi-minor axes of $r_{\mathrm{MAJ}}\times r_{\mathrm{MIN}}$ =
0$\farcs$42$\times$0$\farcs$36 and major axis P.A. of 25∘ (N through E). For
the 315 GHz circumstellar ring, an elliptical annulus with inner
$r_{\mathrm{MAJ}}$ and $r_{\mathrm{MIN}}$ of 0$\farcs$45 $\times$ 0$\farcs$42,
outer $r_{\mathrm{MAJ}}$ and $r_{\mathrm{MIN}}$ of 1$\farcs$45 $\times$
1$\farcs$35, and 85∘ P.A. from N was chosen. The extents of the regions were
selected independently across the different bands to best match the features
in each image, but only vary slightly. For comparison with previous lower
spatial resolution observations, the total system emission is also calculated
by summing emission within an ellipse defined by the outer ring extent above.
This total system integration includes the contribution from the gap between
the ring and the ejecta.

The Cycle 2 images, aside from Band 9, exhibit a slight decrease in integrated
flux density for radii just beyond the outer edge of the ring, which could be
due to undersampled flux or errors in calibration or deconvolution. We take
the RMS of flux densities in background pixels from a large annulus beyond the
ring to estimate the level of these effects, and the resulting uncertainty
contribution is typically of order a few percent of the flux density.

As the reconstruction of the images may propagate systematic as well as random
noise in the background, we have used our images to make a series of
measurements using the same aperture as for the source. The distribution of
these measurements is roughly gaussian, and thus we adopt the RMS of this
distribution as our aperture error, $\sigma_{\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$\farcs$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 $\sim$15%, dominated by the random-
position aperture uncertainty). We include an additional 10% uncertainty for
systematic (calibration) error in Bands 6 & 7, and 20% in Band–9333ALMA Cycle
2 Technical Handbook,
https://almascience.eso.org/documents-and-tools/cycle-2/alma-technical-
handbook/. Finally, we considered the possibility that we may be missing
diffuse emission from cold dust within the SN structure due to over-resolving
an extended source. To address this issue, we simulated observations for
multiple synthetic sources resembling SN 1987A but with varying extended
ellipse components, and found that this effect is at a level below the ALMA
systematic uncertainties. The reader is referred to Appendix C for more
details.

The ejecta, ring and total system flux densities and uncertainties are shown
in Fig. 5 as gold circles, green rings, and purple diamonds, respectively, and
their values are listed in Table 3. Previous measurements are also shown in
Fig. 5 for reference. Preliminary Cycle 6 flux densities (Matsuura et al., in
preparation) from 11,522 days after the explosion, are included here. The
ejecta flux density is 1.6 mJy at 252.4 GHz, and 1.7 mJy at 254.3 GHz. The
total system flux densities at these frequencies are 17.9 mJy and 18.1 mJy,
respectively. The uncertainty on each of these flux density measurements is
estimated as 0.4 mJy.

We note that our Cycle 2 total flux densities in Bands 6 and 7 are
systematically lower than the ALMA Cycle 0 and 6 flux densities, though they
agree within the error bars. They are typically around 50% lower than the
equivalent levels from Cycles 0 and 6, therefore we include an additional 50%
to their positive systematic uncertainties. Potential causes, as noted in §
2.1, include decorrelation from poor weather or a mis-scaling of the flux
calibrator. The integrated 350 and 360 GHz ejecta flux densities are
particularly low, either due to inherently low flux at this epoch, or due to
data quality. As mentioned in § 2.1, weather affected the phase stability of
observations between 346–362 GHz, which can result in decorrelation and
therefore reduced flux recovery. As these measurements are less reliable, they
have been omitted from the remainder of this study. We note that the
systematic offset will not have affected the analysis of the resolved dust
distribution discussed in Section 3, since the offset is not seen in the Band
9 data where the dust peak (the blob) was identified.

The literature values for the total SN 1987A flux densities at various
wavelengths are shown as grey hexagons, and represent the overall spectral
energy distribution (SED) of the system. The total emission at 1.4 GHz, 18
GHz, 44 GHz (Zanardo et al., 2013), 9 GHz (Ng et al., 2013), 94 GHz (Lakićević
et al., 2012), and 102 GHz (Zanardo et al., 2014) is dominated by synchrotron
emission from the ring. The total emission at 213, 345, & 672 GHz (Zanardo et
al., 2014), on the other hand, gradually consists of a higher and higher
fraction of thermal emission until that is dominant in the submm and the FIR.
As the synchrotron brightness increases in time (Staveley-Smith et al., 2014;
Cendes et al., 2018), these literature flux densities were scaled to their
levels at day 9280 by Zanardo et al. (2014) to match the average epoch of the
ALMA cycle 0 observations. The details of the ejecta and ring portions of the
SED will be discussed in turn in the following two sections.

###  4.2 Modified Blackbody Fits to the Ejecta Dust Emission

####  4.2.1 Description of the modified blackbody fits

Table 4:  Modified Blackbody Fits

Dust | Reference | Grain Density | Mdust | T | $\kappa_{850}$ | Good Fit
---|---|---|---|---|---|---
|  | (g cm-3) | (Mdust ) | (K) | (m2 kg-1) |
amC (ACH2 sample), Mie 0.1$\mu$m | Zubko et al. (1996) | 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. (1996) | 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. (1996) | 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 (1991) | 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. (1998) | 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 (1984) | 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 (1993) | 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. (2003) | 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. (2003) | 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. (1994) | 3.2 | 3.46 ${}^{+0.24}_{-0.22}$ | 17.7 ${}^{+0.2}_{-0.2}$ | 0.034 | Y
Silicate/Carbon Composite CDE – fC=0.18 | Dwek & Arendt (2015) | 2.95 | 0.38 ${}^{+0.02}_{-0.02}$ | 21.1 ${}^{+0.2}_{-0.2}$ | 0.270 | N
Silicate – LMC Average | Weingartner & Draine (2001) | $\cdots$ | 2.49 ${}^{+0.15}_{-0.14}$ | 18.0 ${}^{+0.2}_{-0.2}$ | 0.047 | Y
Silicate – 30K Average | Demyk et al. (2017) | $\cdots$ | 0.27 ${}^{+0.02}_{-0.02}$ | 18.5 ${}^{+0.2}_{-0.2}$ | 0.265 | Y
Silicate – Composite Aggregate | Semenov et al. (2003) | $\cdots$ | 35.06 ${}^{+2.18}_{-2.18}$ | 21.4 ${}^{+0.2}_{-0.2}$ | 0.003 | Y
Silicate – Porous Multilayer Spheres | Semenov et al. (2003) | $\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. (2011) | $\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. (2011) | $\cdots$ | 0.60 ${}^{+0.04}_{-0.04}$ | 18.3 ${}^{+0.2}_{-0.2}$ | 0.184 | Y
Silicate – Naked Grains, n${}_{H}=10^{5}$ | Ossenkopf & Henning (1994) | $\cdots$ | 0.64 ${}^{+0.04}_{-0.04}$ | 20.9 ${}^{+0.2}_{-0.2}$ | 0.163 | N
Silicate – Thin Ice Mantles, n${}_{H}=10^{5}$ | Ossenkopf & Henning (1994) | $\cdots$ | 0.74 ${}^{+0.04}_{-0.04}$ | 19.1 ${}^{+0.2}_{-0.2}$ | 0.142 | Y
SiC, Mie 0.1$\mu$m | Pegourie (1988) | 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 (1996) | 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. (1995) | 5.7 | 0.28 ${}^{+0.01}_{-0.02}$ | 28.2 ${}^{+0.4}_{-0.3}$ | 0.259 | N
SiO2, Mie 0.1$\mu$m | Henning & Mutschke (1997) | 2.196 | 4.41 ${}^{+0.39}_{-0.33}$ | 17.5 ${}^{+0.2}_{-0.2}$ | 0.022 | N
TiO2, Mie 0.1$\mu$m | Posch et al. (2003) | 3.78 | 81.77 ${}^{+5.10}_{-4.49}$ | 17.7 ${}^{+0.2}_{-0.2}$ | 0.001 | Y
Al2O3 “Compact” sample, Mie 0.1$\mu$m | Begemann et al. (1997) | 3.2 | 0.90 ${}^{+0.06}_{-0.06}$ | 19.0 ${}^{+0.2}_{-0.2}$ | 0.112 | Y
NaAlSi2O6, Mie 0.1$\mu$m | Mutschke et al. (1998) | 2.4 | 0.10 ${}^{+0.01}_{-0.01}$ | 23.1 ${}^{+0.3}_{-0.3}$ | 0.982 | Y
MgAl2O4, Mie 0.1$\mu$m | Fabian et al. (2001) | 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 (1996) | 7.87 | 3.97 ${}^{+0.28}_{-0.23}$ | 19.9 ${}^{+0.2}_{-0.2}$ | 0.025 | Y

Note. — Mass and temperature fits for greybodies for a selection of 28 of the
dust models discussed in the text. amC: amorphous carbon. Fit values and
uncertainties are from bootstrap resampling of the data within their error
bars, and are determined from the 50th, 84th, and 16th percentiles of the
distributions of fits from 1000 samplings of the observed flux densities. For
$\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. (2004) (and Rouleau & Martin (1991) for their AC1 sample), graphite and
SiC from Laor & Draine (1993), PAHs from (Li & Draine, 2001), stoichiometric
varieties of olivines and pyroxenes from Henning & Stognienko (1996),
silicate/carbon composite CDE with carbon volume filling factor fC=18% from
Dwek & Arendt (2015), jadeite from Mutschke et al. (1998), and in general from
the Jena optical constants database444https://www.astro.uni-
jena.de/Laboratory/OCDB/index.html (Henning et al., 1999). $\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$.

Fig. 5 displays the mm to FIR SED. In this figure, the brown asterisks show
the FIR flux densities measured by Herschel for the total SN 1987A
(unresolved) system, and the gold circles show the mm flux densities from this
work measured with ALMA for the resolved ejecta. The shape of the SED shows
that the ejecta emission arises from thermal (dust) radiation, all the way
into the mm, confirming the results of Zanardo et al. (2014) and Matsuura et
al. (2015).

The next step is to fit the thermal dust emission using dust models. In order
to cover the peak of the thermal emission, we use the Herschel flux densities
from Matsuura et al. (2015) in our model fits since they are measuring the
emission from the ejecta dust, albeit unresolved. Two Herschel flux densities
are treated as upper limits: 70 $\mu$m, as it is possibly contaminated by warm
ring dust (Matsuura et al., 2019) and/or [O I] 63 $\mu$m emission; and 500
$\mu$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 $\sim$1300 days before the ALMA cycle 2 data, and the FIR emission could
potentially vary over time. The heating source of the ejecta was suggested to
be primarily from 44Ti decay, which has an estimated lifetime of 85 years
(Ahmad et al., 2006; Jerkstrand et al., 2011, also see later Matsuura et al.
(2011)). The predicted decrease in this decay energy between the 2012 Herschel
and 2015 ALMA observations is 4.2%. Assuming the FIR luminosity decreased by
this amount between the 2012 and 2015 epochs, the reduction in the temperature
of a 20K blackbody would be $\sim$0.2K, translating into individual Herschel
flux density decrements of 2–5%. This is several times smaller than the
uncertainties on the PACS and SPIRE flux densities. Therefore, if the ejecta
heating is dominated by 44Ti decay, the use of the Herschel flux densities
with the latest ALMA measurements is valid. An alternative additional heating
source will be discussed in Section 6.3.

We tested the robustness of the SED results using three common parameter
estimation techniques: using a maximum likelihood estimation (MLE) with
uncertainties determined by bootstrap resampling; bayesian estimation with
Markov Chain Monte Carlo (MCMC) posterior distributions, using the emcee
package (Foreman-Mackey et al., 2013); and finally by checking with ordinary
least squares (OLS) regression. The OLS, MLE, and MCMC routines all yield
consistent fits for a given dust emission profile. In order to take into
account the systematic offset in our Band 6&7 flux densities from other cycles
(see § 4.1), we determine our best fits and uncertainties from resampling of
the flux densities within their error bars. The best fit and uncertainties are
taken to be the 50th, 16th, and 84th percentiles of the distributions from
1000 samplings.

The modified blackbody (modBB) function we use follows the form
$S_{\nu}(\lambda)=M_{dust}\kappa_{\rm{abs}}(\lambda)B_{\nu}(\lambda)/d^{2}$,
where $S_{\nu}(\lambda)$ is flux density, $\kappa_{\rm{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_{\rm{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_{\rm{abs}}(\lambda)$ is defined as $\kappa=3/4\,\,Q\rho a$, where $Q$
is the absorption efficiency of the dust, which can be calculated from optical
constants with Mie theory. $\kappa_{\rm{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_{\rm{abs}}(\lambda)$ at
wavelength $\lambda_{0}$ and $\beta$ is the power law emissivity index of
$\kappa_{\rm{abs}}(\lambda)$).

Fitting a dust greybody to the FIR-mm SED with the power law approximation to
$\kappa_{\rm{abs}}$ gives fit parameters of $\beta=2.05^{+0.11}_{-0.10}$,
${M_{\mathrm{d}ust}}$ $=1.53^{+0.13}_{-0.13}$ $M_{\odot}$, and
${T_{\mathrm{d}}}$ $=17.83^{+0.60}_{-0.57}$ K if
$\kappa_{0}(850\mu\rm{m})$=0.07 m2 kg-1 (using the empirical measurement of
$\kappa_{\rm{abs}}$ assuming the fraction of metals locked in dust is constant
across the local ISM, James et al., 2002). As this is a significant amount of
dust, here we also fit the SED using a wide variety of compositions and the
full characterization of $\kappa_{\rm{abs}}(\lambda)$, following Indebetouw et
al. (2014) and Matsuura et al. (2015).

We perform dust greybody fits to the FIR-mm SED using a wide variety of
$\kappa_{\rm{abs}}(\lambda)$ profiles directly obtained from the literature or
calculated using Mie theory; in total we used 134 $\kappa_{\rm{abs}}(\lambda)$
profiles. These include: Weingartner & Draine (2001) LMC average (as an
approximation to the conditions near SN 1987A in the LMC), Demyk et al. (2017)
amorphous silicate samples at 30K, Ormel et al. (2011) calculations of bare
and icy silicate+graphite grains, Ossenkopf & Henning (1994) bare and icy
mantle grains (protostellar core coagulation models), and using Mie scattering
calculations for amorphous carbon (amC, Zubko et al., 1996; Rouleau & Martin,
1991), graphite (Draine & Lee, 1984), ‘cosmic silicates’ (amorphous FeMgSiO4,
Jaeger et al., 1994), other silicates including enstatite and forsterite from
Henning & Stognienko (1996), pure iron (Henning & Stognienko, 1996), and 68
profiles from the Jena database (Henning et al., 1999) of minerals including
silicates, amorphous carbons, carbides, oxides, and sulfides. Many of the
models in the Jena database only extend to 500 $\mu$m or 1000 $\mu$m. In these
cases we extrapolate to our longest ALMA wavelength of 1329 $\mu$m in log
space (as lines, as most models follow power laws in the submm–mm), from the
last 300–500 $\mu$m of each curve to ensure a smooth continuation of the
general trend in each. We also consider a “continuous distribution of
ellipsoids” (CDE) model for a composite of carbons and silicates (Dwek &
Arendt, 2015), with a carbon volume filling factor fC of 18%.

For Mie theory calculations, we adopt a grain radius of $a=0.01$ $\mu$m for
PAHs, and $a=0.1$ $\mu$m for all other models. For most dust models, Mie-
derived FIR $\kappa_{\rm{abs}}$ curves are nearly identical for grain radii of
$a<5$ $\mu$m. A representative sample of 28 dust models was selected from this
list to span a wide variety of dust types. Selected grain types include Zubko
et al. (1996) amorphous carbons (amCs), Jäger et al. (2003) amorphous
silicates, and iron from Henning & Stognienko (1996), among others. The fitted
masses and temperatures for this sample are listed in Table 4.

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

####  4.2.2 Results of the dust fits

Fig. 6 shows 28 of the resulting dust fits to the FIR-submm SED (top panel of
the figure) using the various dust emissivity curves (bottom panel of the
figure) listed in Table 4. These have been grouped into similar dust
compositions: amC, graphite, PAHs, silicates, oxides, carbides, sulphides, and
pure iron to show the qualitative differences for fits using the various
composition types. (Individual models are not labeled; for quantitative
measures of specific models, we refer the reader to Table 4.) Fig. 6
illustrates that vastly different dust properties can still translate to
relatively similar fits to the observed flux densities. All models that result
in good fits to the SED have temperatures between $\sim$18–23K (Table 4); nine
of the 28 dust varieties listed here fail our $\chi^{2}$ criteria for a ‘good’
fit.

Models that represent amorphous carbons and graphites tend to have higher
$\kappa_{\rm{abs}}$ values and thus return lower inferred dust masses,
$\mathrm{M_{dust}}\sim$0.4–0.8$M_{\odot}$(Table 4). Amorphous silicates from
Demyk et al. (2017), which have more than a factor of 10 larger $\kappa_{\rm
850}$ compared to those from Semenov et al. (2003), also yield a relatively
moderate mass of $\sim$0.27$M_{\odot}$. The silicate models of Ormel et al.
(2011) and Ossenkopf & Henning (1994) also give moderate dust masses, and
represent grains that originate from coagulation processes as proposed in
dense molecular clouds of the ISM. These aggregated grains, some of which are
coated with ice, can be as large as 100 $\mu$m in radius, and the mass
absorption coefficients increase in the FIR and mm. As a consequence, the
inferred dust masses for these grain types are $\sim$0.5–0.74 $M_{\odot}$ for
Ormel et al. (2011) and Ossenkopf & Henning (1994). Because the dust grains
required for ISM coagulation are calculated to form over timescales of tens of
thousands to millions of years, it is unclear whether such icy grains can be
formed in the SN ejecta in such a short timescale as $<$30 years, however we
present the results of their fits here for comparison. Models of larger
composite silicate grains, such as those of Jaeger et al. (1994) or Jäger et
al. (2003), require even larger masses of $\sim$1–4 $M_{\odot}$ to fit the
observed SED. The CDE composite model of Dwek & Arendt (2015) results in a
moderate mass of 0.38$M_{\odot}$, due to the increased emissivity of the
carbon inclusion.

Dust varieties that generally satisfy our criteria for a good fit include
several amorphous carbon and silicate models (amorphous pyroxene and olivine
varieties), graphite, PAHs, and alumina. Varieties that tend to fit the data
poorly in the optically thin limit include FeO, FeS, SiC, SiO2, organics, and
water ice.

The largest source of uncertainty in the observed dust mass is the choice of
dust emission profile555We did not attempt to fit a two-temperature component
modBB as the SED shape is narrow, which suggests only one component is
necessary (Matsuura et al., 2011, 2015; Mattsson et al., 2015)., specifically
the value of the mass absorption coefficient ($\kappa_{\rm{abs}}$) in the
submm. We attempted to investigate spatial variations in the fitted dust
parameters across the ALMA maps, however the limiting beam size translates to
only 2–3 independent elements across the ejecta, and we found that the
differences in these flux ratios are smaller than their uncertainties.

We compare the inferred dust masses from different dust absorption
coefficients (Fig. 7). The inferred masses very clearly follow an inverse
linear relation based on the submm $\kappa_{\rm{abs}}$ value where
${M_{\mathrm{d}ust}[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.

Figure 6:  Top: Modified blackbody fits to the Herschel and ALMA observations
of the ejecta continuum in SN 1987A, using the 28 dust varieties and
parameters listed in Table 4. The Herschel 70µm and 500µm flux densities are
used as upper limits. Most dust models give reasonable fits to the observed
flux densities. The colors denote the mineralogy of the dust models: amorphous
carbon models are shown in red, graphite in orange, PAH in yellow, silicates
and silicate composites in blue, oxides (including FeO, TiO2, Al2O3,
NaAlSi2O6, and MgAl2O4) in light blue, carbide (SiC) in pink, sulfide (FeS) in
green, and pure iron in light purple. Bottom: Dust
$\kappa_{\rm{abs}}(\lambda)$ curves used for the modified blackbody fits,
showing the large variation in emission profiles.  Figure 7:  Fitted dust mass
as a function of $\kappa_{850}$ from various dust models. The colors for the
different dust varieties are the same as in Fig. 6. The fitted dust mass
closely follows a linear inverse trend over 3 orders of magnitude:
M${}_{\mathrm{dust}}\propto{\kappa_{850}}^{-1}$.

###  4.3 The Ring

The ring flux densities are shown in Fig. 5, along with previous measurements
of the ring from ALMA and the Australia Telescope Compact Array (ATCA): 1.4
GHz, 18 GHz, 44 GHz (Zanardo et al., 2013); 9 GHz (Ng et al., 2013); 94 GHz
(Lakićević et al., 2012). As the frequency increases towards the submm and
FIR, the contribution of the thermal ejecta emission to the total emission
becomes more significant, so the ring emission at day 9280 was estimated from
the total flux densities at 102, 213, 345, & 672 GHz by scaling and
subtracting an ejecta model component in Fourier space based on the Band–9
ejecta flux density at each frequency (Zanardo et al., 2014).

The Cycle 2 ring flux densities exhibit more scatter and are lower than the
ATCA values by $\sim 30$%. The integrated ring emission follows a non-thermal
power law profile of the form $S_{\nu}\propto\nu^{\alpha}$. The spectral index
$\alpha$ was previously estimated from ATCA radio (Zanardo et al., 2014) and
ALMA Cycle 0 data (Indebetouw et al., 2014) to be $-0.73\pm 0.02$ at day 9280
(Zanardo et al., 2014). An updated fit, using the ALMA Cycle 2 ring flux
densities and the more recent 1–9 GHz measurements from Callingham et al.
(2016) from day 9686, results in a slightly lower $\alpha=-0.70\pm 0.06$.

The radio ring emission has steadily increased due to the synchrotron-
producing electrons being accelerated by expanding shockwaves (Zanardo et al.,
2010; Staveley-Smith et al., 2014). Recently, Cendes et al. (2018) have fit
the radio emission of the (2D) ring and (3D) torus emission models across many
epochs as a power law of the form S($\nu,t$)=K$\nu^{\alpha}(t-t_{0})^{\beta}$,
where $\alpha$ is the spectral index of the emission across the spectrum,
$\beta$ is the power law slope, and K is the offset constant for the given
model. The ALMA ring flux densities are generally in excellent agreement
(within 5%) with the Cendes et al. (2018) model prediction for day 10402 where
K=1.5$\pm$0.1, $\alpha$=0.70 and $\beta$=0.59$\pm$0.02. We find no evidence of
a contribution from dust in the ring to the mm wavelength flux densities (see,
e.g., Bouchet et al. 2006; Matsuura et al. 2019).

##  5 Analysis and results of molecular lines

Figure 8:  Molecular line ratios of SiO. The line emission in each channel was
converted to flux [W m-2] integrated over the 300 km s-1 bin before division.
Channel velocity centers are LSRK; for reference, the SN 1987A systematic
velocity is 287 km s-1. Left: SiO $J$=7 $\\!\rightarrow\\!$ 6 / SiO $J$=5
$\\!\rightarrow\\!$ 4. Note that while the 300 km s-1 window covers the line
centers, the SiO $J$=7 $\\!\rightarrow\\!$ 6 line has poor coverage there due
to that frequency falling at the edges of the two observed tunings. Right: SiO
$J$=6 $\\!\rightarrow\\!$ 5 / SiO $J$=5 $\\!\rightarrow\\!$ 4. Intensity maps
were convolved to the SiO $J$=7 $\\!\rightarrow\\!$ 6 beam size before
dividing. Cyan contours are 679 GHz dust 3$\sigma$ and 5$\sigma$ levels. The
apparently high ratios at the edges of the map features are primarily due to
reduced S/N in those outer regions.

###  5.1 Analysis of Molecular Lines Using RADEX

In the previous Section, we modelled the integrated SED of the dust emission
under the assumption of uniform temperature and density within the ejecta.
However, in §3, we saw that the SiO $J$=5 $\\!\rightarrow\\!$ 4 and CO $J$=2
$\\!\rightarrow\\!$ 1 images exhibit a hole (Fig. 3) where the dust emission
peaks, and the SiO line ratio indicates lower SiO $J$=5 $\\!\rightarrow\\!$ 4
brightness with respect to SiO $J$=6 $\\!\rightarrow\\!$ 5 and SiO $J$=7
$\\!\rightarrow\\!$ 6 (seen in the SiO line ratio maps in Fig. 8). In order to
understand the spatial distribution of the line ratios and intensities
qualitatively, we use the non-LTE (Local Thermal Equilibrium) line radiative
transfer code radex (van der Tak et al., 2007).

####  5.1.1 Description of the procedure adopted for SiO

This subsection describes the analysis of the SiO lines and intensities in
detail. The analysis for CO was carried out in a similar manner, thus it is
described only briefly in that subsection.

radex calculates the molecular line intensities, using the escape
probabilities from (Osterbrock, 1989), and the uniform sphere method for the
gas distribution was chosen for this calculation. The code involves
calculations of level populations, using the Einstein coefficients and
collisional cross sections of molecular lines assembled by the LAMBDA database
(Schöier et al., 2005), and we use H2–SiO collisional cross sections based on
the calculations by Dayou & Balança (2006). In the ISM, H2 is widely assumed
to be the dominant collisional partner in molecular clouds; however, that is
probably not the case for the ejecta of SNe. As a consequence of a series of
nuclear burning processes, the progenitor star’s core will have built up
layers of different newly synthesised elements, with hydrogen being depleted.
In the two layers containing abundant Si, the major elements are O and S
(e.g., Woosley, 1988), and the collisional partner of SiO is likely to be O2
and SiS. This would potentially change the collisional cross section by a
factor of 1–10, depending on the transitions (Matsuura et al., 2017).

One of 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 (1999). Previous analyses of SN
1987A with lower angular resolution suggested the range of 2.5–45 %
(Kamenetzky et al., 2013; Matsuura et al., 2017); thus, we assume that the
filling factor of 1–50 % is a reasonable range. In this analysis, we adopted a
column density grid 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. (2017) 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 $\\!\rightarrow\\!$ 4 and SiO $J$=6
$\\!\rightarrow\\!$ 5 were convolved to match the lowest spatial resolution –
that of the SiO $J$=7 $\\!\rightarrow\\!$ 6 beam – with a uniform pixel width
of 0$\farcs$015\. The convolved and regridded line ratio maps, made on a
channel-by-channel basis, are displayed in Fig. 8. Intensities were averaged
over 5$\times$5 pixels in order to increase the signal to noise ratio in the
SiO $J$=7 $\\!\rightarrow\\!$ 6 image. As the minor-axis FWHM of the beam is
0$\farcs$17 for SiO $J$=7 $\\!\rightarrow\\!$ 6, there is a small loss of
spatial information by this averaging.

One caveat: the continuum subtraction for SiO $J$=6 $\\!\rightarrow\\!$ 5 was
performed on the final imaged data cube, whereas for SiO $J$=5
$\\!\rightarrow\\!$ 4 and SiO $J$=7 $\\!\rightarrow\\!$ 6 the continuum
subtraction was done in $uv$ space before imaging (imcontsub vs. uvcontsub in
Casa). Continuum subtraction in $uv$ space is generally considered preferable
if the continuum dominates the line emission, and the difference could
slightly affect the ratios, but at high S/N as in the case of SiO $J$=6
$\\!\rightarrow\\!$ 5 the two methods should give similar results. As
discussed in § 3, there was no appreciable difference between the two methods
for the CO $J$=6 $\\!\rightarrow\\!$ 5 data.

The uncertainties in the measured SiO lines and line intensities are dominated
by calibration uncertainties, and we adopted 10 % of the line intensities for
this analysis. The intensity uncertainties that were measured as a fluctuation
of the ‘blank’ sky level are 3 % for SiO $J$=6 $\\!\rightarrow\\!$ 5 and SiO
$J$=5 $\\!\rightarrow\\!$ 4, and 5 % for SiO $J$=7 $\\!\rightarrow\\!$ 6 after
5-pixel averaging, so that the uncertainties based on the observations are
smaller than the systematic uncertainties from calibration errors. Because the
dominant uncertainties are systematic, these propagate in an asymmetric way:
i.e., if the actual line intensity of SiO $J$=6 $\\!\rightarrow\\!$ 5 were
higher than the measured value, the line ratio of SiO $J$=6
$\\!\rightarrow\\!$ 5/SiO $J$=5 $\\!\rightarrow\\!$ 4 would increase, but the
ratio of SiO $J$=7 $\\!\rightarrow\\!$ 6/SiO $J$=6 $\\!\rightarrow\\!$ 5 would
decrease. The Cycle 2 flux densities being systematically lower than Cycle 0
as discussed in § 4.1 is not an issue for this analysis, as the important
result here is that the relative change of temperature and density within the
ejecta can explain the difference in molecular line ratios.

####  5.1.2 Results of SiO analysis

Figure 9:  Locations of regions A and B used for the radex analysis of the SiO
transitions. The background image is the SiO $J$=5 $\\!\rightarrow\\!$ 4 line
at 0 km s-1 (LSRK) channel, which clearly shows the molecular hole at region
A. Region B is a representative region of the general SiO emitting ejecta.

Figure 10:  (a)–(c) Results of radex modeling of the SiO line ratios at region
A, representing the molecular “hole”. These calculations assume a column
density of $N_{\rm SiO}=10^{16}$ cm-2. (d) The combination of all three best-
fit curves. For each line, the solid curve denotes the best fit locus of the
grid of kinetic temperature and H2 density values to the observed line ratios,
while the dashed and dotted curves represent the 1 $\sigma$ background
uncertainties and 1 $\sigma$ systematic uncertainties, respectively.

Figure 11:  The SiO line ratios that match the measured ALMA values at the SiO
hole. The same as Fig. 10 (d) but for a column density $N_{\rm SiO}=10^{17}$
cm-2.

Figure 12:  The same SiO calculations as in Fig.10 (d) but for region B with
column density $N_{\rm SiO}=10^{16}$ cm-2.

We selected two representative spatial regions, shown in Fig. 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 $\\!\rightarrow\\!$ 4 and CO $J$=2
$\\!\rightarrow\\!$ 1 which we call region A, and the other is representative
of the neighboring ‘typical’ SiO ejecta, named region B. The actual hole
itself has almost negligible SiO $J$=5 $\\!\rightarrow\\!$ 4 line intensity,
so we chose the nearest possible point for region A. All ratios in this
analysis were determined from the 0 km s-1 LSRK channel, which has good S/N
for all three SiO lines, and the hole is clearly identified.

In region A, which has a hole in SiO $J$=5 $\\!\rightarrow\\!$ 4, we found
that a column density of $N_{\rm SiO}=10^{16}$ cm-2 is reasonable. For $N_{\rm
SiO}=10^{15}$ 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}$ cm-2 can match the measurements, but beyond $10^{18}$
cm-2 the predicted and measured line intensities do not match.

Fig. 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}$ cm-2. The collisional partner for the radex calculations
is H2; see the last paragraph of this section for discussion of this point. In
Fig. 10 (a)–(c), the colored contours show the calculated ratios, the black
lines show the measured SiO line ratios, and the 1 $\sigma$ background
uncertainties and 1 $\sigma$ systematic uncertainties are indicated with
dashed and dotted lines, respectively. Fig. 10 (d) compares the possible
ranges of the kinetic temperature ($T_{\rm kin}$) and collision partner
density ($n_{\rm coll}$) from the three SiO line ratios. The uncertainty of
the SiO $J$=7 $\\!\rightarrow\\!$ 6/SiO $J$=5 $\\!\rightarrow\\!$ 4 ratio has
the tightest constraint, so only the uncertainty of this ratio is plotted in
Fig 10 (d). Fig 10 (d) shows that the ratios of all three SiO lines can be
matched with a very similar set of $T_{\rm kin}$ and $n_{\rm coll}>10^{7}$
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}$ cm-3, the second possibility is a non-LTE
condition with a range of ($10^{6}<n_{\rm coll}<10^{7}$ 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}$
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_{\rm SiO}=10^{16}$ cm-2 to $N_{\rm
SiO}=10^{17}$ cm-2, the SiO ratios change much more gradually as a function of
the kinetic temperature ($T_{\rm kin}$) and collision partner density ($n_{\rm
coll}$) (Fig. 11), expanding the feasible parameter space. This is the effect
of higher optical depth. The sets of $T_{\rm kin}$ and $n_{\rm coll}$ from the
three different calculated SiO ratios do not overlap in Fig. 11, but it is
still possible to consider a wide range of $T_{\rm kin}$ at $n_{\rm
coll}>10^{7}$ cm-3, from 19–22 K within the 1 $\sigma$ uncertainty. In the
non-LTE range, the required $n_{\rm coll}>10^{7}$ cm-3 is more or less
comparable to that for $N_{\rm SiO}=10^{16}$ 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_{\rm SiO}=10^{16}$ cm-2 case.

At region B, the physical parameters are slightly different from those at the
hole (region A). Fig. 12 shows the possible parameter space which fits the SiO
ratios for region B; only $N_{\rm SiO}=10^{16}$ cm-2 gives a feasible range.
The most plausible temperature is $T_{\rm kin}$=18 K for LTE conditions
($n_{\rm coll}>10^{7}$ cm-3), and an alternative possibility is non-LTE with
$3\times 10^{5}<n_{\rm coll}<10^{7}$ cm-3 and $18<T_{\rm kin}<$200 K. In the
optically thin regime, the filling factor and the column density are inversely
correlated, thus the accuracy of the filling factor is limited by our column
density grid.

In summary, the difference in the modelled line ratios and intensities near
the SiO $J$=5 $\\!\rightarrow\\!$ 4 and $J$=6 $\\!\rightarrow\\!$ 5 hole
(region A) with respect to the ‘representative’ ejecta region B can be
explained in the following three ways. The first possibility is that the hole
region has a higher temperature (35 K) than the surrounding locations (19–22
K) with both having LTE conditions – i.e., high density of the collisional
partner. The second possibility also requires LTE conditions, but instead of
high temperature, the hole region has a higher column density in a much
smaller area, represented by a small beam filling factor. The third
possibility is that the entire area is non-LTE, i.e., a lower density of the
collisional partner, but with the hole having a factor of a few to a few tens
higher density of the collisional partner than the surrounding region. These
three explanations are not exclusive to each other; a mixture of these three
scenarios is possible.

The uncertainties involved with this analysis arise from uncertainties in the
collisional cross section. The collisional partner is most likely not H2, but
rather other molecules such as O2 or SiS, according to chemical models
(Sarangi & Cherchneff, 2013). Therefore, instead of higher H2 density at the
hole the collisional partner may be different in the hole and in region B –
for example, region B may be dominated by collisions with O2, whereas in the
hole the dominant partner could be SiS. However, as will be discussed in §
6.3.1, hydrodynamic simulations predict that such spatial distributions for
the Si and O atoms are unlikely, therefore, our conclusion that there is a
higher temperature, column density or density at region A is still valid even
with this uncertainty.

####  5.1.3 CO analysis and results

Figure 13:  Top: Regions chosen for comparing the CO $J$=2 $\\!\rightarrow\\!$
1 and CO $J$=6 $\\!\rightarrow\\!$ 5 fluxes and line profiles; each region
spans one CO $J$=6 $\\!\rightarrow\\!$ 5 beam. Regions of particular interest
are labeled 1–9. The shaded greyscale contours illustrate the CO $J$=2
$\\!\rightarrow\\!$ 1 emission (lower left panel of Fig. 3). Middle: A
qualitative comparison of the stacked CO $J$=2 $\\!\rightarrow\\!$ 1 (blue)
and CO $J$=6 $\\!\rightarrow\\!$ 5 (red) spectra in each region with a
representative scale shown in the bottom left corner. Bottom: A zoomed in view
of the spectra for regions 1–9 shown in the middle panel. The vertical dashed
line is the systemic velocity of SN 1987A, 287 km s-1 (LSRK).

Although the CO $J$=6 $\\!\rightarrow\\!$ 5 line does not have enough S/N for
a quantitative analysis on a pixel-by-pixel basis, we can sum the spectra in
independent (single beam-width) regions to aid our analysis. The top panel of
Fig. 13 shows the location of 20 regions selected across the CO-emitting
ejecta, overlaid on top of the CO $J$=2 $\\!\rightarrow\\!$ 1 emission. Nine
regions, each the size of the CO $J$=6 $\\!\rightarrow\\!$ 5 beam, are
highlighted as areas of interest, potentially probing different conditions
(numbered 1–9 in rows from left to right). The middle panel of Fig. 13
compares the summed spectra of CO $J$=6 $\\!\rightarrow\\!$ 5 and CO $J$=2
$\\!\rightarrow\\!$ 1 (with the latter convolved to the CO $J$=6
$\\!\rightarrow\\!$ 5 beam before integration) for these 9 regions showing
their location with respect to the CO $J$=2 $\\!\rightarrow\\!$ 1 emission.
The spectra are in units of mJy, having been spatially integrated over each
region.

The bottom panel of Fig. 13 provides a zoomed in view of these spectra for a
more detailed comparison. Across the CO ejecta, we see that the majority of
the line profiles in the CO $J$=2 $\\!\rightarrow\\!$ 1 and CO $J$=6
$\\!\rightarrow\\!$ 5 transitions are similar to each other in the scaled
spectra (see for example regions 6, 7, and 9 in Fig. 13). However, in regions
1 and 8, the CO $J$=6 $\\!\rightarrow\\!$ 5 profile is suppressed at negative
velocities with respect to the CO $J$=2 $\\!\rightarrow\\!$ 1 line. The CO
$J$=6 $\\!\rightarrow\\!$ 5 emission is also suppressed with respect to CO
$J$=2 $\\!\rightarrow\\!$ 1 in region 3, though this is across the entire
velocity profile. At the location of the CO molecular hole (regions 4 and 5),
we see strong CO $J$=6 $\\!\rightarrow\\!$ 5 emission at velocities of
–1000–+1000 km s-1 with respect to CO $J$=2 $\\!\rightarrow\\!$ 1 (bottom
panel Fig. 13); this continues to neighboring region 2. We note that regions A
and B from the SiO analysis fall within the CO map regions 4 and 8,
respectively, but they are not interchangeable. They were selected
independently, and are centered at different locations and serve different
purposes (region B is representative of the general SiO emitting properties
across the molecular ejecta based on the SiO line ratios, whereas Fig. 13
demonstrates that region 8 has very different CO gas properties to most of the
other regions).

Using this information, we carry out an analysis with radex similar to that
for SiO. Fig.14 shows the resulting radex calculations of the CO $J$=6
$\\!\rightarrow\\!$ 5/2 $\\!\rightarrow\\!$ 1 ratios. The three black and
white lines show the curves for flux ratio values of 38, 20 and 3,
respectively, corresponding to the higher, intermediate and lower ends of the
line ratios observed across the 20 regions. We note that the units of the flux
densities used to derive the line ratios in this (CO) Section are in Jy, and
are therefore different to the W m-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 analysis666The spectra and channel
maps in Figs. 13 and A.2 are spectral density units (mJy and mJy per beam). In
order to compare the line ratios in the integrated fluxes in $\rm W\,m^{-2}$,
the flux densities in mJy per velocity channel units need to be multiplied by
a factor of $\partial f$/$\partial$v=–$f_{0}$/c to account for the change from
integrating in velocity v to frequency $f$. For the CO 6 $\\!\rightarrow\\!$
5/2 $\\!\rightarrow\\!$ 1 ratio, multiply by
$f_{0,\mathrm{CO65}}$/$f_{0,\mathrm{CO21}}$ $\sim$ 3..

As demonstrated originally in Fig. 3 in Kamenetzky et al. (2013), the CO $J$=6
$\\!\rightarrow\\!$ 5 line is sensitive to temperature change. We propose here
that the CO $J$=6 $\\!\rightarrow\\!$ 5 suppression with respect to CO $J$=2
$\\!\rightarrow\\!$ 1 indicates the gas is at a lower temperature in regions 1
and 8 (where we also see a peak in the dust emission, § 3, Fig. 4) compared to
the surrounding regions. In these regions, the blue wing of the CO $J$=6
$\\!\rightarrow\\!$ 5 emission is lower, thus if the dust and CO $J$=6
$\\!\rightarrow\\!$ 5 are spatially coincident in regions 1 and 8, this could
imply that the CO and dust originate from a discrete region on the near side
of the ejecta, though this is speculative as we have no velocity information
on the dust. Due to the low S/N of the CO $J$=6 $\\!\rightarrow\\!$ 5 line, we
cannot specify the exact excitation temperatures in these regions. However,
the CO excitation temperature is higher near the CO $J$=2 $\\!\rightarrow\\!$
1 hole (regions 4 and 5) at velocities of $-$1000–1000 km s-1.

Figure 14:  Ratios of CO $J$=6 $\\!\rightarrow\\!$ 5 and CO $J$=2
$\\!\rightarrow\\!$ 1 flux densities in Jy units derived using radex for the
optically thin case. The black and two white lines show ratio values of 38, 20
and 3, respectively indicating high, intermediate and low CO 6–5/2–1 end
ratios observed in the 20 regions in Fig. 13.

##  6 Discussion

###  6.1 Dust Properties from the Integrated SED

Any dust model that produces a mass higher than the total abundance of metals
formed in the SN ejecta is clearly unphysical. Here we discuss whether the
observed, or predicted, metal yields formed in the ejecta of SN 1987A can be
used to rule out some of the dust varieties and compositions that produce good
fits to the SED (Table 4). Inferred dust masses of several tens of solar
masses, as in the case of TiO2 (Posch et al., 2003) which requires 82
$M_{\odot}$ of dust, or MgAl2O4 Fabian et al. (2001) which requires 122
$M_{\odot}$ for example, are clearly untenable as they are larger than the
progenitor star mass (18–20 $M_{\odot}$, Woosley, 1988). 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. (2013), the total mass of the key limiting metals C, N, Ne, Mg, Si, S,
and Fe is 0.77 $M_{\odot}$, with 1.21$M_{\odot}$ of oxygen. Focusing on
individual dust varieties, the limits can be further differentiated. The
Nomoto et al. (2013) 18 $M_{\odot}$ model predicts 0.149$M_{\odot}$ of carbon,
putting a limit on the mass of graphite and amorphous carbon grains, as well
as a rough limit for PAH varieties since their masses are dominated by carbon.
The carbonaceous grain model that gives both a good fit to the SED and
produces the nearest fitted mass to the predicted carbon yield is the amC (AC1
sample) from Rouleau & Martin (1991), with a dust mass of 0.43$M_{\odot}$. For
silicate dust, good fits to the SED produce masses of $0.3-0.7$ $M_{\odot}$,
though the yields from Nomoto et al. (2013) would result in a maximum combined
silicate metal mass (and therefore dust mass, limited by the available Mg) of
$<0.4$ $M_{\odot}$.

The iron yield from the Nomoto et al. (2013) 18 $M_{\odot}$ nucleosynthesis
model is 0.079 $M_{\odot}$ (56Fe only) or 0.085 $M_{\odot}$ (including all
isotopes); this is orders of magnitude less than the inferred pure iron dust
model (Henning & Stognienko, 1996), which requires 3.97 $M_{\odot}$ of dust to
fit the SED. The Woosley & Weaver (1995) 15 $M_{\odot}$ and 18 $M_{\odot}$
$Z=0.1Z_{\odot}$ progenitor models predict iron masses ranging from 0.14–0.20
$M_{\odot}$, where roughly half of the iron originates from 56Ni decaying to
Fe. This is an order of magnitude less mass than the fit requires. We can
therefore rule out a scenario where iron-rich grains in SN 1987A are producing
the bulk of the thermal emission due to either not fitting the SED (in the
case of FeO and FeS, Henning & Stognienko 1996; Henning et al. 1995) or
resulting in unrealistically high dust masses for the pure iron model.

Limits from the Nomoto et al. (2013) 18 $M_{\odot}$ explosive synthesis model
for some other common dust varieties include 0.337 $M_{\odot}$ for forsterite
(Mg2SiO4), 0.481 $M_{\odot}$ for enstatite (MgSiO3), 0.014 $M_{\odot}$ for
alumina (Al2O3), and 0.373 $M_{\odot}$ for silica (SiO2) – assuming 100% of
all isotopes are locked in dust grains. These are all notably lower than the
dust masses resulting from the modBB fits – e.g., 4.0 $M_{\odot}$ for
forsterite, 4.1 $M_{\odot}$ for enstatite and 0.9 $M_{\odot}$ for alumina
(Table 4). Given the discussion above, one can place a rough upper limit on
the total mass of dust that could form in the ejecta of SN 1987A of $<1.5$
$M_{\odot}$ given the available metal budget. From this it is possible to
immediately rule out a further seven dust varieties listed in Table 4 as
producing unphysical dust masses.

Many of the fits to the SED of SN 1987A require more mass in dust grains than
the mass of available metals for the corresponding species. This can be
explained if the SED is made up of a mixture of several species each
contributing to the overall dust budget, as originally proposed in Matsuura et
al. (2015). For example, locking all available mass into a combination of
C+MgSiO3+FeS would give a total metal mass, and an upper limit on the total
dust mass, of 0.72 $M_{\odot}$. However, taking the yields of the Woosley
(1988) 18 $M_{\odot}$ 0.1$Z_{\odot}$ progenitor model, for example, would give
0.55 $M_{\odot}$ of total metals available for dust formation for this
combination of species, or $\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.,
2017 model) and amorphous carbons as also concluded by Matsuura et al. (2015),
limiting the total dust mass to potentially $<$0.7$M_{\odot}$.

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

Our spatially resolved images (Fig. 3) show that the dust distribution is
clumpy and asymmetric. Comparing the spatial distribution of the dust with the
CO $J$=6 $\\!\rightarrow\\!$ 5 reveals a weak anti-correlation with the
integrated CO $J$=6 $\\!\rightarrow\\!$ 5 and dust distributions, while there
is little spatial correlation between the dust and SiO images. We suggested
earlier that this anti-correlation occurs because CO $J$=6 $\\!\rightarrow\\!$
5 is suppressed compared with CO $J$=2 $\\!\rightarrow\\!$ 1 in the dust
bright regions, indicating that the excitation temperature of CO is lower than
in other neighboring regions. This provides new information about the
chemistry and physics involved in the formation of dust.

Note that we see an exception to this in one region (region 4 in our CO
analysis and roughly corresponding to region A in our SiO radex analysis).
Here, the hole in SiO and CO $J$=2 $\\!\rightarrow\\!$ 1 coincides with the
dust peak, with relatively strong CO $J$=6 $\\!\rightarrow\\!$ 5 emission
observed at $-$1000–1000 km s-1 with respect to CO $J$=2 $\\!\rightarrow\\!$
1, and this strong CO $J$=6 $\\!\rightarrow\\!$ 5 continuing to its
neighboring region (region 5 in Fig. 13). We will discuss this region
separately in § 6.3.1.

The SN chemistry after the explosion inherits a series of nuclear synthesis
processes at the stellar core prior to and during the SN explosion (e.g.,
Sarangi & Cherchneff, 2013). The outermost region is the H-envelope, followed
by the He shell, which can contain more carbon atoms than oxygen atoms
(Woosley & Weaver, 1995; Rauscher et al., 2002). The He shell can also form
CO, as it enriches with C and O (e.g., Sarangi & Cherchneff, 2013). The inner
region is roughly represented by an O+Ne zone, which can also contain C,
followed by an O+Mg+S+Si zone, and finally with a 56Ni 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 $\\!\rightarrow\\!$ 5 spatial
distributions as the result of both CO and dust components having originated
from the He-envelope or O+Ne nuclear burning zone containing both C and O
prior to the explosions. Thus we propose that the dust grains formed in this
region could be carbonaceous.

Our suggestion of carbonaceous dust contradicts predictions of supernova dust
formed in chemical models. Sarangi & Cherchneff (2015) predicts that SiO
molecules formed in the O+Mg+S+Si zone, and SiO molecules directly condense
into silicate dust. The majority of nuclear zones have more O atoms than C
atoms, and the formation of CO blocks the formation of graphite or amorphous
carbon. Deneault et al. (2006) and Clayton (2011) previously discussed the
formation of carbonaceous grains by dissociating CO via highly energetic
electrons, making unbound carbon available for carbonaceous dust formation.
However, the Clayton (2011) calculations involved only a few reaction rates
involving C and CO. In contrast, Sarangi & Cherchneff (2013) and Sarangi &
Cherchneff (2015) included more extensive chemical networks, as well as
dissociation of molecules by energetic electrons, and found this resulted in
very few carbonaceous dust grains ($6\times 10^{-3}$ $\rm M_{\odot}$ of
carbonaceous dust out of total of 0.04 $\rm M_{\odot}$ dust mass for a 19 $\rm
M_{\odot}$ star). One simple explanation is that Sarangi & Cherchneff (2013)
and Sarangi & Cherchneff (2015) modelled the SN chemistry only up to day 1500,
and the dust composition might have changed since then. However, it still
requires CO to be dissociated in order to remove the blockage of carbonaceous
dust formation, 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., 2017), which was largely under-predicted in the chemical
model of Sarangi & Cherchneff (2013), this suggests some tensions exist in
molecular chemistry models of SNe. However, this tension between our proposed
C-rich grain formation and the lack of C-rich grains “grown” in the SN dust
models could be alleviated if macroscopic mixing is efficient enough to allow
Si and C dust to form in the same regions. Indeed, some 3D hydrodynamical
models (e.g., B15, N20; see Utrobin et al., 2019) suggest that $\sim$30–70% of
the Si can be mixed out of the central regions into the C shell (A.
Wongwathanarat, private communication).

An alternative explanation for the anti-correlation between dust and CO
suggested in this work is that carbonaceous dust could originate from the He
layer while the CO gas is restricted to the C+O core, possibly resulting in a
projected anti-correlation between dust and CO. However, this does not explain
why the gas temperature is lower at the dust emitting region (though, the
presence of dust can lower the gas temperature, see below). In this scenario,
the CO might have originated from two different layers of nuclear burning
zones.

Our ALMA observations have also shown that regions of bright dust emission
have lower CO excitation temperature than other regions. There are two
possibilities to explain this. First, more dust may have led to cooler gas
temperatures due to radiative cooling via dust emission, i.e., the cooler gas
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 Fig. 3, for example, in comparison with
the 315 GHz image, but the surface brightness is lower in those locations. The
CO-dust anti-correlation is not as obvious in these fainter, extended regions,
therefore we do not claim that the entire dust content of the SN 1987A ejecta
is carbonaceous in nature; some of the dust components within the system could
be associated with silicates as well.

The high resolution ALMA observations in this work show that the dust
distribution in the ejecta of SN 1987A is clumpy even at small scales (as also
seen in the CO and SiO ejecta shown originally by Abellán et al., 2017). Aside
from SN 1987A, we know from the knots and filaments observed across the
Galactic SNR Cassiopeia A (e.g., Milisavljevic & Fesen, 2015) that the gas
ejected in SN explosions can be clumpy. Interestingly, Sarangi & Cherchneff
(2015) included clumpiness in their chemical models and showed that clumpy
gas, compared with smoothly distributed gas, provides density enhancements in
the ejecta resulting in larger grain sizes for SN-formed dust.

###  6.3 Interpreting the Bright Point Source (the Blob)

####  6.3.1 The blob: predictions from hydrodynamical models

The bright dust peak (the blob) is observed at the location of the hole in the
CO $J$=2 $\\!\rightarrow\\!$ 1 and SiO $J$=5 $\\!\rightarrow\\!$ 4 line-
emitting ejecta, and we also see emission located in the hole in the SiO $J$=7
$\\!\rightarrow\\!$ 6 transition. To explain this additional SiO emission, our
analysis of the SiO line ratios gives three possibilities: a higher SiO
temperature, a high column density with high area filling factor, and high
density of the collisional partner. The higher intensity of the emission of
the dust blob compared to its surroundings could be explained by higher dust
(‘column’) density or by higher dust temperatures inside the blob. For both
dust and SiO molecules, the two key physical parameters are temperature and
density, where the latter is also associated with the column density.

We first discuss the possibility of enhanced density (of dust, as well as SiO
and CO) in the blob based on model predictions. To form the SiO and dust, a
significant amount of Si must be present in the blob region. The dust could
also be carbon based, such that instead of Si, C could also be enhanced.
However, recent hydrodynamical simulations (Gabler et al. in preparation)
compared with ALMA observations (Fig. 4 in Abellán et al. 2017) show, that in
the explosions of three out of the four simulated progenitor models, the final
density distributions of Si and C (each multiplied 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. (2016) (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
$\\!\rightarrow\\!$ 4 or CO $J$=2 $\\!\rightarrow\\!$ 1 (Fig. 4 in Abellán et
al. 2017). However, in the hydrodynamic simulations, one model based on the
progenitor model from Shigeyama & Nomoto (1990) (model N20 in Wongwathanarat
et al. 2015, Abellán et al. 2017) leads to a weaker reverse shock and, hence,
a weaker self-reflected shock. The latter then is not able to significantly
accelerate the central ejecta. Therefore, the central region of this explosion
model still has similar densities of Si and C compared to that of the fastest
moving Si- and C-rich ejecta. While the N20 model may offer the possibility of
a high central density, the light curve (Utrobin et al., 2015, 2019) disfavors
that model.

####  6.3.2 The blob: gas and dust heated by the compact source

To explain the higher SiO and CO gas temperature and increased brightness for
the observed properties of the blob, we suggest two possibilities: (i) gas
heated by a compact object or (ii) a clump heated by 44Ti decay. Here we argue
that the most probable explanation for the detected dust blob is that the
innermost part of dust and gas is heated by radiation from the compact object,
with an early development of a pulsar wind nebula.

44Ti was synthesized at the time of the SN explosion, and its decay energy is
the main source of the heating of the inner ejecta (Jerkstrand et al., 2011;
Larsson et al., 2016). SN explosion models show that 44Ti is located more or
less at a similar radial extent as 56Ni and 28Si (Woosley & Weaver, 1995;
Wongwathanarat et al., 2015, 2017) with almost identical bulk velocities,
though the modeled distribution of 44Ti 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 44Ti ejecta is redshifted,
suggesting that the bulk is moving away from the observer (Boggs et al.,
2015). Since the predicted distribution of 44Ti shows qualitatively similar
properties as Fig. 4 in Abellán et al. (2017), because it is subject to
effectively the same hydrodynamical history (Jerkstrand et al., in
preparation), it is unlikely that gas heated by 44Ti decay would be more
centrally distributed and also be co-located with a small blob of gas and dust
at one small region in the very inner ejecta.

Therefore, we suggest that the central blob is due to warmer ejecta, and the
most likely source of the heating energy is from the compact object (i.e.,
possibility (i)). The dust and gas in the blob could be directly heated by
X-rays from the surface of the compact object (Alp et al., 2018a), or the dust
could be heated by synchrotron radiation generated by the compact object.

In the latter case, previously it was shown that synchrotron radiation from
neutron stars can heat up the dust grains in pulsar wind nebulae, as seen in
Galactic SNRs including the Crab Nebula (Temim et al., 2006; Gomez et al.,
2012b; Owen & Barlow, 2015), G54.1$+$0.3 (Temim et al., 2017; Rho et al.,
2018), G11.2$-$0.3, G21.5$-$0.9 and G29.7$-$0.3 (Chawner et al., 2018). In the
Crab Nebula, the radiation from the pulsar wind can partially dissociate gas,
contributing to the formation of the molecules OH+ and 36ArH+ (Barlow et al.,
2013). SN 1987A may currently be undergoing this phase, i.e., just beginning
to develop a pulsar wind nebula in the innermost region. If the compact source
is a black hole (Brown et al., 1992; Blum, & Kushnir, 2016), instead of a
neutron star, jets from the black hole can also heat up gas and dust locally
(e.g., Russell et al., 2006).

Instead of synchrotron, direct thermal radiation from the neutron star can
also heat the gas and dust locally (Alp et al., 2018a). The SN 1987A compact
object is still surrounded by a dense metal-rich gas. Metals absorb and
scatter (soft) X-ray, UV and optical light efficiently, so that it would be
challenging to detect light from the compact source directly (McCray &
Fransson, 2016; Alp et al., 2018b). The absorbed energy is re-processed into
longer wavelengths, and eventually ends up heating dust and molecules.
However, this picture does not consider any radio emission produced directly
by the compact object.

From the neutron star kick inferred by the distribution of intermediate-mass
elements (e.g., Katsuda et al., 2018) and the redshifted 44Ti spectrum (Boggs
et al., 2015), the direction of the compact object is predicted to be moving
toward us and extending along the north-east direction in the sky projection
(e.g., Janka et al., 2017). The location of this blob is consistent with that
prediction. There is an offset between the dust blob and the estimated
location of the progenitor star as derived by Alp et al. (2018a):
$\alpha=05^{\mathrm{h}}35^{\mathrm{m}}27^{\mathrm{m}}.9875$,
$\delta=-69^{\circ}16^{\prime}11^{\prime\prime}.107$ (ICRS). The brightest
pixel in the 679 GHz emission is offset 72 mas to the east and 44 mas to the
north of the position of the progenitor star from Alp et al. (2018a). This
offset is $\sim$3–5 times the total alignment uncertainty, however, it is
still within the range of the neutron star kick. Zanardo et al. (2014)
proposed that the compact object may have traveled 20-80 mas from the site of
the SN, towards the west, in comparison to the circular radius of 100 mas used
in Alp et al. (2018a). The peak of the dust blob from our work falls within
this range (though it is at the upper end, at approximately 700 km s-1 in our
data, assuming the SN 1987A position of Alp et al. (2018a)).

Moreover, the location of the compact object and the brightest part of the
pulsar wind in the Crab Nebula are also not coincident (Weisskopf et al.,
2000; Gomez et al., 2012b). Although the Crab’s pulsar is powerful enough to
affect its environment dynamically whereas the compact object in SN 1987A
would be at an earlier evolutionary stage, nevertheless, a misalignment
between the estimated location of the progenitor star and the location of the
dust blob is not unprecedented.

The approximate temperature of the dust of the blob can be estimated from the
679 GHz flux density. The peak is a factor of $\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. (1996)
or to 26 K from 18 K for Jäger et al. (2003) dust grains would be required to
explain the peak).

We estimate the dust peak flux density as S679GHz = 3–5 mJy, depending on how
a 2D gaussian is placed on the 679 GHz map; the peak pixel flux densities are
at S/N$\sim$7 (above the RMS level), and the uncertainty is caused by the
source being blended with other nearby features. This flux density contains
contamination from the underlying continuum emission, which can contribute
about 2 mJy, thus, the estimated flux density of the compact source is of
order 1–2 mJy. Figure 15 shows our estimated 679 GHz flux density range of 1–2
mJy (orange bar), along with spectra for the Crab Nebula and its central
pulsar (Bühler & Blandford, 2014), scaled to the distance of SN 1987A, for
comparison. The estimated blob flux density falls between the total spectrum
of the pulsar wind nebula and the sole pulsar spectrum, using the Crab as a
template. Using our 1–2 mJy range and assuming a dust model, one can estimate
the millimeter spectrum and the luminosity of the compact source. For the
Zubko et al. (1996) ACAR model and a temperature of 33 K, we obtain the light
orange shaded region in Fig 15. Alp et al. (2018a) derived upper limits to the
flux densities of the compact object, shown here as blue triangles. At 213
GHz, our ACAR modBB curve gives Sν 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 Lbol,dust=40–90 L⊙ for
679 GHz flux densities of 1–2 mJy. This is an order of magnitude estimate,
incorporating uncertainties in the flux density measurements and temperature
estimate. This luminosity range is an upper limit assuming that the compact
source heats dust from 0 K to 33 K. However, the compact source is most likely
additional heating, on top of 44Ti-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⊙.

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 44Ti; (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.

Figure 15:  Limits on the luminosity of the compact object in SN 1987A. The
flux densities for the entire ejecta are shown by the yellow curve. The
estimated 679 GHz compact object flux density range (1–2 mJy) is denoted by
the vertical orange bar. The light orange shaded region shows the Zubko et al.
(1996) ACAR model for this flux density range, and assumes it corresponds to a
temperature of 33K. The flux densities assuming this ACAR model and 1–2 mJy at
679 GHz are consistent with the limits from Alp et al. (2018a) (blue-grey
triangles). For comparison, the spectra of the Crab Nebula and the Crab Pulsar
(purple and red curves, respectively; Bühler & Blandford, 2014, and references
therein) are shown scaled to the distance of the LMC.

###  6.4 Dust as a Source of Extinction

The precise origin of the H$\alpha$ hole, whether from a physical lack of
material (due to the reverse shock or a simple void due to the expansion of
ejecta), or illumination of the ring X-rays brightening the outer rim of the
ejecta, or from dust extinction, has been discussed for many years (e.g.,
McCray, 2003; Larsson et al., 2011; Fransson et al., 2013; Larsson et al.,
2016, 2019). A simple estimate of the extinction can be made by assuming the
dust fills a spherical shell with uniform density, and calculating the optical
depth using the mass extinction coefficient $\kappa_{\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\,ds$. For
shell edges corresponding to the features in the images
(0$\farcs$05–0$\farcs$2, or $\sim$0.4–1.5$\times 10^{15}$m), the Zubko et al.
(1996) ACAR model and its corresponding dust mass fit give an opacity at 6563Å
of $\tau_{\mathrm{H}\alpha}\sim$ 560\. Silicate dust varieties tend to have
lower $\kappa_{\mathrm{ext,\lambda}}$ at a given wavelength – the amorphous
forsterite model of Jäger et al. (2003), for example, gives
$\tau_{\mathrm{H}\alpha}\sim$ 400\. $\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$m, the optical depth of ACAR dust at H$\alpha$ is 10. Even this
level is optically thick, meaning that dust extinction can play a non-
negligible role in the observed H$\alpha$ distribution.

##  7 Summary

We have observed SN 1987A with ALMA in Cycle 2, 10352–10441 days after the
explosion, in Bands 7 and 9. This follows on from the first ALMA results from
Cycle 0 in Kamenetzky et al. (2013); Indebetouw et al. (2014); Zanardo et al.
(2014) and Matsuura et al. (2015), and complements the molecular line studies
of Matsuura et al. (2017) and Abellán et al. (2017) from Cycles 2 and 3. In
this paper we describe the observations, data reduction, calibration, and
photometry in the ALMA bands at the highest angular resolution to date for the
continuum of SN 1987A.

  * •

We find that the dust emission in the ejecta is clumpy and asymmetric, fitting
within the H$\alpha$ “key hole”, with a peak of emission that we name “the
blob”. The dust ejecta region is smaller in scale than the cool, clumpy CO and
SiO ejecta regions. Dust, in the amounts we fit here in a simple uniform
spherical shell geometry, is optically thick at optical wavelengths, with
$\tau\sim$500\. Dust extinction could be a factor in the appearance of the
H$\alpha$ hole.

  * •

We see an anti-correlation between the CO $J$=6 $\\!\rightarrow\\!$ 5 emission
and dust, and this anti-correlation is not seen when comparing dust and CO
$J$=2 $\\!\rightarrow\\!$ 1. Our radex analysis suggests this is the result of
a lower CO gas temperature where the dust emission is stronger, compared to
the surrounding ejecta, hinting that the dust may be C-rich and may have
formed due to dissociation of CO, contrary to chemical predictions.

  * •

We observe a dust peak (the blob) at the location of the molecular hole
detected at lower CO and SiO transitions, and the higher CO and SiO line
transitions are stronger in that location than the emission in the lower
transitions. We suggest that this is the result of warm gas and dust at the
location of the blob, and discuss the possibility that this could be due to
slow moving reverse shock material originating from the explosion (the self-
reflected shock), heating from a high concentration of radioactive decay, or
the compact source. The most likely scenario is an indirect detection of the
compact source.

  * •

We fit the spectral energy distribution of the dust emission from the ejecta
with modified blackbody profiles. The derived dust masses and temperatures
depend on the submm emissivity of the dust, which is not very well determined.
Temperatures from the fits are generally between 18–23K. Amorphous carbon and
graphite models have dust with high emissivity and so yield the lowest dust
masses, around $0.4-1.6~{}{\rm M}_{\odot}$. Silicates return higher dust mass
estimates in the range $0.6-4~{}{\rm M}_{\odot}$. “Typical ISM” dust varieties
from the Milky Way and nearby galaxies give masses between $1-2~{}{\rm
M}_{\odot}$. Taking the total mass of available metals excluding oxygen
predicted to be ejected by a SN with progenitor mass appropriate for SN 1987A,
we rule out several grain models and compositions for the ejecta dust. We
revise the possible range of dust masses in the ejecta to $0.2-0.4$$M_{\odot}$
for carbon or Si grains, or a total of $<0.7$$M_{\odot}$ for a mixture of
grain species. A mixture of dust species, including silicates and carbonaceous
grains, seems necessary to reconcile the continuum SED, nucleosynthesis model
yields, and the molecular line analysis.

Acknowledgments: We acknowledge Arka Sarangi, Loretta Dunne, Steve Maddox, and
the anonymous reviewer for interesting and useful discussions. PJC and HLG
acknowledge support from the European Research Council (ERC) in the form of
Consolidator Grant CosmicDust (ERC-2014-CoG-647939, PI H L Gomez). MM
acknowledges support from the UK STFC Ernest Rutherford fellowship
(ST/L003597/1). Research at Garching was supported by the ERC through grant
ERC-AdG No. 341157-COCO2CASA, and by the Deutsche Forschungsgemeinschaft
through Sonderforschungsbereich SFB 1258 “Neutrinos and Dark 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. Atacama Large Millimeter Array

## References

  * Abellán et al. (2017) Abellán, F. J., Indebetouw, R., Marcaide, J. M., et al. 2017, ApJ, 842, L24
  * Ahmad et al. (2006) Ahmad, I., Greene, J. P., Moore, E. F., et al. 2006, Phys. Rev. C, 74, 065803
  * Alp et al. (2018a) Alp, D., Larsson, J., Fransson, C., et al. 2018, ApJ, 864, 174
  * Alp et al. (2018b) Alp, D., Larsson, J., Fransson, C., et al. 2018, ApJ, 864, 175
  * Arshakian & Ossenkopf (2016) Arshakian, T. G., & Ossenkopf, V. 2016, A&A, 585, A98
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Barlow et al. (2013) Barlow, M. J., Swinyard, B. M., Owen, P. J., et al. 2013, Science, 342, 1343
  * Begemann et al. (1997) Begemann, B., Dorschner, J., Henning, T., et al. 1997, ApJ, 476, 199
  * Blum, & Kushnir (2016) Blum, K., & Kushnir, D. 2016, ApJ, 828, 31
  * Boggs et al. (2015) Boggs, S. E., Harrison, F. A., Miyasaka, H., et al. 2015, Science, 348, 670
  * Bouchet et al. (1991) Bouchet, P., Danziger, I. J., & Lucy, L. B. 1991, AJ, 102, 1135
  * Bouchet et al. (2006) Bouchet, P., Dwek, E., Danziger, J., et al. 2006, ApJ, 650, 212
  * Brown et al. (1992) Brown, G. E., Bruenn, S. W., & Wheeler, J. C. 1992, Comments on Astrophysics, 16, 153
  * Bühler & Blandford (2014) Bühler, R., & Blandford, R. 2014, Reports on Progress in Physics, 77, 066901
  * Burrows (1988) Burrows, A. 1988, ApJ, 334, 891
  * Burrows et al. (2000) Burrows, D. N., Michael, E., Hwang, U., et al. 2000, ApJ, 543, L149
  * Callingham et al. (2016) Callingham, J. R., Gaensler, B. M., Zanardo, G., et al. 2016, MNRAS, 462, 290
  * Cendes et al. (2018) Cendes, Y., Gaensler, B. M., Ng, C.-Y., et al. 2018, ApJ, 867, 65
  * Chawner et al. (2018) Chawner, H., Marsh, K., Matsuura, M., et al. 2018, MNRAS,
  * Cherchneff & Dwek (2009) Cherchneff, I., & Dwek, E. 2009, ApJ, 703, 642
  * Cherchneff & Dwek (2010) Cherchneff, I., & Dwek, E. 2010, ApJ, 713, 1
  * Cigan (2019) Cigan, P. 2019, MultiColorFits: Colorize and combine multiple fits images for visually aesthetic scientific plots, ascl:1909.002
  * Clayton (2011) Clayton, D. D. 2011, New A Rev., 55, 155
  * Danziger et al. (1989) Danziger, I. J., Gouiffes, C., Bouchet, P., et al. 1989, IAU Circ.4746, 1
  * Dayou & Balança (2006) Dayou, F., & Balança, C. 2006, A&A, 459, 297
  * De Looze et al. (2017) De Looze, I., Barlow, M. J., Swinyard, B. M., et al. 2017, MNRAS, 465, 3309
  * Demyk et al. (2017) Demyk, K., Meny, C., Leroux, H., et al. 2017, A&A, 606, A50
  * Deneault et al. (2006) Deneault, E. A.-N., Clayton, D. D., & Meyer, B. S. 2006, ApJ, 638, 234
  * Draine & Lee (1984) Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89
  * Dwek et al. (2010) Dwek, E., Arendt, R. G., Bouchet, P., et al. 2010, ApJ, 722, 425
  * Dwek et al. (2014) Dwek, E., Staguhn, J., Arendt, R. G., et al. 2014, ApJ, 788, L30
  * Dwek & Arendt (2015) Dwek, E., & Arendt, R. G. 2015, ApJ, 810, 75