© ESO 2022

1 Introduction

The detection of H I with the 21 cm line (²S_1∕2, F = 1–0; f₀ = 1420.4058 MHz) is a milestone in the study of interstellar matter in our Galaxy (Ewen & Purcell 1951). Immediately after, astronomers tried to investigate star formations in atomic clouds, only to realize that the atomic clouds have too high a temperature and too low a density to form stars. In the 1960s, molecular clouds were identified to be the birth place of stars. Since hydrogen is the major component of molecular clouds, the syntheses of molecular hydrogen (H₂) from H atoms is critical for star formation. The measurement of [H]/[H₂] is very important for us to understand the chemical state of a targeted source (e.g., Cameron 1962; Gould & Salpeter 1963; Dickey et al. 1983). However, the observation of H I and H₂ in molecular clouds is very difficult. H₂ is a non-polarized molecule lacking permanent electric dipole moment, and thus it cannot be detected in the microwave band. Meanwhile, the crowded velocity components contributed by background and foreground emitting sources make it also hard to identify H I emission associated with any exact cloud component in the Galaxy (Li & Goldsmith 2003).

H I narrow-line self-absorption (HINSA) was demonstrated to be a good tracer of molecular clouds (Li & Goldsmith 2003, and references therein). H I self-absorption (HISA) occurs if cold atomic hydrogen is in front of a warmer emission background (e.g., Knapp 1974; Baker & Burton 1979; Wang et al. 2020). HISA is common but diverse in gas outside the solar circle. The majority of HISA features have no obvious ¹²CO emission counterparts (Gibson et al. 2000; Gibson 2010). Kavars et al. (2005) searched for HISA features within the Southern Galactic Plane Survey (SGPS), finding H I number densities of a few particles per cubic centimeter. The origins of HISA can be varied, including the cold neutral medium with low temperatures (15–35 K, Heiles 2001) and the “missing link” clouds in status between atomic clouds and molecular clouds (Kavars et al. 2005), as well as the cold H I component within the molecular clouds (Wilson & Minn 1977). As a special case of HISA, HINSA was defined as self-absorption feature with corresponding CO emission and line width comparable or smaller than that of CO (Li & Goldsmith 2003). The cold H I traced by HINSA has been shown to be tightly associated with molecular components. HINSA is an efficient method to detect cold H I mixed with molecular hydrogen H₂. Besides CO,other molecular tracers such as OH also have central velocity tightly correlated with that of HINSA (Li & Goldsmith 2003), although no clear correlation was found for the nonthermal velocity dispersion between OH and HINSA (Tang et al. 2021).

Since it was first detected in our Galaxy, the HINSA science has been developed substantially with the observations of the Arecibo 305 m Telescope (Goldsmith & Li 2005; Krčo et al. 2008; Krčo & Goldsmith 2010; Tang et al. 2016; Zuo et al. 2018). In the optically selected nearby dark cores (Lee & Myers 1999), the volume densities of H atoms (n(H)) appear to be slightly higher than the steady-state value resultant from the competition between H I formation and destruction (Li & Goldsmith 2003). Guided by ¹³CO J = 1−0, the HINSA detection rate (the fraction of ¹³CO emitting clouds where HINSA is detected) can be as high as ~80% (Krčo et al. 2008).

The Five-hundred-meter Aperture Spherical radio Telescope (FAST; Nan et al. 2011; Li et al. 2018) located in the southwest ofChina is the most sensitive radio telescope in the L band. In the FAST era, HINSA science has a new opportunity for great advancement (Heiles 2019). Tang et al. (2020) conducted a pilot H I 21 cm spectra survey toward 17 Planck Galactic cold clumps (PGCCs) using the FAST in single-point mode. In that survey, 58% of PGCCs have detections of HINSA features, with the detection rate slightly smaller than the rates in the Taurus and Perseus regions (77%; Li & Goldsmith 2003) and molecular cores (80%; Krčo & Goldsmith 2010). This deviation may arise from limited sample and different background H I emission or reflect the different evolution statuses of these sources.

The Planck satellite (Tauber et al. 2010; Planck Collaboration I 2011) provided an all-sky survey in the submillimeter-to-millimeter range with unprecedented sensitivity. PGCCs (Planck Collaboration XXIII 2011; Planck Collaboration XVIII 2016) were released as one of the products. They are mostly cold quiescent samples but are not as dense as low-mass cores, implying that they are at very early evolution stages (Wu et al. 2012; Liu et al. 2012, 2019; Juvela et al. 2015; Zhang et al. 2020; Xu et al. 2021). The detection of HINSA features in early cold molecular clouds such as PGCCs (Tang et al. 2020) reveals the potential to improve our understanding about the transition from atomic cloud to molecular cloud.

In the present work, we observed twelve PGCCs and five comparison objects using the L-band 19-beam receiver employed on the FAST. The 19-beam receiver (Li et al. 2018) can help us quickly map the target sources. Seven of the samples are observed in the snapshot mode, which is made up of four sets of 19-beam tracking mode observations and consists of 76 pointings. We introduce the basic information of the sample in Sect. 2. The FAST H I observations are described in Sect. 3.

Another feature of this work is that the sources we observed are carefully chosen as the carbon chain molecule (CCM) abundant clouds or clumps (Sect. 2). The observed PGCCs are mainly chosen from PGCCs with HC₃N emission. PGCCs with CCM production regions are chosen because they tend to contain abundant H I, since the hydrogen atom is one of the main components of hydrocarbon molecules and N-bearing CCMs (e.g., Taniguchi et al. 2016). The initial searching for carbon chain molecule (CCM) production regions began in 1970s (Turner 1971). Since then, CCMs were detected toward molecular clouds in virious environments, including the gas envelopes warmed up by young stellar objects (Suzuki et al. 1992; Sakai et al. 2008) or affected by young stellar shocks (Wu et al. 2019a,b). However, it is not clear which gas component in a PGCC is associated with the CCM production regions, the warm component accompanying young stellar objects, or the cold component embedded in the cold and dense central region. Searching for HINSA toward PGCCs could help us understand this issue. Another important reason for choosing a sample of CCM production sources is that we may have a higher detection rate of HINSA toward these sources.

Besides observations, this work also focuses on the HINSA extracting methods. The HINSA features are extracted by basically following the method of Krčo et al. (2008), denoted as “method 1” in this work. However, we point out the caveats of method 1 and give an analytic formula to quantify the requirement of high signal-to-noise (S/N) if method 1 is adopted. An improved HINSA-extracting method is proposed in this work, denoted as “method 2”. The S/N requirement for method 2 is much lower than that for method 1. Method 2 is sensitive to HINSA spectrum with high velocity resolution (thus a broad line width relative to channel width) and low S/N or HINSA features with multiple velocity components. The algorithms and the S/N thresholds of these two HINSA-extracting methods are described in Sect. 4.

Applying the improved method to the observed H I 21 cm spectra, HINSA features can be extracted in 14 of the 17 observed sources with a detection rate of 82%. The detection rate of HINSA in PGCCs is 11/12 (90%), higher than that in the sample of Tang et al. (2020). The fitting results of HINSA, HC₃N J = 2–1 and CO spectra as well as the HINSA images are presented in Sect. 5. Combining the information of HINSA, HC₃N, and CO emission as well as the dust continuum, the present work aims to improve knowledge of the H I abundances in PGCCs, theevolutionary statuses of PGCCs, and the excitation mechanism of HINSA in CO and CCM emission regions. In Sect. 6, we discuss the morphologies, abundances, and excitations of H I, CO, and HC₃N in PGCCs. We summarize this paper in Sect. 7.

2 Sample

This work focuses on the HINSA observations toward PGCCs and the HINSA extracting methods. Besides 12 PGCCs, our sample consists of five comparison objects chosen from literature. They are also CCM abundant sources, including one starless core, L1521B, as well as four star forming sources. The star forming sources are chosen from an outflow catalog (Wu et al. 2004 and references therein).

Among the 12 PGCCs, nine have been observed in HC₃N J = 2–1 using the Tian Ma radio telescope (TMRT)¹, conducted as part of the survey searching for Ku-band CCM emission in PGCCs led by Yuefang Wu. The main line of HC₃N J = 2–1 is F = 3–2 (18 196.3104 MHz). Eight of the 12 PGCCs have valid detections of HC₃N J = 2–1.

All of the 12 PGCCs have single-point observations in J = 1–0 of CO and its isotopomers by Wu et al. (2012), using the Purple Mountain Observatory (PMO) 13.7 m telescope. Among them, 11 were mapped in these CO transitions by Wu et al. (2012). The analysis in this work focuses on PGCCs. Comparing the PGCCs and the comparison objects can lead to a clearer conclusion about the effects of the core statuses and star formation activities on the morphologies of HINSA features. The H I spectra from sources of different types can also help us to test the robustness of the HINSA extracting algorithm. However, we do not try to give any quantitative analysis about those comparison objects.

The names and coordinates of sources in our sample are listed in the first three columns of Table 1. The distribution of the sources in the Galactic plane is shown in Fig. 1. The background of Fig. 1 shows the dust temperature decomposed by Planck Collaboration XLVIII (2016) from the Planck continuum images (Planck Collaboration I 2016; Planck Collaboration XLVIII 2016)² using the Generalized Needlet Internal Linear Combination (GNILC; Remazeilles et al. 2011) component-separation method. The PGCCs tend to be located at the margins of the emission regions of the Planck 353 μm continuum (Fig. 1).

