Very Deep inside the SN 1987A Core Ejecta: Molecular Structures Seen in 3D

Supernova 1987A in the Large Magellanic Cloud (LMC) has provided an excellent laboratory for supernova physics. Its short distance (50 kpc) has provided the very rare opportunity to spatially resolve the supernova as it evolves (currently the outer shock is greater than $1\buildrel{\prime\prime}\over{.} 5$ or ∼1 lt-yr in diameter). The environment of the ejecta of SN 1987A consists of a thin dense circumstellar ring, the equatorial ring, that is tilted 43° from the line of sight (Jakobsen et al. 1991; Tziamtzis et al. 2011). The first evidence of interaction of the outermost ejecta with the equatorial ring appeared in 1995 (Sonneborn et al. 1998; Lawrence et al. 2000), and was followed by the emergence of hotspots along the equatorial ring in the following years. These shocks and the reverse shock have provided much of the energy for the recent electromagnetic display from the supernova, as the decay of the radioactive energy sources produced in explosive nucleosynthesis provides ever decreasing amounts of energy.

At present, the ejecta emission at optical wavelengths is dominated by gas illuminated by the X-ray radiation from the ring (Larsson et al. 2011; Fransson et al. 2013). Our view of the inner debris where most of the nucleosynthesis products reside is highly obscured at optical/near-infrared by a massive cloud of dust grains (Matsuura et al. 2015; McCray & Fransson 2016). By contrast, millimeter molecular emission originates in the cold innermost ejecta and is free of the obscuration, thus being particularly powerful for understanding the details of elemental distribution and mixing.

Elements in the ejecta material are not simply radially stratified as might be expected from stellar evolution. Theory suggests that instabilities mixed the stellar nucleosynthetic products into a clumpy 3D structure. These began during the neutrino-driven phase of the explosion by convective overturn (Burrows et al. 1995) and/or standing accretion shock instability (SASI; Blondin et al. 2003), and grew by Rayleigh–Taylor instabilities as the shock propagated outward through the progenitor during the first $\sim {10}^{4}\,{\rm{s}}$ (Ebisuzaki et al. 1989; Benz & Thielemann 1990). During the next days, radioactive decay of ⁵⁶Ni $\to$ ⁵⁶Co $\to$ ⁵⁶Fe heated clumps rich in those heavy elements, causing some additional expansion of those clumps. If the process was sufficiently important, a large-scale reduction of heavy products from the central regions and lowest velocities may have occurred ("nickel heating"; Woosley et al. 1988; Herant & Benz 1991, 1992; Basko 1994). After the first few days, this structure was frozen in homologous free expansion, a record that can be analyzed long afterward.

Early detection of hard X-rays and the smooth shape of the light curve provided observational evidence that freshly synthesized material from the core of the supernova had found its way to the outer, less optically thick, regions of the envelope (Arnett et al. 1989; McCray 1993). Observations of the structure of atomic emission lines (Haas et al. 1990; Spyromilio et al. 1990) also suggested that such macroscopic mixing of the ejecta had taken place. Detailed spectral modeling (McCray 1993 and references therein), however, placed limits on the microscopic elemental mixing. Adiabatic expansion and radiative cooling caused temperatures to drop rapidly and allowed molecules to form where different atoms coexist. The molecular excitation temperature of 20–170 K (Kamenetzky et al. 2013; Matsuura et al. 2017) reflects the balance between this cooling and heating due to the gamma-rays and leaking positrons from the ⁴⁴Ti decay and possible X-rays and UV emission from the ring. Vibrational emission of CO and SiO was observed within the first 2 years (Rank et al. 1988; Spyromilio et al. 1988; Roche et al. 1991) and provided further constraints on the amount of mixing. However, none of these observations spatially resolved the heterogeneous ejecta. Modest angular resolution (∼500 mas) Atacama Large Millimeter/submillimeter Array (ALMA) data (Kamenetzky et al. 2013) revealed the presence of bright rotational emission from cold ($\lt 150\,{\rm{K}}$) relatively slow moving (∼2000 km s⁻¹) carbon monoxide CO and silicon monoxide SiO.

