# Counterfactual Explanation for Regression via Disentanglement in Latent
Space

Xuan Zhao  Innovation & digitale Transformation
SCHUFA Holding AG
Germany
<EMAIL_ADDRESS>Klaus Broelemann  Innovation & digitale Transformation
SCHUFA Holding AG
Germany
<EMAIL_ADDRESS>Gjergji Kasneci  Responsible Data Science
Technical University of Munich
Germany
<EMAIL_ADDRESS>
###### Abstract

Counterfactual Explanations (CEs) help address the question: How can the
factors that influence the prediction of a predictive model be changed to
achieve a more favorable outcome from a user’s perspective? Thus, they bear
the potential to guide the user’s interaction with AI systems since they
represent easy-to-understand explanations. To be applicable, CEs need to be
realistic and actionable. In the literature, various methods have been
proposed to generate CEs. However, the majority of research on CEs focuses on
classification problems where questions like “What should I do to get my
rejected loan approved?” are raised. In practice, answering questions like
“What should I do to increase my salary?” are of a more regressive nature. In
this paper, we introduce a novel method to generate CEs for a pre-trained
regressor by first disentangling the label-relevant from the label-irrelevant
dimensions in the latent space. CEs are then generated by combining the label-
irrelevant dimensions and the predefined output. The intuition behind this
approach is that the ideal counterfactual search should focus on the label-
irrelevant characteristics of the input and suggest changes toward target-
relevant characteristics. Searching in the latent space could help achieve
this goal. We show that our method maintains the characteristics of the query
sample during the counterfactual search. In various experiments, we
demonstrate that the proposed method is competitive based on different quality
measures on image and tabular datasets in regression problem settings. It
efficiently returns results closer to the original data manifold compared to
three state-of-the-art methods, which is essential for realistic high-
dimensional machine learning applications. Our code will be made available as
an open-source package upon the publication of this work.

###### Index Terms:

XAI, Counterfactual Explaination, Regression

##  I Introduction

Machine learning models have recently become ubiquitous in our society, and
the demand for explainability is on the rise. Especially in high-stakes
applications like healthcare, finance, employment, etc., explainability of the
model’s behavior is a prime prerequisite for establishing trust in machine
learning-based decision systems. Researchers have developed various techniques
to explain the relationships between the input and output of a model. While
models with relatively simple designs (e.g., logistic regression, decision
trees, rule fit algorithm) can be interpreted straightforwardly, more complex
models could be analyzed based on simpler surrogate models. These surrogate
models imitate the more complex model in a post-hoc style locally and provide
explanations in terms of local feature importance (i.e., for a given input).
Another important line of research aims to answer the following question in a
classification setting: “How can the input be changed to achieve a prediction
representing the favored class instead of the unfavored one?”. The exploration
of outcomes in alternative, similar yet non-occurring worlds is called
counterfactual analysis [1]. With a pre-trained and fixed model, the only way
for the model to produce a different output is by altering the input. To this
end, counterfactual explanations (CEs) provide prescriptive suggestions on the
features of the query sample (i.e., input) that have to change (and also by
how much they need to change) to achieve the desired outcome. We require this
change to be minimal (i.e., associated with low costs) and actionable (i.e.,
realistic and feasible). With a clear semantic implication and a common
logical grounding, CEs are generally easy to understand by end-users [2]. In
fact, in the literature, the most direct potential use of CEs is to provide
advice and guidance to end-users.

To generate practical CEs, various requirements are proposed in the
literature: CEs should be actionable, sparse, valid, proximate, and
computationally efficient to generate [3]. Different methods have been
proposed that aim to meet these requirements. It is also clear that there
might be a trade-off between these requirements. However, all methods face the
following problems:

  * •

The high computational costs of the constrained search for adequate CEs,
especially in the case of high-dimensional input space.

  * •

The generation of out-of-sample CEs, which leads to suggestions that are not
actionable (i.e., unlikely or even impossible to be achieved since they do not
correspond to the training data distribution).

  * •

Uncertainties of the user about the implementation of the suggested CEs.

Nonetheless, it is well understood that an actionable CE should be close to
the data manifold and suggest meaningful (i.e., realistic and easily
achievable) changes to a query sample. While most of the existing work on CEs
focuses on classification models, in many practical scenarios, it is quite
essential to answer questions in a regressive way, e.g., “what should I do to
improve my credit score?” Since the model outputs are discrete for
classification, the notion of “change” for a model’s output is unambiguous –
when the decision boundary is crossed. In contrast, regression models have
continuous outputs which can change under arbitrarily small perturbations [4].
In this situation, the existing search algorithms designed for classification
are unsuitable for regression problem settings. A naive way to transform
regression to classification is simply binning the target labels. However,
such a transformation would not address the counterfactual search within the
binned intervals. To address the problems mentioned above, we develop a novel
method to generate CEs for a pre-trained regressor by first disentangling the
label-relevant from the label-irrelevant dimensions with the help of a special
adversarial design (see details in Section III) in the latent space. CEs are
then generated by combining the label-irrelevant dimensions and the predefined
output. We show that our method maintains the characteristics of the query
sample during the counterfactual search. In the problem setting of Figure 1,
we assume the counterfactual question is “what should be changed if we want
the original handwritten digit images with 0 rotational degree to become
images with 45 rotational degree?” According to our method, the label-
irrelevant characteristics like writing styles remain intact when the
rotational angles (label) change during the search. Our method only requires
access to the predictions of a previously trained regressor that we aim to
explain and its training dataset. We demonstrate that it is computationally
less expensive and more robust than other state-of-the-art methods in Section
IV. Our main contributions are: (1) a model-agnostic framework for finding
actionable CEs that lie on the data manifold with low computational expense
for regression problem settings, (2) a novel strategy for manipulating the
latent space to enable efficient counterfactual search, (3) a comparison of
methods across image and tabular datasets.

Figure 1: Gradual changes of digits from rotational 0 to 45 degrees (We
consider a regression task of estimating rotation angles of MINIST images), as
expected, show that the label-irrelevant characteristics like writing styles
do not change during the search.

The remainder of this paper is structured as follows: First, we briefly review
the existing methods of counterfactual explanation in Section II. Then, in
Section III, we propose an adversarial disentanglement approach for computing
CEs. In Section IV, we present the evaluation results of our method against
three state-of-the-art methods. Finally, Section V concludes the paper.

##  II Background and Related Work

CE methods fall into two main categories: model-agnostic or model-specific.
The model-specific methods usually require full access to the specific model,
which needs to be explained, in order to exploit its internal structure. On
the other hand, model-agnostic methods work for arbitrary machine learning
(ML) models and often only require access to the predictions of an already
fitted model that needs to be explained. In this paper, we focus on the model-
agnostic situation. We will present some relevant work below.

###  II-A Generating Counterfactuals by Perturbing the Original Input Space

One large branch of the literature generates counterfactuals by perturbing the
original input feature space. Inverse classification [5] maintains sparsity by
partitioning the features into immutable and mutable features and imposing
budgetary constraints on the mutable features. A sampling approach is proposed
by Laugel et al. [6] using the growing spheres method to traverse the input
space for CEs. Gradient descent is utilized in the input space to find
contrastive explanations [7], which are separated into pertinent positives and
negatives. The authors include an autoencoder loss to keep the explanations
within the data manifold. Gradient descent methods are improved with the
introduction of prototypes, guiding the gradient descent towards the average
value of the target class by averaging the representations of the training set
in the latent space [8]. CE generation with reinforcement learning [9], on the
other hand, proposes a fast, model-agnostic strategy by replacing the usual
optimization procedure with a learnable process. Moreover, GRACE [10] is
designed for neural networks on tabular data and combines the concepts of
contrastive explanations with interventions by performing constrained gradient
descent and adding a loss that measures information gain to keep the resulting
explanations sparse. Many of the methods mentioned above already include
generative models to obtain CEs on the data manifold. However, perturbing
inputs without proper regularization can generate unconvincing and infeasible
CEs [11] that resemble adversarial samples. Practical desiderata for the
generated CEs, such as actionability, sparsity, validity, proximity, and
computational efficiency, are included in the methods as constraints of the
optimization design, which might lead to complicated optimization, especially
when the input is high-dimensional.

###  II-B Generating Counterfactuals by Perturbing the Latent Space

Latent space perturbation methods utilize generative and probabilistic models
in algorithm design to ensure that CEs have a high probability under the data
distribution $p(X)$. StylEx [12] includes the classifier in the generative
model and manipulates the latent space to visualize the counterfactual search.
ExplainGAN [13] is a method for finding CEs for images by training multiple
autoencoders and using the signal from the classifier and discriminators to
inform the learned representations. Growing sphere search is performed in the
latent space of a conditional variational autoencoder to generate
counterfactuals [14]. Along similar lines, Balasubramanian et al. [15] employ
gradient descent in the latent space of a variational autoencoder with
regularization terms. Sharpshooter [16] searches for counterfactuals by linear
interpolation in the latent space with aid from two separate autoencoders.
Xspells [17] searches through the latent space with the help of a latent
decision tree for CEs for short text classification. However, searching in the
latent space does not guarantee sparsity in the input space, as we show in
Section IV.

Figure 2: Diagram of the data flow in the algorithm. During Training: an input
sample $x$, which receives a regression prediction $\hat{y}$ from the pre-
trained regressor $f$. The sample $x$ is embedded into a latent space learned
by the Encoderu with an Adversarial Regressor that helps remove the
information related to the regression label $\hat{y}$. A Discriminator aims to
discriminate the real sample $x$ from the reconstructed $x^{\prime}$ by the
Decoder. During Generating: with a query sample $x_{q}$ with a regression
label $\hat{y}_{q}$ (e.g., a rotation of 10 degrees), the label-irrelevant
embedding $z_{u}$ is produced by running $x_{q}$ through the fixed Encoderu
from the Training Step. Then a linear interpolation is performed by combining
the predefined output $\hat{y}_{t}$ (e.g., an output of 45 degrees) and
$\hat{y}_{q}$ to get $y_{CE}$. $y_{CE}$ and $z_{u}$ are sent through the fixed
Decoder from the Training Step to get a potential CE. The potential CE goes
through the pre-trained regressor $f$ to get a regression score. This process
stops until the predefined output $\hat{y}_{t}$ is reached.

##  III An Adversarial Disentanglement Approach

###  III-A Problem Setting

In this paper, we consider generating CEs in a regression setting. We are
given a query sample (i.e., a designated input instance) for which we seek an
actionable counterfactual explanation toward a predefined output.

Formally, the training dataset of the pre-trained regressor
$\hat{y_{i}}=f(x_{i})$ is given by $\mathcal{D}=(x_{i},y_{i})_{i=1}^{K}$,
where $y_{i}\in\mathbb{R}$ is a continuous label and $\hat{y_{i}}$ is the
regression output/prediction. $x_{q}$ is a query sample with an output
$\hat{y}_{q}$ from the regressor $f$. $\hat{y}_{t}$ is the predefined output.
Hence, CEs are needed for the query sample $x_{q}$ under the regressor $f$. We
consider the CEs for regression defined as below.

Definition 1 During the CE search, if the regression output just reaches the
_predefined output_ $\hat{y}_{t}$, $X_{CE}$ is the CE for the regressor $f$
with a specific predefined output $\hat{y}_{t}$: $f(X_{CE})=\hat{y}_{t}$.

###  III-B Desiderata for CE Search

Our approach is driven by the following main desiderata:

Desideratum 1 Intrinsic and label-irrelevant characteristics of the input
query should be retained during the search process, i.e., the favoured label
should be reached based on intrinsic properties of the input and with low
effort.

Desideratum 2 The generation step, should return highly realistic, i.e.,
feasible/actionable CEs.

Desideratum 3 The CEs should be retrieved efficiently also in high-dimensional
settings to ensure the practical application to real-life settings.

###  III-C Learning Components

