# H i Narrow-Line Self-Absorptions Toward the High-Mass Star-Forming Region
G176.51+00.20

Yingjie Li  Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing
210008, China  Ye Xu  Purple Mountain Observatory, Chinese Academy of
Sciences, Nanjing 210008, China  Jin-Long Xu  National Astronomical
Observatories, Chinese Academy of Sciences, Beijing 100101, China  CAS Key
Laboratory of FAST, National Astronomical Observatories, Chinese Academy of
Sciences, Beijing 100101, China  Dejian Liu  Purple Mountain Observatory,
Chinese Academy of Sciences, Nanjing 210008, China  University of Science and
Technology of China, Hefei, Anhui 230026, China  Jingjing Li  Purple Mountain
Observatory, Chinese Academy of Sciences, Nanjing 210008, China  Zehao Lin
Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210008,
China  University of Science and Technology of China, Chinese Academy of
Sciences, Hefei, Anhui 230026, China  Peng Jiang  National Astronomical
Observatories, Chinese Academy of Sciences, Beijing 100101, China  CAS Key
Laboratory of FAST, National Astronomical Observatories, Chinese Academy of
Sciences, Beijing 100101, China  Shuaibo Bian  Purple Mountain Observatory,
Chinese Academy of Sciences, Nanjing 210008, China  University of Science and
Technology of China, Hefei, Anhui 230026, China  Chaojie Hao  Purple Mountain
Observatory, Chinese Academy of Sciences, Nanjing 210008, China  University of
Science and Technology of China, Hefei, Anhui 230026, China  Xiuhui Chen
College of Mathematics and Physics, Hunan University of Arts and Science,
Changde, Hunan 415300, China

###### Abstract

Using the Five-hundred-meter Aperture Spherical radio Telescope (FAST) 19-beam
tracking observational mode, high sensitivity and high-velocity resolution H i
spectral lines have been observed toward the high-mass star-forming region
G176.51+00.20. This is a pilot study of searching for H i narrow-line self-
absorption (HINSA) toward high-mass star-forming regions where bipolar
molecular outflows have been detected. This work is confined to the central
seven beams of FAST. Two HINSA components are detected in all seven beams,
which correspond to a strong CO emission region (SCER; with a velocity of
$\sim$ $-$18 km s-1) and a weak CO emission region (WCER; with a velocity of
$\sim$ $-$3 km s-1). The SCER detected in Beam 3 is probably more suitably
classified as a WCER. In the SCER, the HINSA is probably associated with the
molecular material traced by the CO. The fractional abundance of HINSA ranges
from $\sim 1.1\times 10^{-3}$ to $\sim 2.6\times 10^{-2}$. Moreover, the
abundance of HINSA in Beam 1 is lower than that in the surrounding beams
(i.e., Beams 2 and 4–7). This possible ring could be caused by ionization of H
i or relatively rapid conversion from H i to H2 in the higher-density inner
region. In the WCER (including Beam 3 in the SCER), the HINSA is probably not
associated with CO clouds, but with CO-dark or CO-faint gas.

ISM: jets and outflows - ISM: Molecules - stars: formation – ISM: abundances –
ISM: kinematics and dynamics

##  1 Introduction

Stars form in molecular clouds, and it is important for star formation to
synthetize molecular H2 from atomic H. However, it is difficult to directly
extract individual components from H i emission and correlate them with
molecular clouds (Li & Goldsmith, 2003; Liu et al., 2022). In such conditions,
H i self-absorption (HISA) lines have become a very useful tracer of H i. HISA
occurs in environments where cold H i gas is in the front of a warmer emission
background (e.g. Knapp, 1974); such conditions are ubiquitous in the Milky Way
(e.g., Heeschen, 1955; Gibson et al., 2000; Wang et al., 2020). HISA is
associated with both CO-dark clouds and molecular clouds with strong CO
emission (e.g., Hasegawa et al., 1983; Gibson et al., 2000; Gibson, 2010). H i
narrow-line self-absorption (HINSA), a special case of HISA, is usually
associated with CO clouds, and its linewidth is comparable or smaller than
that of CO (see Li & Goldsmith, 2003; Goldsmith & Li, 2005; Tang et al., 2016;
Wang et al., 2020). HINSA has been found to be an excellent tracer of
molecular clouds (Li & Goldsmith, 2003), and this tight correlation has been
confirmed by numerous studies (e.g., Krčo & Goldsmith, 2010; Tang et al.,
2020; Liu et al., 2022).

Numerous surveys have provided abundant knowledge about the environment and
physical properties of the regions producing HINSA. For instance, the
detection rate of HINSA is $\sim$ 77% (Li & Goldsmith, 2003) in optically
selected dark clouds (Lee & Myers, 1999), $\sim$ 80% (Krčo & Goldsmith, 2010)
in molecular cores (i.e., the Lynds dark cloud, which has a peak optical
extinction of 6 mag or more; see the catalog of dark clouds in Lynds, 1962),
and $\sim$ 58% (Tang et al., 2020) to $\sim$ 92% (Liu et al., 2022) in Planck
Galactic cold clumps (PGCCs; cold and quiescent clumps in very early evolution
stages of star formation; see Planck Collaboration et al., 2011; Wu et al.,
2012; Planck Collaboration et al., 2016). These results indicate that HINSA is
associated with molecular clouds, especially molecular cores and/or clumps.
These papers also constrained the abundance of HINSA, which ranges from $\sim$
10-4 to $\sim$ 10-2 (e.g., Krčo & Goldsmith, 2010; Tang et al., 2020; Liu et
al., 2022); the relationship of the central velocity, linewidth, and column
density between HINSA and the molecular gas traced by 12CO, 13CO, OH, etc.
(e.g., Li & Goldsmith, 2003; Tang et al., 2020; Liu et al., 2022) and the
atomic gas traced by CI (e.g., Li & Goldsmith, 2003); and even a ring of
enhanced HINSA abundance inside of a dark molecular cloud (i.e., in B227; see
Zuo et al., 2018).

Star-forming regions characterized by molecular outflows have also been
incorporated into samples arising from the search for HINSA (e.g., Liu et al.,
2022). In this study, the corresponding detection rate of HINSA was only 25%
(i.e., one in four sources characterized by molecular outflows). As such, one
of our goals is to enlarge the sample of known star-forming regions
characterized by molecular outflows. A good choice is the nine regions where
outflowing gases traced by 12CO, 13CO, HCO+, and CS have been recently
detected with high sensitivity (i.e., main beam rms noise of tens of mK) with
the 13.7 m millimeter telescope of the Purple Mountain Observatory in Delingha
by Liu et al. (2021). In one of the nine regions, i.e., G176.51+00.20, a
neutral stellar wind traced by atomic hydrogen (H i wind) has also been
detected (Y. J., Li et al., 2022, ApJ, accepted). Therefore, studying HINSA
toward G176.51+00.20 not only would be a bridge to connect the atomic and
molecular gas therein, but may also be helpful for studying the relationship
between H i winds and molecular outflows.

G176.51+00.20, an active high-mass star-forming region, is located 1.8 kpc
from Earth (Moffat et al., 1979; Snell et al., 1988). The dense NH3 core in
the center of this region has become synonymous with this region (e.g.,
Torrelles et al., 1992a; Chen et al., 2003; Jiang et al., 2013). The central
engine of this massive star-forming region was identified as a zero-age main-
sequence B3 star by using a 3.6 cm map produced from Very Large Array
observations (Torrelles et al., 1992b). A more detailed description of this
region was provided by Dewangan (2019), where the center of the region of
study in this work was marked as “H ii region” in their figure 1(a).

The remainder of this paper is organized as follows. In Section 2, we describe
the data used in this work. Section 3 presents the results of the HINSA survey
toward G176.51+00.20, and the physical properties of the HINSA features and
those of the corresponding molecular lines. In Section 4, we discuss the
association between the HINSA features and the molecular clouds, and estimate
the abundance of the HINSA features, and summarize our results.

##  2 Data

###  2.1 H i Observations with FAST

The most sensitive ground-based single-dish radio telescope is the Five-
hundred-meter Aperture Spherical radio Telescope111https://fast.bao.ac.cn/,
which located in Guizhou Province of southwest China (Nan, 2006; Nan et al.,
2011; Qian et al., 2020). By using the 19-beam receiver equipped on FAST (see
Li et al., 2018; Jiang et al., 2019, 2020), we observed highly sensitive and
high-velocity resolution H i spectra toward G176.51+00.20. The mean main beam
rms noise for the 19 obtained spectra is $\sim$ 7 mK @ 0.1 km s-1 resolution.
The observations were conducted on August 19th and 20th, 2021, with total
integration time is 335 minutes. The half-power beam width (HPBW) is $\sim$
2.9′ at 1.4 GHz, and the pointing error is $\sim$ 0.2′ (see Jiang et al.,
2019, 2020).

###  2.2 Summary of Molecular Lines Observed Toward G176.51+00.20

Five high-sensitivity molecular lines were observed toward G176.51+00.20 (with
a cell of $\sim$ $15^{\prime}\times 15^{\prime}$) by using the 13.7 m
telescope: 12CO ($J=1-0$)(115.271 GHz), 13CO ($J=1-0$)(110.201 GHz), and C18O
($J=1-0$)(109.782 GHz) observed from 2019 July to November; HCO+ ($J=1-0$)
(89.189 GHz) observed from 2019 November to 2020 February; and CS
($J=2-1$)(97.981 GHz) observed in 2020 May (for more details, see tables 2–3
of Liu et al., 2021). The velocity resolution of these lines ranges from 0.159
to 0.212 km s-1 (the corresponding channel width is 61 kHz), and the HPBW
ranges from 49′′ to 61′′. The lines were observed with the on-the-fly mode,
and gridded to 30′′ pixels. The rms noises range from $\sim$ 60.9 to $\sim$
15.3 mK (for more details see Liu et al., 2021).

