The messy death of a multiple star system and the resulting planetary nebula as observed by JWST

Planetary nebulae (PNe) are the ejected envelopes of intermediate-mass ($\sim$1–8 M_⊙) stars that have recently terminated their asymptotic giant branch (AGB) stage of evolution. Moving outwards from the hot pre-white dwarf star ($T\sim 10^{5}$ K) that is the progeny of the AGB star, the structure of a canonical quasi-spherical PN consists of a hot, sparse, wind-heated bubble ($T\sim 10^{7}$K) surrounded by a dense shell of displaced, ionised AGB gas ($T\sim 10^{4}$ K), which in turn may still be surrounded by “pristine,” cold ($T\sim 10^{2}$ K), molecule- and dust-rich AGB ejecta. On the other hand, if the progenitor star interacted with a companion(s) during its post-main sequence evolution, we would expect departures from spherical symmetry, perhaps including spiral structures and arcs (e.g., Mastrodemos1999, ; Mohamed2012, ; Maercker2012, ), the presence of a dense, molecule-rich torus (e.g., Santander-Garcia2017, ), one or more pairs of polar lobes formed by fast, collimated outflows and jets (e.g., Sahai1998, ; Sahai11, ), and/or a dusty, circumbinary disk VanWinckel2003 . The type of interaction depends on the orbital radius, and ranges from common envelope evolution for close binaries Ivanova2013 , to accretion disks and gravitational focussing of the wind for wider systems Mastrodemos1998 ; Mohamed2007 ; deValBorro2009 , to displacement of the central star from the geometric centre of the nebula for the widest systems Soker1999b .

The first Hubble Space Telescope (HST) images of PNe revealed a breathtaking new world of details and far more complex structures than had been gleaned from ground-based images (e.g., Balick1998, ; Sahai98, ). The superb spatial resolution of HST, combined with high-resolution, kinematic mapping, enabled the construction of detailed 3D, morpho-kinematic models, which, together with hydrodynamic models (e.g., Sabbadin2006, ; Steffen2006, ), started to connect our understanding of the evolution of the structures and kinematics of PNe with their possible binary star origins (e.g., Balick2002, ; DeMarco2009, ; Jones2017, ).

The James Webb Space Telescope (JWST), with its superb sensitivity and high spatial resolution from near- to mid-IR, is now poised to enable a leap of similar magnitude in our understanding of PNe. This journey began when JWST released near-IR and mid-IR images of just one PN, NGC 3132, as part of its ERO program. NGC 3132 is a nearby ($D\sim 750$ pc), molecule-rich Sahai1990 ; Kastner1996 , ring-like PN, long known to harbour a visual binary comprising the central (progenitor) star and an A star companion. In this paper we show that the JWST ERO images contain multiple, new lines of evidence that NGC 3132 is the recent product of a hierarchical multiple progenitor stellar system, which has experienced both indirect and direct interactions involving one or more components. Such binary interactions have taken on new importance in the era of gravitational wave detectors (LIGO Abramovici1992 , LISA AmeroSeoane2017 ) and ambitious transient surveys Ivezic2008 . Indeed, PNe like NGC 3132 offer unique insight into the formation pathways of the close, single and double degenerate binaries that are eventual gravitational wave sources and (perhaps) type Ia supernova progenitors (SantanderGarcia2015, ; Chiotellis2020, ; Cikota2017, ).

The first striking discovery of JWST is the presence of the dusty disk around the ultra-hot central star. This indicates that JWST can accurately detect dusty disks lighter than Ceres, as far as $\sim$700 pc away. For our PN, the presence of such a disk orbiting the PN central star favours a close binary interaction, where the companion either merged with the primary star, or is still in orbit but is undetected (mass $<$ 0.2 M_⊙; based on an unresolved or barely resolved, equal-brightness companion); in either case, the companion has donated a substantial fraction of its angular momentum to the gas Huang1963 ; Soberman1997 . Observationally, such disks around PN central stars, though rare, appear to be by and large associated with known or strongly suspected binarity Clayton2014 and may be related to circumbinary disks detected around other classes of post-AGB binary stars vanWinckel2009 .

