We study Individual Fairness (IF) for Bayesian neural networks (BNNs).
Specifically, we consider the $\epsilon$-$\delta$-individual fairness notion, which requires that, for any pair of input points that are $\epsilon$-similar according to a given similarity metrics, the output of the BNN is within a given tolerance $\delta>0.$
We leverage bounds on statistical sampling over the input space and the relationship between adversarial robustness and individual fairness to derive a framework for the systematic estimation of $\epsilon$-$\delta$-IF, designing Fair-FGSM and  Fair-PGD as global, fairness-aware extensions to gradient-based attacks for BNNs.
We empirically study IF of a variety of approximately inferred BNNs with different architectures on fairness benchmarks, and compare against deterministic models learnt using frequentist techniques. Interestingly, we find that BNNs trained by means of approximate Bayesian inference consistently tend to be markedly more individually fair than their deterministic counterparts.

§ INTRODUCTION

Deep learning models have achieved state-of-the-art performance in a wide variety of tasks [Goodfellow et al., 2016].

Several calls for caution have, however, recently been raised about their deployment in tasks where fairness is of concern [Barocas and Selbst, 2016]. In fact, Neural Networks (NNs) have been found to reinforce negative biases from sensitive datasets [Bolukbasi et al., 2016], discriminating against individuals on the basis of attributes such as gender or race.
To address this, research efforts have been directed at both measuring the fairness of NNs, their de-biased training, as well as defining precise notions of fairness.

Given a BNN $f$ and a similarity metric between individuals $d_{fair}$, which encodes a task-dependent notion of similarity [Ilvento, 2020],  Individual Fairness (IF) enforces that all pairs of similar individuals in the input space get treated similarly by $f$ [Dwork et al., 2012].  As opposed to the statistical nature of group fairness [Mehrabi et al., 2021],  IF aims at computing worst-case bias measures over a model input space.
Consequently, albeit being defined on the full input space and over a fairness similarity metric, because of its worst-case nature, IF has been linked to adversarial robustness [Yeom and Fredrikson, 2021].
As it has recently been show that Bayesian Neural Networks (BNNs) have a tendency to be less fragile to adversarial attacks than their frequentist counter-parts [Carbone et al., 2020], it is natural to wonder whether approximate Bayesian inference may also have a positive impact over the IF of a neural network.
However, to the best of our knowledge, no work has been conducted along these lines of inquire.

In this paper, we investigate the IF of BNNs and empirically evaluate it on various benchmarks.
While exact computations of IF in BNNs is infeasible due to their non-convexity, we exploit the relationship between IF and adversarial robustness [Yeom and Fredrikson, 2021, Benussi et al., 2022]  to develop a framework for the adaptation of adversarial attack methods for IF.
In particular, we explicitly instantiate Fair-FGSM and Fair-PGD as extensions of their corresponding adversarial attacks [Goodfellow et al., 2015, Athalye et al., 2018] by employing gradient steps modifications and projections specific to $d_{fair}$ metrics commonly used in the fairness literature.
Furthermore, while attack methods estimate IF locally  around a given individual, we use concentration inequalities [Boucheron et al., 2004] to statistically bound the worst-case expected IF over all pairs of individuals.

We perform an empirical evaluation of IF in BNNs on a selection of benchmarks, including the Adult dataset [Dua and Graff, 2017] and Folktables [Ding et al., 2021], and on a variety of different architectural, approximate Bayesian inference, and similarity metric parameters.
We compare the results obtained by BNNs with those obtained by deterministic NNs and deep ensembles.
We find that BNNs consistently outperform their deterministic counterparts in terms of individual fairness. That is, albeit still learning certain biases, all things being equal,  BNNs show a tendency to be fairer than their deterministic counter-parts.
We finish the paper with an empirical analysis  inquiring this interesting property and a discussion on why approximate Bayesian inference may lead to fairer prediction models.

§ INDIVIDUAL FAIRNESS FOR BNNS

We consider a NN $f^w : X \rightarrow \mathbb{R}^m$ trained on a dataset  $\mathcal{D} = \{ (x_i,y_i)\}_{i=1}^{N} $ composed by $N$ points, where $w$ is the aggregated vector of weights and biases, $X \subset \mathbb{R}^n$ is the input space, and $n$ and $m$ are respectively the input and output dimensions.
For simplicity, we focus on the case of binary classification, however, our results can be trivially extended to the regression and multi-class case.
In a Bayesian setting, $w$ is sampled from a random variable $\mathbf{w}$ whose posterior distribution is obtained by application the Bayes rule:
$ p(w|\mathcal{D}) \propto p(\mathcal{D}|w)p(w) $,
where  $p(w)$ is the prior and $p(\mathcal{D}|w)$ is the likelihood.
The posterior predictive distribution $\pi$ over an input $x$ is then obtained  by averaging the posterior distribution over a sigmoid function $\sigma$, i.e.,
$\pi(x) = \int \sigma(f^w(x))p(w|\mathcal{D})dw .$

To define individual fairness for BNNs, we adopt the $\epsilon$-$\delta$-IF [John et al., 2020], which depends on a similarity metric $d_{fair}: X  \times X \rightarrow \mathbb{R}_{\geq 0 }  $  encapsulating task-dependent information about sensitive attributes.