For the twelve PGCCs in our sample, their distances were calculated using the Bayesian distance calculator³ (Reid et al. 2016) based on their coordinates and velocities. The Bayesian calculator constrains the probability density function (PDF) of distance based on four types of information: kinematic distance (KD), the spiral arm model (SA), Galactic latitude (GL), and parallax source (PS). The velocity dispersion has been taken into account to give the uncertainty of kinematic distance. Their angular size (D_dust) and dust temperatures (T_ECC) can be quoted from the Planck early cold core (ECC) catalog (Planck Collaboration XXIII 2011), which is the pioneer sample of PGCCs. For the other five sources, their distances are adopted from literature (see Table 1). Their D_dust can also be obtained if the Hershcel continuum data are available. The distance, D_dust, and T_ECC are listed in the fourth to sixth columns of Table 1.

3 Observation and archive data

3.1 H I observations using the FAST

The FAST⁴ is a ground-based radio telescope built in the Guizhou province of southwest China (Nan et al. 2011), and it is the most sensitive single-dish telescope. The tracking accuracy is about 0.2′ (Jiang et al. 2019). For zenith angle θ_ZA < 26.4°, the illuminated aperture is 300 meters. A 19-beam receiver in the L band is employed. The full width at half maximum (FWHM) beam size is about 3′ around 1420 GHz for θ_ZA < 26.4°, which is about half of the angular spacing between two neighboring beams, 5.87′. Our observations(2019a-020-S) were conducted on August 5 and 11, 2019. The on-source time was ten minutes for each source. Our targets all have θ_ZA smaller than 20° during the observations, as listed in Table 1. It takes approximately ten minutes or less to change the targeted sources. The spectral mode with 1024 k channels evenly spaced between 1 and 1.5 GHz was adopted for our observations, corresponding to a channel spacing of ~0.1 km s⁻¹ for H I linewith a rest frequency of ~1420 GHz.

Seven ofour sources were observed using the tracking mode with the central beam (beam 1) always pointing on the center of the target, while the other ten sources were observed with the snapshot mode (see Table 1). The snapshot mode is made up of four sets of tracking mode observations. In snapshot mode, the targeted locations of the central beams between two successive tracking mode observations are slightly shifted (see Fig. 2). The map obtained in snapshot mode is spatially half-Nyquist sampled.

The power spectra are recorded at every 0.1 s. For calibration, a high-level noise (10 K) was injected lasting one second for every two seconds. The system temperature is about 20 K for θ_ZA < 20°. For a point source, the beam efficiency is constant, with a value of ~0.6 for θ_ZA < 26.4°. It will decrease to ~0.5 if θ_ZA increases to about 36° (Jiang et al. 2019). For an extended source, the beam efficiency should be adopted as 0.85, instead. Figure 3 shows the comparison between the H I 21 cm spectrum of G173.3-16A1 obtained by the FAST and that from the Effelsberg-Bonn HI Survey (EBHIS; Winkel et al. 2016). They are consistent with each other if a beam efficiency of 0.85 is adopted for the FAST observation.

Binned with a channel spacing 0.1 km s⁻¹, the H I 21 cm spectra we obtained have rms noises of 30–40 mK for sources observed in the tracking mode and 60–80 mK for sources observed in the snapshot mode.

Fig. 1 Galactic longitude-latitude positions of the 17 sources in our sample. Yellow and cyan hexagons represent the 12 PGCCs and the other five sources, respectively. The size of the hexagon represents the covering area of the 19-beam receiver of the FAST. The dot within each hexagon shows targeted locations of the central beam (Beam 1; Fig. 2). The background image shows Planck dust temperature distribution. The black contours with levels of 3, 7, and 11 Myr sr−1 represent the Planck 850 μm continuum. The red contours represent the Hα emission (Finkbeiner 2003) in units of R (106∕4π photons cm−2 s−1 sr−1), with levels from 15 to 215 stepped by 40.

Fig. 2 Configuration of the 19-beam receiver of FAST (Li et al. 2018). Each snapshot mode observation consists of four single-point-mode observations. Equatorial coordinate system is adopted for the tracking mode, while Galactic coordinate system is adopted for the snapshot mode.

3.2 HC₃N J = 2–1 observations of nine PGCCs

Among the 12 PGCCs, nine were observed in HC₃N J = 2–1 using the Tian Ma radio telescope (TMRT). Observations were conducted in July of 2017 and August of 2018.

The TMRT is a 65-m diameter fully steerable radio telescope located in Shanghai (Li et al. 2016). The pointing accuracy is better than 10′′. The main-beam efficiency is 0.60 in the Ku band from 12 to 18 GHz (Wang et al. 2015). The front end is a cryogenically cooled receiver covering a frequency range of 11.5–18.5 GHz. The digital back-end system (DIBAS Bussa & VEGAS Development Team 2012) is employed, which supports a variety of observing modes for molecular line observations. Mode 22 was adopted for our observation.In this mode, each of the three banks (Bank A, Bank B, and Bank C) has eight subbands with a bandwidth of 23.4 MHz and 16384 channels. The center frequency of each subband is tunable to an accuracy of 10 kHz. The calibration uncertainty is 3% (Wang et al. 2015). The Ku band observations cover the transition of HC₃N J = 2–1 (18.196226 GHz). The velocity resolution is 0.023 km s⁻¹ near 18 GHz.

Fig. 3 Blue and orange solid line are the H I 21 cm spectra obtained by the FAST beam 1 and beam 8, respectively. The spectra werecalibrated adopting a beam efficiency of 0.85 for extended source. The red dotted lines are the spectra from the EBHIS (Winkel et al. 2016). The black solid lines are the FAST H I 21 cm spectra smoothed to have a velocity resolution similar to that of the EBHIS spectrum. To avoid overlapping, the red dotted lines and the black solid lines have been moved up by 10 and 20 K for beam 1 and beam 8, respectively.

3.3 Archive data of CO and dense gas tracers

The spectral of the J = 1 − 0 transitions of ¹²CO as well as its isotopomers ¹³CO and C¹⁸O were extracted from Wu et al. (2012). Following up the releasing of PGCCs, Wu et al. (2012) observed the three CO lines towards 673 PGCCs using the PMO 13.7 m telescope⁵. The morphologies and dynamic properties of PGCCs are well understood thanks to the CO mapping observations (Wu et al. 2012; Liu et al. 2012, 2013; Meng et al. 2013; Zhang et al. 2016, 2020). All 12 PGCCs have single-point observations of these three CO lines. Among them, all except for G174.08-13.2 have been mapped in these lines. The CO spectra can help to identify HINSA by providing us the initial guesses of the velocity and line width of the HINSA feature (see Sects. 4 and 5.2).

Observations of emission lines of dense gas tracers toward PGCCs are rare. Among the sources in this work, two PGCCs, G165.6-09A1 and G174.4-15A1, have been mapped in C₂H N = 1–0 (87.316925 GHz) and N₂H⁺ J = 1–0 (93.173397 GHz) using the PMO 13.7 m telescope (Liu et al. 2019). The G192.2-11A2 has been mapped in HCO⁺ J = 1–0 (89.188523 GHz) using the Instituto de Radioastronomía Milimétrica (IRAM) 30 m telescope⁶ (project ID: Delta03-13). The observation of HCO⁺ J = 1−0 in G192.2-11A2 was conducted on May 3 of 2014, with a velocity resolution of 0.16 km s⁻¹.

4 Methods of extracting HINSA

The H I 21 cm spectrum toward a molecular cloud is a mixture of foreground and background emission and absorption features originated from cold H I component in that molecular cloud. It is not easy to extract the absorption components from an H I 21 cm spectrum, since the shape of the spectrum without absorption features (T_H i) is unknown. A reference point with the same foreground and background H I emission with the target source but showing no absorption features is difficult to find.

Combining Eqs. (1) and (5) from Li & Goldsmith (2003), the observed continuum-subtracted H I spectrum brightness temperature along the HINSA sight line is (1)

Here, τ is the optical depth of the cold HI component τ_b = p(τ_f + τ_b), τ_f and τ_b denote the foreground and background H I optical depths, and T_c is the continuum brightness temperature. T_ex, T_b, and T_f are the excitation temperatures of the cold H I component, the background warm H I component and the foreground warm H I component, respectively. T_H I can then be expressed asT_r|_τ=0.

If T_f = T_b and both τ_f ≪ 1 and τ_b ≪ 1, the absorbed component of H I spectrum can be expressed as in Eq. (8) of Li & Goldsmith (2003): (2)

Equation (2) leads to the expression of T_H I as (3)

To avoid overfitting, parameters p, T_ex, τ_f, and T_c should be known parameters or fixed as reasonable values.

The excitation temperature (T_ex) of HINSA is assumed to be the value of molecular gas ~10 K, which is much smaller than T_r. T_c should be smaller thanT_ex and T_r, if it is adopted as the typical value of the L-band intensity of the Galactic background emission ~ 3.5 K (Li & Goldsmith 2003). We extracted the L-band background continuum intensities for our sources from the Bonn-Stockert survey (Reich 1982; Reich & Reich 1986), which gives a mean value of 3.9 K with a standard deviation of 0.2 K. Neglecting T_c is reasonable and may lead to an overestimation of the optical depth of HINSA by ~10%.