An interacting binary scenario is reinforced by the shape of the ionised cavity, which represents the inner, most recent mass-loss phase, when the already hot central star emitted a fast, tenuous wind. Pairing the JWST images with spatially resolved spectroscopy we constructed a 3D visualisation of this cavity (see Morpho-kinematic modelling in Supplementary Material). In Figure 3 we show that this inner cavity is inferred to be an expanding prolate ellipsoid with its long axis tilted at approximately 30 $\deg$ to the line of sight. Its surface is not smooth and presents instead a number of protuberances, most of which can be paired via axes passing through, or very near the central star. Prolate cavities such as these, with misaligned structures, are common in PN and are likely sculpted by jets from interacting binaries in the earlier, pre-PN phase of the nebula Sahai2000 , with additional details added during the interaction between the AGB wind and post-AGB fast wind and via the process of PN ionisation.

The numerous protuberances clearly evident in the 3D reconstruction could arise from ionised gas breaking out of the inner cavity through an uneven outer shell. The apparent pairing between these protuberances may argue instead for the presence of intermittent and toppling jets (AkashiSoker2021, ). To generate jets over such a wide range of axes, an interacting binary is not enough, and one would have to conjecture that the central star is or was a member of not just a close binary, but of an interacting triple system (BearSoker2017, ). Recent studies of interactions in triple systems (Hamers2022, ; Glanz2021, ) also argue for the possibility of interactions yielding complex ejecta.

Outside the ionised ellipsoid, one encounters material ejected earlier in the star’s history. The AGB mass loss, at rates of up to $\sim$10^-5 M_⊙yr^-1 and speeds of $\sim$10 km s^-1 over a $\sim$10⁵ yr timescale (HofnerOlofsson2018, ), generates an enormous, expanding envelope of molecular gas and dust. The H₂ halo imaged by JWST constitutes the most recently ejected (inner) region of this AGB envelope. The spikes observed in the halo (Figure 1, right panel) show that the inner cavity is very porous, though less so near the minor axis where the cavity edges are brightest, densest, and least fractured.

The JWST images motivated 2D hydrodynamic simulations to replicate these flocculent structures. In Figure 4 we see two time snapshots towards the end of a simulation where an inner, faster wind from the heating central star and its ionising radiation, plough into the dense AGB (halo) material (see Methods, Section Hydrodynamic modelling). The fragmentation that happens at the interface of the swept-up material also creates the variable opacity needed to shield some of the wind material from ionising radiation, which then quickly recombines and allows the formation of molecules. Non ionising radiation leaks more readily because the opacity above 913 Å is lower. These photons produce florescence of H₂.

In Figure 4 we see two time snapshots towards the end of the simulation. In the first panel we see a set of approximately radial spikes, but 200 years later those straight and thin spikes evolve to thicker and sometimes curved ones. In the right column of Figure 4 two different parts of the nebula exhibit thinner and straighter spikes (top-right panel) or thicker, bent ones (bottom-right panel). Although the entire nebula was ejected and ionized over a short time interval, there can be a delay in the evolution of a given spike in a specific part of the nebula, related to the local opacity in the swept-up shell. Figure 4 suggests that differences of only $\sim$ 200 years in the timescales of mass ejection and/or the progress of illumination along specific directions can explain the marked differences observed in the flocculent structure around the nebula.

The successful modelling of illumination percolating unevenly into the molecular halo (Figure 4) motivated a further geometric model of the halo, presented in Figure 5. This Figure compares the extended H₂ structures as imaged by JWST with a model consisting of two thick, concentric, unbroken but clumpy, shells of material that are illuminated by the central star through a porous ellipsoid representing the boundary of the ionised cavity, with reduced opacity in the polar regions. As a result of the uneven illumination the distribution of H₂ material appears fragmented and is generally brighter toward the polar regions (and suppressed along the equatorial plane) of the central ellipsoidal, ionised region. The distribution seen in the JWST H₂ images could be reproduced more closely by altering the opacity of the inner ellipsoid. Fly-though movies of the 3D reconstructions of both the inner ellipsoid (Figures 3) and the outer H₂ halo (Figure 5) can be found following the links.

The arches in the JWST images, are not smeared as is typical of those seen in projection (e.g., Balick2012, ), but are instead sharp. This possibly indicates that these arches are on or near the plane of the sky, indicating that the orbit of the companion at $\sim$40-60 AU is closely aligned to the waist of the inner ellipsoid. This companion cannot partake in the formation of the disk around the central star, though it may play a secondary role in the shaping of other PN structures. It is also unlikely to have launched strong jets because at such distance the accretion rate would be very low. As such, this would be an additional companion to the inner binary (or triple), making it a tertiary (or quaternary) companion.

