FAST search for circumstellar atomic hydrogen. II. Is BD+30$\arcdeg$3639 an interacting planetary nebula?

As the descendant of low- and intermediate-mass stars (0.8 to 8 ${\rm M}_{\sun}$), planetary nebulae (PNe) are important tools for understanding the history of stellar winds and the life cycle of materials in a Galaxy (Kwitter & Henry, 2022). However, extensive PN observations show that the total mass of the ionized material and the central star is only up to 1.5 M_☉ (e.g., Kimura et al., 2012), much lower than the theoretical upper limit. This is the so-called ‘PN missing mass problem’ (Kwok, 1994). CO and H₂ observations indicate that PNe have molecular envelopes with a mass from $10^{-2}$ to a few M_☉ (Huggins et al., 1996; De Marco et al., 2022), while infrared observations suggest that PNe have low dust-to-gas mass ratios ($\sim 0.01$, Phillips, 2007; Andriantsaralaza et al., 2020), suggesting that the missing mass of PNe is not in either molecules or dust grains. Furthermore, only a modest number of PNe have detected so-called ‘AGB’ haloes where at least some of this missing mass is located (see, e.g., Corradi et al., 2003), so the overall problem remains. Another potentially substantial nebular component is atomic hydrogen (H I), especially in ionization-bound PNe (see, e.g., Kwitter & Henry, 2022, for a review).

H I has been detected through the 21 cm emission in several circumstellar shells surrounding evolved stars ($\geq$1 pc; Gérard & Le Bertre, 2006), demonstrating that the H I emission could serve as a tracer of the early history of PN mass loss. However, such observations are very challenging, in part because of strong contamination from interstellar H I emission. To our best knowledge, only five PNe have been reported as having a detection of the H I emission line (Taylor & Pottasch, 1987; Taylor et al., 1990; Rodríguez et al., 2000, 2002; Gérard & Le Bertre, 2006). Taking this into account, we launched a project to search for H I in PNe (Ouyang et al., 2022, Paper I hereafter) using the new Five-hundred-meter Aperture Spherical radio Telescope (FAST; Nan et al., 2011), which is the most sensitive telescope at L-band frequencies of relevance here (Jiang et al., 2020).

In Paper I we reported a FAST detection of H I absorption from the young PN IC 4997. This paper is devoted to investigate the dynamic interaction between PN and the interstellar medium (ISM) through the observation of H I emission toward another young (600-800 year kinematic age) PN BD+30°3639 (BD+30 hereafter). Further details and multi-wavelength imagery, including from the Hubble Space Telescope (HST) and near infrared and optical spectroscopy of BD+30 can be found in the Hong Kong/AAO/Strasbourg/H$\alpha$ PNe database (HASH, Parker et al., 2016) as catalogue entry 512. It is known that interaction between PN ejecta and the surrounding ISM can distort the morphology of circumstellar envelopes and cause a very extensive spatial distribution of both PN ejecta and swept up materials (e.g., Villaver et al., 2003; Sabin et al., 2010). Therefore, such a study might help to understand the ‘PN missing mass problem’.

The PN-ISM interaction has been an intensively studied topic (e.g. Gurzadian, 1969). Based on a thin-shell approximation and snowplow model, Smith (1976) found that the central star of interacting PN have a certain displacement with respect to nebular center. The simple analysis of Borkowski et al. (1990) showed that the PN-ISM interaction could lead to an asymmetrical surface brightness of PN. Soker et al. (1991) suggested that the PN-ISM interaction could be summarized into three phrases as briefly outlined below. Initially, when nebular density is greater than that of the surrounding ISM, the PN is in a free-expanding phase, typically at velocities of $\sim 25$–$30$ km s^-1. With this nebular expansion and without further material injection, the PN density naturally decreases and eventually reaches a value comparable to that of the ambient, surrounding ISM. The upstream of the PN shell is compressed so that the surface brightness of the interacting region increases—examples include PNe like PFP1 (HASH ID 2522, Pierce et al., 2004). This happens for already very evolved PN and there are few clear examples in the literature. As the PN density further drops, the PN shell becomes significantly distorted and the central star of the PN (CSPN) appears to be located closer to the upstream edge. Subsequently, the CSPN would move out of the PN shell, while the PN shell itself eventually dissipates and is completely mixed with the ISM. One putative example of this is Sh 2-174 (HASH ID 643, Tweedy & Napiwotzki, 1994) but this object is now classified as ionised ISM in HASH. The 30-year old analysis of Soker et al. (1991) assumed that the PN shell formed before the PN-ISM interaction. However, observations show that an interaction may occur at various stages of PNe evolution (Rauch et al., 2000).