We design a CE generation method based on two main steps: a Training Step and
a Generating Step. It only requires access to the training dataset and
prediction of the pre-trained regressor $f$ that we aim to explain through
appropriate CEs. For our Training Step, we form a new training dataset
$\mathcal{D}_{t}=(x_{i},\hat{y}_{i})_{i=1}^{K}$ by replacing original $y$ in
$\mathcal{D}$ with $\hat{y}$. We then train (1) an Encoderu that captures
label-irrelevant embedding $z_{u}$ of the original input $x$, (2) an
Adversarial Regressor helps remove label-relevant information from the
embedding $z_{u}$, (3) a Decoder that takes $z_{u}$ and $\hat{y}$ to output
reconstruction $x^{\prime}$ and (4) a Discriminator aims to discriminate the
original input $x$ from the reconstruction $x^{\prime}$ (see Figure 2). In the
Generating Step, we combine the label-irrelevant dimensions $z_{u}$ and the
predefined output $\hat{y}_{t}$ and project them back to the original space to
generate a relatively less expensive and more realistic CE. The algorithm flow
is shown in Figure 2. In the following part of this subsection, we will
describe each component of the Training Step in Figure 2. The details of
Generating Step are described in Section III-D.

Encoderu We notice that generative models usually try to generate various
unseen samples. However, the generation of CEs needs to satisfy different
criteria, as mentioned in Section I. In our problem setting, to satisfy
_Desideratum 1_ which requires maintaining the characteristics of the query
sample during the search, it is intuitive to adopt disentanglement methods in
the latent space of an autoencoder.

Adversarial Regressor Inspired by a Two-Step Disentanglement Method [18], we
introduce an Adversarial Regressor to ensure that the embedding $z_{u}$
captured by the Encoderu is regression label-irrelevant. The motivation for us
to utilize this component is that a CE should maintain the query sample’s
characteristics according to _Desideratum 1_. It is inspired by GANs, where
the discriminator gradually loses the capacity to tell the generated data from
the real data during the training phase. While GANs are usually used to
improve the quality of generated output – telling fake from real, the
adversarial component encourages the Encoderu to dismiss information about the
labels, leading to disentanglement. With the Adversarial Regressor, we could
guarantee that the $z_{u}$ and $\hat{y}$ are disentangled and independent,
which prepares for the interpolation in the Generating Step.
$\hat{y}^{\prime}$ is the regression label of $z_{u}$ through the Adversarial
Regressor. The adversarial regression loss is shown in Equation 1.

$\mathcal{L}_{Adv}=\frac{1}{N}\sum_{i=1}^{N}||\hat{y}_{i}-\hat{y}^{\prime}_{i}||_{2}^{2}$
(1)

Decoder As indicated in Figure 2, we can see that the Decoder takes regression
label $\hat{y}$ and the embedding $z_{u}$ as inputs. In the Training Step,
with the training going on, we aim to reconstruct a close sample to the query
sample. We show in Figure 4(a) how the label is added to the architecture of
the Decoder. $x^{\prime}_{i}$ is the reconstruction after the Decoder. The
reconstruction loss of the Encoderu and Decoder is shown in Equation 2.

$\mathcal{L}_{Rec}=\frac{1}{N}\sum_{i=1}^{N}||x_{i}-x^{\prime}_{i}||_{2}^{2}$
(2)

Discriminator To guarantee the reconstruction $x^{\prime}$ is realistic and
close the training data distribution to obtain actionable CEs (_Desideratum
2_), we borrow the ideas of GANs again. Existing conditional GANs (cGANs) are
mainly used for class labels, and the discriminator’s loss function is shown
in Equation 3. Inspired by the CcGAN architecture introduced in recent work
[19], we aim to use the discriminator in CcGAN to assure the quality of
reconstruction $x^{\prime}$. In CcGAN, a reformulation of the loss function of
cGAN based on kernel density estimates is introduced, which leads to two
improved empirical discriminator losses for regression labels based on the
cGAN for classification – the hard vicinal discriminator loss (HVDL) and the
soft vicinal discriminator loss (SVDL). After checking the reported
performance in the paper of CcGAN, we find that the performance between the
two loss functions is similar.

Hence, we adopt HVDL by plugging $\hat{p}^{HVE}_{r}(x,y)$ and
$\hat{p}^{HVE}_{g}(x,y)$ in Equations 4 and 5 back to Equation 3 to get a
discrimination loss $\mathcal{L}_{D}$. We also replace all labels $y$ in
Equations 3, 4 and 5 by regression labels $\hat{y}$ (e.g., $y_{i}^{r}$ to
$\hat{y}_{i}^{r}$ and $y_{i}^{g}$ to $\hat{y}_{i}^{g}$) because we perform the
Training Step on the new training dataset
$\mathcal{D}_{t}=(x_{i},\hat{y}_{i})_{i=1}^{K}$ instead of the original
dataset $\mathcal{D}$. How the regression label is input into the network of
the discriminator is shown in Figure 4(b).

$x_{i}^{r}$ and $x_{i}^{g}$ are the original/real and the reconstructed/fake
$i^{\prime}$th sample; $y_{i}^{r}$ and $y_{i}^{g}$ are the corresponding
regression labels of $x_{i}^{r}$ and $x_{i}^{g}$; $k$ and $\sigma$ are two
positive hyper-parameters, $C_{1}$ and $C_{2}$ are two normalizing factors
ensuring $\hat{p}^{HVE}_{r}$ and $\hat{p}^{HVE}_{g}$ valid probability density
functions; $N^{r}_{y,k}$ is the number of the labels $y_{i}^{r}$ satisfying
$\lvert y-y^{r}_{i}\lvert\leqslant{k}$; $N^{g}_{y,k}$ is the number of the
labels $y_{i}^{g}$ satisfying $\lvert y-y^{g}_{i}\lvert\leqslant{k}$. The
terms in the first square brackets of $\hat{p}^{HVE}_{r}$ and
$\hat{p}^{HVE}_{g}$ show that we estimate the marginal label distributions
($p_{r}(y)$ and $p_{g}(y)$) through kernel density estimates. The terms in the
second square brackets are based on the assumption that a small perturbation
around $y$ leads to negligible changes to $p_{r}(x|y)$ and $p_{g}(x|y)$.

$\displaystyle\begin{aligned} \mathllap{\mathcal{L}_{D}}&=-\mathds{E}_{y\sim
p_{r}(y)}[\mathds{E}_{x\sim p_{r}({x|y})}[\textrm{log}D(x,y)]]\\\
&\quad-\mathds{E}_{y\sim p_{g}(y)}[\mathds{E}_{x\sim
p_{g}({x|y})}[1-\textrm{log}D(x,y)]]\\\
&=-\int{\textrm{log}(D(x,y))p_{r}(x,y)dxdy}\\\
&\quad-\int{\textrm{log}(1-D(x,y))p_{g}(x,y)dxdy}\end{aligned}$ (3)
$\displaystyle\begin{aligned}
\mathllap{\hat{p}^{HVE}_{r}(x,y)}&=C_{1}\left[\dfrac{1}{N^{r}}\sum_{j=1}^{N^{r}}\textrm{exp}\left(-\dfrac{(y-y_{j}^{r})^{2}}{2\sigma^{2}}\right)\right]\\\
&\quad\left[\dfrac{1}{N^{r}_{y,k}}\sum_{i=1}^{N^{r}}\mathbb{1}_{\\{\lvert
y-y^{r}_{i}\lvert\leqslant{k}\\}}\delta(x-x^{r}_{i})\right]\end{aligned}$ (4)
$\displaystyle\begin{aligned}
\mathllap{\hat{p}^{HVE}_{g}(x,y)}&=C_{2}\left[\dfrac{1}{N^{g}}\sum_{j=1}^{N^{g}}\textrm{exp}\left(-\dfrac{(y-y_{j}^{g})^{2}}{2\sigma^{2}}\right)\right]\\\
&\quad\left[\dfrac{1}{N^{g}_{y,k}}\sum_{i=1}^{N^{g}}\mathbb{1}_{\\{\lvert
y-y^{g}_{i}\lvert\leqslant{k}\\}}\delta(x-x^{g}_{i})\right]\end{aligned}$ (5)

Summary The configuration of the network in the Training Step is composed of
three network branches: first, the Adversarial Regressor is trained to
minimize the adversarial regression loss $\mathcal{L}_{Adv}$ in Equation 1 –
it is trained to regress $z_{u}$ to $\hat{y}$. Second, the autoencoder network
is trained to minimize the total loss $\mathcal{L}$, the sum of three terms as
shown in Equation 6: (i) the reconstruction loss $\mathcal{L}_{Rec}$ as shown
in Equation 2, (ii) minus the adversarial regression loss $\mathcal{L}_{Adv}$
in Equation 1, and (iii) minus the discrimination loss $\mathcal{L}_{D}$ in
Equation 3 (with adjustment). Third, the discriminator network is trained to
discriminate the reconstructed from the real samples by minimizing the
discrimination loss $\mathcal{L}_{D}$.

$\mathcal{L}=\mathcal{L}_{Rec}-\lambda_{Adv}\mathcal{L}_{Adv}-\lambda_{D}\mathcal{L}_{D}$
(6)

###  III-D Algorithm

Algorithm 2 and 2 show the pseudo-code of the Training and Generating Steps in
Figure 2, which is used to generate the CEs in Section IV.

Algorithm 1 Training Step of our proposed architecture

1:$\psi$, $\phi$ , $\omega$ and $\sigma$ the initial parameters of Encoderu,
Decoder, Adversarial Regressor and Discriminator; $n$ the number of
iterations; $\lambda_{Adv}$ the weight of regularization term
$\mathcal{L}_{Adv}$

2:while not converged do

3: for $i=0$ to $n$ do

4: Sample $\\{x,\hat{y}\\}$ a batch from dataset $\mathcal{D}$.

5:
$\omega\stackrel{{\scriptstyle+}}{{\leftarrow}}\bigtriangledown_{\omega}\mathcal{L}_{Adv}$

6:
$\sigma\stackrel{{\scriptstyle+}}{{\leftarrow}}\bigtriangledown_{\sigma}\mathcal{L}_{D}$

7: Sample $\\{x,\hat{y}\\}$ a batch from dataset $\mathcal{D}$.

8:
$\mathcal{L}\leftarrow\mathcal{L}_{Rec}-\lambda_{Adv}\mathcal{L}_{Adv}-\lambda_{D}\mathcal{L}_{D}$

9:
$\psi,\phi\stackrel{{\scriptstyle+}}{{\leftarrow}}\bigtriangledown_{\psi,\phi}\mathcal{L}$

10: end for

11:end while

Algorithm 2 Generating Step of our proposed architecture

1:$S$ samples $\alpha=\alpha_{s=1}^{S}$ in $(0,1]$; $f$ regressor; $x_{q}$ the
query sample; $\hat{y}_{q}$ regression label of the query sample;
$\hat{y}_{t}$ predefined regression output; $tol$ tolerance; $x_{s}$ potential
CE

2:for $\alpha_{s}$ in $\alpha$ do

3: $z_{u}\leftarrow Encoder\textsubscript{u}(x_{q})$

4: $\hat{y}_{CE}=(1-\alpha)\hat{y}_{q}+\alpha\hat{y}_{t}$

5: $x_{s}\leftarrow Decoder(z_{u},\hat{y}_{s})$

6: if $|f(x_{s})-\hat{y}_{CE}|<tol$ then

7: $x_{CE}=x_{s}$

8: end if

9:end for

10:return $x_{CE}$

For our Training Step, we sample a batch from $\mathcal{D}_{t}$ to update the
parameters of the Adversarial Regressor and Discriminator ($\omega$ and
$\sigma$). Then we sample a batch from $\mathcal{D}_{t}$ to update the
parameters: $\psi$ of the Encoderu and $\phi$ of the Decoder under the loss
function introduced by Equation 6. We iterate this procedure until convergence
of the loss functions is reached. At the end of the Training Step: we should
get (1) an Encoderu that could extract the regression label-irrelevant
dimensions $z_{u}$ from the query sample, (2) an Adversarial Regressor that
could not tell the regression label $\hat{y}$ of a given input $x$, (3) a
Decoder that takes the input of the encoded $z_{u}$ and the regression label
$\hat{y}$ to reconstruct a realistic sample $x^{\prime}$, and (4) a
Discriminator which fails to tell the reconstructed samples $x^{\prime}$ from
the real input $x$.