The visual A-type companion would then be a fourth (fifth) member of the group, an almost complete bystander from the point of view of interaction and shaping, but critically important for this study: Its well measured mass, and slight evolved status, constrained the initial mass of the central star: (2.86$\pm$0.06) M_⊙.

To reconstruct the events that lead to the demise of the progenitor of NGC3132, the PN acts like a murder scene. The A-type companion, could not have partaken to the interaction that unravelled the AGB star, but was (and is) certainly present. A second companion at 40-60AU left an indelible trail of its presence in the form of arcs, but was not close enough to generate the dusty disk, nor shape the ionised cavity, implying that there must have been at least another accomplice. This points the finger at a close-by companion, that is either avoiding detection, or has perished in the interaction (merged). If the numerous protuberances seen in the ionised cavity come in pairs, then tumbling jet axes would be needed and this would point the finger to the presence of a second, close companion (Hamers2022, ; Glanz2021, ), which would make the system a quintet. Even ignoring the putative second, close companion, we can state with good degree of certainty that the system is at least a quartet. Systems of four or five stars orbiting within a few $\times$1000 AU are not impossibly rare for primary stars in the progenitor mass range of interest here (e.g., HD 104237; Feigelson2003, ); indeed, present estimates indicate that 50% or more of stars of 2-3 M_⊙ are in multiple systems, and of order 2% of A-type stars have four companions (Duchene2013, ).

JWST is at the starting gate of its promise as an astrophysical pathfinder. With complementary radio, interferometric and time resolved observations, it can find the temporal signatures of active convective mass ejection from the surfaces of AGB stars and the subsequent gravitational influence of companion stars in dynamically- and thermally-complex outflows. Thus JWST offers the potential to intimately connect the histories of PNe and the role of close stellar companions to studies of chemical evolution, nebular shaping and binary interactions for the next century.

The inner, ionised cavity of NGC 3132 is elliptical in shape, with a major axis of $\sim$40 arcsec (0.15 pc) and an electron density of $n\sim 10^{3}$ cm^-3. The ionization structure and abundances were the subject of a recent study by MonrealIbero2020 . The nebula is also known to be molecule-rich Sahai1990 ; it is among the brightest PNe in near-IR H₂ emission Storey1984 ; Kastner1996 .

A bright A2 V star is found near the centre of the PN, but is too cool to be the ionizing star; the actual PN progenitor is much fainter and is located $\sim$1.7 arcsec to the South-West of the A star Kohoutek1977 ; Ciardullo1989 . The A2 V star has the same radial velocity and extinction as the PN, and its proper motion ($\mu_{\alpha}=-7.747$ mas/yr $\sigma_{\mu_{\alpha}}=0.026$; and $\mu_{\delta}=-0.125$ mas/yr $\sigma_{\mu_{\delta}}=0.031$) agrees with that of the central star ($\mu_{\alpha}=-7.677$ mas/yr $\sigma_{\mu_{\alpha}}=0.235$; and $\mu_{\delta}=0.197$ mas/yr $\sigma_{\mu_{\delta}}=0.275$), demonstrating that the PN progenitor and A-type companion constitute a comoving visual binary. The distance to NGC 3132 is obtained from Gaia DR3 measurements of this visual binary. No Gaia DR3 radial velocity is available for the optically faint central star (the PN progenitor). However, the brighter (A-type) visual companion and the PN have the same radial velocity: $(-11.4\pm 1.6)$ km s^-1 for the A star from Gaia, and $(-10\pm 3)$ km s^-1 for the PN from Meatheringham1988 . The A star and PN central star also have compatible Gaia DR3 proper motions (within 1.5$\sigma$).

The brighter, A-type star has a Gaia DR3 geometric distance (median of geometric distance posterior) of 754 pc, with lower and upper 1$\sigma$-like confidence intervals (16th and 87th percentiles of the posterior) of 18 pc and 15 pc respectively (Bailer-Jones2021, ). The fainter central star has a Gaia DR3 geometric distance of 2124.7 pc, with lower and upper 1$\sigma$-like bars of 559.1 pc, and 1464.5 pc. The quality flags of the astrometric solution for this star are not optimal, most likely due to the vicinity of the much brighter A-star; in particular, the goodness-of-fit along the scan is 16.9, while it should be close to unity. We therefore adopt the Gaia DR3 distance to the central star’s visual A-type companion, $754^{+15}_{-18}$ pc, as the distance to the PN.