Each of our sources has a short distance d compared with the galactics scale. Because of the short distance and low Galactic latitude b (Table 1), the height from the Galactic plane dsin(b) is small compared with the scale height of the Galactic H I layer (~350 pc; Lockman 1984; Nakanishi & Sofue 2003). Most of the H I emitters along the line of sight should be located behind the absorption source. We believe that p should be close to 1 for most of our targets, and we adopt p = 1 hereafter. We note that, if p < 1, this assumption may lead to underestimations of the HINSA optical depth and column density. In the optical thin limit, the obtained optical depth of HINSA would be approximately pτ_real, where τ_real is the real optical depth.

If p = 1, τ_f = 0, and T_ex and T_c are ignored, Eq. (3) can be simplified as

Despite its simple form, Eq. (5) does not deviate much from Eq. (3) for HINSA, which has low optical depth.

The key step in extracting HINSA is to make reasonable estimations of the optical depths of cold H I (τ) in Eqs. (3) or (5), and to obtain the corresponding T_H I with the highest probability to represent the real unabsorbed spectrum. The HINSA extracting methods (method 1 and method 2) described below are optimized for HINSA with T_ex ≪ T_r and T_c ≪ T_r. Equation (5) helps to estimate the S/N requirement of method 1 and shed light on the motivation of method 2. Great care is required if these methods are implemented to extract absorption features of HISA clouds with warmer gas, significant foreground HI emission, or significant background continuum emission.

4.1 Method 1

HINSA usually shows significant features in the second derivative representation of T_r (Krčo et al. 2008). It is assumed the τ can be expressed as (multi-)Gaussian function: (6)

where τ_i, σ_i, and v_i represent the optical depth, velocity dispersion, and central velocity of the ith component, respectively. Free parameters (τ_i, σ_i, v_i) are requiredto be fit during the extraction of HINSA features. The initial trial values of (τ_i, σ_i, v_i) and the fixed parameter T_ex can be estimated from emission lines of other species, especially those supposed to have space distributions similar to the cold H I in molecular cloud.

The Gaussian fitting tries to obtain a smooth curve, which can be expressed as a straight line or other analytic function. However, the T_H I is unknown and cannot be assumed as a linear function of velocity or frequency. Thus, the Gaussian fitting cannot be directly applied to extract HINSA. Instead, the fitting process can be conducted under a relaxed judgment of the fit curve (T_H I); that is, T_H I should look smoothest.

T_H I can be obtained through minimizing the sum of the square of second derivative of T_r, denoted as in this work. This method (method 1) was first highlighted by Krčo et al. (2008) to extract HINSA. Here, is expressed as (7)

where v_u –v_l is the frequency or velocity range that we are interested in. The in Eq. (7) represents the numerical version of the second derivative of T_H I: (8)

where Δ_ch is the channel spacing of the spectrum.

The basic idea of this method is to find a set of parameters (τ_i, σ_i, v_i) that will reasonably fill the dips of T_r through Eq. (3). The T_H I with minimal value is the smoothest T_H I that can represent the unabsorbed spectrum.

However, this method (also denoted method 1) has theoretical defects and is noise sensitive. We give an analytic formula to quantify the S/N threshold of method 1 in Sect. 4.1. An improved method (method 2) is proposed and described in Sect. 4.2.

4.1.1 Shortcomings of method 1

T_H I usually shows a complex profile. To estimate the S/N of the HINSA detection in the second derivative representation (), T_H I is approximated by a Gaussian function with a line width of Δ_w and a peak intensity of T_peak. Here, T_peak and Δ_w = N_chΔ_ch are adopted as the typical height and typical width of T_H I, respectively. S/N can then be estimated from Eq. (8), and the second derivative of a Gaussian function (e.g., Eq. (2) of Krčo et al. 2008), and expressed as (9)

The term requires a very high S/N of the spectrum. Another problem is that, when T_r is multiplied by a factor of g during the fitting process, the noise will also be amplified by the same factor. The value of this amplification factor g can be estimated from Eq. (5) as (10)

In theory, the original T_H I will not be fully recovered through this method, since the value contributed by the noise will also be amplified. The deviation between the original T_H I and the fit T_H I is systematic, and it can be reduced but never completely eliminated by an improving S/N.

4.1.2 S/N threshold of method 1

Here, we further explain why method 1 can only partly extract the HINSA feature especially when S/N is low, and we try to quantitatively assess the performances of method 1 when applied to spectra with different S/Ns. We considered that there is a linear spectrum that has an absorption feature with a Gaussian-like optical depth: (11)

The value corresponding to a test value of τ₀ (denoted as τ_0,t) can be expressed as (12)

The fit value of τ₀ (denoted as τ_0,f) can be obtained by solving the following equation: (13)

Here, we assume that τ_0,f is small and thus e^2τ_0,f ~ 1. A recovering factor (rf) is defined as (14)

which represents the proportion of the HINSA extracted through method 1. Equation (13) leads to (15)

where σ(T_noise) is the rms of the noise. The S/N threshold for method 1 (S/N^c) is defined as the value of T_ab∕σ(T_noise) when rf = 3∕4, which leads to (16)

The recovering factor can then be expressed as (17)

If S/N(T_ab) is smaller than S/N^c, rf should be small, and Eq. (17) is only statistically valid. Thus, it is not permitted to obtain the optical depth of HINSA through dividing the fit value by rf. In such a case, method 1 either can only partly extract the HINSA feature or totally fails.

Taking the channel spacing of our spectra Δ_ch as 0.1 km s⁻¹, the absorption line width as 1 km s⁻¹ (corresponding to N_ch = 10), and τ = 0.1, Eq. (16) leads to S∕N^c > 30. Keep in mind that S∕N^c is defined on T_ab scale, and it will be multiplied by a factor of ten if defined on a T_r scale. Such a high S/N cannot always be achieved in observations.

4.2 Method 2

The main difficulty of method 1 is that it recovers T_H I by multiplying T_r by a factor of ~e^τ, which leads to the multiplication of noise simultaneously. To inhibit the increasing of noise, we replace the multiplication on T_r in Eq. (5) (or Eq. (4)) by an addition operation. After a simple manipulation, Eq. (5) can be rewritten as (18)

where means the fit T_H I. This equation highlights the possibility of keeping the noise component in T_H I constant during the fitting process and equal to that in T_r if the T_H I term on the right side is replaced by its smoothed version T_H I,smooth. It is reasonable since we expect that T_H I is smooth and can be approached by low-order polynomial.

More specifically, the fitting process of method 2 consists of five steps, including

• 1.

  Calculating T_H I through Eq. (3) for each given test value of free parameters.

• 2.

  Obtaining T_HI,smooth, a smooth version of T_H I; for example, by n_th-order polynomial fitting.

• 3.

  Updating T_H I using Eq. (18).

• 4.

  Calculating the value.

• 5.

  Repeating the above steps to find the minimal value and its corresponding T_H I.

4.2.1 Advantage of method 2

Since most of our HINSA features are weak, it is safe to assume that τ₀ ≪ 1, which would produce a near-Gaussian absorption dip. Method 2 is to recover T_H I by adding a smoothed near-Gaussian component to T_r. It splits the multiplying term e^τT_HI in Eq. (5) into two terms: one keeps the noise and details, the other is smoothed. What we try to do is to balance the resolution and S/N. In method 2, the noise will not increase much, and the velocity resolution is kept since the smoothing is not conducted on T_r, but on T_H I which is expected to be smooth. In this work, we used polynomial fitting to obtain a smooth version of , but other methods such as convolution could also be used.

Once is obtained, the absorbed component can be calculated through (19)

where represents the fit T_H I, and T_ab should be apositive quantity. If the HINSA is optically thin, the noise intensity of T_ab given by method 1 would be small with a value of approximately τT_noise. If method 2 is adopted, the calculated T_ab should look even smoother since the amplification of noise in has been restrained.

4.2.2 S/N threshold of method 2

We now quantitatively assess the performance of method 2 and compare it with method 1. Method 2 works when (20)

which leads to (21)

Here, δ(X) represents the standard deviation of X. It is importantto note that T_ab is a fit curvethat does not contain noise. If Eq. (21) is satisfied, we say that the absorption feature against the noise has a significance level higher than five sigma when extracted using method 2. Equation (21) can be satisfied when S/N(T_ab) is higher than the S/N threshold for method 2, denoted as S/N^c;n. To derive S/N^c;n from Eq. (21), we used (22)

and the approximation (23)

Equation (21) can then be approximated as (24)

which leads to (25)

From Eq. (25), we can say that method 2 is analogous to the Gaussian fitting, since both of them have sensitivities proportional to the square root of the line widths. In contrast to method 1, method 2 favors extracting features with relatively large line widths resulting from multiple absorption components or high velocity resolution. The S/Ns of the spectra required by method 2 are much smaller than those by method 1.

4.3 Comparison between method 1 and method 2

Method 2 is an updated version of method 1. It overcomes the downsides of method 1 by introducing a intermediate step of smoothing, and thus has less strict requirement on S/N. If the T_H I spectral structure is too complex, method 2 may lose resolution.

If the absorption dip deviates from Gaussian shape, the fit unabsorbed spectrum () may be larger or smaller thanthe real T_H I. For method 1, the fit tends to be smaller than T_H I, because a solution of larger than T_H I will be more severely penalized by the amplification of noise. For method 2, the fit has no system bias of being larger or smaller than the real value. Thus, method 2 is less predictable, and method 1 seems to be more robust when the absorption dip deviates much from Gaussian shape. However, the results of the two methods in such cases are both somewhat unreliable.