We have now observed SN 1987A with ALMA in the CO $J=2-1$, SiO $J=5-4$, and $J=6-5$ rotational transition lines using the long baseline, high angular resolution mode of the telescope. These observations are of the highest resolution ever taken for SN 1987A and enable us to create and compare the first 3D maps of molecular emission from deep inside of the remnant with state-of-the-art explosion models.

ALMA observations of SN 1987A at 1.3 mm (Band 6, 211–275 GHz) were performed in two different epochs: cycle 2 low angular resolution images were made in 2014 September 2 (A001/X10e/X140), and cycle 3 high angular resolution images were obtained from 2015 November 1 to 15 (A001/X1ee/X620). Combining data from both cycles improves the image quality. Both data sets were calibrated against quasi-stellar objects.

For the cycle 2 data set, J0519−4546 (05:19:49.72, −45:46:43.85; 0.75 Jy at 234 GHz) was the absolute flux calibrator and J0635−7516 (06:35:46.51, −75:16:16.82; 0.68 Jy at 234 GHz) was the phase calibrator. The cycle 3 high-resolution data consisted of two separate observations for each spectral setup. The ¹²C¹⁶O 2−1 (230.538 GHz) and ²⁸Si¹⁶O 5−4 (217.105 GHz) images include data using J0519−4546 (0.75 Jy at 224 GHz) as the absolute flux calibrator and J0601−7036 (06:01:11.25, −70:36:08.79; 0.70 Jy at 224 GHz) as the phase calibrator. The SiO 6−5 (260.518 GHz) image includes data using either J0519−4546 (0.63 Jy at 253 GHz) or J0334−4008 (03:34:13.65, −40:08:25.10; 0.44 Jy at 253 GHz) as the absolute flux calibrator and J0601−7036 (0.58 Jy at 253 GHz) as the phase calibrator.

Combining the cycle 2 and 3 data sets results in baselines between 34 m (25 kλ) and 16,196 m (12,600 kλ). We used the Common Astronomy Software Application (casa.nrao.edu) to process the interferometric data into a cube with 50 mas angular resolution and 100 km s⁻¹ spectral bins (for spherical expanding ejecta at the time of these observations, 100 km s⁻¹ corresponds to a path length of ∼9 × 10¹⁵ cm or 12 mas). For imaging and deconvolution we used the task tclean with different multiscale scales depending on the resulting beam sizes. For CO 2−1 and SiO 5−4 images, we adopted $\mathrm{scales}=[0,7,21]$, and for the SiO 6−5 image we adopted $\mathrm{scales}=[0,5,15]$. In both cases, we used a 6 mas pixel size.

To determine the 3D distribution, we converted velocity into angular size on the sky by assuming the distance to the LMC to be 50 kpc. We adopted a spherical free expansion, and for the age of the supernova we took the time elapsed since the supernova explosion and the date of the observations (∼28.7 years). The assumption of free expansion can be tested in the spatially summed image. Using the full width at zero intensity (FWZI) ∼ 4000 km s⁻¹, and a corresponding angular size of the velocity-integrated image of ∼480 mas, we obtain a time of 28.5 years since the explosion. There is no evidence of deceleration of the supernova ejecta, and therefore the free expansion approximation is valid.

In Figure 1, we present the CO 2−1 and SiO 5−4 emission framed by the circumstellar equatorial ring. Phenomenologically, the image provides strong clues to the structure of the inner ejecta of the supernova. The CO and SiO distributions are spatially distinct: both have a toroidal or shell-like distribution around the center, but most of the CO emission presents a maximum extension larger than SiO emission (±1700 km s⁻¹ versus ±1300 km s⁻¹; or $\pm 1.53\times {10}^{17}$ cm versus $\pm \,1.17\times {10}^{17}$ cm, approximately).

Zoom In Zoom Out Reset image size

The 3D distribution of CO 2−1 and SiO 5−4 is shown in Figures 2 and 3. Detailed examination reveals that the CO emission forms a torus-like shape perpendicular to the equatorial ring, clear evidence of asymmetry in the explosion. By contrast, SiO is clumpier and distributed in a broken shell rather than a torus. Within a region $\leqslant 1500$ km s⁻¹ ($1.35\times {10}^{17}$ cm) from the center of the explosion, we find that 25% of the clumps have peaks brighter than $1.9\times$ and $1.6\times$ the mean CO and SiO intensity, respectively. Translating the peak brightness of clumps into peak density of clumps requires modeling the non-LTE excitation of each clump, which we cannot do with only one emission line. However, if we use the average kinetic temperature and collider density derived for the entire nebula (Matsuura et al. 2017), then a peak brightness $3\times$ brighter than average would translate into an SiO density $5\times$ higher than average, for the typical clump size.