The clumpiness of NGC 3132 in H₂ emission links this nebula to other molecule-rich PNe, such as the Helix Nebula (NGC 7293, Odell_etal_2004_HelixKnots ; Meixner_etal_2005 ; matsuura_etal_2007 ; Matsuura_etal_2009 ), Ring Nebula (NGC 6720, Kastner1994 ), and the hourglass-shaped (bipolar) nebula NGC 2346 (Manchado2015, ), in which the molecular emission seems to be associated with dense knots that are embedded in or surround the ionised gas. The origin of such H₂ knots in PNe — as overdensities in the former AGB wind, vs. formation in situ following recombination of H, as the central star enters the cooling track — remains an open question (Fang2018, ). In contrast to the Helix Nebula, there is little evidence for cometary tails emanating from the knots in the inner regions of NGC 3132. However, NGC 3132’s system of approximately radially-directed H₂ spikes external to the main H₂-bright ring system has close analogues in, e.g., the Ring and Dumbbell Nebulae Kastner1994 ; Kastner1996 .

Some H₂ knots in NGC 3132 are seen in absorption against the bright background nebular emission. This extinction is apparent not only in optical (HST) images but also, surprisingly, even in the JWST NIRCam near-infrared images (see Supplementary Figure 4). We measured the extinction at 1.87 $\mu$m for two knots seen in absorption against the (Pa$\alpha$) nebula background: the largest knot on the west side (coordinates 10:07:00.4, $-$40:26:08.8), and one of the darkest on the east side (10:07:02.5, $-$40:26:00.3). The diameters of these knots are $\sim$0.36 arcsec and $\sim$0.15 arcsec, while their extinction is $\sim$0.57 mag and $\sim$0.25 mag (at 1.87 $\mu$m), respectively; using the dust extinction law $A(\lambda)/A(V)$ from Cardelli:1989p1977 , the corresponding values of $A(V)$ are 3.9 and 1.7 mag assuming $R_{V}=3.1$. We then estimate the hydrogen column densities $N(\rm{H})$ from these extinction measurements, and convert to the hydrogen density $n(\rm{H})$ of the knot by assuming that the knot diameters are roughly equivalent to their depths along the line of sight. Using the conversion between $A(V)$ and $N(\rm{H})$ from Bohlin.1978 , where H is the combination of H⁰, H⁺ and H₂, the estimated column densities are $N(\rm{H})=7.3\times 10^{21}$ cm^-2 and $N(\rm{H})=3.2\times 10^{21}$ cm^-2, respectively. For the adopted distance of 754 pc, the estimated densities are $n(\rm{H})\sim 2\times 10^{6}$ cm^-3 for both knots. These densities suggest knot masses of $10^{-5}$ M_⊙, similar to the typical knot (“globule”) masses found in the Helix Nebula (Andriantsaralaza2020, ).

The critical density of excitation of the 2.12 $\mu$m H₂ 1–0 S(1) line at a kinetic temperature of 2000 K is 9$\times$10⁵ cm^-3 (Bourlot1999, ), if the collision partner is H. The critical density is higher for the 1–0 S(1) line than for the 0–0 S(9) 4.69 $\mu$m H₂ line (6$\times$10⁴ cm^-3; Wolniewicz_etal_1998, ; Bourlot1999, ). Hence, the excitation of H₂ should be nearly thermal if the gas temperature is sufficiently high, with the caveat that both critical densities are higher if the primary collision partner is H₂ rather than H.

To ascertain whether the mid-IR source associated with the PN central star is extended, we measured the JWST instrumental point spread function (PSF), using Gaussian fitting of field stars. We measured Gaussian FWHMs of 0.29, 0.40, 0.44 and 0.58 arcsec at 7.7, 11.3, 12.8 and 18 $\mu$m, respectively. We also measured two compact, slightly resolved galaxies in the field.