##  3 Data Analysis and Results

Because the molecular line observations are confined to the central $\sim$
$15^{\prime}\times 15^{\prime}$ region centered at (R.A., Decl.) (J2000) =
(05h37m52s.1, 32∘00′03′′.0), the HINSA survey is also limited to the
corresponding regions; i.e., central seven beams (Beams 1–7; see Figure 1(a)).
We followed the method of Liu et al. (2022) to extract HINSA features. To
extract HINSA, it was required to understand the physical properties of the
emission regions, at least the central velocity, $v_{\mathrm{c,m}}$, velocity
dispersion, $\sigma_{\mathrm{m}}$, and excitation temperature,
$T_{\mathrm{ex}}$ (i.e., the initial conditions to extract HINSA; see Liu et
al., 2022). These three physical parameters were obtained by using the high-
sensitivity molecular line observations of the high-mass star-forming region
G176.51+00.20 made by Liu et al. (2021).

###  3.1 Overview of the Region with High-Sensitivity Molecular-Line
Observations

From the study of molecular outflows toward G176.51+00.20 by Liu et al.
(2021), the velocity of the most conspicuous components is $\sim$ $-$18 km
s-1, which corresponds to a high-mass star-forming region with a strong
bipolar outflow (see Snell et al., 1988; Zhang et al., 2005; Liu et al.,
2021). In fact, there is also a weak component with a velocity of $\sim$ $-3$
km s-1 in the region close to the edge of the observed region. Throughout this
work, we have divided the observed region into two parts: a strong CO emission
region (hereafter SCER) with a velocity of $\sim$ $-$18 km s-1 and a weak CO
emission region (WCER) with a velocity of $\sim$ $-$3 km s-1.

Figure 1 presents an integrated intensity map of the WCER (left) and the SCER
(right). In the WCER, the 13CO emission is weak and no C18O emission is
detected. Therefore, only the integrated intensity map of 12CO, with a
velocity range of [$-$5, $-$1] km s-1, is presented here. There are three
peaks in the CO map of the WCER; i.e., Peaks 1–3 (see the labels in Figure
1(a)). In the SCER, there is strong 12CO (blue contours), 13CO (red contours),
and C18O (background map) emission, where the velocity range of the maps is
[$-$30, $-$10] km s-1. Because the various emission components are all
connected, and it is difficult to distinguish each of them clearly, we marked
the position centered at the peak intensity of C18O as Peak 4. We note that
this large emission region covers, in fact, almost the entire observed region
(including the regions covered by all seven FAST beams; see below).

(a) WCER (b) SCER

Figure 1: Integrated intensity map of the WCER (left) and SCER (right). Left:
the background map is the integrated intensity map of 12CO in the range of
[$-5$, $-1$] km s-1. Peaks 1–3 label the peaks in the WCER (i.e.,
corresponding to the intensity peaks in the 12CO map), and the red circles and
numbers present the positions and indices of the central seven beams of FAST.
Right: the background is an integrated intensity map of C18O in the range of
[$-30$, $-10$] km s-1, where the blue and red contours show the integrated
intensity map of 12CO and 13CO in the same velocity range as the C18O map,
respectively. Their levels are 0.3, 0.4, 0.5, 0.6, 0.7, and 0.8 times the
corresponding peak intensity. The color bars in the left and right panels are
in units of $\mathrm{(K\;km\;s^{-1})^{1/2}}$.

In the WCER, because the line emission of 13CO is weak, and lines from other
molecules (i.e., C18O, HCO+ and CS) are not detected, we only conducted
Gaussian fits for the 12CO spectra. The values of $v_{\mathrm{{c,m}}}$ and
$\sigma_{\mathrm{m}}$ are given based on those Gaussian fits. The peak
brightness temperatures of the 12CO spectra,
$T^{\mathrm{{}^{12}CO}}_{\mathrm{peak}}$, were used to calculate the value of
$T_{\mathrm{ex}}$ (e.g., Garden et al., 1991) as

$T_{\mathrm{ex}}=\frac{5.53}{\displaystyle\ln\left(1+\frac{5.53}{T^{\mathrm{{}^{12}CO}}_{\mathrm{peak}}+0.82}\right)}.$
(1)

The results of $v_{\mathrm{{c,m}}}$, $\sigma_{\mathrm{m}}$, and
$T_{\mathrm{ex}}$, with 1$\sigma$ errors, are $-$3.92 $\pm$ 0.01, 0.35 $\pm$
0.01 km s-1, and 7.18 $\pm$ 0.18 K for Peak 1, and $-$2.20 $\pm$ 0.01, 0.44
$\pm$ 0.01 km s-1, and 4.76 $\pm$ 0.16 K for Peak 2, and $-$2.54 $\pm$ 0.04,
0.42 $\pm$ 0.04 km s-1, and 5.08 $\pm$ 0.31 K for Peak 3, respectively. In the
SCER, the value of $T_{\mathrm{ex}}$ was also calculated based on
$T^{\mathrm{{}^{12}CO}}_{\mathrm{peak}}$. The values of $v_{\mathrm{{c,m}}}$
and $\sigma_{\mathrm{m}}$ are based on the Gaussian fit of the 13CO spectrum.
The values of $v_{\mathrm{{c,m}}}$, $\sigma_{\mathrm{m}}$, and
$T_{\mathrm{ex}}$ are $-$18.53 $\pm$ 0.01, 1.13 $\pm$ 0.01 km s-1 and 20.30
$\pm$ 0.08 K, respectively. These values were used to extract the HINSA
features.

###  3.2 HINSA Features

To extract as many HINSA features as possible, we assumed that, for the
spectrum of each FAST beam, HINSA features corresponding to each of the four
components of the molecular clouds (i.e., Peaks 1–4) exist. The specific
description of the method used to extract the HINSA feature is presented in
Section A.

We extracted fourteen HINSA features in total from the seven beams (see Figure
2). Their velocity, $V_{\mathrm{LSR}}$, peak intensity, $T_{\mathrm{ab}}$,
optical depth, $\tau_{0}$, and velocity dispersion, $\sigma_{\mathrm{HI}}$,
with 1$\sigma$ errors are listed in Table 1.

The fitted values of $\tau_{0}$ and $\sigma_{\mathrm{HI}}$ were used to
calculate the column density of the HINSA features, $N_{\mathrm{HI}}$ (see Li
& Goldsmith, 2003) as

$N_{\mathrm{HI}}=1.95\times 10^{18}\tau_{0}T_{\mathrm{ex}}\Delta V,$ (2)

where $\Delta V=\sqrt{8\ln(2)}\sigma_{\mathrm{HI}}$. The results of
$N_{\mathrm{HI}}$ and the corresponding $T_{\mathrm{ex}}$ with 1$\sigma$
errors are listed in Table 1. The values of $N_{\mathrm{HI}}$ range from
$\sim$ 5.3 to 33.1 $\times$ 1018 cm-2 in the SCER, and from $\sim$ 0.5 to 3.4
$\times$ 1018 cm-2 in the WCER.

(a) Beam 1 (b) Beam 2 (c) Beam 3 (d) Beam 4 (e) Beam 5 (f) Beam 6 (g) Beam 7

Figure 2: HINSA features of Beams 1–7, where the green lines (marked as
$T_{\mathrm{r}}$) show the observed spectrum, the red lines (marked as
$T_{\mathrm{HI}}$) the recovered spectrum of H i, and the blue lines (marked
as $T_{\mathrm{ab}}$) the fitted spectrum of the HINSA features. Table 1:
Physical Parameters of the HINSA Features