A more comprehensive hydrodynamic simulation of these types of interaction was performed by Wareing et al. (2007), taking into account a wide range of PN moving velocities, mass loss rates, and ISM density. This work shows that the PN-ISM interaction can be divided into four stages: (i) the PN expands freely; (ii) a bow shock is formed in the interacting region, and the density and surface brightness increase in the upstream of CSPN movement direction; (iii) PN geometric center moves downstream with respect to CSPN; (iv) the PN shell is completely stripped and the CSPN moves outside the PN shell. In truth there are very few examples available for study of PNe that have gone through stages (iii) and (iv).

BD+30 ($\alpha_{\rm J2000}=19^{h}34^{m}45.2^{s}$, $\delta_{\rm J2000}=30\arcdeg 30\arcmin 59.0\arcsec$) is one of the brightest young PNe in the solar neighborhood. It was used as one of four PN examples to demonstrate the value of the so-called distance mapping technique (DMT) based on an expansion parallax method (Gómez-Gordillo et al., 2020). They determined the distance to BD+30 of $D=1.3$ kpc, while the distance calculated from Gaia Data Release 3 (Chornay & Walton, 2021) parallax is $1.62\pm 0.08$ kpc. The Surface-Brightness radius relation statistical technique (Frew et al., 2016), arguably the best such distance technique available for the general PN population in the absence of Gaia distances, provides a distance of $D=2.22\pm 0.63$ kpc, in agreement with Gaia values within the errors. We adopt $D=1.5$ kpc in this paper.

The radial velocity of the PN with respect to the local standard of rest (LSR) is $V_{\rm LSR}=-13.0\pm 0.8$ km s^-1 (Schneider et al., 1983). BD+30 is a relatively nearby, very high surface brightness PNe saturated in IPHAS (INT Photometric H$\alpha$ Survey of the Northern Galactic Plane), PanSTARRS, and $\mathit{WISE}$ (Wide-field Infrared Survey Explorer) imagery. There is clear detection of radio and extended X-ray emission (Leahy et al., 2000). It has a Wolf-Rayet CSPN [WC9] spectrum (Mendez et al., 1991) indicating close-in, fast winds. The optical morphology of BD+30 taken from the available HST imagery is as a boxy-ellipsoid. It is classified as morphological class Eamrs according to the ERBIAS-sparm classification scheme of HASH (refer to Parker, 2022, for details). The optical luminosity is dominated by the bright CSPN (Li et al., 2002).

H₂ was detected in the outer layer of the optical ring, which has a diffuse and highly uneven spatial distribution (Shupe et al., 1998). CO is approximately symmetrically distributed at the north and south ends (Bachiller et al., 2000). We refer the reader to Freeman & Kastner (2016) for a 3D multi-wavelength reconstruction of BD+30. Taylor et al. (1990) reported the discovery of H I emission from BD+30 using the Westerbork Synthesis Radio Telescope (WSRT), which appears to be spatially unresolved and has a minor southeastward extension. Among the five PNe with detected H I emission, BD+30 can be reached by FAST at the most optimum zenith angle. BD+30 has a moderate Galactic disk elevation, which compromises between minimizing the contamination from the ISM H I emission and maximizing the ISM density for the occurrence of substantial PN-ISM interaction. Interferometric observations generally lack zero-spacing data, so emission sources above a certain angular size are resolved out and undetectable. This problem becomes more serious when considering issues such as the missing mass problem in the PN and the interaction between the PN and the ISM. By using a single-dish telescope with a large aperture, such as FAST, it is expected that radio emissions missed by conventional interferometry will be detected.