We then repeated the procedure for the central star. No fit was possible at 18 $\mu$m, due to saturation (see Supplementary Figure 5). At 7.7 $\mu$m the central star is on the edge of the diffraction spike of the A star, and only an upper limit on FWHM could be obtained. However, measurements of the PN central star image in the 11.3 and 12.8 $\mu$m filters gave consistent results, with measured FWHMs of 0.55 and 0.60 arcsec, significantly larger than the respective PSFs and comparable to the two field galaxies. Gaussian deconvolution using the PSF yields deconvolved FWHM values for the central star of $\leq 0.3$ ($\leq 230$ AU) at 7 $\mu$m, and 0.4 arcsec (300 AU) at 11.3 and 12.8 $\mu$m. The extent of the central star at 18 $\mu$m  is $\gtrsim 0.9$ arcsec in diameter (see Supplemenraty Figures 5 and 7).

We determined the mass of the A-star companion using version 1.2 of the PARSEC isochrones (Marigo2017, ) for solar metallicity, taken as $Z=0.0152$. We used $M_{\rm bol,0}=(0.34\pm 0.25)$ mag and the GAIA DR3 spectroscopic temperature $T_{\rm eff}=(9200\pm 200)$ K, where the errors are conservative. The star is confirmed to be beginning to turn off the main sequence, in a phase where the luminosity of $(57\pm 15)\,\rm L_{\odot}$ increases by 0.1% per Myr and the temperature decreases by 7 K per Myr (see Supplementary Figure 6). The isochrones yield an age of $(5.3\pm 0.3)\times 10^{8}\,$yr and a mass of $M_{\rm A2V}=(2.40\pm 0.15)\,\rm M_{\odot}$. The central star of the PN is evolving on the same isochrone, but from a more massive star as it has evolved further. We use the same isochrones to determine the initial mass of a star on the thermal-pulsing AGB, the phase where the central star ejected the envelope. This gives an initial mass for the central star of $M_{\rm i}=(2.86\pm 0.06)\,\rm M_{\odot}$. We have carried out the same isochrone fitting using an alternative stellar evolutionary model (the DARTMOUTH code; Dotter2008 ). Both the A2V star mass and the mass of the progenitor of the central star decrease by 0.15 M_⊙.

The final, CS, mass for such a star is 0.66 M_⊙. However, we have shown that such a star would show a high C/O$\sim$2, while the presence of silicate features in the Spitzer spectrum indicate that C/O$\lesssim$1. To reconcile the mass and the abundances we conjecture that the evolution was interrupted by the binary interaction that formed the disk, when the core mass was 0.61 M_⊙. With such a mass the evolutionary time to the current position on the HR diagram is in better agreement with the age of the nebula. This mass is also in better agreement with that derived from the photoionisation model (0.58$\pm$0.03) M_⊙.

The stratified ionisation and excitation structure of NGC 3132 is evident in Fig. 1, wherein the bright rim of ionized gas, as traced by [S iii] and Br$\alpha$ emission, lies nestled inside the peak H₂  emission. However, significant ionised hydrogen and high-excitation plasma — traced by [Ne ii] and [S iii] emission in the MIRI F1280W and F1800W filter images, respectively — is observed beyond the bright inner, elliptical ring.

We constructed a three-dimensional photoionisation model using the code Mocassin Ercolano2003 . To constrain the model we used the Multi Unit Spectroscopic Explorer (MUSE) emission line maps and absolute H$\beta$ flux of MonrealIbero20 , the optical integrated line fluxes from Tsamis2004 , the IR line fluxes from Mata_etal_2016 , as well as the velocity-position data obtained from the high-resolution scanning Fabry-Perot interferometer, SAM-FP, mounted on the SOAR telescope adaptive module. The observations were taken under photometric conditions. The seeing during the observations was 0.7 arcsec for the [N ii] observations to 0.9 arcsec for the H$\alpha$ one. The FWHM of a Ne calibration lamp lines was 0.586 Å or 26.8 km s^-1, which corresponds to a spectral resolution of about 11 200 at H$\alpha$.

We determined the density structure by fitting the emission line maps to the SAM-FP images of [N ii] $\lambda$6584 and H$\alpha$, using a distance of 754 kpc. The model adopts as free parameters the temperature and luminosity of the ionising source, and the elemental abundance of the gas component (assumed constant throughout the nebula); we assumed that no dust is mixed in the gas. For the ionising source we use the NLTE model atmospheres of central stars of planetary nebulae from Rauch2000 .

