Evidence for the Disruption of a Planetary System During the Formation of the Helix Nebula

Our understanding of the afterlives of planetary systems, once the host star has evolved off the main sequence, is predominantly based on modeling (e.g., Villaver & Livio 2009; Mustill & Villaver 2012; Villaver et al. 2014; Veras et al. 2017; Mustill et al. 2018; Veras & Fuller 2019; Maldonado et al. 2020; Veras & Heng 2020). Observing the architectures of planetary systems around white dwarfs would critically constrain these models. We have various lines of indirect evidence that such systems exist, including infrared excess emission (e.g., Jura et al. 2007; Xu & Jura 2012; Xu et al. 2015) and contamination of white dwarf spectra with metal absorption lines (e.g., Jura 2006; Xu et al. 2013, 2014). Spitzer mid-infrared surveys of hot white dwarfs found ≃20% of the objects located in planetary nebula exhibit infrared excess (Bilíková et al. 2012).

White dwarf absorption spectra are one of the few ways in which the composition of planetesimals around stars can be determined, yielding important clues to the composition of planetary companions to other stars (Jura & Young 2014). Surviving multiple planets bound to the stellar host are believed to play an important role in scattering planetesimals into the white dwarf's Roche limit producing the observed excesses and contamination (e.g., Debes & Sigurdsson 2002; Maldonado et al. 2021). The recent discovery of Jovian-mass exoplanets around white dwarfs has strengthened the case that the post-main-sequence survival of planetary system is indeed possible, and the mechanisms proposed to place material close to the white dwarf are plausible (Gänsicke et al. 2019; Vanderburg et al. 2020; Blackman et al. 2021; Scaringi et al. 2022).