The observations were performed in the tracking mode of FAST with the 19-beam (M01–M19, see Jiang et al., 2020, for the arrangement of the beams) receiver on 2021 August 9. The integration time was 40.5 minutes. The central beam (M01), which has a half-power beamwidth (HPBW) of about 2.82$\arcmin$ at 1420 MHz, was consistently oriented toward the target during the observation. The standard deviation of the pointing accuracy is 7.9$\arcsec$. The 19-beam observation mode enables us to sample the sky over a region of $20\arcmin\times 20\arcmin$ around the target at 1.42 GHz. The L-band (1.05–1.45 GHz) spectra were recorded by a backend with 1048576 channels, resulting in a frequency resolution of 476.84 Hz and a velocity resolution of 0.01 km s^-1. For the flux calibration, a high intensity noise of about 12 K is injected periodically. The antenna temperature was obtained by

where $T_{\rm cal}$ is the injected noise diode temperature, $P^{\rm cal}_{\rm on}$ and $P^{\rm cal}_{\rm off}$ are the power values when the diode is on and off, respectively. The antenna temperature was converted to flux density by $S_{\rm A}=T_{\rm A}/G$, where $G$ is the gain given by $G=\eta*25.6$ K Jy^-1. The aperture efficiency $\eta$ depends on the zenith angle at the time of the observation (see Jiang et al., 2020, for details). For our observations, the $G$ values of the 19 beams range from 13.3 to 15.9 K Jy^-1. The linear polarization channels XX and YY were separately calibrated and a baseline was subtracted. Standing waves may be generated between the dish and the feed cabin, complicating the baseline (Jiang et al., 2020). Therefore, we used a sinusoidal function to fit the baseline. Nevertheless, the standing wave has a low amplitude (see Figure 1 for an example), and barely affects our results. The XX and YY spectra then were co-added to obtain the final spectra. Based on the line-free spectral regions, we obtained a root mean square noise of the final spectrum of 1.05 mJy.

Figure 2 displays the spectra of the 19 beams. The overall spectral shapes are quite similar among the different beams, suggesting that the H I spectra are dominated by interstellar emission. The observed lines are clearly grouped into two main velocity ranges, $-110$ to $-50$ km s^-1 and $-40$ to $+40$ km s^-1. The first negative velocity group (in terms of LSR) has a sharp blue wing and a broad red tail down to low flux levels, while the latter is much stronger and does not exhibit a broad red tail. Both major velocity groups are composed of several sub-peaks with different relative intensities for different beams.

Assuming that our tentative detection of the H I absorption (Fig. 3c) is real and using the same method described in Paper I, we can estimate the expansion velocity, mass, and column density of the neutral shell traced by the absorption at $-13.29$ km s^-1. The expansion velocity along the line of sight is derived from $V_{\rm LRS}-V_{\rm m}+{\rm FWHM}/2=5.5$ km s^-1. An optical spectroscopic study has shown that the N⁺ and O⁺⁺ ionized regions have an expansion velocity of 28 km s^-1 and 35.5 km s^-1, respectively (Bryce & Mellema, 1999). It is well known that BD+30 contains a central hot-bubble (e.g. Heller et al., 2018) also revealed as extended X-ray emission (Leahy et al., 2000). The ionization and shock heating may results in a radial pressure gradient, causing an acceleration of the ionized gas. Another potential cause for the low H I velocity is that the proposed interaction with the ISM has caused significant deceleration of a much more extended neutral shell, indicating the presence of a stalled astrosphere. Under a thin-shell assumption, the angular radius of the H I envelope is about equal to that of the ionized region ($4^{\prime\prime}$.93, Bojičić et al., 2021). The excitation temperature of H I in the PN is poorly known, but it is reasonable to assume an intermediate value of 100 K (see Altschuler et al., 1986, and Paper I for a discussion). From the absorption feature at $-13.29$ km s^-1, the column density and mass of the H I are thus determined to be N${}_{\rm abs1}=4.8\times 10^{20}$ cm^-2 and M${}_{\rm abs1}=3.1\times 10^{-2}$ M_☉, respectively. From the continuum flux density we estimate the mass of ionized gas to be $3.7\times 10^{-2}$ M_☉. Consequently, the mass ratio of atomic to ionized gas is about 0.8, which is close to the value in another young PN IC 4997 (Paper I). It should be emphasized that the atomic/ionized gas mass ratio is independent of the distance and nebular size under the thin-shell assumption. The number of neutral hydrogen atoms in BD+30 is estimated as $4\times 10^{55}$.

