Neutral Stellar Winds Toward the High-Mass Star-Forming Region G176.51+00.20

We observed the high-mass star-forming region G176.51+00.20 using the Five-hundred-meter Aperture Spherical radio Telescope (FAST) with the 19-beam tracking observational mode. This is a pilot work of searching for neutral stellar winds traced by atomic hydrogen (i.e., HI winds) using the high sensitivity HI line toward high-mass star-forming regions where bipolar molecular outflows have been detected with high sensitivity by Liu et al. HI wind was detected in this work only in Beam 1. We find here that, similar to low-mass star formation, no matter how large the inclination is, the HI wind is likely sufficiently strong to drive a molecular outflow. We also find that the abundance of HI in the HI wind is consistent with that of the HI narrow-line self-absorption (HINSA) in the same beam (i.e., Beam 1). This implies that there is probably an internal relationship between HI winds and HINSA. This result also reinforces the assertion that HI winds and detected molecular outflows are associated with each other.

One of the important goals of contemporary astrophysics is to understand star formation (e.g., Tan et al. 2014;Bally 2016). Relative to low-mass star formation, the formation mechanism(s) of high-mass stars remain poorly understood (Shu et al. 1987;McKee & Ostriker 2007;Tan et al. 2014), as too are the driving scenario(s) of massive molecular outflows, an essential phase of early high-mass star formation (Arce et al. 2007).
Early studies, especially in the 1980s and 1990s, attempted to ascertain the driving source and mechanism of molecular outflows, mostly for low-mass stars.  and Strom et al. (1986) found that the mass-loss rate of ionized gas is too low to drive CO lobes, and suggested that the stellar winds which drive the CO lobes are largely comprised of neutral gas. Lizano et al. (1988) carefully compared the neutral stellar wind traced by atomic hydrogen (H i wind hereafter) detected in HH 7-11 and the corresponding CO/HCO + outflows, and found that the H i wind was strong enough to drive the molecular outflows. A similar conclusion was drawn by carefully comparing the H i winds in both HH 7-11 and L1551 with the corresponding molecular outflows (see Giovanardi et al. 1992). In such studies, which were performed with the 305 m radio telescope of the Arecibo Observatory, the root mean square (rms) noise was typically ∼ 3.5-8.5 mK @ 3.9 km s −1 with beam size of 3 .2-3 .5. A model based on the detected H i wind in HH 7-11 was developed by Lizano et al. (1988) to explain the entrainment of ambient molecular gas (i.e., driving the molecular outflow) by the neutral wind. These studies suggested that H i winds are a promising driving source of molecular outflows. However, is such an assertion also true for massive stars?
To answer this question, samples of the molecular outflows have been assembled. For example, main beam rms noise of dozens of mK) toward nine high-mass star-forming regions using the 13.7 m millimeter telescope of the Purple Mountain Observatory in Delingha (see more details in Liu et al. 2021). Accordingly, an rms noise of ∼ 3 mK @ 1.0 km s −1 can be achieved with about six hour integrations and a beam size of ∼ 2 .9 by using the powerful 21 cm H i line detected by the most sensitive ground-based, single-dish Five-hundred-meter Aperture Spherical radio Telescope 1 (FAST; Nan 2006;Nan et al. 2011;Jiang et al. 2020). The use of FAST will enabled us to conduct broad investigations of stellar winds toward star-forming regions where molecular outflows have been detected. As such, we have carried out a pilot survey of the H i wind toward detected molecular outflows with highly sensitive observations. The selected object in this work is one of the nine high-mass star-forming regions in Liu et al. (2021), G176.51+00.20.
G176.51+00.20, also known as AFGL 5157 and IRAS 05345+3157, is located 1.8 kpc from Earth (Moffat et al. 1979;Snell et al. 1988). It is a complex region containing multi-generational star formation (e.g., Chen et al. 2003;Dewangan 2019). The main research object in this work is a highmass star-forming region at a relatively young phase, where H 2 O masers and an H ii region have been detected (Torrelles et al. 1992a). This region is marked as an "H ii region" in figure 1(a) in Dewangan (2019). There is a dense NH 3 core in the center of this region, which has become synonymous with this region (e.g., Torrelles et al. 1992a;Chen et al. 2003;Jiang et al. 2013). By constructing a 3.6 cm map with the Very Large Array, the excitation source of this region was identified as a zero-age main-sequence B3 star (Torrelles et al. 1992b). The bolometric luminosity of this excitation source was determined as 1.7 × 10 3 L via an 8-1200 µm SED fit (see Molinari et al. 2008). The age of this star-forming region was suggested to be ∼ 2 × 10 5 yr by investigating the near-infrared H 2 line emission using the 1.88 m telescope of Okayama Astronomical Observatory, Japan (see Chen et al. 2003).
The famous bipolar outflow (traced by 12 CO) detected by Snell et al. (1988) is centered on the dense NH 3 core (Molinari et al. 2002;Zhang et al. 2005;Dewangan 2019), suggesting that the excitation source of this bipolar outflow is probably the zero-age main-sequence B3 star (see Torrelles et al. 1992b). Liu et al. (2021) confirmed this bipolar outflow with high-sensitivity 12 CO, 13 CO, C 18 O, HCO + , and CS line emission observations, and broadened the blue lobe of this outflow (traced by 12 CO) by a factor of ∼ 2. This makes the dense NH 3 core more suitable to investigate whether its H i wind is strong enough to drive the molecular outflows.
The remainder of this paper is organized as follows. In Section 2, we describe the data used in this work. Section 3 analyses the H i wind associated with the molecular outflows. In Section 4, discussion of whether the H i wind is strong enough to drive the molecular outflows is presented.
Finally, Section 5 gives a summary and the main conclusions of this work.