The brightest SiO emission is from a single blob located off-center below the equatorial plane (i.e., the side located closer to Earth). In both species, the emission joins into a more continuous torus or shell at about 50% of the peak of emission, and the central deficit of emission becomes fully enclosed at about 30% of the peak of emission (see Figure 3).

Zoom In Zoom Out Reset image size

The observed low emission at the central region agrees with predictions from numerical models according to which different mechanisms (e.g., ⁵⁶Ni-bubble effect, second outward shock; Ertl et al. 2016) would accelerate the innermost material outward. The zones with the lowest emission ($\lt 20 \%$ of peak of emission) have similar sizes as the brightest clumps ($\gt 70 \%$ of peak of emission) and a similar filling fraction of $\sim 25 \%$. The parts of the ejecta with the least SiO emission are preferentially found around the center ($r\leqslant 0\buildrel{\prime\prime}\over{.} 1$ or $r\leqslant 7.5\times {10}^{16}$ cm), whereas holes of emission in CO span a range of radii similar to the range spanned by the brightest CO clumps. We have examined the distribution of these zones of minimal emission and find that both CO and SiO holes are roughly perpendicular to the CO torus. One possibility could be that these regions are filled with material of different composition, like the heavy-element-dominated ⁵⁶Ni fingers predicted by some models.

The structure of the molecular emission is not aligned with the emission of 1.644 μm [Si i]+[Fe ii], which is concentrated in two asymmetric lobes fairly close to the plane of the ring and is brightest at larger radii than both CO and SiO (Larsson et al. 2016). Differences could be due to chemistry, molecular dissociation by positrons, different excitation mechanisms (thermal versus non-thermal), or dust obscuration. Figures 2 and 3 also reveal the first direct evidence of non-spherical instabilities, as SiO extends to greater radial distances (velocities) than CO in some directions. Considering emission brighter than 30% of the peak of emission, SiO has a greater extent than CO in 25% of radial directions, and in a few directions, extends up to 500 km s⁻¹ greater radial velocities.

We compare our data to 3D numerical models of neutrino-driven core-collapse supernovae obtained with the finite-volume Eulerian multifluid hydrodynamics code PROMETHEUS (Utrobin et al. 2015; Wongwathanarat et al. 2015) that employs an approximate, gray neutrino transport. All models considered are single-star models. To follow the simulations from 11 to 15 ms after core-bounce to about 150 days we also include a description of the ⁵⁶Ni heating. The additional heating increases the outward velocities of the innermost ¹²C and ²⁸Si by about 20%, so while it is not likely the only cause of the observed central evacuation, this "nickel heating" effect does contribute.

The progenitor mass for SN 1987A is in the range of 14–20 ${M}_{\odot }$ (Arnett et al. 1989; McCray 1993; Smartt 2009). We study four different non-rotating, non-magnetic progenitors: two 15 ${M}_{\odot }$ red supergiants (RSGs), model W15 (Woosley & Weaver 1995) and model L15 (Limongi et al. 2000), and two blue supergiants (BSGs) tailored to represent SN 1987A progenitors: a 20 ${M}_{\odot }$ model N20 (Shigeyama & Nomoto 1990) and a 15 ${M}_{\odot }$ model B15 (Woosley et al. 1988). The BSG models have a shallower density profile at the H/He interface than the RSG models and differ significantly in the C+O core mass (${\rm{N}}20\sim 3.8\,{M}_{\odot }$ versus ${\rm{B}}15\sim 1.7\,{M}_{\odot }$).

One-dimensional models of molecule formation in supernovae have been published (Lepp et al. 1990; Liu & Dalgarno 1996; Sarangi & Cherchneff 2013; Sluder et al. 2016), but there are no models that combine higher-dimensional hydrodynamic evolution of the ejecta with chemical processes. The emitted line radiation also depends on molecular excitation. This deep in the ejecta, it is reasonable to assume that excitation and heating are predominantly due to ⁴⁴Ti decay. CO and SiO emission may therefore reflect the ⁴⁴Ti abundance distribution in addition to the CO and SiO abundance distributions. We have not yet detected variations in the CO or SiO excitation temperature that would result from inhomogeneous excitation, but such variations may be detectable with future analysis of ALMA data.

