FAST Search for Circumstellar Atomic Hydrogen. I. The Young Planetary Nebula IC 4997

Using the Five-hundred-meter Aperture Spherical radio Telescope in Guizhou, China, we detect the 21 cm neutral atomic hydrogen absorption in the young planetary nebula IC 4997. The absorption arises from a shell that is also associated with Na i D lines. The H i shell has a mass of 1.46 × 10−2 M ⊙ and a dynamic age of 990 yr. The column density of H i is estimated to be 7.1 × 1020 cm−2, which can be well explained in terms of a photodissociation region around the ionized nebula, limited by the self-shielding of H2. We find that the atomic-to-ionized hydrogen ratio is 0.6, suggesting that H i substantially contributes to the overall nebular mass.


INTRODUCTION
Planetary nebulae (PNe) are the ionized, ejected remnants of the stellar winds from stars of low to intermediate mass (0.8 to 8 M ) at the end of their lives.They serve as important tools for understanding stellar evolution and the life cycle of dispersed materials in a Galaxy.However, observations show that the total mass of the central star and ionized nebula of up to 1.5 M (Kimura et al. 2012) is significantly lower than the theoretical upper limit, indicating not all of the mass ejected by the star is seen in the PN.Indeed, the total mass usually estimated for the main ionised PN shell itself is only ∼ 0.1 M (see, e.g., Villaver et al. 2003).The connection between the birth mass of PN central stars and the mass left when they die (the initial-final mass relation) is thus poorly known, (leading to the so-called 'PN missing mass problem', see Kwok 1994).Such gaps in our knowledge on the abundant PN progenitor stars, which dominate the recycling of stellar material, severely impact our ability to model the evolution and chemical enrichment of any galaxy.The asymptotic giant branch (AGB) halos (see, e.g., Corradi et al. 2003) are thought to be the reservoirs for much of this missing mass.Combining the ionized gas masses derived from radio observations with the calculated molecular masses, Huggins et al. (1996) obtained nebular masses ranging from 10 −3 to 1 M and found that the median value of molecular mass is 0.031 M in 44 PNe.Many PNe are ionization bounded, where part of the matter lost by the stellar wind during the pre-PN stage has not yet been ionized by the central star.A substantial part of their envelope may remain neutral (see, e.g., Kwitter & Henry 2022, for a review).These outer regions of PNe show dust emission, but the abundance of refractory elements limits the dust mass to no more than 1% of the total mass of the nebula.Infrared observations suggest a lower dust-to-gas mass (Phillips 2007;Andriantsaralaza et al. 2020), and hence most of the missing PN mass cannot be in dust grains.More recently, studies of molecules and ions in PNe which have evolved from common envelopes of binary stars have also shown that the missing mass is not present in these forms in these objects (Sarkar & Sahai 2021;Santander-García et al. 2022).
However, there remains a vital circumstellar component that has not been sufficiently investigated and that could significantly contribute to the overall mass-loss budget: atomic hydrogen.
The presence of atomic hydrogen can be indicated through the spin-flip transition at 21 cm wavelength in the radio regime.Unfortunately, the detection of circumstellar atomic hydrogen is hampered by the omnipresence of 21 cm emission from the ambient interstellar medium.Rodríguez & Moran (1982) were the first to detect H I absorption in the young PN NGC 6302 using the Very Large Array after earlier unsuccessful attempts (Thompson & Colvin 1970;Zuckerman et al. 1980), and soon thereafter, Altschuler et al. (1986, A86 hereafter) detected circumstellar H I absorption towards IC 4997 using the Arecibo telescope.Subsequent investigations confirmed the presence of circumstellar H I absorption in five PNe and circumstellar H I emission in one PN (Taylor et al. 1990).To date, definite circumstellar H I absorption has been detected in nine PNe (Rodríguez & Moran 1982;Altschuler et al. 1986;Schneider et al. 1987;Taylor et al. 1990), and five PNe have been detected with circumstellar H I emission (Taylor & Pottasch 1987;Taylor et al. 1990;Rodríguez et al. 2000Rodríguez et al. , 2002;;Gérard & Le Bertre 2006).
The number of PNe with circumstellar H I detections is extreme small, compared to the total number of verified Galactic PNe (∼3800; e.g.Acker et al. 1992;Kohoutek 2001;Drew et al. 2005;Parker et al. 2006;Miszalski et al. 2008;Parker et al. 2016).To improve our knowledge on the mass and dynamics of the atomic gas budget in PNe, we have begun a program to search for H I feature arising from PNe using the Five-hundred-meter Aperture Spherical radio Telescope (FAST; Nan et al. 2011), the most sensitive single dish telescope in L-band available today (Jiang et al. 2020).In this paper, we report our first results for IC 4997 (l = 58 • .33,b = −10 • .98), a compact, young, rapidly evolving, well-studied, high-surface brightness PN with a prominent central star.