H I OBSERVATIONS WITH FAST
FAST is located in Guizhou Province of southwest China (Nan et al. 2011). It is equipped with a 19-beam receiver (the frequency range is 1.0-1.5 GHz with a bandwidth of 500 MHz) and has dual linear polarizations (i.e., XX and YY; see Li et al. 2018;Jiang et al. 2019Jiang et al. , 2020. The spectral resolution is ∼ 477 Hz, corresponding to a velocity resolution of ∼ 0.1 km s −1 at 1.4 GHz. The half-power beam width (HPBW) is ∼ 2.9 at 1.4 GHz, and the pointing error is ∼ 0.2 (see Jiang et al. 2019Jiang et al. , 2020. Observations toward G176.51+00.20 were conducted on August 19th and 20th, 2021. The observational mode was 19-beam tracking, and the total integration time was 335 minutes with a sampling rate of one second. A new algorithm is used to obtain a flat baseline. In this algorithm, there are three main steps: (1) a polynomial fitting over the velocity range of [−7000, 7000] km s −1 ; (2) removing the standing waves by using fast Fourier transforms; (3) calibrating the unsmooth parts in the baseline by extreme envelope curves (for more details of the pipeline used to obtain the flat baseline, see Liu et al. 2022). For Beams 2, 3, 5, 6, and 8-19, the spectra of two linear polarizations (i.e., XX and YY) were dealt with separately and then averaged. However, because of the bad baseline, only the YY polarization spectrum from Beam 1 and the XX spectra from Beams 4 and 7 were retained.
The mean rms noise of the nineteen spectra is ∼ 7 mK @ 0.1 km s −1 (see the rms noise of each spectrum in Table 1). Figure 1 presents the positions of the central seven beams superposed on CO molecular maps, and Figure 2 shows the spectrum of the central seven beams and the H i narrow-line self-absorption lines (HINSA, T ab ) in each beam (see ).

Identification of the H i Wind
We defined the excess brightness temperature (EBT, hereafter) as the spectrum resulting from the subtraction of the average spectrum of the off-positions from the spectrum of the on-position (i.e., the target beam). Similar to the work of Lizano et al. (1988), the six beams located immediately around the target beam (i.e., which acted as the on-position) were used as the off-positions. The baseline of the EBT was flat, i.e., the rms noise was ∼ 1.8-2.3 mK @ 2.0 km s −1 (i.e., smooth over twenty channels), except for Beams 3 and 4, whose rms noises were, respectively, 3.3 and 6.1 mK @ 2.0 km s −1 . We also defined the standard deviation (STD) of EBT (SEBT) to evaluate the variation among the spectrum from different off-positions. For instance, for Beam 1, we first separately calculated the EBT of Beam 1 relative to its six surrounding beams, and then computed the STD of these six EBTs, which is the SEBT of Beam 1.

H i Wind in Beam 1
From Figure 3, it can be seen that Beam 1 presents an evident high-velocity red wing whose maximum velocity exceeds ∼ 120 km s −1 . The SEBT is also small in the high-velocity red wing, indicating this wing may well be true, rather than the changes or fluctuations of the emission from one beam area to another. This high-velocity red wing was crudely fitted by a "triangular fit" (see e.g., Lizano et al. 1988). Assuming this high-velocity line wing is bipolar, the symmetric line of the high-velocity red wing with respect to the line center, v c,w , is also presented as a red oblique line (see the red isosceles triangle of the high-velocity bipolar wing in Figure 3(a)). Therein, v c,w is taken as the velocity of the corresponding HINSA (see Figure 2 and , which is indicated as a red vertical line in Figure 3(a). This symmetric line is referred to as the presumed high-velocity blue wing, although there is no blue wing present in Figure 3(a). However, we find that the velocity of the bulge in Beam 7 (see Figure 3(g)) overlaps with that of the presumed high-velocity blue wing in Beam 1. Therefore, a further analysis of the high-velocity line wing in Beam 1 was conducted where we took different beams as the off-positions (see Figure 4). five cases are referred to as Normal, W-E, NE-SW, NW-SE, and Lack Beam 7. The high-velocity red wing remains basically unchanged no matter which beams are used as off-positions (see Figure 4 and quantitative comparison in Section 3.2), and the high-velocity blue wing is also present as long as the off-positions do not contain Beam 7. Therefore, Beam 1 probably presents a high-velocity bipolar wing, but the high-velocity blue wing is contaminated by other components. In addition, there is a small pit between the "triangular fit" of the EBT for the cases of W-E, NE-SW, and NW-SE (see the left green oblique line in Figures 4(a)-(c)) and the presumed high velocity blue wing (see red oblique line). This small pit likely causes the green oblique line to deviate from the red oblique line, which implies that this high-velocity bipolar wing is probably symmetrical; therefore, we identify it as a bipolar H i wind. Because the blue lobe is contaminated by other components, and we presume that this H i wind is a symmetrical bipolar H i wind, we applied the fitting result of the red lobe to the blue lobe.   The bulge in the blue line wing of Beam 2 is, in fact, present in only one polarization spectrum (i.e., YY). Therefore, this blue line wing is highly likely a fake high-velocity line wing. Figure 3 also presents other high-velocity features, e.g., wings with negative intensities in right side of EBTs in Beams 2 and 7, and wings with positive intensities in right side of EBTs in Beam 6, etc. These high-velocity features are more likely to stem from the changes or fluctuations of the emission from one beam area to another, because these features are located where SEBTs are large.

Physical Parameters of the H i Wind
We calculated the mass, M HI,wind , momentum, P HI,wind , maximum velocity, V max,wind , and average velocity, V mean,wind , of the red lobe of the H i wind ("triangular fit") corresponding to different offpositions, i.e., the case of Normal (see red line in Figure 3 Table 2). This reinforces the conclusion that the red lobe of the H i wind is real, and the case of Normal can reflect the physical properties of the red lobe of the H i wind. Therefore, the calculations below are based on the case of Normal.
The "triangular fit" of the red lobe of the H i wind produced T wind (v) = (14.94 ± 1.32) + (−0.15 ± correction are, respectively, 120.6 ± 16.3 and 40.2 ± 5.4 km s −1 , where V mean,wind = V max,wind /3 for the triangular profile. V max,wind (and also V mean,wind ) are 0.7 and 0.8 times that in HH 7-11 and in L1551 (both are located in low-mass star-forming regions; see Lizano et al. 1988;Giovanardi et al. 1992), respectively. The column density of the red lobe of the H i wind, N HI,wind , is (1.94 ± 0.30) × 10 18 cm −2 (see Table 3). Assuming that the H i wind occupies the same spatial location as the molecular outflows, and both of them are angularly unresolved (e.g., Lizano et al. 1988), N HI,wind is two orders of magnitude smaller than the column density of the outflow traced by 12 CO, and three orders of magnitude smaller than that traced by 13 CO, HCO + , and CS (see the column density of the molecular outflows in Liu et al. 2021).
13 CO is proved to be a great tracer of the molecular column density, and the fractional HINSA abundance (where the molecular column density is traced by 13 CO) is ∼ 1.1 × 10 −3 .
We also calculated the abundance of H i in the H i wind, X HI,wind , as: where N H 2 ,outflow is the column density of the molecular outflow (that is traced by 12 CO, 13 CO, HCO + , or CS; see Liu et al. 2021). The results of X HI,wind are listed in Table 3. X HI,wind , determined by the molecular outflow traced by 13 CO, is ∼ 1.0 × 10 −3 (only for the red lobe), which is consistent with the fractional HINSA abundance.
This consistency reveals an internal correlation between the H i wind and HINSA. It also supports the conclusion of the mixture of the H i wind and the molecular outflow because HINSA is probably mixed with the gas in the cold, well-shielded regions of the molecular clouds (see Li & Goldsmith 2003;. However, observations with higher resolution are required to confirm the correlation between the H i wind and HINSA, and a larger sample is essential to rule out that the observed consistency here is only a coincidence. We also calculated the mass, M HI,wind , momentum, P HI,wind , and kinetic energy, E HI,wind , of the H i wind (combine the values of the two lobes) without inclination correction (see the calculations in Section A and the values in Table 3). The mass and momentum over the overall flow history are larger than the values estimated above (see Lizano et al. 1988). The values of M HI,wind , P HI,wind , and E HI,wind are larger than that in HH 7-11 (see Lizano et al. 1988;Giovanardi et al. 1992) by a factor of ∼ 7, ∼ 5, and ∼ 4, respectively, and larger than that in L1551 (see Giovanardi et al. 1992) by a factor of ∼ 26, ∼ 20, and ∼ 6, respectively.

Dynamical Comparison Between the H i Wind and the Molecular Outflow
To determine whether the H i wind is strong enough to drive the molecular outflow, analyses of two dimensions were conducted and are presented here; i.e., the mass loss of the excitation star and the lifetime of the H i wind.

Mass Loss of the Excitation Star
If we assume that the molecular outflow is entrained by the H i wind originally ejected from the star, then the mass loss, ∆M * ,m , from the excitation star over the entire history of mass loss is M * ,m = P m,outflow /V max,wind (see Lizano et al. 1988), where P m,outflow is the momentum of the molecular outflow (see Liu et al. 2021). Both the blue and red lobes of the molecular outflow are included in the calculation. The total ∆M * ,m , which combines contributions from the 12 CO, 13 CO, HCO + , and CS outflows, is 6.26 ± 0.80 M (see the contribution from each molecular outflow in Table 3). Because Liu et al. (2021) did not give the error of P m,outflow , the error of ∆M * ,m is only from V max,wind and no error is presented in Table 3. It is worth noting that the total ∆M * ,m is probably overestimated because the values of P m,outflow are corrected for the abundance of the corresponding tracer (see Liu et al. 2021).
Similar to the work in Lizano et al. (1988), we calculated the atomic mass excesses, ∆M * ,ame , of the line core of Beam 1 relative to the line core of the six beams tightly around Beam 1, assuming that the low-velocity material is associated with the source and fills the beams of FAST. We tested different velocity spans, v s , as line core (see Figure 5 Beam 1, and these atomic mass excesses may be an indicator of accumulated atomic mass loss from the excitation star, ∆M * ,a (e.g., Lizano et al. 1988). The atomic mass excesses of Beam 1 relative to its surrounding beams with a velocity range of |v c,w − v| ≤ 20 km s −1 are listed in Table 3, which are set to be two times of the corresponding atomic mass excesses with a velocity range of v c,w − v ≤ 20 km s −1 . The average value of ∆M * ,ame over the five atomic mass excesses (i.e., except for that between Beams 1 and 3, where the atomic mass excesses are significantly different from the other atomic mass excesses) is 6.0 ± 3.5 M (i.e., ∆M * ,a ∼ 6.0 ± 3.5 M ), where 3.5 M is the STD of the five atomic mass excesses. The error of atomic mass excess inherited from the original spectrum is far less than 3.5 M , and therefore no error of atomic mass excess is presented in Table 3. ∆M * ,a ∼ ∆M * ,m , and the outflowing atomic gas would undergo deceleration when it entrained molecular outflow gas, indicating that the velocity associated with ∆M * ,a may be not less than that associated with ∆M * ,m and H i wind is likely strong enough to drive the molecular outflow.

Lifetime of the H i Wind
The lifetime of the H i wind flow, i.e., the total duration of the flow, t, without inclination correction was estimated according to Lizano et al. (1988). The distance from the H i wind to the central excitation star, l HI,wind , was set as the outflow length, l outflow , of the molecular outflow (which, here, corresponds to the red lobe of 12 CO; see Liu et al. 2021), i.e., l HI,wind ∼ 0.8 pc. The crossing time of the H i wind, t a , is then t a = l HI,wind /V mean,wind = (1.9 ± 0.2) ×10 4 yr, where the uncertainty comes only from V mean,wind because no error of l outflow was given by Liu et al. (2021). Therefore, the stellar mass-loss rate,Ṁ a , isṀ a = M HI,wind /t a = (12.2 ± 3.6) × 10 −6 M yr −1 . This value (without inclination correction) is larger than that in HH 7-11 by a factor of ∼ 2 and larger than that in L1551 by a factor of ∼ 10 (see Lizano et al. 1988;Giovanardi et al. 1992). The lifetime of the H i wind is t = ∆M * ,a /Ṁ a = (5.3 ± 3.5) × 10 5 yr, where the accumulated atomic mass, ∆M * ,a , is set to 6.0 ± 3.5 M (see above). All the physical properties above are listed in Table 3. The dynamical timescale of the red lobe of the molecular outflow, t outflow , as an indicator of the age of the flow, ranges from 1-5 ×10 4 yr traced by four species (i.e., 12 CO, 13 CO, HCO + and CS; see Liu et al. 2021). t is about one order of magnitude larger than t outflow , indicating that the H i wind is likely strong enough to drive the molecular outflow. Observations with higher spatial resolution of both H i and molecules are required to further confirm this conclusion.
To further investigate the influence of the inclination, Figure 6 presents both t and t outflow (including those for the 12 CO, 13 CO, HCO + , and CS outflows) under different inclinations. It can be seen that the lifetime of the H i wind is much larger than the timescale of the outflow traced by all the four molecules, no matter what inclination is considered. This suggests that the inclination does not change the inequality of t t outflow .

Physical Properties with an Inclination of 69 •
We next calculated the physical properties of the H i wind assuming that an inclination of 69 • (see Table 4). The maximum velocity of the H i wind changes to ∼ 336.2 ± 44.6 km s −1 . This value is larger than that for HH 7-11 (assuming an inclination of 22 • ; see Herbig & Jones 1983;Lizano et al. 1988) by a factor of ∼ 2, and is ∼ 60% of that for L1551 (assuming an inclination of 75 • ; see Snell & Schloerb 1985;Giovanardi et al. 1992). The value of the momentum of the H i wind after inclination correction is 22.66 M km s −1 , which is larger than that in HH 7-11 by a factor of ∼ 12 (see Lizano et al. 1988;Giovanardi et al. 1992), and a factor of ∼ 15 for that in L1551 (see Giovanardi et al. 1992). The values of the wind's kinetic energy, i.e., ∼ 2.6 × 10 47 erg, and the stellar mass-loss rate, i.e., ∼ 3.2 × 10 −5 M yr −1 , are larger than those in HH 7-11 by factors of ∼ 28 and ∼ 10, and larger than those in L1551 by factors of ∼ 3 and ∼ 7, respectively (see these values in HH 7-11 and L1551 in Lizano et al. 1988;Giovanardi et al. 1992).
The kinetic luminosity of the H i wind, L HI,wind = 1 2Ṁ a V 2 max,wind , is 16.0 ± 10.1 L without inclination correction, which is much less than the bolometric luminosity, L * , of the excitation source (i.e., Note-See description of each quantity in Table 3. 1.7 × 10 3 L ; see Molinari et al. 2008). However, the kinetic luminosity of the H i wind after inclination correction is 322.0 L , i.e., ∼ 19% of L * . The ratio of L HI,wind /L * after inclination correction is comparable to that in HH 7-11 (i.e., 20%) and in L1551 (i.e., 127%, see Giovanardi et al. 1992), while the ratio of L HI,wind /L * before inclination correction is only ∼ 1%. This result implies that an inclination correction is necessary when we calculate the physical properties of the H i wind.

SUMMARY AND CONCLUSIONS
We obtained high-sensitivity H i spectra with an average rms noise of ∼ 7 mK @ 0.1 km s −1 toward the high-mass star-forming region G176.51+00.20 using FAST with a 19-beam tracking observational mode. We searched for H i wind in the central seven beams (i.e., Beams 1-7), but detected it only in Beam 1. The main results are as follows: 1. The abundance of H i in the H i wind is consistent with that of the HINSA (where the molecular column density is traced by 13 CO). This indicates that there probably exists an internal correlation between the H i wind and HINSA; such a correlation would enhance the argument of the association between the H i wind and the molecular outflows.
2. The H i wind is likely strong enough to drive the molecular outflows. This conclusion is not affected by the value of the inclination.
3. The mass, momentum, and kinetic energy of the H i wind and the associated mass-loss rate of the excitation star in G176.51+00.20, a high-mass star-forming region, are about one orders of magnitude larger than those in low-mass star-forming regions (i.e., HH 7-11 and L1551).

APPENDIX
A. MASS AND MOMENTUM OF THE H I WIND A crude "triangular fit" was adopted to obtain the relationship between the brightness temperature and the radial velocity (e.g., Lizano et al. 1988) as where a and b are constants with units of K and K (km s −1 ) −1 , respectively. The result of the "triangular fit" is plotted as the red solid line in Figure 3(a). Considering that the H i wind may be contaminated or affected by the ambient environment, and manifesting as an asymmetry between the blue and red lobes, the blue and red lobes were fitted separately with the line center, v c,w , as the velocity of HINSA (see , see also Figure 2 and the red vertical line in Figure 3(a)).
T wind can be used to calculate the column density of the H i wind as: N HI,wind = 1.82 × 10 18 T wind dv, where H i emission is assumed to be optically thin (Dickey & Benson 1982;Lizano et al. 1988;Saha et al. 2018).
The maximum velocity of the H i wind, V max,flow , reads: and the mean velocity of the H i wind is V mean,flow = V max,flow /3 for the triangular profile. Considering that the beam size of FAST is 2 .9 and the distance to the H i wind is 1.8 kpc, the mass, M HI,wind , and momentum, P HI,wind , of the H i wind without inclination correction are: and P HI,wind = C 0 |v − v c,w |T wind dv, respectively, where the constant C 0 is: where we assume that the hydrogen mass abundance is 0.74 (see Garden et al. 1991) and the H i wind emission is from a point source (e.g., Lizano et al. 1988). This is reasonable because the value of V max,flow is 120.6 ± 16.3 km s −1 (see Lizano et al. 1988;Frank et al. 2014, and references therein).
The value of C 0 is larger than that in Equation (A6) by a factor of 2 ln 2 for an extended source with a uniform brightness temperature over the beam. The kinetic energy, E HI,wind , of the H i wind is (see Giovanardi et al. 1992): E HI,wind = 1 2 M HI,wind V 2 max,flow .
This work made use of the data from FAST. FAST is a Chinese national mega-science facility, operated by National Astronomical Observatories, Chinese Academy of Sciences. We would like to thank the anonymous referee for the helpful comments and suggestions that helped to improve the paper. This work was sponsored by the Natural Science Foundation of Jiangsu Province (grant No.