In our Generating Step, we are motivated to do a linear interpolation because
it is relatively less computationally expensive than other existing
perturbation methods (_Desideratum 3_). We first pass a query sample $x_{q}$
with an unfavorable regression output $\hat{y}_{q}$ through the Encoderu to
obtain the sample’s label irrelevant embedding $z_{u}$. Then we generate CEs
by increasing or decreasing the query regression label $\hat{y}_{q}$ according
to the task (increase or decrease the label toward the predefined output
$\hat{y}_{t}$) using linear interpolation in latent space
$y_{CE}=(1-\alpha)\hat{y}_{q}+\alpha\hat{y}_{t}$ with $\alpha$ in $(0,1]$.
$y_{CE}$ and $z_{u}$ are then projected back to the original space through the
Decoder to get the potential CE $x_{s}$, and the pre-trained regressor $f$
assesses the regression score of the potential CE $x_{s}$. The counterfactual
search stops if $x_{s}$ satisfies the predefined regression output (within a
user-specified tolerance $tol$). The search is performed by sampling between
$\hat{y}_{q}$ and $\hat{y}_{t}$ with a finite number of $\alpha$, which is
quite efficient compared to other perturbing methods.

##  IV Experiments and Evaluation

We compare our method to three other counterfactual methods introduced in
Section II: Gradient Descent Method improved with Prototypes (Prototype) [8],
CE Generation with Reinforcement Learning (RL) [9] and Gradient Descent in the
Latent Space of a VAE (GDL) [15] on three datasets: MNIST [20] (image), Car
Price Prediction [21] (tabular), and House Sales in King County [22]
(tabular). The reason why we choose Prototype and RL is that they are both
designed not only for image datasets but also for tabular datasets. Besides,
they both use autoencoders to learn the representation to remain close to the
data distribution, similar to our design, while they use other constraints to
ensure desiderata like sparsity. We use the package ALIBI [23] for Prototype
and RL to stay close to the original design. Since the application of
Prototype and RL is for classification in ALIBI, we bin the continuous target
labels into 20 classes. GDL is a relatively simple baseline, but it is similar
to our method because both operate in latent space. We implement GDL from
scratch since the source code is not available. For the experiments in this
Section, the predefined output $\hat{y}_{t}$ is set to be
$\hat{y}_{t}-\hat{y}_{q}=0.2$ (after normalization) for further comparison
among different CE methods. Our test are run on an Intel(r) Core(TM) i7-8700
CPU. The networks in the experiments are built based on Pytorch [24].

###  IV-A Measures for Comparison

We evaluate the quality of CEs by the measures taken from literature [3, 16,
25]: (i) counterfactual generation time, time required to find a CE for a
given query sample (ii) validity, the percentage success in generating CEs
that the predefined outputs by the users are reached (iii) proximity,
distance($L_{2}$ norm) from query samples to CEs in original space (iv)
sparsity, $L_{1}$ norm of the change vector in original space (v)
reconstruction loss \- a measure of the CE being in sample. We pass a CE
through the autoencoder and measure its loss. A smaller loss indicates closer
to the original data distribution because the autoencoder is trained on the
same training dataset as the pre-trained regressor. Further details about the
metrics used in the experiments can be found in the Appendix (-A).

Figure 3: Original and counterfactual samples for MNIST showing from top to
bottom: original query images and CEs found via Our Method, GDL, Prototype and
RL. Note Our Method exhibits more of the characteristics of the predefined
output (further rotation) than the other methods. Our Method and GDL remain
(perturbation in latent space) more realistic than the Prototype and RL
(perturbation in original space). Additionally, only Our Method captures the
higher-level concept like rotation.

###  IV-B Datasets, Training and Evaluation

(a) The label input method for the Decoder

(b) The label input method for the Discriminator

Figure 4: The label $y$ input methods of the Decoder and Discriminator for
MNIST (for tabular datasets, multi-layer perceptrons are used)

MNIST The MNIST dataset [26] is an extensive database of handwritten digits.
The original MNIST dataset’s target label is suitable for the multi-
classification problem. Here, we appeal to its visualization and randomly
rotate the original MNIST handwritten digits from 0 to 45 degrees since it is
not easy to find a more proper image dataset with a regression task. We
generate 20000 handwritten digit images with rotational angles randomly from 0
to 45 degrees as training and 5000 as testing sets. We adjust the problem to a
regression predicting the rotational angle of the images. MNIST provides a
naturally intuitive visualization of the result – we could show a gradual
change from a certain query rotational angle to the predefined output of
rotational angle requested by the users.

For MNIST, we use a CNN-based network (Figure 4). We train the model for 400
epochs by stochastic gradient descent using the Adam optimizer and a batch
size of 200. Our model uses two hyperparameters, $\lambda_{Adv}$ and
$\lambda_{D}$, for regularizing the reconstruction loss function for the
network. These weigh the effect of the adversarial losses from both the
Adversarial Regressor and Discriminator. We varied $\lambda_{Adv}$ between 0.1
and 5 during training. During the training process, we want the Encoderu and
Adversarial Regressor ‘grow’ together. We found $\lambda_{Adv}=0.5$ to give us
the highest harmonic training process of the Encoderu, Adversarial Regressor
and Decoder. We use this value to report our results. Changing $\lambda_{D}$
while keeping $\lambda_{Adv}$ fixed at 0.5 does not result in any significant
changes, as these changes were largely overshadowed by the gains from changing
$\lambda_{Adv}$. Hence, we set $\lambda_{D}=1$ for our experiments. The
details of the hyper-parameters of the Discriminator are shown in the Appendix
(-C).

The CEs generated by our method shown in Figure 1 exhibit a combination of
characteristics from the query sample (e.g., the thickness of the stroke,
writing styles) and from the predefined output (e.g., a different rotational
angle). Our _Desideratum 1_ of keeping the characteristics of the query sample
which are not related to the predefined output (mentioned in Section III-B) is
reached. Figure 3 also shows that only our method captures the higher-level
characteristics like rotation. Table I (mean±SD) shows the quality of
counterfactual explanations measured by the above metrics. Our method
outperforms the other methods in counterfactual generation time,
reconstruction and validity but not always in proximity. Our Method and GDL
remain (perturbation in latent space) more realistic than Prototype and RL
(perturbation in original space). It also indicates that for high-dimensional
applications, it is important to avoid the situation of adversarial samples
mentioned in Section II, especially for methods perturbing the original space.

TABLE I: Summary of metrics for MNIST

Method | time(s) | reconstruction | sparsity | validity(%) | proximity
---|---|---|---|---|---
Our Method | 0.044$\pm$0.005 | 0.008$\pm$0.001 | 2.343$\pm$0.627 | 90.3 | 0.896$\pm$0.245
GDL | 5.803$\pm$0.224 | 0.520$\pm$0.135 | 0.090$\pm$0.002 | 84.2 | 2.170$\pm$0.835
Prototype | 7.345$\pm$0.784 | 4.463$\pm$0.563 | 0.014$\pm$0.005 | 75.3 | 1.075$\pm$0.235
RL | 1.804$\pm$0.112 | 5.453$\pm$0.673 | 0.012$\pm$0.008 | 65.3 | 0.875$\pm$0.132

TABLE II: Summary of metrics for Car Price Prediction

Method | time(s) | reconstruction | sparsity | validity(%) | proximity
---|---|---|---|---|---
Our Method | 0.007$\pm$0.001 | 0.012$\pm$0.003 | 2.119 $\pm$0.857 | 89.5 | 0.193 $\pm$0.075
GDL | 1.568 $\pm$0.047 | 0.247 $\pm$0.048 | 2.121 $\pm$0.085 | 70.9 | 0.172$\pm$0.068
Prototype | 1.345 $\pm$0.095 | 0.926 $\pm$0.175 | 0.014$\pm$0.005 | 58.7 | 1.075$\pm$0.055
RL | 0.085 $\pm$0.012 | 0.857$\pm$0.058 | 0.015$\pm$0.026 | 60.5 | 0.834$\pm$0.047

TABLE III: Summary of metrics for House Sales in King County

Method | time(s) | reconstruction | sparsity | validity(%) | proximity
---|---|---|---|---|---
Our Method | 0.008$\pm$0.002 | 0.247$\pm$0.155 | 1.713 $\pm$0.973 | 91.2 | 3.534$\pm$0.892
GDL | 3.465 $\pm$0.765 | 0.475 $\pm$0.406 | 1.739 $\pm$0.895 | 71.2 | 1.953$\pm$0.876
Prototype | 2.164$\pm$1.973 | 0.859 $\pm$0.675 | 0.871$\pm$0.105 | 67.5 | 1.402$\pm$0.379
RL | 0.095 $\pm$0.006 | 0.870 $\pm$0.501 | 0.375$\pm$0.159 | 85.3 | 0.964$\pm$0.105

Car Price Prediction The Car Price Prediction dataset contains features about
used cars which could be used for price prediction. For simplicity, we treat
all the features as mutable. We train a regressor on the selling prices with
Age of the cars, Number of Kilometres the car is driven, Number of seats as
continuous features, and Fuel type of car (petrol/diesel/CNG/LPG/electric) as
categorical features.

For tabular datasets, more prepossessing is usually performed compared to
image datasets. We normalize the continuous features and use one-hot encoding
to deal with the categorical features. We train the model for 300 epochs by
stochastic gradient descent using the Adam optimizer and a batch size of 100.
$\lambda_{Adv}$ and $\lambda_{D}$ are set at 0.5 and 0.5 for the best
performance. For tabular datasets, we use a multi-layer perceptron-based
network. At the end of our training step, as we expected, the Adversarial
Regressor gradually outputs the mean of the prediction label $\hat{y}$
(predicted normalized car prices), which means the Encoderu extracts the label
irrelevant information of the input $x$ as designed. The Adversarial Regressor
could only ‘guess’ an output.

Based on the comparison results shown in Table II, our method outperforms in
terms of time, validity, and reconstruction dimensions, while Prototype/RL is
stronger in sparsity and GDL in proximity.

House Sales in King County This dataset contains house sale prices for King
County, Seattle. It includes homes sold between May 2014 and May 2015. This
tabular dataset contains information about the number of bedrooms, age of the
house, number of bathrooms, square feet of living area, grade of the house,
floor, and other features. We trained a regressor to predict the house price
using five continuous and one categorical feature. From Table III, we find
that, similar to the Car Price Prediction dataset, our method is much faster
at generating CEs and excels in realism, which means it can provide more
actionable CEs as suggestions. The training details of this dataset are
similar to the Car Price Prediction dataset, and more details are included in
the Appendix.

##  V Conclusion, Limitations and Future Work

In this work, we presented a novel model-agnostic algorithm for finding CEs
via disentanglement and interpolation in the latent space for regression
models. Our method implements a framework that combines label-irrelevant
dimensions and the desired predefined output to generate actionable CEs for a
given query sample. We demonstrated the advantages and disadvantages of our
method by comparing it to three state-of-the-art methods on three different
datasets (MNIST, Car Price Prediction, and House Sales in King County). The
evaluation results suggest that our method provides more valid CEs that are
closer to the original data distribution in a more efficient way. We also
provide a visual path of how the features change while moving towards a more
favorable predefined output. Based on the presented comparison, we see further
opportunities for our work to focus on improving the latent space structure to
achieve better representations. Since the search is in the latent space, the
sparsity of our approach is low in the original input space. However, this
aligns with recent findings [27] that CEs from dense regions of the data
distribution tend to have higher validity and are less sparse than those from
sparse regions. Future work, however, should focus on improving the sparsity
of the CEs generated by generative models.

##  VI Acknowledgement

This work has received funding from the European Union’s Horizon 2020 research
and innovation programme under Marie Sklodowska-Curie Actions (grant agreement
number 860630) for the project “NoBIAS - Artificial Intelligence without
Bias”.