The Helix Nebula (NGC 7293) is a polypolar, rather than bipolar, planetary nebula viewed at a pole-on orientation (O'Dell et al. 2004). Polypolar planetary nebulae exhibit structure associated with multiple, episodic bipolar outflow events at different orientations due to precession of the bipolar outflow (e.g., Kaler & Aller 1974; Manchado et al. 1996; López et al. 1998; Hsia et al. 2019); around 20% of young planetary nebulae are found to be polypolar (Sahai et al. 2011). Lying at a distance of 201 ± 3 pc, the Helix Nebula is one of the closest such systems to the Sun (Gaia Collaboration et al. 2016, 2018). Its central white dwarf (WD 2226-210) has a candidate low-mass binary companion in a few days' orbit based on Hα line emission (Gruendl et al. 2001) and photometric observations of periodic variability by the Transiting Exoplanet Survey Satellite (TESS; Aller et al. 2020). The most likely explanation for the variability is found to be irradiation of a substellar- or planetary-mass body that with a period of a few days that must have gone through common-envelope evolution.

WD 2226-210 is a young, DAO-type white dwarf with an effective temperature of around T_⋆ = 103,600 ± 5500 K, $\mathrm{log}g=7.0\,\pm \,0.2$ (Napiwotzki 1999; Traulsen et al. 2005), and an estimated age of ${10.6}_{-1.2}^{+2.3}$ kyr based on the expansion velocity of the nebula (O'Dell et al. 2002). More recent results record a slightly lower temperature T_⋆ = 94,640 ± 3349 K and similar surface gravity (Gianninas et al. 2011), but these are consistent with each other. For the purposes of this work we have adopted the values of Napiwotzki (1999) and Traulsen et al. (2005) to model the white dwarf spectrum. It is a powerful source of X-ray emission, providing a direct constraint on the density of material along the line of sight (Montez et al. 2015).

Binary post-asymptotic giant branch (post-AGB) stars, thought to be the progenitors of bipolar and polypolar planetary nebulae, are surrounded by gas- and dust-rich disks that resemble protoplanetary disks around pre-main-sequence stars (e.g., De Ruyter et al. 2006; Hillen et al. 2017). These post-AGB disks are sites of dust grain growth (Molster et al. 1999; Scicluna et al. 2020) and potentially planet(esimal) formation and growth (Bear & Soker 2014; Völschow et al. 2014). In combination, these properties make the Helix Nebula an excellent candidate to study the post-main-sequence fate of circumstellar material.

WD 2226-210 was resolved from the surrounding nebula at infrared wavelengths by Spitzer, revealing excess emission consistent with a circumstellar dust disk (Su et al. 2007). Subsequent observations at far-infrared wavelengths by Herschel supported this finding (Van de Steene et al. 2015), although those observations were focused on the wider nebula and not the white dwarf specifically. The presence of circumstellar emission from the white dwarf points to either the existence of remnant planetesimals from the main-sequence progenitor system (either asteroids or comets), a debris disk evolved to the post-main-sequence phase (Bonsor & Wyatt 2010; Veras & Heng 2020) or perhaps a secondary disk formed from stellar ejecta during its post-AGB phase if the system is indeed a binary (Kluska et al. 2022).

In this article we focus our efforts on deducing the origin and nature of the unresolved excess emission observed around the central white dwarf of the Helix Nebula. This is part of our wider effort to observationally constrain the post-main-sequence evolution of planetary systems (see also Ertel et al. 2019; Scicluna et al. 2020). The remainder of the paper proceeds as follows. In Sections 2 and 3, we summarize the observations and describe the methodology used to interpret the multiwavelength imaging and spectroscopy. Next, in Section 4, we present the results of multiwavelength image analysis and radiative transfer modeling. We then discuss the implications for the origin of the observed excess in Section 5, before giving our conclusions in Section 6.

Here we provide details of the new observations of the Helix Nebula we have obtained with the Stratospheric Observatory for Infrared Astronomy (SOFIA) and the Atacama Large Millimeter/submillimeter Array (ALMA), used to determine the brightness and spatial extent of the white dwarf's infrared excess emission. We present the imaging observations most relevant to this work, spanning mid-infrared to millimeter wavelengths, in Figure 1. A summary of the flux density measurements used to scale the white dwarf photosphere model and to fit the disk spectral energy distribution (SED), combining these new data with existing archival measurements, are provided in Table 1. In this table we break down the excess flux measurements of the target into three categories: "Total," "Extended," and "Compact." "Total" flux measurements are a combination of the observed flux (or upper limits) associated with both a point-like source and extended emission associated with the white dwarf position. "Compact" flux measurements denote detection of a point source (or 3σ upper limit) associated with the white dwarf position. "Extended" measurements for Herschel are related to the presence of a broad feature centered to the NW of the white dwarf with an angular extent of 28'' × 24'', while the "Extended" measurement for ALMA references the 3σ limit on the presence of a face-on, spatially resolved disk 3 beams in diameter. The reduction and analysis summarized here is not a complete and detailed overview of the process. The exact reduction steps used can be inferred from the scripts available, along with the reduced data, analysis scripts, and models, in a GitHub repository. ¹²

Zoom In Zoom Out Reset image size

SOFIA/HAWC+ Band A (53 μm) observations centered on WD 2226-210 were taken on 2017 October 17 for program 05_0054 (PI: S. Ertel). The observation was downloaded as a level 3 data product from the SOFIA archive. ¹³ The Lissajous map covers a roughly rectangular area 15 × 2' around the white dwarf, with an rms of 37 mJy beam⁻¹ for a total integration time of 27 minutes. The instrument's point-spread function (PSF) is 3'' FWHM, equivalent to a spatial resolution of 600 au at the distance of the Helix Nebula.

Deep Herschel/PACS small-map observations of the region around the Helix Nebula's white dwarf (PID: OT1_ksu_2; PI: K. Su) were taken using the 70 and 160 μm wave band combination. These unpublished, archival observations were obtained as level 2.5 (pipeline-reduced, mosaicked) data products from the Herschel Science Archive ¹⁴ (Verdugo 2019). For the analysis presented here we adopt the maps generated using the JScanamorphos algorithm as most reliable for the retention of extended structure. The map has nonuniform coverage, centered on the source, with a rectangular extent of $3^{\prime} \times 1\buildrel{\,\prime}\over{.} 5$ and an rms of 3 mJy within the central 1' radius of the map where coverage is highest and approximately uniform. The instrument PSF at 70 μm is 5.4'' FWHM, equivalent to a spatial resolution of 1000 au at the distance of the Helix Nebula.

ALMA Band 6 observations of WD 2226-210 were obtained in Cycle 3 (PID: 2015.1.00762.S; PI: S. Ertel) and Cycle 4 (PID: 2016.1.00608.S; PI: S. Ertel). These data sets were obtained from the ESO ALMA Science Archive. ¹⁵ In both cases the setup consists of four spectral windows providing a total of 7.5 GHz bandwidth to study the target emission. For the Cycle 3 (C3) observations two windows measure the continuum while the third and fourth windows cover the ¹²CO (2-1) and ²⁹SiO (5–4) lines. For the Cycle 4 (C4) observations, three windows lie over the continuum while the fourth covers the ¹²CO (2–1) line. Calibration and reduction of the ALMA observations were carried out in Common Astronomy Software Applications (CASA) using the appropriate CASA version (4.5.3 in C3 and 4.7.6 in C4). Image reconstruction was carried out using scripts supplied by the observatory, adopting Briggs weighting as a compromise between signal-to-noise and resolution. Both data sets had comparable spatial resolution, with a beam FWHM 025 along the major axis (50 au at 200 pc), sufficient to spatially resolve a solar-system-scale disk around the white dwarf. In the C3 data, the measured continuum sensitivity was 9 μJy with 3σ line sensitivities of 1.38 mJykm^-1s^-1 for ¹²CO (2–1) and 1.53 mJykm^-1s^-1 for ²⁹SiO (5–4). In the C4 data, the continuum sensitivity was 11 μJy with a 3σ line sensitivity of 1.74 mJykm^-1s^-1 for ¹²CO (2–1).

In summary, we have obtained new images of the central regions of the Helix Nebula centered on the white dwarf with SOFIA/HAWC+ at 53 μm and ALMA in Band 6 (≃1 300 μm). We complement these observations with previous archival imaging from Spitzer, presented in Su et al. (2007), and Herschel, presented in Van de Steene et al. (2015). The Spitzer data show a compact source at infrared wavelengths (from 8 to 70 μm) centered on the white dwarf position. The Herschel far-infrared observations (70 and 160 μm) from Van de Steene et al. likewise show compact emission associated with the white dwarf, but no emission is detected at submillimeter wavelengths (250 to 500 μm). The archival Herschel far-infrared observations (also 70 and 160 μm), which are deeper than those of Van de Steene et al., show compact emission associated with the white dwarf and its environment consistent with the Spitzer/Multiband Imaging Photometer for Spitzer (MIPS) 70 μm observation. The new SOFIA/High-resolution Airborne Wide-band Camera (HAWC+) 53 μm observation has a higher angular resolution than the previous Spitzer and Herschel 70 μm observations, but has much lower sensitivity and is a nondetection for the expected excess emission. The ALMA millimeter-wavelength observations are likewise nondetections in both continuum and line emission. In the next section we will combine the available detections and upper limits to model the excess emission and infer its likely origin.

Here we summarize the modeling approach used to constrain the spatial extent and density distribution of the circumstellar dust around the white dwarf and determine the dust properties based on the spatial constraints from multiwavelength imaging and source SED. Once we have an understanding of the dust emission, we then apply a radiative transfer model to determine the minimum dust grain size and size distribution.

3.1. Image Analysis and Modified Blackbody Fitting

Initially, we examine the assembled imaging observations of the Helix Nebula to determine the degree of excess emission as a function of wavelength and its spatial distribution; these observations are presented in Figure 1. We measure the extent of the extended disk using imaging observations at 24, 54, 70, and 1300 μm. From the 24 μm Spitzer observations of Su et al. (2007), the spatial extent of the proposed circumstellar disk is 30 to 100 au, spanning an angular diameter of 025 to 1''; this component is therefore unlikely to be spatially resolved in the SOFIA and Herschel far-infrared observations. We therefore consider that the system may comprise two components: a compact (point-like) source and an extended component. The compact component is modeled as a point source using a wavelength-appropriate PSF model, while the extended component is modeled as a two-dimensional (2D) Gaussian profile. The convolved model is then subtracted from the observation and significant residuals identified in the image.

We are aware that in the ALMA observations the point-like component may itself be spatially resolved, and therefore the disk orientation (surface brightness) becomes important as a constraint at millimeter wavelengths. A broad extended component may remain undetected by ALMA due to either low surface brightness or insensitivity to the appropriate angular scale due to the interferometer configuration. Furthermore, we also consider the presence of (spatially resolved) gas emission at millimeter wavelengths from the system. Two common species, carbon monoxide (CO) and silicon monoxide (SiO), were covered by spectral windows in the ALMA observations. These lines are diagnostic of emission from post-AGB envelopes and icy planetesimal belts. Both CO and SiO emission might be expected from a remnant post-AGB disk, whereas volatile rich planetesimals may leave a detectable CO emission in a debris disk. The geometry of the line emission would further inform our understanding of the origin of the excess emission.

We then proceed to model the SED. We first model the emission as a combination of modified blackbodies specific to the number of spatial components in the system obtained in the previous step (e.g., Wyatt 2008). The modified blackbodies are defined by temperature T_bb, break wavelength λ₀, and exponent β, such that

$$\begin{eqnarray}F(\lambda )\propto \left\{\begin{array}{ll}B(\lambda ,{T}_{{bb}}), & \lambda \leqslant {\lambda }_{0}\\ B(\lambda ,{T}_{{bb}})\times {\left({\lambda }_{0}/\lambda \right)}^{\beta }, & \lambda \gt {\lambda }_{0}\end{array}\right.,\end{eqnarray} \tag{ 1 }$$

with the following additional constraints:

• 1.  
  The rising Spitzer/InfraRed Spectrograph (IRS) spectrum at mid-infrared wavelengths is entirely accounted for by the compact component.
• 2.  
  The temperature of the compact component, being strongly constrained by the Spitzer/IRS spectrum (Su et al. 2007), is fixed as 102 K, and the temperature of the extended component must be lower than the compact component.
• 3.  
  The break wavelength (λ₀) for both components must be at least the shortest wavelength for which a detection exists (i.e., 24 μm for the compact component and 70 μm for the extended component).
• 4.  
  The beta exponents for both components must lie between 0 and 4, which covers the range observed for debris disks. The beta exponent may be different for the two components.

We determine the best-fit parameters for all components simultaneously. To do this we create 10,000 sets of values drawn from the appropriate ranges for each parameter (T_bb , λ₀, β) and scale the resultant modified blackbodies to the measured flux density of the system. Models that violated the upper limits for either component were discarded from the final ensemble. For the remainder that were consistent with the available data, we calculated the best-fit values for the two components in combination through error-weighted least-squares fitting, and the mean and standard deviation of the valid parameter sets in isolation.

3.2. Radiative Transfer Modeling

Here we present our findings regarding the structure and emission properties of the Helix white dwarf's excess emission, based on the multiwavelength imaging and photometry data sets. We first consider the extent of the emission region, combining archival mid- and far-infrared imaging from Spitzer and Herschel with new far-infrared SOFIA/HAWC+ and millimeter ALMA data. We then determine the spectral shape of the excess emission, informed by infrared to millimeter-wavelength photometry and upper limits. We then combine these constraints using a radiative transfer model to infer the dust grain properties necessary for consistency between them.

4.1. Imaging and Spectral Energy Distribution

Mid- and far-infrared excess emission close to the Helix central white dwarf was first discovered with Spitzer, revealing spatially unresolved emission located at the stellar position that was seen at wavelengths from 8 to 70 μm (Su et al. 2007); see Figure 1. The emission could be clearly disentangled from the bulk of the nebulas emission due to its compactness (diameter <2'' at 200 pc, based on the 8 μm image). Fitting the excess SED revealed it to have a temperature of ≥86 K and a corresponding size of several tens of au, comparable to the extent of the Edgeworth–Kuiper Belt, but it was not possible to constrain the dust properties from the limited wavelength range of the available data.

The SOFIA/HAWC+ Band A (54 μm) observation does not have any evidence of a point-like source at the white dwarf position, as shown in Figure 1. At first glance, this is seemingly inconsistent with the previously inferred dust temperature and morphology of the Spitzer imaging observations and photometry, including the rising Spitzer/IRS spectrum up to 35 μm. Since the Spitzer/MIPS and SOFIA/HAWC+ maps have comparable resolution (FWHM ≃ 6''), the nondetection in the longer-wavelength observation places a strong constraint on the total brightness of the excess emission, despite its limited sensitivity (rms 37 mJy). We can therefore infer that within 600 au of the white dwarf the total emission must be under 120 mJy.

In contrast, archival Herschel/PACS 70 μm imaging data (PI: K. Su) surprisingly show extended emission associated with the white dwarf's position. That emission can be decomposed into two components, which we have fitted using a PSF model for the compact component, with a flux density of 36 ± 4.5 mJy, and a 2D Gaussian for the extended emission, with a flux density of 222 ± 12 mJy, as shown in Figure 2. The brightness of both components is consistent with the nondetection in the SOFIA/HAWC+ 54 μm image due to the much greater sensitivity of Herschel/PACS at 70 μm. Furthermore, the level of excess emission from the extended component and its nondetection in submillimeter Herschel/Spectral and Photometric Imaging Receiver (SPIRE) maps of the wider nebula (Van de Steene et al. 2015) rules out the presence of a substantial mass of cold dust at larger distances from the white dwarf, which would have been undetected in the SOFIA/HAWC+ map due to dilution across multiple beams.

Zoom In Zoom Out Reset image size

The ALMA observations of the Helix white dwarf are also nondetections, both in continuum and line measurements. The largest angular scale for the array configuration was about 4'', which rendered the ALMA observations blind to the extended component seen in the Herschel/PACS map. However, these maps are deep, with a continuum rms around 10 μJy. The compact component of the emission should have been easily recoverable in the observations given its presumed extent based on the shorter-wavelength maps, if it followed a spectral slope consistent with debris disk emission (i.e., including large grains up to millimeter sizes). Therefore, the ALMA observations provide a strong constraint on the surface brightness of the excess emission at millimeter wavelengths. The rms sensitivities of the line observations were 0.38 mJykm^-1s^-1 for CO (2–1) and 0.51 mJykm^-1s^-1 for SiO (5–4). Assuming a line width of 10 km s^–1 (typical for debris disks), optically thin emission, and local thermodynamic equilibrium for the gas (T_gas = 100 K), we derive 5σ upper limits to the integrated line fluxes of 2.3 × 10⁻¹⁷ W m² for SiO (5–4) and 3.1 × 10⁻¹⁷ W m² for SiO (5–4). These limits are equivalent to gas masses of 1.2 × 10⁻⁴ M_⊕ in CO and 2.7 × 10⁻² M_⊕ in SiO. The nondetection of CO or SiO line emission from the white dwarf rules out the presence of substantial gas mass, contrary to expectations if the excess emission originated from a remnant post-AGB disk.

From the SED presented in Figure 3 we see that the shape of the compact component (purple) is strongly constrained by both the Spitzer/IRS spectrum and the Herschel/PACS photometry. In conjunction, these observations confine the shape of the excess to being a sharp rise followed by a sharp falloff. The SOFIA HAWC+ upper limit at 54 μm further constrains the total brightness of the excess such that the peak of emission must occur before that wavelength. However, the nondetection at (sub)millimeter wavelengths of the compact component leaves its spectral slope unconstrained, with the ALMA measurements offering only weak constraints on the shape if we assume the dust emission region would be unresolved (consistent with a compact Kuiper Belt analog around the white dwarf). By contrast, the extended component (blue) is only weakly constrained with a single measurement in the Herschel/PACS 70 μm map. Attributing the detected excess at shorter wavelengths to the compact component alone provides some restriction on the temperature of the extended component, while the weak upper limits in the submillimeter and millimeter from Herschel SPIRE and ALMA enable us to rule out a massive cold component to the total emission. The best-fit parameters of the modified blackbody models used to interpret the SED are given in Table 2.

Zoom In Zoom Out Reset image size

Table 2. Summary of SED Model Fitting

Parameter	Compact	Extended
Temperature (K)	102	50 ± 6
λ0 (μm)	25 ± 7	96 ± 29
β	1.4 ± 0.7	1.5 ± 0.2

Note. The temperature of the compact component was held fixed at 102 K.

Download table as:  ASCIITypeset image

The available data offer no real constraint on the nature of the extended component seen in the Herschel 70 μm map with weak upper limits either side of that to constrain its behavior. Multiple interpretations could fit the evidence equally well, such as it being a blowout halo from the compact component, part of the wider Helix Nebula coincidentally superimposed upon the white dwarf, or a diffuse remnant of the dispersing post-AGB envelope. However, the compact component extent is relatively well constrained to within a few 100 au of the white dwarf, and its SED is surprisingly well constrained with a steep rise, quick turnover, and steep falloff. In the next subsection we perform a further analysis of the compact component, applying these constraints in a radiative transfer model.

4.2. Radiative Transfer Modeling

We used standard debris disk modeling tools to fit the far-infrared SED of the compact excess (Ertel et al. 2012) and simulate spatially resolved observations (Ertel et al. 2011). This is appropriate to determine the spatial and grain size distribution of any accumulation of optically thin circumstellar dust (the fractional luminosity of the Helix excess suggests it is optically thin). For each model fit to the SED, a spatially resolved ALMA image of the emission was simulated and compared to the ALMA upper limit. We used simple power-law radial dust and grain size distributions and realistic dust emissivities commonly employed in debris disk modeling. A broad parameter space of the spatial and grain size distribution was explored, as illustrated in Figure 4, and the exact model parameters for each scenario are summarized in Table 3.

Zoom In Zoom Out Reset image size

Table 3. Summary of Model Parameters Used to Calculate SEDs Presented in Figure 4

Parameter	Model A	Model B	Model C	Model D
Disk inner radius, Rin (au)	60	30	30	30
Disk outer radius, Rout (au)	80	50	90	90
Exponent of surface brightness profile, α	0.0	−1.0	−1.0	−1.0
Minimum grain size, ${a}_{\min }$ (μm)	2	60	60	10
Maximum grain size, ${a}_{\max }$ (μm)	20	1000	1000	1000
Exponent of size distribution, γ	3.5	3.5	3.5	3.5
Mdust (×10−3 M⊕)	8	10	10	10

Download table as:  ASCIITypeset image

The lower limit of the grains size derived from SED modeling of debris disks around luminous main-sequence stars is similar to the blowout size (i.e., the minimum size for which dust grains can remain in a bound orbit around the host star under the influence of radiation pressure) of the dust grains inferred from the stellar properties (see Section 4.9.1 of Horner et al. 2020). For the Helix white dwarf (T_WD =103,600 ± 5500 K, L_WD = 67.6 L_⊙), this size is 60 μm. Such dust grain sizes are not consistent with the steep slope of the far-infrared excess emission toward longer wavelengths as constrained by the Herschel data. The Herschel data also rule out that the compact component is too extended to be detected by ALMA. Instead, the ALMA nondetection provides strong constraints on the excess brightness at millimeter wavelengths.

The steep slope of the SED can only be explained by the presence of small, micron-sized grains that dominate the emission but have a very low emissivity at long wavelengths (Ertel et al. 2012; Marshall et al. 2016). Such a grain size is similar to that of typical dust grains released by evaporating comets in the solar system (Mason et al. 2001). These grains cannot be on stable orbits around the Helix central star due to their removal by radiation pressure and Poynting–Robertson drag, and must be replenished continuously. The removal of the dust by radiation pressure could also explain the broad extended emission seen in the Herschel data, the radial distribution of which is consistent with such a scenario.

Here we consider the probability, based on the lines of evidence presented above, that the observed excess emission originates from a leftover planetesimal belt, a dissipating cometary cloud, or the remnants of the post-AGB circumbinary disk.

The white dwarf central star WD 2226-210 has a well-determined mass of 0.60 ± 0.02 M_⊙ (Benedict et al. 2009). Late-type binary companions as cool as late-M stars are not detected (Ciardullo et al. 1999) and the photometric variability with 2.77 days period detected using the TESS light curves excludes a secondary companion with a mass in the range 0.16 M_⊙ ≤ M_⋆ ≤ 2.5 M_⊙ (Aller et al. 2020). At this period the now-substellar secondary must have gone through common-envelope evolution. During the system's post-main-sequence evolution, any circumbinary planets or planetesimal belts would have been disrupted by the effect of mass loss from the primary, causing orbits to expand outward if located at a distance that avoids tidal engulfment (Villaver & Livio 2007, 2009). Assuming a standard initial-to-final mass relation for the primary, the Helix Nebula white dwarf must have a 1.5 M_⊙ stellar progenitor (Benedict et al. 2009). This assumes that the secondary star did not affect the mass-loss rate evolution, which is largely unknown if the system has experienced evolution through a common-envelope phase. The engulfment scenario is complex; whether or not a planet is engulfed by the expanding progenitor as it evolves from the main sequence is dependent not only on the progenitor mass and planetary semimajor axis, but also mass of the planetary companion. Different simulations demonstrate engulfment by a 1–2 M_⊙ progenitor for a gas giant planet with an initial orbital semimajor axis up to 4 au (for a 5 M_Jup planet). The orbital distance of the most distant planet that would be engulfed by this star is between 1.7 au (for a terrestrial planet) and 3 au (for a Jupiter-like planet; Villaver & Livio 2009; Mustill & Villaver 2012). Note, however, that planetary systems that host more than one planet are expected to undergo instabilities following post-AGB mass loss causing planets/planetesimals to be sent to orbits close to the WD tidal disruption radius (Mustill et al. 2018; Maldonado et al. 2021).

The nondetection of significant gas emission (either CO or SiO) in the ALMA observations of WD 2226-210 places strict constraints on the evolution timescale for any post-AGB binary disk in this system. The Helix Nebula has an estimated total mass of 0.9 M_⊙, with an ionized gas mass of 0.3 M_⊙ (Henry et al. 1999) and 0.6 M_⊙ in molecular gas (Young et al. 1999). This is a substantial amount of molecular gas, but there is no strong evidence for a remnant disk of material around the central white dwarf. Most of the H₂ is distributed at large angular scales beyond 13 up to 51 from the central white dwarf (Matsuura et al. 2009).

Observations of post-AGB disk systems have revealed disk masses ranging from 8 × 10⁻⁴ up to 10⁻² M_⊙ (Bujarrabal et al. 2018; Gallardo Cava et al. 2021; Kluska et al. 2022). These disks have short lifetimes, believed to be around 10⁴ but perhaps up to 10⁵ yr (Oomen et al. 2019); the nondetection obtained here is therefore consistent with the shortest expected timescales for the dissipation of the post-AGB circumbinary disks. As a hot, young white dwarf WD 2226-210 exhibits substantial radiation pressure, with a minimum dust grain size of ${a}_{min}\simeq 60\,\mu m$. Dust growth in the post-AGB circumbinary disk could create grains large enough to remain bound to the stellar remnant and withstand this radiation pressure, but dust grains of this size are too large to satisfactorily fit the observed SED. We therefore discard the idea of a post-AGB circumbinary disk as the potential origin for the excess emission.

The size distribution of dust grains in debris disks, produced by the collisional attrition of planetesimals, typically starts around 1 to 10 μm, at a few times the radiation blowout limit for the host star luminosity (Pawellek et al. 2014; Pawellek & Krivov 2015). The size distribution extends up to millimeter- or centimeter-sized pebbles, with a (sub)millimeter size distribution exponent q between 3 and 4 (MacGregor et al. 2016; Marshall et al. 2017; Norfolk et al. 2021). We find that the spectral slope is constrained by the ALMA upper limits to the top end of this range for compact configurations consistent with a planetesimal belt, particularly in the edge-on case. For more spatially diffuse emission, the limit on the spectral slope is relaxed, allowing for smaller q values, but then the architecture ceases to be consistent with a planetesimal belt and begins to look more like a shell (e.g., Matrà et al. 2018; Marshall et al. 2021). The nondetection of millimeter-wavelength emission from the system therefore counts against the planetesimal belt hypothesis for the origin of the observed excess.

The incidence of debris disks around main-sequence stars of comparable masses to the progenitor star of the Helix Nebula is around 20% to 30% (Thureau et al. 2014; Holland et al. 2017; Sibthorpe et al. 2018). Many young white dwarfs have infrared excesses associated with substantial dust- and gas-rich disks, with an incidence of a few percent. As the white dwarf cools the brightness of an attendant debris disk likewise drops, especially at infrared wavelengths, confounding detection of any dusty excess for systems beyond a few million years into the post-main-sequence phase of their evolution. Models of the post-main-sequence of A-star debris disks establish that, with current sensitivity constraints, one such system within 200 pc would be detectable by current facilities (Bonsor & Wyatt 2010). The Helix Nebula excess identified by Su et al. (2007) is the only example of such a system found to date, and is consistent with the expectations of Bonsor & Wyatt (2010).

The key differences between old white dwarfs and the Helix Nebula's white dwarf are its youth and its much higher luminosity. Around old white dwarfs relatively massive debris disks can persist, and are still undetectable because the dust is not heated to any noteworthy degree; this is not the case for any belt around the Helix Nebula's white dwarf (Bonsor et al. 2011). So it is much harder to hide a belt there, in particular if it is massive enough to produce the huge influx of material into the inner system over the age of the nebula through scattering by a planet. Potentially, we just do not see a planetesimal belt's excess because it is on the Rayleigh–Jeans tail of the warmer dust's emission and we have few, weak constraints on the total emission from the system at long wavelengths. However, the slope of that emission is already too steep to explain for that dust in a normal debris disk, hence we need dust grains too small for a steady-state debris disk, so it would be hard to hide much cool dust in there.

Steep submillimeter slopes have been observed at far-infrared or millimeter wavelengths for a handful of debris disk systems (Ertel et al. 2012; Marshall et al. 2016). The nature of the dust in these steep SED debris disk systems is still uncertain. However, in each case the minimum dust grain sizes inferred for those debris disk systems is consistent with the blowout grain size expected from radiation pressure (Burns et al. 1979), which is explicitly not the case here.

The main constraint that speaks against a debris disk as the origin for the observed excess around the Helix Nebula white dwarf is the small grain size required to fit the shape of the SED from 12 to 70 μm. The sharp rise and fall of the SED can only be fit satisfactorily with dust grains of a minimum size far below (by a factor ≃10) the minimum grain size expected from radiation pressure arguments. Such grains could not persist around the white dwarf due to its luminosity, so must continually be created, or delivered, by some mechanism absent of a remnant planetesimal belt. The size distribution of solar system cometary dust grains suggests a population of primarily submicron-sized dust grains with a narrow size distribution (Mason et al. 2001). This would result in a narrow, steeply declining SED as has been observed here. A cometary origin for the dust grains in the system is therefore consistent with the continuum observations.

However, the cometary cloud scenario might be expected to produce substantial amounts of volatiles such as water or CO from volatile outgassing from the planetesimals (as seen in several debris disks; Greaves et al. 2016; Kral et al. 2017; Marino et al. 2020; Rebollido et al. 2022). The water ice line (150 K) lies at 28 au. This is passingly consistent with the inner edge of the dust emission region (35 au) as inferred from the Spitzer/IRS spectrum. The CO ice line (30 K) lies at 740 au from WD 2226-210 (T_eff = 110,000 K, L_WD = 67 L_⊙), well beyond the proposed location of the dust emission. Furthermore the dissociation timescale for CO in this environment would be very short, such that the amount of CO which could be produced by cometary outgassing, given the observed dust mass of 0.13 M_⊕ (and assuming a volatile composition similar to solar system comets of 0.4%–30%; Mumma & Charnley 2011), would be rapidly photoionized and dissipated before it could build up into detectable levels around the white dwarf.

Remnant planetesimals from the main-sequence system on more eccentric orbits will undergo a collisional cascade; the timescale for this evolution is dependent on the initial mass in planetesimals, the radial location of the planetesimal belt, and their distribution of eccentricities and inclinations. We can estimate the sublimation timescale for a planetesimal being irradiated by the white dwarf. Assuming a typical comet has a radius of 100 km, a perihelion distance of 30 au, and an eccentricity of 0.5, we can calculate the sublimation timescale for a 100 km comet, typical for solar system Edgeworth–Kuiper Belt bodies, to be <1000 yr (see Equation (11) of Stone et al. 2015). This suggests that the observed dust emission is not the product of in situ destruction of remnant planetesimals.

Alternatively, the observed dust emission is the result of cometary bodies on eccentric orbits that are returning to the inner parts of the system after being kicked onto high-eccentricity orbits during the post-main-sequence evolution. A timescale of 10⁵ yr between post-main-sequence evolution and the current state of the Helix Nebula is comparable to the period of comets originating in the solar system's Oort cloud and might point to the potential origin of this material. Given the observed dust mass of 0.13 M_⊕ we can equate this to ≃500 × 10⁶ bodies with a mass equivalent to that of the Hale–Bopp comet, likely representative of an Oort cloud comet or mid-sized Kuiper Belt object. Considering that the vast bulk of cometary bodies is icy matter rather than dust alone, this may be considered a conservative estimate of the number required. Further, assuming that the dust we currently observe has been deposited over the past 10⁵ yr without any loss due to radiation pressure, etc., therefore requires 5000 Hale–Bopp equivalent comets per year being completely disrupted around the Helix Nebula white dwarf. Again, this is a conservative assumption as many cometary bodies will not be completely broken up. The observed dust must be deposited near-instantaneously and not built up over the duration by comets on relatively long periods as blowout will remove dust which is generated mostly at periastron. If the cometary orbits are too short, thereby generating dust more frequently, then the comets will not survive the full 10⁵ yr to produce the excess visible today. The disruption of a massive Kuiper Belt during the post-main-sequence evolution of the system, driving planetesimals to high orbital eccentricities, which are now still entering for only the first time, matches all the observable properties of the system. This leads to the preferred explanation of excess emission originating from a disrupted planetary system.

If these are indeed cometary bodies returning to the inner parts of the Helix Nebula system, detecting the volatiles released by those evaporating comets would be one avenue to test the scenario we have proposed here (e.g., Rebollido et al. 2020). However, this may be challenging against the bright and complex structure of the nebula. JWST medium-resolution spectrometer spectra of the white dwarf and its environment may yield evidence of dust features (e.g., crystalline or amorphous silicate), but there was no evidence for any spectral features in the relatively high signal-to-noise IRS spectrum. However, the higher angular resolution of JWST may help disentangle the circum-white dwarf environment from the nebula. Additionally, Extremely Large Telescope/Mid-infrared E-ELT Imager and Spectrograph N-band imaging of the circum-white dwarf environment could yield a detection of its emission, as the dust is deposited in a fairly well-defined region around the white dwarf.

We have modeled the emission from the vicinity of the Helix Nebula's white dwarf (WD 2226-210) at infrared to millimeter wavelengths. We separate the emission into two physical components associated with the white dwarf for the modeling, based on structure seen in Herschel/PACS observations at far-infrared wavelengths. Upper limits at far-infrared and millimeter wavelengths obtained, respectively, with SOFIA/HAWC+ and ALMA provide constraints on the spatial and emission properties of the dust associated with these physical components.

The observed dust emission is inconsistent with either a remnant of the post-AGB envelope or the in situ collisional destruction, or sublimation, of remnant planetesimals surviving from the main-sequence system. The steep rise at mid-infrared wavelengths and absence of detectable submillimeter emission create a spectral shape of the excess emission inconsistent with the grain sizes expected from dust produced by planetesimals in a collisional equilibrium as would be expected in a classical debris disk scenario. Instead, we conclude that small dust grains must be present and that they need to be replenished at high rates due to the intense radiation field of the central white dwarf.

We propose a scenario where the dust is released by heavy cometary activity around the white dwarf, demonstrating that the longevity of cometary bodies against sublimation necessitates a continual injection of dust by several thousand comets per year. This leads us to conclude that the excess emission observed around the central white dwarf in the Helix Nebula originates from the disruption of a massive Kuiper Belt analog. This would likely have happened when a putative planetary system was destabilized during the star's post-main-sequence evolution.

We thank the referee for their thoughtful review. J.P.M. thanks Alexander Mustill for his insight regarding the evolution of planetesimal belts in post-main-sequence binary systems.

This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. This research has made use of NASA's Astrophysics Data System.

J.P.M. and F.K. acknowledge support by the Ministry of Science and Technology of Taiwan under grant No. MOST107-2119-M-001-031-MY3, and Academia Sinica under grant No. AS-IA-106-M03. J.P.M. acknowledges support by the Ministry of Science and Technology of Taiwan under grant No. MOST109-2112-M-001-036-MY3.

This work was also partly supported by the Spanish program Unidad de Excelencia María de Maeztu CEX2020-001058-M, financed by MCIN/AEI/10.13039/501100011033.

E.V. acknowledges support from the "On the rocks II project" funded by the Spanish Ministerio de Ciencia, Innovaciòn y Universidades under grant PGC2018-101950-B-I00.

D.K. acknowledges the support of the Australian Research Council (ARC) Discovery Early Career Research Award (DECRA; grant No. DE190100813) and the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013.

The work was partially based on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA). SOFIA is jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NNA17BF53C, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 0901 to the University of Stuttgart. Financial support for this work was provided by NASA through award SOF_05-0054_Ertel issued by USRA.

This paper makes use of the following ALMA data sets: ADS/JAO.ALMA#2015.1.00762.S and #2016.1.00608.S. ALMA is a partnership of ESO (representing its member states), the NSF (USA) and NINS (Japan), together with the NRC (Canada) and 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.

This work benefited from the FEARLESS collaboration (FatE and AfteRLife of Evolved solar systems; PI: S. Ertel). The authors particularly acknowledge helpful discussions with Siyi Xu, whose input guided the formative stages of the project, and in obtaining the data.

Facilities: ALMA - Atacama Large Millimeter Array, Herschel - European Space Agency's Herschel space observatory, SOFIA - Stratospheric Observatory For Infrared Astronomy, Spitzer - Spitzer Space Telescope satellite.

Software: This paper has made use of the Python packages astropy (Astropy Collaboration et al. 2013, 2018), NumPy (Harris et al. 2020), matplotlib (Hunter 2007).