OBSERVATION AND DATA REDUCTION
The observations were performed using the tracking mode of the FAST 19-beam receiver (see the upper panel of Figure 1) on August 10, 2021.The half-power beamwidth (HPBW) of the center beam is about 2.82 at 1420 MHz.The HPBWs of the outer beams is slightly larger than that of the central beam, but the deviation is not more than 0.2 .The standard deviation of pointing accuracy is within 7.9 (Jiang et al. 2020).The central beam (M01) was pointed at the target.As IC 4997 has a size smaller than the beam size, the other beams sample the adjacent off-source positions.The backend of the spec (F) was used to record the spectral line.This mode records the full frequency range (1.05-1.45GHz) of the L-band observation of FAST with 1048576 channels, giving a frequency resolution of 476.84 Hz or a velocity resolution of 0.10 km s −1 .A high intensity noise of about 12 K was injected periodically during about 50 minutes observation times for the flux calibration.The integration time was 3000 seconds.
The frequency-time 'waterfall' plots were examined to remove the radio interference that occasionally appears in the spectra recorded for all the beam.The integrated spectra were then obtained for each beam, as shown in the lower panel of Figure 1.Considering the potential ripple effect, we fitted the baselines using a sinusoidal function, although a linear fitting does not make any practical difference.The interstellar 21 cm emission in the spectra has at least two discrete velocity components, which were not resolved in the observations performed by A86 due to their much poorer spectral resolution.We integrated the spectra over three different time spans.For each beam, there is practically no difference in the different integrated spectra, demonstrating the stability of the observations.
In order to investigate the spatial distributions of interstellar 21 cm emission, we derived the flux density integrated over the velocity range from −50 km s −1 to 80 km s −1 in the local standard of rest (LSR) for each beam.The upper panel of Figure 1 shows that the spatial variation of the interstellar 21 cm line radiation is smooth around IC 4997, providing the potential to reliably subtract the interstellar contamination from the on-source spectrum.The average of the spectra obtained with the six beams (M02-M07) surrounding the central beams was taken as the off-source spectrum.Considering the potential ripple effect, we fit the the baseline of the 'ON minus OFF' spectrum using a sinusoidal function although a linear fitting does not make any practical difference.The 'ON minus OFF' spectrum is shown in the upper panel of Figure 2, which has been binned over ten adjacent frequency channels to improve the signal-to-noise ratio, resulting in a spectral resolution of 1 km s −1 .