Beams | Index | $V_{\mathrm{LSR}}$ | $\sigma_{\mathrm{HI}}$ | $\tau_{\mathrm{HI}}$ | $T_{\mathrm{ex}}$ | $T_{\mathrm{ab}}$ | $N_{\mathrm{HI}}$
---|---|---|---|---|---|---|---
|  | (km s-1) | (km s-1) |  | (K) | (K) | (1018 cm-2)
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8)
Beam 1 | 1 |  $-$18.42 $\pm$ 0.02 | 1.14 $\pm$ 0.02 | 0.31 $\pm$ 0.01 | 20.30 $\pm$ 0.08 | 18.97 $\pm$ 0.008 | 33.1 $\pm$ 1.1
| 2 |  $-$3.47 $\pm$ 0.02 | 0.85 $\pm$ 0.02 | 0.12 $\pm$ 0.01 | 7.18 $\pm$ 0.18 | 11.89 $\pm$ 0.008 | 3.4 $\pm$ 0.3
Beam 2 | 1 |  $-$17.65 $\pm$ 0.02 | 0.78 $\pm$ 0.02 | 0.07 $\pm$ 0.01 | 20.30 $\pm$ 0.08 | 4.33 $\pm$ 0.006 | 5.3 $\pm$ 0.7
| 2 |  $-$3.83 $\pm$ 0.02 | 0.53 $\pm$ 0.02 | 0.03 $\pm$ 0.01 | 7.18 $\pm$ 0.18 | 2.57 $\pm$ 0.006 | 0.5 $\pm$ 0.2
Beam 3 | 1 |  $-$16.65 $\pm$ 0.02 | 1.46 $\pm$ 0.02 | 0.16 $\pm$ 0.01 | 20.30 $\pm$ 0.08 | 11.25 $\pm$ 0.007 | 21.5 $\pm$ 1.4
| 2 |  $-$2.35 $\pm$ 0.09 | 0.80 $\pm$ 0.09 | 0.04 $\pm$ 0.01 | 5.08 $\pm$ 0.31 | 4.36 $\pm$ 0.007 | 0.8 $\pm$ 0.2
Beam 4 | 1 |  $-$17.63 $\pm$ 0.02 | 0.96 $\pm$ 0.02 | 0.28 $\pm$ 0.01 | 20.30 $\pm$ 0.08 | 17.47 $\pm$ 0.011 | 25.4 $\pm$ 0.9
| 2 |  $-$2.57 $\pm$ 0.09 | 0.78 $\pm$ 0.09 | 0.11 $\pm$ 0.01 | 5.08 $\pm$ 0.31 | 10.88 $\pm$ 0.011 | 2.0 $\pm$ 0.2
Beam 5 | 1 |  $-$18.19 $\pm$ 0.02 | 1.15 $\pm$ 0.02 | 0.21 $\pm$ 0.01 | 20.30 $\pm$ 0.08 | 13.91 $\pm$ 0.006 | 22.8 $\pm$ 1.1
| 2 |  $-$3.23 $\pm$ 0.09 | 0.77 $\pm$ 0.09 | 0.13 $\pm$ 0.01 | 5.08 $\pm$ 0.31 | 13.01 $\pm$ 0.006 | 2.3 $\pm$ 0.3
Beam 6 | 1 |  $-$18.38 $\pm$ 0.02 | 1.24 $\pm$ 0.02 | 0.18 $\pm$ 0.01 | 20.30 $\pm$ 0.08 | 10.06 $\pm$ 0.006 | 20.5 $\pm$ 1.2
| 2 |  $-$2.41 $\pm$ 0.02 | 0.80 $\pm$ 0.02 | 0.10 $\pm$ 0.01 | 4.76 $\pm$ 0.16 | 9.61 $\pm$ 0.006 | 1.7 $\pm$ 0.2
Beam 7 | 1 |  $-$17.19 $\pm$ 0.02 | 0.76 $\pm$ 0.02 | 0.08 $\pm$ 0.01 | 20.30 $\pm$ 0.08 | 4.94 $\pm$ 0.008 | 5.7 $\pm$ 0.7
| 2 |  $-$2.14 $\pm$ 0.02 | 0.83 $\pm$ 0.02 | 0.14 $\pm$ 0.01 | 4.76 $\pm$ 0.16 | 12.23 $\pm$ 0.008 | 2.5 $\pm$ 0.2

Note. — Columns 3-–8 present the central velocity, velocity dispersion,
optical depth, excitation temperature, peak intensity, and column density of
HINSA features.

###  3.3 Physical Properties of the Molecular Line Emission

To make a comparison between the HINSA features and the molecular line
emission, we resampled the molecular data observed by Liu et al. (2021) to
pixel sizes of $3^{\prime}\times 3^{\prime}$, and aligned pixel center of
molecular emission with the corresponding FAST beam center.

Figure 3 presents the molecular line emission after resampling, corresponding
to each beam as well as the corresponding HINSA features, i.e.,
$T_{\mathrm{ab}}$. 12CO emission with a velocity of $\sim$ $-$18 km s-1 (i.e.,
the SCER) is detected in all seven beams, indicating that the detected 12CO
emission in the SCER covers almost the entire observed region of the molecular
lines. The line centers and velocity dispersions of the molecular lines are
from the Gaussian fits to the 12CO, 13CO, and/or C18O spectra. We also
calculated the column density of all observed molecular species, i.e., 12CO,
13CO, C18O, HCO+, and CS when these were detected (see Section B). The
velocity center, velocity dispersion, and column density as well as the
corresponding $T_{\mathrm{ex}}$, are listed in Table A1.

(a) Beam 1 (b) Beam 2 (c) Beam 3 (d) Beam 4 (e) Beam 5 (f) Beam 6 (g) Beam 7

Figure 3: HINSA features in Beams 1–7 and the molecular line spectra of the
corresponding regions. The red, green, blue, black, and yellow lines show the
spectra of 12CO, 13CO, C18O, HCO+, and CS, respectively, and the purple lines
represent the HINSA features. We fitted the profiles of both 12CO and 13CO
with a Gaussian profile for each component, except for spectra that showed
strong self-absorption or were contaminated by other components. For Beam 6,
only the line emission of C18O was fitted because there were strong self-
absorption in the spectra of 12CO and 13CO.

##  4 Discussion and and Conclusions

###  4.1 Association Between the HINSA Features and the Molecular Clouds

The association between the HINSA features and the molecular line emissions
can be seen in Figure 3. In the SCER (shown as the first entry for each beam
in Table 1), the HINSA features are probably associated with the molecular
line emission. In the WCER (shown as the second entry in Table 1), no
molecular line emission associated with HINSA features is detected toward
Beams 1, 3, 4, and 6. The only detected 12CO emission, corresponding to Beams
2, 5, and 7, is also weak. The HINSA features in the WCER are probably
associated with CO-dark or CO-faint gas (e.g., Wolfire et al., 2010; Tang et
al., 2016; Liu et al., 2022). This is an issue worthy of further observations.

In addition, a further quantitative comparison was conducted regarding the
relationship of the center velocity and the non-thermal velocity dispersion
between the HINSA features and the molecular gases. In the SCER, the
differences between the center velocity of the HINSA features (see Table 1)
and those of the 12CO, 13CO, and/or C18O lines are within the velocity
resolution of the CO spectra (i.e., $\sim$ 0.2 km s-1). This indicates that
the HINSA features are associated with the molecular clouds traced by the CO.
The non-thermal velocity dispersion of the HINSA features,
$\sigma_{\mathrm{nt,HI}}=\sqrt{\sigma_{\mathrm{HI}}^{2}-kT_{\mathrm{ex}}/m_{\mathrm{H}}}$,
is smaller than that of the 12CO lines,
$\sigma_{\mathrm{nt,^{12}CO}}=\sqrt{\sigma_{\mathrm{{}^{12}CO}}^{2}-kT_{\mathrm{ex}}/m_{\mathrm{{}^{12}CO}}}$,
and is similar to that of the 13CO lines,
$\sigma_{\mathrm{nt,^{13}CO}}=\sqrt{\sigma_{\mathrm{{}^{13}CO}}^{2}-kT_{\mathrm{ex}}/m_{\mathrm{{}^{13}CO}}}$,
except for that in Beam 3, where the value of $\sigma_{\mathrm{nt,HI}}$ is
larger than that of $\sigma_{\mathrm{nt,^{12}CO}}$ by $\sim$ 20% (see
$\sigma_{\mathrm{HI}}$ in Table 1, and $\sigma_{\mathrm{{}^{12}CO}}$ and
$\sigma_{\mathrm{{}^{13}CO}}$ in Table A1). This result is in agreement with
that of Li & Goldsmith (2003), indicating that the HINSA in the SCER is
probably mixed with the gas in cold and well-shielded regions of the molecular
clouds.

In the WCER, the association between the HINSA features and the CO cloud
appears to be poor. First, the difference between the center velocity of the
HINSA features (see Table 1) and that of 12CO is larger than in the case of
the SCER, and the maximum difference reaches $\sim$ 0.8 km s-1 (i.e., in Beam
7). Second, the value of $\sigma_{\mathrm{nt,HI}}$ is larger than that of
$\sigma_{\mathrm{nt,^{12}CO}}$. No 13CO emission is detected toward the
regions corresponding to all the seven beams after resampling. In fact, we
also do not detect any 13CO emission in Beam 3 in the SCER, and the value of
$\sigma_{\mathrm{nt,HI}}$ in this region is larger than that of
$\sigma_{\mathrm{nt,^{12}CO}}$. Therefore, it is probably more accurate to
classify Beam 3 in the SCER as a WCER component. The discussion below is based
on this reclassification. We also find that the value of
$\sigma_{\mathrm{HI}}$ in the WCER (including Beam 3 in the SCER) is
comparable to some of the values found by Tang et al. (2016) for CO-dark
clouds. These facts indicate that the HINSA toward the WCER is probably not
mixed with CO clouds, but probably associated with more diffuse CO-dark or CO-
faint clouds. Further observations of CI, CII, and/or OH may be helpful in
clarifying the relationship between the HINSA features and the CO-dark or CO-
faint gas (e.g., Li & Goldsmith, 2003; Wolfire et al., 2010; Tang et al.,
2016) toward the WCER.

###  4.2 Abundance of H i the HINSA Features

The fractional abundance of HINSA, $X_{\mathrm{HI}}$, based on the column
density is valuable. In this case
$X_{\mathrm{HI}}=N_{\mathrm{HI}}/N_{\mathrm{H}}$, where $N_{\mathrm{HI}}$ is
the column density of the HINSA, and
$N_{\mathrm{H}}=N_{\mathrm{HI}}+2N_{\mathrm{H_{2}}}$. $N_{\mathrm{H_{2}}}$ is
the column density of H2 in the molecular cloud, which can be traced by 12CO,
13CO, C18O, HCO+, and CS, i.e., $N^{\mathrm{{}^{12}CO}}_{\mathrm{H_{2}}}$,
$N^{\mathrm{{}^{13}CO}}_{\mathrm{H_{2}}}$,
$N^{\mathrm{C^{18}O}}_{\mathrm{H_{2}}}$,
$N^{\mathrm{HCO^{+}}}_{\mathrm{H_{2}}}$, and
$N^{\mathrm{CS}}_{\mathrm{H_{2}}}$ (see Table A1).