Nevertheless, in order to focus on the size scales and spatial distributions, we compare the CO and SiO emission distributions directly to the square root of the product of the ¹²C and ¹⁶O, and of ²⁸Si and ¹⁶O modeled atomic density distributions. We also analyzed ¹²C and ²⁸Si individually, with and without density ceilings at which CO and SiO become optically thick. If the yields are high (a large fraction of the available atoms remains in gas-phase molecules), then the molecular densities will be high enough to make the observed CO and SiO lines optically thick.

We tested this in our structure comparisons by introducing an upper cutoff to the model volume density. For a rotational transition from level J at energy E_u to level $J-1$ at energy E_l, the Sobolev large velocity gradient optical depth is related to the molecular density n, rotational partition function Q, and excitation temperature T_x. For a diatomic rotator, like CO or SiO, the energy levels are ${E}_{J}={BJ}(J+1)$, where B is the rotational constant and ${g}_{J}=2J+1$ is the statistical weight. As long as T_x is not too low compared to B, the partition function can be approximated by $Q\sim T/B$, with $B/{k}_{B}=2.766\,{\rm{K}}$. Since the kinetic temperature is $T\sim 50\,{\rm{K}}$ (Matsuura et al. 2017), the approximation is very reasonable. Assuming the same T_x for all levels,

$$\begin{eqnarray}&&\tau \approx \displaystyle \frac{{\lambda }^{3}{A}_{{ul}}{g}_{u}t}{8\pi Q}({e}^{-{E}_{l}/{T}_{x}}-{e}^{-{E}_{u}/{T}_{x}})n.\end{eqnarray} \tag{ 1 }$$

For the CO 2−1 transition, the Einstein coefficient ${A}_{21}\,=6.9\times {10}^{-7}$ s⁻¹. At $T\sim 50\,{\rm{K}}$, and at the time of the observation, τ_CO $\approx \,2.7\times {10}^{-3}$ n_CO. Optical depth unity therefore corresponds to a CO density of n_CO ∼ 374 cm⁻³. This is within a factor of order unity of the maximum model densities, so does not affect our conclusions about relative structure size distributions. The same calculation for the SiO 5−4 transition, with ${A}_{54}=5.2\times {10}^{-4}$ s⁻¹, yields τ_SiO ≈4.2 n_SiO. Thus, optical depth unity corresponds to a SiO density of n_SiO ∼ 0.24 cm⁻³, a rather more severe threshold compared to model densities. From lower angular resolution ALMA observations, the total CO and SiO masses are a few $\times {10}^{-2}$ M_⊙ and a few $\times {10}^{-4}$ M_⊙, respectively, with uncertainties of a factor of several (Matsuura et al. 2017). These masses account for $\lt 20 \%$ and $\lt 0.1 \%$ of the expected C and Si yields (Woosley & Weaver 1995). Models suggest that formation of silicate dust is very efficient in this environment, leaving only a small fraction of ²⁸Si currently in molecular SiO (Sarangi & Cherchneff 2013). Thus, until more complete 3D models including dust formation exist, we compare the modeled ²⁸Si × ¹⁶O distribution without a cutoff with observed SiO.

The results obtained for the clump sizes and distributions are robust to the different approaches. The main difference is that ²⁸Si × ¹⁶O has a more pronounced central hole than ²⁸Si alone, more similar to the SiO data. Due to the fact that ¹⁶O extends further out than ²⁸Si, their intersection has more of a hole. Instead, the differences between ¹²C × ¹⁶O and ¹²C alone are less significant, as expected since these nucleosynthetic products start out at more similar radii before the instabilities take place. Figure 4 shows the phenomenological comparison of the data and the four models homologously scaled from 150 days to 29 years and convolved to the observational resolution.

Zoom In Zoom Out Reset image size