RESULTS
IC 4997 has a V LSR value of −49.8 ± 1.3 km s −1 (Schneider et al. 1983).This velocity is marked in the upper panel of Figure 2, which clearly reveals a dip just blueward.We do not see a circumstellar H I emission feature, which if present should have a LSR velocity of −49.8 km s −1 .The subtraction of the off-source spectrum is imperfect redward of the target velocity because of the strong interstellar emission.After closely examining the waterfall plots (Figure 3), we confirm that the dip is stably present in the data recorded by beam M01, and is not detected in the other beams.To further investigate the circumstellar origin of the dip, we derive two OFF-source spectra by averaging the spectra obtained for the west and east beams (M02 and M05) and those for the north and south beams (M03, M04, M06, and M07), respectively.The difference between the two OFF-source spectra is shown in the middle panel of Figure 2, where the dip is invisible and the baseline nearby the dip is flat.We compute the deviation (σ v ) of the individual spectra obtained with beams M02-M07 with respect to their average (i.e. the off-source spectrum described in Section 2), as shown in the lower panel of Figure 2. σ v has a value of ∼ 0.8 mJy at the velocity of the dip, which is much smaller than the peak depth of the dip, and can be taken as the flux uncertainty introduced by the subtraction of OFF-source spectrum.We thus unambiguously confirm the conclusion of A86 that the dip originates from the 21 cm absorption on the near side of the expanding H I envelope around IC 4997.The σ v value however is relatively large at −50 < V LSR < 80 km s −1 so that features at this velocity range shown in the upper panel of Figure 2 are more likely to be spurious.We can see σ v ∼ 2 mJy at v = −49.8km s −1 .This in part explains the nondetection of the circumstellar H I emission feature, whose peak flux density has an upper limit of 3σ v = 6 mJy.
The high-resolution and high-sensitivity FAST spectrum allows us to precisely estimate the parameters of the circumstellar 21 cm absorption.By fitting a Voigt function to the absorption feature, we obtain the LSR velocity of the peak absorption V HI = −63.27± 0.07 km s −1 , the full width at half maximum (FWHM) of 12.74 ± 0.11 km s −1 , the maximum absorption flux density of −8.07 ± 0.07 mJy, and the integrated flux density of −136.7 ± 1.5 mJy km s −1 .As pointed out by A86, such a V HI value together with the fact that IC 4997 is in the first quadrant of the Galactic plane implies that the absorption is unlikely to arise from a foreground cloud situated along the line of sight, because otherwise the distance to this PN or the velocity of the cloud would be unrealistic high.Assuming the continuum flux density at 21 cm to be 37 ± 2 mJy (Isaacman 1984), we deduce the integration of the optical depth over the velocity τ dv = 3.70 ± 0.06 km s −1 , which is about 1.8 times higher than that derived by A86.
The measurement of the FWHM places a loose upper limit on the gas temperature as T < 3600 K. Given the relatively low critical density (< 1 cm −3 ) for the H I transition, T is expected to be close to the spin-excitation temperature (T ex ).As discussed by A86, T ex is unlikely to be smaller than 50 K.To facilitate the comparison of results, we follow A86 and simply adopt T ≡ T ex = 100 K in the following calculations.After deconvolution for the thermal broadening, the H I line width caused by nebular dynamics is determined to be ∆V HI = (FWHM 2 − 5.5kT /m H ) 0.5 = 12.6 km s −1 , where m H is the mass of the hydrogen atom, and k is the Boltzmann constant.If ignoring the projection effect (i.e.assuming a plane-parallel slab) and turbulent broadening and assuming that the velocity increases with radius, ∆V HI represents the velocity difference between the inner and outer H I layers.The maxima expansion velocity of the H I shell is approximately determined through V exp (HI) = V LSR − V HI + ∆V HI /2 = 19.8km s −1 , which is higher than the expansion velocity of the inner ionized region (∼ 14.5 km s −1 ; Sabbadin 1984).If assuming T = 1000 K, ∆V HI and V exp (HI) would be scaled down by a factor of 1.16 and 1.05, respectively.Rao et al. (2020) compare the heliocentric radial velocities of several atomic and molecular lines, (see their Figure 6).For the sake of comparison, we convert V HI into the heliocentric radial velocity and obtain −80.70 km s −1 , which is intermediate between the heliocentric radial velocities of [Fe II] and H 2 emission lines.The Na I D lines show two absorption components at −84 and −97 km s −1 with the former stronger than the latter (Rao et al. 2020).Our H I measurement coincides with the velocity of the stronger Na D component, while the weaker velocity component is not detected in the H I absorption spectrum.From the report of Rao et al. (2020), the stronger Na D component appears to have a width of ∆V NaI ∼ 8 km s −1 .Considering the projection effect (see Figure 4 and the next section), ∆V NaI may more accurately represent the velocity difference between the inner and outer neutral layer compared to ∆V HI .
Following the approach of A86, we estimate the mass of circumstellar atomic hydrogen by where Γ is a geometric factor to calibrate the total mass of H I in the envelope relative to the gas detected in absorption, D is the distance to IC 4997, and θ is the angular radius of the H I envelope.Note that the mass of atomic hydrogen at the far side of the nebular is already included in this equation through a spherical symmetry assumption.If the atomic hydrogen is located within a thin shell, we could assume Γ ≈ 1, or otherwise M HI would be underestimated.
The uncertainty is dominated by the distance measurement.Various distances to IC 4997 have been reported in the literature.Recently, Frew et al. (2016) report a value of 4.85 kpc from their statistical distance estimator based on a robust surface brightness-radius relation, which is more than twice as large as that used in A86.An even larger distance of 7.12 +3.83 −1.85 kpc is reported by Chornay & Walton (2021) based on the Gaia's Early Data Release 3 (Gaia Collaboration et al. 2021).We thus adopt D = 5 kpc, as a compromise between these recently reported distances and given the large errors.The angular radius θ is assumed to be 0.86 (Bojičić et al. 2021), which is only slightly larger than the value used in A86 (0.8 ).Thus we estimate the radius of the atomic hydrogen envelope to be r(HI) = 0.02 pc, which indicates a dynamic age of t(HI) = r(HI)/V exp (HI) = 990 yr.
From the above we obtain M HI = 1.46 × 10 −2 M .This is nearly 7.7 times higher than that derived by A86.Even after correcting for the effects of different distances and θ values used by us and A86, our calculations still give a H I mass that is 1.8 times higher.Using the formula given by Schneider et al. (1987) but assuming the same D and θ values as above, we derive the mass of ionized hydrogen M HII = 2.46 × 10 −2 M .As a result, the mass ratio between atomic and ionized gas (M HI /M HII ) is about 0.6, suggesting that atomic hydrogen has a significant contribution on nebular mass.Note that the M HI /M HII ratio is distance independent.
Under the optically thin assumption, the H I emission line can be used to determine M HI by where F ν dv is the integrated flux density of the H I emission line.Assuming that the H I emission line has a Gaussian profile with the same FWHM of the absorption line and based on the measured upper limit of the peak flux density, we have F ν dv < 81.5 mJy km s −1 , resulting in M HI < 0.48 M .This is fully compatible with that obtained using the absorption line.Comparing this with Equations 1 shows T ex < 1950 K.
If assuming that the atomic hydrogen is primarily injected from the central star, we can estimate the mass loss rate ( Ṁ ).As observations clearly reveal an outwardly accelerating stellar wind (Rao et al. 2020), we simply assume a constant acceleration v ∝ t.The dynamic age t(HI) then is half of the age of the H I shell.Taking V exp as the velocity of the outer H I layer, we roughly estimate the duration time of H I injection to be ∆t(HI) = 2t(HI)∆V NaI /V exp = 800 yr.The thickness of the H I shell is estimated to be ∆r(HI) = r(HI)[4t 2 (HI) − (2t(HI) − ∆t(HI)) 2 ]/4t 2 (HI) = 0.01 pc.Consequently, the mass loss rate of ṀHI = M HI /∆t(HI) is derived to be 2 × 10 −5 M yr −1 .After correcting for the helium abundance (Ruiz-Escobedo & Peña 2022), we obtain Ṁ = 3 × 10 −5 M yr −1 .It should be noted that the value is accurate only to an order of magnitude since the nebular dynamics are presumably more complex than is assumed here.The derived Ṁ is a typical value for a PN during the 'superwind' epoch at the tip of the AGB (Decin et al. 2019).If the stellar ejecta are largely molecular, which are subsequently dissociated by UV irradiation to enhance the abundance of circumstellar H I, the obtained Ṁ value is still valid, but the duration time is that of injecting the H I precursor (H 2 ).This possibility will be discussed in the following section.