In practice, method 1 is applicable only when S/N > S/N^c (Sect. 5.2). It confirms the validity of the threshold of the S/N (S/N^c) defined in Eq. (16), which describes the requirement of method 1. It also highlights the usefulness of method 2. For molecular clouds such as PGCCs, S/N ~S/N^c requires a FAST on-source integration time of ~10 min, where the S/N is that of the absorption feature in the observed H I spectrum. For FAST drift-mode observations (Li et al. 2018), the effective integration time is ~40 s (Zhang et al. 2019). Method 2 provides a possibility to extract HINSA features from spectra with such short integration time.

5 Result

5.1 CO parameters

Among the 12 PGCCs in our sample, all except for G174.08-13.2 have been mapped by Wu et al. (2012) in J = 1–0 of ¹²CO and its isotopomers (Sect. 3.3). To match the H I 21 cm spectra, the J = 1–0 spectra of ¹²CO and its isotopomers of the 11 mapped PGCCs were smoothed to have an angular resolution of 3 arcmin. Figure 4 shows the CO spectra extracted from the locations of the central beam of the H I observations. Those spectra are fit to obtain line parameters, which can help us to extract HINSA features from the H I 21 cm spectra (Sect. 5.2). For G174.08-13.2, the CO parameters derived from its single-point observations by Wu et al. (2012) are adopted.

Assuming ¹²CO J =1−0 is optically thick (τ ≫ 1) and the beam filling factor is a unit, the excitation temperature can be expressed as (Kwok 2007) (26)

where is the peak brightness temperature of ¹²CO J = 1–0, the background temperature (T_bg) can be adopted as 2.73 K, and (27)

Among the 11 sources observed by Wu et al. (2012), nine have single-peak spectra of ¹³CO J = 1–0. We fit them with a single Gaussian profile. The other three sources, G170.88-10.92, G172.8-14A1, and G174.4-15A1, have double-peak spectra of ¹³CO J = 1–0, and fittings of double Gaussian profiles are applied. Gaussian fitting is also applied to the spectra of ¹²CO J = 1–0 and C¹⁸O J = 1–0, with the number of velocity components same as that of ¹³CO J = 1–0.

The column densities of ¹³CO and C¹⁸O are calculated following the method described in Wu et al. (2012). Adopting the typical abundance ratios X[H₂]/X[¹³CO] = 89 × 10⁴ (McCutcheon et al. 1980; Pineda et al. 2013) and X[H₂]/X[C¹⁸O] = 7 × 10⁶ (Frerking et al. 1982), the H₂ column densities (H₂) and (H₂) can be calculated. The peak brightness temperatures of ¹²CO J = 1–0, the Gaussian parameters of ¹³CO J = 1–0 and C¹⁸O J = 1–0, and the gas column densities derived form the CO emission are listed in Table 2.

Fig. 4 Spectra of J = 1–0 of 12CO, C13O, and C18O (Wu et al. 2012) are shown from top to bottom in each panel.

5.2 Extracting HINSA

Both method 1 and method 2 are applied to extract HINSA features in our spectra. The free parameters to be fit are τ_i, v_i, and σ_i with i = 0, 1,..., m − 1 (Sect. 4), and m is the number of velocity components determined from CO emission or quoted from literature. For all sources, τ_i is initially set as 0.1. T_ex is fixed as 10 K, and it will only introduce a small value of uncertainty (Sect. 5.2.4).

For PGCCs, initially, v_i is adopted as the central velocity of ¹³CO J = 1–0 (Table 2), and σ_i is adopted as (28)

Here, the ΔV represents the FWHM of the spectral line, and the correction for the different molecular weight between ¹³CO and H I has been conducted.

For comparison objects, the initial guesses of central velocities are adopted as the system velocities quoted from literature. The system velocities are 6.6 km s⁻¹ for L1489 (Wu et al. 2019a), 6.4 km s⁻¹ for L1521B (Hirota et al. 2004), 9.0 km s⁻¹ for IRAS 05413-0104 (Zhang et al. 2021), 10.0 km s⁻¹ for HH25 MMS (Gibb & Davis 1997; Wu et al. 2004), and 0.5 km s⁻¹ for CB34 (Codella & Scappini 2003). The initial values of line widths are set as 1 km s⁻¹. We used the L-BFGS-B algorithm implemented by the Python package Scipy⁷ to minimize and obtain corresponding fit parameters.

5.2.1 Performances of method 2 toward individual objects

The fitting results of G174.4-15A1 are shown in Fig. 5 as examples to present the different performances of the two methods. For G174.4-15A1, there are two gas components with velocities traced by CO emission close to each other (Table 2). Method 1 can partly extract HINSA corresponding to the bluer velocity component (V_LSR ~ 5 km s⁻¹), but it fails to extract HINSA corresponding to the redder velocity component (V_LSR ~ 6.8 km s⁻¹). This result is expected by Eq. (16). For the bluer component, S/N(T_ab) is similar to S/N^c. However, for the redder component, S/N(T_ab) is smaller than S/N^c (see Table 3). If method 2 is adopted, both of the two HINSA features can be extracted. This example confirms that method 2 can extract absorption features from H I spectra with low S/N and blended absorption components.

To test the robustness of method 2, we also compared the fitting result of the two methods toward the spectra of G170.88-10.92. The H I spectrum of G170.88-10.92 shows a single dip with a high S/N. The fit results of the two methods are very close to each other (see Fig. 6). For the center beam, the optical depths fit by the two methods have a deviation of less than 10 percent. The performances of the two methods are similar to each other when the S/N is high. The fitting result of method 2 should be reliable. This is confirmed by the similarity between the distributions of CO emission lines and H I absorption features for both G174.4-15A1 and G170.88-10.92 (see Sect. 5.2.3).

Fig. 5 HINSA features of G174.4-15A1 extracted by method 1 (left) and method 2 (right). Method 1 fails to fit the HINSA component because of the low S/N and the blending of two absorption components with similar velocities. In each panel, the gray line is the observed spectrum Tr, the blue line is the recovered TH I, and the red line is Tab. The Tab shown here has been multiplied by a factor of five.

5.2.2 Performances of method 2 toward central beams

Method 2 was adopted to obtain the HINSA parameters for all the H I 21 cm spectra. For each source, the observed spectrum of the central beam (beam 1) as well as the extracted HINSA are shown in Fig. 7. Among the 12 PGCCs, ten have HINSA features on the central beam spectra that can be successfully extracted using method 2. The exceptions are G160.51-17.07 and G165.6-09A1. It corresponds to a HINSA detection rate of 83%.

The fit parameters of HINSA, including the optical depth, central velocity, line width of HINSA and the height of the absorption dip (T_ab) are listed in Table 3. S/N(T_ab), defined as the ratio between the peak intensity of T_ab fit through method 2 and the noise intensity T_noise, is also listedin Table 3.

For a central beam spectrum that can be successfully fit using method 2, the S/N threshold of method 1 (S/N^c) is calculated through Eq. (16). To obtain S/N^c, the values ofτ₀ and ΔV fit through method 2 are used, and N_ch is adopted as ΔV∕Δ_ch. Most of our spectra have S/N(T_ab) comparable to or even larger than S/N^c. The sources whose HINSA features can be partly extracted using method 1 are labeled in Table 3. From Fig. 8, it is clear that spectra with successful fittings by method 1 all have S/N(T_ab) > S/N^c. This impliesthat the S/N threshold of method 1 given by Eq. (16) is reasonable, and it is essential to apply method 2 to extract HINSA features from the H I spectra of this present work.

If method 1 is applied, HINSA features can only be extracted in seven PGCCs (Table 3), corresponding to an extracting rate of 58%. This value is similar to the value of Tang et al. (2020) but lower than the HINSA detection rate of method 2.

Fig. 6 Same as Fig. 5, but for G170.88-10.92.

Fig. 7 H I 21 cm spectra at the center beam. In each panel, the gray line is the observed spectrum Tr, the blue line is the recovered TH I using method 2, and the red line is Tab. Values of Tab on the y-axis range from 0 to 20 K for all sources except G174.08-13.2, which uses 0–60 K. The orange and blue squares represent the sources that can and cannot be successfully fit using method 1. The blue line is not shown for spectrum with fit τ smaller than 0.01.

5.2.3 HINSA maps

The maps of the integrated intensity of the optical depth of HINSA features extracted by method 2 are overlaid by the contours of CO emission. The sources observed in snapshot mode and tracking mode are shown in Figs. 9 and 10, respectively. Comparing the emission regions of HINSA and molecular lines can help us to confirm the detections of HINSA. It can also help us to study the correlation between the distribution of cold H I and other species.

We find that HINSA features can be extracted in the north-west margin of the map of G165.6-09A1. This makes G160.51-17.07 the only observed PGCC that does not show any HINSA feature. The detection rate of HINSA features in PGCCs could be as high as 11/12 (~90%) using method 2 if H I spectral of all beams are taken into consideration.

Among five comparison objects, only L1251B, HH25 MMS, and CB34 have HINSA features, with a detection rate lower than the value in PGCCs. The morphologies and environments of HINSA in PGCCs and comparison objects are further discussed in Sect. 6.2.

5.2.4 Uncertainties of HINSA parameters

The T_ex of CO is close to 10 K (Sect. 6.1). The T_ECC of all sources except G165.6-09A1 have values of about 10 ± 2 K. Thus, it is reasonable to adopt the excitation temperature of HINSA as 10 K (Sect. 5.2). For T_r ~ 50 K and a uncertainty of excitation temperature σ(T_ex) = 2 K, the relative error of the fit optical depth can be estimated as (29)