Given $\epsilon \geq 0$ and $\delta \geq 0$,  we say that a BNN $f^{\mathbf{w}}$ is $\epsilon$-$\delta$-individually fair iff $
\forall x', x'' \in X s.t.\ d_{fair}(x',x'') \leq \epsilon  \implies \vert \pi(x') - \pi(x'') \vert \leq \delta,$
that is, if the predictive distribution of $\epsilon$-similar individuals are $\delta$-close to each other.

Notice that while related to notion of adversarial robustness for BNNs [Cardelli et al., 2019, Wicker et al., 2020], Definition <ref> has some key differences.
First, rather than being a local definition specific to a test point and its neighbourhood, individual fairness looks at the worst-case behaviour over the whole input space, i.e., simultaneously looking at neighbourhoods around all the input points in $X$.
Furthermore, while adversarial attacks are generally developed for an $\ell_p$ metric, individual fairness is built around a task-specific similarity metric $d_{fair}$.
Intuitively, $d_{fair}$ needs to encode the desired notion of similarity between individuals [Dwork et al., 2012].

While a number of metrics have been discussed in the literature [Ilvento, 2020], we focus on the following ones, which are widely used and
can be automatically learnt from data [Yurochkin et al., 2020]:[While not directly investigated in this paper, our fairness attacks can be generalised to similarity metrics built over embeddings [Ruoss et al., 2020], by attacking the embeddings as well.]

* Weighted $\mathbf{\ell_p}$. In this case $d_{\text{fair}}(x',x'')$ is defined as a weighted version of an $\ell_p$ metric, i.e.$d_{\text{fair}}(x',x'') = \sqrt[p]{\sum_{i=1}^n \theta_i |x'_i - x''_i|^p}  $, where the  weights $\theta_i$ can be set accordingly to their correlation with the sensitive attribute.

* Mahalanobis distance. In this case we have $d_{\text{fair}}(x',x'')  = \sqrt{(x'-x'')^T S^{-1} (x'-x'')}$, for a given positive semi-definite (SPD) matrix $S$. Intuitively, $S$ accounts for the intra-correlation of features to capture latent dependencies w.r.t. the sensitive features.

§ FAIRNESS ATTACKS

Definition <ref> could be reformulated according to the following optimization problem

\begin{align}\label{eq:optimisation}
\delta^* = \max_{x' \in X} {\max_{\substack{x'' \in X\\ d_{fair}(x',x'') \leq \epsilon  }} \vert \pi(x') - \pi(x'') \vert}.
\end{align}

Checking IF is then equivalent to check whether $\delta^* \leq \delta $ or not.
In Eqn (<ref>), the inner maximization problem finds the point in a specific $d_{fair}$-neighbourhood that maximizes the variation in the predictive posterior, while the outer one considers the global part of the definition.
We proceed by solving the inner-optimisation problem with gradient-based techniques, and rely on statistical methods to solve with high confidence the outer problem.

Inner Problem
Because of the non-convexity of NN architectures, exact computation of the inner maximization problem in Eqn (<ref>) is generally infeasible.
We proceed by adapting gradient-based optimization techniques
commonly used to detect adversarial attacks to our setting. In particular, in the case of the Fast Gradient Sign Method (FGSM) [Goodfellow et al., 2015], given the log-loss $L$ associated to the likelihood function (e.g., the binary entropy for the sigmoid), and a point $x$ with associated label $y$, we obtain that:

$$ x_{x',\max} \approx x' +\eta \cdot  \text{sgn}( \mathbb{E}_{w \sim p(w|\mathcal{D})}[\nabla_x L(  f^w(x),y  )  ]).$$

$\eta \in \mathbb{R}^n$ scales the attack strength according to  $\epsilon$ and $d_{fair}$. It is the first ingredient that needs to be adapted for $\epsilon$-$\delta$-IF.
In the case of the weighted $\ell_p$ metric, $\eta$ can be simply defined as the vector with entries $\eta_i = \epsilon/\sqrt{\theta_i}$.
For the Mahalanobis distance, given $x'$ and $\epsilon$, $d_{fair}(x',x'') \leq \epsilon$ describes a hyper-ellipsoid centred in $x'$.
We can hence set $\eta$ to be the vector whose generic component is scaled according to the ellipsoid axis, i.e., $\eta_i = \epsilon \sqrt{S_{ii}}$.
Furthermore, while one is guaranteed to remain inside the $d_{fair}$-neighbourhood in the $\ell_p$ case, for the Mahalanobis metric one needs to check whether the attack is inside the $\epsilon$-$d_{fair}$ ball or not. In the latter case, we set up a local line search problem to project $x_{x',\max}$ back to a point inside the ball.
Other gradient-based attacks such as PGD [Madry et al., 2018] or CW [Carlini and Wagner, 2017] can be similarly adapted by taking consecutive gradients. Similarly, other more general attacks could be employed [Yuan et al., 2020, Däubener et al., 2020] as well as training techniques [Liu et al., 2019, Wicker et al., 2021]. In Section <ref>, we will consider both FGSM and PGD, which we denote with Fair-FGSM and Fair-PGD to distinguish them from the adversarial attack counter-parts.

Outer Problem
The final step needed for the computation $\epsilon$-$\delta$-IF is to check whether:

\begin{align}\label{eq:verification}
\max_{x' \in X}\big( \pi(x_{x',\max})- \pi(x') \big)  \leq \delta.
\end{align}

Again, exact computation of Eqn (<ref>) is infeasible as it would require solving a non-convex optimisation problem over a possibly large input space. We relax the problem and instead compute the probability that a point sapled from $X$ is fair. In order to do that we iid sample $n_{\text{samples}}$ points from $X$ and build the following empirical estimator:
$\hat{p}=\sum_{i=1}^{n_{\text{samples}}} \frac{\phi(x_i)}{n_{\text{samples}}}, \quad \text{where }  \phi(x_i)=\begin{cases} 1 &\text{if }\pi(x_{x_i,\max}) - \pi(x_i)  \leq \delta \\
0 & \text{otherwise}
\end{cases}.$
Note that $\phi(x_i)$ is a binary variable. Consequently, $\hat{p}$ approximate the probability that a point sampled from $X$ satisfies the inner condition in Eqn (<ref>). Lemma <ref> below is a straightforward application of the Chernoff bound [Hellman and Raviv, 1970] and guarantees that if $n_{\text{samples}}$ is large enough, then the approximation error introduced by the empirical estimator can be made arbitrarily small with high confidence.

Assume that $n_{\text{samples}}$ points are sampled iid from $X$ according to a (possibly unknown) probability measure $P.$
Call $ p=\mathbb{E}_{x'\sim P}[ \phi(x')]. $
Then for any $\gamma,\theta >0$ such that  $ n_{samples} > \frac{1}{2\theta^2} \log \left( \frac{2}{\gamma} \right),$ it holds that
$ P( |\hat{p} - p| > \theta\big) \leq \gamma. $

Therefore, in order to approximate the solution of the optimisation problem in Equation <ref> with statistical confidences $\theta$ and $\gamma$  we compute the smallest $\delta$ such that $\hat{p}=1$ with $n_{samples}$ points randomly sampled from the input space.

The overall resulting methodology is summarised in the pseudocode in Algorithm <ref>. Need to add some explanation of the algorithm if you decide to keep it. Personal opinion is that it does not add much on what we already say, so I would remove it and eventually put it in the Appendix

$\epsilon$-$\delta$-IF estimation for BNNs Might remove
Input: $X$ – input space, $\pi$ – posterior predictive distribution, $d_{fair}$ – similarity metric, $\epsilon$, $\delta$ – IF parameters, $\theta, \gamma$ – Chernoff bound parameters

Output: $\hat{p}$ – IF estimation
[1]

$n_{samples} \gets$ number of samples by Chernoff bounds (<ref>)
$n, k \gets$ 0, 0
$ n <  \lceil n_{samples} \rceil$
$x' \gets$ sample uniformly from $X$
$x'_{bias} \gets x'$; $\eta \gets$ build step size for $\epsilon$-$d_{fair}$ balls; $l \gets 0$
$ l <  steps$
$x'_{bias} \gets x'_{bias} + \frac{\eta}{steps} \mathbb{E}_{w \sim p(w|\mathcal{D})}[\nabla_x L(  f^w(x'_{bias}),y'  ) $
Project $x'_{bias}$ back to $\{ x | d_{fair}(x',x) \leq \epsilon \}$
$l \gets l+1$
$\mathsf{SAT}\gets$ check $|\pi(x') - \pi(x'_{bias})|  \leq \delta$
$k \gets k + \mathsf{SAT}$; $n \gets n+1$
$\hat{p} \gets$ $k/n$

return $\hat{p}$

§ RESULTS

Heatmaps for the estimations of individual fairness ($\delta$) across various architectures, inference methods (dNN, VI or HMC), datasets, and similarity metrics. The first two rows give results for the Adult dataset (respectively Fair-FGSM with $\ell_p$, and Fair-PGD with Mahalanobis). The last two for the Folktables dataset.

We employ Fair-FGSM and Fair-PGD to evaluate and compare the IF of deterministic Neural Networks (dNNs) and Bayesian Neural Networks (BNNs).
We approximate the BNNs' posterior either with VI (in particular using variational online Gauss-Newton [Khan et al., 2018]) or with Hamiltonian Monte Carlo (HMC) [Neal, 2012].
The NNs are trained to solve the taks associated to the Adult [Dua and Graff, 2017] and the Folktables [Ding et al., 2021] datasets, where we take gender as the sensitive.[Code to run the experiments can be found at <https://github.com/alicedoherty/bayesian-fairness>.]

Parametric Analysis
We evaluate how IF is affected by architectural parameters.
Namely, across the two similarity metrics, we begin by fixing $\epsilon=0.1$ and analysing how $\delta$ is affected by the numbers of hidden layers and of neurons per layer, in ranges typically used for these datasets [Yurochkin et al., 2020, John et al., 2020, Yeom and Fredrikson, 2021, Benussi et al., 2022]. The results for these analyses are plotted in the heatmaps in Figure <ref>.
The three columns of heatmaps respectively show the IF estimations for models learned with dNN, VI, and HMC.
The first two rows give the results for the Adult dataset, respectively for Fair-FGSM with the weighted $\ell_p$ metric and Fair-PGD with the Mahalanobis similarity metric. The last two rows show the analogous results on the Folks dataset.

We immediately observe how, in the overwhelming majority of the cases, the models learnt by HMC are markedly fairer (i.e., they have lower values of $\delta$) than their deterministic counter-parts.
The same generally applies also to VI models, except for the IF obtained on the Adult dataset by means of Fair-PGD for the Mahalanobis metric (second row, second column).
This combined observations almost perfectly mirrors the behaviour of approximate Bayesian inference in adversarial robustness settings [Carbone et al., 2020], albeit here over the whole input space, and around $d_{fair}$-neighbourhoods.

Interestingly, both adversarial robustness and overfitting generally worsen as the model becomes deeper. Instead, we observe that more often than not this pattern is reversed in IF. We posit that while the NNs are indeed becoming more fragile to attacks as they get deeper, it is also less reliant on each specific input feature, rather it builds on a non-trivial representation of them – therefore having a tendency to perform better on IF tasks. On the other hand, as expected, when the number of neurons per layer increases, we get less fair models.

In the left plot of Figure <ref> we analyse how individual fairness is affected by $\epsilon$, i.e., the maximum dissimilarity allowed between two individuals. We give a selection of the results for an architecture with 3 layers and 16 neurons (thick line) or 32 neurons (dashed line) trained on the Adult dataset. Of course, we observe that as $\epsilon$ increases $\delta$ increases as well, as a greater $\epsilon$ implies that we are looking at the models' behaviour on increasingly more dissimilar pairs of individuals.
Interestingly, we find that BNNs' tendency to be more fair is consistent across the different values of epsilon.

Analysis of Posterior Predictive
We experimentally inquire on the reason behind BNNs' improved individual fairness.
To this end, for a given BNN posterior predictive distribution, we solve the inner problem of Section <ref> over 100 randomly sampled test points.
Rather than doing this over the full posterior, we randomly sample 15 times a number $k$ of weights realisation from the posterior, and compute the predictive posterior over them.
We investigate values of $k$ ranging from $1$ (i.e., a single deterministic-like network extracted from the BNN) to $50$.

We also compare BNNsr, with an averaging of deterministic neural networks in the form of a Deep Ensemble (DE) [Lakshminarayanan et al., 2017].
The results for this analysis are given in the right plot of Figure <ref>.
Surprisingly, even when using only one posterior sample from the BNN distribution, the models learnt with Bayesian approximation method are already significantly more fair than with the frequentist DE technique.
Furthermore, especially in the case of HMC, the more samples are drawn from the posterior distribution at prediction time, the fairer the model becomes – with a $50\%$ bias reduction when going from $1$ to $50$.
On the other hand, notice that averaging  does not help with DE, whose fairness estimation remains almost constant throughout the plot.

Left: Fairness for different maximum distances $\epsilon$ for deterministic and Bayesian NNs with 3 layers and either 16 (thick lines) or 32 (dashed lines) neurons per layer. Right: Changes in empirical distribution of fairness across deterministic and Bayesian predictors computed over an increasing number of posterior samples.

§ CONCLUSION

In this paper we extended adversarial attack techniques to study individual fairness (IF) for BNNs. On a set of experiments we empirically showed that BNNs tend to be intrinsically fairer than their deterministic counter-parts.  Furthermore, we showed how the particular approximate inference method has an impact on IF, with more approximate methods being less fair compared to approaches such as HMC.
Finally, we have empirically shown how increased fairness is likely due to a combination of Bayesian training and Bayesian averaging, which may have a beneficial effect in reducing the magnitude of the gradients.

[Athalye et al., 2018]
Anish Athalye, Nicholas Carlini, and David Wagner.
Obfuscated gradients give a false sense of security: Circumventing
defenses to adversarial examples.
In International conference on machine learning, pages
274–283. PMLR, 2018.

[Barocas and Selbst, 2016]
Solon Barocas and Andrew D Selbst.
Big data's disparate impact.
Calif. L. Rev., 104:0 671, 2016.

[Benussi et al., 2022]
Elias Benussi, Andrea Patane, Matthew Wicker, Luca Laurenti, and Marta
Individual fairness guarantees for neural networks.
International Joint Conference on Artificial Intelligence,

[Bolukbasi et al., 2016]
Tolga Bolukbasi, Kai-Wei Chang, James Y Zou, Venkatesh Saligrama, and Adam T
Man is to computer programmer as woman is to homemaker? debiasing
word embeddings.
Advances in neural information processing systems, 29, 2016.

[Boucheron et al., 2004]
Stéphane Boucheron, Gábor Lugosi, and Olivier Bousquet.
Concentration inequalities.
Advanced Lectures on Machine Learning: ML Summer Schools 2003,
Canberra, Australia, February 2-14, 2003, Tübingen, Germany, August 4-16,
2003, Revised Lectures, pages 208–240, 2004.

[Carbone et al., 2020]
Ginevra Carbone, Matthew Wicker, Luca Laurenti, Andrea Patane, Luca Bortolussi,
and Guido Sanguinetti.
Robustness of bayesian neural networks to gradient-based attacks.
Advances in Neural Information Processing Systems,
33:0 15602–15613, 2020.

[Cardelli et al., 2019]
Luca Cardelli, Marta Kwiatkowska, Luca Laurenti, Nicola Paoletti, Andrea
Patane, and Matthew Wicker.
Statistical guarantees for the robustness of bayesian neural
International Joint Conference on Artificial Intelligence,

[Carlini and Wagner, 2017]
Nicholas Carlini and David Wagner.
Magnet and" efficient defenses against adversarial attacks" are not
robust to adversarial examples.
arXiv preprint arXiv:1711.08478, 2017.

[Däubener et al., 2020]
Sina Däubener, Joel Frank, Thorsten Holz, and Asja Fischer.
Efficient calculation of adversarial examples for bayesian neural
In Third Symposium on Advances in Approximate Bayesian
Inference, 2020.

[Ding et al., 2021]
Frances Ding, Moritz Hardt, John Miller, and Ludwig Schmidt.
Retiring adult: New datasets for fair machine learning.
Advances in Neural Information Processing Systems,
34:0 6478–6490, 2021.

[Dua and Graff, 2017]
Dheeru Dua and Casey Graff.
UCI machine learning repository, 2017.
URL <http://archive.ics.uci.edu/ml>.

[Dwork et al., 2012]
Cynthia Dwork, Moritz Hardt, Toniann Pitassi, Omer Reingold, and Richard Zemel.
Fairness through awareness.
In Proceedings of the 3rd innovations in theoretical computer
science conference, pages 214–226, 2012.

[Goodfellow et al., 2016]
Ian Goodfellow, Yoshua Bengio, Aaron Courville, and Yoshua Bengio.
Deep learning, volume 1.
MIT Press, 2016.

[Goodfellow et al., 2015]
Ian J Goodfellow, Jonathon Shlens, and Christian Szegedy.
Explaining and harnessing adversarial examples.
International Conference on Learning Representations, 2015.

[Hellman and Raviv, 1970]
Martin Hellman and Josef Raviv.
Probability of error, equivocation, and the chernoff bound.
IEEE Transactions on Information Theory, 160
(4):0 368–372, 1970.

[Ilvento, 2020]
Christina Ilvento.
Metric learning for individual fairness.
In 1st Symposium on Foundations of Responsible Computing (FORC
2020). Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2020.

[John et al., 2020]
Philips George John, Deepak Vijaykeerthy, and Diptikalyan Saha.
Verifying individual fairness in machine learning models.
In Conference on Uncertainty in Artificial Intelligence, pages
749–758. PMLR, 2020.

[Khan et al., 2018]
Mohammad Khan, Didrik Nielsen, Voot Tangkaratt, Wu Lin, Yarin Gal, and Akash
Fast and scalable bayesian deep learning by weight-perturbation in
In International conference on machine learning, pages
2611–2620. PMLR, 2018.

[Lakshminarayanan et al., 2017]
Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell.
Simple and scalable predictive uncertainty estimation using deep
Advances in neural information processing systems, 30, 2017.

[Liu et al., 2019]
Xuanqing Liu, Yao Li, Chongruo Wu, and Cho-Jui Hsieh.
Adv-bnn: Improved adversarial defense through robust bayesian neural
International Conference on Learning Representations, 2019.

[Madry et al., 2018]
Aleksander Madry, Aleksandar Makelov, Ludwig Schmidt, Dimitris Tsipras, and
Adrian Vladu.
Towards deep learning models resistant to adversarial attacks.
International Conference on Learning Representations, 2018.

[Mehrabi et al., 2021]
Ninareh Mehrabi, Fred Morstatter, Nripsuta Saxena, Kristina Lerman, and Aram
A survey on bias and fairness in machine learning.
ACM Computing Surveys (CSUR), 540 (6):0 1–35,

[Neal, 2012]
Radford M Neal.
Bayesian learning for neural networks, volume 118.
Springer Science & Business Media, 2012.

[Ruoss et al., 2020]
Anian Ruoss, Mislav Balunovic, Marc Fischer, and Martin Vechev.
Learning certified individually fair representations.
Advances in Neural Information Processing Systems,
33:0 7584–7596, 2020.

[Wicker et al., 2020]
Matthew Wicker, Luca Laurenti, Andrea Patane, and Marta Kwiatkowska.
Probabilistic safety for bayesian neural networks.
In Conference on Uncertainty in Artificial Intelligence, pages
1198–1207. PMLR, 2020.

[Wicker et al., 2021]
Matthew Wicker, Luca Laurenti, Andrea Patane, Zhuotong Chen, Zheng Zhang, and
Marta Kwiatkowska.
Bayesian inference with certifiable adversarial robustness.
In International Conference on Artificial Intelligence and
Statistics, pages 2431–2439. PMLR, 2021.

[Yeom and Fredrikson, 2021]
Samuel Yeom and Matt Fredrikson.
Individual fairness revisited: transferring techniques from
adversarial robustness.
In Proceedings of the Twenty-Ninth International Conference on
International Joint Conferences on Artificial Intelligence, pages 437–443,

[Yuan et al., 2020]
Matthew Yuan, Matthew Wicker, and Luca Laurenti.
Gradient-free adversarial attacks for bayesian neural networks.
Third Symposium on Advances in Approximate Bayesian Inference,

[Yurochkin et al., 2020]
Mikhail Yurochkin, Amanda Bower, and Yuekai Sun.
Training individually fair ml models with sensitive subspace
International Conference on Learning Representations, 2020.


* Bow [1997] Numerous review articles appear in Parity and Time Reversal Violation in Compound States and Related Topics (edited by N. Auerbach and J. D. Bowman, World Scientific, 1997) https://www.worldscientific.com/doi/pdf/10.1142/3227 .
  * Bowman and Gudkov [2014] J. D. Bowman and V. Gudkov, Phys. Rev. C 90, 065503 (2014).
  * Ryndin [1968] R. M. Ryndin, $3^{\text{rd}}$ Proceedings of the LNPI Winter School (in Russian) (Leningrad, USSR, 1968).
  * Bilenky  _et al._ [1964] S. M. Bilenky, L. I. Lapidus, and R. M. Ryndin, Usp. Fiz. Nauk 84, 243 (1964).
  * Bilenky  _et al._ [1969] S. M. Bilenky, L. I. Lapidus, and R. M. Ryndin, Phys. Usp. 11, 512 (1969).
  * Kabir [1988] P. K. Kabir, Phys. Rev. D 37, 1856 (1988).
  * Lamoreaux and Golub [1994] S. K. Lamoreaux and R. Golub, Phys. Rev. D 50, 5632 (1994).
  * Masuda  _et al._ [1992] Y. Masuda _et al._ , Proceedings of the WEIN ‘92 (edited by T. D. Vylov World Scientific, Singapore, 1992).
  * Skoy [1996] V. R. Skoy, Phys. Rev. D 53, 4070 (1996).
  * Lukashevich  _et al._ [2011] V. V. Lukashevich, A. V. Aldushchenkov, and D. Dallman, Phys. Rev. C 83, 035501 (2011).
  * Masuda [1998] Y. Masuda, Nuclear Physics A 629, 479 (1998), quark Lepton Nuclear Physics.
  * Masuda [2000] Y. Masuda, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 440, 632 (2000).
  * Kabir [1989] P. Kabir, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 284, 63 (1989).
  * Baryshevsky and Podgoretsky [1965] V. Baryshevsky and M. Podgoretsky, Sov. Phys JETP 20, 704 (1965).
  * Abragam and Goldman [1982] A. Abragam and M. Goldman, Nuclear magnetism: Order and Disorder, International Series of Monographs on Physics (Oxford University Press, Oxford, UK, 1982).
  * Masuda [1996] Y. Masuda, Parity and Time Reversal Violation in Compound Nuclear States and Related Topics (edited by N. Auerbach and J. D. Bowman, World Scientific, Singapore, 1996) p. 83.
  * Serebrov [1996] A. P. Serebrov, Parity and Time Reversal Violation in Compound Nuclear States (edited by N. Auerbach and J. D. Bowman, World Scientific, Singapore, 1996) p. 327.
  * Flambaum [1993] V. V. Flambaum, Physica Scripta T46, 198 (1993).
  * Flambaum and Vorov [1993] V. V. Flambaum and O. K. Vorov, Phys. Rev. Lett. 70, 4051 (1993).
  * Flambaum and Izrailev [1997a] V. V. Flambaum and F. M. Izrailev, Phys. Rev. E 56, 5144 (1997a).
  * Flambaum  _et al._ [2015] V. V. Flambaum, M. G. Kozlov, and G. F. Gribakin, Phys. Rev. A 91, 052704 (2015).
  * Flambaum and Vorov [1994] V. V. Flambaum and O. K. Vorov, Phys. Rev. C 49, 1827 (1994).
  * Flambaum and Vorov [1995a] V. V. Flambaum and O. K. Vorov, Phys. Rev. C 51, 1521 (1995a).
  * Flambaum and Vorov [1995b] V. V. Flambaum and O. K. Vorov, Phys. Rev. C 51, 2914 (1995b).
  * Flambaum and Gribakin [2000] V. Flambaum and G. Gribakin, Philosophical Magazine Part B 80, 2143 (2000).
  * Flambaum  _et al._ [1994] V. V. Flambaum, A. A. Gribakina, G. F. Gribakin, and M. G. Kozlov, Phys. Rev. A 50, 267 (1994).
  * Dzuba  _et al._ [2017] V. A. Dzuba, J. C. Berengut, C. Harabati, and V. V. Flambaum, Phys. Rev. A 95, 012503 (2017).
  * Flambaum  _et al._ [1998] V. Flambaum, A. Gribakina, and G. Gribakin, Physical Review A 58, 230 (1998).
  * Flambaum  _et al._ [2002] V. V. Flambaum, A. A. Gribakina, G. F. Gribakin, and C. Harabati, Phys. Rev. A 66, 012713 (2002).
  * Harabati  _et al._ [2017] C. Harabati, J. C. Berengut, V. V. Flambaum, and V. A. Dzuba, Journal of Physics B: Atomic, Molecular and Optical Physics 50, 125004 (2017).
  * Frankle  _et al._ [1992] C. M. Frankle, J. D. Bowman, J. E. Bush, P. P. J. Delheij, C. R. Gould, D. G. Haase, J. N. Knudson, G. E. Mitchell, S. Penttilä, H. Postma, N. R. Roberson, S. J. Seestrom, J. J. Szymanski, S. H. Yoo, V. W. Yuan, and X. Zhu, Phys. Rev. C 46,
778 (1992).

  * Stephenson  _et al._ [1998] S. L. Stephenson, J. D. Bowman, B. E. Crawford, P. P. J. Delheij, C. M. Frankle, M. Iinuma, J. N. Knudson, L. Y. Lowie, A. Masaike, Y. Matsuda, G. E. Mitchell, S. I. Penttilä, H. Postma, N. R. Roberson, S. J. Seestrom, E. I. Sharapov, Y.-F. Yen, and V. W. Yuan, Phys. Rev. C 58, 1236 (1998).
  * Desplanques and Noguera [1996] B. Desplanques and S. Noguera, Nuclear Physics, A 598, 139 (1996).
  * Blons  _et al._ [1984] J. Blons, C. Mazur, D. Paya, M. Ribrag, and H. Weigmann, Nucl. Phys. A 414, 1 (1984).
  * Blons [1989] J. Blons, Nucl. Phys. A 502, 121c (1989).
  * Flambaum and Sushkov [1980] V. Flambaum and O. Sushkov, Physics Letters B 94, 277 (1980).
  * Sushkov and Flambaum [1981a] O. Sushkov and V. Flambaum, Yad. Fiz. 33(3) 629, 1981. Sov. J. Nucl. Phys. 33, 329 (1981a).
  * Sushkov and Flambaum [1981b] O. Sushkov and V. Flambaum, Yad. Fiz. 33(1) 59, 1981. Sov. J. Nucl. Phys. 33, 31 (1981b).
  * Flambaum and Zelevinsky [1995] V. Flambaum and V. Zelevinsky, Physics Letters B 350, 8–12 (1995).
  * Auerbach  _et al._ [1995] N. Auerbach, J. D. Bowman, and V. Spevak, Phys. Rev. Lett. 74, 2638 (1995).
  * Spevak and Auerbach [1995] V. Spevak and N. Auerbach, Physics Letters B 359, 254 (1995).
  * Butler [2016] P. A. Butler, Journal of Physics G: Nuclear and Particle Physics 43, 073002 (2016).
  * Cao  _et al._ [2020] Y. Cao, S. E. Agbemava, A. V. Afanasjev, W. Nazarewicz, and E. Olsen, Phys. Rev. C 102, 024311 (2020).
  * Butler [2020] P. Butler, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 476, 20200202 (2020).
  * Nazarewicz [1994] W. Nazarewicz, Nuclear Physics A 574, 27 (1994).
  * Flambaum and Feldmeier [2020] V. V. Flambaum and H. Feldmeier, Phys. Rev. C 101, 015502 (2020).
  * [67] The NUDAT program for nuclear data, https://www.nndc.bnl.gov/nudat2/.
  * Flambaum and Dzuba [2020] V. V. Flambaum and V. A. Dzuba, Phys. Rev. A 101, 042504 (2020).
  * Lane and Pendlebury [1960] A. Lane and E. Pendlebury, Nuclear Physics 15, 39 (1960).
  * Flambaum [1992] V. V. Flambaum, Phys. Rev. C 45, 437 (1992).
  * Hussein  _et al._ [1995] M. S. Hussein, A. K. Kerman, and C.-Y. Lin, Z. Phys. A 351, 301 (1995).
  * Feshbach  _et al._ [1996] H. Feshbach, M. S. Hussein, and A. K. Kerman, Parity and Time Reversal Violation in Compound Nuclear States and Related Topics (edited by N. Auerbach and J. D. Bowman, World Scientific, Singapore, 1996) p. 157.
  * Auerbach [1992] N. Auerbach, Phys. Rev. C 45, R514 (1992).
  * Auerbach and Bowman [1992] N. Auerbach and J. D. Bowman, Phys. Rev. C 46, 2582 (1992).
  * Auerbach and Spevak [1994] N. Auerbach and V. Spevak, Phys. Rev. C 50, 1456 (1994).
  * Bowman  _et al._ [1992] J. D. Bowman, G. T. Garvey, C. R. Gould, A. C. Hayes, and M. B. Johnson, Phys. Rev. Lett. 68, 780 (1992).
  * Kerman  _et al._ [1963] A. K. Kerman, L. S. Rodberg, and J. E. Young, Phys. Rev. Lett. 11, 422 (1963).
  * Hussein  _et al._ [2016] M. S. Hussein, B. V. Carlson, and A. K. Kerman, Acta Physica Polonica. Series B 47, 391 (2016).
  * Gudkov and Shimizu [2018] V. Gudkov and H. M. Shimizu, Phys. Rev. C 97 (2018).
  * Flambaum and Gribakin [1994] V. V. Flambaum and G. F. Gribakin, Phys. Rev. C 50, 3122 (1994).
  * Bowman  _et al._ [1990] J. D. Bowman, C. D. Bowman, J. E. Bush, P. P. J. Delheij, C. M. Frankle, C. R. Gould, D. G. Haase, J. Knudson, G. E. Mitchell, S. Penttila, H. Postma, N. R. Roberson, S. J. Seestrom, J. J. Szymanski, V. W. Yuan, and X. Zhu, Phys. Rev. Lett. 65, 1192 (1990).
  * Bohr and Mottelson [1998] A. Bohr and B. R. Mottelson,  _Nuclear Structure_ (World Scientific Publishing Company, 1998).
  * Flambaum  _et al._ [1984a] V. Flambaum, I. Khriplovich, and O. Sushkov,
Phys. Lett. B 146, 367 (1984a).

  * Desplanques  _et al._ [1980] B. Desplanques, J. F. Donoghue, and B. R. Holstein, Annals of Physics 124, 449 (1980).
  * Flambaum and Murray [1997] V. V. Flambaum and D. W. Murray, Phys. Rev. C 56, 1641 (1997).
  * Noguera and Desplanques [1986] S. Noguera and B. Desplanques, Nuclear Physics A 457, 189 (1986).
  * Dubovik and Zenkin [1986] V. Dubovik and S. Zenkin, Annals of Physics 172, 100 (1986).
  * Feldman  _et al._ [1991] G. B. Feldman, G. A. Crawford, J. Dubach, and B. R. Holstein, Phys. Rev. C 43, 863 (1991).
  * Wasem [2012] J. Wasem, Phys. Rev. C 85, 022501 (2012).
  * Blyth  _et al._ [2018] D. Blyth _et al._ (NPDGamma Collaboration), Phys. Rev. Lett. 121, 242002 (2018).
  * Haxton and Holstein [2013] W. C. Haxton and B. R. Holstein, Progress in Particle and Nuclear Physics 71, 185 (2013), fundamental Symmetries in the Era of the LHC.
  * Flambaum  _et al._ [1984b] V. V. Flambaum, I. B. Khriplovich, and O. P. Sushkov,  _On the possibility to study P- and T-odd nuclear forces in atomic and molecular experiments_, Tech. Rep. (Akad. Nauk Novosibirsk. Inst. Yarn. Fiz., Novosibirsk, 1984).
  * Flambaum  _et al._ [1986] V. Flambaum, I. Khriplovich, and O. Sushkov, Nuclear Physics A 449, 750 (1986).
  * Flambaum  _et al._ [2014] V. V. Flambaum, D. DeMille, and M. G. Kozlov, Phys. Rev. Lett. 113, 103003 (2014).
  * Haxton and Henley [1983] W. C. Haxton and E. M. Henley, Phys. Rev. Lett. 51, 1937 (1983).
  * Khriplovich and Korkin [2000] I. Khriplovich and R. Korkin, Nuclear Physics A 665, 365 (2000).
  * Dmitriev and Sen’kov [2003] V. F. Dmitriev and R. A. Sen’kov, Physics of Atomic Nuclei 66, 1940–1945 (2003).
  * Yamanaka  _et al._ [2017] N. Yamanaka, B. Sahoo, N. Yoshinaga, T. Sato, K. Asahi, and B. Das, The European Physical Journal A 53, 1 (2017).
  * de Vries  _et al._ [2015] J. de Vries, E. Mereghetti, and A. Walker-Loud, Phys. Rev. C 92, 045201 (2015).
  * Pospelov and Ritz [2005] M. Pospelov and A. Ritz, Annals of Physics 318, 119 (2005), special Issue.
  * Palos [2018] L. B. Palos, (NOPTREX Collaboration), talk at the Thirteenth Conference on the Intersections of Particle and Nuclear Physics (unpublished) (2018).
  * Moody and Wilczek [1984] J. E. Moody and F. Wilczek, Phys. Rev. D 30, 130 (1984).
  * Marsh [2016] D. J. E. Marsh, Physics Reports 643, 1 (2016), arXiv:1510.07633 [astro-ph.CO] .
  * Graham  _et al._ [2015] P. W. Graham, D. E. Kaplan, and S. Rajendran, Phys. Rev. Lett. 115, 221801 (2015).
  * Gupta  _et al._ [2016] R. S. Gupta, Z. Komargodski, G. Perez, _et al._ , J. High Energ. Phys. 02 , 166 (2016).
  * Flacke  _et al._ [2017] T. Flacke, C. Frugiuele, E. Fuchs, R. S. Gupta, and G. Perez, Journal of High Energy Physics 6, 1 (2017).
  * Stadnik  _et al._ [2018] Y. V. Stadnik, V. A. Dzuba, and V. V. Flambaum, Phys. Rev. Lett. 120, 013202 (2018).
  * Graner  _et al._ [2016] B. Graner, Y. Chen, E. G. Lindahl, and B. R. Heckel, Phys. Rev. Lett. 116, 161601 (2016).
  * Swallows  _et al._ [2013] M. D. Swallows, T. H. Loftus, W. C. Griffith, B. R. Heckel, E. N. Fortson, and M. V. Romalis, Phys. Rev. A 87, 012102 (2013).
  * Mantry  _et al._ [2014] S. Mantry, M. Pitschmann, and M. J. Ramsey-Musolf, Physical Review D 90, 10.1103/physrevd.90.054016 (2014).
  * Yuan  _et al._ [1991] V. W. Yuan, C. D. Bowman, J. D. Bowman, J. E. Bush, P. P. J. Delheij, C. M. Frankle, C. R. Gould, D. G. Haase, J. N. Knudson, G. E. Mitchell, S. Penttilä, H. Postma, N. R. Roberson, S. J. Seestrom, J. J. Szymanski, and X. Zhu, Phys. Rev. C 44, 2187 (1991).
  * Sharapov  _et al._ [2000] E. I. Sharapov, J. D. Bowman, B. E. Crawford, P. P. J. Delheij, C. M. Frankle, M. Iinuma, J. N. Knudson, L. Y. Lowie, J. E. Lynch, A. Masaike, Y. Matsuda, G. E. Mitchell, S. I. Penttilä, H. Postma, N. R. Roberson, S. J. Seestrom, S. L. Stephenson, Y.-F. Yen, and V. W. Yuan, Phys. Rev. C 61, 025501 (2000).
  * Auerbach  _et al._ [1996] N. Auerbach, V. V. Flambaum, and V. Spevak, Phys. Rev. Lett. 76, 4316 (1996).
  * Spevak  _et al._ [1997] V. Spevak, N. Auerbach, and V. V. Flambaum, Phys. Rev. C 56, 1357 (1997).
  * Kadmenskii  _et al._ [1982] S. G. Kadmenskii, V. P. Markushev, and V. I. Furman, Soviet Journal of Nuclear Physics 35(2) , 166 (1982).
  * Kadmenskii  _et al._ [1984] S. G. Kadmenskii, V. P. Markushev, Y. P. Popov, and V. I. Furman, Soviet Journal of Nuclear Physics 39(1) , 4 (1984).
  * Chishti  _et al._ [2020] M. Chishti, D. O’Donnell, G. Battaglia, M. Bowry, D. Jaroszynski, B. Singh, M. Scheck, P. Spagnoletti, and J. Smith, Nature Physics 16, 853 (2020).
  * Flambaum  _et al._ [1996] V. V. Flambaum, G. F. Gribakin, and F. M. Izrailev, Phys. Rev. E 53, 5729 (1996).
  * Horoi  _et al._ [1995] M. Horoi, V. Zelevinsky, and B. A. Brown, Phys. Rev. Lett. 74, 5194 (1995).
  * Zelevinsky  _et al._ [1995] V. Zelevinsky, M. Horoi, and B. A. Brown, Physics Letters B 350, 141 (1995).
  * Zaretsky and Sirotkin [1983] D. Zaretsky and V. Sirotkin, Yad. Fiz. 37(1) 45, 1983. Sov. J. Nucl. Phys. 37, 361 (1983).
  * Koonin  _et al._ [1992] S. E. Koonin, C. W. Johnson, and P. Vogel, Phys. Rev. Lett. 69, 1163 (1992).
  * Carlson and Hussein [1993] B. V. Carlson and M. S. Hussein, Phys. Rev. C 47, 376 (1993).
  * Lewenkopf and Weidenmüller [1992] C. H. Lewenkopf and H. A. Weidenmüller, Phys. Rev. C 46, 2601 (1992).
  * Flambaum and Izrailev [1997b] V. V. Flambaum and F. M. Izrailev, Phys. Rev. E 55, R13 (1997b).


20XX Vol. X No. XX, 000–000 11institutetext: Purple Mountain Observatory and
Key Laboratory of Radio Astronomy, Chinese Academy of Sciences, 10 Yuanhua
Road, Nanjing 210033<EMAIL_ADDRESS>
22institutetext: University of Chinese Academy of Sciences, 19A Yuquan Road,
Shijingshan District, Beijing 100049, China
Received 2019/10/22; accepted 2019/12/10

# The Column Density Structure of Orion A Depicted by N-PDF

Yuehui Ma 11 2 2 Hongchi Wang 11 Chong Li 11 2 2 Ji Yang 11

###### Abstract

We have conducted a large-field simultaneous survey of 12CO, 13CO, and C18O
$J=1-0$ emission toward the Orion A giant molecular cloud (GMC) with a sky
coverage of $\sim$ 4.4 deg2 using the PMO-13.7 m millimeter-wavelength
telescope. We use the probability distribution function of the column density
(N-PDF) to investigate the distribution of molecular hydrogen in the Orion A
GMC. The H2 column density, derived from the 13CO emission, of the GMC is
dominated by log-normal distribution in the range from $\sim$4$\times 10^{21}$
to $\sim$1.5$\times 10^{23}$ cm-2 with excesses both at the low-density and
high-density ends. The excess of the low-density end is possibly caused by an
extended and low-temperature ($\sim$10 K) component with velocities in the
range of 5$-$8 km s-1. Compared with the northern sub-regions, the southern
sub-regions of the Orion A GMC contain less gas with column density in
$N_{H_{2}}>1.25\times 10^{22}\ \rm{cm}^{-2}$. The dispersions of the N-PDFs of
the sub-regions are found to correlate with the evolutionary stages of the
clouds across the Orion A GMC. The structure hierarchy of Orion A GMC is
explored with the DENDROGRAM algorithm, and it is found that the GMC is
composed of two branches. All structures except one in the tree have virial
parameters less than 2, indicating self-gravity is important on the spatial
scales from $\sim$0.3 to $\sim$4 pc. Although power-laws and departures from
log-normal distributions are found at the high-density end of N-PDFs of active
star-forming regions, the N-PDFs of structures in the Orion A GMC are
predominantly log-normal on scales from R$\sim$0.4 to 4 pc.

###### keywords:

ISM: clouds — ISM: individual objects (Orion A) — ISM: structure — stars:
formation — surveys — turbulence

##  1 Introduction

Stars form in the cold and dense molecular clouds. The cloud structure is
affected by various physical processes, such as turbulence, cloud-cloud
collisions, and feedbacks from massive stars. The probability distribution
function (PDF) is a simple but potent tool to depict the structure of
molecular clouds, therefore, providing an efffective means to investigate the
various physical processes that influence the structures of molecular clouds.
Theoretical studies (Vazquez-Semadeni 1994; Padoan et al. 1997; Klessen 2000)
of supersonic turbulence suggest that the volume density PDFs ($\rho$-PDFs) of
molecular clouds exhibit log-normal shapes when self-gravity is not important.
This characteristic has been attributed to the influence of some steady and
multiple independent dynamical events that shape the molecular clouds (McKee &
Ostriker 2007). When star formation activities occur, the $\rho$-PDF is
strongly affected by gravity in addition to turbulence, resulting in excess
above the log-normal distribution at the high density end (Kainulainen et al.
2014; Federrath & Klessen 2013). $\rho$-PDF is also proved to be a reliable
indicator for the star formation rate and efficiency of molecular clouds in
different star formation models (Elmegreen 2008; Padoan et al. 2014). However,
it is challenging to constrain the underlying $\rho$-PDFs of molecular clouds
with observations because of the projection effect along the line-of-sight.
Instead, a lot of theoretical and observational studies focused on the column
density PDFs (N-PDFs). Simulations have suggested an evolutionary trend of
molecular clouds from the turbulent dominated stages to the gravity dominated
star-forming stages (Ballesteros-Paredes et al. 2011; Kritsuk et al. 2011).
The gravity first broadens the log-normal N-PDF at the beginning and then
produces a power-law tail in the N-PDF (Ballesteros-Paredes et al. 2011). The
slope of the power-law tail develops over time and gets shallower as the star
formation efficiency (SFE) increases.

Observationally, three tracers of column density are widely used, which are
the near infrared dust extinction, the far-infrared dust emission, and the
optically thin molecular line emission. Comparatively, dust extinction probes
the largest dynamic range of column density and less relies on models (Goodman
et al. 2009a). However, the gas tracer can provide us the kinematic
information of the molecular clouds, which helps to disentangle the structure
of molecular clouds. Log-normal shaped N-PDFs derived with the three tracers
have been observed in the Perseus cloud (Goodman et al. 2009a), and power-law
N-PDFs have been observed in active star-forming regions (Kainulainen et al.
2009; Schneider et al. 2013; Lombardi et al. 2015). In the observations
mentioned above, quiescent clouds often exhibit log-normal N-PDFs in the whole
dynamic range or only have modest excesses at the high density ends. In
contrast, clouds active in star formation have N-PDFs with prominent non-log-
normal or power-law high-density tails arising from the dense star-forming
regions.

As the nearest massive star-forming region (d = 414 pc, Menten et al. 2007),
the Orion A GMC is notable for its complex hierarchical filamentary structure
and the intense star formation activities (Bally 2008). There have been plenty
of large-scale surveys toward the Orion A GMC at multi-wavelengths (Bally et
al. 1987; Sakamoto et al. 1994; Johnstone & Bally 1999; Shimajiri et al. 2011;
Megeath et al. 2012; Fischer et al. 2013; Ripple et al. 2013; Stutz et al.
2013; Berné et al. 2014; Shimajiri et al. 2014; Nishimura et al. 2015; Stutz &
Kainulainen 2015; Großschedl et al. 2018; Hacar et al. 2018; Kong et al. 2018;
Suri et al. 2019). The GMC is mainly composed of a large and dense integral-
shaped-filament (ISF, Bally et al. 1987) in the northern part that contains
star-forming regions OMC 1-4, and a relatively extended and less dense tail in
the southern part that covers the LDN 1641 region. Kong et al. (2018) have
made a brief review of the surveys that covers the ISF in their table 1. The
broad dynamic range of column densities revealed by the previous studies makes
the GMC an ideal case for N-PDF analysis. The N-PDF of Orion A GMC has been
derived with multi-wavelength observations (Kainulainen et al. 2009; Lombardi
et al. 2015; Stutz & Kainulainen 2015; Berné et al. 2014), most of them
presenting a N-PDF of power-law shape. However, most of the previous N-PDF
studies of Orion A are limited to the ISF region except for the $Herschel$
observations (Stutz & Kainulainen 2015) (S15 hereafter). S15 found that all
the N-PDFs in different subregions of the Orion A GMC exhibit power-law
distributions, and the power-law exponent of the N-PDF correlates with the
fraction of the Class 0 protostars. For a panoramic view of the Orion A GMC,
it is valuable to study the N-PDF of the entire Orion A GMC using molecular
line emission and make comparison with N-PDFs derived using dust emission.
Furthermore, for the investigation of the relation between N-PDF forms and the
underlying physics, it is necessary to find out the dependence of the N-PDF
forms and the role of self-gravity on spatial scales. There are several
automatic algorithms that can identify structures in molecular clouds, like
CLUMPFIND (Williams et al. 1994), GAUSSCLUMPS (Stutzki & Guesten 1990),
FELLWALKER (Berry 2015), CPROPS (Rosolowsky & Leroy 2006), and DENDROGRAM
(Rosolowsky et al. 2008). Among the algorithms, DENDROGRAM can reveal the
hierarchical nature of molecular clouds, and it is relatively insensitive to
the input parameters.

In this work, we present a detailed analysis of N-PDF of the Orion A GMC using
the 12CO and 13CO $J=1-0$ emission line data from a large-scale ($\sim$ 4.4
deg2) survey toward the GMC. The N-PDFs of the hierarchical structures and the
role of gravity on different spatial scales in the Orion A GMC have been
studied using the DENDROGRAM technique. In Section 2, we describe the
observations and the data reduction processes. The main results are presented
in Section 3. The N-PDF properties and the importance of gravity in the GMC on
different spatial scales are discussed in Section 4. We summarize our results
and make conclusions in Section 5.

##  2 Observations and data reduction

###  2.1 Observations

The observations were made using the PMO-13.7 m millimeter-wavelength
telescope located at Delingha in China during 2011 June as a pilot survey of
the Milky Way Imaging Scroll Painting project
111http://www.radioast.nsdc.cn/mwisp.php (Su et al. 2019). The PMO-13.7m
telescope is equipped with a nine-beam Superconducting Spectroscopic Array
Receiver (SSAR) that works in the 85-115 GHz frequency band (Shan et al.
2012). With a carefully selected local oscillator (LO) frequency, the 12CO,
13CO, and C18O $J=1-0$ line emission data can be obtained simultaneously with
this instrument. The upper sideband contains the 12CO $J=1-0$ and the lower
sideband contains the 13CO and C18O $J=1-0$ line emission. For each sideband,
a Superconductor-Insulator-Superconductor (SIS) mixer works as the front end
of the receiver. The back-end of the receiver is a Fast Fourier Transform
Spectrometer (FFTS) with 16,384 frequency channels and a total bandwidth of 1
GHz, providing a velocity resolution of 0.17 km s-1 at 110 GHz. At 110 GHz and
115 GHz, the half-power beam width (HPBW) of the telescope is about 52
$\arcsec$ and 50$\arcsec$, respectively.

The observations were conducted in the position-switch on-the-fly (OTF) mode.
The Orion A GMC is divided into twelve 30$\arcmin\times$30$\arcmin$ cells. For
each cell two scannings were made along the directions of right ascension and
declination, respectively. The scanning rate of the observation is 50$\arcsec$
per second and the dump time is 0.3 s. The antenna temperature is calibrated
according to $\rm{T_{MB}=T^{*}_{A}/\eta_{MB}}$ during the data reduction
processes. According to the status report of the PMO-13.7 m telescope, the
main beam efficiencies are 0.44 and 0.48 at the 12CO $J=1-0$ and the 13CO
$J=1-0$ wavelengths, respectively. We use the GILDAS/CLASS software package to
reduce the data, which includes the subtraction of a second order baseline
from each spectrum and the re-griding of the raw data. At the positions where
multiple spectra are taken at different times, a combination of the spectra is
made using the reciprocal square of the RMS noise level of each spectrum as
the weight. The final reduced data cube covers a sky area of $\sim$4.4 deg2,
which contains the majority CO emission of the GMC as revealed by Bally et al.
(1987). However, the spectra at the edge of the surveyed area have low signal-
to-noise ratios and they are removed in the analysis of this work. The trimmed
edge area, indicated with grey lines in Figure 1, accounts for about twenty
percent of the scanned area. The final median RMS noise level of the $\sim$3.5
deg2 effective region is 0.61 K per channel at 115 GHz and 0.37 K per channel
at 110 GHz.

###  2.2 Method of Structure Identification

In this work, we use the DENDROGRAM algorithm to extract “tree” structures
from the H2 column density map and investigate the properties of the N-PDFs of
the hierarchical structures and the importance of self-gravity at various
spatial scales in Orion A in Section 4. The threshold and the difference
between two separate column density peaks for the identification of “leaf” of
the tree structure are set to be 2.0 and 3.0 times the median noise level of
N${}_{H_{2}}$ (see text in Section 3.1), respectively. The minimum pixel
number of a “leaf” structure is set to 200, which physically corresponds to a
clump of the size of $\sim$0.8$\times$0.8 pc2 at the cloud distance. The
outputs of the DENDROGRAM provide masks of the identified structures.

We evaluate the importance of self-gravity in a molecular structure using the
virial parameter $\alpha_{vir}$, which measures the ratio between the total
kinetic energy and the gravitational energy. The widely used definition of
$\alpha_{vir}$ is $5\sigma_{v}^{2}R/GM_{c}$, where $\sigma_{v}$ is the one-
dimensional velocity dispersion, G is the gravitational constant, $R$ and $M$
are the effective radius and mass of the structure, respectively. The critical
value of $\alpha_{vir}$ for a marginally bound structure is 2 (Kauffmann et
al. 2013). The angular size of a structure is defined as the geometric mean of
the FWHMs of its major and minor axes which can be obtained in the output
table of the DENDROGRAM algorithm, and the LTE mass can be obtained using the
total column density within the periphery of the structure. The velocity
dispersion of a molecular structure is obtained by fitting a gaussian function
to the average spectrum of the 13CO emission.

Figure 1: Color-coded image of integrated intensity of the Orion A GMC, with
the 12CO emission in the velocity range from $-$10 to 20 km s-1 in blue, the
13CO emission from 0 to 18 km s-1 in green, and the C18O emission from 0 to 14
km s-1 in red. The grey lines outline the noisy edge of the surveyed area that
has been trimmed in the analysis. Active star-forming regions in the GMC are
indicated with red letters.

##  3 Results

Figure 1 presents the composite intensity map of the 12CO, 13CO, and C18O
$J=1-0$ emission of the Orion A GMC. The 12CO and 13CO emission in the region
are much more extended than the C18O emission, which is mostly concentrated in
the active star-forming regions indicated in Figure 1. Compared with the
emission of the 12CO and 13CO, the C18O emission traces positions with high
column density, therefore, only providing a limited dynamical range for the
N-PDF study. In the present work, we use the 12CO and 13CO data for the N-PDF
analysis and the analysis of the C18O data is deferred to a future paper.

###  3.1 Calculation of the Column Density

The 12CO emission is usually optically thick in typical molecular cloud
environments, while the 13CO emission is often optically thin. Providing the
molecular clouds are under the local thermodynamic equilibrium (LTE)
conditions and the molecules are uniformly excited, the column density of the
clouds can be derived with the 12CO and 13CO spectra.

Firstly, the excitation temperature can be obtained with the peak brightness
temperature of the optically thick 12CO line (Li et al. 2018),

$T_{ex}=5.532[ln(1+\dfrac{5.532}{T_{peak}+0.819})]^{-1}$ (1)

where $T_{peak}$ is the peak brightness temperature of the 12CO $J=1-0$ line.
The optical depth and the column density of the 13CO $J=1-0$ emission can be
derived through (Pineda et al. 2010; Li et al. 2018)

$\tau^{13}_{v}=-ln\left\\{1-\dfrac{T_{MB}(^{13}CO)}{5.29[J(T_{ex})-0.164]}\right\\}$
(2)

$N_{{}^{13}CO}=2.42\times
10^{14}\frac{T_{ex}+0.88}{1-e^{-5.29/T_{ex}}}\int{\tau^{13}_{v}dv}.$ (3)

where $T_{ex}$ is the excitation temperature, and
$J(T_{ex})=[exp(5.29/T_{ex})-1]^{-1}$. According to Pineda et al. (2010), when
the optical depth of the 13CO emission is small, $\tau_{v}^{13}<1$, the
integral item $T_{ex}\int{\tau}dv$ can be approximated to
$\tau_{0}/(1-e^{-\tau_{0}})\int T_{mb}\ dv$, where $\tau_{0}$ is the peak
optical depth of the 13CO emission line. Providing that the abundance ratios
[12C]$/$[13C] $=$ 69 (Wilson 1999) and [CO]$/$[H] $\sim$ 10-4 (Solomon &
Klemperer 1972; Herbst & Klemperer 1973), the [$H_{2}/^{13}$CO] ratio is
expected to be $7\times 10^{5}$, the column density of molecular hydrogen can
be obtained with the following formula (Li et al. 2018)

$N_{H_{2}}=1.694\times
10^{20}\frac{\tau_{0}}{1-e^{-\tau_{0}}}\frac{1+0.88/T_{ex}}{1-e^{-5.29/T_{ex}}}\int{T_{MB}(^{13}CO)dv},$
(4)

The peak brightness temperature of 12CO emission and the integrated intensity
of 13CO emission are calculated in the velocity range from 0 to 18 km s-1.
Observationally, molecular clouds are defined as contiguous structures above a
fixed brightness temperature in the position-position-velocity space. The
molecular line emission from the clouds should have a certain width of
velocity dispersion ($\sigma_{v}\sim 1$ km s-1). To avoid missing weak signals
and to ensure an adequate signal-to-noise ratio, we only use the spectra that
have at least five contiguous velocity channels with intensities higher than 3
$\sigma_{rms}$ for ${}^{12}CO$ for the calculation of excitation temperature,
and 1.5$\sigma_{rms}$ for ${}^{13}CO$ for the calculation of column density.
According to the above criteria, the signal-to-noise ratio of selected 13CO
spectra is above 3.35, which is sufficiently high to exclude the noisy
channels. The three-sigma detection limit of the H2 column density is
$\sim$1.5$\times$1021 cm-2 under the criteria.

###  3.2 Basic Properties of the Orion A GMC

Figure 2(a) shows the spatial distribution of the excitation temperature
derived using Eq. (1). The excitation temperatures of Orion A range from 5 to
107 K. In the ISF and L 1641 N regions, most of the excitation temperatures
are above 25 K. In the OMC 1 region, which lies behind the Orion Nebula, the
temperatures are $\sim$ 50 K on average and reach the maximum 107 K at the
position of Orion-KL. Comparatively, the southern part of the GMC is much
colder except for the south-eastern tail. The distribution of the 13CO optical
depth is shown in Figure 2(b). The optical depths, $\tau^{13}$, range from
0.04 and to 1.35 with a median value of 0.37, and only 23 pixels have optical
depths greater than 1, indicating the 13CO emission in the region is mostly
optically thin.

(a)

(b)

Figure 2: (a) Excitation temperature derived from the 12CO emission. (b)
Distribution of peak optical depth of the 13CO emission.

(a)

(b)

Figure 3: (a) Spatial distribution of H2 column density of the Orion A GMC.
The black contours and the white contours at the middle of ISF correspond to
the lower and higher critical value 4$\times$1021 cm-2 and 1.5$\times$1023
cm-2, respectively, at which the N-PDF departures from lognormal distribution.
(b) N-PDF of H2 column density of the Orion A GMC. The purple dashed line
delineates the fitted Gaussian function of $s=ln(N/<N>)$ within the fitting
range from $N_{H_{2}}=4.4\times 10^{21}$ cm-2 to $N_{H_{2}}=1.4\times 10^{23}$
cm-2, which is indicated with light blue shadow. The vertical and horizontal
dashed lines mark the median detection limit and the level at which the number
of pixels in a bin is ten. The fitted parameters $\mu$ and $\sigma$ of the
Gaussian function are indicated in purple below the horizontal dashed line.
The grey dashed lines are the N-PDFs of the dendrogram structures 5 (narrow)
and 11 (wide), respectively, and the blue dashed line is the summation N-PDF
of the two structures. For structures 5 and 11, see discussions in Section
4.2.

###  3.3 N-PDF of the Entire Orion A GMC

####  3.3.1 Overall properties of the H2 column density distribution

The spatial distribution of the H2 column density of Orion A GMC is displayed
in Figure 3(a). The corresponding pixel-by-pixel statistics is given in 3(b).
The median detection limit based on the selection criteria described in
Section 3.1 is 1.5$\times$1021 cm-2. As shown in Figure 3(a), the northern ISF
and the L 1641 N regions contain more high-density gas than the southern
regions. In Figure 3(b), we totally use 40592 pixels for the statistics and
the mean H2 column density is about 1.5$\times$1022 cm-2, which is much higher
than the threshold of star formation of 6.3$\times$1021 cm-2 (Johnstone et al.
2004; Lada et al. 2010; Kainulainen et al. 2014). The highest column density
is associated with the Orion-KL region, up to a value of
$\sim$3$\times$1023cm-2. The column densities derived with CO line emission in
the Orion A GMC in this work are higher than that derived with the dust tracer
(Stutz & Kainulainen 2015) and are also higher than the column densities
derived with CO in other nearby star-forming regions, like the Perseus
molecular cloud (Goodman et al. 2009a; Kainulainen et al. 2009).

Mathematically, if a variable N follows a log-normal distribution, its natural
logarithm distribution is normal. The log-normal probability density function
has the form

$p(N)=\frac{1}{N}\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(lnN-\mu)^{2}}{2\sigma^{2}}}$
(5)

where $\mu$ and $\sigma$ are the mean and the dispersion of $lnN$,
respectively. In the histogram in Figure 3(b), the H2 column density N is
normalized by its mean, i.e. the statistical variable is $s=ln(N/<N>)$. As
shown in Figure 3(b), $s$ can be well fitted with a Gaussian function rather
than a power-law function in the column density range from 4$\times$1021 cm-2
to 1.5$\times$1023 cm-2. In the log-log space, the log-normal function is a
parabolic curve, as shown by the purple dashed line in Figure 3(b). Beyond the
fitting range, the N-PDF deviates from the log-normal distribution at both the
low-density and the high-density ends. The low-density deviation mostly
corresponds to the diffuse gas at the edge of the cloud, which is outside the
black contour in Figure 3(a), while the high-density excess is associated with
the very active star-forming regions Orion-KL and Orion-Bar.

####  3.3.2 The low-density excess above the log-normal distribution

Figure 4: Channel maps of the 12CO emission with channel width of 1 km s-1 in
the velocity range from 0 to 16 km s-1. The colors start from 1.5
$\sigma_{RMS}$. The black contour in each panel corresponds to N${}_{H_{2}}$ =
4$\times$1021cm-2. The central velocity of each panel is given in the upper-
left corner of each panel. Figure 5: N-PDFs of the two parts of Orion A GMC
divided by the declination of $\delta=-5^{\circ}57\arcmin 00\arcsec$, overlaid
with the fitted log-normal functions as indicated in purple dashed lines. The
fitting range in each panel is indicated with light blue shadow. Figure 6:
Position-velocity diagrams of the 12CO emission from declination of -4.5∘ to
-8.5∘ in steps of 0.25∘. The colors indicate the integrated intensity in units
of K arcdeg from 2 sigma to the maximum of each panel. The declination
interval of each panel is shown in white in the upper-right corner.

In Figure 3(b), the low-density deviation from the log-normal distribution of
the N-PDF is quite obvious and is different from the previous results of the
Orion A region using either dust or gas tracers (Berné et al. 2014; Stutz &
Kainulainen 2015). We present the channel maps of the 12CO emission in Figure
4, and that of 13CO emission in Figure 14, to investigate which velocity
components contribute to the low-density end of the N-PDF. The
4$\times$1021cm-2 contour is overlaid in each panel of Figures 4 and 14.
Outside the 4$\times$1021cm-2 contour, extended diffuse gas is present in the
velocity channels from 5.5 to 8.5 km s-1, more apparent in the 12CO channel
maps than the 13CO channel maps. From Figure 4, it can be seen that the
diffuse gas mainly exists in the region below $\delta=-5^{\circ}57\arcmin
00\arcsec$. Dividing the Orion A GMC by the declination of
$\delta=-5^{\circ}57\arcmin 00\arcsec$, we separately calculated the N-PDFs of
the northern and southern parts of the GMC. The results are presented in
Figure 5. The N-PDFs of the two parts are quite different. The excesses above
log-normal distribution at the low-density end of the N-PDF is more obvious in
the southern portion than in the northern portion of the GMC.

Figure 7: Histograms of the 12CO peak brightness temperature of the EC
counterpart (blue), the main ridge (green), and the entire observed region
(pink) of Orion A GMC. The blue vertical line and the black number on that
line show the median temperature of the EC counterpart.

Figure 6 presents the position-velocity diagrams of the 12CO emission at every
declination interval of 0.25∘. The P-V diagrams of the 13CO emission are given
in Figure 15 in the Appendix. From Figures. 6 and 15, it can be seen that
there is an apparent velocity gradient from north to south for the main cloud,
which is represented by the red color in the panels. However, it can also be
seen from the 7th and 8th panels that there is a weak component at 6 km s-1,
represented by the narrow white strips at the left side of the two panels.
This 6 km s-1 component is different from the main cloud in terms of the
velocity gradient. This component also can be seen in other panels from
$-$6°to $-$7.75°, although it has similar velocity to the main cloud in
regions below $\delta=-$6.5°, therefore cannot be easily distinguished from
the main cloud.

The 6 km s-1 diffuse component is most likely the counterpart of the “Extended
Component (EC)”, defined by Sakamoto et al. (1997). The EC has been fully
observed by Nishimura et al. (2015) and they suggest that the EC has a typical
velocity of $\sim$ 6 km s-1 and is widely distributed to the east of the main
ridge of Orion A GMC. We present the histograms of the 12CO peak brightness
temperature for the EC counterpart (N${}_{H_{2}}<4\times 10^{21}$ cm-2), the
main ridge (N${}_{H_{2}}>4\times 10^{21}$ cm-2), and the entire observed
region of Orion A GMC in Figure 7. The 12CO emission of the EC counterpart has
peak temperatures in the range from 7 to 13 K, which accounts for the low-
temperature excess in the histogram of the entire GMC.

###  3.4 N-PDFs in Sub-regions

####  3.4.1 Calculation of N-PDFs with Star-formation Activities

S15 have studied the column density distribution in Orion A using the
_Herschel_ observations of the dust emission. They used a joint protostar
catalog that contains the PACS bright red sources (PBRS) (Stutz et al. 2013)
and the Class 0, I and Flat-spectrum sources in the Herschel Orion Protostar
Survey (HOPS) (Furlan et al. 2016) to divide the Orion A GMC into eight
0.5°$\times$0.5°sub-regions according to the protostar distribution. They
found that the N-PDF of hydrogen nucleus (NH-PDF) derived from dust emission
has power-law distribution in each sub-region and the index of the power-law
distribution varies with the number fraction of the Class 0 protostars. For
direct comparison, we use the first seven sub-regions in S15 to calculate the
NH-PDFs except for sub-region 3 that is not fully covered by our observation.
The seven sub-regions are indicated with squares in Figure 8(a) and their
corresponding NH-PDFs are presented in Figure 8(b). The physical parameters of
each sub-region, such as the averaged H2 column density, excitation
temperature, line width, LTE mass, number of young stellar objects (YSO), and
star formation efficiency (SFE), are given in Table 1. In addition to the
Class 0 and Class I catalog used by S15, we also included the disk dominated
pre-main-sequence stars (Class II sources) from Megeath et al. (2012, 2016) to
calculate the SFE in each sub-region. The Class II sources are marked with
green dots in Figure 8(a).

(a)

(b)

Figure 8: (a) H2 column density map overlaid with sub-regions 1-7 in S15. The
number of each sub-region is indicated in the upper-right corner of the
squares. The magenta pluses and the green dots show the positions of the
protostars used in S15, and the disk dominated pre-main-sequence stars from
the Spitzer catalog (Megeath et al. 2012, 2016), respectively. (b) Hydrogen
nucleus NH(2N${}_{H_{2}}$)-PDFs of sub-regions 1-7 and the summation N-PDF of
all the sub-regions.

Table 1: Physical parameters of each sub-region

Region | $\overline{N}_{H_{2}}$ | $\mu_{s}$ | $\sigma_{s}$ |  Tex |  Tdust | $\sigma_{v}$ |  MLTE |  n${\rm{{}_{YSO}}^{\prime}}$ |  nYSO | SFE | Mass Fraction
---|---|---|---|---|---|---|---|---|---|---|---
|  ($\times$1022 cm-2) |  |  | (K) | (K) |  (km s-1) |  (M⊙) |  |  |  |
1 | 2.6 | -0.33 | 0.85 | 29.5 | 27.3 | 0.83 | 5638 | 44 | 420 |  3.59$\%$ |  92.38$\%$
2 | 2.0 | -0.28 | 0.76 | 24.7 | 23.0 | 1.03 | 5583 | 25 | 408 |  3.52$\%$ |  88.14$\%$
3 | 2.4 | -0.18 | 0.65 | 29.2 | 18.2 | 1.07 | 6967 | 22 | 166 |  1.18$\%$ |  92.58$\%$
4 | 1.8 | -0.17 | 0.59 | 25.2 | 17.8 | 0.94 | 5265 | 17 | 74 |  0.70$\%$ |  83.34$\%$
5 | 1.2 | -0.10 | 0.47 | 20.4 | 16.4 | 0.99 | 3265 | 31 | 127 |  1.91$\%$ |  68.93$\%$
6 | 1.2 | -0.12 | 0.51 | 18.1 | 15.4 | 0.96 | 3444 | 26 | 98 |  1.40$\%$ |  72.14$\%$
7 | 1.3 | -0.06 | 0.37 | 19.4 | 15.0 | 0.83 | 3634 | 48 | 125 |  1.69$\%$ |  66.97$\%$

Columns 2-4 give the average H2 column density, the mean and dispersion of the
natural logarithm of the normalized H2 column density (N${}_{H_{2}}>4\times
10^{21}$ cm-2) of each sub-region. Columns 5-8 give the median excitation
temperature, median dust temperature from the _Herschel-Planck_ temperature
map (Lombardi et al. 2014), median velocity dispersion of the 13CO emission,
and the LTE mass of each sub-region. The last four columns give the total
number of protostars (Class 0, I and Flat-spectrum sources), the total number
of all YSOs (protostars and Class II sources), the SFE, and the dense gas
(N${}_{H_{2}}>$1.25$\times$1022 cm-2) mass fraction of each sub-region. The
information of the protostars are from the PBRS and HOPS catalogs (Stutz et
al. 2013; Furlan et al. 2016), while the information of Class II sources is
from the Spitzer catalog (Megeath et al. 2012, 2016). The numbers of the
protostars and YSOs are counted within the boundaries of each sub-region,
which are the boxes in Figure 8. The SFE is derived using the assumption that
the average YSO mass is 0.5 M⊙.

(a)

(b)

Figure 9: (a) Relationship between the mean and the dispersion of the
logarithmic normalized column densities in the sub-regions and (b) the same
relationship but excluding the positions where N${}_{H_{2}}$ is lower than
4$\times$1021 cm-2. The solid diagonal in each panel represents the relation
$\mu_{s}=-\sigma_{s}^{2}/2$.

As shown in Figure 8(b), although the NH-PDFs of the sub-regions cannot be
well fitted with log-normal or power-law functions, it can be seen that the
width of the NH-PDF gradually increases from south to north, which resembles
the trend of the increasingly flattened power-law NH-PDFs derived using the
dust emission (fig. 2 in S15). However, unlike the dust-derived NH-PDFs that
can be fitted well by power-law distributions at the high density end, the
gas-derived NH-PDFs show quick decline at the high density end, but have
relatively extended log-normal-like low-density tails. There is a trend of
increasing column density from the south to the north sub-regions. The
southernmost sub-regions 5$-$7 have median and maximum H column densities of
2.5$\times$1022 and 7.9$\times$1022 cm-2, respectively, which are a factor of
about 1.5 lower than the northern sub-regions. As listed in the last column of
Table 1, only $\sim$70$\%$ of the H2 mass is distributed at positions with H2
column density above 1.25$\times$1022 cm-2 in sub-regions 5$-$7, while this
ratio is $\sim$90$\%$ in the northernmost three sub-regions. This difference
in column density between the south and north sub-regions may be in part
caused by the different extents of CO depletion in these sub-regions.
According to chemical evolutionary models of molecular clouds (Bergin et al.
1995), the gas-phase CO abundance is strongly affected by the depletion effect
in cold (T${}_{dust}\sim$10 K even up to 20 K), dense (n${}_{H_{2}}>$ 104
cm-3), and well-shielded (A${}_{v}\sim$10 mag, corresponds to N${}_{H_{2}}=$
9.4$\times$1021 cm-2) (Bohlin et al. 1978) regions of molecular clouds, and it
is dominated by thermal evaporation in regions with dust temperatures higher
than 22 K. It can be seen from Table 1 that the dust temperatures in the three
southern sub-regions are less than 17 K, indicating that the southern sub-
regions are significantly colder than the northern sub-regions. Ripple et al.
(2013) have studied the relationship between the 13CO column density and the
extinction in Orion A and found that the 13CO column density in the southern
portion of the GMC, which contains sub-regions 5$-$7 and part of sub-region 4
in this work, no longer increases with extinction when A${}_{v}\sim$ 10 mag.
Compared with sub-regions 1$-$4, the lack of high column density gas in sub-
regions 5$-$7 is consistent with the scenario that the 13CO abundance is
reduced by depletion in the cold and well-shielded southern sub-regions.

Whether an N-PDF is intrinsically a log-normal distribution can be tested with
the relationship $\mu_{s}=-\sigma_{s}^{2}/2$ resulting from the normalization
of the column density (Goodman et al. 2009a) (see details in the Appendix).
The relationship between the mean and the dispersion of the normalized
logarithmic column densities of sub-regions 1$-$7 is shown in Figure 9(a). The
diagonal in Figure 9 represents the relationship $\mu_{s}=-\sigma_{s}^{2}/2$
and intrinsic log-normal distributions should fall on this line. As shown in
Figure 9(a), the N-PDFs of the sub-regions more or less deviate from the mean-
dispersion relationship for log-normal functions, and the N-PDFs of the
southern sub-regions are closer to the line than the northern sub-regions. As
discussed in section 3.3.1, the low-density excess in the N-PDF of the GMC
probably originates from the EC counterpart of the Orion A GMC. The
$\mu_{s}-\sigma_{s}$ relationship when the low-density is excluded shows that
the N-PDFs of all sub-regions distribute much closer to the log-normal line
(Figure 9(b)), suggesting intrinsic log-normal distributions of the H2 column
density in these sub-regions.

####  3.4.2 The evolutionary trend of the N-PDF

(a)

Figure 10: Relation between the dispersion of the natural logarithmic
normalized H2 column density, $\sigma_{s}$, and the fraction of protostars
across the sub-regions. The color and black symbols distinguish between the
two cases where the calculation of $\sigma_{s}$ includes (color) and does not
include (black) the positions of N${}_{H_{2}}<$4$\times$1021 cm-2. The dashed
lines give the linear fit to the data points. The Pearson correlation
coefficient between $\sigma_{s}$ and the protostar fraction is $-$0.90 (color)
and $-$0.93 (black).

In Table. 1, we summarize the total number of the YSOs and the total number of
protostars (Class 0, I, and the Flat-spectrum sources) of each sub-region. The
protostars are at an earlier evolutionary stage than the Class II YSOs (Dunham
et al. 2014). Figure 10(a) shows the relationship between the protostar
faction and $\sigma_{s}$ for the sub-regions. If we assume that the Orion A
GMC is undergoing star formation with a constant star formation rate, the
fraction of protostars is related to the evolutionary status of the molecular
cloud. Higher protostar fraction indicates an earlier evolutionary stage.
Figure 10(a) shows that the protostar fraction is anti-correlated with
$\sigma_{s}$, indicating that the column density structure of the molecular
clouds in the Orion A region is coupled with the evolutionary stages of star
formation, with such a trend that a later evolutionary stage corresponds to a
broader column density dispersion.

##  4 Hierarchical Column Density Structure

###  4.1 Hierarchical “Tree” and the Importance of Self-gravity on Different
Spatial Scales

Figure 11: Resulting “tree” form the DENDROGRAM algorithm, where the values of
the x axis have no physical meaning, and different colors correspond to
separate “trunk” structures of the “tree”.

The dendrogram of the tree structure of the Orion A GMC is given in Figure 11.
The observed region contains two “trunk” structures, as shown in Figure 11.
The right trunk is the main cloud, while the left trunk is the small cloud in
the top-left corner of Figure 3(a). To investigate the properties of
structures at various scales, we only use the right trunk in the following
analysis. There are 23 individual structures in the right trunk, among which
twelve are “leaves”, ten are “branches”, and one is the trunk itself. The
sizes of the leaves span a large range and there is no favorable scale for the
leaf structure. The leaves contain both large filaments and small clumps. The
“genealogy” diagrams of the spatial distribution of all the structures are
given in Figure 12. As shown in Figure 12, the ridge of the Orion A GMC is
composed of two main structures, branch 11 that contains the entire ISF and
the L1641 N regions, and branch 5 that contains the southern tail of Orion A
GMC. Comparatively, the density structures are richer in the northern branch
than in the southern branch, and there are more levels in branch 11 than
branch 5. Downward from branch 3, the northern and southern structures of the
GMC merge into one cloud.

The significance of self-gravity for the dynamics of dense cores on scales of
$\sim$0.1 pc is well-known (Larson 2003; di Francesco et al. 2007). However,
turbulent simulations without self-gravity successfully yield density and
velocity distributions of molecular clouds that resemble observations
(Vázquez-Semadeni & García 2001; Padoan & Nordlund 2002; Federrath et al.
2009). Goodman et al. (2009b) have made a “dendrogram” analysis of the L 1448
cloud using 13CO observations and found that self-gravity is important on the
full range of scales of the structures in that cloud, from subparsec to $\sim$
1 pc. We present the genealogy diagram of the derived virial parameters and
N-PDFs of the structures identified in this work in Figure 13. The panels in
Figure 13 that correspond to gravitationally bound structures are filled with
light blue. Except for leaf 1 that corresponds to a clump at the southernmost
end of the GMC, all structures, on scales from 0.3 pc (leaf 20) to 4.3 pc
(trunk 2), are gravitationally bound. This result suggests that self-gravity
is also important on much larger scales, i.e., on the clump- to cloud- scales.
Our result is consistent with the results of Goodman et al. (2009b) and
further highlights the necessity of including self-gravity in simulations of
molecular clouds on pc scales.

###  4.2 Properties of N-PDFs on Different Spatial Scales

The N-PDFs of the dendrogram structures shown in Figure 13 are uniformly
normalized by the mean column density of the Orion A GMC, which is
1.5$\times$1022 cm-2. We have performed log-normal or power-law fittings on
the N-PDFs as appropriate. The results are indicated with dotted and dashed
lines in Figure 13. We note that for most of the dendrogram structures in
Figure 13, the adopted uniform mean column density for the normalization of
N-PDF is lower than the actual mean column density of the individual
structures, which translates the N-PDFs shown in Figure 13 to the right along
the $\ln N/<N>$ axis and therefore increases the value of the fitted $\mu$
parameter. Due to this reason, most of the fitted $\mu$ parameters have
positive values.

The leaves of branch 11 can be divided into three categories, leaves 18 and 12
with pure log-normal N-PDFs, leaves 21 and 14 that have log-normal N-PDFs with
high-density excesses, and leaves 20 and 15 that have power-law N-PDFs at the
high-density end. As shown in Figure 12, leaf 18 is located to the south of
the ISF structure, and corresponds to the center of sub-region 3 in Section
3.4. A relatively low number of protostars are associated with leaf 18,
indicating inactive star formation within leaf 18. Leaf 12 is a relatively
faint filament located at the southernmost end of branch 11, and this filament
is also quiescent in star formation as indicated by the few associated Class
II YSOs in Figure 8. Leaves 21 and 14 correspond to the OMC-1/2/3 regions,
which contain a cluster of intermediate-mass protostars and the NGC 1999
region. The high-density excesses in leaves 21 and 14 are concentrated in
regions with s$\sim$ 2.3$-$3.2 ($N_{H_{2}}\sim 1.5-3.7\times 10^{23}$ cm-2)
and s$\sim$ 1.2$-$1.5 ( $N_{H_{2}}\sim 5.0-6.7\times 10^{22}$ cm-2),
respectively. The high-density regions account for 12.8% of the mass and 4.2%
of the area for leaf 21. For leaf 14, the high-density mass and area fractions
are 23.7% and 16.6%, respectively. Leaves 20 and 15 correspond to OMC-4, which
includes a dense star-forming region and the H II region L 1641 N. Leaf 20 is
possibly on a later evolutionary stage than leaf 21, as indicated by the pure
power-law PDF and the lower protostar fraction as discussed in Section 3.4.2.
The N-PDFs of branches 19, 17, 16, 13, and 11 can be well-fitted with log-
normal functions below $s=2.3$. The high-density excesses in the N-PDFs of
these branches are caused by the pixels of density above $N_{H_{2}}\sim
1.5\times 10^{23}$ cm-2 which is confined to the Orion KL region in leaf 21
and is shown with the white contour in Figure 3a. The high-density pixels with
$s>2.3$ only account for a small mass and area proportion of the branch N-PDFs
below leaf 21. The mass proportion decreases from 12.8% of leaf 21 to 2.6% of
trunk 2, while the corresponding area proportion decreases from 4.2% to 0.2%.
Therefore, the majority of the structures of branch 11 are dominated by log-
normal N-PDFs and departures from log-normal N-PDFs at high-densities only
occur in active star-forming regions.

Except for leaf 7 that corresponds to the H II region L 1641 S, all leaves in
branch 5 have log-normal N-PDFs without high-density excesses. The N-PDF of
leaf 4 is well fitted with double log-normal functions, indicating the
possible existence of mixed components along the line-of-sight. We can see
from the last two panels in Figure 6 that there are two separate velocity
components at the location of leaf 4, one of which is in the range from 5 to 8
km s-1, the same as the EC counterpart. However, it is most likely that the
velocity separation in the leaf 4 region is a local dynamic property and not
caused by the EC gas. The N-PDFs of branches 9 and 5 are fitted with log-
normal functions, while the N-PDF of branch 6 that contains the H II region L
1641 S is fitted with power-law at its high density tail.

Generally, the $N_{H_{2}}$ dynamical ranges of the structures above branch 5
are smaller than those of the structures above branch 11. From leaf 21 to
branch 11, the $\mu$ of the fitted log-normal function decreases from
$\sim$1.5 to $\sim$0.3, and the $\sigma_{s}$ increases from $\sim$0.3 to
$\sim$0.70. From leaf 8 to branch 5, the $\mu$ and the $\sigma_{s}$ of the
fitted log-normal functions vary gradually from $\sim$0.4 to $\sim$0.1 and
$\sim$0.1 to $\sim$0.3, respectively. The ranges of $\mu$ and $\sigma_{s}$ of
structures in branch 5 are much smaller than those in branch 11, indicating
the southern filament (branch 5) of the Orion A GMC is at an earlier
evolutionary stage than the northern filament (branch 11), which is consistent
with the results in Section 3.4. We also plot the fitted log-normal N-PDFs of
branches 5 and 11, and their summation in Figure 3. The N-PDFs of the two
filaments account for the high-density part of the column density distribution
of Orion A, supporting that the Orion A molecular cloud is mainly composed of
these two filaments.

Below branches 11 and 5, there are two leaves and three branches. Leaves 1 and
22 correspond to a southernmost and a northernmost clump, respectively, both
inactive in star formation. The N-PDF of leaf 1 is fitted with a log-normal
function while the N-PDF of leaf 22 has an irregular shape. Branch 3 is mainly
the combination of filaments 5 and 11, and is approximately the main ridge of
the Orion A GMC. The N-PDFs of branches 3, 0, and 2 are fitted with log-normal
functions and getting closer to the overall N-PDF of the whole cloud
successively. The high-density excesses in the N-PDFs of these three branches
are caused by leaf 21, while the low-density excess of branch 2 is caused by
the EC gas, as discussed in Section 3.3.2.

In summary, regardless the sub-structure is in the northern or southern part
of the GMC, power-laws and deviations from log-normal distributions at the
high-density tail of N-PDFs only exist in active star-forming regions. The
majority of the structures in Orion A GMC have log-normal dominated N-PDFs on
scales from $R\sim 0.4$ to 4 pc. The predominant existence of log-normal
N-PDFs in structures on such a broad range of scales suggests a dominant role
of turbulence on these spatial scales.

Figure 12: “Genealogy” diagram of the sub-structures in the dendrogram in
Figure 11 and their spatial distribution. The number in each panel is the
index of the structure in the dendrogram. The green line shows the
relationship between the structures. Figure 13: Same configuration as Figure
12, but with N-PDF rather than the spatial distribution of each structure
being presented in each panel. The log-normal and power-law fit to the N-PDFs
are indicated with the purple dotted line and the blue dashed line,
respectively. The fitted log-normal parameters are given in each panel. The
virial parameter of each structure is shown in the upper-left corner in each
panel. The panels corresponding to structures with $\alpha<2$ are filled with
light blue.

##  5 Summary and Conclusion

We have conducted a large-field survey of 12CO, 13CO, and C18O $J=1-0$
emission toward the Orion A GMC with a sky coverage of $\sim$ 4.4 deg2. Using
the 12CO and 13CO data, we have investigated the N-PDFs of H2 column density
of the GMC and compared the NH-PDFs with previous results from _Herschel_
observations in seven sub-regions. The importance of self-gravity and
properties of N-PDFs on different spatial scales are studied using the
DENDROGRAM method. The main results are summarized as follows.

  1. 1.

The H2 column density N-PDF for the entire GMC is fitted with a log-normal
function in the range from $\sim$4$\times 10^{21}$ to $\sim$1.5$\times
10^{23}$ cm-2, with excesses at both the low-density and the high-density
ends. The excess at the high-density end corresponds to the Orion KL region,
and the excess at the low-density end is possibly caused by an extended and
low-temperature component ($\sim$10 K) with velocities in the range of 5$-$8
km s-1.

  2. 2.

To compare with the results from _Herschel_ observations, we divided the Orion
A GMC into seven sub-regions. The N-PDFs of the sub-regions are irregular in
shape, but are log-normal distributed at the high column density end. Compared
with the northern sub-regions, the three southern sub-regions have less high-
density gas, which may be in part caused by the CO depletion effect. The H2
column density structure is coupled with the evolutionary stage of clouds in
Orion A GMC. Broader column density dispersions correspond to later
evolutionary stages.

  3. 3.

In terms of structure hierarchy, Orion A GMC is mainly composed of two
filamentary structures located in the north and south, respectively. All the
structures in the dendrogram tree except leaf 1 are gravitationally bound
across spatial scales from 0.3 to 4.3 pc. Although power-laws and departures
from log-normal distributions exist in N-PDFs of structures in five active
star-forming regions mostly with scales $<0.6$ pc, the N-PDFs of structures in
the Orion A GMC are predominantly log-normal on scales from $R\sim 0.4$ to 4
pc.

###### Acknowledgements.

We thank the PMO-13.7 m telescope staffs for their supports during the
observation. MWISP project is supported by National Key R$\&$D Program of
China under grant 2017YFA0402701 and Key Research Program of Frontier Sciences
of CAS under grant QYZDJ-SSW-SLH047. Y. Ma acknowledges supports by NSFC
grants 11503086 and 11503087. This work makes use of the SIMBAD database,
operated at CDS, Strasbourg, France.

##  Appendix A The relationship between the mean and the dispersion of the
logarithmic normalized column density

As the definition of the log-normal distribution introduced in Eq. (5) in
Section 3.3, the mathematical expectation of N is

$\centering
E(N)=\int_{0}^{+\infty}Nf(N)dN=\int_{0}^{+\infty}\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(lnN-\mu)^{2}}{2\sigma^{2}}}dN.\@add@centering$
(6)