DISCUSSION
The total number of H atoms in the atomic shell is 1.7 × 10 55 .Under the thin-shell approximation at r(HI) = 0.02 pc this gives a column density of N HI = 7.1 × 10 20 cm −2 .H I can be formed by photodissociation of the H 2 molecules in the original stellar wind, by photons in the 912-1100 Å range.The dissociation is limited by self shielding of H 2 where the absorption lines become optically thick, and by dust shielding where the far-UV photons suffer heavy extinction.Using the formula provided by Sternberg et al. (2014) and assuming a solar metallicity, a gas density of n = 100 cm −3 , a standard interstellar radiation field (Draine 1978), we find that the self shielding of H 2 yields a N HI value of 6.0 × 10 20 cm −2 .However, the radiation field near a hot central star of a PN is considerably stronger than the standard interstellar radiation field.The column density to the H 2 dissociation front scales as χ/n HI where χ is the far-UV flux relative to the standard field.For a stellar blackbody temperature of T eff = 55, 000 K (Miranda et al. 2022) and an assumed luminosity of 7 × 10 3 L , χ is of order 5 × 10 3 .To obtain the measured N HI requires a density of n HI ∼ 6 × 10 5 cm −3 .For comparison, the electron density of the ionized region varies from n e ∼ 3 × 10 4 to ∼ 10 6 cm −3 (Hyung et al. 1994;Ruiz-Escobedo & Peña 2022), encompassing the derived density of the H I region.The H I shell thickness is roughly δr(HI) ≈ N HI /n HI = 4 × 10 −4 pc, much lower than r(HI), thus justifying the thin-shell approximation.Note that here we have assumed a molecular instead of an atomic wind as in Section 3, and δr(HI) represents the mean free path of UV photons.Thus it is not surprising to see δr(HI) < ∆r(HI).
Dust shielding becomes important for an extinction of A V > 2 mag.The Balmer decrement suggests that the ionized region has an extinction of A V ∼ 0.7 mag (Burlak & Esipov 2010).Under the assumption of typical interstellar dust, the relation between A V and hydrogen column density of Güver & Özel (2009) suggests that the H I shell contributes approximately A V = 0.14 mag to the extinction.However, IC 4997 is a carbon-rich PN with an unusually high dustto-gas ratio (Lenzuni et al. 1989).Therefore, the circumstellar dust may contribute a higher extinction than expected for interstellar silicates.Nevertheless, the agreement with the density of the ionized region suggests that the size of the atomic shell is set by self shielding of H 2 , and that internal extinction plays only a minor role.
Based on this calculation, we predict that under the thin-shell approximation, the atomic layer surrounding a PN will have a N HI which depends mainly on the stellar temperature and luminosity, and on the inner radius and density of the neutral shell: where L * is the stellar luminosity in solar units, T eff is the effective stellar (blackbody) temperature in K, r in is the inner radius of the atomic shell in units of pc, and n HI is the H I density in units of cm −3 .The photodissociation rate is k ∼ 4 × 10 −11 χ s −1 .For χ = 5000, full dissociation therefore takes around 100 years.This is less than the age of the shell, noting that the star will have been hot enough for dissociation to occur during only part of the time since the ejection of the shell.The H I shell is therefore likely in equilibrium with the current radiation field.The H 2 dissociation front is defined by where the dissociation rate (after self-shielding) equals the H 2 formation rate which for our density and assuming normal dust is ∼ 10 −11 s −1 .
The agreement between the derived density of the atomic shell and the ionized region, and that between the H I expansion velocity and the velocity field of the ionized region indicate that the H I shell is closely related to the ionized region: it is a continuous shell encompassing an ionization front and a dissociation front.
IC 4997 is a very bright, well studied PN that is characterized by variability in magnitude, emission line ratios, electron temperature and density, line profile, and radio continuum flux density over a timescale of several decades (see Miranda et al. 2022, and references therein).The variable physical conditions may have led to variable abundance of H I. The difference between our results and those of A86 may be due to a substantial increase of the H I abundance during the past 36 years.An increasing number of UV photons (larger χ) can increase the H I column density over that time period, based on the photodissociation time scale of H 2 calculated above, In contrast, the equivalent widths of the Na I D absorption lines have decreased by more than half between 1989 and 2014 (Rao et al. 2020), suggesting a more complex picture.Rao et al. (2020) show that the Na I D lines trace two shells, one coincident in velocity with H I and the other with a higher expansion velocity.This second shell has a slightly weaker Na I D absorption depth, and is visible in CH as well.These optical lines are seen in absorption against the central star and therefore trace a smaller region than the radio data which is seen in absorption against the ionized nebula.Comparing the Na I D and H I absorption features, the Na column density derived for the velocity component of H I is N NaI = 9.4 × 10 11 cm −2 .This value is not inconsistent with that obtained by Welty et al. (1994).The faster Na I D component, with N NaI = 4.7 × 10 11 cm −2 , has no H I counterpart.This discrepancy is most easily explained if this component is not a full shell but a clump in front of the star which does not cover most of the radio-emitting nebula.The proposed nebular structure is sketched in Figure 4, where the neutral clump and shell generate two Na I D absorption components of comparable intensity, while the H I absorption dominantly originates from the shell.It is possible that there exist more neutral clumps that do not manifest themselves because of the absence of background continuum emission.As shown in Figure 4, the projection effect of the expanding shell may significantly contribute to the H I line broadening, while the Na I D line profile is essentially free from the projection effect.If the radial velocity of the neutral shell is constant and thus the line width is only due to the projection effect, the stronger Na I D component will peak at the blue edge of the H I line.This is clearly inconsistent with the observation, and validates again the radial acceleration model.
The estimation of M HI is significantly affected by the uncertainty of T ex , as one can tell from Equation 1.Based on photoionization modelling using Cloudy (Ferland et al. 2017), we find that the H I region has an electron temperature of 500 and 300 K under the assumption of n = 10 5 and 4 × 10 5 cm −3 , respectively.This suggests that T ex is most likely in the range from 100-500 K.If taking T ex = 500 K, we would obtain M HI = 7.3 × 10 −2 M .Therefore, the H I absorption measurement indicates that the atomic shell accounts for about 37-75% of the mass of the known nebula (ionized plus atomic).It is a significant nebular component.However, the total mass of M HII+HI =(3.9-9.8)×10−2 M is still very low compared to the total mass expected to be lost in a superwind.We suggest that the second component seen in Na I D is due to a clump in the wind, and does not add significant mass.Although the thickness of the H I shell set by H 2 self shielding suggests the presence of molecular material outside of this shell, the lack of a CO detection does not indicate a massive molecular envelope.Lenzuni et al. (1989) obtained a dust-to-gas ratio of 0.02, yielding a dust mass of (7.8-19.6)×10−4 M .Assuming that in reality the dust-to-gas ratio is a more typical 0.01 (limited by the fraction of refractory elements), the undetected molecular mass may contain as much mass as the ionized plus atomic region.That still leaves the total mass of the nebula relatively small.
It follows that, for this object at least, the 'missing mass' does not appear to be in the nebula.