which is comparable with the values transferred from the data errors (Table 3).

It is important to note that the H I column density traced by HINSA is the value for cold H I mixed within a molecular cloud. If there is no background emission with T_H I > T_ex, there will be no HINSA, and this may lead to the underestimations of the column density of cold H I and the detection rate of HINSA.

Although method 2 can significantly reduce the requirement for a S/N, the HINSA features may still be difficult to extract if T_H I is too complex. If the value of the fit optical depth is too low, it would not be easy to judge if the extracted absorption feature is a real HINSA signal. The lower limit of τ to be safely extracted is discussed in Sect. 5.2.5.

Fig. 8 Distribution of S/Nc and S/N(Tab) for the H I spectra of the central beams (Table 3). Here, S/N(Tab) is the ratio between the intensity of absorption dip revealed by method 2 and the spectra noise, and S/Nc is the signal-to-noise threshold required by method 1 (Sect. 5.2.2). The orange and blue squares represent the sources that can and cannot be successfully fit using method 1, respectively. The black line represents x = y. As predicted, the black line is the separation line between the two cases.

5.2.5 Lower limit of τ to be safely extracted

Assuming both T_H I and T_ab have Gaussian shapes with their central velocities aligned, and their line widths are 10 and 1 km s⁻¹, respectively, T_r will only have a dip-like structure when (30)

which leads to (31)

Only HINSA with τ larger than 0.01 can be safely extracted. If the fit optical depth τ_0,f has value smaller than 0.01, the fitting results should be considered as unreliable. In such case, the fit absorption feature may be a misleading signal introduced by the irregularity of the background H I emission. If the fit optical depth is smaller than or equal to 0.01, it would be listed in Table 3 as 0.01, which means that the fitting result is not reliable.

5.2.6 Column density of H I traced by HINSA

The column density of H I can be calculated through (Li & Goldsmith 2003) (32)

The total gas column densities (N(H₂)) are adopted as the values derived from CO spectra (Sect. 5.1) to calculate the abundances of cold H I.

The values of column densities and abundances of H I are also listed in Table 3. We provide discussions concerning the H I abundances of PGCCs in Sect. 6.3.

5.3 HC₃N

Nine PGCCs have been observed in HC₃N J = 2–1 using the TMRT 65-m telescope (Sect. 3.2). To compare the intensities and line widths of HINSA and HC₃N in PGCCs, the spectra of HC₃N J = 2–1 are shown in Fig. 11. The 1sigma noise levels are about 45 mK when the spectra are binned to have a channel width of ~0.05 km s⁻¹. All of those PGCCs, except for G178.98-06.7, have a valid detection of HC₃N J = 2–1. Although the CO spectra of several PGCCs such as G159.2-20A1, G172.8-14A1, and G174.4-15A1 have double velocity components with separations of ~1 km s⁻¹, all of the detected HC₃N J = 2–1 spectra show single velocity components.

The HC₃N J = 2–1 spectra are fit using the hyperfine structure (HFS) fitting procedure provided by Gildas/CLASS⁸. If the fit optical depth of the spectrum tends to be smaller than 0.1, the fitting procedure will give a fitting result with optical depth adopted as 0.1. The fit parameters include the peak intensity (T_peak), central velocity (V_LSR), line width (ΔV), and the optical depth (τ) of the main line of HC₃N J = 2–1. These parameters are listed in Table 4.

G174.06-15A1 and G174.08-13.2 have peak brightness temperature of HC₃N J = 2–1 larger than 1 K. Their optical depths are 0.4 ± 0.1 and 0.5 ± 0.1, respectively. Assuming that the beam filling factor f is unit and the background temperature T_bg = 2.73 K, the excitation temperatures of HC₃N J = 2–1 would be derived as 6.4 K and 7.7 K for G174.06-15A1 and G174.08-13.2, respectively. For G165.6-09A1 and G172.8-14A1, the fit optical depths are larger than 1. The reason may be that the S/Ns are not high enough for these two sources. The HC₃N J = 2–1 of G159.2-20A1, G173.3-16A1 and G174.08-13.2 are optically thin with optical depths smaller than 0.1.

Under local thermodynamic equilibrium (LTE) assumption, the equation to calculate the molecular column density is (Garden et al. 1991; Mangum & Shirley 2015) (33)

where B, μ, and Q are the rotational constant, the permanent dipole moment, and the partition function, respectively. These transition parameters can be obtained from “Splatalogue⁹”. The beaming factors may be much smaller than one for most sources, except for G174.06-15A1 and G174.08-13.2. The excitation temperatures are adopted as 5 K to calculate the beam-averaged column densities of HC₃N. The column densities of HC₃N are also listed in Table 4. For G165.6-09A1 and G172.8-14A1, if the excitation temperature given by HFS fittings is adopted, the value of the calculated column densities will only change by less than 5 percent. The average abundance of HC₃N is 2.2 × 10⁻⁹ with a standard error of 0.8 × 10⁻⁹.

Fig. 9 Backgrounds: strengths of HINSA measured in the snapshot mode, extracted using method 2. Black contours: 13CO J = 1–0 emissions, which have been smoothed to the angular resolution of H I observations. Big black boxes: areas covered by CO data. Blue crosses: sampled points of the FAST, and red circle: center beam. Small red boxes and triangles: Class I/II and Class III YSOs. In panels a and c, the white contours represent the emission of C2H N = 1–0 within the white rectangles (Liu et al. 2019). In panel f, the small black rectangle represents the mapped region of the IRAM HCO+ J = 1–0 and C18O J = 2–1 observations, which are shown in Fig. 14.

Fig. 10 Backgrounds: strengths of HINSA measured in the tracking mode. Symbols have the same meanings with those in Fig. 9.

6 Discussions

6.1 CO and dust emission

The left panel of Fig. 12 shows the correlation between the H₂ gas column densities derived from ¹³CO (N^13CO(H₂)) and C¹⁸O ((H₂)). N^13CO(H₂) and N^C18O(H₂) are consistent with each other within a factor of 3. For sources with low column densities ((H₂) < 10²²), the gas column densities derived from C¹⁸O tend to be smaller than those derived from ¹³CO. This may result from the low S/Ns for C¹⁸O J = 1–0 when the gas column densities are small. Another explanation is that C¹⁸O is more easily photo-dissociated by the external radiation field compared with ¹³CO, especially for sources with low column density and extinction, such as PGCCs (Shimajiri et al. 2014; Schultheis et al. 2014; Shetty et al. 2011).

Planck Collaboration I (2016); Planck Collaboration XLVIII (2016) decomposed the all-sky maps of foreground dust column density and temperature. For sources in this present work, we extracted the values of dust temperatures (T_dust) from the temperature map and the gas column densities (N^dust(H₂)) from the dust column density map, according to their coordinates. Here, a gas-to-dust mass ratio of 100 is assumed to calculate the gas column density (N^dust(H₂)) from the dust column density. T_dust and N^dust(H₂) are listed in the seventh and eighth columns of Table 1, respectively. The right panel of Fig. 12 shows the correlation between N^dust(H₂) and (H₂). The gas column densities derived from dust emissions and CO spectra are consistent with each other, and their ratios vary within a factor of three.

From ¹³CO intensity maps, which have been smoothed to have angular resolutions of 3 arcmin, we measured the angular diameters of the ¹³CO emission regions (D_CO). D_CO is comparable or even larger than D_dust, the diameter of the Planck dust emission region (see Sect. 3). This means that the CO emissions, if smoothed to have similar angular resolutions to those of the Planck continuum, will look as extended as the Planck dust emission regions. CO emission is coupled well with the dust. The CO cores and subclumps extracted from CO intensity maps may be parts of the continuous CO emission structures with column densities locally enhanced. This is in contrast to the case of C₂H. In G165.6-09A1, G174.06-15A1, and G173.3-16A1, the emission regions of C₂H (also have been smoothed to have angular resolutions of 3 arcmin) are much more compact than CO emission regions (Figs. 9 and 10).

Figure 13 shows the correlation between CO J = 1–0 excitation temperatures (T_ex(CO)) and dust temperatures T_ECC and T_dust (see Sect. 3 and Table 1). T_ex(CO) is estimated from the peak brightness temperatures of ¹²CO J = 1–0 (T₁₂) through Eq. (26), assuming the ¹²CO J = 1–0 lines are optically thick and the beam filling factors are equal to a unit. T_ex (CO) tends to be larger than the dust temperatures of the PGCCs derived from the SED fittings of the continuum fluxes (T_ECC), but smaller than the values of T_dust extracted from the decomposed dust maps (Planck Collaboration I 2016; Planck Collaboration XLVIII 2016). The mean values of T_ex(CO), T_ECC, and T_dust (in unit of K) are 12.6(0.8), 10.8(0.7), and 17.3(0.7), respectively. The numbers in parentheses are the values of corresponding standard errors. The temperatures derived here are generally similar to those measured for infrared drak clouds (Xie et al. 2021). T_ECC should be the temperatures of the coldest dust components, since large-scale emissions were already filtered out before the application of SED fitting. On the contrary, the contribution of the extended warmer dust components makes the value of T_dust higher. Part of CO emission may come from more diffuse region besides the regions associated with the coldest dust components. It can also explain the extended morphology of CO emission. This reflects the nature of PGCC, the coldest and quiescent molecular cloud with a mixture of relatively diffuse and dense components. This is also consistent with the high detection rate of HINSA especially in HC₃N-harboring PGCCs (Sect. 6.5).

