# On the Three-Dimensional Structure of Local Molecular Clouds

Catherine Zucker Center for Astrophysics $\mid$ Harvard & Smithsonian, 60
Garden St., Cambridge, MA, USA 02138  Alyssa Goodman Center for Astrophysics
$\mid$ Harvard & Smithsonian, 60 Garden St., Cambridge, MA, USA 02138  João
Alves University of Vienna, Department of Astrophysics, Türkenschanzstraße 17,
1180 Vienna, Austria  Shmuel Bialy Center for Astrophysics $\mid$ Harvard &
Smithsonian, 60 Garden St., Cambridge, MA, USA 02138  Eric W. Koch University
of Alberta, Department of Physics, 4-183 CCIS, Edmonton AB T6G 2E1, Canada
Center for Astrophysics $\mid$ Harvard & Smithsonian, 60 Garden St.,
Cambridge, MA, USA 02138  Joshua S. Speagle Department of Statistical
Sciences, University of Toronto, Toronto, ON M5S 3G3, Canada  David A. Dunlap
Department of Astronomy & Astrophysics, University of Toronto, Toronto, ON M5S
3H4, Canada  Dunlap Institute for Astronomy & Astrophysics, University of
Toronto, Toronto, ON M5S 3H4, Canada  Center for Astrophysics $\mid$ Harvard &
Smithsonian, 60 Garden St., Cambridge, MA, USA 02138  Michael M. Foley Center
for Astrophysics $\mid$ Harvard & Smithsonian, 60 Garden St., Cambridge, MA,
USA 02138  Douglas Finkbeiner Center for Astrophysics $\mid$ Harvard &
Smithsonian, 60 Garden St., Cambridge, MA, USA 02138  Reimar Leike Max Planck
Institute for Astrophysics, Karl-Schwarzschildstraße 1, 85748 Garching,
Germany  Ludwig-Maximilians-Universität, Geschwister-Scholl Platz 1, 80539
Munich, Germany  Torsten En${\ss}$lin Max Planck Institute for Astrophysics,
Karl-Schwarzschildstraße 1, 85748 Garching, Germany  Ludwig-Maximilians-
Universität, Geschwister-Scholl Platz 1, 80539 Munich, Germany  Joshua E.G.
Peek Space Telescope Science Institute, 3700 San Martin Dr, Baltimore, MD
21218, USA  Department of Physics & Astronomy, Johns Hopkins University,
Baltimore, MD 21218, USA  Gordian Edenhofer Max Planck Institute for
Astrophysics, Karl-Schwarzschildstraße 1, 85748 Garching, Germany  Ludwig-
Maximilians-Universität, Geschwister-Scholl Platz 1, 80539 Munich, Germany

###### Abstract

We leverage the 1 pc spatial resolution of the Leike et al. (2020) 3D dust map
to characterize the three-dimensional structure of nearby molecular clouds
($d\lesssim 400$ pc). We start by “skeletonizing” the clouds in 3D volume
density space to determine their “spines,” which we project on the sky to
constrain cloud distances with $\approx 1\%$ uncertainty. For each cloud, we
determine an average radial volume density profile around its 3D spine and fit
the profiles using Gaussian and Plummer functions. The radial volume density
profiles are well-described by a two-component Gaussian function, consistent
with clouds having broad, lower-density outer envelopes and narrow, higher-
density inner layers. The ratio of the outer to inner envelope widths is
$\approx 3:1$. We hypothesize that these two components may be tracing a
transition between atomic and diffuse molecular gas or between the unstable
and cold neutral medium. Plummer-like models can also provide a good fit, with
molecular clouds exhibiting shallow power-law wings with density, $n$, falling
off like $n^{-2}$ at large radii. Using Bayesian model selection, we find that
parameterizing the clouds’ profiles using a single Gaussian is disfavored. We
compare our results with 2D dust extinction maps, finding that the 3D dust
recovers the total cloud mass from integrated approaches with fidelity,
deviating only at higher levels of extinction ($A_{V}\gtrsim 2-3$ mag). The 3D
cloud structure described here will enable comparisons with synthetic clouds
generated in simulations, offering unprecedented insight into the origins and
fates of molecular clouds in the interstellar medium.

molecular clouds, solar neighborhood, interstellar reddening, astronomy data
visualization

††journal: ApJ††software: astropy (Astropy Collaboration et al., 2018), bokeh
(Bokeh Development Team, 2018), dynesty (Speagle, 2020), glue (Robitaille et
al., 2017), numpy (Harris et al., 2020), PyVista (Sullivan & Kaszynski, 2019),
scipy (Virtanen et al., 2020), FilFinder (Koch & Rosolowsky, 2015), scikit-
image (van der Walt et al., 2014), skan (Nunez-Iglesias et al., 2018), RadFil
(Zucker & Chen, 2018)

##  1 Introduction

Much of our understanding of the structure and properties of molecular clouds
is obtained via either two-dimensional (2D) observations of dust in emission
or extinction (see André et al., 2010; Lombardi, 2009; Planck Collaboration et
al., 2011) or three-dimensional (3D) spectral-line observations of gas (see
Dame et al., 2001), where the third dimension is not distance, but radial
velocity. Until very recently, any insights into the true 3D spatial structure
of molecular clouds were available only in numerical simulations, and were
limited by the extent to which these simulations could approximate the
observed properties of real clouds. Even then, detailed comparisons between
observations and simulations could, again, only be made either in 2D projected
space, or via spectral line maps, where the correspondence of intensity
features in position-position-velocity (p-p-v) space to “real” density
features in position-position-position (p-p-p) space has been shown to be
fraught with complications (Beaumont et al., 2013).

Interstellar dust reddens starlight, so by simultaneously modeling the
distance and extinction to large numbers of stars, it is possible to infer a
self-consistent 3D distribution of stars and dust that produces the observed
distribution of stellar colors. Prior to the advent of Gaia (Gaia
Collaboration et al., 2018) in 2018, accurate (e.g. astrometric) measurements
of stars’ distances independent from measures of their colors were scarce, so
transforming the 2D sky into a 3D physical picture of the distribution of
stars and clouds that form them was extraordinarily challenging statistically.
Several path-breaking 3D dust maps were created prior to Gaia (Marshall et
al., 2006; Sale et al., 2014; Green et al., 2015), but Gaia’s astrometric
constraints on stellar distances have allowed for significant gains in the
spatial resolution of 3D dust maps. And the fidelity of 3D maps has improved
even more as a variety of photometric catalogs, which offer billions of
stellar color measurements needed to constrain extinction in fits, have been
crossmatched with Gaia, (see e.g. Green et al., 2019; Lallement et al., 2019;
Rezaei Kh. et al., 2018; Chen et al., 2019; Hottier et al., 2020; Leike &
Enßlin, 2019). Customization of the 3D dust mapping framework, where distance
resolution can be enhanced at the expense of spatial resolution, has yielded
highly accurate distance measurements for nearby molecular clouds, with
typical distance uncertainties of $\approx 5\%$ (Zucker et al., 2019, 2020;
Chen et al., 2020; Yan et al., 2019), a factor of five improvement over
typical distance uncertainties in the pre-Gaia era.

While these new catalogs of molecular cloud distances have shown good
agreement with maser distances (Zucker et al., 2020), revealed new Galactic-
scale features (e.g. the Radcliffe Wave; Alves et al., 2020), and even mapped
distance gradients within an individual molecular cloud (Zucker et al., 2019,
2018b), up until now, catalogs have not had the resolution to resolve the
detailed 3D structure (e.g. thickness, orientation) of clouds on scales
comparable to numerical simulations. Most existing approaches group batches of
stars into pixels on the sky, and then fit stellar distance and extinction
measurements in each pixel to obtain the distribution of dust along the line
of sight. A 3D map can then be created by combining 1D line-of-sight
information for all 2D sky pixels. However, the line-of-sight dust
distributions suffer from a distortion effect, where the angular resolution on
the sky is always superior to the distance resolution along the line of sight,
resulting in dust clouds appearing to extend radially outward from the Sun in
an unphysical manner.

Recently Leike et al. (2020) presented a new 3D dust map of the solar
neighborhood out to a distance of $\approx 400$ pc with a distance resolution
of 1 pc. Leike et al. (2020) overcome the “fingers of god” effect by directly
modeling the dust density in 3D x-y-z cartesian space. Their method combines
metric Gaussian variational inference (Knollmüller & Enßlin, 2019) with
Gaussian Processes (e.g. Rasmussen & Williams, 2006) in the context of
information field theory (see Enßlin et al., 2009; Enßlin, 2019; Arras et al.,
2019; Enßlin & Frommert, 2011). The integrated extinction for a star known to
lie in a given 3D voxel with some uncertainty provides information on the
integrated dust density along that given line of sight, acting as an “anchor”
point that can constrain the smooth 3D distribution of dust. By inferring the
spatial correlation of the dust on different scales, one can connect these
anchor points, and infer the dust distribution at arbitrary points in space,
including along lines of sight where no stars are observed. Specifically,
Leike et al. (2020) model the dust as a log-normal Gaussian process, while
simultaneously inferring the correlation kernel of the process, corresponding
to the physical spatial correlation power spectrum of the dust. In this way,
Leike et al. (2020) use information on the distance and extinction to 5
million stars in the vicinity of the Sun (obtained from the Starhorse
algorithm; see Anders et al., 2019) to reconstruct the 3D spatial distribution
of dust on a fixed cartesian grid with a spatial resolution of 1 pc. The Leike
et al. (2020) results offer the first clear 3D spatial view of the local
interstellar medium.111The methodology is similar to that summarized in Rezaei
Kh. et al. (2017, 2018) and implemented in Rezaei Kh. et al. (2020) to obtain
a 3D cartesian view of the Orion A molecular cloud at a spatial resolution of
approximately 10 pc. We will compare the Leike et al. (2020) results to the
Rezaei Kh. et al. (2020) results in Orion A in §5.

In this work, we utilize the high spatial resolution of the Leike et al.
(2020) 3D dust map to analyze the 3D spatial structure and thicknesses of
famous nearby star forming regions for the first time. Specifically, we extend
two existing algorithms geared towards 2D data — FilFinder (Koch & Rosolowsky,
2015) and RadFil (Zucker & Chen, 2018) — into 3D. We apply the FilFinder3D
algorithm to determine the topological skeletons of nearby clouds. We then use
RadFil3D to determine their average radial volume density profiles, and
characterize the shapes of these profiles via both Gaussian and Plummer
fitting. We then validate our results by creating synthetic 2D dust extinction
maps from the 3D density distributions, and comparing these maps to
independent results (derived from the NICEST algorithm; Lombardi, 2009;
Lombardi et al., 2011) often used in star formation studies of nearby
molecular clouds.

In §2, we present the datasets used in this work, including the Leike et al.
(2020) dust map and the NICEST 2D extinction maps. In §3, we detail the
methodology we use to derive the topological spines of the clouds, and to
construct radial volume density profiles for each region. We also discuss the
derivation of our 2D extinction maps. In §4, we present our results in the
context of a custom 3D visualization environment, including the x-y-z-l-b-d
values (the Heliocentric Galactic Cartesian coordinates, the Galactic
coordinates, and the distance) for each point in the spine of every cloud, and
their thicknesses derived from the radial profile fitting. We also compare the
2D dust extinction maps derived from the projected Leike et al. (2020) 3D dust
distribution with NICEST 2D extinction maps, compute cloud masses from each
method, and quantify the dynamic range in density that the Leike et al. (2020)
dust map is able to probe. In §6, we discuss the implications of our results,
including potential evidence of a phase transition in our radial density
profiles. Finally, we conclude in §7.

##  2 Data

###  2.1 Leike et al. (2020) 3D Dust Map

Leike et al. (2020) compute the 3D dust distribution on a cartesian grid
centered on the Sun, extending $x$ = [-370 pc, 370 pc], $y$ = [-370 pc, 370
pc], and $z$ = [-270 pc, 270 pc], where $x$, $y$, and $z$ are the Heliocentric
Galactic Cartesian coordinates. As such, the distance range of the map varies
from $d=370$ pc in the midplane, up to $d\approx$ 590 pc at the corners of the
grid. This distance range encompasses about a dozen well-studied nearby star-
forming regions, and in this work we target the Chamaeleon, Ophiuchus, Orion
B, Orion A, $\lambda$ Orionis, Lupus, Taurus, Perseus, Musca, Pipe, and
Cepheus molecular clouds for our radial profile analyses. We also provide the
topological skeleton for one additional cloud — Corona Australis — but do not
fit for its width, since it does not meet the density threshold required for
inclusion in our study (see §3.1). We emphasize that three clouds – Orion A,
Orion B, and $\lambda$ Orionis – lie at the very edge of the Leike et al.
(2020) grid, and the results for these clouds are subject to additional
scrutiny, as discussed further in §4.

####  2.1.1 Converting to Gas Volume Density

The Leike et al. (2020) 3D dust map is provided in terms of the differential
extinction, specifically the optical depth in the Gaia G-band, $\tau_{G}$, per
one parsec. Leike et al. (2020) generate twelves samples of their 3D dust map,
and we adopt the mean of the samples in our main analysis. To ensure that
uncertainties in the underlying dust reconstruction are not affecting our
results, we will additionally repeat our radial profile analysis using all
twelve realizations of the dust map, as discussed further in the
Appendix.222For more details, see the public release of the Leike et al.
(2020) map available on Zenodo (Leike et al., 2020). The mean dust map is
available in the mean_std.h5 file, while the twelve samples are available in
samples.h5.