CONCLUSIONS
Circumstellar H I absorption is clearly detected in the radio spectrum of IC 4997, the first such detection from FAST.This feature originates from a shell associated with one of the velocity components of the Na I D lines.Adopting an updated distance of 5 kpc, we obtain M HI = 1.46 × 10 −2 M and N HI = 7.1 × 10 20 cm −2 .The H I shell was ejected at least 990 yr ago during the superwind phase.The neutral hydrogen in this PN is comparable to the ionized nebula in mass, and may have increased over the past 36 years.Nevertheless, the total nebular mass is still far less than that expected for the superwind.This study demonstrates the feasibility of using FAST to search for atomic hydrogen in a large sample of PNe and to provide new insights into the so-called 'PN missing mass problem'.

Figure 1 .Figure 2 .
Figure 1.(a) The layout of the 19 beams (M01-M09) with the central beam (M01) pointed to IC 4997.The color represents the integrated flux density of the interstellar 21 cm emission feature.(b) The spectra obtained with each beam.The scales of all panel are the same.Different colors represent the spectra integrated over different time spans, which actually superpose each other.

Figure 4 .
Figure 4.The schematic sketch (not drawn to scale) of the nebular structure.The central filled circle represents the central star, which is surrounded from the inside out by ionized, neutral, and molecular shells.A neutral clump lies within the molecular shell.The radial arrows indicate the outwardly accelerating expansion.The purple arrows represent VLSR.The regions where the H I and Na I absorption features, as shown in Figure 2(a), are formed, are denoted with red and green, respectively.