Fig. 11 HC3N spectra binned with a velocity resolution of 0.05 km s−1. The red lines show the results of HFS fitting.

Fig. 12 Relations between the column densities of H2 derived from CO and dust. Upper: correlation between the column densities of H2 derived from13CO (N13CO(H2)) and C18O (NC18O(H2)) at the positions of the FAST beams for 11 PGCCs (G174.08-13.2 is not included because of lacking CO mapping observations). Data of different sources are shown in different colors. Lower: correlation between the H2 column densities at the map centers derived from dust emissions (Ndust(H2)) and C18O (NC18O(H2)).

Fig. 13 Comparison between CO J = 1–0 excitation temperatures and Planck dust temperatures. The black line represents y = x.

6.2 Morphologies and environments of HINSA

In the scales of source sizes, the morphologies of HINSA emission are similar to those of ¹³CO for most PGCCs detected with HINSA features (Figs. 9 and 10). For G165.6-09A1 and G174.4-15A1, the ¹³CO J = 1–0 emission and the HINSA features are weak at the map centers. This is not caused by pointing bias. In panel a and panel c of Fig. 9, the blue contours represent the emission of C₂H N = 1–0 (Liu et al. 2019). At the map center of G165.6-09A1, emissions of C₂H N = 1−0 were detected. Figure 14 shows the HCO⁺ J = 1−0 and C¹⁸O J = 2−1 emission toward the map center of G174.4-15A1. The HCO⁺ J = 1−0 and C¹⁸O J = 2−1 emission tends to be enhanced in the southwest corner (Fig. 14), in constract with the CO J = 1−0 emission extending to the north (panel f of Fig. 9). These results show that HINSA features are more tightly correlated with the emission of CO than that of dense gas tracers.

However, there are also several exceptions to the similarity between the emission regions of CO and HINSA. For example, in G159.2-20A1 (see the first panel of Fig. 10), this correlation is weak. For a molecular cloud, its large-scale environment and the star formation activities within it may have an influence on the existence and distribution of HINSA.

HINSA features are not detected (τ₀ ≲ 0.01) at the central beam toward G165.6-09A1 and G160.51-17.07, and marginally detected (τ₀ ~ 0.3) toward G159.2-20A1. There are no central-beam HINSAs detected in G165.6-09A1 and G159.2-20A1, while G160.51-17.07 has no detected HINSA anywhere. The three sources (G165.6-09A1, G160.51-17.07, and G159.2-20A1) all have T_dust larger than 19 K, the largest compared with other PGCCs (Table 1). All PGCCs are located at the margin of large-scale dust emission regions traced by the Planck 353 GHz continuum (Fig. 1). G160.51-17.07, the only PGCC showing no HINSA feature, is located within an H_α emission region represented by red contours in Fig. 1. These results all demonstrate the importance of the environments for the detection of HINSA in PGCCs and the local deviation between HINSA and CO emission regions.

Star formation activities also have significant influence on the detections of HINSA. Among 12 PGCCs, 11 have detections of HINSA, and the detection rate is 90%. HINSA features are also detected in the starless cores, L1521B. In four star formation regions, L1489E, IRAS 05413-0104, CB34, and HH25 MMs, HINSA features are only detected in CB34 and HH25 MMS. The HINSA features detected in HH25 MMS are not very reliable. HH25 MMS was observed in the tracking mode instead of the snapshot mode. There are only several (<5) beams have H I spectra with HINSA features. In the H I spectrum of HH25 MMS (Fig. 7), besides the dip with velocity similar to those of CO lines (~10 km s⁻¹), there is another dip with redder velocity (~13 km s⁻¹) where no corresponding CO emission lines can be seen. These two dips may be misleadingly produced by a H I emission peak between them. HINSA features tend to be inhibited by star formation activities.

Inhibitions of star formation activities on HINSA features can also be seen in PGCCs. Based on the WISE data, Marton et al. (2016) presented an all-sky catalog of 133 980 Class I/II and 608 606 Class III young stellar object (YSO) candidates. All sources in our sample have WISE YSO associations within the fields of view, except for G172.8-14A1 (see Figs. 9 and 10). The YSOs tend to locate at the margin of the HINSA detection regions, especially for sources G174.4-15A1, G175.34-10.8, G192.2-11A2, and G170.88-10.92. In G165.6-09A1, the map center is associated with four YSOs, and no HINSA feature is detected. However, at the northwest margin of G165.6-09A1, there are strong HINSA features (τ_HINSA ~ 0.1) but no YSO association. Since we do not know the velocities of YSOs, the YSOs matched based on angular separations may not be truly associated with PGCCs in 3D space. The relations between the distribution of YSOs and HINSA confirm that WISE YSOs are spatially related to the PGCCs, as pointed out by Marton et al. (2016) according to the YSO surface distributions around PGCCs.

Overall, HINSA tends not to be detected in regions associated with warm dust emission and background H_α emission. In the region associated with young stellar objects or dense gas emission, the HINSA feature also tends to be inhibited. Besides these effects, the varied abundances of H I in molecular clouds (Sect. 6.3) and the different excitation conditions between HINSA and CO (Sect. 6.4) may also lead to the different distributions between the HINSA and CO emission in some sources.

Fig. 14 Cutout region of G192.2-11A2 represented by the black box in Fig. 9f. The background and black contours present the emissions of C18O J = 2–1 and HCO+ J = 1–0 measured with the IRAM 30 m telescope (Xu et al., in prep.), respectively. The contours have levels from 0.5 to 0.9 stepped by 0.1 in unit of K km s−1. The red circle shows the location of the central beam of the H I observation. The yellow dashed circle and solid circle represent the beam size of HCO+ J = 1–0 and C18O J = 2–1, respectively. The white dashed lines are the grid lines of the Galactic coordinate.

6.3 H I abundances and evolution status

The H I traced by HINSA is assumed to be well mixed with molecular gas (this issue is further discussed in Sect. 6.4). Figure 15 shows the correlation between the column densities of H I (N(H I)) derived from HINSA and the column densities of H₂ (N(H₂)) derived from CO spectra. N(H I) and N(H₂) are positively correlated, especially for data at the map centers. The derived abundances of H I are approximately 3 × 10⁻⁴, varied by a factor of ~3. This confirms the similarity between CO and HINSA.

If a molecular cloud has a typical distance d = 400 pc, angular diameter Θ = 5′, and H I column density N(H I) = 2 × 10¹⁸ cm⁻² (Table 3), the typical volume density of H I can be estimated as (34)

This value is consistent with the value at the steady state (n^st (H I)) given by Li & Goldsmith (2003). If the beam dilution effect and the relatively large uncertainty of N(H I) (Fig. 15) are considered, the volume density of H I may exceed the n^st (H I). In some work such as Tang et al. (2020) some sources with large H I abundances (~10⁻²) were reported, and those sources were thought to be in transition phase between atomic and molecular states.

Fig. 15 Correlation between the column densities of H I and those of H2 derived fromC18O J = 1–0 emission. The orange squares represent the data of the central beams. The blue dots represent the data of other beams. The red square shows the data at the peak position of the HINSA map of G165.6-09A1.

6.3.1 Production efficiency of H₂ from H I

The steady-state value of H I abundance reflects the balance between the destruction rate of H₂ by cosmic rays and the formation rate of H₂ on the grain surfaces. However, the formation efficiency of H₂ relies on the temperature and grain types. When the dust temperature is low (typically several K depending on the binding energies of hydrogen), the grain will be covered by a monolayer of H I or H₂, which prevents H I further sticking onto the grain surface (Hollenbach & Salpeter 1971). When the dust temperature is higher (> 15 K in astrophysical environments; Katz et al. 1999), the formation efficiency of H₂ is small since the depleted H I will be thermally desorbed before reacting with another H I. Thus, H₂ will only be efficiently produced within a narrow temperature range (Hollenbach & Salpeter 1971; Katz et al. 1999). The H₂ formation at higher temperature can be efficient if chemisorption site is taken into consideration (Cazaux & Tielens 2002), and this has a much higher binding energy (~10⁴ K) than the value for the physisorption site (~500 K) and provides a shield for H I when the dust temperature is high (>20 K).

The assumption that a physisorption site occupied by one H I will not further absorb another H I is not necessary. When the dust temperature is low, there may be several H I in a single physisorption site. The bound energy of H-H is 4.52 eV (Duley 1996). If the reaction heat of a forming H₂ is able to peel all H atoms in that site from grain surface, only the empty site is then truly capable of absorbing H I in gas phase. Under this assumption, the formation efficiency of H₂ is not altered, and absorption and synthesis of other elements and molecules on grain surface at low dust temperatures will not be weakened by the occupation of sites by hydrogens. The production efficiency of H₂ depends on the grain properties.

6.3.2 Evolution status

The sources with a relatively higher abundance of H I may still be in steady state if their dust temperatures and grain properties are not suitable for H₂ formation. There is no HINSA in G160.51-17.07, which has a relatively high dust temperature (19.2 K) but not a weak CO emission (Fig. 4). However, we cannot rule out that G160.51-17.07 is in transition phase, considering its relatively lower surface density of dust compared with other PGCCs (Table 1).