where $N$, $\sigma$, and $\mu$ has the same meaning as in Eq.(5). If we let
$t=(ln\ N-\mu)/(\sqrt{2}\sigma)$, then $N=exp(\sqrt{2}\sigma t+\mu)$ and
$dN=\sqrt{2}\sigma exp(\sqrt{2}\sigma t+\mu)\ dt$. The expectation can be
written as

$E(N)=\int_{-\infty}^{+\infty}\frac{1}{\sqrt{\pi}}e^{-t^{2}+\sqrt{2}\sigma
t+\mu}dt=\frac{1}{\sqrt{\pi}}\int_{-\infty}^{+\infty}e^{-(t-\frac{\sqrt{2}}{2}\sigma)^{2}+\mu+\frac{\sigma^{2}}{2}}dt.$
(7)

Again, after make the substitution of $u=t-\sqrt{2}/2\sigma$ and make use of
the integral formula $\int_{-\infty}^{+\infty}e^{ax^{2}}dx=\sqrt{\pi/a}$ the
above equation is finally converted to

$E(N)=\frac{1}{\sqrt{\pi}}e^{\mu+\frac{\sigma^{2}}{2}}\int_{-\infty}^{+\infty}e^{-u^{2}}du=e^{\mu+\frac{\sigma^{2}}{2}}.$
(8)