There are several steps needed to convert the data values of the mean map —
referred to as the extinction density $s_{x}$ in Leike et al. (2020) — to gas
volume density. First, we use the relationship between
$\frac{A_{\lambda}}{N_{\rm H}}$ (the ratio of the extinction $A$ at some
wavelength $\lambda$ over the hydrogen column density $N_{\rm H}$) and inverse
wavelength from Draine (2009, see Figure 1). Assuming a mean wavelength for
the Gaia G-band of 673 nm ($\lambda^{-1}=1.49\;\micron^{-1}$; Jordi et al.,
2010), we obtain $\frac{A_{G}}{N_{\rm H}}=4\times 10^{-22}\rm\;cm^{2}\;mag$
from Draine (2009).333Draine (2009) is assuming an $R_{V}$ = 3.1. Thus,
variation in $R_{V}$ could induce changes in our adopted conversion to gas
volume density. Given the known relationship between extinction $A_{\rm G}$
and optical depth ($A_{G}=1.086\times\tau_{G}$) we obtain $N_{\rm
H}=2.71\times 10^{21}\;\tau_{G}\;\rm cm^{-2}$. Finally, given that the Leike
et al. (2020) map is provided in the form of $\frac{\tau_{G}}{1\;\rm{pc}}$
(rather than the usual $\tau_{G}$), we obtain $n_{\rm H}=880\;{\rm
cm^{-3}}\times s_{x}$. Thus, to convert the Leike et al. (2020) 3D dust map to
volume density of hydrogen, we simply multiply the entire volume by a factor
of 880. We emphasize that this volume density is the total volume density of
hydrogen nuclei, independent of phase ($n=n_{\rm HI}+2n_{\rm H_{2}}$) and is
used to derive all results reported in §4. Volume density cubes for each
cloud, with WCS information in a cartesian reference frame, are available on
the Harvard Dataverse (see https://doi.org/10.7910/DVN/IADP7W).

####  2.1.2 Systematic Effects

While Leike et al. (2020) provide their results on a 1 pc grid, they argue
that the functional resolution of the map is 2 pc, and only features which
manifest at spatial scales greater than 2 pc should be considered reliable.
Accordingly, for the main analysis, we exclude the inner 2 pc of each dust
cloud’s radial volume density profile when inferring their structure, while
reporting results including the inner 2 pc in the Appendix.

By performing tests on synthetic mock data, Leike & Enßlin (2019); Leike et
al. (2020) also find that there are a small number of outliers in their 3D
dust distribution, in which the differential dust extinction values in certain
voxels is significantly larger than predicted, to a degree that cannot be
captured by the reported uncertainties. While Leike & Enßlin (2019) and Leike
et al. (2020) use different underlying stellar data for their statistical
reconstruction, the methodology is similar, and we find the same high-
differential-extinction outliers in the newer Leike et al. (2020) map. In this
analysis, we will be averaging over many voxels and will find that our results
are insensitive to these outliers.

We emphasize that the 3D volume densities reported here are likely lower
limits. We observe inhomogeneities in column density maps of nearby molecular
clouds at spatial scales smaller than 2 pc, the effective spatial resolution
of the Leike et al. (2020) 3D dust map (see e.g. Zari et al., 2016). Thus,
much of the dense gas in these clouds is unresolved, not just due to the map’s
spatial resolution, but due to the inability of the 3D dust mapping technique
to probe high density gas, stemming from the lack of stellar data constraining
the 3D dust reconstruction in this regime.

Finally, one additional consideration is the effect that YSOs have on the 3D
dust reconstruction in Leike et al. (2020) and whether this would bias any 3D
volume density estimates. The stellar catalog underpinning Leike et al. (2020)
is based on the StarHorse algorithm, presented in Anders et al. (2019). Anders
et al. (2019) does not employ any pre-main sequence models when determining
stellar properties. Thus, as argued in Anders et al. (2019), any young stars
detected with Gaia and included in the StarHorse sample will have poor fits,
triggering a reliability flag. Since Leike et al. (2020) only use high-quality
stellar data with no questionable reliability flags set, there should be very
little contamination from YSOs in the 3D dust map.

###  2.2 2D NICEST Extinction Maps

To validate our 3D cloud structure results (see §3) and quantify any “missing
mass” in the 3D maps, we compute 2D maps of the $K$ band extinction, $A_{K}$,
using the NICEST algorithm (Lombardi, 2009, A Near-Infrared Color Excess
Method Tailored for Small-Scale Structures) for each cloud. Goodman et al.
(2009) showed that dust extinction estimates based on the near-infrared color
excess are relatively free from bias and provide more robust estimates of the
column density, at higher dynamic range, than gas-based approaches. NICEST,
and its predecessor NICER (Lombardi & Alves, 2001), rely on the fact that the
difference between the observed and intrinsic colors of a star (background to
a cloud) provide information on the cloud’s extinction along the line of sight
to that star. By measuring the unreddened colors of stars in a nearby, low
extinction control field, it is possible to model the color excess on a star-
by-star basis. Because there is a small dispersion in the intrinsic near-
infrared colors of stars (and photometric errors), the technique requires
averaging the extinction estimates for many individual stars in batches across
the cloud, and pixel-by-pixel extinction estimates are computed as spatial
averages via the adoption of a Gaussian smoothing kernel. Compared to NICER,
NICEST removes contamination by foreground stars, which by their nature do not
probe the cloud’s extinction, and better characterizes small-scale
inhomogeneities (i.e. filaments) present in the dust.

Like Lombardi et al. (2011), to compute the extinction, we leverage the near-
infrared $J$, $H$, and $K$ band colors of millions of stars detected in 2MASS
(Skrutskie et al., 2006) for each field. Using the 2MASS photometry, NICEST
creates two sets of stellar samples: main sequence stars and YSOs, defined
using their respective K-band magnitude distributions. In practice, given a
star, the algorithm computes the intrinsic K-band magnitude (i.e. removes the
extinction) and compares this with the K-band number counts of field stars
(typically an exponential distribution) and YSOs (typically a Gaussian-like
curve). This comparison is used to assign to each star a probability of being
a YSO. The functionality to produce custom cloud-by-cloud NICEST extinction
maps is available online (Lombardi & Alves, 2001; Lombardi, 2009)444See
http://interstellarclouds.fisica.unimi.it/html/index.html and we adopt the
default map making and cluster parameters. The default option is to include
all stars in the extinction calculation, including YSOs. To check that YSOs
are not affecting our mass calculations, we compute two masses for Orion A
both with and without including YSOs in the sample. We find that including the
YSOs in the extinction calculation only increases the mass by 2%. Since Orion
A contains the largest number of YSOs in the solar neighborhood, our decision
to compute the extinction maps using all stars should have no effect on our
results. The configuration files used to produce the NICEST maps are available
online at the Harvard Dataverse; the NICEST maps for all clouds targeted in
this study are likewise available online (see
https://doi.org/10.7910/DVN/WONVUH for both).

##  3 Methods

###  3.1 FilFinder 3D

To determine the topological skeletons of nearby molecular clouds, we employ a
modified 3D version of the FilFinder package from Koch & Rosolowsky (2015).
FilFinder was originally developed in 2D to constrain the filamentary
structure in $\rm H_{2}$ column density maps obtained as part of the Herschel
Gould Belt Survey (André et al., 2010). FilFinder3D uses similar morphological
operations to the existing 2D version, but now employs the skan package for
efficient operations on large morphological skeletons (Nunez-Iglesias et al.,
2018). A full description of the 3D algorithm will be presented in Koch et al.
(2021, in preparation).

To compute the skeletons, we start by defining the same density threshold for
all clouds, equivalent to $n=35\;\rm cm^{-3}$, based on the volume density
distribution computed in §2.1.555Corona Australis has no significant extended
density features above a level of $n=35\;\rm cm^{-3}$, so we skeletonize it at
a level of $n=5\;\rm cm^{-3}$ to provide distance information. It is included
in Table 1 but not in Tables 2 and 3. Hierarchical structure finding
algorithms for feature identification do exist for 3D volumes, the most common
being dendrograms (Rosolowsky et al., 2008), which are frequently applied to
spectral-line cubes of CO emission in position-position-velocity space (Rice
et al., 2016). However, it is well known that dendrograms are not able to
capture filamentary structure, with rarely any structures exceeding an aspect
ratio of $2:1$ (see Rice et al., 2016). Since we find that nearby star-forming
regions are indeed filamentary, even at lower densities, we forgo dendrograms
for this analysis and defer to future work. Before settling on a fixed
threshold of $n=\rm 35\;cm\;^{-3}$, we test a variety of thresholds, ranging
from $20-50\;\rm cm^{-3}$. Thresholds significantly higher than $n=\rm
35\;cm\;^{-3}$ were too patchy to skeletonize, while thresholds significantly
lower than $n=\rm 35\;cm\;^{-3}$ meant that the peak radial densities were
often offset from the main spine.666In previous analyses based on 2D dust
emission maps (Zucker et al., 2018a), we implemented a “shift” option in the
RadFil package to account for dense clumps offset from the main spine, since
FilFinder does not take into the density range inside each mask when defining
the spine. The “shift” option corrects for this by shifting each cut so that
its peak density is centered at a radial distance of 0 pc. However, we forgo
shifting in this work due to the small handful of noisy outliers at very high
differential extinction, which would bias the density results at small radii.
Our choice of density threshold is validated in §4, since the spines defined
in 3D tend to align very well with the spines of filaments seen in 2D dust
emission from Planck (Planck Collaboration et al., 2011), and the radial
density profiles peak near a radial distance of zero. Overall we find that all
clouds are well-defined at $n=\rm 35\;cm\;^{-3}$ and we can produce the most
unbiased widths possible by not tweaking feature identification on a cloud-by-
cloud basis.

After defining this initial $n=\rm 35\;cm\;^{-3}$ threshold, we
morphologically close the 3D mask by performing a dilation, followed by an
erosion. This procedure removes small voids in the mask, rendering it suitable
for skeletonization.777More details on mathematical morphology, including the
erosion, dilation, and closing of features in 2D or 3D space, can be found in
the scikit-image documentation, which we use to perform these operations (van
der Walt et al., 2014). We then skeletonize the 3D mask using the Medial Axis
Transform, convert the skeleton to a graph, and run a longest path algorithm
on the resulting skeleton to obtain the main trunk of each density feature.
While we only use the spine computed from the longest path algorithm to derive
the widths, we also preserve and report the full skeleton without pruning
applied, in order to retain distance information on all branches. In this
scheme, nearby molecular clouds can either be defined by a single, or
multiple, trunks, depending on whether they are connected in density space
above a threshold of $n=35\;\rm cm^{-3}$. Thus, this skeletonization process
is performed between one and four times, depending on how many isolated 3D
features there are after performing the initial thresholding. For clouds like
Perseus, which manifest as a single feature in density space, we identify only
a single spine. However, clouds like Taurus are actually composed of several
distinct features in density space. The skeletonization and longest path
algorithms are individually run on each isolated mask, to produce a set of
trunks defining the cloud.888The two exceptions are the Musca filament and the
Pipe nebula. We find that Musca is connected to the Chamaeleon cloud in 3D
density space at a level of $n=35\;\rm cm^{-3}$. The same is true for the Pipe
nebula, which is connected to Ophiuchus. Because these clouds have always been
treated as individual entities in the literature, and obtaining widths for the
Pipe nebula and Musca (as opposed to the Pipe-Ophiuchus, and Musca-Chamaeleon
complexes) is beneficial for comparison with prior work, we mask out the Pipe
nebula when obtaining widths for Ophiuchus, and the Musca cloud when obtaining
widths for Chamaeleon (and vice versa). In reality, the complexes are
connected by only very thin tendrils, so this procedure does not affect the
morphology of the derived masks, or the resulting widths.

Finally, we emphasize again that the FilFinder algorithm is designed to work
on elongated, filamentary structures, with aspect ratios of at least $3:1$. We
find that the vast majority of our masks tend to be elongated with aspect
ratios of at least $3:1$, but there are a few sub-features in clouds which are
better described as ellipsoids than filaments, some of which may be unresolved
in the Leike et al. (2020) map. Ellipsoids do not dominate the set of density
features, either across the sample, or within an individual cloud, so this has
a negligible effect on the results. The volume density distributions
(https://doi.org/10.7910/DVN/IADP7W), 3D masks
(https://doi.org/10.7910/DVN/NMRLP1), and 3D spines
(https://doi.org/10.7910/DVN/YQYBRD) are available online with full WCS
information at the Harvard Dataverse. We refer interested readers there to
explore our masking and skeletonization results on a cloud-by-cloud basis in
more detail.

###  3.2 RadFil 3D

####  3.2.1 Building Profiles

To characterize the thicknesses of local clouds, we employ the radial profile
builder RadFil (Zucker & Chen, 2018). RadFil was designed in 2D to build
radial column density profiles using filament spines determined by the
FilFinder package (Zucker & Chen, 2018; Koch & Rosolowsky, 2015; Zucker et
al., 2018a). RadFil builds profiles by taking radial cuts across the spine.
The same logic for constructing a radial profile in 2D can easily be extended
into 3D. We take the 3D skeleton returned by FilFinder3D, and fit a B-spline
to obtain a smoothed continuous version of the spine and the first derivative
along its length. Using functionality within the Python package PyVista
(Sullivan & Kaszynski, 2019), we generate a slice through the volume density
cube at each point along the spine given the normal vector obtained from the
first derivative. Then, for each point in the slice, we determine its radial
distance from the spine point and its volume density via linear interpolation.
On a slice by slice basis, we compute the median density in forty radial
distance bins, spaced half a parsec apart, out to 20 pc. Given the full set of
profiles, an average profile is obtained by computing the median volume
density in each radial distance bin. Before taking the average, we remove the
small minority of profiles in which the peak radial density is less than the
density threshold used to define the mask ($n=35\;\rm cm^{-3}$). For clouds
composed of multiple isolated density features, the average is taken using
slices across all density features, to compute a representative width for each
cloud.

This averaged radial profile is used to perform the Gaussian and Plummer
fitting described in §3.2.2. Because of the 1 pc voxel size of the 3D dust
map, and the fact that Leike et al. (2020) recommend that scales less than 2
pc are suspect, we remove the inner 2 pc of each averaged profile when fitting
the profiles. However, in the Appendix, we report results for the full profile
fits (including the inner 2 pc) and find that our choice to remove radial
distances between 0 and 2 pc does not significantly affect our results.

####  3.2.2 Fitting Profiles

We fit three functions to the averaged radial volume density profile computed
for each cloud in §3.2.1 — a single-component Gaussian, a two-component
Gaussian, and a Plummer function. One important consideration is the modeling
of the scatter in the averaged profile. In our case, the scatter in the volume
density within a bin at fixed radial distance is not a true measurement error,
but rather encompasses real physical variation in the radial density profile
taken at different slices along the spine of a cloud. For this reason, we
choose to model the scatter on the density as part of our fitting procedure,
introducing it as an additional free parameter in all three of our models. We
perform our inference in a Bayesian framework, adopting a simple Gaussian
likelihood (see e.g. Hogg et al., 2010) and flat priors on all our parameters.
None of our parameters are prior-dominated, so this choice of flat prior has
no effect on our results. For computational expediency, we minimize the
negative log-likelihood, in lieu of maximizing the likelihood. Our Gaussian
log-likelihood is of the form:

$\ln\,p(n\big{|}\theta,\sigma_{n}^{2})=-\frac{1}{2}\sum_{i}\left[\frac{[n_{i}-n_{\theta,i}]^{2}}{\sigma_{n}^{2}}+\ln(2\pi\sigma_{n}^{2})\right]$
(1)

where $n_{i}$ is the measured density in the $i$th radial distance bin,
$n_{\theta,i}$ is the modeled density in the $i$th radial distance bin
(defined by the free parameters $\theta$), and $\sigma_{n}^{2}$ is the square
of the scatter in the density across all bins that we also infer as part of
our modeling.

For the single-component Gaussian fit, our model for the volume density of
hydrogen nuclei $n$ as a function of radial distance $r$ is of the form:

$n(r)=a\,\exp{\left(\frac{{-r}^{2}}{2\sigma^{2}}\right)}$ (2)

Our free parameters $\theta$ are the amplitude ($a$) and standard deviation
($\sigma$). The mean of the Gaussian is fixed to zero. We adopt a flat prior
between $0-80$ $\rm cm^{-3}$ for $a$, between $0-15$ pc for $\sigma$ and
between $0-100$ $\rm cm^{-3}$ for $\sigma_{n}^{2}$.

To better model the significant tails of the distribution, we additionally fit
a two-component Gaussian, again with the mean of each component fixed to zero:

$n(r)=a_{1}\,\exp{\left(\frac{{-r}^{2}}{2\sigma_{1}^{2}}\right)}+a_{2}\,\exp{\left(\frac{{-r}^{2}}{2\sigma_{2}^{2}}\right)}$
(3)

where the amplitudes ($a_{1},a_{2}$) and their standard deviations
($\sigma_{1},\sigma_{2}$) are the free parameters in the modeling. We adopt a
flat prior between $20-80$ $\rm cm^{-3}$ for $a_{1}$, between $0-6$ pc for
$\sigma_{1}$, between $2-20$ $\rm cm^{-3}$ for $a_{2}$, between $6-25$ pc for
$\sigma_{2}$, and between $0-100$ $\rm cm^{-3}$ for $\sigma_{n}^{2}$. The
offset in the prior range between $a_{1}$ and $a_{2}$ is adopted to avoid
degeneracies in the modeling.

Finally, we fit a Plummer function, which we parameterize as:

$n(r)=\frac{n_{0}}{[{1+({\frac{r}{R_{\mathrm{flat}}})^{2}}]}^{\;\frac{p}{2}}}$
(4)

where our free parameters are $n_{0}$ (peak profile height),
$R_{\mathrm{flat}}$ (the flattening radius), and $p$ (index of the density
profile). Equation 4 is the same ubiquitous Plummer function applied to column
density maps of smaller-scale filaments (widths a few tenths of a parsec,
lengths a few parsecs; Arzoumanian et al., 2019; Panopoulou et al., 2017) from
the Herschel Gould Belt Survey (André et al., 2010), except we have
transformed the model from column density space back to volume density space
(see Equation 1 of Arzoumanian et al., 2011). We adopt a flat prior between
$20-170$ $\rm cm^{-3}$ for $n_{0}$, between $0-15$ pc for $R_{\mathrm{flat}}$,
between $0-6$ for $p$, and between $0-100$ $\rm cm^{-3}$ for $\sigma_{n}^{2}$.

We sample for all model parameters using the nested sampling code dynesty
(Speagle, 2020). We adopt the default parameters of the nested sampler, with a
stopping criterion of dlogz = 0.1. We utilize the logarithm of the evidence
$\ln(Z)$ extracted from the dynesty chains to compute a Bayes Factor for the
purpose of model comparison, as discussed further in §4.

###  3.3 Creating 2D Dust Extinction Maps from 3D Dust Data

To compare the Leike et al. (2020) 3D dust map with traditional 2D approaches
(see §2.2), we create projected 2D dust extinction maps from the 3D volume
density distributions used in the construction of the radial profiles. To do
so, we utilize the FitsOffAxisProjection functionality from the Python package
yt (Turk et al., 2011). The FITSOffAxisProjection function integrates the
volume density $n$ along a line of sight $\hat{l}$, where $\hat{l}$ is
determined by the vector between the Sun and the median cartesian coordinates
of the points in a cloud’s skeleton, $(x_{\rm cen},y_{\rm cen},z_{\rm cen})$.
The resulting projection corresponds to the total hydrogen column density
distribution $N$ = $N_{\rm HI}+2N_{\rm H_{2}}$ (in $\rm{cm^{-2}}$) as seen
from the Sun. In order to compare with the NICEST maps of the $K$ band
extinction $A_{\rm K}$ (from §2.2), we convert from total hydrogen column
density to $A_{K}$ by assuming a constant conversion factor from Lada et al.
(2009) of $\frac{N({\rm HI})+2N({\rm H_{\rm 2}})}{A_{K}}=1.67\times
10^{22}\rm\;cm^{-2}\;mag^{-1}$.

To assign sky coordinates to each projection, we convert from the skeleton’s
central coordinates $(x_{\rm cen},y_{\rm cen},z_{\rm cen})$ to sky coordinates
$(l_{\rm cen},b_{\rm cen},d_{\rm cen})$, where $l_{\rm cen}$ and $b_{\rm cen}$
become the central Galactic Longitude and Latitude of the projected image on
the plane of the sky, and $d_{\rm cen}$ is the distance between the Sun and
the central point of the cloud’s skeleton. The pixel scale of the image is
computed by determining the angular size spanned by a 1 pc voxel (the voxel
size of the original 3D dust grid) at the central distance of the cloud’s
skeleton $d_{\rm cen}$. These projected $A_{K}$ maps are available for
download at the Harvard Dataverse (see https://doi.org/10.7910/DVN/GXXKHD).

###  3.4 Mass Calculation

To quantify the fraction of cloud mass being recovered with 3D dust mapping,
we compute cloud masses using both traditional NICEST 2D dust extinction maps
(from §2.2) and the projected $A_{K}$ maps we derive from the 3D dust in §3.3.
In both cases, we adopt a consistent threshold of $A_{K}$=0.1 mag. This
threshold is the same used in previous literature to calculate cloud masses
from the NICEST 2D dust extinction maps (see e.g. Lada et al., 2010).999Like
Lada et al. (2010), we also subtract a constant background of $A_{K}=0.15$ mag
from the Pipe nebula, due to its position close to the plane and towards the
Galactic center. Furthermore, because the 2D dust extinction maps based on
Leike et al. (2020) are computed using only dust local to the cloud, rather
than integrated along the line of sight, we limit our comparison to latitudes
$|b|>5^{\circ}$. In reality, this latitude cut only masks out a small section
of the Pipe and Lupus clouds. The masses are calculated inside the cloud
boundaries by assuming the same mass surface density relation as Zari et al.
(2016). This relation is equivalent to $\frac{\Sigma_{\rm gas}}{A_{K}}=\mu
m_{p}\beta_{K}=183\;M_{\sun}\;\rm pc^{-2}\;mag^{-1}$, where $\mu=1.37$ is the
mean molecular weight, $m_{p}=1.67\times 10^{-24}$ g is the mass of the
proton, and $\beta_{K}$ is the same gas-to-dust ratio that we adopt above,
equivalent to $1.67\times 10^{22}\rm\;cm^{-2}\;mag^{-1}$. To ensure a fair
comparison, the NICEST maps from §2.2 are convolved to the same resolution as
the Leike et al. (2020) based extinction maps. In both cases, the same
distance $d_{cen}$ — the average distance to a cloud’s skeleton — is adopted.

Table 1: Skeletonization Results

Cloud | l | b | $d_{min}$ | $d_{med}$ | $d_{max}$ | $x_{min}$ | $x_{max}$ | $y_{min}$ | $y_{max}$ | $z_{min}$ | $z_{max}$ | $N_{components}$ | $l_{skeleton}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
| ∘ | ∘ | pc | pc | pc | pc | pc | pc | pc | pc | pc |  | pc
Chamaeleon | 299.0 | -15.9 | 173 | 183 | 190 | 69 | 99 | -159 | -148 | -57 | -35 | 1 | 68
Ophiuchus | 355.4 | 17.0 | 115 | 136 | 150 | 104 | 144 | -26 | 5 | 24 | 48 | 3 | 74
Lupus | 341.7 | 5.1 | 155 | 165 | 198 | 139 | 196 | -73 | -28 | 4 | 47 | 3 | 80
Taurus | 172.6 | -15.6 | 131 | 148 | 168 | -166 | -122 | 6 | 44 | -45 | -26 | 3 | 79
Perseus | 158.9 | -20.2 | 279 | 295 | 301 | -270 | -236 | 95 | 100 | -116 | -81 | 1 | 57
Musca | 301.0 | -9.7 | 171 | 172 | 173 | 86 | 91 | -146 | -145 | -36 | -22 | 1 | 16
Pipe | 2.1 | 4.1 | 147 | 152 | 163 | 146 | 163 | -8 | 7 | 6 | 18 | 1 | 27
Cepheus | 112.1 | 15.9 | 337 | 345 | 370 | -133 | -65 | 300 | 323 | 81 | 133 | 4 | 52
Corona AustralisaaCorona Australis is too diffuse to meet the minimum density threshold to be included in the width analysis for this study. We skeletonize the cloud at a much lower threshold ($n=5\;\rm cm^{-3}$) to provide distance information, but do not include it in subsequent analysis in Tables 2 and 3. | 6.5 | -24.6 | 136 | 161 | 179 | 132 | 162 | -2 | 24 | -77 | -33 | 1 | 86
Orion A* | 212.5 | -19.6 | 391 | 405 | 445 | -349 | -315 | -241 | -177 | -149 | -118 | 3 | 52
Orion B* | 205.5 | -14.0 | 397 | 404 | 406 | -362 | -332 | -190 | -157 | -114 | -76 | 1 | 69
$\lambda$ Orionis* | 194.3 | -14.6 | 375 | 386 | 397 | -367 | -356 | -129 | -73 | -109 | -73 | 3 | 62

**footnotetext: Clouds below the horizontal line – Orion A, Orion B, and
$\lambda$ Orionis – should be treated with caution, as they lie at the very
edge of the Leike et al. (2020) 3D dust grid and are subjected to additional
biases. Part of Orion A, toward the Orion Nebula Cluster, is also devoid of
dust in the Leike et al. (2020) 3D dust map despite being highly extinguished.
See §5 for full details.

Note. — Properties of the 3D skeletons for local molecular clouds. (1) Name of
the cloud. (2-3) Median Galactic longitude $l$ and latitude $b$ of the points
defining the cloud’s skeleton. (4 - 6) Minimum, median, and maximum distance
of the points defining the skeleton. (7-8) Minimum and maximum extent of the
cloud’s skeleton in the Heliocentric Galactic cartesian $x$ direction. (9-10)
Minimum and maximum extent of the cloud’s skeleton in the Heliocentric
Galactic cartesian $y$ direction. (11-12) Minimum and maximum extent of the
cloud’s skeleton in the Heliocentric Galactic cartesian $z$ direction. (13)
Number of skeletons per cloud, determined by the number of significant
features above a density threshold of $n=\rm 35\;cm^{-3}$. (14) Total length
of the set of skeletons in (13), computed using only each skeleton’s main
trunk via FilFinder’s longest path algorithm. At the Harvard Dataverse, we
provide tables in FITS format with the ($x,y,z,l,b,d$) values of each point in
the spine (for both the longest path skeleton summarized here, and the full,
unpruned skeleton). We also provide a machine readable version of this table
(see https://doi.org/10.7910/DVN/CAVMAQ). Interactive figures highlighting
these skeletonization results are available at
https://faun.rc.fas.harvard.edu/czucker/Paper_Figures/3D_Cloud_Topologies/gallery.html.

Figure 1: Skeletonization results for the Perseus Molecular Cloud shown in 2D
(left panel) and 3D (right panel). In both panels, the scatter points (colored
by distance over the range $d=279-301$ pc) delineate the skeleton of the
cloud, determined in 3D Cartesian space and projected back on the plane of the
sky. In the 3D version of the figure we show the 3D volume density of the
cloud $n$, where the cyan and white distributions indicate lower thresholds of
$n=\rm 4\;cm^{-3}$ and $n=\rm 35\;cm^{-3}$ respectively. In the plane of the
sky version, we show the Planck reddening map (Planck Collaboration et al.,
2011). An interactive version of this figure is available at
https://faun.rc.fas.harvard.edu/czucker/Paper_Figures/3D_Cloud_Topologies/perseus_topology/perseus.html
or in the online version of the published article. Comparable interactive
figures for all clouds in the sample are available at
https://faun.rc.fas.harvard.edu/czucker/Paper_Figures/3D_Cloud_Topologies/gallery.html.
In the online version, the cyan and white boundaries (toggleable on/off)
correspond to the FWHM of the inner and outer envelopes of the cloud, derived
from the two-component Gaussian fitting results in §4.  Figure 2: Same as in
Figure 1 but for the Chamaeleon Molecular Cloud, with the spine (colored
points) spanning a distance range of $d=173-190$ pc. Opaque spine points
constitute the skeleton’s main spine, while semi-transparent points are part
of the full, unpruned skeleton. An interactive version of this figure is
available at
https://faun.rc.fas.harvard.edu/czucker/Paper_Figures/3D_Cloud_Topologies/chamaeleon_topology/chamaeleon.html
or in the online version of the published article. Comparable interactive
figures for all clouds in the sample are available at
https://faun.rc.fas.harvard.edu/czucker/Paper_Figures/3D_Cloud_Topologies/gallery.html.

##  4 Results

###  4.1 Skeletonization

In Table 1 we summarize the properties of the skeletons, determined using
FilFinder3D’s longest path algorithm, for all clouds in our sample. This
includes the minimum, median, and maximum distance to each cloud’s skeleton,
along with the extent of the cloud’s skeleton along the Heliocentric Galactic
cartesian $x$, $y$, and $z$ directions. We also include the number of skeletal
components (skeletonized at a density threshold above $n=\rm 35\;cm^{-3}$) and
the total length of these components. The length is the total length of denser
features (computed using only the skeleton returned by the longest path
algorithm) so it does not include features branching off the main trunk. The
properties of both the longest path skeletons, and the full unpruned
skeletons, are available online in machine readable format on the Harvard
Dataverse (see https://doi.org/10.7910/DVN/YQYBRD). The skeletonization
results for two clouds (Perseus and Chamaeleon) are shown in Figures 1 and 2,
respectively. These figures are interactive in their online version; it is
possible to manipulate the skeleton in 3D, or hover over any point in the 2D
skeleton to obtain its plane-of-the-sky coordinates and distance. Comparable
interactive figures for the rest of the clouds in the sample are available at
https://faun.rc.fas.harvard.edu/czucker/Paper_Figures/3D_Cloud_Topologies/gallery.html
or in the online version of the published article.

We have separated Table 1 into two parts, representing “complete” (top rows)
and “incomplete” (bottom rows) clouds. Incomplete clouds – Orion A, Orion B,
and $\lambda$ Orionis — are those which lie at the very edge of the Leike et
al. (2020) 3D dust map, which extends to 370 pc along the cartesian $x$ and
$y$ directions. While these three clouds are still capable of being
skeletonized, it is likely that some fraction of each cloud lies beyond the
map boundary, which would bias the skeletons to be mildly closer in distance
than they actually are. This incompleteness will be further validated in §4.4,
where projected 2D extinction maps based on the Leike et al. (2020) 3D dust
distribution indicate that Orion A, Orion B, and $\lambda$ Orionis are the
only three clouds in the sample whose masses are significantly discrepant from
counterpart masses calculated out to “infinity” (that is, using the NICEST
algorithm; Lombardi, 2009). This discrepancy indicates that there is missing
mass present beyond the boundaries of the map.

=1.5in

Table 2: Radial Volume Density Fitting Results

| Two-Component Gaussian | Plummer | Single Gaussian
---|---|---|---
Cloud | $\sigma_{1}$ | $a_{1}$ | $\sigma_{2}$ | $a_{2}$ | $\sigma_{n,\rm G2}^{2}$ | $\frac{\sigma_{2}}{\sigma_{1}}$ | $lnZ_{\rm G2}$ | $n_{0}$ | $R_{\rm flat}$ | $p$ | $\sigma_{n,\rm P}^{2}$ | $lnZ_{P}$ | $\sigma$ | $a$ | $\sigma_{n,G1}^{2}$ | $lnZ_{\rm G1}$
| pc | cm-3 | pc | cm-3 | cm-3 |  |  | cm-3 | pc |  | cm-3 |  | pc | cm-3 | cm-3 |
Chamaeleon | $2.9_{-0.1}^{+0.1}$ | $31.3_{-0.6}^{+0.5}$ | $10.2_{-0.3}^{+0.4}$ | $8.7_{-0.5}^{+0.5}$ | $0.2_{-0.0}^{+0.0}$ | 3.5 | -2.9 | $47.8_{-1.2}^{+1.5}$ | $3.1_{-0.2}^{+0.1}$ | $1.8_{-0.1}^{+0.1}$ | $0.1_{-0.0}^{+0.0}$ | -1.5 | $4.8_{-0.3}^{+0.3}$ | $30.8_{-2.0}^{+2.1}$ | $6.1_{-1.3}^{+1.9}$ | -58.0
Ophiuchus | $2.7_{-0.2}^{+0.2}$ | $31.1_{-1.5}^{+1.6}$ | $10.1_{-0.5}^{+0.6}$ | $12.6_{-1.1}^{+1.1}$ | $0.8_{-0.2}^{+0.2}$ | 3.8 | -27.2 | $56.0_{-2.9}^{+3.5}$ | $2.3_{-0.2}^{+0.2}$ | $1.4_{-0.0}^{+0.0}$ | $0.3_{-0.1}^{+0.1}$ | -11.9 | $5.8_{-0.4}^{+0.5}$ | $29.1_{-2.0}^{+2.0}$ | $9.5_{-2.0}^{+2.6}$ | -64.7
Lupus | $2.6_{-0.2}^{+0.3}$ | $30.6_{-1.8}^{+1.8}$ | $8.1_{-0.5}^{+0.6}$ | $15.5_{-2.0}^{+1.9}$ | $0.9_{-0.2}^{+0.3}$ | 3.1 | -29.2 | $48.6_{-1.9}^{+2.3}$ | $3.6_{-0.3}^{+0.3}$ | $1.9_{-0.1}^{+0.1}$ | $0.5_{-0.1}^{+0.1}$ | -19.3 | $5.2_{-0.3}^{+0.3}$ | $32.7_{-1.7}^{+1.8}$ | $6.2_{-1.3}^{+1.7}$ | -58.1
Taurus | $2.6_{-0.1}^{+0.1}$ | $30.9_{-1.1}^{+1.2}$ | $11.9_{-0.5}^{+0.6}$ | $9.2_{-0.5}^{+0.5}$ | $0.4_{-0.1}^{+0.1}$ | 4.6 | -17.4 | $72.8_{-11.1}^{+22.7}$ | $1.2_{-0.3}^{+0.3}$ | $1.2_{-0.0}^{+0.0}$ | $0.4_{-0.1}^{+0.1}$ | -17.3 | $5.8_{-0.6}^{+0.6}$ | $24.3_{-2.2}^{+2.4}$ | $10.3_{-2.3}^{+3.3}$ | -65.7
Perseus | $4.1_{-0.1}^{+0.1}$ | $34.3_{-0.9}^{+0.9}$ | $11.9_{-0.6}^{+0.8}$ | $10.9_{-1.0}^{+1.0}$ | $0.3_{-0.1}^{+0.1}$ | 2.9 | -15.0 | $47.8_{-0.9}^{+0.9}$ | $6.1_{-0.3}^{+0.3}$ | $2.4_{-0.1}^{+0.1}$ | $0.5_{-0.1}^{+0.1}$ | -20.7 | $6.0_{-0.2}^{+0.2}$ | $39.0_{-1.5}^{+1.6}$ | $6.9_{-1.4}^{+2.1}$ | -60.4
MuscaaaAt small radial distances, Musca is largely unresolved in the Leike et al. (2020) map, and at large radial distances, its structure is consistent with being prior-dominated, so we urge caution when interpreting the results. | $3.1_{-0.1}^{+0.1}$ | $31.8_{-1.5}^{+1.5}$ | $17.9_{-3.0}^{+4.1}$ | $4.9_{-0.6}^{+0.7}$ | $1.3_{-0.3}^{+0.4}$ | 5.8 | -34.4 | $48.3_{-5.5}^{+8.4}$ | $2.5_{-0.5}^{+0.6}$ | $1.6_{-0.2}^{+0.2}$ | $2.4_{-0.5}^{+0.7}$ | -44.5 | $4.4_{-0.4}^{+0.4}$ | $29.8_{-2.9}^{+3.0}$ | $8.9_{-2.0}^{+2.7}$ | -63.8
Pipe | $4.1_{-0.4}^{+0.4}$ | $25.5_{-2.8}^{+2.7}$ | $9.7_{-1.2}^{+2.5}$ | $9.0_{-3.4}^{+3.3}$ | $1.5_{-0.3}^{+0.5}$ | 2.4 | -36.7 | $35.5_{-1.4}^{+1.5}$ | $7.7_{-1.0}^{+1.2}$ | $3.4_{-0.5}^{+0.6}$ | $1.4_{-0.3}^{+0.4}$ | -35.5 | $5.7_{-0.2}^{+0.2}$ | $30.6_{-1.1}^{+1.2}$ | $3.5_{-0.7}^{+1.1}$ | -49.3
Cepheus | $2.5_{-0.1}^{+0.1}$ | $46.6_{-1.6}^{+1.7}$ | $10.9_{-0.6}^{+0.7}$ | $10.1_{-0.7}^{+0.7}$ | $0.7_{-0.2}^{+0.2}$ | 4.4 | -26.4 | $116.9_{-17.5}^{+23.7}$ | $1.3_{-0.2}^{+0.2}$ | $1.4_{-0.0}^{+0.0}$ | $0.5_{-0.1}^{+0.2}$ | -21.6 | $4.3_{-0.4}^{+0.4}$ | $39.1_{-3.7}^{+3.9}$ | $14.3_{-3.0}^{+4.2}$ | -71.4
Orion A* | $2.9_{-0.1}^{+0.1}$ | $52.4_{-1.8}^{+1.8}$ | $8.4_{-0.8}^{+1.0}$ | $8.8_{-1.7}^{+1.9}$ | $1.3_{-0.3}^{+0.4}$ | 2.9 | -35.0 | $68.3_{-3.0}^{+3.4}$ | $4.4_{-0.4}^{+0.5}$ | $3.1_{-0.3}^{+0.3}$ | $1.5_{-0.3}^{+0.5}$ | -37.2 | $3.7_{-0.1}^{+0.1}$ | $53.1_{-2.4}^{+2.5}$ | $5.8_{-1.3}^{+1.8}$ | -57.5
Orion B* | $4.9_{-0.1}^{+0.1}$ | $39.2_{-1.1}^{+1.1}$ | $14.6_{-1.3}^{+1.9}$ | $9.0_{-1.2}^{+1.3}$ | $0.3_{-0.1}^{+0.1}$ | 3.0 | -13.6 | $50.3_{-0.8}^{+0.8}$ | $8.1_{-0.4}^{+0.5}$ | $3.0_{-0.2}^{+0.2}$ | $0.6_{-0.1}^{+0.2}$ | -23.6 | $6.5_{-0.2}^{+0.2}$ | $43.8_{-1.3}^{+1.4}$ | $5.9_{-1.2}^{+1.8}$ | -58.0
$\lambda$ Orionis* | $3.7_{-0.1}^{+0.1}$ | $27.2_{-0.5}^{+0.5}$ | $12.8_{-0.4}^{+0.5}$ | $9.3_{-0.4}^{+0.4}$ | $0.1_{-0.0}^{+0.0}$ | 3.4 | -0.7 | $40.2_{-1.0}^{+1.1}$ | $4.3_{-0.2}^{+0.3}$ | $1.7_{-0.1}^{+0.1}$ | $0.3_{-0.1}^{+0.1}$ | -14.0 | $6.3_{-0.4}^{+0.3}$ | $29.1_{-1.6}^{+1.6}$ | $7.2_{-1.5}^{+2.0}$ | -60.2

**footnotetext: Clouds below the horizontal line – Orion A, Orion B, and
$\lambda$ Orionis – should be treated with caution, as they lie at the very
edge of the Leike et al. (2020) 3D dust grid and are subjected to additional
biases. Part of Orion A, toward the Orion Nebula Cluster, is also devoid of
dust in the Leike et al. (2020) 3D dust map despite being highly extinguished.
See §5 for full details.

Note. — Cloud thickness results, obtained via radial volume density profile
fitting of each cloud’s 3D density distribution. (1) Name of the cloud. (2 -
6) Best-fit parameters for the two-component Gaussian fit, including the
standard deviation and amplitude of the narrower Gaussian ($\sigma_{1}$,
$a_{1}$), the standard deviation and amplitude of the broader Gaussian
($\sigma_{2}$, $a_{2}$), and the square of the modeled scatter in the density
measurements ($\sigma_{n,\rm G2}^{2}$). (7) Ratio of the outer to inner
Gaussian standard deviations. (8) Logarithm of the evidence of the two-
component Gaussian fit. (9-12) Best-fit parameters of the Plummer fit,
including the peak amplitude $n_{0}$, the flattening radius $R_{\rm flat}$,
the power index of the density profile $p$, and the square of the modeled
scatter in the density measurements ($\sigma_{n,\rm P}^{2}$). (13) Logarithm
of the evidence of the Plummer fit. (14-16) Standard deviation $\sigma$,
amplitude $a$, and the square of the modeled scatter in the density
measurements ($\sigma_{n,\rm G1}^{2}$) for the single component Gaussian fit.
(17) Logarithm of the evidence of the single-component Gaussian fit. We report
the single-component Gaussian results here for completeness only, as it is a
poor fit to the data particularly at large radial distances. A machine
readable version of this table is available at
https://doi.org/10.7910/DVN/QKYR3G.

Figure 3: Radial volume density profile for the Perseus Molecular Cloud (top)
and the Chamaeleon Molecular Cloud (bottom). Because the Leike et al. (2020)
3D dust map is more reliable at scales $>2$ pc and is not capable of probing
the highest gas densities, we exclude the inner 2 pc (coincident with the
yellow bar) from the fit. In each panel the gray dots indicate the profiles
for individual slices along the spine, while the blue dots indicate the median
of all slices. We fit three functions to the averaged radial volume density
profile, whose fits are shown in solid red: a two-component Gaussian function
(left), a Plummer function (middle), and a single-component Gaussian function
(right). In the left panel, we additionally overlay the decomposition of the
two Gaussians via the dashed white and cyan lines, corresponding to the
“inner” (black-white dotted) and “outer” (black-cyan dotted) components.
Identical plots for the rest of the sample are available in the Appendix.
Figure 4: Summary of the radial volume density profile fitting results for the
two-component Gaussian function, including the standard deviation $\sigma_{1}$
of the inner Gaussian, the standard deviation $\sigma_{2}$ of the outer
Gaussian, the amplitude of the inner Gaussian $a_{1}$, and the amplitude of
the outer Gaussian $a_{2}$. Each point represents a single local cloud in the
sample.  Figure 5: Summary of the radial volume density profile fitting
results for the Plummer function, including the distribution of the flattening
radius $R_{\mathrm{flat}}$, the peak density $n_{0}$, and the index of the
density profile $p$. Each point represents a single local cloud in the sample.
We observe a strong degeneracy in the Plummer modeling, particularly between
the $R_{\mathrm{flat}}$ and $p$ parameters.

###  4.2 Width Results

In Table 2, we summarize the results of our radial volume density fitting
procedure. A machine readable version of Table 2 is available at the Harvard
Dataverse (https://doi.org/10.7910/DVN/QKYR3G). For each cloud and fit (two-
component Gaussian, Plummer, single-component Gaussian) we list the average
and 1$\sigma$ uncertainty on each parameter. In each case, we obtain a set of
samples derived from the dynesty nested sampling procedure. The average value
reported is the median (50th percentile) of the samples. The upper and lower
bounds are the difference between the 50th and 84th percentiles of the
samples, and the 16th and 50th percentiles of the samples, respectively.
Cornerplots showing the marginalized posterior distributions for each cloud
and type of fit are available online at the Harvard Dataverse
(https://doi.org/10.7910/DVN/YVQC7X). The radial volume density profiles for
two clouds — Perseus and Chamaeleon — are shown in Figure 3, and similar
figures are available for the rest of the clouds in the Appendix. The data
behind these figures, including both the median and full ensemble of radial
profiles computed for each cloud, are available at the Harvard Dataverse (see
https://doi.org/10.7910/DVN/JBRWHB).

We summarize the width results in Figure 4 and Figure 5. We exclude the
single-component Gaussian from the summary due to its poor Bayes Factor
relative to the two-component Gaussian and Plummer profiles (as justified in
§4.2.1). Figure 4 shows a cornerplot of the best-fit results for the two-
component Gaussian model parameters, while Figure 5 shows a cornerplot of the
best-fit results for the Plummer model parameters. We find that for our two-
component Gaussian results, inner widths $\sigma_{1}$ span 2.5 to 4.9 pc, with
an average of $\approx$ 3.0 pc. The outer widths span 8 and 18 pc, with an
average of $\approx$ 11 pc. We find a relatively tight spread in the ratio of
the widths of the outer to inner envelopes ($\frac{\sigma_{2}}{\sigma_{1}}$),
with an average value of 3.4 and a scatter of only 0.9. The amplitude $a_{1}$
of the inner Gaussian spans $n=$ 25 - 53 $\rm cm^{-3}$, about three times as
large as the amplitude $a_{2}$ of the outer Gaussian, which spans $n=$ 5 to 15
$\rm cm^{-3}$. The widths of our Plummer fits $R_{\rm flat}$ span 1.2 - 8.1 pc
with a typical value of $R_{\rm flat}=3.6$ pc. The power-law indices $p$ of
the Plummer fit range from 1.2 and 3.4, with a median value of 1.8. The peak
density parameter $n_{0}$ ranges from $35-117\;\rm cm^{-3}$, with a median of
$\rm 49\;cm^{-3}$. When the inner 2 pc of the profile is excluded, the Plummer
fits often have a higher value of $n_{0}$ than measured, suggesting, as we
will argue in §6, that much of the densest structure in the 3D dust map is
unresolved.

####  4.2.1 Model Comparison

In addition to the best-fit parameters, we evaluate the performance of the
models in light of the data. From Bayes’ Theorem the posterior probability of
a model $\mathcal{M}$ given the data $\mathcal{D}$ is proportional to:

$p(\mathcal{M}|\mathcal{D})\propto p(\mathcal{D}|\mathcal{M})p(\mathcal{M})$
(5)

where $p(\mathcal{D}|\mathcal{M})$ is the likelihood of the data given the
model and $p(\mathcal{M})$ is the prior probability of the model. To directly
compare two models, $\mathcal{M}_{1}$ versus $\mathcal{M}_{2}$, we are
interested in the ratio of their posterior probabilities (or posterior odds),
given as:

$\frac{p(\mathcal{M}_{1}|\mathcal{D})}{p(\mathcal{M}_{2}|\mathcal{D})}=B_{12}\frac{p(\mathcal{M}_{1})}{p(\mathcal{M}_{2})}$
(6)

with the Bayes Factor, $\mathcal{B}_{12}$, equal to the ratio of the marginal
likelihood (also known as the evidence $\mathcal{Z}$) of the two models:

$B_{12}=\frac{p(\mathcal{D}|\mathcal{M}_{1})}{p(\mathcal{D}|\mathcal{M}_{2})}=\frac{\mathcal{Z}_{1}}{\mathcal{Z}_{2}}$
(7)

Assuming equal prior odds ($\frac{p(\mathcal{M}_{1})}{p(\mathcal{M}_{2})}=1$),
the posterior odds is simply the Bayes Factor. The logarithm of the evidence
$\ln(\mathcal{Z})$ is conveniently returned for each model fit by the nested
sampling package dynesty and is summarized in Table 2. The logarithm of the
Bayes Factor is then given as:

$\ln(\mathcal{B}_{12})=\ln(\mathcal{Z}_{1})-\ln(\mathcal{Z}_{2})$ (8)

Computing $\ln(\mathcal{B}_{12})$ for the two-component Gaussian and Plummer
Models, respectively, with regards to the single-component Gaussian model, we
find that the two-component Gaussian has on average $\ln(\mathcal{B}_{12})=44$
over the single-component Gaussian, whereas the Plummer model has on average
$\ln(\mathcal{B}_{12})=39$ over the single-component Gaussian. Based on the
empirical “Jeffrey’s” scale for evaluating the relative strength of two models
(Jeffreys, 1939), $\ln(\mathcal{B}_{12})>5$ represents decisive evidence in
favor of a given model. Because the Bayes Factor has a built in penalty for
model complexity, the decisive evidence in favor of the two-component Gaussian
and Plummer models is not due to the increased parameterization of the fit.
Based on these Bayes Factor calculation, we can conclude that the single-
component Gaussian fit is heavily disfavored with respect to both the two-
component and Plummer fits. When comparing $\ln(\mathcal{B}_{12})$ for the
two-component Gaussian and Plummer models, we find that the Plummer is
decisively favored for two clouds, while the two-component Gaussian is
decisively favored for four clouds, with the remaining clouds only having
moderate evidence in favor of either model. Ruling out the single-component
Gaussian model, we will discuss the theoretical implications of both the two-
component Gaussian and Plummer profiles in §6.

###  4.3 Potential Additional Sources of Fitting Uncertainty

While the formal uncertainties on all parameters are reported in Table 2,
there are potential additional systematic effects, which we discuss here.

####  4.3.1 Degeneracy of the Plummer Profile

While the Bayes Factor indicates that the Plummer function is sometimes
preferred over the two-component Gaussian function, there are clear caveats to
using a Plummer fit that should not go unsaid. In particular, it is well-known
— both from observations and simulations — that two parameters in the Plummer
function — $p$ and $R_{\mathrm{flat}}$ are highly degenerate. We see this
trend clearly in Figure 5, where $p$ and $R_{\mathrm{flat}}$ are strongly
correlated. In particular, as $R_{\mathrm{flat}}$ increases so does the value
of $p$. This covariance has also been shown in synthetic filaments on much
smaller scales (lengths $\approx$ a few parsecs, widths $\approx$ a few tenths
of a parsec) in AREPO simulations from Smith et al. (2014), where the Plummer
function was fit to the radial column density distribution. We likewise see
this in observations of the radial column density distributions of large-scale
filaments (lengths $\approx$ tens of parsecs, widths $\approx$ 1 pc) (Zucker
et al., 2018a), where $p$ and $R_{\mathrm{flat}}$ show strong correlation.
Accordingly, we urge caution when intepreting the Plummer profiles, as
discussed in §6.2.

####  4.3.2 Intra-Cut Variation in the Radial Profiles

In producing the averaged radial profile for fitting, we are not explicitly
taking into account the variation in the radial profiles along the spine of
each cloud. To summarize the variation, we have chosen to model the variable
$\sigma_{n}^{2}$, which represents the scatter in the density profile values,
capturing how much variation in the density we observe in a radial distance
bin across all cuts. In general, for our better fitting models (Plummer and
two-component Gaussian), we find that our fits prefer a small value for
$\sigma_{n}^{2}$, equivalent to a few particles per cubic centimeter in
density $n$. This trend is also evident in Figures 1 and 2. We do see a large
scatter from cut to cut at small radial distances ($r<2$ pc), but these
distances are excluded from the fit. For the remainder of the profile the
scatter is significantly reduced, and is likely driving the low inferred value
for $\sigma_{n}^{2}$. From this, we conclude that while inter-cut variation in
the radial volume density profile may have some modest affect on our inferred
value for the inner width $\sigma_{1}$ of the two-component Gaussian, it
likely has less effect on the outerwidth $\sigma_{2}$ of the two-component
Gaussian or the $R_{\mathrm{flat}}$ of our Plummer fits.

####  4.3.3 Bias in Derived Widths as a Function of Cloud Distance

In order to confirm that we are probing intrinsic variation in the cloud
widths, rather than the resolution of the underlying 3D dust map, it is
necessary to verify that our derived cloud widths are not dependent on the
distance of the cloud. If this dependence exists, we expect clouds at larger
distances to have larger widths, due to degradation in the quality of Gaia
parallax measurements underpinning the Leike et al. (2020) 3D dust map at
larger distances. To examine this potential correlation, we plot all width
parameters in Table 2 as a function of cloud distance, which we show in
Appendix B and Figure 12. We find no correlation between cloud distance and
cloud width for any width parameter, suggesting that our properties are being
driven by the intrinsic widths of clouds, rather than by distance biases.

####  4.3.4 Including the Inner 2 pc of the Radial Profiles When Fitting

In §4, we report results where radial distances between 0 and 2 pc (closest to
the spine) are excluded from the radial profiles, due to combination of the 3D
dust map’s 1 pc voxel width and the caveat from Leike et al. (2020) that
scales less than 2 pc are unreliable. To ensure that our choice to exclude the
inner 2 pc does not significantly affect our results, we refit the radial
profiles using all radial distances between 0 and 20 pc, and those results are
reported in Appendix D. Overall, we find that including the inner 2 pc has a
modest effect on our results, changing the widths on average by $\lesssim
20\%$ across all models. However, the quality of the fit degrades, validating
our choice to remove the inner 2 pc of each radial profile in our main
analysis.

####  4.3.5 Uncertainties in the Underlying 3D Dust Distribution

All results reported in this section are derived from the mean 3D dust
distribution reported in Leike et al. (2020). However, Leike et al. (2020)
also report twelve different realizations of their dust distribution, which
are used to derive the mean map. Since the different realizations reflect the
underlying uncertainty in the 3D dust reconstruction, it is necessary to
quantify the effect this uncertainty has on our results. Accordingly, in
Appendix E, we implement a more sophisticated fitting procedure, which
simultaneously takes into account all twelve realizations of the 3D dust map
when computing the model parameters. We find that overall, the uncertainty in
the 3D dust map does not dominate our results, particularly for the two-
component Gaussian fit, with best-fit parameters again typically changing
$\lesssim 20\%$ across all clouds and all model fits. See Appendix E for full
details.

####  4.3.6 Relationship Between 3D Dust Correlation Kernel and Profile Shapes

Finally, to construct their 3D dust map, Leike et al. (2020) model the dust as
a log-normal Gaussian process, while simultaneously inferring the correlation
kernel of the process, corresponding to the physical spatial correlation power
spectrum of the dust. This correlation kernel acts as a prior on the dust
extinction density101010The dust extinction density in the context of the
Leike et al. (2020) is the optical depth in the Gaia $G$ band per one parsec..
The adoption of a log-normal Gaussian process prior was primarily chosen for
statistical reasons, to ensure that the dust density is always positive, and
to allow the dust density to vary by orders of magnitude. While this kernel is
inferred from the data, it is necessary to quantify the similarity between the
shapes of the radial profiles and the shape of the inferred kernel. If the
shapes of the kernel and the profiles differ, it is likely that the data in
the vicinity of the cloud directly determines the profile shapes. However, if
the shapes of the kernel and profiles look similar, it is possible that the
kernel is being imprinted on the clouds’ structure.

In Appendix C and Figure 13 therein we compare the structure of our radial
profiles to the structure of the Leike et al. (2020) kernel inferred in the 3D
dust map’s reconstruction. We find that at radial distances between 0 and 10
pc, the cloud profiles deviate significantly from the kernel, indicating that
their shapes are being driven by data local to the cloud. This deviation
provides evidence that we are probing the intrinsic structure of the cloud at
smaller radial distances. However, at large radii ($r>10$ pc), the shapes of
the kernel and the cloud profiles become increasingly similar, particularly
for the Musca cloud. This similarity indicates that the cloud structure at
lower densities may be driven by the prior on the dust extinction density.
However, since the kernel is inferred using data inclusive of the cloud’s
structure (in addition to structure beyond each cloud of interest) the outer
envelopes may still be physically meaningful. Accordingly, while we will
discuss the physical significance of the profiles at lower densities in §6, we
caution against overinterpretation of the tails of the radial volume density
until the lower density outer envelopes can be better validated in future
work.

Table 3: Mass and Extinction Results

Cloud | ${\rm Mass}_{\rm NICEST}$ | ${\rm Mass}_{\rm Leike}$ | $\rm\frac{Mass_{\rm NICEST}}{Mass_{Leike}}$ | Maximum $A_{K\rm,NICEST}$ | Maximum $A_{K\rm,Leike}$
---|---|---|---|---|---
| $M_{\sun}$ | $M_{\sun}$ |  | mag | mag
Chamaeleon | 4.9e+03 | 4.4e+03 | 1.1 | 1.19 | 0.29
Ophiuchus | 9.4e+03 | 7.9e+03 | 1.2 | 1.99 | 0.31
Lupus | 1.0e+04 | 8.9e+03 | 1.2 | 1.55 | 0.28
Taurus | 1.6e+04 | 1.5e+04 | 1.0 | 0.90 | 0.38
Perseus | 1.6e+04 | 1.3e+04 | 1.2 | 1.40 | 0.28
Musca | 5.9e+02 | 4.9e+02 | 1.2 | 0.47 | 0.19
Pipe | 8.5e+02 | 5.2e+02 | 1.6 | 1.29 | 0.23
Cepheus | 1.2e+04 | 1.0e+04 | 1.2 | 1.15 | 0.33
Orion A* | 3.2e+04 | 9.4e+03 | 3.4 | 1.48 | 0.34
Orion B* | 2.9e+04 | 1.7e+04 | 1.7 | 1.34 | 0.33
$\lambda$ Orionis* | 1.2e+04 | 8.6e+02 | 13.7 | 1.36 | 0.22

**footnotetext: Clouds below the horizontal line – Orion A, Orion B, and
$\lambda$ Orionis – should be treated with caution, as they lie at the very
edge of the Leike et al. (2020) 3D dust grid and are subjected to additional
biases. Part of Orion A, toward the Orion Nebula Cluster, is also devoid of
dust in the Leike et al. (2020) 3D dust map despite being highly extinguished.
See §5 for full details.

Note. — Summary of mass and extinction results for local molecular clouds. (1)
Cloud Name (2) Mass derived from a traditional 2D extinction map computed
using the NICEST algorithm from Lombardi (2009). (3) Mass derived for the
cloud from the Leike et al. (2020) 3D dust map, which has been projected back
on the plane of the sky. In both cases, the mass is computed above an
$A_{K}=0.1$ mag cloud boundary. (4) Ratio of the mass from columns (2) and
(3). In (5) Maximum $A_{K}$ inside the cloud boundary based on the NICEST
algorithm. (6) Maximum $A_{K}$ inside the cloud boundary for the 2D extinction
map derived from the 3D dust data. We provide a machine readable version of
this table at the Harvard Dataverse (https://doi.org/10.7910/DVN/EIPHPR).

###  4.4 2D Extinction Maps and Masses

In this section, we compare the cloud masses derived from 3D dust to those
from traditional 2D dust approaches, with the goal of quantifying any
discrepancies. In Figure 6, for the Perseus and Chamaeleon molecular clouds,
we show a comparison between traditional 2D extinction maps (from the NICEST
algorithm; Lombardi, 2009, see §2.2) and 2D extinction maps derived by
projecting 3D dust data (from Leike et al., 2020, see §3.3) on the plane of
the sky. Both maps have been convolved to the same (lower) angular resolution
of the projected 3D dust maps, given the average distance to the cloud shown
in Table 1. In the same figures, we overlay for both NICEST (red) and 3D dust
(blue) the $A_{K}=0.1$ mag contours used to derive cloud masses in Table 3. As
will become more clear in Figure 7, despite being morphologically very similar
at the $A_{K}=0.1$ mag level, the NICEST maps are able to probe a much larger
dynamic range of extinction in comparison to the Leike et al. (2020) derived
results. In particular, the NICEST maps show much more structure at the
highest extinction levels, revealing clumps, cores, and filaments that are
invisible in the Leike et al. (2020) projected version.

This lack of structure at high extinction is not surprising, given how the
Leike et al. (2020) 3D dust map is produced. In order to probe highly
extinguished regions (see Figure 6), any extinction-based dust mapping method
must be able to see stars through dust. In contrast to NICEST, which relies
only on near-infrared photometry, the Leike et al. (2020) analysis require
Gaia optical photometry with relatively high signal-to-noise, which precludes
stars behind dense regions from being included in the analysis. Another
possible source of discrepancy between the Leike et al. (2020) projected 3D
dust and traditional 2D dust maps might stem from the calibration of the
likelihood function used to construct the Leike et al. (2020) 3D dust map.
Specifically, the likelihood used in the construction of the Leike et al.
(2020) 3D dust map is calibrated using stars in lines of sight with
approximately zero extinction, which could introduce biases in regions of
higher extinction.

This missing high extinction information in the Leike et al. (2020) 3D dust
map is more apparent in Figure 7, which compares $A_{K}$ PDFs for the
inherently 2D NICEST technique (Lombardi, 2009) and the 2D projection of the
3D maps (Leike et al., 2020), across the full cloud sample. The histograms are
computed by stacking every pixel above an $A_{K}=0.1$ mag for every cloud. The
$A_{K}$ PDF for the Leike et al. (2020) method shows more extinction than
NICEST at low values ($A_{K}\approx 0.1$ mag), and much less (essentially no)
extinction above an $A_{K}\approx 0.3$ mag. While the shapes of the two PDFs
are $\approx$ similar at the low extinction end ($A_{K}<0.15$ mag), only the
NICEST PDF shows a tail extending to high extinction ($A_{K}\approx 1.0$ mag).
This high extinction tail manifesting in the NICEST maps is consistent with
the maximum $A_{K}$ we report on a cloud-by-cloud basis in Table 3. For the
NICEST maps, the maximum $A_{K}$ toward each cloud is $\approx 0.2-0.3$ mag,
while for NICEST it is a factor of a few higher, with peak $A_{K}=0.9-1.5$ mag
inside each cloud boundary. So, these comparisons suggest that the Leike et
al. (2020) map is a faithful representation of projected dust column density
up to $A_{K}=0.2-0.3$ mag, equivalent to $A_{V}=2-3$ mag.

How does reduced sensitivity at high extinction affect masses derived from the
Leike et al. (2020) map? Table 3 lists masses for clouds (above their
$A_{K}=0.1$ mag boundaries; see §3.4) computed using both methods and the
ratio of those masses. In general, the ratio is near unity because mass at
high extinction makes up a small fraction of clouds’ total mass (see also e.g.
Lada et al., 2010). In all cases, NICEST masses are larger, but not by a
significant amount: the ratio of the NICEST mass to the 3D dust mass ranges
from 1.0 to 1.6, with an average of 1.2, for the “complete clouds,” listed
above the line in Table 3. Discrepancies are, as should be expected,
significant for the “incomplete” clouds (below the line, Orion A, Orion B, and
$\lambda$ Orionis) which lie at the edge of the Leike et al. (2020) 3D dust
map. While the mass ratio for Orion B is only modestly higher than the rest of
the sample ($\rm\frac{Mass_{NICEST}}{Mass_{Leike}}=1.7$), the NICEST mass for
Orion A is $3.4\times$ larger than its 3D dust derived mass, and for $\lambda$
Orionis it is $13.7\times$ as large. This high ratio indicates that while
Orion B is largely complete, significant amounts of mass are missing for Orion
A and $\lambda$ Orionis, not only because the Leike et al. (2020) method is
not sensitive enough at high extinction, but also likely because some fraction
of the mass lies beyond the Leike et al. (2020) grid.

Figure 6: A comparison of the $K$ band extinction ($A_{K}$) maps derived from
the NICEST algorithm (Lombardi, 2009) (red) and from the Leike et al. (2020)
3D dust map projected back on the plane of the sky (blue) for the Perseus
Molecular Cloud (left) and Chamaeleon Molecular Cloud (right). The contours
indicate the $A_{K}=0.1$ mag cloud boundaries which are used to derive the
clouds masses listed in Table 3.  Figure 7: A comparison of the PDFs between
the $K$ band extinction ($A_{K}$) derived from the 2D NICEST algorithm
(Lombardi, 2009) (red) and from the Leike et al. (2020) 3D dust map projected
back on the plane of the sky (blue). We only consider structure above the
$A_{K}=0.1$ mag threshold used to derive the cloud masses in Table 3. We
extract the extinction distribution above this $A_{K}=0.1$ mag boundary for
each cloud, and then stack the per-cloud results to obtain a representative
histogram for each method. The histograms are normalized such that their total
area is equal to one. While similar at the $A_{K}\approx 0.1$ mag level,
traditional 2D dust extinction maps show a larger dynamic range, with
increased sensitivity at higher extinction ($A_{K}\gtrsim 0.2-0.3$ mag).

##  5 Cloud-by-Cloud Exposition

In this section, we summarize the skeletonization and radial profile fitting
results of §4 on a cloud-by-cloud basis, and in light of complementary results
for each cloud from the literature. A gallery of 3D interactive figures for
the full sample is available here.

###  5.1 Chamaeleon

As shown in Figure 2 (interactive), we find that the skeleton of the
Chamaeleon Molecular Cloud is located at distances between 173 and 190 pc,
with a median distance of 183 pc. Curiously, and as apparent in Figure 2, we
find that all the Chamaeleon clouds – Chamaleon I, II, and III — are connected
via a hub-like structure, with a void in the center. The Chamaeleon complex is
also physically connected to the Musca dark cloud at a distance of 172 pc.
Musca is also shown in the 3D dust map in Figure 2 but skeletonized separately
in this work.

Based on two-component Gaussian fitting results, Chamaeleon has an inner
envelope full-width half max (FWHM$=2.35\sigma$) of 6.8 pc, and an outer
envelope FWHM of 24 pc, with a ratio of the outer to inner width equal to 3.5.

###  5.2 Ophiuchus

As shown online (interactive), we find that the skeleton of Ophiuchus is
located at distances between 115 and 150 pc. The median distance is 136 pc,
which is also the distance of L1688, consistent with previous distance
estimates based on 3D dust mapping (from Zucker et al., 2019, 2020) and from
VLBI parallax measurements of YSOs (from Ortiz-León et al., 2017). There is a
10 pc distance gradient across the B44 streamer, which extends from $\rho$
Ophiuchus out to a distance of $d\approx 145$ pc. The farthest component is
the “arc” of Ophiuchus, lying at $d\approx 150$ pc, while the nearest
component is the Ophiuchus north complex, lying between $d\approx 115-125$ pc.
The effect of feedback from OB stars (e.g. $\sigma$ Sco, $\delta$ Sco) on the
3D distribution of dust, in particular the formation of the elongated
streamers B44 and B45, will be discussed in detail in Alves et al. (2021, in
prep.). While B44 appears as part of the main spine of Ophiuchus, B45 appears
as a skeletal branch which is pruned prior to the width analysis. We find that
the B45 streamer extends from B44, rather than from the core L1688 region.
With a typical width of $\approx 1\;\rm pc$ we believe this is an artifact of
B45 being only marginally resolved in the Leike et al. (2020) 3D dust, rather
than a reflection of the true morphology of the B45 streamer.

Based on the two-component Gaussian fitting results, Ophiuchus has an inner
envelope FWHM of 6.3 pc, and an outer envelope FWHM of 24 pc, with a ratio of
the outer to inner width equal to 3.8.

###  5.3 Lupus

As shown online (interactive), the skeleton of the Lupus molecular cloud lies
between 155 and 198 pc, with a median distance of 165 pc. On a region-by-
region basis, we find Lupus 9 to be the farthest cloud in the complex, at a
distance of $\approx 197$ pc. Lupus 3-8 are all connected within the same
complex in 3D space (at the level of our thresholding algorithm, $n=\rm
35\;cm^{-3}$) at a distance of $\approx 163$ pc. Lupus 1 lies the closest in
distance, at $\approx 155$ pc, though the Lupus 1 line of sight suffers from
confusion. By confusion, we mean that there is evidence of multiple components
along the line of sight, in that there is a component immediately to the south
of Lupus 1 (appearing to connect to Lupus 1 on the plane of the sky) but which
lies at a much farther distance of 190 pc. Our distances to Lupus 1, 3, 4, 5,
and 6 are broadly consistent with those from Galli et al. (2020b), whose
distances range from $150-160$ pc based on Gaia DR2 astrometry towards young
stars associated with the complex. Curiously, we find no prominent feature at
higher gas densities ($n>\rm 35\;cm^{-3}$) associated with Lupus 2. However,
Galli et al. (2020b) place Lupus 2 at a distance of 158 pc, suggesting it is
indeed associated with the Lupus 1, 3, 4, 5, and 6 regions.

Based on the two-component Gaussian fitting results, Lupus has an inner
envelope FWHM of 6.1 pc, and an outer envelope FWHM of 19 pc, with a ratio of
the outer to inner width equal to 3.1.

###  5.4 Taurus

As shown online (interactive), the skeleton of the Taurus Molecular Cloud
extends between 131 and 168 pc. When viewed top down, it is evident that
Taurus is comprised of dust in two different layers, separated on average by
about 10-15 pc. This layering is consistent with several previous studies,
which find evidence of significant depth effects (see e.g. Galli et al., 2018,
2019; Zucker et al., 2019) over this same distance range. The furthest
component of Taurus, at a distance of 168 pc, is an elongated component of the
cloud extending out from L1538 towards $(l,b)=(177.5^{\circ},-8.5^{\circ})$.

Based on the two-component Gaussian fitting results, Taurus has an inner
envelope FWHM of 6.1 pc, and an outer envelope FWHM of 28 pc, with a ratio of
the outer to inner width equal to 4.6.

###  5.5 Perseus

As shown in Figure 1 (interactive), the skeleton of the Perseus Molecular
Cloud extends between 279 to 301 pc. It exhibits a prominent $\approx 20$ pc
distance gradient. This distance gradient is very similar to that found
independently in Zucker et al. (2018b), which also used 3D dust mapping to
obtain distances to star-forming regions across Perseus, finding them to be
located between 279 and 302 pc. This gradient is also broadly consistent with
the VLBI parallax results towards young stars in Perseus from Ortiz-León et
al. (2018), though Ortiz-León et al. (2018) find IC348 to be at a distance of
320 pc, about 20 pc farther away than we find either in this work or in Zucker
et al. (2018b). Upon further examination of the distribution of the 3D dust
and young stars in this region, we find that neither the Zucker et al. (2018b)
nor Ortiz-León et al. (2018) distance to IC348 is inaccurate. Rather, it
appears that the IC348 stellar cluster has actually migrated beyond the
distance of the dust cloud from which it formed, which results in the
discrepancy between the dust- and star-based approaches. This migration is
consistent with the estimated age of IC348 of $\approx 2-6$ Myr (Muench et
al., 2007). No such discrepancy is found for the younger NGC1333 cluster. When
considering solely the 3D dust distribution, neglecting information from YSOs,
no component of Perseus lies beyond 300 pc.

Based on the two-component Gaussian fitting results, Perseus has an inner
envelope FWHM of 9.6 pc and an outer envelope FWHM of 28 pc, with a ratio of
the outer to inner width equal to 2.9.

###  5.6 Musca

As shown online (interactive), the Musca dark cloud is located between
distance of 171 and 173 pc. We find no evidence that Musca is a sheet viewed
edge-on, but rather is a true filament in 3D space (c.f. Tritsis & Tassis,
2018). It is physically connected to the Chamaeleon Molecular Cloud complex at
a level of $n=35\;\rm cm^{-3}$ but is skeletonized separately in this work to
obtain an independent width.

Based on the two-component Gaussian fitting results, Musca has an inner
envelope FWHM of 7.3 pc, and an outer envelope FWHM of 42 pc, with a ratio of
the outer to inner width equal to 5.8. Musca is a clear outlier in the
distribution of its outer Gaussian profile, which is significantly larger and
more extended than the rest of the sample. At large radial distances, the
shape of Musca’s radial profile shows the highest correspondence to that of
the underlying kernel inferred during the reconstruction of the Leike et al.
(2020) 3D dust map (see §4.3.6 and Appendix C), indicating that this highly
extended envelope is likely prior-dominated. In addition to the outer envelope
being prior-dominated, we also believe that the inner envelope of Musca –
manifesting as an infrared dark cloud – is unresolved in our map, as 2D dust
results targeting the infrared dark component of the cloud find widths as
small as a few tenths of a parsec (Cox et al., 2016). As a result, we urge
extreme caution when interpreting these width results.

###  5.7 Pipe

As shown online (interactive), the Pipe dark cloud is located between a
distance of 147 and 163 pc, with a 15 pc distance gradient. The stem of the
Pipe nebula, near B59, is located at a distance of $150$ pc, and the bowl,
extending down to latitudes of $b=2.4^{\circ}$, lies as far as 163 pc. About
halfway down the stem, a thin tendril of gas connects the Pipe nebula with the
Ophiuchus cloud complex, meeting an extension of Ophiuchus that reaches beyond
the B44 streamer. However, the Pipe nebula is skeletonized separately from
Ophiuchus in this work, in order to obtain an independent width. The distance
to B59 in the stem of the Pipe nebula is smaller than, but still consistent
with, independent Gaia parallax measurements of YSOs in the complex (163 pc,
from Dzib et al., 2018).

Based on the two-component Gaussian fitting results, the Pipe nebula has an
inner envelope FWHM of 9.6 pc and an outer envelope FWHM of 23 pc, with a
ratio of the outer to inner width equal to 2.4. This ratio of the outer to
inner width is the smallest of all clouds in the sample.

###  5.8 Cepheus

As shown online (interactive), the skeleton of the Cepheus molecular cloud is
located between distances of 337 and 370 pc, with a median distance of 345 pc.
Our skeletonization only includes the northern section of Cepheus, as the
southern section (closer to the Galactic plane) actually lies at distances
$>900$ pc (Grenier et al., 1989; Zucker et al., 2019). The median skeleton
distance is in excellent agreement with the average Cepheus “near” distance of
352 pc, based on an independent 3D dust mapping pipeline from Zucker et al.
(2019). Our Cepheus skeletonization includes four isolated density features at
a level above our chosen threshold ($n=35\;\rm cm^{-3}$), the highest of any
cloud in the sample. These features include the sub-regions L1148/L1552/L1155
($d=338$ pc), L1172 ($d=340$ pc), L1228 ($d=368$ pc), L1241 ($d=349$ pc), and
L1247/L1251 ($d=340$ pc).

Based on the two-component Gaussian fitting results, Cepheus has an inner
envelope FWHM of 5.9 pc, and an outer envelope FWHM of 26 pc, with a ratio of
the outer to inner width equal to 4.4.

###  5.9 Corona Australis

As shown online (interactive), the skeleton of the Corona Australis cloud is
located between distances of 136 and 179 pc, with a median distance of 161 pc.
Corona Australis is one of the fluffiest local clouds in the Leike et al.
(2020) map and has no significant extended features above a density threshold
of $n=35\;\rm cm^{-3}$ so it is skeletonized at a much lower level ($n=5\;\rm
cm^{-3}$) to provide distance information. We find that the Corona Australis
cloud exhibits a prominent distance gradient, with the head of the cloud about
20-30 pc closer than the tail. The most extinguished part of Corona Australis,
which manifests as a dark cloud at optical wavelengths, lies at a distance of
$d=144$ pc, which is in good agreement with the results of Galli et al.
(2020a) (149 pc), and Zucker et al. (2019) (151 pc). The reason we determine a
higher overall median spine distance is due to the existence of Corona
Australis’ diffuse tail, which extends as far as 179 pc at a latitude of
$b=-26^{\circ}$.

###  5.10 Orion A

As shown online (interactive), we determine that the skeleton of the Orion A
cloud spans distances of 391 to 445 pc, with a median distance of 405 pc. The
3D dust distribution towards Orion A is peculiar, due to the fact that there
is a large piece of the cloud missing towards the Orion Nebula Cluster — the
densest part of the entire complex. This lack of dust is not physical, and
represents the inability of the Leike et al. (2020) to fully probe high gas
densities. As a result, the head (toward the Orion Nebula Cluster) and tail of
Orion A (toward L1647) appear not to connect as they should in 3D space.
Analysis of Orion’s YSO distribution (Großschedl et al., 2018) and other high
resolution 3D dust maps of Orion A (Rezaei Kh. et al., 2020) clearly show that
the head and tail are physically connected. The missing ONC region could
explain in part the significant discrepancy between the cloud mass derived
from 2D dust extinction versus from the projected 3D dust distribution (see
§4.4).

In general, we find good agreement with previous studies towards regions
probed by the dust map, as well as a similar range of distances. For example,
Großschedl et al. (2018) find that the bulk of the Orion A cloud lies between
distances of $d\approx 380-460$ pc, which is not inconsistent with our spread
of $391-445$ pc. However, there is no visible connection between the farthest
point in Orion A (L1647-S at a distance of 445 pc) and the main Orion A cloud
consisting of L1641, while this connection does appear in other approaches
(Großschedl et al., 2018; Rezaei Kh. et al., 2020). This lack of connection at
the tail could be physical, or it could be indicative of another missing mass
situation, like toward the Orion Nebula Cluster. In addition to mass missing
within the grid, it is also possible that there is mass missing beyond the
grid, since Orion A lies at the boundary of the Leike et al. (2020) dust mass.
Thus, our skeletonization and width analysis results for Orion A should be
treated with more caution than for other “complete” clouds. The most likely
effect of Orion A lying at the edge of the grid is that we will slightly
underestimate the distance to the complex.

Based on the two-component Gaussian fitting results, we find that Orion A has
an inner envelope FWHM of 6.8 pc and an outer envelope FWHM of 20 pc, with a
ratio of the outer to inner width equal to 2.9.

###  5.11 Orion B

As shown online (interactive), we determine that the skeleton of the Orion B
cloud spans distances of 397 to 406 pc, with a median distance of 404 pc. Our
median distance is consistent with Kounkel et al. (2018), which leverages Gaia
DR2 parallax measurements toward young stars in the complex to obtain a
distance of $407\pm 4$ pc. Unlike Orion A, we find no evidence of a strong
distance gradient, and also find that Orion B consists of a single density
feature with a length of 69 pc. However, since Orion B again lies at the very
edge of the Leike et al. (2020) dust map, it is possible that parts of Orion B
could still be missing, which would bias our distance estimates lower compared
to a cloud that lies well within the map boundaries. Thus, like our
skeletonization and width analysis results for Orion A, we caution that our
results for Orion B should be treated with more caution than for other
“complete” clouds.

Based on the two-component Gaussian fitting results, Orion B has an inner
envelope FWHM equal of 11.4 pc and an outer envelope FWHM of 34 pc, with a
ratio of the outer to inner width equal to 3.0.

###  5.12 $\lambda$ Orionis

Finally, as shown online (interactive), we determine that the skeleton of the
$\lambda$ Orionis cloud spans distances of 375 to 397 pc, with a median
distance of 386 pc. A clear ring-like morphology is observed in 2D dust
extinction maps, which Dolan & Mathieu (2002) argue formed from ejected
material driven by a supernovae blast at the center of the ring about 1 Myr
ago.

We observe the same ring-like morphology in 3D. However, the skeletonization,
and resulting widths, should be considered with extreme caution, as the
skeleton lies only a few parsecs from the boundary of the Leike et al. (2020)
3D dust map, so the skeletonization is likely to underestimate the cloud’s
distance. There is clearly mass missing in the $\lambda$ Orionis region, given
that the mass for the cloud based on 2D dust maps (from NICEST) is a factor of
fourteen higher than what we calculate using just the 3D dust volume local to
the cloud (see §4.4 and Table 3). Consistent with this missing mass scenario,
average distances to the $\lambda$ Orionis region tend to be higher in other
studies (e.g. 402 pc from Zucker et al., 2020). The average distance to
$\lambda$ Orionis from Zucker et al. (2020) is also consistent with
independent distances of young stars, for which Kounkel et al. (2018) obtain
$d=404$ pc.

Based on the two-component Gaussian fitting results, $\lambda$ Orionis has an
inner envelope FWHM of 8.7 pc and an outer envelope FWHM equal of 30 pc, with
a ratio of the outer to inner width equal to 3.4. Despite only measuring the
cloud width using voxels in front of $\lambda$ Orionis (due to the Leike et
al. (2020) grid ending directly beyond the spine), these results are
consistent with results for the “complete” clouds.

##  6 Discussion

In Section 4, we have shown that the radial volume density profiles measured
in 3D for local molecular clouds are well-described by two functions: a two-
component Gaussian function and a Plummer function. A single-component
Gaussian is universally a poor fit to the data. We have also shown that the
Leike et al. (2020) technique is capturing a majority of the cloud mass for
all “complete” clouds in the sample, deviating from 2D integrated dust
extinction measurements (from the NICEST methodology; Lombardi, 2009; Lombardi
et al., 2011) only at higher extinction ($A_{K}\gtrsim 0.3$ mag). In this
section, we discuss the potential implications of our results, and provide an
initial comparison of the 2D and 3D shapes of the radial profiles for a subset
of the sample, as a basis for future work.

###  6.1 Two-Component Gaussian — Probing a Phase Transition?

In terms of the two-component Gaussian model, one possible explanation for the
shape of the radial density profiles is that we are probing gas in different
phases, where the transition between the outer and inner Gaussian profiles
represents a shift in temperature, density, and/or chemical composition of the
gas. In the density range consistent with our radial profiles, possible
transitions include a chemical transition between atomic and molecular
hydrogen (the HI-to-$\rm H_{2}$ transition), or a thermal transition between
the unstable neutral medium (UNM) and cold neutral medium (CNM), driven by the
temperature dependence of the atomic cooling curve. While the data are
suggestive, we do not yet have enough detailed evidence to prove the link
between the shape of the radial volume density profiles and a phase
transition. Thus, in what follows, we only seek to present a plausibility
argument.

####  6.1.1 Chemical Transition Between Atomic and Molecular Gas

The first explanation for the bimodality in the radial profiles is that the
broader, lower amplitude “outer” Gaussian probes the shape of the atomic HI
gas envelope, while the narrower higher amplitude “inner” Gaussian probes the
molecular phase. If so, the noticeable break in volume density coincident with
the start of the inner Gaussian would represent a chemical transition of the
gas, delineating the radial distance in 3D space at which point molecular
hydrogen becomes self-shielded. Models suggest that the transition to self-
shielding — at which point molecular hydrogen can persist — can be either
sharp or gradual, depending on the strength of the UV radiation field relative
to the gas density. However, once a critical column density of $\rm H_{2}$ has
been produced, the fraction of hydrogen in molecular form should rise rapidly
(Bialy & Sternberg, 2016; Sternberg et al., 2014; Krumholz et al., 2008). If
this transition is occurring in local molecular clouds — as predicted by
analytic theory — it might explain why we find the two-component Gaussian
model to be so stongly favored over a single-component Gaussian model.

We have known for a few decades that HI envelopes around molecular clouds are
common (see e.g. Moriarty-Schieven et al., 1997; Williams & Maddalena, 1996;
Williams et al., 1995). These HI envelopes are significantly extended with
respect to the molecular gas and may be either the remnants of atomic clouds
which condensed to form GMCs, or photodissociated gas on the outskirts of the
molecular cloud. More likely, the atomic envelope arises from a combination of
both phenomena. In addition to being observed in solar neighborhood clouds
(e.g. in Perseus and Taurus; see Stanimirovic et al., 2014; Pascucci et al.,
2015; Lee et al., 2015) and toward the inner galaxy (Beuther et al., 2020;
Bialy et al., 2017), atomic envelopes have also been observed at much lower
resolution in the nearby Large Magellanic Cloud (Fukui et al., 2009). Recent
numerical simulations characterizing the properties of highly resolved
synthetic molecular clouds (e.g. Duarte-Cabral & Dobbs, 2016; Seifried et al.,
2017) also note the presence of large extended HI envelopes.

We find that the typical amplitudes of the outer Gaussian profiles are
consistent with the typical densities of the atomic envelopes seen in
observations, mainly via HI emission on the plane of the sky. In this study,
the peak densities of the outer envelopes — characterized by the amplitude
parameter $a_{2}$ — span $n=5-15\;\rm cm^{-3}$, with an average amplitude of
$n=9\;\rm cm^{-3}$, falling off to only a few particles per cubic centimeter
at large radii. This range is consistent with HI densities inferred for atomic
envelopes around molecular clouds in both observations and simulations.
Looking at a sample of 123 HI envelopes in the Large Magellanic Cloud, Fukui
et al. (2009) derive a typical envelope density of $\rm 10\;cm^{-3}$,
calculated by dividing each cloud’s HI column density by its approximate size
on the plane of the sky.111111We note, however, that the HI-to-H2 transition
strongly depends on the dust-to-gas ratio, and the Large Magellanic Cloud has
lower metallicity than Milky Way (Van der Swaelmen et al., 2013), so any
direct comparison between the two systems should be treated with caution.
Using GALFA HI and IRAS data, Bialy et al. (2015) also calculate the atomic
hydrogen volume density towards the Perseus molecular cloud, finding that the
atomic densities vary between $n=2-11\;\rm cm^{-3}$. This range of densities
is similar to the range of densities captured by the outer Gaussian in the
radial volume density profile we measure for Perseus, spanning $n=12\;\rm
cm^{-3}$ toward the spine, down to $n\approx 2\;\rm cm^{-3}$ toward the
outskirts of the cloud.

In terms of the inner Gaussian, we find typical peak densities — parameterized
by $a_{1}$ — to lie between $n\approx 25-52\;\rm cm^{-3}$. While the amplitude
for the averaged profile is of this order, a significant fraction of cuts show
peak densities with $n=100\;\rm cm^{-3}$, with some extending as high as
$n=200\;\rm cm^{-3}$. While these average volume densities are not high enough
to be observationally probed by CO-bright molecular gas manifesting at
$n\gtrsim 200\;\rm cm^{-3}$ (see e.g. Busch et al., 2019), they are consistent
with densities of so-called “CO-dark” gas, a phase associated with warmer,
more diffuse molecular hydrogen which is not detectable in CO or HI. While
invisible in both CO and HI, $\gamma$ ray observations — capable of tracing a
cloud’s gas mass independent of its chemical state — show that diffuse
molecular gas surrounds all nearby CO-detected molecular clouds targeted in
this study. This CO-dark phase effectively acts as an observational bridge
between clouds’ extended atomic envelopes and their embedded CO-bright cores
(Grenier et al., 2005), and would mark the true chemical transition between
atomic and molecular gas. Tang et al. (2016) use [CII] emission to trace CO-
dark diffuse molecular gas in the Galactic plane, finding that this gas
manifests at a range of densities (from $n\approx 60-1000\;\rm cm^{-3}$), with
an average density of $n\approx 200\;\rm cm^{-3}$.121212While [CII] emission
is commonly used as a tracer of CO dark gas, we note that a recent study by
Hall et al. (2020) towards the local Perseus Molecular Cloud concludes that a
[CII] envelope can also arise from the CNM phase, without a need to invoke CO-
dark molecular gas. The average density from Tang et al. (2016) is a factor of
a few larger than our average peak density $a_{1}$, but the lower estimates
from Tang et al. (2016) remain consistent with our results. Also consistent
with these observations, analysis of the CO-dark diffuse molecular gas for
synthetic clouds from the SILCC-Zoom project (Seifried et al., 2019) find
number densities between $10-1000\;\rm cm^{-3}$, with an average density of
CO-dark gas on the order of $100\;\rm cm^{-3}$ — only a factor of two higher
than what we estimate here for the inner Gaussian. Given that the volume
densities derived from the 3D dust data — particularly for the inner Gaussian
component — should be lower limits (see also §2.1.2), our results are
generally in agreement with literature estimates for CO-dark molecular gas
from both observations and simulations.

Finally, there is a straightforward test to disprove this HI-to-$\rm H_{2}$
transition hypothesis, involving a comparison between the radial extent of CO
emission projected on the plane of the sky and radial extent of the inner
Gaussian component derived here. The detection of CO emission significantly
beyond the boundaries of the inner Gaussian profile would be in direct
contradiction to the proposed origin, since no molecular gas should lie beyond
the inner Gaussian. To test this scenario, we build CO radial profiles for a
subset of the clouds in the sample, as described in detailed in Appendix F. We
find that the CO profiles shows widths equal to or marginally smaller than
that of the inner Gaussian profiles, consistent with the scenario that the
inner Gaussian profiles could be tracing CO dark gas.

####  6.1.2 Thermal Transition Between Unstable and Cold Neutral Medium

While a diffuse molecular gas component offers one explanation for the shape
of the inner Gaussian profile, no $\rm H_{2}$ component necessarily needs to
be invoked to explain the clouds’ density structure. It is well established
that the atomic neutral medium consists of two thermally stable states at
significantly different densities and temperatures — the warm neutral medium
(WNM) at $n\approx 0.2-0.9\;\rm cm^{-3}$ and the cold neutral medium (CNM) at
$n\approx 7-70\;\rm cm^{-3}$ (Field et al., 1969; Wolfire et al., 2003; Murray
et al., 2018b). This two-phase structure is a natural consequence of the
density and temperature dependence of the cooling curve at play in the diffuse
interstellar medium. While the peak densities we find for the inner Gaussian
($a_{1}=25-52\;\rm cm^{-3}$) are consistent with lower-density estimates of
CO-dark diffuse molecular gas, they are also in strong agreement with the
densities of the CNM. However, the peak densities of the outer envelopes are
much larger than expected for the WNM. In traditional two-phase models, HI at
intermediate temperatures should be expected to quickly evolve into one of the
stable phases. However, ample evidence has recently suggested that a
significant fraction of Galactic HI lies in the thermally unstable phase,
driven by dynamical processes in the interstellar medium (e.g. turbulence;
c.f. Heiles & Troland, 2003; Roy et al., 2013; Murray et al., 2018b). Thus, if
the inner envelope is probing the CNM, then the outer envelopes (with peak
densities $a_{2}=5-15\;\rm cm^{-3}$) may probe the thermally unstable HI lying
at densities between the WNM and CNM with $n\approx 1-10\;\rm cm^{-3}$. The
transition observed in the two-component Gaussian profiles would then be
driven not by a change in the chemical state of the gas, but by a change in
the thermal phase of the atomic neutral medium.

Consistent with this scenario, Murray et al. (2018b) analyzed HI spectra
toward 57 lines of sight taken from the 21-SPONGE survey to infer the fraction
of HI mass in the thermally unstable regime, finding that 20% is in the UNM
phase, as opposed to 28% in the CNM phase (with the remaining mass in the WNM
phase). Murray et al. (2018b) further find no significant difference in the
UNM fraction as a function of latitude, with both high ($|b|>10^{\circ}$) and
low latitude lines of sight showing similar mass distributions. Given its
ubiquity, both from Murray et al. (2018a) and complementary studies (Roy et
al., 2013; Heiles & Troland, 2003), the UNM-to-CNM transition represents a
viable alternative to the HI-to-$\rm H_{2}$ transition for explaining the two-
component Gaussian shape of the radial profiles. We note, however, that while
the Murray et al. (2018b) lines of sight span a wide range in latitudes, none
of them intersect the clouds targeted in this work.

Nevertheless, to broadly verify the consistency of our results with those of
Murray et al. (2018b), we can compare for each cloud the relative amount of
mass contained in the outer and inner envelope. If our results are consistent
with Murray et al. (2018b), we would expect the ratio of the outer and inner
envelope masses to be $\frac{M_{2}}{M_{1}}=\frac{M_{\rm UNM}}{M_{\rm
CNM}}=\frac{20\%_{\rm UNM}}{28\%_{\rm CNM}}=0.71$. We compute this ratio in
two ways. First, we compute the mass of the inner envelope, $M_{1}$, and outer
envelope $M_{2}$ by taking the definite integral of the Gaussian density
profiles (shown via the black-white dotted and cyan-black dotted curves in
e.g. Figure 1) over the volume, assuming a cylindrical geometry:

$M_{1}=2\pi l\,\mu
m_{p}\int^{+\infty}_{-\infty}a_{1}e^{\frac{-r^{2}}{2\sigma_{1}^{2}}}\,r\,dr$

$M_{2}=2\pi l\,\mu
m_{p}\int^{+\infty}_{-\infty}a_{2}e^{\frac{-r^{2}}{2\sigma_{2}^{2}}}\,r\,dr$

where the parameters $a_{1}$, $\sigma_{1}$, $a_{2}$, $\sigma_{2}$ are the same
as those reported in Table 2, $\mu=1.37$ is the mean molecular weight,
corrected for the helium abundance, $m_{p}$ is the mass of a proton, and $l$
is the length of the cloud reported in Table 1. Computing $M_{1}$ and $M_{2}$
numerically and taking the ratio across the sample we obtain a median ratio
for $\frac{M_{2}}{M_{1}}=2.9$, with a spread of 1.3, which is higher than the
predicted Murray et al. (2018b) results of 0.72. Given that we are likely
underestimating the densities in potentially CNM-associated regions, one
possible explanation is that while the typical extents of the inner and outer
components could be correct, their amplitudes, particularly for the inner
Gaussian, may not be. Thus, a second, more extreme way to test this is to
assume that the inner envelope is entirely composed of CNM, and UNM only
exists beyond the edge of the inner envelope. To compute this, we create
cylindrical representations of the inner and outer envelopes, where the width
of each 3D cylinder is equivalent to the FWHM derived from the $\sigma_{1}$
and $\sigma_{2}$ values reported in Table 2 (FWHM $=2.355\times\sigma$). These
cylinders are toggle-able on and off in the gallery of interactive cloud
figures. The cylinders define a mask in 3D volume density space over which the
masses are computed:

$M_{\rm cylinder_{inner}}=\Sigma_{i}\;\mu m_{p}\times n_{i,\rm inner}\;dV$

$M_{\rm cylinder_{outer}}=\Sigma_{i}\;\mu m_{p}\times n_{i,\rm outer}\;dV$

Here, $dV$ is the volume element of an individual voxel $i$ in the 3D volume
density map, equal to 1 pc3. The $\Sigma_{i}\;n_{i,\rm inner}$ and
$\Sigma_{i}\;n_{i,\rm outer}$ terms represent the sum over the set of voxels
(in gas density space $n$) which lie within the 3D inner and outer cylindrical
regions, respectively. Under the assumption that no UNM lies within the bounds
of the inner envelope, the resulting mass ratio would be:

$\frac{M_{2}}{M_{1}}=\frac{M_{\rm cylinder_{outer}}-M_{\rm
cylinder_{inner}}}{M_{\rm cylinder_{inner}}}$

As expected, we calculate a lower value of $\frac{M_{2}}{M_{1}}$, with a
median of 2.1 and a spread of 1.0. In reality the true mass ratio is likely
somewhere between the two examples given here, and additional follow-up work
will be needed to rule out the possibility of a UNM-to-CNM transition.

####  6.1.3 Constant Ratio Between Outer and Inner Cloud Envelopes

Regardless of whether the transition between the outer and inner envelopes is
due to a a) chemical transition between HI and $\rm H_{2}$ gas b) a thermal
transition between the UNM and the CNM; or neither, we find that the ratio
between the widths of the outer and inner envelopes across the sample is
roughly 3.4, with a relatively small scatter of $\approx 1$. This ratio could
be physically significant, with the potential to constrain the formation
timescale of $\rm H_{2}$, or the rate at which gas flows between different
thermal phases of the diffuse interstellar medium (e.g. between the UNM and
CNM). More targeted analytic modeling (e.g. in the style of Sternberg et al.,
2014; Bialy & Sternberg, 2016; Krumholz et al., 2009), further
contextualization with CO and HI spectral-line maps, and comparisons to
simulations that track HI, $\rm H_{2}$ and CO formation and destruction in a
time-dependent chemical network (e.g. the Cloud Factory Simulations; Smith et
al., 2020) should shed more light on this behavior.

Figure 8: Comparison between 2D radial column density profiles and 3D radial
volume density profiles for the Perseus (left), Chamaeleon (middle) and Taurus
(right) molecular clouds. The thick blue lines show 3D volume density derived
from the Leike et al. (2020) 3D dust map. The thin blue lines show 2D column
density profiles derived by projecting the Leike et al. (2020) 3D dust map on
to the plane of the sky (§3.3). The red lines show 2D column density from the
NICEST map (Lombardi, 2009). The green lines show profiles derived from 2D
dust emission maps (Schlegel et al., 1998). All maps have been scaled to a
peak value of 1.0, to facilitate a comparison of the shapes of the volume and
column density profiles. While the 2D and 3D results largely agree at radial
distances $<5$ pc, significant deviations are observed at larger radial
distances, particularly for the Chamaeleon and Taurus clouds.

###  6.2 An Isothermal Self-Gravitating Cylinder?

For the past ten years, Plummer models (originally describing the stellar
surface density of globular clusters; Plummer, 1911) have often been used to
characterize the radial profiles of filaments extracted from $\rm H_{2}$
column density maps as part of the Herschel Gould Belt Survey (André et al.,
2010). These studies (Arzoumanian et al., 2011, 2019; Palmeirim et al., 2013;
Cox et al., 2016) adopt an idealized Plummer-like analytic model for the
(column) density profile of a cylindrical filament, featuring a denser, flat
inner portion and an outer portion that falls off like a power-law at larger
radii:

$n(r)=\frac{n_{0}}{\left[1+\left({r/R_{\rm flat}}\right)^{2}\right]^{p/2}}\
\longrightarrow\Sigma(r)=A_{p}\frac{N_{\rm 0}}{\left[1+\left({r/R_{\rm
flat}}\right)^{2}\right]^{\frac{p-1}{2}}},\ \ $ (9)

where $n(r)$ describes the volume density profile adopted throughout this
work, and $\Sigma(r)$ is the corresponding column density profile employed in
the literature, with $N_{0}$ the peak column density, $R_{\mathrm{flat}}$ the
flattening radius, $p$ the index of the density profile, and $A_{p}$ is a
constant factor that accounts for the filament’s inclination angle to the
plane of the sky. In column density studies, the inclination of the filament
to the plane of the sky is assumed to be zero, which is not the case in our
volume density analysis. Following the formalism outlined in Ostriker (1964)
and Equation 9, the density structure of an isothermal, self-gravitating
gaseous cylinder in hydrostatic equilibrium corresponds to a special case with
a power-law index $p=4$ or $n\propto r^{-4}$ at large radii ($r>>R_{\rm
flat}$). Analyzing the Herschel column density distribution of filaments,
various works (Arzoumanian et al., 2011, 2019; Palmeirim et al., 2013; Cox et
al., 2016) have found that these filaments are inconsistent with the $p=4$
scenario, instead preferring much shallower power-law indices ($p=1.5-2.5$).

The Herschel filaments analyzed in previous work lie within the local
molecular clouds we analyze here, persisting at much higher column densities
($N_{\rm H_{2}}>10^{22}\;\rm cm^{-2}$), and on much smaller scales (widths
$\approx$ a few tenths of a parsec, lengths $\approx$ a few parsecs). As such,
they are unresolved in the Leike et al. (2020) 3D dust maps we analyze.
Nevertheless, the same formalism can be applied to the much larger molecular
clouds these dense filaments inhabit, since the analytic models are scale-
free, and we have shown here that local molecular clouds are well-described as
cylindrical filaments, even at lower densities ($n<100\;\rm cm^{-3}$), where
their aspect ratios are often in excess of $3:1$ (see lengths and widths in
Tables 1 and 2). Thus, while we target a very different regime – about two
orders of magnitude lower in density, and an order of magnitude larger spatial
scale – the comparison to the classical $p=4$ Ostriker (1964) model for an
isothermal, self-gravitating filament in hydrostatic equilibrium, is still
relevant.

We find that the local, large-scale, molecular clouds studied here have
shallower radial density profiles, with typical power-law indices $p$ between
$1.5-2.5$, than would an isothermal self-gravitating cylinder in hydrostatic
equilibrium. This shallower power-law, while inconsistent with isothermal
self-gravitating filaments in the absence of a magnetic field, could be
consistent with some models of magnetized filaments (Fiege & Pudritz, 2000).
Previous studies of the structure of the magnetic field in solar neighborhood
clouds find a strong coupling between the gas column density and the magnetic
field orientation, which is consistent with the idea that magnetic fields play
an important role in structuring the interstellar medium in and around local
molecular clouds (see e.g. Soler, 2019).

However, given that the parameters of the Plummer profile are highly
degenerate (see discussion in §4.2), and driven by behavior at the tail of the
distribution (at very low densities), we are hesitant to assign any physical
significance to the shape of the Plummer profile, other than broadly noting
the inconsistency with the Ostriker (1964) result. Similar to the two-
component Gaussian results, future comparison with numerical simulations
should allow for more robust interpretation of the physical significance, if
any, of the shape of the Plummer profiles.

###  6.3 Comparison between 3D Volume Density and 2D Column Density Radial
Profiles

The ultimate goal of this work is to lay the foundation for a robust
characterization and interpretation of the structure of molecular clouds. To
do so will require the contextualization of these newly presented 3D results
with extant 2D approaches. Traditionally, there have been three ways of
measuring the column density structure of molecular clouds: near-infrared
extinction mapping (e.g. the NICEST method, §2.2; Lombardi, 2009), thermal
dust emission mapping in the far-infrared (e.g. the SFD or Planck dust maps;
Schlegel et al., 1998; Planck Collaboration et al., 2016), and spectral-line
mapping of CO isotopologues. A systematic analysis of and comparison to the 2D
radial profiles of molecular clouds is beyond the scope of this work. So, as a
precursor to a more expansive study, we offer in Figure 8 an initial
comparison of our 3D radial volume density profile results with traditional 2D
column density results for three clouds in the sample (Perseus, Chamaeleon,
and Taurus). We compare the 3D radial profiles with three measurements of the
2D radial profile: one obtained from near-infrared extinction using the NICEST
maps (see §2.2), one obtained from far-infrared dust emission maps (Schlegel
et al., 1998), and one obtained by projecting the 3D dust map from Leike et
al. (2020) on to the plane of the sky, using the same column density maps from
§3.3. Each curve has been smoothed with a kernel roughly equal to 2 pc (the
effective spatial resolution of the Leike et al. (2020) dust map), to
alleviate noise in the profiles at large radial distances and facilitate
comparison between the different tracers.131313In all three cases, the radial
profiles are obtained identically to our 3D results presented in §3.2.1,
except the threshold for mask (and thus skeleton) creation was defined in 2D
column density space as opposed to 3D volume density space. All 2D maps are
convolved to the same resolution. The difference in “zero point” emission and
extinction levels in traditional 2D approaches precludes the application of an
identical column density threshold for all three column density maps. To avoid
biases in our analysis, we define thresholds using the total area of the cloud
masks, rather than absolute column density. That is, we adopt the same
$A_{K}=0.1$ mag threshold (converted to $N_{\rm H}$) used in the §3.4 NICEST
results to define the optimal mask, and set the thresholds for the 2D dust
emission (Schlegel et al., 1998) and projected Leike et al. (2020) 2D results
so they occupy approximately the same area. By doing so, all masks are defined
using a narrow range in column density, spanning $\rm N=1-3\times
10^{21}\;cm^{-2}$. A tutorial for how to build 2D radial column density
profiles with the RadFil package is available on Zenodo (Zucker & Chen, 2018).
We largely defer comparison to the 2D structure of clouds obtained from
velocity-integrated spectral-line maps to future work, given the need for
careful accounting of opacity effects and the presence of CO dark gas. In
order to facilitate a comparison of the shapes of the profiles, we scale the
radial volume density profile and the three radial column density profiles to
a peak value of one.

In contrast to the 3D radial profiles, the 2D radial profiles show a much
larger variation in structure both within and between clouds. In general, we
find good agreement between 2D and 3D results at small radial distances ($r<5$
pc), particularly for the Perseus and Taurus clouds. However, at radial
distances $>5$ pc, we observe a flattening of the 2D radial profiles that does
not manifest in 3D. These larger radial distances are beyond the $A_{K}=0.1$
mag boundary used to validate the 3D-dust-based cloud masses in §4.4. We
observe flattening at large radii even in the projected dust map from Leike et
al. (2020) (thin blue trace), which makes use of the same data as the 3D
radial profile (thick blue trace), but simply integrates along the line of
sight. The flattening seen (only) in 2D suggests two underlying causes. First,
diffuse dust unassociated with (but potentially in the same general vicinity
of) the cloud could be contaminating the inferred radial profiles at large
radial distances. Second, projection effects intrinsic to the cloud itself
(e.g. due to the cloud’s orientation with respect to the observer) could be
driving this variation. In either or both scenarios, discrepancy between 2D
and 3D cloud structure will have profound implications for interpreting the
shapes of the column density PDFs of molecular clouds (c.f. Alves et al.,
2017; Chen et al., 2019). These deviations also underline the difficulty of
validating the atomic outer envelopes of clouds and the tails of the Plummer
profiles (§6.2) using traditional 2D dust maps. Disentangling these effects
will require careful calibration of 2D dust extinction and emission maps;
careful supplementation of the 2D dust information with opacity-corrected
spectral line maps of CO, HI, and [CII], and contextualization of both types
of data in light of synthetic observations of local cloud analogs.

##  7 Conclusion

The rise of Gaia has transformed our understanding not just of the 3D
distribution of stars, but also of the 3D distribution of dust that reddens
their colors. The recent release of the Leike et al. (2020) 3D dust map, which
provides 1 pc spatial resolution out to a distance of several hundred parsecs,
offers an unprecedented opportunity to measure and describe the shape and
internal structure of local molecular clouds on scales previously obtainable
only in simulations. Converting from differential dust extinction in the Gaia
G-band to volume density using a wavelength-dependent extinction curve, we
characterize the shapes and thicknesses of famous local clouds in 3D gas
density space for the first time. Local clouds appear filamentary, even at gas
densities $n<100\;\rm cm^{-3}$. Thus, we start by extending the FilFinder
algorithm – originally applied to 2D maps of dust emission (Koch & Rosolowsky,
2015) — into 3D, to determine the skeletons of nearby clouds, equivalent to a
one-voxel-wide representation of their 3D structure. We then extend the RadFil
algorithm (Zucker & Chen, 2018) to measure the radial volume density profiles
with respect to the cloud spines. Our results are as follows:

  * •

By projecting each cloud’s skeleton – determined in 3D $(x,y,z)$ cartesian
space — back to sky coordinates $(l,b,d)$, we constrain the distances to
famous local molecular clouds. The 1 pc uncertainty on the distances allows
for complete characterization of distance gradients and the presence of
multiple components.

  * •

We model each cloud’s radial volume density profile with three functions: a
two-component Gaussian, a single-component Gaussian, and a Plummer-like
function. Comparing the relative quality of models, we find that a single-
component Gaussian function is universally disfavored, with significantly
worse odds compared with the two-component Gaussian and Plummer models.

  * •

For the two-component Gaussian model, we find widths $\sigma_{2}$ of the
“outer” Gaussian between $8-18$ pc, with amplitudes between $n=5-15$ $\rm
cm^{-3}$. For the “inner” Gaussian, we find widths $\sigma_{1}$ of $2.5-4.9$
pc, with amplitudes between $n=25-52$ $\rm cm^{-3}$. The typical ratio of the
outer to inner widths is $3.4:1$, with relatively small scatter ($<1$). We
hypothesize that the boundary might represent a chemical transition between HI
and $\rm H_{2}$ gas in local clouds, with the inner regions probing warm,
diffuse molecular gas, and the outer regions tracing the extended atomic
envelopes. Another plausible scenario includes a thermal transition between
the unstable neutral medium (UNM) and the cold neutral medium (CNM). Future
comparison with numerical simulations (e.g. Smith et al., 2020) will be needed
to verify these hypotheses.

  * •

We produce synthetic 2D dust extinction maps ($A_{K}$) by projecting the Leike
et al. (2020) 3D dust distribution on the plane of the sky. We compare the 2D
extinction maps derived from 3D dust with traditional integrated 2D dust
extinction maps (from the NICEST algorithm; Lombardi, 2009). Computing cloud
masses, we determine that the Leike et al. (2020) map is recovering a majority
of each cloud’s mass. However, the 3D dust is insensitive to structure above
$A_{K}=0.3$ mag ($A_{V}=3$ mag)

  * •

All data underpinning this work are publicly available online at the Harvard
Dataverse (https://dataverse.harvard.edu/dataverse/3D_cloud_structure/). An
online gallery of interactive figures illustrating our cloud results is also
available online here.

Ultimately, these results are only the beginning of what is possible through
the combination of Gaia astrometry, large photometric surveys, innovative
computational and statistical techniques, and exploratory 3D data
visualization. The quality of the Leike et al. (2020) 3D dust map, and of
similar maps extending out many more kiloparsecs (Green et al., 2019; Zucker
et al., 2020), relies on accurate measurements of the distance and integrated
extinction to individual stars. Not only will future Gaia data releases
produce increasingly better estimates of stellar distance, but upcoming
photometric surveys like LSST will allow for improved constraints on the
integrated extinction. Future 3D dust maps, built on much larger catalogs of
stellar properties (e.g. the brutus algorithm; Speagle et al. 2021, in prep.)
should allow for improved models of the 3D dust distribution, probing not only
clouds at larger distances, but also with improved sensitivity to higher
densities within molecular clouds.

We would like to thank our anonymous referee, whose thoughtful reading of our
manuscript improved the quality of this work, particularly with regards to the
physical interpretation of the radial profile results. We would like to thank
Charles Lada, Thomas Dame, Robert Benjamin, and Gina Panopoulou for helpful
feedback and discussions which improved the quality of this manuscript. The
visualization, exploration, and interpretation of data presented in this work
was made possible using the glue visualization software, supported under NSF
grant numbers OAC-1739657 and CDS&E:AAG-1908419.
D.P.F. and C.Z. acknowledge support by NSF grant AST-1614941, “Exploring the
Galaxy: 3-Dimensional Structure and Stellar Streams.”
We thank the developers of the X-Toolkit (Haehn et al., 2014), which enabled
the fast volume rendering of our results in the browser, in Figures 1 and 2.

## References

  * Alves et al. (2017) Alves, J., Lombardi, M., & Lada, C. J. 2017, A&A, 606, L2, doi: 10.1051/0004-6361/201731436
  * Alves et al. (2020) Alves, J., Zucker, C., Goodman, A. A., et al. 2020, Nature, 578, 237, doi: 10.1038/s41586-019-1874-z
  * Anders et al. (2019) Anders, F., Khalatyan, A., Chiappini, C., et al. 2019, A&A, 628, A94, doi: 10.1051/0004-6361/201935765
  * André et al. (2010) André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102, doi: 10.1051/0004-6361/201014666
  * Arras et al. (2019) Arras, P., Frank, P., Leike, R., Westermann, R., & Enßlin, T. A. 2019, A&A, 627, A134, doi: 10.1051/0004-6361/201935555
  * Arzoumanian et al. (2011) Arzoumanian, D., André, P., Didelon, P., et al. 2011, A&A, 529, L6, doi: 10.1051/0004-6361/201116596
  * Arzoumanian et al. (2019) Arzoumanian, D., André, P., Könyves, V., et al. 2019, A&A, 621, A42, doi: 10.1051/0004-6361/201832725
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipócz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Beaumont et al. (2013) Beaumont, C. N., Offner, S. S. R., Shetty, R., Glover, S. C. O., & Goodman, A. A. 2013, ApJ, 777, 173, doi: 10.1088/0004-637X/777/2/173
  * Beuther et al. (2020) Beuther, H., Wang, Y., Soler, J., et al. 2020, A&A, 638, A44, doi: 10.1051/0004-6361/202037950
  * Bialy et al. (2017) Bialy, S., Bihr, S., Beuther, H., Henning, T., & Sternberg, A. 2017, ApJ, 835, 126, doi: 10.3847/1538-4357/835/2/126
  * Bialy & Sternberg (2016) Bialy, S., & Sternberg, A. 2016, ApJ, 822, 83, doi: 10.3847/0004-637X/822/2/83
  * Bialy et al. (2015) Bialy, S., Sternberg, A., Lee, M.-Y., Le Petit, F., & Roueff, E. 2015, ApJ, 809, 122, doi: 10.1088/0004-637X/809/2/122
  * Bokeh Development Team (2018) Bokeh Development Team. 2018, Bokeh: Python library for interactive visualization.  https://bokeh.pydata.org/en/latest/
  * Busch et al. (2019) Busch, M. P., Allen, R. J., Engelke, P. D., et al. 2019, The Astrophysical Journal, 883, 158, doi: 10.3847/1538-4357/ab3a4b
  * Chen et al. (2019) Chen, B. Q., Huang, Y., Yuan, H. B., et al. 2019, MNRAS, 483, 4277, doi: 10.1093/mnras/sty3341
  * Chen et al. (2020) Chen, B. Q., Li, G. X., Yuan, H. B., et al. 2020, MNRAS, 493, 351, doi: 10.1093/mnras/staa235
  * Cox et al. (2016) Cox, N. L. J., Arzoumanian, D., André, P., et al. 2016, A&A, 590, A110, doi: 10.1051/0004-6361/201527068
  * Dame et al. (2001) Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792, doi: 10.1086/318388
  * Dolan & Mathieu (2002) Dolan, C. J., & Mathieu, R. D. 2002, AJ, 123, 387, doi: 10.1086/324631
  * Draine (2009) Draine, B. T. 2009, in Astronomical Society of the Pacific Conference Series, Vol. 414, Cosmic Dust - Near and Far, ed. T. Henning, E. Grün, & J. Steinacker, 453.  https://arxiv.org/abs/0903.1658
  * Duarte-Cabral & Dobbs (2016) Duarte-Cabral, A., & Dobbs, C. L. 2016, MNRAS, 458, 3667, doi: 10.1093/mnras/stw469
  * Dzib et al. (2018) Dzib, S. A., Loinard, L., Ortiz-León, G. N., Rodríguez, L. F., & Galli, P. A. B. 2018, ApJ, 867, 151, doi: 10.3847/1538-4357/aae687
  * Enßlin (2019) Enßlin, T. A. 2019, Annalen der Physik, 531, 1800127, doi: 10.1002/andp.201800127
  * Enßlin & Frommert (2011) Enßlin, T. A., & Frommert, M. 2011, Phys. Rev. D, 83, 105014, doi: 10.1103/PhysRevD.83.105014
  * Enßlin et al. (2009) Enßlin, T. A., Frommert, M., & Kitaura, F. S. 2009, Phys. Rev. D, 80, 105005, doi: 10.1103/PhysRevD.80.105005
  * Fiege & Pudritz (2000) Fiege, J. D., & Pudritz, R. E. 2000, MNRAS, 311, 85, doi: 10.1046/j.1365-8711.2000.03066.x
  * Field et al. (1969) Field, G. B., Goldsmith, D. W., & Habing, H. J. 1969, ApJ, 155, L149, doi: 10.1086/180324
  * Fukui et al. (2009) Fukui, Y., Kawamura, A., Wong, T., et al. 2009, ApJ, 705, 144, doi: 10.1088/0004-637X/705/1/144
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051