The derived value of n(H I) for PGCCs may also be underestimated. It is not sure if the column densities of H I are evenly distributed within regions shown HINSA features. The beam dilution effect may be not ignorable, just like the case of CO emission.

6.4 Does the HINSA we detected trace molecular gas?

We assumed that HINSA mainly traces H I that is well mixed with molecular gas in Sect. 6.3, based on the fact that the morphologies of HINSA emissions are overall correlated to those of ¹³CO in PGCCs (Sect. 6.2). However, it is still not clear if the H I absorption contributed by the H I component around molecular gas cloud is also important. The spatial resolution of the H I observations of this present work prevents us from digging deeply into this issue. However, from the discussion below about the excitation status of HINSA against background H I emission, we try to illustrate that H I traced by HINSA should be a gas component of molecular cloud and confined within and around CO emission regions.

6.4.1 Collisional excitation of HINSA

The H-H spin-exchange collision is the main mechanism of the collisional excitation of H I if the abundance of H atom overwhelm that of electron (Field 1958; Furlanetto & Furlanetto 2007). The coefficient of H-H spin-exchange collision tabulated in Col. (4) of Table 2 of of Zygelman (2005) can be empirically fit as (35)

for 10 < T_k < 300 K.

In optically thin conditions, the excitation temperature of HINSA can be estimated through (Purcell & Field 1956; Li & Goldsmith 2003) (36)

in which (37)

T_k is the kinematic temperature, T_background is the brightness temperature of background emission, n_H I is the volume density of hydrogen atom (the major collision partner to excitate the H I 21 cm line), and A_1,0 is the spontaneous decay rate 2.85 × 10⁻¹⁵ s⁻¹ (Wild 1952).

If T_k = 10 K is adopted, Eqs. (35) and (37) lead to y ~ 6.9n_H I. Here, n_H I is in units of cm⁻³. For cold H I within a molecular cloud, if we further adopt n_H I as the typical volume density of H I given by Eq. (34) or the steady-state value, the value of y would be larger than seven. If HINSA is not fully excited, the density of H I given by Eq. (34) will be underestimated, which leads to an even larger value of y. Thus, within a molecular cloud, the density of H I in steady state should be larger than the critical density of HINSA. If a molecular cloud is in transition phase or has grain properties not suitable for H₂ formation (Sect. 6.3), its density of H I would be even higher. It is safe to assume that cold H I within molecular cloud is collisionally excited.

Fig. 16 Blue line: average H I spectrum of the observed sources in this work. Orange line: estimated bymultiplying the blue line by sin(14°). Green line: evaluated from the full-sky H I survey HI4PI based on EBHIS and GASS (HI4PI Collaboration 2016).

6.4.2 External excitation of HI 21 cm line

For a diffuse gas component around a molecular cloud with high abundance but low H atom density, the H I 21 cm line may be mainly excited by the background emission. The L-band continuum brightness temperature (T_Lband) is on average smaller than4 K (Sect. 4). This is consistent with the value estimated adopting the intensity of L-band standard interstellar radiation field (ISRF) as ~0.8 K (Winnberg et al. 1980; Li & Goldsmith 2003) and that of the cosmic microwave background as 2.73 K. If T_L band is adopted as the intensity of external emission, the absorption feature contributed by diffuse H I component may not be ignorable.

We note that the L-band external emission is mainly contributed by , the brightness temperature of the averaged background H I emission: (38)

The sources we observed all have |V_LSR|≲ 10 km s⁻¹ and Galactic latitude |b|≲ 20°. The mean value of |b| is 14°. The lower limit of the intensity of can be estimated through multiplying the average spectrum of observed sources by sin (14°), as shown in Fig. 16. The peak intensity of is expected to be ~14 K. This is confirmed by the H I spectrum averaged over the all-sky spectra from the HI4PI survey (HI4PI Collaboration 2016), which gives brightness temperatures of ~20, ~10, and ~5 K for V_LSR = 0, 10, and 20 km s⁻¹, respectively (Fig. 16). These values are larger than T_L band.

For a molecular cloud with short distance to the Earth and thus similar compared with that seen from the Earth, T_background will be dominated by the average brightness temperature of H I 21 cm spectra if |V_LSR| is small (|V_LSR| < 20 km s⁻¹). If is small, the warmer and diffuser H I around the molecular cloud may be excited by T_L band. Then, the H I in the diffuse gas may also contribute absorption features similar to HINSA and make us overestimate the column density of cold H I associated with molecular gas. If is not ignorable, the diffuse H I will not contribute much absorption, even if it is not collisionally excited. Then, we can safely assume that the HINSA features we extracted are tightly associated with molecular gas. This effect inhibits the possible absorption contributed by the outer diffuse H I-rich shell of a molecular cloud, and it confines HINSA within molecular cloud.

Fig. 17 Comparisons between the velocity dispersions and column densities of HC3N and H I. Left: cumulative distributions of velocity dispersions. Solid lines and dashed lines represent σ and σNT, respectively (Sect. 6.5). Right: correlation between the column densities of H I and those of HC3N. Green and orange lines represent the linear fitting with and without the two data in blue shadowing, respectively. For data in the blue shadowing, the N(H I) parameters are not reliable with regard to the fit optical depth of HINSA τ < 0.01 (Sect. 5.2.5).

6.4.3 HINSA and CO emission region

A molecular cloud may contain both CO-bright gas and CO-dark gas (e.g. Bolatto et al. 2013). In molecular region with lower gas density, the H I abundance should be higher with the volume density of H I unchanging. There may also be HINSA features produced in CO-dark regions. However, HINSA features can be restrained because the cloud tends to be warmer in low-density regions. CO is important for the cooling of molecular clouds (Goldsmith & Langer 1978). For a molecular cloud with n(H₂) = 10³ cm⁻³, T_k = 10 K, and abundance of CO X_CO = 2 × 10⁻⁴, the cooling rate contributed by CO emission is (Whitworth & Jaffa 2018) (39)

The heatingrate contributed by cosmic rays and X-rays is (Glassgold & Langer 1973) (40)

If n(H₂) is smaller than 10³ cm⁻³, CO cooling is not effective enough, and the gas will be warmer and unfavorable to HINSA. The density threshold 10³ cm⁻³ is also the lower limit for the excitations of CO transitions, and this is another reason why HINSA features are confined around CO emission regions. PGCCs may be the coldest and earliest samples showing HINSA features. Overall, the non-small and CO cooling both help to confine HINSA features to regions within and around CO emission kernels.

6.5 HINSA and HC₃N

The velocity dispersions (σ) can be derived from the line widths (ΔV) through (41)

The nonthermal velocity dispersions (σ^NT) can be calculated through (42)

where m_X represents the mass of molecule X. The left panel of Fig. 17 shows the cumulative distributions of σ and σ^NT for HINSA, ¹³CO, C¹⁸O and HC₃N. σ_{HC_3N} is similar to σ_¹³CO, and tends to be larger than σ_C¹⁸O. However, the nonthermal velocity dispersion traced by HINSA is comparable with that of C¹⁸O. The broadening effect introduced by the possible thick optical depth of ¹³CO J = 1–0 cannot fully explain the larger σ^NT traced by ¹³CO compared with those traced by C¹⁸O. An optical depth of 2 contributes to line widths smaller than 30% (Phillips et al. 1979). ¹³CO may trace more extended gas component compared with C¹⁸O because ¹³CO J = 1–0 has larger optical depth compared with C¹⁸O J = 1–0 and is easier to be excited. HINSA and C¹⁸O J = 1–0 are optically thin in molecular clouds, and originate from similar regions. However, the spectra of HC₃N are quite narrow, with line widths (<0.1 km s⁻¹ ; Fig. 17) smaller than those of HINSA and ¹³CO.

Previous observations show that HC₃N is abundant in some but not all starless cores (Suzuki et al. 1992; Wu et al. 2019b). In low-mass, star formation regions, HC₃N will be enhanced again under the heating of YSOs in early stages (Sakai et al. 2008, 2009), but it may be destroyed by radiation from YSOs (Liu et al. 2021). In molecular cloud harboring young stars, it is also correlated with outflows and shocks (Wu et al. 2019b; Yu et al. 2019). The narrow line widths of HC₃N J = 2–1 imply that the HC₃N emission may originate from more compact gas component within PGCCs, instead of the outer regions dissociated by external radiation fields or shocked by turbulent flows. Although its CO emission is extended (Sect. 6.1), the PGCC could harbor dense and compact regions. The narrow line widths of HC₃N confirm that these PGCCs are in a quiescent state. Carbon chain molecule abundant PGCCs have both cold and condensed structures at the center regions traced by HC₃N and extended gas components at the outer regions traced by CO. This makes them good targets to observe and study HINSA.

The right panel of Fig. 17 shows the correlation between the column density of H I derived from HINSA and the column densities of HC₃N. There is a weak positive correlation between N(HC₃N) and N(H I), and the abundance ratio N(HC₃N)/N(H I) is ~10⁻⁵. The PGCCs with the strongest HC₃N emission tend to have more remarkable HINSA features. However, there is not enough evidence to assert that HINSA in sources with weak HC₃N emission will be suppressed. The G178.98-06.7 is one such source with no detection of HC₃N but with strong HINSA. This is not a surprise since not all PGCCs are CCM abundant, and the detection rate of HINSA is considerablein PGCCs with no previous observations of CCM emission (Tang et al. 2020). The prior of the presence of HC₃N raises the detection rate of HINSA close to ~90% in PGCCs of this work, although the sample size is still limited.