## References

  * [1] P. Menzies and H. Beebee, “Counterfactual Theories of Causation,” in _The Stanford Encyclopedia of Philosophy_ , winter 2020 ed., E. N. Zalta, Ed. Metaphysics Research Lab, Stanford University, 2020.
  * [2] C. Fernandez, F. Provost, and X. Han, “Full Text: Explaining Data-Driven Decisions made by AI Systems: The Counterfactual Approach,” https://onikle.com/articles/48299, Jan. 2020.
  * [3] S. Verma, J. Dickerson, and K. Hines, “Counterfactual Explanations for Machine Learning: A Review,” _arXiv:2010.10596 [cs, stat]_ , Oct. 2020.
  * [4] T. Spooner, D. Dervovic, J. Long, J. Shepard, J. Chen, and D. Magazzeni, “Counterfactual Explanations for Arbitrary Regression Models,” Jun. 2021\.
  * [5] M. T. Lash, Q. Lin, W. N. Street, and J. G. Robinson, “A budget-constrained inverse classification framework for smooth classifiers,” _arXiv:1605.09068 [cs, stat]_ , Jun. 2017.
  * [6] T. Laugel, M.-J. Lesot, C. Marsala, X. Renard, and M. Detyniecki, “Inverse Classification for Comparison-based Interpretability in Machine Learning,” _arXiv:1712.08443 [cs, stat]_ , Dec. 2017.
  * [7] A. Dhurandhar, P.-Y. Chen, R. Luss, C.-C. Tu, P. Ting, K. Shanmugam, and P. Das, “Explanations based on the missing: Towards contrastive explanations with pertinent negatives,” in _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , ser. NIPS’18. Red Hook, NY, USA: Curran Associates Inc., Dec. 2018, pp. 590–601.
  * [8] A. Van Looveren and J. Klaise, “Interpretable Counterfactual Explanations Guided by Prototypes,” _arXiv:1907.02584 [cs, stat]_ , Feb. 2020.
  * [9] R.-F. Samoilescu, A. Van Looveren, and J. Klaise, “Model-agnostic and Scalable Counterfactual Explanations via Reinforcement Learning,” Jun. 2021.
  * [10] T. Le, S. Wang, and D. Lee, “GRACE: Generating Concise and Informative Contrastive Sample to Explain Neural Network Model’s Prediction,” _arXiv:1911.02042 [cs, stat]_ , Oct. 2020.
  * [11] Y. Goyal, Z. Wu, J. Ernst, D. Batra, D. Parikh, and S. Lee, “Counterfactual Visual Explanations,” _arXiv:1904.07451 [cs, stat]_ , Jun. 2019.
  * [12] O. Lang, Y. Gandelsman, M. Yarom, Y. I. Wald, G. Elidan, A. Hassidim, B. Freeman, P. Isola, A. Globerson, M. Irani, and I. Mosseri, “Explaining in Style: Training a GAN to explain a classifier in StyleSpace,” in _Proc. ICCV 2021_ , 2021.
  * [13] P. Samangouei, A. Saeedi, L. Nakagawa, and N. Silberman, “ExplainGAN: Model Explanation via Decision Boundary Crossing Transformations,” in _ECCV (10)_ , Jan. 2018.
  * [14] M. Pawelczyk, J. Haug, K. Broelemann, and G. Kasneci, “Learning Model-Agnostic Counterfactual Explanations for Tabular Data,” _Proceedings of The Web Conference 2020_ , pp. 3126–3132, Apr. 2020.
  * [15] R. Balasubramanian, S. Sharpe, B. Barr, J. Wittenbach, and C. B. Bruss, “Latent-CF: A Simple Baseline for Reverse Counterfactual Explanations,” _arXiv:2012.09301 [cs]_ , Jun. 2021.
  * [16] B. Barr, M. R. Harrington, S. Sharpe, and C. B. Bruss, “Counterfactual Explanations via Latent Space Projection and Interpolation,” _arXiv:2112.00890 [cs]_ , Dec. 2021.
  * [17] O. Lampridis, L. State, R. Guidotti, and S. Ruggieri, “Explaining short text classification with diverse synthetic exemplars and counter-exemplars,” _Machine Learning_ , 2022.
  * [18] N. Hadad, L. Wolf, and M. Shahar, “A Two-Step Disentanglement Method,” _arXiv:1709.00199 [cs, stat]_ , Jan. 2020.
  * [19] X. Ding, Y. Wang, Z. Xu, W. J. Welch, and Z. J. Wang, “Continuous Conditional Generative Adversarial Networks: Novel Empirical Losses and Label Input Mechanisms,” Aug. 2021.
  * [20] C. C. Yann LeCun and C. Burges, “MNIST handwritten digit database,” http://yann.lecun.com/exdb/mnist/, 2010.
  * [21] V. Thackerm, “Car Price Prediction,” https://kaggle.com/code/vatsalthacker22/car-price-prediction, 2021.
  * [22] harlfoxem, “House Sales in King County, USA,” https://www.kaggle.com/datasets/harlfoxem/housesalesprediction, 2016.
  * [23] J. Klaise, A. Van Looveren, G. Vacanti, and A. Coca, “Alibi Explain: Algorithms for Explaining Machine Learning Models,” pp. 1–7, Jun. 2021\.
  * [24] A. Paszke _et al._ , “Pytorch: An imperative style, high-performance deep learning library,” in _NeurIPS_ , 2019, pp. 8024–8035.
  * [25] S. Wachter, B. Mittelstadt, and C. Russell, “Counterfactual Explanations without Opening the Black Box: Automated Decisions and the GDPR,” Mar. 2018.
  * [26] Y. LeCun and C. Cortes, “MNIST handwritten digit database,” 2010. [Online]. Available: http://yann.lecun.com/exdb/mnist/
  * [27] M. Pawelczyk, K. Broelemann, and G. Kasneci, “On counterfactual explanations under predictive multiplicity,” in _Conference on Uncertainty in Artificial Intelligence_. PMLR, 2020, pp. 809–818.

###  -A Metrics for Comparison

1) conterfactual generation time - time required to find a counterfactual for
a query sample. Low computational expense is one of the largest strengths of
our method. We report the average time spent generating counterfactual
explanations for a query sample. Our approach is particularly fast compared to
other methods because generating a CE requires a search through a relatively
lower dimension and projection through the Decoder for no more than a fixed
number of interpolated points (search stops when criteria are met). By
comparison, iterative search based on gradients or perturbing in the input
space can be quite expensive if the search distance is large, the learning
rate is low, or the mutable dimensions are high. We do not include the time of
training the Auto-encoder-based architecture with an Adversarial Regressor and
a Discriminator because even though this training time is not trivial, it is a
one-time cost and is not related to the scalability of CEs generation for a
query sample. Ideally, our method triumphs even more, when the query sample is
high dimensional.  2) validity - the percentage of success in generating
counterfactuals that indeed meet the stopping search criteria. The goal for
generating CEs is that they are actually CEs, which requires them to be
regressed as the target scoring. We evaluate validity by including the pre-
trained regressor in the algorithm. Since linear interpolation ends when
certain criteria are met for the regressor, we could see, in general, that our
method performs better at this dimension. Our method uses linear interpolation
to find a proper CE, while the baselines are designed with optimization for
each query sample. When the optimization is designed with complicated
constraints, it is intuitive that the validity is lower since valid solutions
might not be available.  3) proximity - a measure of distance from query
sample to counterfactual in original space. Proximity is usually used as a
constraint to promote the performance of counterfactual generation [3]. This
metric measures how close a counterfactual is to the query sample it is
generated from. Usually this measure is calculated by the $L2$ norm in input
space.  4) sparsity - a measure of features changed. Sparsity is a common
metric in the counterfactual literature [3]. Sparsity is a measure of how
sparse the change vector is. We calculate $L1$ norm of the change vector in
input space).  5) reconstruction loss - a measure of the CE being close to the
data manifold. We can measure the closeness of a CE to its original dataset by
passing a counterfactual through the autoencoder and measuring its
reconstruction loss. During the training process of the autoencoder, we try to
minimize the reconstruction loss. Intuitively, a sample with a smaller loss
should be more in-sample and closer to the original data distribution. An
unseen sample should have a larger reconstruction loss to the autoencoder.

###  -B Further Thoughts about the Comparison Results

We could see that our method performs better on three dimensions above:
counterfactual generation time, validity and reconstruction loss due to its
design. We design the architecture based on our three desiderata and in the
end, through the experiment, we show that the desiderata are achieved. Our
method performs linear interpolation in the latent space instead of
optimization in the original space, it has the highest validity, and also
because the linear interpolation happens in latent space, we cannot guarantee
the sparsity and proximity in the original space, which is stated in Section
II.

###  -C Experiment Details of MNIST

The regression model used to regress the modified handwritten digits into 0 to
45 rotational angles is a CNN-based neural network with convolutional layers,
fully connected layers and ReLU activations in each layer. The model is
trained on standardized features with stochastic gradient descent using Adam
optimizer for 100 epochs with batch size 200.

The Encoder, Adversarial Regressor and Decoder of our design for MNIST are
also CNN-based neural networks with convolutional layers, fully connected
layers and ReLU activation in each layer. The number of latent dimension is
25.

As for the hyper-parameter tuning, please refer to Figure 5 About the
Discriminator, we follow the instructions from the original CcGAN paper, and
set the $\sigma=0.047$ and $k=0.015$.

Figure 5: Pareto plot of $-L_{adv}$ versus $-L_{d}$ for the hyperparameter
search for the MNIST dataset.

###  -D Experiment Details of Car Price Prediction

There is more preprocessing for tabular datasets. To stabilize the training,
we remove 10 percent extreme instances (where the prices is either too high or
too low). Than we normalize the continuous features and encode the categorical
features with one-hot encoding. The regression model used to regress the
features of the cars to their prices is is a 3 layer feedforward neural
network and ReLU activations in each layer. The model is trained on
standardized features with stochastic gradient descent using Adam optimizer
for 100 epochs with batch size 100.

The Encoder, Adversarial Regressor and Decoder of our design for Car Price
Prediction are also feedforward neural network and ReLU activations in each
layer. The number of latent dimension is 2.

As for the hyper-parameter tuning of the Discriminator, we follow the
instructions from the original CcGAN paper, and set the $\sigma=0.035$ and
$k=0.004$.

###  -E Experiment Details of House Sales

The preprocessing of House Sales dataset is similar to Car Price Prediction.
The regression model used to regress the features of the cars to their prices
is is a 3 layer feedforward neural network and ReLU activations in each layer.
The model is trained on standardized features with stochastic gradient descent
using Adam optimizer for 100 epochs with batch size 100.

The Encoder, Adversarial Regressor and Decoder of our design for Car Price
Prediction are also feedforward neural network and ReLU activations in each
layer. The number of latent dimension is 2.


# Pilot Hi Survey of Planck Galactic Cold Clumps with FAST

Ningyu Tang11affiliation: National Astronomical Observatories, CAS, Beijing
100012, China<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>22affiliation: CAS Key Laboratory of FAST, National Astronomical
Observatories, Chinese Academy of Sciences, Beijing 100101, China , Di
Li11affiliation: National Astronomical Observatories, CAS, Beijing 100012,
China<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>22affiliation: CAS
Key Laboratory of FAST, National Astronomical Observatories, Chinese Academy
of Sciences, Beijing 100101, China 33affiliation: School of Astronomy and
Space Science, University of Chinese Academy of Sciences, Beijing 100049,
China , Pei Zuo44affiliation: Kavli Institute for Astronomy and Astrophysics,
Peking University, 100871 Beijing, PR China , Lei Qian11affiliation: National
Astronomical Observatories, CAS, Beijing 100012, China<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>22affiliation: CAS Key Laboratory of FAST,
National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100101, China , Tie Liu55affiliation: Shanghai Astronomical Observatory,
Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, People’s
Republic of China 66affiliation: East Asian Observatory, 660 North A’ohoku
Place, Hilo, HI 96720, USA , Yuefang Wu77affiliation: Department of Astronomy,
School of Physics, Peking University, 100871 Beijing, PR China , Marko
Kr$\rm\check{c}$o11affiliation: National Astronomical Observatories, CAS,
Beijing 100012, China<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>22affiliation: CAS Key Laboratory of FAST, National Astronomical
Observatories, Chinese Academy of Sciences, Beijing 100101, China , Mengting
Liu11affiliation: National Astronomical Observatories, CAS, Beijing 100012,
China<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>22affiliation: CAS
Key Laboratory of FAST, National Astronomical Observatories, Chinese Academy
of Sciences, Beijing 100101, China , Youling Yue11affiliation: National
Astronomical Observatories, CAS, Beijing 100012, China<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>22affiliation: CAS Key Laboratory of FAST,
National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100101, China , Yan Zhu11affiliation: National Astronomical Observatories,
CAS, Beijing 100012, China<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>22affiliation: CAS Key Laboratory of FAST, National
Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China
, Hongfei Liu11affiliation: National Astronomical Observatories, CAS, Beijing
100012, China<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>22affiliation: CAS Key Laboratory of FAST, National Astronomical
Observatories, Chinese Academy of Sciences, Beijing 100101, China , Dongjun
Yu11affiliation: National Astronomical Observatories, CAS, Beijing 100012,
China<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>22affiliation: CAS
Key Laboratory of FAST, National Astronomical Observatories, Chinese Academy
of Sciences, Beijing 100101, China , Jinghai Sun11affiliation: National
Astronomical Observatories, CAS, Beijing 100012, China<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>22affiliation: CAS Key Laboratory of FAST,
National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100101, China , Peng Jiang11affiliation: National Astronomical Observatories,
CAS, Beijing 100012, China<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>22affiliation: CAS Key Laboratory of FAST, National
Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China
, Gaofeng Pan11affiliation: National Astronomical Observatories, CAS, Beijing
100012, China<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>22affiliation: CAS Key Laboratory of FAST, National Astronomical
Observatories, Chinese Academy of Sciences, Beijing 100101, China , Hui
Li11affiliation: National Astronomical Observatories, CAS, Beijing 100012,
China<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>22affiliation: CAS
Key Laboratory of FAST, National Astronomical Observatories, Chinese Academy
of Sciences, Beijing 100101, China , Hengqian Gan11affiliation: National
Astronomical Observatories, CAS, Beijing 100012, China<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>22affiliation: CAS Key Laboratory of FAST,
National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100101, China , Rui Yao11affiliation: National Astronomical Observatories,
CAS, Beijing 100012, China<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>22affiliation: CAS Key Laboratory of FAST, National
Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China
, Shu Liu11affiliation: National Astronomical Observatories, CAS, Beijing
100012, China<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>22affiliation: CAS Key Laboratory of FAST, National Astronomical
Observatories, Chinese Academy of Sciences, Beijing 100101, China , and FAST
Collaboration