We compared the values of $N^{\mathrm{{}^{12}CO}}_{\mathrm{H_{2}}}$ and
$N^{\mathrm{{}^{13}CO}}_{\mathrm{H_{2}}}$ within all regions where we
simultaneously detected 12CO and 13CO line emission (i.e., Beams 1, 2, and 4–7
of the SCER). The results show that the value of
$N^{\mathrm{{}^{13}CO}}_{\mathrm{H_{2}}}$ is larger than that of
$N^{\mathrm{{}^{12}CO}}_{\mathrm{H_{2}}}$ by a factor of $\lesssim$ 3,
indicating that $N^{\mathrm{{}^{13}CO}}_{\mathrm{H_{2}}}$ is a better
indicator of the column density of the molecular clouds. The values of
$N^{\mathrm{{}^{12}CO}}_{\mathrm{H_{2}}}$ are probably underestimated due to
optical depth effects on the 12CO spectra. In addition, the species tracing
higher densities, i.e., C18O, HCO+, and CS, are detected only in Beams 1 and
6. In Beam 1, the ratios of the column density traced by 13CO,
$N^{\mathrm{{}^{13}CO}}_{\mathrm{H_{2}}}$, to those traced by the three dense
gas tracers, $N^{\mathrm{C^{18}O}}_{\mathrm{H_{2}}}$,
$N^{\mathrm{HCO^{+}}}_{\mathrm{H_{2}}}$, and
$N^{\mathrm{CS}}_{\mathrm{H_{2}}}$, are $\sim$ 0.92, $\sim$ 1.48, and $\sim$
0.78, respectively, indicating that $N^{\mathrm{{}^{13}CO}}_{\mathrm{H_{2}}}$
is a good measure of the column density of the molecular cloud. In Beam 6, the
dilution of the dense gas tracers is relatively stronger (see Figure 1(b)).
Nevertheless,
$N^{\mathrm{{}^{13}CO}}_{\mathrm{H_{2}}}/N^{\mathrm{C^{18}O}}_{\mathrm{H_{2}}}\sim
1.2$ also indicates that 13CO is a good tracer of the molecular cloud.
Therefore, we used $N^{\mathrm{{}^{13}CO}}_{\mathrm{H_{2}}}$ to represent the
column density of molecular clouds in the SCER.

(a) $N^{\mathrm{{}^{12}CO}}_{\mathrm{H_{2}}}$ .VS.
$N^{\mathrm{{}^{13}CO}}_{\mathrm{H_{2}}}$ (b)
$N^{\mathrm{{}^{13}CO}}_{\mathrm{H_{2}}}$ (or
$N^{\mathrm{{}^{12}CO}}_{\mathrm{H_{2}}}$) .VS. $N_{\mathrm{HI}}$

Figure 4: Left: comparison of $N^{\mathrm{{}^{12}CO}}_{\mathrm{H_{2}}}$ and
$N^{\mathrm{{}^{13}CO}}_{\mathrm{H_{2}}}$ in the SCER, where 13CO is not
detected in Beam 3. Right: comparison of the column density of HINSA,
$N_{\mathrm{HI}}$, and of the total gas, $N_{\mathrm{H}}$. The red points
represent the results corresponding to the SCER (except for Beam 3), where the
column density of the molecular clouds is traced by 13CO. The blue points
present the results corresponding to the WCER, where the column density of the
molecular clouds is traced by 12CO. The green points indicate the result
corresponding to Beam 3 in the SCER (which may be more suitably classified as
a WCER component), where the column density of the molecular clouds is traced
by 12CO.

The abundance of HINSA, $X_{\mathrm{HI}}$, in most cases of the SCER (i.e.,
Beams 1 and 4–7), ranges from $\sim$ $1.1\times 10^{-3}$ to $\sim$ $4.8\times
10^{-3}$ (see the red points in Figure 4(b)). The one exception is the value
of $X_{\mathrm{HI}}$ is $\sim$ $2.6\times 10^{-2}$ in Beam 2, where the value
of $N^{\mathrm{{}^{13}CO}}_{\mathrm{H_{2}}}$ is one order of magnitude smaller
than those in the other five beams (i.e., Beams 1 and 4–7), and the relative
error of $N^{\mathrm{{}^{13}CO}}_{\mathrm{H_{2}}}$ is large (see Figure 4(b)
and Table A1). The abundances of HINSA in Beams 1 and 4–7 are similar to those
toward optically selected dark clouds and molecular cores (i.e., with an
average value of $\sim$ 1.5 $\times$ 10-3, Li & Goldsmith, 2003; Krčo &
Goldsmith, 2010). The HINSA abundances in the SCER (including Beam 2) fall
into the HINSA abundance range toward PGCC sources (i.e., from 5.1 $\times$
10-4 to 1.3 $\times$ 10-2) of Tang et al. (2020), but larger than the HINSA
abundances derived from other PGCC sources (i.e., $\sim$ 3 $\times$ 10-4
varied by a factor of $\sim$ 3) by Liu et al. (2022).

The HINSA abundance toward Beam 1 is consistent with the abundance of H i in
the H i wind (i.e., the ratio of the column density of the H i wind to that of
the total outflowing gas) in the same beam. This fact indicates that HINSA is
probably a physical bridge that connects H i winds and molecular outflows.

We also find that the value of $X_{\mathrm{HI}}$ in Beam 1 is smaller than the
values obtained from the surrounding beams (i.e., Beams 2–7; see Figure 1(a)).
This structure seems to be a ring, but unlike that observed by Zuo et al.
(2018), the abundance of HINSA, $X_{\mathrm{HI}}$, diminishes toward the
center of the cloud. One of the reason for this could be that the central
high-mass zero-age main sequence B3 star (see Torrelles et al., 1992a) ionizes
the H i detected in Beam 1 in the SCER, which is located in an H ii region
(see Dewangan, 2019). This then decreases the abundance of H i. The ionization
can be studied using radio recombination lines (e.g., Zhang et al., 2014,
2021).

The mean densities of H i, $\bar{n}_{\mathrm{HI}}$, and total proton,
$\bar{n}_{\mathrm{H}}$, in Beam 1 are $\sim$ 7.2 and $\sim$ 6.1 $\times$ 103
cm-3, respectively, assuming that the cloud diameter is $\sim$ 1.5 pc,
corresponding to the diameter being the HPBW of FAST. We assume that the
density profile has a constant density core surrounded by an envelope of the
form

$n(r)=\frac{2n_{0}}{1+(r/r_{0})^{2}},$ (3)

where $n_{0}$ and $r_{0}$ (set to 1017 cm) are the central density and the
radius of the constant density core. In this case, $n_{0}\sim
9.4\,\bar{n}_{\mathrm{H}}$ would be larger than that of model 1 in Goldsmith
et al. (2007), indicating that the density is high enough to initiate the
conversion of H i $\rightarrow$ H2. Following the model developed by Goldsmith
et al. (2007), the time that has elapsed since the material was UV irradiated
is $\sim$ 2 $\times$ 105 yr (see Section C), during which H i could convert to
H2. This time is consistent with the timescale of the star formation in this
region (i.e., $\sim$ 2 $\times$ 105 yr) derived from near-infrared H2 knots by
Chen et al. (2003). The model of Goldsmith et al. (2007) also implies another
possible reason why $X_{\mathrm{HI}}$ diminishes toward the center of the
cloud, i.e., the conversion to molecular form occurs more rapidly in the
higher-density inner region.

In the WCER, a comparison of $N_{\mathrm{HI}}$ and $N_{\mathrm{H}}$ (only for
the regions where CO emission is detected) is plotted as blue and green points
in Figure 4(b). The abundance of the HINSA features, $X_{\mathrm{HI}}$, is
$\sim$ 0.04, $\sim$ 0.14, and $\sim$ 0.17 in Beams 2, 4, and 7, respectively.
The value of $X_{\mathrm{HI}}$ is $\sim$ 0.14 in Beam 3 in the SCER where 13CO
emission is not detected. The high abundances are most likely caused by an
underestimation of the column density of the molecular clouds, because the
associated molecular clouds (including the region in Beam 3 in the SCER) are
probably CO-dark or CO-faint gas (e.g., Bolatto et al., 2013; Tang et al.,
2016).

The comparison of the HINSA features and the molecular line emission indicates
that the HINSA is probably mixed with the gas in cold and well-shielded CO
clouds in which 13CO emission is detected. There are HINSA features that have
no detected 13CO and/or 12CO emission counterparts, meaning they are probably
associated with CO-dark or CO-faint gas. In addition, HINSA is a promising
physical bridge to connect the H i winds with molecular outflows. There is a
possible ring where the abundance of HINSA diminishes toward the center of the
cloud.

##  Appendix A Extracting HINSA

The method used here to extract HINSA is based on the methodology developed by
Liu et al. (2022). The first step of this method is to obtain the initial
values of the optical depth, $\tau_{i}$, velocity dispersion, $\sigma_{i}$,
and central velocity, $v_{i}$. $\tau_{i}$ is initially set to 0.1. Obtaining
the initial values of $v_{i}$ and $\sigma_{i}$ requires knowledge derived from
the corresponding high-sensitivity CO spectrum, including the central
velocity, $v_{\mathrm{{c,m}}}$, velocity dispersion, $\sigma_{\mathrm{m}}$,
and excitation temperature, $T_{\mathrm{ex}}$.To extract HINSA as much as
possible, we assume that every component of CO emission (i.e., Peaks 1–4)
corresponds to a HINSA feature. The initial value of $v_{i}$ equals
$v_{\mathrm{{c,m}}}$, and the initial value of $\sigma_{i}$ reads as