7 Summary

Twelve PGCCs, one starless core L1521B, and four star forming sources were searched for H I narrow-line self-absorption (HINSA) features using the FAST. Combining the data of H I 21 cm lines, CO J = 1–0, HC₃N J = 2–1, and dust continuum, we studied the morphologies, abundances, and excitation conditions of HINSA, CO, and HC₃N in PGCCs and star forming sources. We also investigated the relations between these parameters and the environments.

The main results are as follows:

• 1.

  HINSA features are detected in 11 of 12 PGCCs, corresponding to a detection rate of 90%. Eight of the detected ones have HC₃N emissions.

• 2.

  We improved the method of extracting HINSA features (method 2) based on the method of Krčo et al. (2008) (method 1). The S/N thresholds of the two methods are deduced. Method 2 has a much lower S/N requirement compared to method 1. Method 2 performs better when extracting HINSA features with low S/Ns, and large line widths resulting from blended components or high velocity resolution.

• 3.

  The line widths and velocities of CO lines are adopted as the initial values of the fitting parameters of HINSA. Applying the improved method, HINSA features are extracted in 11 of 12 PGCCs, the starless core L1251B, and two of four star forming sources. If method 1 was applied, HINSA features can only be extracted in L1251B and seven PGCCs.

• 4.

  The optical depths and line widths of HINSA, and the column densities of H I were obtained. The abundances of cold H I in PGCCs are approximately 3 × 10⁻⁴, varied bya factor of ~3.

• 5.

  HINSA tends not to be detected in regions associated with warm dust emission, background ionizing radiations, stellar objects, and dense gas emission.

• 6.

  The distribution of HINSA is similar to that of CO emission. External excitation by the background H I emission and the CO cooling both helps to confine HINSA features to regions within and around CO emission kernels.

• 7.

  CO emission looks as extended as the Planck dust continuum. Part of the CO emission may come from more diffuse regions besides the regions associated with the cold dust or dense gas.

• 8.

  HFS fittings are applied to HC₃N J = 2–1. The column densities of HC₃N are calculated. The average abundance of HC₃N is 2.2(±0.8) × 10⁻⁹. The nonthermal dispersions traced by HINSA are comparable with those of C¹⁸O, and larger than those of HC₃N. HC₃N emissions may originate from the condensed regions within PGCCs. CCM-abundant PGCCs provide a good sample to study HINSA.

Guided by CCM emissions, HINSA can be a good tool to detect and analyze the extended neutral gas components in PGCCs. The sample size of this work is limited. It would be helpful to search more PGCCs in the future, with the supports of the good sensitivity of the FAST and the improved HINSA extracting method.

Acknowledgements

This project was supported by the National Key R&D Program of China No. 2017YFA0402600, and the NSFC No. 12033005, 11433008, 11725313, 11373009, 11503035 and 11573036. N.-Y.T. is supported by National Key R&D Program of China No. 2018YFE0202900 and NSFC No. 11803051. K.W. acknowledges support by the National Science Foundation of China (12041305, 11973013), and the High-performance Computing Platform of Peking University through the instrumental analysis fund of Peking University (0000057511). We are grateful to the help of the staff of the FAST during the observations and data conduction. We show warm thanks to the anonymous referee for providing many constructive argues and suggestions which make the content of this work more substantial and clear.

References

All Tables

All Figures

Fig. 1 Galactic longitude-latitude positions of the 17 sources in our sample. Yellow and cyan hexagons represent the 12 PGCCs and the other five sources, respectively. The size of the hexagon represents the covering area of the 19-beam receiver of the FAST. The dot within each hexagon shows targeted locations of the central beam (Beam 1; Fig. 2). The background image shows Planck dust temperature distribution. The black contours with levels of 3, 7, and 11 Myr sr−1 represent the Planck 850 μm continuum. The red contours represent the Hα emission (Finkbeiner 2003) in units of R (106∕4π photons cm−2 s−1 sr−1), with levels from 15 to 215 stepped by 40.

In the text

	Fig. 2 Configuration of the 19-beam receiver of FAST (Li et al. 2018). Each snapshot mode observation consists of four single-point-mode observations. Equatorial coordinate system is adopted for the tracking mode, while Galactic coordinate system is adopted for the snapshot mode.
In the text

Fig. 3 Blue and orange solid line are the H I 21 cm spectra obtained by the FAST beam 1 and beam 8, respectively. The spectra werecalibrated adopting a beam efficiency of 0.85 for extended source. The red dotted lines are the spectra from the EBHIS (Winkel et al. 2016). The black solid lines are the FAST H I 21 cm spectra smoothed to have a velocity resolution similar to that of the EBHIS spectrum. To avoid overlapping, the red dotted lines and the black solid lines have been moved up by 10 and 20 K for beam 1 and beam 8, respectively.

In the text

	Fig. 4 Spectra of J = 1–0 of 12CO, C13O, and C18O (Wu et al. 2012) are shown from top to bottom in each panel.
In the text

Fig. 5 HINSA features of G174.4-15A1 extracted by method 1 (left) and method 2 (right). Method 1 fails to fit the HINSA component because of the low S/N and the blending of two absorption components with similar velocities. In each panel, the gray line is the observed spectrum Tr, the blue line is the recovered TH I, and the red line is Tab. The Tab shown here has been multiplied by a factor of five.

In the text

	Fig. 6 Same as Fig. 5, but for G170.88-10.92.
In the text

Fig. 7 H I 21 cm spectra at the center beam. In each panel, the gray line is the observed spectrum Tr, the blue line is the recovered TH I using method 2, and the red line is Tab. Values of Tab on the y-axis range from 0 to 20 K for all sources except G174.08-13.2, which uses 0–60 K. The orange and blue squares represent the sources that can and cannot be successfully fit using method 1. The blue line is not shown for spectrum with fit τ smaller than 0.01.

In the text

Fig. 8 Distribution of S/Nc and S/N(Tab) for the H I spectra of the central beams (Table 3). Here, S/N(Tab) is the ratio between the intensity of absorption dip revealed by method 2 and the spectra noise, and S/Nc is the signal-to-noise threshold required by method 1 (Sect. 5.2.2). The orange and blue squares represent the sources that can and cannot be successfully fit using method 1, respectively. The black line represents x = y. As predicted, the black line is the separation line between the two cases.

In the text

Fig. 9 Backgrounds: strengths of HINSA measured in the snapshot mode, extracted using method 2. Black contours: 13CO J = 1–0 emissions, which have been smoothed to the angular resolution of H I observations. Big black boxes: areas covered by CO data. Blue crosses: sampled points of the FAST, and red circle: center beam. Small red boxes and triangles: Class I/II and Class III YSOs. In panels a and c, the white contours represent the emission of C2H N = 1–0 within the white rectangles (Liu et al. 2019). In panel f, the small black rectangle represents the mapped region of the IRAM HCO+ J = 1–0 and C18O J = 2–1 observations, which are shown in Fig. 14.

In the text

	Fig. 10 Backgrounds: strengths of HINSA measured in the tracking mode. Symbols have the same meanings with those in Fig. 9.
In the text

	Fig. 11 HC3N spectra binned with a velocity resolution of 0.05 km s−1. The red lines show the results of HFS fitting.
In the text

Fig. 12 Relations between the column densities of H2 derived from CO and dust. Upper: correlation between the column densities of H2 derived from13CO (N13CO(H2)) and C18O (NC18O(H2)) at the positions of the FAST beams for 11 PGCCs (G174.08-13.2 is not included because of lacking CO mapping observations). Data of different sources are shown in different colors. Lower: correlation between the H2 column densities at the map centers derived from dust emissions (Ndust(H2)) and C18O (NC18O(H2)).

In the text

	Fig. 13 Comparison between CO J = 1–0 excitation temperatures and Planck dust temperatures. The black line represents y = x.
In the text

Fig. 14 Cutout region of G192.2-11A2 represented by the black box in Fig. 9f. The background and black contours present the emissions of C18O J = 2–1 and HCO+ J = 1–0 measured with the IRAM 30 m telescope (Xu et al., in prep.), respectively. The contours have levels from 0.5 to 0.9 stepped by 0.1 in unit of K km s−1. The red circle shows the location of the central beam of the H I observation. The yellow dashed circle and solid circle represent the beam size of HCO+ J = 1–0 and C18O J = 2–1, respectively. The white dashed lines are the grid lines of the Galactic coordinate.

In the text

	Fig. 15 Correlation between the column densities of H I and those of H2 derived fromC18O J = 1–0 emission. The orange squares represent the data of the central beams. The blue dots represent the data of other beams. The red square shows the data at the peak position of the HINSA map of G165.6-09A1.
In the text

	Fig. 16 Blue line: average H I spectrum of the observed sources in this work. Orange line: estimated bymultiplying the blue line by sin(14°). Green line: evaluated from the full-sky H I survey HI4PI based on EBHIS and GASS (HI4PI Collaboration 2016).
In the text

Fig. 17 Comparisons between the velocity dispersions and column densities of HC3N and H I. Left: cumulative distributions of velocity dispersions. Solid lines and dashed lines represent σ and σNT, respectively (Sect. 6.5). Right: correlation between the column densities of H I and those of HC3N. Green and orange lines represent the linear fitting with and without the two data in blue shadowing, respectively. For data in the blue shadowing, the N(H I) parameters are not reliable with regard to the fit optical depth of HINSA τ < 0.01 (Sect. 5.2.5).

In the text