###### Abstract

We present a pilot Hi survey of 17 Planck Galactic Cold Clumps (PGCCs) with
the Five-hundred-meter Aperture Spherical radio Telescope (FAST). Hi Narrow
Self-Absorption (HINSA) is an effective method to detect cold Hi being mixed
with molecular hydrogen H2 and improves our understanding of the atomic to
molecular transition in the interstellar medium. HINSA was found in 58% PGCCs
that we observed. The column density of HINSA was found to have an
intermediate correlation with that of 13CO, following $\rm
log(N(HINSA))=(0.52\pm 0.26)log(N_{{}^{13}CO})+(10\pm 4.1)$. Hi abundance
relative to total hydrogen [Hi ]/[H] has an average value of $4.4\times
10^{-3}$, which is about 2.8 times of the average value of previous HINSA
surveys toward molecular clouds. For clouds with total column density
N$\rm{}_{H}>5\times 10^{20}$ cm-2, an inverse correlation between HINSA
abundance and total hydrogen column density is found, confirming the depletion
of cold Hi gas during molecular gas formation in more massive clouds.
Nonthermal line width of 13CO is about 0-0.5 km s-1 larger than that of HINSA.
One possible explanation of narrower nonthermal width of HINSA is that HINSA
region is smaller than that of 13CO. Based on an analytic model of H2
formation and H2 dissociation by cosmic ray, we found the cloud ages to be
within 106.7-107.0 yr for five sources.

###### Subject headings:

ISM: clouds — ISM: evolution — ISM: molecules.

##  1\. Introduction

Cold molecular cores/clumps is the cradle of star formation. Statistical
properties of prestellar cores/clumps are beneficial for understanding general
evolution progress of star formation. The Planck Catalogue of Galactic cold
clumps (PGCCs) released by Planck satellite (Planck Collaboration et al. 2016)
contains 13188 sources with dust temperature ranging from 5.8 to 20 K. These
cold clumps, identified based on the 5’ beam of Planck, contain low-mass cores
to large molecular clouds, providing us a chance of investigating star
formation in various scales, from 0.1 to tens of pc.

Spectral observations toward PGCCs provide extra kinematic and chemical
properties compared to multi-band continuum observations. Molecular line
surveys have been taken in CO (1-0) (Wu et al. 2012; Liu et al. 2012, 2016,
2018), C2H (N=1-0) and N2H+ (J=1-0) (Liu et al. 2019), HCO+ and HCN (J=1-0)
(Yuan et al. 2016). These molecular transitions are collisionally excited with
critical density $n_{crit}\gtrsim$ 103 cm-3, making them tracers of dense gas.
With multiple tracers, the age of a cloud can be modeled through chemical
reaction networks to reproduce the observed line ratio, e.g., [C2H]/[N2H+]
(Pan et al. 2017; Liu et al. 2019). Such estimate, however, is hardly
quantitative calculation and model dependent, due to the complexities of
reaction networks and incomplete knowledge of key parameters, such as
densities.

Hi could be a tracer of cloud age too. When a cloud cools and evolves, the
abundance of Hi in the cloud decreases with time. Cosmic ray dissociation
maintains a small content of Hi even in the central, dense region of molecular
clouds, where Hi is shielded from being dissociated by UV photons and cools
through collision with H2 (Solomon & Werner 1971; Goldsmith & Li 2005). The
cold Hi appears as absorption against its bright Hi background and correlates
with molecular emission. It is called Hi Narrow Self-Absorption (HINSA), and
has been widely detected toward dark clouds in the Taurus/Perseus region and
nearby molecular clouds (Li & Goldsmith 2003) and molecular clouds (Krčo &
Goldsmith 2010). HINSA is a special case of Hi Self-Absorption (HISA), which
is common and has been introduced as an effective tool to study Galactic Hi
distribution (e.g., Gibson et al. 2000; Dickey et al. 2003; Kavars et al.
2003). HISA usually arises in cold neutral medium (CNM) from different
mechanisms, e.g., temperature fluctuations, thus may not match with molecular
emission. For the case of HINSA, it is closely associated with molecular
clouds in terms of sky position, velocity, and non-thermal line width. With
HINSA analysis, a “ring" of enhanced Hi abundance was clearly seen toward dark
cloud B227 (Zuo et al. 2018).

With an instantaneous effective area of 300 m diameter, the Five-hundred-meter
Aperture Spherical radio Telescope (FAST, Nan et al. 2011) located in
Southern-west of China is the most sensitive radio telescope in Hi band. FAST
began its commissioning after Sep 25, 2016. Its commissioning progress was
reviewed in Jiang et al. (2019). The 19 beam receiver was installed on FAST in
mid May of 2018. Spectral backend for the 19 beam receiver was installed in
late July, 2018.

In this paper, we introduced pilot spectral observations with 19 beam receiver
of FAST. The paper is structured as follows: Observation information is
described in Section 2. Analysis of deriving physical parameters is shown in
Section 3. The results and discussion are presented in Section 4. We summarize
our conclusion in Section 5.

##  2\. Observations and Data Archive

###  2.1. Observations of Hi

We selected our sources from the catalogue of Wu et al. (2012) and its
followup survey projects, in which more than 674 PGCCs were observed in J=1-0
transition of CO and its isotopes. Two criteria were applied: (1) position of
the source is within FAST sky coverage but is out of sky coverage of Arecibo
305m telescope; (2) strong CO emission is detected toward the source. Position
distribution of observed 17 PGCCs is shown in Fig. 1.

The FAST observations were conducted between 2018 December 16 and 2019 January
2. We observed the Hi 1420.4058 MHz transition with the central beam of the
19-beam receiver on FAST, which corresponds to a beam width (Full Width at
Half Maximum; FWHM) of 2.9 arcmin. The pointing error is $\sim$ 0.15 arcmin at
1420 MHz (Jiang et al. 2019).

Single pointing observations of PGCCs were taken with the narrow mode of
spectral backend. This mode has 65536 channels in 31.25 MHz bandwidth,
resulting in velocity resolution of 0.1 km s-1 at 1.420 GHz. Each source was
integrated for 10 minutes with a sampling rate of 1.0066 s. For intensity
calibration, noise signal with amplitude of $\sim$ 11 K was injected under a
period of 2.0132 s, which is synchronized with sampling rate.

After intensity calibration, spectral baseline was removed by applying a
linear fitting of line-free range within [-300,300] km s-1. During
observations, system temperature was around 18 K. Reduced root mean square
(rms) of the final spectrum is $\sim$35 mK, which is consistent with
theoretical sensitivity calculation. FAST’s current setup does not include a
local oscillator (LO), and therefore neither doppler tracking nor frequency
switching are available.

In order to test the validity of our data reduction, the PGCC source
G089.29+03.99 located within the sky coverage of Arecibo telescope was
observed. The observed spectrum with FAST is expressed in antenna temperature
($T\rm_{A}$) while the data of GALFA-HI (Peek et al. 2011, 2018) is expressed
in brightness temperature ($T\rm_{B}$). We define a ratio $R$ to convert
measured $T\rm_{A}$ to $T\rm_{B}$, $R=$$T\rm_{B}/$$T\rm_{A}$. A value of
$R=1.15$ is used for comparing observation results of G089.29+03.99. As shown
in Fig. 2, we have consistent spectrum for FAST and Arecibo.

At Hi band, telescope efficiency of FAST is almost a constant value of
$\sim$0.6 when zenith angle (ZA) is less than 26.4∘ and decreases linearly
when ZA is greater than 26.4∘ (Jiang et al. 2019). This implies a $R$ value of
1.67 within $\rm ZA<26.4^{\circ}$. This value is larger than required
conversion value of 1.15 between FAST and Arecibo toward G089.29+03.99. Three
main effects may contribute to this difference.

  1. 1.

The GALFA-HI map was not sampled at the center position precisely of the PGCC
source G098.29+03.99, while we pointed directly at the source.

  2. 2.

The flux for the GALFA-HI map was bootstrapped using the Leiden/Argentine/Bonn
(LAB) survey (Kalberla et al. 2005). This means that while on a large scale
the flux of the map is very accurate, there may be small-scale fluctuations.

  3. 3.

The GALFA map was partially corrected for stray radiation, while we did not.

We applied $R=1.67$ and corrected its dependence of ZA to derive $T\rm_{B}$
from $T\rm_{A}$.

Figure 1.— Spatial distribution of observed 17 PGCCs. Sources with and without
HINSA detection are labeled with cyan circles and red triangles. Sky coverage
of FAST within zenith angle (ZA) of 26.4 and 40 degree are shown with yellow
and grey color. The Galactic plane and Galactic center (GC) are shown with a
blue solid line and a red circle, respectively.  Figure 2.—  Hi spectra for
PGCC source G089.03-41.2. Blue and red color represents data from GALFA-HI and
FAST, respectively. The amplitude of FAST data has been multiplied by a factor
of 1.15.

###  2.2. CO and Dust Archive Data

The 12CO(1-0) and 13CO(1-0) data were obtained under single pointing mode for
investigating dynamic evolution of PGCCs with PMO 13.7m and TRAO-14m
telescopes (Wu et al. 2012; Liu et al. 2016). Angular resolution of both PMO
13.7m and TRAO-14m is $\sim$1 arcmin (FWHM) at 115 GHz. The backend of PMO
13.7m provides bandwidth of 1 GHz and frequency resolution of 61 kHz (0.16 km
s-1 for 12CO and 0.17 km s-1 for 13CO). While the backend of TRAO-14m provides
bandwidth of 125 MHz and frequency resolution of 15 kHz. This frequency
resolution leads to velocity resolution of 0.04 km s-1 at 115.271 GHz. The
spectral resolution of TRAO-14m were smoothed to that of PMO 13.7m.

The dust temperature and dust reddening toward each source were derived from
Planck Collaboration et al. (2016), in which the data were smoothed into FWHM
of 5 arcmin.

##  3\. Analysis

We start analysis by identifying molecular components in velocity based on
multiple gaussian fitting of molecular spectrum. There are two factors: 12CO
profile may be affected by self-absorption due to its large optical depth,
making it difficult to distinguish component numbers; C18O has a relatively
low abundance compared to 12CO. Its emission is weak and may be ‘missed’ under
our observational sensitivity. Considering these factors, we identify velocity
components with 13CO spectral data, whose emission is relatively optically
thinner than 12CO and strong enough to be detected. A total of 23 components
are identified toward 17 PGCCs. Each component is fitted with a Gaussian
function.