Figure 5 presents a comparison between our spectrum and the [N II]$\lambda 6548$ emission line observed by Bryce & Mellema (1999). The [N II] line traces the velocity of the outermost ionized gas, which can be accelerated by the pressure from the ionized gas. It exhibits two peaks arising from the front and back sides of the expanding ionized shell; the blue one has a velocity coinciding with the dip seen in our H I spectrum between $-40$ and $-50$ km s^-1. This leads us to the speculation that this dip is associated a H I absorption tracing the outer boundary of the ionized nebula, although the identification is not certain because the velocity is unfortunately similar to that of the interstellar H I. A further support for this speculation comes from the wide-aperture observation of Neiner et al. (2000) who show that [O I] and [S II] emission lines have a similar width suggesting that the velocity remains similar across the ionization front. However, even if validated, the absorption over the velocity range from $-40$ to $-50$ km s^-1 has an equivalent width comparable to that of the absorption at $-13.29$ km s^-1, suggesting the mass of the absorption gas to be M${}_{\rm abs2}\sim 10^{-2}$ M_☉, still far less than the missing PN mass.

According to Aller & Hyung (1995), BD+30 has a very low He/H abundance ratio ($\sim 0.02$). Even correcting for helium abundance, the total gas mass is still fairly small compared to the upper theoretical limit of PN mass. Therefore, the observation of H I absorption as in BD+30 cannot solve the ‘PN missing mass problem’. The uncertainties of the mass estimation primarily stem from the assumed radius of the atomic shell. Nevertheless, the observations of Bachiller et al. (2000) demonstrate that the molecular and ionized regions of BD+30 have a similar size, which justifies the thin-shell assumption for the calculation of the H I mass.

Is the missing mass contained in extended nebulosity that has been stripped by the ram pressure of the ISM? To settle this question we need to investigate the interaction between BD+30 and the ISM. BD+30 has a proper motion of $-2.3\pm 0.03$ mas yr^-1 and $-8.8\pm 0.03$ mas yr^-1 in right ascension and declination, respectively (Gaia Collaboration et al., 2021). The values are in good agreement with those given by Kerber et al. (2008) but with lower uncertainties. This corresponds to a velocity of $14.1$ km s^-1 eastward and $62.6$ km s^-1 northward at an adopted distance of 1.5 kpc, and thus the position angle of the proper motion is 192^∘.7 east of north. The heliocentric velocity of BD+30 is $-$31.4 km s^-1 (Schneider et al., 1983). Assuming that the ISM follows the Galactic rotation with no peculiar motion, we can calculate the velocity of BD+30 relative to its ambient ISM ($V_{*}$). Using the Oort constants of $A=15.3$ km s^-1 kpc^-1 and $B=-11.9$ km s^-1 kpc^-1, we obtain the radial and tangential velocities of the local ISM with respect to the position of the sun to be $V_{r}=17.7$ km s^-1 and $V_{t}=-32.4$ km s^-1, respectively. The ISM velocities can be converted to heliocentric values by adopting a solar motion with respect to the LSR of $(U,V,W)_{\odot}=(10.6,10.7,7.6)$ km s^-1 (Reid et al., 2019), where the velocity components $U$, $V$, and $W$ are respectively defined along the radial, orbital, and vertical directions of the Galactic plane. We finally have a $V_{*}$ value of $\sim 50$ km s^-1. The Galactic disk elevation of BD+30 is about 130 pc, where the ISM is largely composed of cold neutral gas (Heiles & Troland, 2003). The large motion velocity and the high ISM density provide favorable conditions for the interaction between stellar wind and ISM, even for such a young PN. However, there is scant evidence in the optical for any obvious ISM-PN interaction probably because of insufficient ultraviolet radiation to excite optical lines in extended regions.