$\sigma_{i}=\sqrt{\sigma_{\mathrm{m}}^{2}+\frac{kT_{\mathrm{ex}}}{m_{\mathrm{H}}}-\frac{kT_{\mathrm{ex}}}{m_{\mathrm{CO}}}},$
(A1)

where $k$ is the Boltzmann constant, $m_{\mathrm{H}}$ the mass of atomic
hydrogen, and $m_{\mathrm{CO}}$ the mass of the 12CO molecule for Peaks 1–3
and of the 13CO molecule for Peak 4.

The optical depth of a HINSA feature can be written (Krčo et al., 2008) as

$\tau(v)=\sum_{m-1}^{i=0}\displaystyle\tau_{i}e^{\displaystyle-\frac{(v-v_{i})^{2}}{2\sigma_{i}^{2}}},$
(A2)

where $\tau_{i}$, $\sigma_{i}$, and $v_{i}$ are the initial values of the
optical depth, velocity dispersion, and central velocity of the $i$-th
component of the HINSA features, respectively.

The H i emission free from self-absorption, $T_{\mathrm{HI}}$, and the self-
absorption of H i, $T_{\mathrm{ab}}$, satisfies (Li & Goldsmith, 2003)

$T_{\mathrm{ab}}=T_{\mathrm{HI}}-T_{\mathrm{r}}=\left[pT_{\mathrm{HI}}+\left(T_{\mathrm{c}}-T_{\mathrm{ex}}\right)\left(1-\tau_{\mathrm{f}}\right)\right](1-e^{-\tau}),$
(A3)

where $p$ is a parameter representing the proportion of the background H i
optical depth, $T_{\mathrm{r}}$ is the observed H i spectrum, $T_{\mathrm{c}}$
is the H i continuum brightness temperature, $T_{\mathrm{ex}}$ is the
excitation temperature (derived from the spectrum of 12CO), and
$\tau_{\mathrm{f}}$ is the foreground H i optical depth.

Following the treatment of Liu et al. (2022), $p$ was set to unity, and then
$\tau_{f}=0$. $T_{\mathrm{ex}}$ and $T_{c}$ were simultaneously ignored.
$T_{\mathrm{HI}}$ then can be expressed as

$T_{\mathrm{HI}}=T_{\mathrm{r}}e^{\tau},$ (A4)

which can be rewritten as

$T_{\mathrm{HI}}=T_{\mathrm{r}}+(1-e^{\tau})T^{f}_{\mathrm{HI,smooth}},$ (A5)

where $T^{f}_{\mathrm{HI,smooth}}$ is the polynomial fit of $T_{\mathrm{HI}}$.

$\tau_{i}$, $\sigma_{i}$, and $v_{i}$ (Krčo et al., 2008) could be obtained by
minimizing $\mathscr{R}$

$\mathscr{R}\cong\sum^{N}_{j=1}(T^{\prime\prime}_{\mathrm{HI}})_{j}^{2}\Delta_{\mathrm{ch}},$
(A6)

where

$T^{\prime\prime}_{\mathrm{HI}}=\frac{T_{\mathrm{HI},j+1}+T_{\mathrm{HI},j-1}-2T_{\mathrm{HI},j}}{\Delta_{\mathrm{ch}}^{2}},$
(A7)

and where $\Delta_{\mathrm{ch}}$ is the channel width of the H i spectrum and
$j$ the index of the channel. Similar to Liu et al. (2022), the L-BFGS-B
algorithm from Scipy 222https://www.scipy.org/(Virtanen et al., 2020) was used
to obtain $\tau_{i}$, $\sigma_{i}$, and $v_{i}$. For more details of the
methodology of extracting HINSA, see Liu et al. (2022). The final values of
the optical depth, velocity dispersion, and central velocity of the HINSA
features are denoted $\tau_{0}$, $\sigma_{\mathrm{HI}}$, and
$V_{\mathrm{LSR}}$, respectively.

##  Appendix B Column Densities of the Molecules

For the tracer of 13CO, under local thermal equilibrium, the column density of
H2, $N^{\mathrm{{}^{13}CO}}_{\mathrm{H_{2}}}$, is (Snell et al., 1984)

$N^{\mathrm{{}^{13}CO}}_{\mathrm{H_{2}}}=2.3\times
10^{19}\frac{T_{\mathrm{ex}}}{e^{-5.29/T_{\mathrm{ex}}}}\int\frac{e^{\tau_{\mathrm{{}^{13}CO}}}}{1-e^{-\tau_{\mathrm{{}^{13}CO}}}}T_{\mathrm{mb}}^{{}^{13}\mathrm{CO}}dv,$
(B1)

where $X_{\mathrm{{}^{13}CO}}$ = [13CO]/[H2] = $2\times 10^{-6}$ is used
(Dickman, 1978). In the third term, the optical depth
($\tau_{\mathrm{{}^{13}CO}}$) correction,
$e^{\tau_{\mathrm{{}^{13}CO}}}/(1-e^{-\tau_{\mathrm{{}^{13}CO}}})$, is
considered only when the values of $T_{\mathrm{mb}}^{{}^{13}\mathrm{CO}}$ are
larger than three times the corresponding rms noise. In this term, the optical
depth of 13CO is

$\tau_{\mathrm{{}^{13}CO}}=-\ln\left[1-\frac{T_{\mathrm{mb}}^{{}^{13}\mathrm{CO}}/5.29}{1/(e^{5.29}/T_{\mathrm{ex}}-1)-0.16}\right].$
(B2)

$\tau_{\mathrm{{}^{13}CO}}$ can be used to calculate the optical depth of
12CO, $\tau_{\mathrm{{}^{12}CO}}$, by multiplying the abundance ratio,
[12CO]/[13CO] $\sim$ 50\. The column density of H2 traced by 12CO is therefore

$N^{\mathrm{{}^{12}CO}}_{\mathrm{H_{2}}}=4.2\times
10^{17}\frac{T_{\mathrm{ex}}}{e^{-5.53/T_{\mathrm{ex}}}}\int\frac{e^{\tau_{\mathrm{{}^{12}CO}}}}{1-e^{-\tau_{\mathrm{{}^{12}CO}}}}T_{\mathrm{mb}}^{{}^{12}\mathrm{CO}}dv,$
(B3)

where $X_{\mathrm{{}^{12}CO}}$ = [12CO]/[H2] = $10^{-4}$ is used in this work
(Snell et al., 1984). If 13CO emission is not detected, the optical depth
correction in the third term is set to unity.

Similarly, the column density of H2 traced by C18O is

$N^{\mathrm{C^{18}O}}_{\mathrm{H_{2}}}=2.4\times
10^{20}\frac{T_{\mathrm{ex}}}{e^{-5.27/T_{\mathrm{ex}}}}\int\frac{e^{\tau_{\mathrm{C^{18}O}}}}{1-e^{-\tau_{\mathrm{C^{18}O}}}}T_{\mathrm{mb}}^{\mathrm{C^{18}O}}dv,$
(B4)

with

$\tau_{\mathrm{C^{18}O}}=-\ln\left[1-\frac{T_{\mathrm{mb}}^{\mathrm{C^{18}O}}/5.27}{1/(e^{5.27}/T_{\mathrm{ex}}-1)-0.17}\right],$
(B5)

where $X_{\mathrm{C^{18}O}}$ = [C18O]/[H2] = $2\times 10^{-7}$ is used (Garden
et al., 1991). For HCO+ (Yang et al., 1991) we have

$N^{\mathrm{HCO^{+}}}_{\mathrm{H_{2}}}=1.87\times
10^{19}\frac{T_{\mathrm{ex}}}{1-e^{-4.28/T_{\mathrm{ex}}}}\int\frac{e^{\tau_{\mathrm{HCO^{+}}}}}{1-e^{-\tau_{\mathrm{HCO^{+}}}}}T_{\mathrm{mb}}^{\mathrm{HCO^{+}}}dv,$
(B6)

with

$\tau_{\mathrm{HCO^{+}}}=-\ln\left[1-\frac{T_{\mathrm{mb}}^{\mathrm{HCO^{+}}}/4.28}{1/(e^{4.28}/T_{\mathrm{ex}}-1)-0.26}\right],$
(B7)

where $X_{\mathrm{HCO^{+}}}$ = [HCO+]/[H2] = $10^{-8}$ is used (Turner et al.,
1997). Finally, for CS we have (Liu et al., 2021)

$N^{\mathrm{CS}}_{\mathrm{H_{2}}}=1.89\times
10^{20}\frac{T_{\mathrm{ex}}}{e^{-4.70/T_{\mathrm{ex}}}}\int\frac{e^{\tau_{\mathrm{CS}}}}{1-e^{-\tau_{\mathrm{CS}}}}T_{\mathrm{mb}}^{\mathrm{CS}}dv,$
(B8)

with

$\tau_{\mathrm{CS}}=-\ln\left[1-\frac{T_{\mathrm{mb}}^{\mathrm{CS}}/4.70}{1/(e^{4.70}/T_{\mathrm{ex}}-1)-0.22}\right],$
(B9)

where $X_{\mathrm{CS}}$ = [CS]/[H2] = 10-9 is used (Tatematsu et al., 1998).
$N^{\mathrm{{}^{12}CO}}_{\mathrm{H_{2}}}$,
$N^{\mathrm{{}^{13}CO}}_{\mathrm{H_{2}}}$,
$N^{\mathrm{C^{18}O}}_{\mathrm{H_{2}}}$,
$N^{\mathrm{HCO^{+}}}_{\mathrm{H_{2}}}$, and
$N^{\mathrm{CS}}_{\mathrm{H_{2}}}$ are listed in Table A1.

Table A1: Physical Parameters Traced by the Molecules