Next, we check if there exists a corresponding Hi absorption feature toward
each 13CO emission component. Hi fluctuation is common in Hi spectra. It may
arise from either kinematic effect of two Hi clouds in different velocity or
absorption of cold Hi cloud against brighter Hi background. A feature of Hi
fluctuation is identified as HINSA only when there is corresponding 12CO or
13CO emission. This method of identifying HINSA through single pointing
spectrum could not absolutely exclude the coincidence between Hi fluctuation
from kinematic effect and 13CO emission. As discussed in Gibson et al. (2000),
geometry inspection of Hi mapping would improve this situation.

For these identified HINSA components, we calculated excitation temperature
and column density of atomic and molecular content. They are shown as follows.

###  3.1. CO Excitation

For all observed PGCCs components, temperature ratio between 12CO and 13CO
emission is around 2. This value is much smaller than Galactic [12C/13C] ratio
($\sim$ 60; Lucas & Liszt 1998) in ISM, leading to a reasonable assumption
that 12CO is optically thick, $\tau_{12}\geq 1$. Assuming a filling factor of
1, the excitation temperature of 12CO equals brightness temperature of 12CO.

The excitation temperature of 13CO, $T\rm_{ex}^{{}^{13}CO}$ is assumed to be
same of $T\rm_{ex}^{{}^{12}CO}$. The optical depth of 13CO, $\tau_{13}$ is
given by

$\tau_{13}=-\rm ln(1-\frac{\textit{T}_{b}^{13}}{\textit{T}_{b}^{12}})\ ,$ (1)

under the assumption of $\tau_{12}\geq 1$.

The column density of 13CO is derived through the following equation (Li 2002;
Qian et al. 2012),

$N_{\rm tot}^{\rm{{}^{13}CO}}=3.70\times 10^{14}\int\frac{T_{\rm
b}^{13}}{K}\frac{d\upsilon}{\rm km\ s^{-1}}f_{u}f_{\rm\tau_{13}}f_{\rm
b}\frac{1}{f_{\rm beam}}\ \rm cm^{-2},$ (2)

in which $T_{b}^{13}$ is brightness temperature of 13, $f_{\tau_{13}}$ is the
optical depth of 13CO,
$f_{\tau_{13}}=\int\tau_{13}d\upsilon/\int(1-e^{-\tau_{13}})d\upsilon$. The
uncertainty of $N_{\rm tot}^{\rm{{}^{13}CO}}$ was calculated by uncertainty
propagation of spectral rms of 12CO and 13CO data.

The above analysis of 13CO are derived under the assumption of local thermal
equilibrium (LTE). Application of non-LTE analysis leads to a consistent
result with that of LTE for PGCC (Liu et al. 2016).

###  3.2. Physical Parameters of HINSA

Li & Goldsmith (2003) proposed a three-compnent model of HINSA, which consists
of warm Hi background and foreground gas with a central cold Hi cloud along
line of sight. In this model, optical depth of background and foreground Hi
gas are $\tau_{b}$ and $\tau_{f}$, respectively. The factor $p$ is defined as
$p$=$\tau_{f}$/($\tau_{f}$+$\tau_{b}$). Kinematic distance of each PGCC
component was derived through Galactic rotation curve (Brand, & Blitz 1993).
We combine kinematic distance and Galactic Hi surface density distribution
(Nakanishi, & Sofue 2003) to estimate $p$ along sightline of the PGCC. This is
different from the formula adopted in Li & Goldsmith (2003), which is
reasonable for cloud with high Galactic latitude but would be invalid for
cloud in the Galactic plane. Details of this calculation could be found in
section 4.2 of Tang et al. (2016). $T\rm_{b}(\upsilon)$ and
$T\rm_{H{\textsc{i}}}(\upsilon)$ represent observed Hi brightness temperature
spectrum and brightness temperature spectrum when central cold Hi cloud is
absent, respectively. Absorption temperature, $T\rm_{ab}(\upsilon)$ =
$T\rm_{H{\textsc{i}}}(\upsilon)-T\rm_{b}(\upsilon)$ could be described with
(Li & Goldsmith 2003) ,

$T{\rm{}_{ab}}(\upsilon)=[pT_{H{\textsc{i}}}(\upsilon)+(T_{c}-T_{x})(1-\tau_{f})](1-e^{-\tau}),$
(3)

where $T\rm_{c}$ is continuum temperature including cosmic microwave
background (CMB, 2.73 K), Galactic synchrotron radiation ($\sim$ 0.7 K), and
continuous emission from other radiation mechanisms. The value of $T\rm_{c}$
is derived by the equation $T\rm_{c}$ = 2.7 +
$T_{c408}$($\nu\rm_{H{\textsc{i}}}$/408 MHz)-2.8, in which $T_{c408}$ is
continuum temperature around 408 MHz (Haslam et al. 1982), and
$\nu\rm_{H{\textsc{i}}}$ is frequency of Hi emission. Derived value of
$T\rm_{c}$ toward our targets ranges from 3.8 to 5.6 K. $T\rm_{x}$ is the
excitation temperature of HINSA component, and depends on kinetic temperature
Tk and effective temperature of radiation field. The gas of observed PGCCs are
in molecular phase, resulting a reasonable assumption that $T\rm_{x}=$
$T\rm_{k}$. The value of $p$ is around 0.9. We follow Li & Goldsmith (2003),
the value of $\tau_{f}$ is assumed as 0.1, which implies a total optical depth
of unity along Galactic sightline. This assumption would not significantly
affect the result. The optical depth of HINSA, $\tau$ is an unknown parameter
to be solved. Assuming a Gaussian broadening, $\tau$ is described in the
formula of
$\tau=\tau_{0}e^{-\frac{(\upsilon-\upsilon_{0})^{2}}{2\sigma^{2}}}$, in which
$\tau_{0}$ is peak optical depth at central velocity.

$T\rm_{HI}$ is recovered Hi spectrum without absorption from the central cold
cloud. There are many methods to derive $T\rm_{HI}$, for instance, averaging
nearby spectra without this absorption cloud (e.g., Gibson et al. 2000) or
fitting observed profile with polynomial (e.g., Li & Goldsmith 2003).
Considering complexity of Hi emission around the Galactic plane, the method of
observing off positions or fitting with polynomial may introduce significant
uncertainty. The method of polynomial fitting depends on choice of polynomial
index.

HINSA and molecular hydrogen are believed to mix well inside the cloud. Total
velocity broadening of each transition,
$\sigma_{tot}=\sqrt{\sigma_{th}^{2}+\sigma_{nt}^{2}}$, in which $\sigma_{th}$
and $\sigma_{nt}$ represent thermal broadening due to collision and non-
thermal broadening including turbulence, magnetic field and other effects.
Krčo et al. (2008) introduced a new method to derive $T\rm_{HI}$ by utilizing
physical information (central velocity, non-thermal linewidth and kinematic
temperature) derived from CO data to constrain fitting of HINSA profile. We
introduce this method for analysis in this paper. For each HINSA component, 4
parameters of HINSA are fitted with the following considerations (Krčo &
Goldsmith 2010):

  1. 1.

Central velocity. A maximum deviation of 10% linewidth from central velocities
of 13CO is allowed.

  2. 2.

Non-thermal linewidth. It is allowed within 50% deviation from fitting 13CO.

  3. 3.

Kinematic temperature, which is derived from 12CO and is kept as a constant
during fitting.

  4. 4.

Optical depth. This is a relative free parameter without priori information.
Initial value of 0.1 with limitation of [0,10] is adopted during fitting.

Since the width of HINSA is narrow, HINSA feature would become apparent in
second derivative of observed spectrum (Kavars et al. 2005). The minimization
of the integrated squared sum of second derivative would be sufficient to
remove HINSA feature (Krčo et al. 2008). We follow Krčo et al. (2008) to
derive $T\rm_{HI}$ by fitting a least minimum of second derivation of the
input spectrum. This method performs perfect for simple spectrum but may lead
to large deviation for complex spectrum. The above constraints of inputting
parameters are introduced to reduce uncertainty during fitting. An example of
HINSA fitting is shown in Fig. 3. Corresponding second derivatives of
$T\rm_{HI}$ and $T\rm_{b}$ profile is shown in Fig. 4. It is clearly seen that
second derivatives at central velocity becomes flat. Column density of HINSA
is calculated by (Li & Goldsmith 2003),

$N({\rm HINSA})=1.95\times 10^{18}\tau_{0}T_{\rm x}\Delta V\ \rm cm^{-2},$ (4)

where $\tau_{0}$, $T\rm_{x}$ and $\Delta V$ are peak optical depth, excitation
temperature, and width of HINSA, respectively. The uncertainty of N(HINSA) is
estimated through uncertainty propagation during fitting $\tau_{0}$,$T\rm_{k}$
and $\Delta V$.

###  3.3. Total Column Density of Hydrogen

We obtained total column density of molecular hydrogen $N\rm_{H_{2}}$ and its
uncertainty from Planck Collaboration et al. (2016). $N\rm_{H_{2}}$ toward
each PGCCs is derived from color extinction E(B-V), which is based on multi-
band continuum emission. Obtained value of $N\rm_{H_{2}}$ consists
contribution from all gas components along line of sight of PGCCs. For PGCCs
with two or more components, we assume $N\rm_{H_{2}}$ of each component is
linearly correlated with its column density of 13CO. The column density of
hydrogen of each component is then calculated with

$N_{H}=N(\rm{HINSA})+2N\rm_{H_{2}}.$ (5)

No estimates of $N\rm_{H_{2}}$ are given toward PGCCs G089.62+02.16 and
G093.16+09.61. $N\rm_{H_{2}}$ of gas components in these two sources are
calculated from N(13CO). Thus we have

$\displaystyle N_{H}$
$\displaystyle=N(\rm{HINSA})+\frac{2N(\rm^{13}CO)}{X{\rm{}_{{}^{13}CO/H_{2}}}}$
(6)
$\displaystyle=N(\rm{HINSA})+\frac{2N(\rm^{13}CO)}{X{\rm{}_{{}^{13}C/^{12}C}}X{\rm{}_{{}^{12}CO/H_{2}}}},$
(7)

where $X{\rm{}_{{}^{13}CO/H_{2}}}$ is 13CO/H2 ratio in a cloud. It contains
two parts, $X{\rm{}_{{}^{13}C/^{12}C}}$ and $X{\rm{}_{{}^{12}CO/H_{2}}}$. They
are 13C/12C and 12CO/H2 ratio, respectively. In previous research, the value
of these two parts are generally adopted as 1/60 and 1$\times 10^{4}$ in
Galactic ISM (e.g., Liu et al. 2016).

We compare validity of calculating $N(\rm H_{2})$ by adopting
$X{\rm{}_{{}^{13}CO/H_{2}}}=1.67\times 10^{-6}$ for observed PGCCs with $N(\rm
H_{2})$ estimation from dust. As shown in Fig. 5, the adoption of
$X{\rm{}_{{}^{13}CO/H_{2}}}=1.67\times 10^{-6}$ would overestimate
$N\rm_{H_{2}}$ by a factor of 6.2. This discrepancy may origin from variation
of 12C and 13C in ISM with different environment. We adopt
$X{\rm{}_{{}^{13}CO/H_{2}}}=(1.05\pm 0.36)\times 10^{-5}$ for calculating
$N\rm_{H}$ in PGCCs G089.62+02.16 and G093.16+09.61.

Figure 3.— A example of HINSA fitting toward G089.75-02.16. $Top$ panel: black
and red solid lines are the observed Hi and the recovered Hi background
spectrum THI, respectively. $Bottom$ panel: The profile of absorption
temperature $T\rm_{ab}$ (blue), 12CO (green) and 13CO spectrum (red).  Figure
4.—  Second derivatives of $T\rm_{HI}$ and $T\rm_{b}$ profile for
G089.75-02.16. The velocity range is [0.5,3.5] km s-1. Dashed line represents
central velocity of the HINSA component.  Figure 5.— The
$N\rm_{H_{2}}$(13CO)/$N_{H_{2}}$(dust) ratio as a function of
$N\rm_{H_{2}}$(dust). $N\rm_{H_{2}}$(13CO) and $N_{H_{2}}$(dust) are estimated
molecular hydrogen from $N(^{13}$CO) with converstion factor of 1.67$\times
10^{-6}$ and dust (Planck collaboration 2016), respectively. Red and green
dashed line represents averaged value and 1 $\sigma$ deviation, respectively.

##  4\. Results and Discussion

