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

TL;DR

This paper presents high-resolution ALMA observations of SN 1987A, revealing detailed 3D molecular structures in the supernova's inner ejecta, providing insights into explosion instabilities and nucleosynthesis.

Contribution

It offers the highest resolution 3D mapping of molecular gas in SN 1987A's core, linking observed structures to progenitor and explosion physics.

Findings

CO and SiO distributions differ but both show a central deficit.

Clumpy molecular structures are observed and quantitatively compared to models.

Results constrain progenitor and explosion physics based on morphology.

Abstract

Most massive stars end their lives in core-collapse supernova explosions and enrich the interstellar medium with explosively nucleosynthesized elements. Following core collapse, the explosion is subject to instabilities as the shock propagates outwards through the progenitor star. Observations of the composition and structure of the innermost regions of a core-collapse supernova provide a direct probe of the instabilities and nucleosynthetic products. SN 1987A in the Large Magellanic Cloud (LMC) is one of very few supernovae for which the inner ejecta can be spatially resolved but are not yet strongly affected by interaction with the surroundings. Our observations of SN 1987A with the Atacama Large Millimeter/submillimeter Array (ALMA) are of the highest resolution to date and reveal the detailed morphology of cold molecular gas in the innermost regions of the remnant. The 3D distributions of carbon and silicon monoxide (CO and SiO) emission differ, but both have a central deficit, or torus-like distribution, possibly a result of radioactive heating during the first weeks ("nickel heating"). The size scales of the clumpy distribution are compared quantitatively to models, demonstrating how progenitor and explosion physics can be constrained.

Equations

$$\tau \approx \frac{{\lambda }^3{A}_{ul}{g}_{u}t}{8\pi Q}\left(e^{-E_l/T_x}-e^{-E_u/T_x}\right)n$$

Full text

ISM: supernova remnants — hydrodynamics — instabilities — supernovae: individual (SN 1987A)

1 Introduction

Supernova 1987A in the 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\farcs 5$ or $\sim 1$ light-year 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 outwards through the progenitor during the first $\sim 10^{4}$ s (Ebisuzaki et al., 1989; Benz & Thielemann, 1990). During the next days, radioactive decay of 56Ni$\rightarrow^{56}$Co$\rightarrow^{56}$Fe heats clumps rich in those heavy elements, causing some additional expansion of those clumps. If the process is sufficiently important, a large-scale reduction of heavy products from the central regions and lowest velocities may occur (“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 afterwards.

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 44Ti 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 ($\sim 500$ mas) ALMA data (Kamenetzky et al., 2013) revealed the presence of bright rotational emission from cold ($<150$ K) relatively slow moving ($\sim 2000$ km s*-1*) 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.

2 Observations and data reduction

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 September 2, 2014 (A001/X10e/X140), and cycle 3 high angular resolution images were obtained from November 1 to 15, 2015 (A001/X1ee/X620). Combining data from both cycles improves the image quality. Both data sets were calibrated against QSOs (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 12C16O $2-1$ (230.538 GHz) and 28Si16O $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$\lambda$) and 16,196 m (12,600 k$\lambda$). 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*-1* spectral bins (for spherical expanding ejecta at the time of these observations, 100 km s*-1* corresponds to a path length of $\sim 9\times 10^{15}$ 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 $scales=[0,7,21]$, and for the SiO $6-5$ image we adopted $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 ($\sim 28.7$ years). The assumption of free expansion can be tested in the spatially summed image. Using the full width at zero intensity (FWZI) $\sim 4000$ km s*-1*, and a corresponding angular size of the velocity-integrated image of $\sim 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.

3 Observational Results

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 ($\pm 1700$ km s*-1* vs. $\pm 1300$ km s*-1*; or $\pm 1.53\times 10^{17}$ cm vs. $\pm 1.17\times 10^{17}$ cm, approximately).

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 $\leq 1500$ km s*-1* ($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).

The observed low emission at the central region agrees with predictions from numerical models according to which different mechanisms (e.g., 56Ni-bubble effect, second outward shock, Ertl et al., 2016) would accelerate the innermost material outwards. The zones with the lowest emission ($<20\%$ of peak of emission) have similar sizes as the brightest clumps ($>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\leq 0\farcs 1$ or $r\leq 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 56Ni fingers predicted by some models.

The structure of the molecular emission is not aligned with the emission of 1.644 $\mu$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*-1* greater radial velocities.

4 Comparison with hydrodynamical models

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-15$ milliseconds after core-bounce to about 150 days we also include a description of the 56Ni heating. The additional heating increases the outward velocities of the innermost 12C and 28Si 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 (N$20\sim 3.8$ $M_{\odot}$ vs. 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 44Ti decay. CO and SiO emission may therefore reflect the 44Ti abundance distribution in addition to the CO and SiO abundance distributions. We have not yet detected variations in the CO or SiO excitation temperature which 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 12C and 16O, and of 28Si and 16O modeled atomic density distributions. We also analyzed 12C and 28Si 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$ K. Since the kinetic temperature is $T\sim 50$ K (Matsuura et al., 2017), the approximation is very reasonable. Assuming the same $T_{x}$ for all levels,

[TABLE]

For the CO $2-1$ transition, the Einstein coefficient $A_{21}=6.9\times 10^{-7}$ s*-1*. At $T\sim 50$ K, and at the time of the observation, $\tau$CO $\approx 2.7\times 10^{-3}$ $n$CO. Optical depth unity therefore corresponds to a CO density of $n$CO$\sim 374$ cm*-3*. 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*-1*, yields $\tau$SiO $\approx 4.2$ $n$SiO. Thus, optical depth unity corresponds to a SiO density of $n$SiO $\sim 0.24$ cm*-3*, 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_{\odot}$ and a few $\times 10^{-4}$ $M_{\odot}$, respectively, with uncertainties of a factor of several (Matsuura et al., 2017). These masses account for $<20\%$ and $<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 28Si currently in molecular SiO (Sarangi & Cherchneff, 2013). Thus, until more complete 3D models including dust formation exist, we compare the modeled 28Si$\times^{16}$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 28Si$\times^{16}$O has a more pronounced central hole than 28Si alone, more similar to the SiO data. Due to the fact that 16O extends further out than 28Si, their intersection has more of a hole. Instead, the differences between 12C$\times^{16}$O and 12C 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.

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.

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 12C$\times^{16}$O clump sizes comparable to the data. Instead, for SiO, the RSG models have a large range of 28Si$\times^{16}$O clump sizes, on average, $\sim 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\farcs 25$ consistent with a toroidal structure of that diameter, while the distribution of SiO shows a minor peak at $\sim 0\farcs 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 Fig. 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\farcs 22$, more similar to the CO data.

5 Conclusions

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 that were 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, instability physics and the duration of the unstable phase affects 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).