The statistical study of Ali et al. (2012) shows that PNe lying close to the Galactic plane are most likely to interact with molecular and cold neutral clouds. Matthews et al. (2015) have demonstrated that a H I observation with a single-dish telescope is more sensitive to detect the diffuse low surface brightness emission surrounding the carbon-rich asymptotic giant branch (AGB) star IRC+10216 compared to the Very Large Array. A H I emission tail formed by the interaction with the ISM has been discovered in an AGB Mira variable (Matthews et al., 2008). It has been shown that the asymmetric distribution of H I emission, the increase of flux in the interaction region, and the formation of bow shock can be perceived as a sign of wind-ISM interaction (Tweedy & Kwitter, 1996; Dgani & Soker, 1998). In Figure 6 we examine the spatial distribution of the integrated intensity of C6 and C7. It is shown that the outer beams detect stronger H I emission than the central beam. The H I emission is obviously enhanced toward the southwest, which is roughly aligned with the direction of the proper motion. This is a typical characteristic of the interaction between a moving PN and the ISM that causes the accumulation of gas upstream of the direction of PN motion. Apart from that, the density gradient of the ISM could also affect the morphology of the circumstellar envelope, and possibly lead to an enhancement of H I emission toward the Galactic disk. This is not supported by our observations (see Figure 6) partly because the density scale height is much larger than the nebular size.

How much is the ISM gas that could be swept up by the motion and expansion of BD+30? Assuming that the stellar winds were launched during the AGB phase and lasted $10^{5}$ years and taking the expansion velocity to be a typical value $V_{\rm exp}=20$ km s^-1, we find that these earlier stellar winds, prior to final PN envelope ejection, can reach up to a radius of 2 pc in the head-on direction, corresponding to an angular radius of $\sim 5$ $\arcmin$. Considering the significant uncertainty of such a crude estimation, this value is not far from the observed extent of H I emission (Figure 6). If BD+30 crosses the ISM at a constant velocity of $V_{*}=50$ km s^-1, its central star has travelled a length of about 5 pc since it joined the AGB phase. A runaway shell expanding at a constant velocity ($V_{\rm exp}<V_{*}$) would take on a raindrop-like shape. A simple geometric calculation shows that the volume of the plowed ISM gas is about 42 pc³. The cold neutral medium has a density of $n_{\rm H}=10$–100 cm³ (Wooden et al., 2004), and thus the mass of the swept ISM gas is estimated to be 10–100 ${\rm M}_{\sun}$. If the central star of BD+30 has been decelerated by dynamic friction (i.e. the initial velocity is larger than 50 km s^-1), the mass of the swept ISM gas will be larger than the above estimated value. Therefore, the mass value estimated for the swept material is compatible with that derived from the intensities of the H I emission lines ($\sim 100$ M_☉). Despite the severe uncertainties in these estimates, there is a reasonable agreement that supports the swept-material interpretation. It should also be noted that the flow surrounding the PN-ISM interface could constrain the nebular expansion and cause a smaller mass of the swept ISM gas. A sophisticated modelling is required to obtain an accurate estimation, which is beyond the scope of this paper. According to Borkowski et al. (1990), the PN-ISM interaction is observable when the nebular density decreases to a critical density, $n_{\rm crit}=[(V_{*}+V_{\rm exp})/{c_{s}}]^{2}n_{\rm H}$, where $c_{s}$ is the isothermal sound speed. Assuming $n_{\rm H}=10$ cm^-3 and $c_{s}=10$ km s^-1, we obtain $n_{\rm crit}=490$ cm^-3. The density of the H I shell is presumably lower than $n_{\rm crit}$, and thus could ensure the occurrence of the PN-ISM interaction.