A total of 17 PGCCs inside FAST sky coverage were surveyed. We have detected
13 HINSA components toward 10 PGCCs, implying a HINSA detection rate of
$\sim$58% in PGCCs. This detection percentage is smaller than that of 77% in
Taurus/Perseus region (Li & Goldsmith 2003, hereafter LG03) and 80% of 48
molecular clouds (Krčo & Goldsmith 2010, hereafter KG10). This deviation may
arise from limited sample number and different strategy of sample selection.

We unbiasedly selected PGCCs without constraint on visual extinction and
presence of Young Stellar Objects (YSOs), which could prevent cold gas from
external and internal heating. A embedded YSO is found to locate inside PGCCs
G108.84-00.81. Though LG03 stated that the presence of YSOs would not
significantly affect the detection rate of HINSA, the statement should be
checked further by including molecular clouds with far distance. In KG10, they
selected sample clouds with peak visual extinction $A\rm_{V}>$ 6 mag and
without YSOs inside.

The spectra of sources with HINSA detection and without HINSA detection are
shown in Appendix A and Appendix B, respectively. There are many reasons that
may be responsible for non-HINSA detection. For example, the excitation
temperature of the HI gas inside a molecular cloud is comparable to that of
background continuum temperature. As shown in Fig. 1, there is no obvious
spatial difference between HINSA and non-HINSA detection. The sample of 17
PGGCs are small. More samples are needed to further investigate this point.

###  4.1. HINSA Column Density and Abundance

N(HINSA) of PGCCs has an average of 6.3$\times$ 1018 cm-2. As a comparision,
average value of N(HINSA) in LG03 and KG10 are 7.2$\times$ $10^{18}$ and
6.3$\times$ 1018 cm-2. These results are consistent though we focus on
different category of molecular cloud.

HINSA abundance, N(HINSA)/NH represents abundance of cold Hi in a cloud. The
average HINSA abundance in PGCCs is $4.4\times 10^{-3}$, which is about 3
times the value derived in both LG03 and KG10. We notice that $N_{H_{2}}$ is
derived from Planck data with 5’ FWHM, which is larger than that of Hi (2.9’)
and CO (1’) data we used. This would lead to overestimation of N(HINSA)/NH. As
an estimation, the peak column density $N\rm_{H_{2}}^{peak}$ of PGCC would be
1.38 times the value of $N\rm_{H_{2}}$ (Planck Collaboration et al. 2016).
Considering this uncertainty, the lower limit of averaged N(HINSA)/NH would be
$3.2\times 10^{-3}$.

The correlation between log($N(\rm HINSA)$) and log($N(\rm^{13}CO)$) is
investigated in Fig. 6. An extra uncertainty of 20% is considered due to flux
calibration of 13CO. As described in section 3.2, the uncertainty of fitting
complex spectrum is difficult to estimate quantitatively and is not included
in the error bar of N(HINSA) of Fig. 6. The correlation can be expressed with

$\rm log(N(HINSA))=a\ log(N_{{}^{13}CO})+b,$ (8)

in which a= 0.52$\pm$0.26 and b=10$\pm$4.1. The Pearson correlation
coefficient is 0.52, indicating an intermediate correlation between HINSA and
13CO.

The uncertainty of $N\rm_{H_{2}}$ from dust could reach 50% the value itself.
It is larger than the uncertainty percent of N(HINSA) and $N(\rm^{13}CO)$. The
uncertainties of HINSA abundance N(HINSA)/NH and total column density NH are
significant and must be included during correlation analysis. When uncertainty
in both X and Y axes are considered, the correlation between N(HINSA)/NH and
NH are fitted with the following equation,

$\rm log(N(HINSA)/N_{H})=c\ log(N_{H})+d.$ (9)

As shown in Fig. 7, there is a trend that HINSA abundance decreases along
increasing hydrogen column density. The fitted parameter c = -1.1$\pm$ 0.32
and d=21$\pm$ 6.7 when all data points are included. The slope of -1.1$\pm$
0.32 could not exclude the possibility that linear decrease of HINSA abundance
is just due to division of NH. The two components, -3.73 km s-1 of
G093.16+09.61 and 4.54 km s-1 of G107.64-09.31 show significant deviation. The
13CO intensity of these two components are low, which may lead to
underestimation of $N$(13CO), and thus underestimation of NH with the
calculation method we used in this paper. When these two components are
excluded, the fitted parameter c = -2.7$\pm$ 1.0 and d= 54$\pm$ 22\. The slope
of -2.7$\pm$ 1.0 indicates that the decreasing trend is physically meaningful.
HINSA column density will decrease in increasing NH. This result implies that
cold Hi would be depleted in denser molecular cloud. Since the two components
with low 13CO emission would affect this analysis, more HINSA samples with
accurate estimation of NH are needed to check this correlation.

###  4.2. Kinematic Environment

We investigated further about kinematic environment of PGCCs. Thermal and
nonthermal effect (e.g., turbulence) are two main factors that lead to line
broadening. Since PGCCs are cold ($\sim$ 5-25 K), the thermal broadening width
is small, $\sim$ 0.3 km s-1. Unlike thermal effect, nonthermal effect does not
depend on kinetic temperature and particle mass. In order to check the
turbulence effect on different species, We compare nonthermal width of HINSA
and CO in Fig. 8. We found that nonthermal width of 13CO is about 0-0.5 km s-2
larger than that of HINSA though spatial resolution of CO ($\sim$ 1 arcmin) is
smaller than that of Hi ($\sim$ 3 arcmin). The smoothing of CO data to same
resolution of Hi is expect to result in larger linewidth of 13CO, enlarging
disparity between CO and Hi linewidth.

In order to check if this result is due to applying constraints of non-thermal
width while fitting HINSA in Section 3.2, a wider limitation was chosen. We
adopted 100% deviation from fitting 13CO. Variation of fitted line width under
new limitation is less than 0.5% compared to old one. Limitation of line width
is used for deriving local extreme solution and does not affect fitting
result.

Line broadening caused by turbulence ($\Delta V$) is found to connect with
cloud size ($L$) through a relation , $\Delta V\sim L^{q}$. The value of q is
0.38 when first proposed by Larson (1981). Further observations toward giant
molecular clouds indicate a value of 0.5 (e.g., Solomon et al. 1987). This
relation implies that smaller regions are less turbulent. The broader line
width of 13CO than HINSA may indicate that HINSA lie in a smaller, interior
region than 13CO.

###  4.3. Transition between Atomic and Molecular Phase

The existence of cold Hi in clumps is due to balance between H2 formation on
dust grain and H2 dissociation by external environment, e.g., illumination
from cosmic ray. HINSA abundance can be an effective indicator of cloud
evolution (Goldsmith & Li 2005). In order to compare with analytic model, we
have to identify volume density of total hydrogen. For PGCCs with both
distance estimation $d$ (including kinematic distance derived rotation curve)
and cloud size in angular resolution $\Delta\theta$, the cloud size is
estimated with $L=d\Delta\theta$. For other PGCCs, we adopt the fitting
formula derived in Heyer & Brunt (2004), $\sigma=0.9(L/1\rm pc)^{0.56}$ km s-1
to calculate the size of PGCCs. Average density is calculated with
$\bar{n}\rm{{}_{H}}=\frac{N\rm{{}_{H}}}{L}$.

We follow the analytic model in Goldsmith & Li (2005). In a cloud with Hi
density of $n\rm_{H{\textsc{i}}}$ and H2 density of $n\rm_{H_{2}}$, its Hi
density fraction
$x_{1}=n\rm_{H{\textsc{i}}}/$$(n\rm_{H{\textsc{i}}}+$$2n\rm_{H_{2}})$ as a
function of time $t$ is given by

$x_{1}(t)=1-\frac{2k^{\prime}n_{0}}{2k^{\prime}n_{0}+\zeta_{H_{2}}}[1-\textrm{exp}(\frac{-t}{\tau\rm_{HI\rightarrow
H_{2}}})],$ (10)

in which $k^{\prime}$ is the formation rate coefficient of H2.
$\zeta\rm_{H_{2}}$ is the ionization rate of H2 by cosmic-ray. The value of
$5.2\times 10^{-17}$ s-1 has been adopted (van der Tak, & van Dishoeck 2000).
Time constant between atomic and molecular phase $\tau\rm_{HI\rightarrow
H_{2}}=$$\frac{1}{2k^{\prime}n_{0}+\zeta\rm_{H_{2}}}$. By adopting
$x_{1}=$$N(\rm HINSA)/N\rm_{H}$ and $n_{0}=\bar{n}\rm{{}_{H}}$, cloud age
could be calculated from Equation 10.

Solution can be derived for only 5 of 13 PGCC components. Fig. 9 shows the
evolution age as a function of abundance of Hi and total proton density. Cloud
age of these 5 clouds lies in the range $10^{6.7}$ yr $\leqslant$
$\tau\rm_{cloud}$ $\leqslant$ $10^{7.0}$ yr. Besides, they locate in a
transition phase between unstable and steady state.

Cloud has density profile distribution with cloud radius, e.g., power law
density profile (Shu 1977), which results in a denser central region. The beam
size of Hi ($\sim$ 2.9 arcmin) and CO ($\sim$ 1 arcmin) observations allows us
to only cover the central dense region of a PGCC clump. The estimation of
cloud density by simply dividing cloud size with column density would
underestimate cloud density, leading to overestimate cloud age. As discussed
in section 4.2, we adopt peak proton density, which is 1.38 times the current
value to estimate a lower limit of cloud age. Lower limits of cloud age for 5
PGCC components is in the range $10^{6.6}$ yr $\leqslant$ $\tau\rm_{cloud}$
$\leqslant$ $10^{6.9}$ yr. These value are about 23% smaller than that shown
in Fig. 9.

We are not able to find available solution for the rest 8 components due to
relatively low Hi abundance under its density. A possible reason is that cloud
size was overestimated, resulting underestimation of total hydrogen density.

Figure 6.— Correlation between column density of HINSA and 13CO. The red
dashed line represents linear fitting of the data.  Figure 7.— The relation
between HINSA abundance $N\rm_{HINSA}$/$N\rm_{H}$ compared to total hydrogen
$N\rm_{H}$. The red and green dashed lines represent fitting for the right 11
data points and all 13 data points.  Figure 8.— Nonthermal width of HINSA and
13CO.  Figure 9.—  Time dependence of Hi abundance for various proton density
(in cm-3). The dashed lines represent different logarithmic densities with
index interval of 0.25. Five sources are plotted with black diamonds. The
cloud age of rest 8 sources is calculated with measured Hi abundance and the
assumption of steady state. They are shown with grey triangles.

##  5\. Summary

In this paper, we have observed Hi spectra toward 17 PGCCs with FAST. Thirteen
HINSA components are found toward 10 sources. The abundance of cold Hi in
these clouds varies from $5.1\times 10^{-4}$ to $1.3\times 10^{-2}$ with an
average value of $4.4\times 10^{-3}$.

We found that HINSA column density in PGCCs has a intermediate correlation
with 13CO column density at logarithmic scale. This correlation follows a
relationship of $\rm log(N(HINSA))=(0.52\pm 0.26)log(N_{{}^{13}CO})+(10\pm
4.1)$.

When all date points are included, we could not exclude the possibility that
decreasing HINSA abundance (N(HINSA)/NH) is due to division of NH. If we only
consider HINSA components with $N_{H}$> $5\times 10^{20}$ cm-2, decreasing
HINSA abundance as a function of $N\rm_{H}$ is physically meaningful. This
confirms depletion of cold Hi for forming molecular gas in more massive PGCCs.

The nonthermal linewidths of 11 of 13 HINSA components are narrower than that
of 13CO. This result may implicate that HINSA located in inner and smaller
region than that of 13CO.

We compare our results with analytic model of cloud evolution. Five of 13
HINSA components satisfy cloud evolution at age range of 106.7\- 107.0 yr.

Due to a limited sample of this pilot survey with FAST, there still exist some
uncertainties in the above results. A larger survey of HINSA features are
necessary for improving these results.

We are grateful to the anonymous referee for his/her constructive suggestions,
which have greatly improved this paper. This work is supported by the National
Natural Science Foundation of China grant No. 11988101, No. 11725313 and No.
11803051, National Key R&D Program of China No. 2017YFA0402600, CAS
International Partnership Program No.114A11KYSB20160008, National Science
Foundation of China (11721303) and the National Key R&D Program of China
(2016YFA0400702), CAS "Light of West China" Program, CAS Interdisciplinary
Innovation Team and Young Researcher Grant of National Astronomical
Observatories, Chinese Academy of Sciences. LQ is supported in part by the
Youth Innovation Promotion Association of CAS (id. 2018075).