We find that a model invoking an unobscured central star with effective temperature $T_{\mathrm{eff}}=110$ kK and luminosity $L=200$ L_⊙ well matches the observational data. However, we find that the present-day central star mass implied by the comparison, between these stellar parameters and the evolutionary tracks of Ventura2018 ($0.58\pm 0.03$ M_⊙) is inconsistent with the (large) initial mass inferred from consideration of the presence of the comoving, wide-separation A-type companion (0.66 M_⊙; see Central star system’s masses). Furthermore, the tracks of Ventura2018 indicate that, for this mass, we would have a post-AGB age of 20 000 yrs, whereas the position-velocity data from the SAM-FP instrument yield an expansion velocity of 25-35 km s^-1  implying a much shorter and inconsistent nebular dynamical age in the range 2200–5700 yrs.

The C/O and N/O abundances of the nebula, as well as the crystalline silicate nature of the dust in the PN, indicate that this object has not undergone hot bottom burning and that it has not undergone sufficient dredge up to have increased the C/O ratio above unity. By the time the 2.86-M_⊙ star reaches the tip of the AGB its C/O ratio is approximately 2. It therefore seems that the mass implied by the initial-to-final mass relation using a main sequence mass of 2.86 M_⊙, is too high. We have two ways to resolve this inconsistency (which may both be operating). The central star is shielded by dust in the circumstellar disk making it appear, to the PN, as a cooler star, and/or the central star mass is actually smaller than 0.66 M_⊙, because the AGB evolution was interrupted by a binary interaction.

If the stellar ascent of the AGB was interrupted, we can determine the upper limit for a mass that would produce a nebula with $C/O\lesssim 1$ and N/O$\sim$0.4. This is 0.61 M_⊙. The time for a star of this mass to move from the AGB to the location on the HR diagram with an approximate temperature and luminosity (110kK, 200 L_⊙) as measured above is $\sim$10 000 yrs. The time-scale of the transition from AGB to post-AGB and PN is tightly connected with the rate at which the envelope is consumed: the results obtained are therefore sensitive to the mass-loss description. The time of 10 000 yr, is based on the classic mass loss rates dictated by Reimers or Blocker Bloecker1995 . This estimate must be considered as an upper limit of the duration of this phase; indeed the recent works on the AGB to post-AGB transition by Kamath2021 and Tosi2022 showed that to reproduce the infrared excess of post-AGB stars in the Galaxy and in the Magellanic Clouds one has to invoke significantly higher mass-loss rates than those based on the aforementioned formulations, something that would reduce the time-scales by a factor of $\sim$5. The timescale of 10 000 yrs is therefore easily reconciled with the observed timescale of 2200-5700 yrs implied by the nebula.

The hydrodynamic simulation used to interpret the fragmentation and radial spikes is a 2-dimensional hydrodynamic simulation using the magneto-hydrodynamic code ZEUS-3D. The computational grid is in spherical coordinates and consists of 800 $\times$ 800 equidistant zones in $r$ and $\theta$ respectively, with an angular extent of $90^{\circ}$. The wind and UV luminosity inputs correspond to a stellar post-AGB model with 0.677 M_⊙  which evolves from an initial 2.5 M_⊙  main sequence star Villaver2002 .

At simulation time 0 yr the star has $T_{\rm eff}=10\,000$ K and the AGB wind ($v=10$ km s^-1, $\dot{M}=10^{-6}$ M_⊙ yr^-1) has a homogeneous distribution outside of the pre-PN. The pre-PN has had 1000 yr of evolution prior to this moment, during which time a wide magnetic jet operated with a velocity $v=230$ km s^-1, and a mass-loss rate $\dot{M}=1.3\times 10^{-7}$ M_⊙ yr^-1; this simulation is taken from Model C6 in GarciaSegura2021 . At this time the star starts emitting a fast tenuous wind with a velocity $v$ from 240 to 14 000 km s^-1 and a mass-loss rate, $\dot{M}$ ranging from $1.06\times 10^{-7}$ to $1.13\times 10^{-10}$ M_⊙ yr^-1 over 4000 yrs that sweeps up the AGB wind material. At the same time (0 yr) the ionisation front propagates into the medium.