Because of the normalization when calculating the N-PDF in this work, the
mathematical expectation of N should be unity which means
$\mu+\sigma^{2}/2=0$, i.e., $\mu=-\sigma^{2}/2$.

##  Appendix B Figures

We present the velocity channel maps, p-v diagrams of the 13CO $J=1-0$
emission in Figs. 14 and 15, respectively.

Figure 14: Same as Figure 4 but for the 13CO emission. Figure 15: Same as
Figure 6 but for the 13CO emission.

## References

  * Ballesteros-Paredes et al. (2011) Ballesteros-Paredes, J., Vázquez-Semadeni, E., Gazol, A., et al. 2011, MNRAS, 416, 1436
  * Bally (2008) Bally, J. 2008, Overview of the Orion Complex, ed. B. Reipurth, Handbook of Star Forming Regions, Volume I, ed. B. Reipurth, 459
  * Bally et al. (1987) Bally, J., Langer, W. D., Stark, A. A., & Wilson, R. W. 1987, ApJ, 312, L45
  * Bergin et al. (1995) Bergin, E. A., Langer, W. D., & Goldsmith, P. F. 1995, ApJ, 441, 222
  * Berné et al. (2014) Berné, O., Marcelino, N., & Cernicharo, J. 2014, ApJ, 795, 13
  * Berry (2015) Berry, D. S. 2015, Astronomy and Computing, 10, 22
  * Bohlin et al. (1978) Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132
  * di Francesco et al. (2007) di Francesco, J., Evans, N. J., I., Caselli, P., et al. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 17
  * Dunham et al. (2014) Dunham, M. M., Stutz, A. M., Allen, L. E., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 195
  * Elmegreen (2008) Elmegreen, B. G. 2008, ApJ, 672, 1006
  * Federrath & Klessen (2013) Federrath, C., & Klessen, R. S. 2013, ApJ, 763, 51
  * Federrath et al. (2009) Federrath, C., Klessen, R. S., & Schmidt, W. 2009, ApJ, 692, 364
  * Fischer et al. (2013) Fischer, W. J., Megeath, S. T., Stutz, A. M., et al. 2013, Astronomische Nachrichten, 334, 53
  * Furlan et al. (2016) Furlan, E., Fischer, W. J., Ali, B., et al. 2016, ApJS, 224, 5
  * Goodman et al. (2009a) Goodman, A. A., Pineda, J. E., & Schnee, S. L. 2009a, ApJ, 692, 91
  * Goodman et al. (2009b) Goodman, A. A., Rosolowsky, E. W., Borkin, M. A., et al. 2009b, Nature, 457, 63
  * Großschedl et al. (2018) Großschedl, J. E., Alves, J., Meingast, S., et al. 2018, A&A, 619, A106
  * Hacar et al. (2018) Hacar, A., Tafalla, M., Forbrich, J., et al. 2018, A&A, 610, A77
  * Herbst & Klemperer (1973) Herbst, E., & Klemperer, W. 1973, ApJ, 185, 505
  * Johnstone & Bally (1999) Johnstone, D., & Bally, J. 1999, ApJ, 510, L49
  * Johnstone et al. (2004) Johnstone, D., Di Francesco, J., & Kirk, H. 2004, ApJ, 611, L45
  * Kainulainen et al. (2009) Kainulainen, J., Beuther, H., Henning, T., & Plume, R. 2009, A&A, 508, L35
  * Kainulainen et al. (2014) Kainulainen, J., Federrath, C., & Henning, T. 2014, Science, 344, 183
  * Kauffmann et al. (2013) Kauffmann, J., Pillai, T., & Goldsmith, P. F. 2013, ApJ, 779, 185
  * Klessen (2000) Klessen, R. S. 2000, ApJ, 535, 869
  * Kong et al. (2018) Kong, S., Arce, H. G., Feddersen, J. R., et al. 2018, ApJS, 236, 25
  * Kritsuk et al. (2011) Kritsuk, A. G., Norman, M. L., & Wagner, R. 2011, ApJ, 727, L20
  * Lada et al. (2010) Lada, C. J., Lombardi, M., & Alves, J. F. 2010, ApJ, 724, 687
  * Larson (2003) Larson, R. B. 2003, Reports on Progress in Physics, 66, 1651
  * Li et al. (2018) Li, C., Wang, H., Zhang, M., et al. 2018, ApJS, 238, 10
  * Lombardi et al. (2015) Lombardi, M., Alves, J., & Lada, C. J. 2015, A&A, 576, L1
  * Lombardi et al. (2014) Lombardi, M., Bouy, H., Alves, J., & Lada, C. J. 2014, A&A, 566, A45
  * McKee & Ostriker (2007) McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565
  * Megeath et al. (2012) Megeath, S. T., Gutermuth, R., Muzerolle, J., et al. 2012, AJ, 144, 192
  * Megeath et al. (2016) Megeath, S. T., Gutermuth, R., Muzerolle, J., et al. 2016, AJ, 151, 5
  * Menten et al. (2007) Menten, K. M., Reid, M. J., Forbrich, J., & Brunthaler, A. 2007, A&A, 474, 515
  * Nishimura et al. (2015) Nishimura, A., Tokuda, K., Kimura, K., et al. 2015, ApJS, 216, 18
  * Padoan et al. (2014) Padoan, P., Federrath, C., Chabrier, G., et al. 2014, 77
  * Padoan & Nordlund (2002) Padoan, P., & Nordlund, Å. 2002, ApJ, 576, 870
  * Padoan et al. (1997) Padoan, P., Nordlund, A., & Jones, B. J. T. 1997, MNRAS, 288, 145
  * Pineda et al. (2010) Pineda, J. L., Goldsmith, P. F., Chapman, N., et al. 2010, ApJ, 721, 686
  * Ripple et al. (2013) Ripple, F., Heyer, M. H., Gutermuth, R., Snell, R. L., & Brunt, C. M. 2013, MNRAS, 431, 1296
  * Rosolowsky & Leroy (2006) Rosolowsky, E., & Leroy, A. 2006, PASP, 118, 590
  * Rosolowsky et al. (2008) Rosolowsky, E. W., Pineda, J. E., Kauffmann, J., & Goodman, A. A. 2008, ApJ, 679, 1338
  * Sakamoto et al. (1997) Sakamoto, S., Hasegawa, T., Hayashi, M., Morino, J.-I., & Sato, K. 1997, ApJ, 481, 302
  * Sakamoto et al. (1994) Sakamoto, S., Hayashi, M., Hasegawa, T., Handa, T., & Oka, T. 1994, ApJ, 425, 641
  * Schneider et al. (2013) Schneider, N., André, P., Könyves, V., et al. 2013, ApJ, 766, L17
  * Shan et al. (2012) Shan, W., Yang, J., Shi, S., et al. 2012, IEEE Transactions on Terahertz Science and Technology, 2, 593
  * Shimajiri et al. (2011) Shimajiri, Y., Kawabe, R., Takakuwa, S., et al. 2011, PASJ, 63, 105
  * Shimajiri et al. (2014) Shimajiri, Y., Kitamura, Y., Saito, M., et al. 2014, A&A, 564, A68
  * Solomon & Klemperer (1972) Solomon, P. M., & Klemperer, W. 1972, ApJ, 178, 389
  * Stutz & Kainulainen (2015) Stutz, A. M., & Kainulainen, J. 2015, A&A, 577, L6
  * Stutz et al. (2013) Stutz, A. M., Tobin, J. J., Stanke, T., et al. 2013, ApJ, 767, 36
  * Stutzki & Guesten (1990) Stutzki, J., & Guesten, R. 1990, ApJ, 356, 513
  * Su et al. (2019) Su, Y., Yang, J., Zhang, S., et al. 2019, ApJS, 240, 9
  * Suri et al. (2019) Suri, S., Sánchez-Monge, Á., Schilke, P., et al. 2019, A&A, 623, A142
  * Vazquez-Semadeni (1994) Vazquez-Semadeni, E. 1994, ApJ, 423, 681
  * Vázquez-Semadeni & García (2001) Vázquez-Semadeni, E., & García, N. 2001, ApJ, 557, 727
  * Williams et al. (1994) Williams, J. P., de Geus, E. J., & Blitz, L. 1994, ApJ, 428, 693
  * Wilson (1999) Wilson, T. L. 1999, Reports on Progress in Physics, 62, 143