## References

  * Brand, & Blitz (1993) Brand, J., & Blitz, L. 1993, A&A, 275, 67
  * Dickey et al. (2003) Dickey, J. M., McClure-Griffiths, N. M., Gaensler, B. M., & Green, A. J. 2003, ApJ, 585, 801
  * 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
  * Haslam et al. (1982) Haslam, C. G. T., Salter, C. J., Stoffel, H., et al. 1982, A&AS, 47, 1
  * Heyer & Brunt (2004) Heyer, M. H., & Brunt, C. M. 2004, ApJ, 615, L45
  * Jiang et al. (2019) Jiang, P., Yue, Y., Gan, H., et al. 2019, Science China Physics, Mechanics, and Astronomy, 62, 959502
  * Kalberla et al. (2005) Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, A&A, 440, 775
  * Kavars et al. (2003) Kavars, D. W., Dickey, J. M., McClure-Griffiths, N. M., Gaensler, B. M., & Green, A. J. 2003, ApJ, 598, 1048
  * Kavars et al. (2005) Kavars, D. W., Dickey, J. M., McClure-Griffiths, N. M., et al. 2005, ApJ, 626, 887
  * 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
  * Larson (1981) Larson, R. B. 1981, MNRAS, 194, 809
  * Li (2002) Li, D. 2002, Ph.D. Thesis
  * Li & Goldsmith (2003) Li, D., & Goldsmith, P. F. 2003, ApJ, 585, 823
  * Liu et al. (2018) Liu, T., Kim, K.-T., Juvela, M., et al. 2018, ApJS, 234, 28
  * Liu et al. (2012) Liu, T., Wu, Y., & Zhang, H. 2012, ApJS, 202, 4
  * Liu et al. (2016) Liu, T., Zhang, Q., Kim, K.-T., et al. 2016, ApJS, 222, 7
  * Liu et al. (2019) Liu, X.-C., Wu, Y., Zhang, C., et al. 2019, A&A, 622, A32
  * Lucas & Liszt (1998) Lucas, R., & Liszt, H. 1998, A&A, 337, 246
  * Nakanishi, & Sofue (2003) Nakanishi, H., & Sofue, Y. 2003, PASJ, 55, 191
  * Nan et al. (2011) Nan, R., Li, D., Jin, C., et al. 2011, International Journal of Modern Physics D, 20, 989
  * Pan et al. (2017) Pan, Z., Li, D., Chang, Q., et al. 2017, ApJ, 836, 194
  * Peek et al. (2018) Peek, J. E. G., Babler, B. L., Zheng, Y., et al. 2018, ApJS, 234, 2
  * Peek et al. (2011) Peek, J. E. G., Heiles, C., Douglas, K. A., et al. 2011, ApJS, 194, 20
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A28
  * Qian et al. (2012) Qian, L., Li, D., & Goldsmith, P. F. 2012, ApJ, 760, 147
  * Shu (1977) Shu, F. H. 1977, ApJ, 214, 488
  * Solomon et al. (1987) Solomon, P. M., Rivolo, A. R., Barrett, J., & Yahil, A. 1987, ApJ, 319, 730
  * Solomon & Werner (1971) Solomon, P. M., & Werner, M. W. 1971, ApJ, 165, 41
  * Tang et al. (2016) Tang, N., Li, D., Heiles, C., et al. 2016, A&A, 593, A42
  * van der Tak, & van Dishoeck (2000) van der Tak, F. F. S., & van Dishoeck, E. F. 2000, A&A, 358, L79
  * Wu et al. (2012) Wu, Y., Liu, T., Meng, F., et al. 2012, ApJ, 756, 76
  * Yuan et al. (2016) Yuan, J., Wu, Y., Liu, T., et al. 2016, ApJ, 820, 37
  * Zuo et al. (2018) Zuo, P., Li, D., Peek, J. E. G., et al. 2018, ApJ, 867, 13

Table 1Derived parameters of observed sources.  Source | RA(J2000) | Dec(J2000) | Comp |  12CO |  13CO | HINSA | N(H)
---|---|---|---|---|---|---|---
$T\rm_{ex}$ |  | Vcen | $\Delta V$ | $N\rm_{{}^{13}CO}$ |  | Vcen | $\Delta V$ | $N\rm_{HINSA}$
(K) |  | (km s-1) | (km s-1) | (cm-2) |  | (km s-1) | (km s-1) | (cm-2) | (cm-2)
G089.29+03.99 | 20:50:43.20 | 50:25:38.7 | 1 | 15.95$\pm$0.31 |  | -4.41$\pm$0.02 | 2.85$\pm$0.04 | 2.18$\pm$0.05 |  | -4.70$\pm$0.02 | 2.48$\pm$0.12 | 4.0$\pm$0.5 | 20.15$\pm$6.32
G089.62+02.16 | 21:00:44.18 | 49:30:29.8 | 1 | 8.47$\pm$0.87 |  | -1.39$\pm$0.44 | 1.85$\pm$0.17 | 0.28$\pm$0.04 |  | -1.21$\pm$0.02 | 1.30$\pm$0.01 | 3.31$\pm$0.07 | 2.63$\pm$0.97
|  |  | 2 | 9.08$\pm$0.86 |  | -0.16$\pm$0.07 | 1.27$\pm$0.05 | 0.65$\pm$0.07 |  | -0.29$\pm$0.02 | 1.35$\pm$0.003 | 16.31$\pm$0.09 | 6.18$\pm$2.25
G089.64-06.84 | 21:37:58.08 | 43:11:7.8 | 1 | 18.07$\pm$0.27 |  | 12.86$\pm$0.01 | 1.33$\pm$0.02 | 1.56$\pm$0.04 |  | 12.77$\pm$0.01 | 0.91$\pm$0.01 | 3.90$\pm$0.14 | 13.01$\pm$7.79
G089.75-02.16 | 21:20:11.99 | 46:38:41.3 | 1 | 20.99$\pm$0.39 |  | 1.86$\pm$0.02 | 1.73$\pm$0.03 | 1.88$\pm$0.05 |  | 1.79$\pm$0.01 | 1.48$\pm$0.01 | 13.88$\pm$0.14 | 13.61$\pm$5.73
G090.76-04.57 | 21:34:00.72 | 45:36:504 | 1 | 17.81$\pm$0.24 |  | 10.13$\pm$0.02 | 1.82$\pm$0.05 | 0.86$\pm$0.03 |  | 10.21$\pm$0.01 | 1.50$\pm$0.02 | 7.21$\pm$0.30 | 5.43$\pm$2.54
G092.79+09.12 | 20:37:48.71 | 56:19:27.1 | 1 | 11.8$\pm$0.5 |  | -2.53$\pm$0.03 | 1.63$\pm$0.07 | 0.81$\pm$0.05 |  | -2.37$\pm$0.01 | 1.43$\pm$0.01 | 14.86$\pm$0.13 | 14.04$\pm$2.93
G093.16+09.61 | 20:36:31.44 | 56:54:48.6 | 1 | 8.43$\pm$0.49 |  | -3.73$\pm$0.10 | 1.40$\pm$0.23 | 0.140$\pm$0.030 |  | -3.61$\pm$0.01 | 1.00$\pm$0.03 | 0.68$\pm$0.04 | 1.34$\pm$0.55
|  |  | 2 | 9.19$\pm$0.49 |  | -1.75$\pm$0.02 | 1.19$\pm$0.04 | 1.28$\pm$0.25 |  | -1.68$\pm$0.01 | 1.21$\pm$0.01 | 4.81$\pm$0.06 | 12.23$\pm$4.88
G107.64-09.31 | 23:14:26.40 | 50:39:44.6 | 1 | 10.9$\pm$0.4 |  | 4.54$\pm$0.22 | 2.20$\pm$0.50 | 0.115$\pm$0.034 |  | 4.72$\pm$0.02 | 1.27$\pm$0.04 | 1.86$\pm$0.13 | 1.16$\pm$0.51
|  |  | 2 | 16.3$\pm$0.4 |  | 6.98$\pm$0.02 | 1.26$\pm$0.03 | 1.19$\pm$0.05 |  | 6.85$\pm$0.02 | 0.79$\pm$0.03 | 1.22$\pm$0.10 | 12.04$\pm$4.25
G116.08-02.40 | 23:56:45.11 | 59:44:20.0 | 1 | 16.52$\pm$0.32 |  | -0.96$\pm$0.01 | 1.67$\pm$0.03 | 1.84$\pm$0.06 |  | -1.00$\pm$0.01 | 1.54$\pm$0.02 | 7.14$\pm$0.23 | 11.22$\pm$5.02
G122.01-07.47 | 00:45:00.02 | 55:23:18.6 | 1 | 17.17$\pm$0.22 |  | -52.18$\pm$0.02 | 1.23$\pm$0.04 | 0.59$\pm$0.03 |  | -52.21$\pm$0.02 | 0.87$\pm$0.04 | 2.7$\pm$0.3 | 6.42$\pm$3.10
G089.36-00.69 | 21:12:24.00 | 47:23:32.3 | 1 | 7.43$\pm$0.28 |  | -24.01$\pm$0.07 | 1.85$\pm$0.17 | 0.155$\pm$0.020 |  | - | - | - | 2.97$\pm$ 0.84
|  |  | 2 | 17.77$\pm$0.27 |  | 2.29$\pm$0.02 | 2.71$\pm$0.03 | 3.19$\pm$0.05 |  | - | - | - | 61.23$\pm$ 17.26
|  |  | 3 | 8.89$\pm$0.28 |  | 4.94$\pm$0.41 | 1.93$\pm$0.56 | 0.052$\pm$0.020 |  | - | - | - | 1.00$\pm$ 0.28
G090.76-04.57 | 21:34:00.72 | 45:36:504 | 2 | 7.31$\pm$0.25 |  | 12.64$\pm$0.06 | 0.47$\pm$0.14 | 0.028$\pm$0.011 |  | - | - | - | 0.18$\pm$0.11
G093.22-04.60 | 21:44:44.88 | 47:13:3.4 | 1 | 17.16$\pm$0.25 |  | 3.60$\pm$0.01 | 1.52$\pm$0.02 | 2.04$\pm$0.04 |  | - | - | - | 32.9$\pm$ 15.8
G093.31-11.68 | 22:09:44.88 | 41:42:30.2 | 1 | 15.32$\pm$0.39 |  | -0.35$\pm$0.02 | 0.83$\pm$0.05 | 0.305$\pm$0.028 |  | - | - | - | 5.47$\pm$1.43
G093.36-12.07 | 22:11:11.52 | 41:25:35.8 | 1 | 8.52$\pm$0.37 |  | 0.24$\pm$0.16 | 0.86$\pm$0.37 | 0.036$\pm$0.020 |  | - | - | - | 4.09$\pm$ 1.84
G101.14-15.29 | 22:54:59.75 | 42:35:51.7 | 1 | 15.45$\pm$0.43 |  | -9.59$\pm$0.01 | 0.97$\pm$0.03 | 0.98$\pm$0.05 |  | - | - | - | 18.82$\pm$ 6.39
G108.84-00.81 | 22:58:48.72 | 58:56:54.6 | 1 | 22.74$\pm$0.45 |  | -49.67$\pm$0.02 | 3.27$\pm$0.04 | 4.73$\pm$0.10 |  | - | - | - | 86.8$\pm$ 64.3
G121.93-07.67 | 00:44:28.34 | 55:11:24.4 | 1 | 17.01$\pm$0.21 |  | -52.07$\pm$0.01 | 1.18$\pm$0.01 | 1.202$\pm$0.024 |  | - | - | - | 35.8$\pm$ 13.2

##  Appendix A Detected HINSA toward PGCCs

Figure 10.— Hi and CO spectra toward other 9 PGCCs with HINSA . Descriptions
are same with that of Fig.3.

##  Appendix B Non-HINSA spectra toward PGCCs

Figure 11.— Hi and CO spectra toward 7 PGCCs without HINSA. The observed Hi
spectrum and CO spectra are shown in the top and bottom panel, respectively.