Different quantitative metrics are sensitive to different aspects of a 3D hierarchical structure. Segmenting the cube into clumps using cprops (Rosolowsky & Leroy 2006) and clumpfind (Williams et al. 1994) algorithms reveals the clump size distribution. On the other hand, the radially averaged 3D Fourier transform, and the distribution of distances between pairs of points above a given intensity threshold (thresholds between 50% and 90% of peak of emission were considered), reveal the range of clump separations and any large-scale structure such as the toroidal arrangement of CO. Figure 5 shows the clump size and clump separation for data and models of CO and SiO. In order to help illustrate the metric, simple geometries are also shown in Figure 5.

Zoom In Zoom Out Reset image size

The results that are most robust to different thresholds and segmentation parameters are the following: none of the models fit both emission lines in all structural metrics, but for CO, the RSG models W15 and L15 have ¹²C × ¹⁶O clump sizes comparable to the data. Instead, for SiO, the RSG models have a large range of ²⁸Si × ¹⁶O clump sizes, on average, ∼35% larger than the typical SiO clump sizes. For both lines, the clump sizes of BSG models N20 and B15 tend to be too small compared to the observations.

The clump separation of CO shows a clear peak at $\sim 0\buildrel{\prime\prime}\over{.} 25$ consistent with a toroidal structure of that diameter, while the distribution of SiO shows a minor peak at $\sim 0\buildrel{\prime\prime}\over{.} 30$. For both CO and SiO, the RSG models do not have such a clear peak, but instead they have a broader range of intermediate scales, reflecting the intermediate-sized structures visible in Figure 4. For SiO, both W15 and L15 models have a large-scale shell or torus comparable in extent to the data. The BSG models are more fragmented into small structures than the observations. Model N20 has small clumps with a wide range of clump-to-clump distances, and a larger shell structure than the data. For CO, the BSG model B15 does not have a range of intermediate distances but shows a strong peak at $\sim 0\buildrel{\prime\prime}\over{.} 22$, more similar to the CO data.

We have used the extraordinary capabilities of ALMA to create unique 3D maps of CO and SiO molecular emission from the inner ejecta of SN 1987A with unprecedented spatial resolution. Our observations show the clumpy mixed structure predicted by models (e.g., Hammer et al. 2010) but not previously well-imaged as well as a clear sign of asymmetry in the supernova explosion. We find that the molecular emission forms a torus-/shell-like shape perpendicular to the equatorial ring, and in some directions, the SiO extends further out than CO from the center of the remnant, a proof that non-spherical instabilities have taken place at the time of the explosion.

From the comparison of these data to hydrodynamical models we conclude that none of the models correctly reproduce neither the radial extent of both lines nor the range of clump sizes. The first discrepancy might be solved by adjusting the explosion energy of the model, but the second is directly related to the progenitor structure and the strongest asymmetries produced during shock revival, or to as of yet unmodeled chemical or excitation effects. However, considering that the models are not particularly tuned to reproduce SN 1987A, the agreement in some features as the overall shape of CO and SiO emission is surprisingly good.

The models compared here are non-rotating, and thus we cannot say how different the signatures of SASI would be. Comparing these data with a larger range of models, also including SASI mechanisms, will constrain explosion physics more precisely than has been possible heretofore using unresolved images and spectra, as different model parameters affect the size scales of clumpy structure in different ways: the explosion energy changes the overall extent, whereas the types and durations of instabilities affect the small-scale clump sizes. Progenitor structure and astrochemistry (fractional molecular yield compared to elemental abundance) also change the range of clump sizes and extent of their envelopes.

This Letter makes use of the following ALMA data: ADS/JAO.ALMA# 2013.1.00280.S and ADS/JAO.ALMA# 2015.1.00631.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC 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. F.J.A. and J.M.M. were partially supported by the Spanish MINECO projects AYA2012-38491-C02-01 and AYA2015-63939-C2-2-P and the Generalitat Valenciana projects PROMETEO/2009/104 and PROMETEOII/2014/057. At Garching, funding by the European Research Council through grant ERC-AdG No. 341157-COCO2CASA is acknowledged. P.C. acknowledges support from the European Research Council (ERC) in the form of Consolidator Grant CosmicDust (ERC-2014-CoG-647939, PI: H. L. Gomez). J.L. acknowledges support from the Knut and Alice Wallenberg foundation. M.M. is supported by an STFC Ernest Rutherford fellowship (ST/L003597/1). S.W. acknowledges support from NASA (NNX14AH34G).