Also shown in Figure 6 is the contour map of the WISE W4 (22 $\mu$m) band. This appears unresolved, as seen in the WISE images in HASH for convenience. It has been shown that thermal emission from cold dust in the mid-infrared is a good tracer for extended PN-ISM structures (e.g. Zhang et al., 2012; Ramos-Larios et al., 2018). However, we do not see any dust emission associated with the H I shell. To further validate our observations, we extracted the data of the Galactic Arecibo L-Band Feed Array H I (GALFA-H I) survey (Peek et al., 2011, 2018). Although the GALFA-H I data have a lower sensitivity than the FAST observations, they cover a contiguous sky area surrounding BD+30. The channel maps of the H I emission is shown in Figure 7. We can see that the maps in the velocity range of $-18$ to $-8$ km s^-1 exhibit a cavity in the position of BD+30 and the emission appears to be stronger in the leading edge of the PN. This is consistent with our observations and may be considered as a possible trail that has be plowed by a moving expending shell in the ISM. Due to the large velocity shear at the PN-ISM interface, a Kelvin-Helmholtz instability is developed and induces turbulent mixing in the wake. If our interpretation is correct, it would mean that the H I shell of the PN has largely merged with the ISM and thus is the key to solve the puzzle of missing mass.

Despite the points mentioned above favoring the PN-ISM interaction, the possibility of the H I features resulting from contamination of a foreground cloud remains. The large swept ISM mass is extraordinary for a PN, but is typical for supernova remnants (SNRs) and Wolf-Rayet stars. We notice that indeed a SNR, G65.2+5.7 ($\alpha_{\rm J2000}=19^{h}33^{m}$, $\delta_{\rm J2000}=30\arcdeg 30\arcmin$), is coincident in position with BD+30. G65.2+5.7 has a diameter of $\sim 4\arcdeg$ and a distance of 800 pc (Xiao et al., 2009), and thus is likely to contribute to the H I absorption and emission. Additional support for this possibility comes from Vincent & Kang (2011) who discovered emission from an expanding H I shell around G65.2+5.7 within the velocity range of $-12$ km s^-1 to 10 km s^-1. This is consistent with the components C6 and C7, and these could be part of the SNR shell. The H I emission found towards BD+30 by Taylor et al. (1990) also falls within this velocity range, but is likely associated with BD+30 as it is spatially unresolved and coincident with BD+30. The BD+30 emission measured by Taylor et al. (1990) may be related to the spatially extended component detected by FAST, supporting some association between the wide H I emission and the PN and the idea that the majority of the H I emission would be swept ISM material. The presence of large scale emission from the SNR complicates the interpretation for velocities redward of the systemic velocity of BD+30.

We have searched for atomic hydrogen surrounding the young PN BD+30 using the FAST 19-beam observations. The new high-quality spectra confirm the existence of H I emission toward this PN, as previously discovered by interferometry (Taylor et al., 1990). However, we find that the H I emission detected by Taylor et al. (1990) is in fact resolved into two individual velocity components and extends to an angular radius of $\sim 12\arcmin$ rather than confined to a compact region. The mass of the H I shell is larger than 100 M_☉, excluding the possibility that it is primarily composed of ejected stellar materials from the host star of BD+30. A reasonable conjecture is that the large H I shell stems from material swept up by the BD+30-ISM interaction processes at earlier stages of evolution although such a large swept ISM mass has never been reported for a PN. The H I emission is enhanced roughly along the leading edge of BD+30, providing further evidence supporting this conjecture.

We tentatively detect circumstellar H I absorption features. The H I absorption is expected at various velocities: systemic if it is located in a stalled asterosphere where it sweeps up interstellar gas or at the blueshifted velocities if just outside the ionized region. The H I gas producing the absorption is estimated to have a mass of $\sim 3\times 10^{-2}$ M_☉.

However, we cannot firmly rule out the possibility that the detected H I absorption and emission may originate from a SNR that overlaps along the line of sight and happens to have similar velocities to BD+30. If this was the case, the number of PNe with detected H I emission would be reduced from five to four. Nevertheless, this possibility is not particularly likely because it does not seem to fit the interferometric observations of Taylor et al. (1990).

Investigating PN-ISM interactions in cases like BD+30 is a promising avenue to solve the ‘PN missing mass problem’. The resultant circumstellar structures are prone to be revealed by the observations of H I emission. Unlike interferometric observations that are insensitive to structures larger than a certain extent, single-dish telescopes have noteworthy benefits in observing extended circumstellar materials. Although the results of this study are somewhat inconclusive, we demonstrate the feasibility of utilizing the H I observations of FAST to efficiently reveal the structures for future PN-ISM interaction studies.