Beams | Index | $V_{\mathrm{{}^{12}CO}}$ | $V_{\mathrm{{}^{13}CO}}$ | $\sigma_{\mathrm{{}^{12}CO}}$ | $\sigma_{\mathrm{{}^{13}CO}}$ | $T_{\mathrm{ex}}$ | $N^{\mathrm{{}^{12}CO}}_{\mathrm{H_{2}}}$ | $N^{\mathrm{{}^{13}CO}}_{\mathrm{H_{2}}}$ | $N^{\mathrm{C^{18}O}}_{\mathrm{H_{2}}}$ | $N^{\mathrm{HCO^{+}}}_{\mathrm{H_{2}}}$ | $N^{\mathrm{CS}}_{\mathrm{H_{2}}}$
---|---|---|---|---|---|---|---|---|---|---|---
|  | (km s-1) | (km s-1) | (km s-1) | (km s-1) | (K) | (1020 cm-2) | (1021 cm-2) | (1021 cm-2) | (1021 cm-2) | (1021 cm-2)
Beam 1 | 1 | … |  $-$18.33 $\pm$ 0.04 | … | 0.93 $\pm$ 0.04 | 20.30 $\pm$ 0.08 | 74.72 $\pm$ 0.73 | 14.62 $\pm$ 0.08 | 15.84 $\pm$ 0.53 | 10.29 $\pm$ 0.13 | 21.26 $\pm$ 0.16
Beam 2 | 1 |  $-$17.93 $\pm$ 0.03 | … | 1.01 $\pm$ 0.03 | … | 20.30 $\pm$ 0.8 | 0.81 $\pm$ 0.05 | 0.10 $\pm$ 0.09 | … | … | …
| 2 |  $-$3.69 $\pm$ 0.09 | … | 0.45 $\pm$ 0.09 | … | 7.18 $\pm$ 0.18 | 0.05 $\pm$ 0.01 | … | … | … | …
Beam 3 | 1 |  $-$16.71 $\pm$ 0.03 | … | 1.14 $\pm$ 0.03 | … | 20.30 $\pm$ 0.08 | 0.68 $\pm$ 0.04 | … | … | … | …
Beam 4 | 1 |  $-$17.55 $\pm$ 0.01 |  $-$17.69 $\pm$ 0.01 | 0.95 $\pm$ 0.01 | 0.74 $\pm$ 0.01 | 20.30 $\pm$ 0.08 | 11.32 $\pm$ 0.21 | 2.61 $\pm$ 0.10 | … | … | …
| 2 |  $-$2.66 $\pm$ 0.02 | … | 0.36 $\pm$ 0.02 | … | 5.08 $\pm$ 0.31 | 0.06 $\pm$ 0.01 | … | … | … | …
Beam 5 | 1 | … |  $-$18.28 $\pm$ 0.02 | … | 0.73 $\pm$ 0.04 | 20.30 $\pm$ 0.08 | 8.32 $\pm$ 0.14 | 2.36 $\pm$ 0.08 | … | … | …
Beam 6 | 1 |  $-$18.35 $\pm$ 0.04$\dagger$$\dagger$footnotemark:  | … | 0.58 $\pm$ 0.04$\dagger$$\dagger$footnotemark:  | … | 20.30 $\pm$ 0.08 | 19.89 $\pm$ 0.27 | 5.23 $\pm$ 0.09 | 4.18 $\pm$ 0.90 | 0.48 $\pm$ 0.14 | 0.65 $\pm$ 0.23
Beam 7 | 1 |  $-$17.22 $\pm$ 0.01 |  $-$17.52 $\pm$ 0.02 | 1.11 $\pm$ 0.01 | 0.74 $\pm$ 0.02 | 20.30 $\pm$0.08 | 3.63 $\pm$ 0.10 | 1.23 $\pm$ 0.09 | … | … | …
| 2 |  $-$1.34 $\pm$ 0.09 | … | 0.72 $\pm$ 0.09 | … | 4.76 $\pm$ 0.16 | 0.06 $\pm$ 0.01 | … | … | … | …

††footnotetext: Determined from the Gaussian fit of the C18O spectrum.

##  Appendix C Transition between Atomic and Molecular Phase

We follow the model developed by Goldsmith et al. (2007) to study the
transition between H i and H2. We only focus on the central constant core, and
assume that the central density is $n_{0}\sim 5.7\times 10^{4}$ cm-3, line
width is 1.14 km s-1, radius is $10^{17}$ cm, and $T=10$ K. The other
parameters follow these in Goldsmith et al. (2007). The equations (3)–(5) in
Goldsmith et al. (2007) would produce the evolution of the density of H i,
$n_{\mathrm{HI}}$, and H2, $n_{\mathrm{H_{2}}}$. Figure 5 shows
$n_{\mathrm{HI}}$ and $X_{\mathrm{HI}}$ as a function of time, where
$X_{\mathrm{HI}}$ set to
$n_{\mathrm{HI}}/(n_{\mathrm{HI}}+2n_{\mathrm{H_{2}}})$ assuming that the
sizes of cold H i and CO cloud are similar (see e.g., Goldsmith et al., 2007).
The time that has elapsed since the material was UV irradiated would be $\sim$
2 $\times$ 105 yr assuming $X_{\mathrm{HI}}\sim 1.1\times 10^{-3}$.

Figure 5: $n_{\mathrm{HI}}$ and $X_{\mathrm{HI}}$ as a function of time. The H
i density in the cloud core drops to its steady state value of $\sim$ 2 cm-3
with a time of $\sim$ $4\times 10^{5}$ yr.

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 all the staff members of the
Delingha 13.7 m radio telescopes for their support making the various
molecular line observations. 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. BK20210999), the Entrepreneurship and Innovation Program of Jiangsu
Province, NSFC grants Nos. 11933011 and 11873019, and the Key Laboratory for
Radio Astronomy, Chinese Academy of Sciences. FAST, PMO 13.7m

## References

  * Bolatto et al. (2013) Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207
  * Chen et al. (2003) Chen, Y., Yao, Y., Yang, J., Zeng, Q., & Sato, S. 2003, A&A, 405, 655
  * Dewangan (2019) Dewangan, L. K. 2019, ApJ, 884, 84
  * Dickman (1978) Dickman, R. L. 1978, ApJS, 37, 407
  * Garden et al. (1991) Garden, R. P., Hayashi, M., Gatley, I., Hasegawa, T., & Kaifu, N. 1991, ApJ, 374, 540
  * Gibson (2010) Gibson, S. J. 2010, in Astronomical Society of the Pacific Conference Series, Vol. 438, The Dynamic Interstellar Medium: A Celebration of the Canadian Galactic Plane Survey, ed. R. Kothes, T. L. Landecker, & A. G. Willis, 111
  * Gibson et al. (2000) Gibson, S. J., Taylor, A. R., Higgs, L. A., & Dewdney, P. E. 2000, ApJ, 540, 851
  * Goldsmith & Li (2005) Goldsmith, P. F., & Li, D. 2005, ApJ, 622, 938
  * Goldsmith et al. (2007) Goldsmith, P. F., Li, D., & Krčo, M. 2007, ApJ, 654, 273
  * Hasegawa et al. (1983) Hasegawa, T., Sato, F., & Fukui, Y. 1983, AJ, 88, 658
  * Heeschen (1955) Heeschen, D. S. 1955, ApJ, 121, 569
  * Jiang et al. (2019) Jiang, P., Yue, Y., Gan, H., et al. 2019, SCPMA, 62, 959502
  * Jiang et al. (2020) Jiang, P., Tang, N.-Y., Hou, L.-G., et al. 2020, RAA, 20, 064
  * Jiang et al. (2013) Jiang, Z.-B., Chen, Z.-W., Wang, Y., et al. 2013, RAA, 13, 695
  * Knapp (1974) Knapp, G. R. 1974, AJ, 79, 527
  * Krčo & Goldsmith (2010) Krčo, M., & Goldsmith, P. F. 2010, ApJ, 724, 1402
  * Krčo et al. (2008) Krčo, M., Goldsmith, P. F., Brown, R. L., & Li, D. 2008, ApJ, 689, 276
  * Lee & Myers (1999) Lee, C. W., & Myers, P. C. 1999, ApJS, 123, 233
  * Li & Goldsmith (2003) Li, D., & Goldsmith, P. F. 2003, ApJ, 585, 823
  * Li et al. (2018) Li, D., Wang, P., Qian, L., et al. 2018, IMMag, 19, 112
  * Liu et al. (2021) Liu, D.-J., Xu, Y., Li, Y.-J., et al. 2021, ApJS, 253, 15
  * Liu et al. (2022) Liu, X., Wu, Y., Zhang, C., et al. 2022, A&A, 658, A140
  * Lynds (1962) Lynds, B. T. 1962, ApJS, 7, 1
  * Moffat et al. (1979) Moffat, A. F. J., Fitzgerald, M. P., & Jackson, P. D. 1979, A&AS, 38, 197
  * Nan (2006) Nan, R. 2006, SCPMA, 49, 129
  * Nan et al. (2011) Nan, R., Li, D., Jin, C., et al. 2011, IJMPD, 20, 989
  * Planck Collaboration et al. (2011) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2011, A&A, 536, A23
  * Planck Collaboration et al. (2016) —. 2016, A&A, 594, A28
  * Qian et al. (2020) Qian, L., Yao, R., Sun, J., et al. 2020, Innov, 1, 100053
  * Snell et al. (1988) Snell, R. L., Huang, Y. L., Dickman, R. L., & Claussen, M. J. 1988, ApJ, 325, 853
  * Snell et al. (1984) Snell, R. L., Scoville, N. Z., Sanders, D. B., & Erickson, N. R. 1984, ApJ, 284, 176
  * Tang et al. (2016) Tang, N., Li, D., Heiles, C., et al. 2016, A&A, 593, A42
  * Tang et al. (2020) Tang, N.-Y., Zuo, P., Li, D., et al. 2020, RAA, 20, 077
  * Tatematsu et al. (1998) Tatematsu, K., Umemoto, T., Heyer, M. H., et al. 1998, ApJS, 118, 517
  * Torrelles et al. (1992a) Torrelles, J. M., Eiroa, C., Mauersberger, R., et al. 1992a, ApJ, 384, 528
  * Torrelles et al. (1992b) Torrelles, J. M., Gomez, J. F., Anglada, G., et al. 1992b, ApJ, 392, 616
  * Turner et al. (1997) Turner, B. E., Pirogov, L., & Minh, Y. C. 1997, ApJ, 483, 235
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261
  * Wang et al. (2020) Wang, Y., Bihr, S., Beuther, H., et al. 2020, A&A, 634, A139
  * Wolfire et al. (2010) Wolfire, M. G., Hollenbach, D., & McKee, C. F. 2010, ApJ, 716, 1191
  * Wu et al. (2012) Wu, Y., Liu, T., Meng, F., et al. 2012, ApJ, 756, 76
  * Yang et al. (1991) Yang, J., Umemoto, T., Iwata, T., & Fukui, Y. 1991, ApJ, 373, 137
  * Zhang et al. (2014) Zhang, C.-P., Wang, J.-J., Xu, J.-L., Wyrowski, F., & Menten, K. M. 2014, ApJ, 784, 107
  * Zhang et al. (2021) Zhang, C.-P., Xu, J.-L., Li, G.-X., et al. 2021, RAA, 21, 209
  * Zhang et al. (2005) Zhang, Q., Hunter, T. R., Brand, J., et al. 2005, ApJ, 625, 864
  * Zuo et al. (2018) Zuo, P., Li, D., Peek, J. E. G., et al. 2018, ApJ, 867, 13


# Model-Driven Synthesis for Programming Tutors

Niek Mulleners  Utrecht UniversityUtrecht<EMAIL_ADDRESS>,  Johan Jeuring
Utrecht UniversityUtrecht<EMAIL_ADDRESS>and  Bastiaan Heeren  Open
University of the NetherlandsHeerlen<EMAIL_ADDRESS>
(2020)

###### Abstract.

When giving automated feedback to a student working on a beginner’s exercise,
many programming tutors run into a completeness problem. On the one hand, we
want a student to experiment freely. On the other hand, we want a student to
write her program in such a way that we can provide constructive feedback. We
propose to investigate how we can overcome this problem by using program
synthesis, which we use to generate correct solutions that closely match a
student program, and give feedback based on the results.

synthesis, programming tutors

††copyright: rightsretained††journalyear: 2020††conference: HATRA 2020: Human
Aspects of Types and Reasoning Assistants; 15–20 November, 2020; Chicago,
IL††booktitle: HATRA 2020: Human Aspects of Types and Reasoning Assistants,
15–20 November, 2020, Chicago, IL††ccs: Applied computing Interactive learning
environments††ccs: Software and its engineering Search-based software
engineering

##  1\. Introduction

Programming tutors help students learning to develop a program. Since the
1960s, more than one hundred of such tools have been developed (Keuning et
al., 2018). Many of these tutors analyse the final program a student submits,
but there are also tutors that analyse a partial, unfinished student program,
and give feedback for such programs. Keuning et al. (2018) note that most of
these tools give solution error feedback and that relatively few tools give
hints to students when they are on their way to a solution. We want to
investigate how we can give hints to students that are developing a program,
at a stage where a program need not be complete yet.

Giving hints to students who are halfway towards a solution is challenging.
The literature describes several techniques to achieve this. Anderson and
Skwarecki (1986) and Gerdes et al. (2017) require a student to follow a path
towards a model solution, where there are usually multiple paths possible.
Price et al. (2016), Piech et al. (2015), and Lazar and Bratko (2014) use data
collected from previous students to generate hints. The disadvantages of these
approaches are that they either restrict what a student can do, or that they
are incomplete in the sense that a student can submit a program for which no
hint can be given. The latter need not be a problem if this happens
infrequently, but current practice shows that this happens regularly. We want
to investigate techniques for giving hints to students who can develop any
program for a problem, and ask for hints at any stage.

We will focus on students working on beginners exercises in a statically typed
functional programming language. Such an exercise is defined by

  * •

a textual description;

  * •

a set of model solutions;

  * •

a set of properties;

  * •

a library of allowed functions.

We allow students to leave parts of their program unimplemented by introducing
holes, denoted by a question mark (`?`). For an unfinished student program, we
want to answer the following questions: is the student still on the right
track? If so, how can we give hints on how to proceed? If not, how can we help
the student to get back on track?

Feldman et al. (2019) use program synthesis in a dynamically typed functional
language to tell students whether they are on the right track, by generating
the part of the program the student still needs to write. Singh et al. (2013)
and Gulwani et al. (2018) apply program synthesis to imperative languages to
give automated feedback on how to repair a program. We want to adapt these
approaches to a statically typed functional language and use the synthesis
results to give constructive feedback. We propose to investigate to what
extent we can use program synthesis to give hints to students on their way to
a solution.

##  2\. An example: sorting algorithm

Throughout this proposal, we use the following exercise as an example.

### Description

Write a function `my_sort :: [Int] -> [Int]` that sorts a list of integers.

### Model solutions

There are many possible ways to implement a sorting algorithm. For simplicity,
we only look at the following model solution, implementing `my_sort` as
insertion sort.

    my_sort = foldr insert []
      where
        insert x [] = [x]
        insert x (y:ys) | x < y     = x:y:ys
                        | otherwise = y:insert x ys

### Properties

The correctness of `my_sort` is expressed in terms of the following
properties, expressing that the result of sorting a list should be a
permutation of the input list and that the result should be in non-descending
order.

    sort_permutes xs = my_sort xs ‘permutes‘ xs
    sort_nondescending xs = nondescending (my_sort xs)

### Library

A suitable prelude containing all functions the students should reasonably be
allowed to use.

##  3\. Program Synthesis

Program synthesis is the process of automatically generating programs based on
their intended behaviour. Rather than writing programs by hand, the programmer
expresses their intention in terms of a specification to which the program
should adhere. Such a specification can range from simple input-output
examples to a formal logical specification. A program synthesizer then tries
to generate a program from this specification, typically by performing a
search, enumerating all possible programs within the program space, until one
is found that adheres to the specification.

To verify the correctness of student solutions, we extensively test them
against a large set of input-output examples, generated from the model
solutions. Therefore, as our synthesizer, we will look at Myth (Osera and
Zdancewic, 2015; Frankle et al., 2016), a state-of-the-art tool in example-
based synthesizers. During the synthesis process, Myth can only introduce
constructs that it can execute, so that the synthesized program can be checked
against the input-output examples. When synthesizing a recursive function,
Myth can “execute” a recursive call if the input to the recursive call is in
the example set, by simply returning the corresponding output. This requires,
however, that the example set is trace complete. That is, for every input-
output pair in the example set, there should also be an input-output pair for
every structurally smaller input, so that the example set contains complete
program traces. This requirement is often considered the main limitation of
Myth. In the presence of model solutions, however, we can simply execute a
model solution at the recursive call, bypassing the need for trace-complete
example sets!

Furthermore, we expect correct solutions to resemble our model solutions. We
can use this insight to guide the synthesis process by giving priority to
constructs appearing in the model solutions. Lee et al. (2018) accelerate
search-based synthesis by performing A* search (Hart et al., 1968), with a
heuristic based on learned probabilistic models. We intend to use a similar
approach to speed up the synthesis, by basing our heuristic on the provided
model solutions.

Unfortunately, Myth cannot be applied directly to our setting, because it
requires local input-output examples on holes, whereas we only have global
input-output examples. As described by Lubin et al. (2020), we can use live
bidirectional evaluation to propagate global input-output examples to local
input-output examples on holes. First, live evaluation (Omar et al., 2019) is
used to partially evaluate a program to a result r by proceeding around the
holes. Live unevaluation (Lubin et al., 2020) then checks this result r
against an input-output example, generating constraints on the holes in the
form of local input-output examples.

##  4\. Finishing programs

In this section, we give an overview of how we intend to apply the described
techniques to complete unfinished student programs and use the results to give
feedback.

###  4.1. Live bidirectional example checking

Before we can perform program synthesis, we apply live bidirectional example
checking to check the student solution against an example set, possibly
returning a set of local input-output examples for every hole.

####  4.1.1. Failure

If live bidirectional example checking fails, this means that it is impossible
to finish a student solution by filling its holes. This results in a
counterexample, which can be checked against the property set to determine
which properties are violated.

##### Finished solution

The student might try to implement `my_sort` as insertion sort, but
incorrectly use the list cons operator instead of the `insert` function.

    my_sort = foldr (:) []

Live evaluation might return the counterexample `my_sort [3,2,1] == [3,2,1]`,
violating the property `sort_nondescending`.

##### Partial solution

Alternatively, the student might give a partial solution with an incorrect
base case, such as the following.

    my_sort = foldr ? [0]

Live evaluation gives the counterexample `my_sort [] == [0]`, which violates
the property `sort_permutes`. Note that, in some cases, the resulting
counterexample might also contain a hole! For example, take the following
partial solution.

    my_sort [] = []
    my_sort (x:xs) = x : ?

Live evaluation might return the counterexample `my_sort [2,1] == 2:?`. Using
unevaluation, we can determine that there is no hole filling for which the
epxression `2:?` evaluates to the expected result `[1,2]`.

##### Recovery

Other than presenting the student the counterexample and showing which
properties her implementation violates, we might want to help her get back on
track. Each counterexample corresponds to a path in the abstract syntax tree
of the program, hinting at what part of the program is incorrect. For example,
the counterexample `my_sort [3,2,1] == [3,2,1]` is introduced by the leaf
`(:)` in the abstract syntax tree of `foldr (:) []`. We replace this leaf with
a hole (resulting in the program `foldr ? []`) and repeat the synthesis
process until we eventually find a correct solution closest to the student’s
attempt.

####  4.1.2. Success

On successful live bidirectional evaluation, we get a set of local input-
output examples on the holes. If the student solution contains no holes, this
set of examples is empty and the exercise is finished. If there are holes, the
generated examples can be used to perform program synthesis.

###  4.2. Program synthesis

Now that we have local input-output examples on the holes in our program, we
can perform program synthesis to determine whether the student is still on
track.

#### Failure

If a program can no longer be finished, we expect our synthesis to fail. This
means that the synthesis has either exhaustively searched the program space or
a timeout occurred. Both of these cases can take quite some time. One way to
skip this process is to check for conflicting examples. Take, for example, the
following attempt at implementing `my_sort` using `map`.

    my_sort = map ?0

For the input `[2,2,1]` we expect the output `[1,2,2]`. As such, for a
function `f` filling the hole we generate the following input-output examples.

    f 2 == 2
    f 2 == 1
    f 1 == 2

These constraints are conflicting, since a function cannot have different
outputs on the same input. As such, we can skip the synthesis process. For
other incorrect programs, it might not be possible to generate conflicting
constraints. Consider the following example:

    my_sort = map ? . zip [0..]

Before mapping, it first combines every value in the input list with a unique
integer. In doing so, the list passed to `map ?` contains only unique values,
so that no conflicting constraints can be generated. At this point, we will
perform program synthesis with the generated constraints, but the process will
fail since `map` cannot change the order of a list.

##### Recovery

To help the student recover, we propose an approach similar to the one
described in section 4.1.1. If the synthesis process fails to fill a hole, we
blame its parent node in the abstract syntax tree of the program and replace
that node with a hole. For example, for the incorrect program `map ? . zip
[0..]`, we replace `map ?` with a hole, resulting in the program `? . zip
[0..]`. Like before, we repeat the synthesis process until we find a correct
solution.

#### Success

If a program can still be finished, our synthesis will either timeout, or it
will return a valid hole filling. In the first case, we cannot tell a student
whether she is on the right track, but we can guarantee that her
implementation is more complex than required and advise another approach. In
the second case, we can give her feedback by using the synthesized programs to
calculate specifications for the holes and present these in a legible way. For
example, take the following partial student program.

    my_sort [] = []
    my_sort (x:xs) = f x (my_sort xs)
      where f y ys = ?

This program can easily be finished by filling the hole with `insert y ys`.
Knowing this, we can generate a specification for the hole, by generating
input-output examples using the hole filling. These input-output examples can
help the student understand the part of the program they are yet to write.
Even more interesting, we might use a technique such as described by Claessen
et al. (2010) to automatically generate properties that should hold on
`insert` and therefore should hold on the hole. Ideally, this would allow us
to inform the student that `ys` is a sorted list and `f` should return a
sorted list containing all elements from `ys` plus `y`.

Note that many existing techniques, such as program transformations (Gerdes et
al., 2017), can already give relevant feedback on the previous example by
comparing it to one of the model solutions. The next example, however, is a
lot harder to reason about automatically.

    my_sort [] = []
    my_sort (x:xs) = foldr ? ? xs

It is far from trivial to automatically prove that this partial program is
still on the right track, but we expect our synthesizer to find the correct
hole fillings (`insert` and `[x]` respectively) with relative ease.

##  5\. Conclusion

To conclude, we believe there is a lot of potential in the application of
program synthesis and live bidirectional evaluation to give meaningful
feedback to students on their way to a solution.

## References

  * (1)
  * Anderson and Skwarecki (1986) John R. Anderson and Edward Skwarecki. 1986.  The automated tutoring of introductory computer programming.  _Commun. ACM_ 29, 9 (1986), 842–849.
  * Claessen et al. (2010) Koen Claessen, Nicholas Smallbone, and John Hughes. 2010\.  QuickSpec: Guessing Formal Specifications Using Testing. In _Tests and Proofs_ , Gordon Fraser and Angelo Gargantini (Eds.). Springer Berlin Heidelberg, Berlin, Heidelberg, 6–21.
  * Feldman et al. (2019) Molly Q Feldman, Yiting Wang, William E. Byrd, François Guimbretière, and Erik Andersen. 2019\.  Towards Answering ”Am I on the Right Track?” Automatically Using Program Synthesis. In _Proceedings of the 2019 ACM SIGPLAN Symposium on SPLASH-E_ (Athens, Greece) _(SPLASH-E 2019)_. Association for Computing Machinery, New York, NY, USA, 13–24.  https://doi.org/10.1145/3358711.3361626
  * Frankle et al. (2016) Jonathan Frankle, Peter-Michael Osera, David Walker, and Steve Zdancewic. 2016.  Example-directed synthesis: a type-theoretic interpretation.  _ACM SIGPLAN Notices_ 51, 1 (April 2016), 802–815.  https://doi.org/10.1145/2914770.2837629
  * Gerdes et al. (2017) Alex Gerdes, Bastiaan Heeren, Johan Jeuring, and L. Thomas van Binsbergen. 2017.  Ask-Elle: an Adaptable Programming Tutor for Haskell Giving Automated Feedback.  _International Journal of Artificial Intelligence in Education_ 27, 1 (01 Mar 2017), 65–100.
  * Gulwani et al. (2018) Sumit Gulwani, Ivan Radiček, and Florian Zuleger. 2018\.  Automated Clustering and Program Repair for Introductory Programming Assignments.  _SIGPLAN Not._ 53, 4 (June 2018), 465–480.  https://doi.org/10.1145/3296979.3192387
  * Hart et al. (1968) P. E. Hart, N. J. Nilsson, and B. Raphael. 1968\.  A Formal Basis for the Heuristic Determination of Minimum Cost Paths.  _IEEE Transactions on Systems Science and Cybernetics_ 4, 2 (1968), 100–107.  https://doi.org/10.1109/TSSC.1968.300136
  * Keuning et al. (2018) Hieke Keuning, Johan Jeuring, and Bastiaan Heeren. 2018\.  A Systematic Literature Review of Automated Feedback Generation for Programming Exercises.  _ACM Trans. Comput. Educ._ 19, 1, Article 3 (Sept. 2018), 43 pages.  https://doi.org/10.1145/3231711
  * Lazar and Bratko (2014) Timotej Lazar and Ivan Bratko. 2014.  Data-Driven Program Synthesis for Hint Generation in Programming Tutors. In _Intelligent Tutoring Systems_ , Stefan Trausan-Matu, Kristy Elizabeth Boyer, Martha Crosby, and Kitty Panourgia (Eds.). Springer International Publishing, Cham, 306–311.
  * Lee et al. (2018) Woosuk Lee, Kihong Heo, Rajeev Alur, and Mayur Naik. 2018\.  Accelerating Search-Based Program Synthesis Using Learned Probabilistic Models.  _SIGPLAN Not._ 53, 4 (June 2018), 436–449.  https://doi.org/10.1145/3296979.3192410
  * Lubin et al. (2020) Justin Lubin, Nick Collins, Cyrus Omar, and Ravi Chugh. 2020\.  Program Sketching with Live Bidirectional Evaluation.  _Proc. ACM Program. Lang._ 4, ICFP, Article 109 (Aug. 2020), 29 pages.  https://doi.org/10.1145/3408991
  * Omar et al. (2019) Cyrus Omar, Ian Voysey, Ravi Chugh, and Matthew A. Hammer. 2019\.  Live Functional Programming with Typed Holes.  _Proc. ACM Program. Lang._ 3, POPL, Article 14 (Jan. 2019), 32 pages.  https://doi.org/10.1145/3290327
  * Osera and Zdancewic (2015) Peter-Michael Osera and Steve Zdancewic. 2015.  Type-and-example-directed program synthesis.  _ACM SIGPLAN Notices_ 50, 6 (Aug. 2015), 619–630.  https://doi.org/10.1145/2813885.2738007
  * Piech et al. (2015) Chris Piech, Mehran Sahami, Jonathan Huang, and Leonidas Guibas. 2015.  Autonomously Generating Hints by Inferring Problem Solving Policies. In _Proceedings of the Second (2015) ACM Conference on Learning @ Scale_ (Vancouver, BC, Canada) _(L@S ’15)_. Association for Computing Machinery, New York, NY, USA, 195–204.  https://doi.org/10.1145/2724660.2724668
  * Price et al. (2016) Thomas W. Price, Yihuan Dong, and Tiffany Barnes. 2016\.  Generating Data-Driven Hints for Open-Ended Programming. In _Proceedings of the 9th International Conference on Educational Data Mining_ , Tiffany Barnes, Min Chi, and Mingyu Feng (Eds.). Raleigh, NC, USA, 191–198.
  * Singh et al. (2013) Rishabh Singh, Sumit Gulwani, and Armando Solar-Lezama. 2013\.  Automated Feedback Generation for Introductory Programming Assignments.  _SIGPLAN Not._ 48, 6 (June 2013), 15–26.  https://doi.org/10.1145/2499370.2462195