# CO mapping of Cygnus-X – volume density distribution

Jonah C. Baade Steward Observatory, The University of Arizona, Tucson AZ 85721
USA  Shuo Kong Steward Observatory, The University of Arizona, Tucson AZ 85721
USA  John H. Bieging Steward Observatory, The University of Arizona, Tucson AZ
85721 USA  Thomas Folkers  Steward Observatory, The University of Arizona,
Tucson AZ 85721 USA

###### Abstract

We present CO(2-1) and 13CO(2-1) maps of the Cygnus-X molecular cloud complex
using the 10m Heinrich Hertz Submillimeter Telescope (SMT). The maps cover the
southern portion of the complex which is strongly impacted by the feedback
from the Cygnus OB2 association. Combining CO(1-0) and 13CO(1-0) maps from the
Nobeyama 45m Cygnus-X CO Survey, we carry out a multi-transition molecular
line analysis with RADEX and derive the volume density of velocity-coherent
gas components. We select those components with a column density in the power-
law tail part of the column density probability distribution function (N-PDF)
and assemble their volume density into a volume density PDF ($\rho$-PDF). The
$\rho$-PDF exhibits a power-law shape in the range of 104.5 cm-3 $\lesssim
n_{\rm H_{2}}\lesssim$ 105.5 cm-3 with a fitted slope of $\alpha=-1.12\pm
0.05$. The slope is shallower than what is predicted by simulations of
rotationally supported structures or those undergoing gravitational collapse.
Applying the same analysis to synthetic observations with feedback may help
identify the cause of the shallow slope. The $\rho$-PDF provides another
useful benchmark for testing models of molecular cloud formation and
evolution.

ISM: clouds - ISM: kinematics and dynamics - ISM: molecules

††software: Python (Oliphant, 2007), SciPy (Jones et al., 2001–), Astropy
(Astropy Collaboration et al., 2013), Numpy (van der Walt et al., 2011),
Matplotlib (Hunter, 2007), SAOImageDS9 (Joye & Mandel, 2003), RADEX (van der
Tak et al., 2007), spectral-cube (Ginsburg et al., 2014)

##  1 Introduction

The density structure of a giant molecular cloud (GMC) is crucial to our
understanding of the star formation in the cloud. Because of shocks from
supersonic turbulence, a GMC typically contains local over-dense regions that
are prone to gravitational collapse. The collapsing loci will see the birth of
the next-generation young stars. Therefore, a census of the cloud density,
especially those exceeding the critical value for collapse, is extremely
useful for estimating star formation properties, including, e.g., the star
formation efficiency.

The density probability distribution function (hereafter $\rho$-PDF) is
particularly helpful for the census. Due to the multiplicative nature of the
density fluctuation, the $\rho$-PDF of a GMC with supersonic turbulence is
well approximated by a lognormal distribution (Vazquez-Semadeni, 1994; Padoan
et al., 1997) for relatively lower densities. The approximation remains valid
when considering chemistry (Glover et al., 2010; Federrath, 2013) and magnetic
fields (Lemaster & Stone, 2008; Padoan & Nordlund, 2011).

For higher densities that are gravitationally unstable, the $\rho$-PDF is
skewed to a power-law shape. Using numerical and analytical methods, Kritsuk
et al. (2011, hereafter K11) showed that the power-law relation originates
from the power-law density profile of spherically symmetric collapsing regions
(also see Girichidis et al., 2014; Jaupart & Chabrier, 2020). For a
spherically symmetric density profile $\rho\propto r^{-n}$, the $\rho$-PDF has
the form of $P(\rho)\propto\rho^{-3/n}$, i.e., a steeper collapsing density
profile gives rise to a shallower power-law $\rho$-PDF. For instance, a
collapsing isothermal sphere with a profile of $\rho\propto r^{-2}$ (e.g.,
Shu, 1977) would show a power-law $\rho$-PDF of $P(\rho)\propto\rho^{-1.5}$.
For a collapsing turbulent core with $\rho\sim r^{-1.5}$ (Murray & Chang,
2015), the $\rho$-PDF is $P(\rho)\propto\rho^{-2}$. K11 also observed a
density profile of $\rho\propto r^{-3}$ for rotationally supported cores at
higher densities ($n_{\rm H_{2}}>5\times 10^{9}$ cm-3) which indicated a
$\rho$-PDF of $P(\rho)\propto\rho^{-1}$. The PDF at the end of the simulation
in K11 showed two power-laws, one being $P(\rho)\propto\rho^{-1.67}$ at
intermediate densities, the other $P(\rho)\propto\rho^{-1}$ at higher
densities. These results were confirmed by the recent simulation work from
Khullar et al. (2021). Including magnetic fields, Collins et al. (2011) found
a $\rho$-PDF power-law slope of -1.64, which is consistent with the power-law
at intermediate densities in K11.

Observationally, the density structure is usually measured by the column
density (N-PDF) which is directly observable along each line-of-sight through,
e.g., dust emission/absorption and molecular line emission (e.g., Kainulainen
et al., 2009; Froebrich & Rowles, 2010; Schneider et al., 2013; Murase et al.,
2023). N-PDFs typically exhibit a lognormal profile at relatively low
extinctions with or without a power-law tail at higher extinctions. The
lognormal profile is believed to be a result of supersonic turbulence (e.g.,
Vazquez-Semadeni, 1994; Federrath et al., 2008; Ballesteros-Paredes et al.,
2011), while the power-law tail can result from pressure confinement (e.g.,
Passot & Vázquez-Semadeni, 1998) and/or gravitational collapse (e.g., Cho &
Kim, 2011; Ward et al., 2014).

Recently, Schneider et al. (2022, hereafter S22) measured N-PDF in 29 Galactic
regions using Herschel dust emission maps. They found that the clouds in
general were well fit with a sum of two lognormal profiles at lower densities
and a sum of two power-law tails at higher densities. The massive clouds tend
to have a shallower second power-law, consistent with the findings in K11
simulations and Schneider et al. (2015) observations. The shallower power-law
is probably due to physical mechanisms that hinder gas concentration,
including, e.g., rotation, magnetic fields, and (massive) stellar feedback.

In this paper, we aim to construct the $\rho$-PDF, instead of the N-PDF, of
the massive star-forming cloud Cygnus-X. While N-PDF is more straightforward
to measure in observations, $\rho$-PDF is more directly linked to the
underlying physics of the density structure. The $\rho$-PDF also provides a
pathway to estimating the star formation efficiency, especially in active
star-forming clouds that are dominated by gravity (Kainulainen et al., 2014).

To measure the $\rho$-PDF, we simply estimate the gas density ($n_{\rm
H_{2}}$) that is responsible for the molecular line emission through detailed
radiative transfer modeling. As will be introduced in §3, we use multi-
transition CO and 13CO lines for the modeling. The assumption is that the
molecular line emission is solely excited through the collisional interaction
with the partner H2. The advantage of the method is that the derived H2 volume
density $n_{\rm H_{2}}$ does not suffer the depletion problem which happens to
molecular column density measurements (e.g., Goodman et al., 2009), as the
freeze-out of H2 ice on dust grains is rather limited (Wakelam et al., 2017),
if there is any (see discussions in Seligman & Laughlin, 2020). The use of
molecular lines also alleviates the problem of projection effects in N-PDF
studies (e.g., Schneider et al., 2016; Wang et al., 2020; Murase et al.,
2023).

The Cygnus-X molecular cloud complex, at a distance of $\sim$1.4 kpc (Rygl et
al., 2012), is an active star-forming region with ample dense gas (total mass
$\sim 10^{6}$ M⊙, Schneider et al., 2006) and strong feedback from clusters of
O/B stars (Reipurth & Schneider, 2008; Wright et al., 2015). It provides an
unprecedented cosmic lab for the study of cloud density distribution in the
presence of the interaction between gravity and massive stellar feedback. As
introduced earlier, a $\rho$-PDF typically contains a power-law profile at
relatively higher densities, which is caused by gravitational collapse.
However, we lack quantitative knowledge about how feedback shapes the
$\rho$-PDF. Qualitatively, for instance, an expanding HII region sweeps the
vicinity of the driving massive star and creates compressed, high-density gas,
potentially flattening the power-law. Stutz & Kainulainen (2015); Schneider et
al. (2015) observed shallower power-law PDF in clouds with stellar feedback.
In this paper, with the measurement of the $\rho$-PDF in Cygnus-X, we aim to
examine how the strong feedback from the O-stars (re-)shapes the density
distribution. We will focus on the power-law part of the $\rho$-PDF at
relatively high densities.

In the following, we first describe our data collection in §2, including our
newly-acquired maps of CO(2-1) and 13CO(2-1) using the Heinrich Hertz
Submillimeter Telescope (SMT) and the publicly available Nobeyama 45m Radio
Observatory (hereafter NRO45) maps of CO(1-0) and 13CO(1-0). Then, in §3 we
describe in detail our multi-line fitting methods for deriving the gas volume
density. Next, we report our findings in §4. Finally, we discuss and conclude
in §5 and §6, respectively.

##  2 Observations and Data Reduction

###  2.1 SMT

Figure 1: Spitzer 8 $\mu$m image of Cygnus-X from the Spitzer Legacy Survey of
the Cygnus-X Region (Beerer et al., 2010; Kryukova et al., 2014). The white
polygon shows the NRO45 mapping area from the Nobeyama 45m Cygnus-X CO Survey
(Yamagishi et al., 2018; Takekoshi et al., 2019). The green polygon shows the
SMT mapping area. Our analysis is done in the overlapping region.

The CO(2-1) and 13CO(2-1) lines were observed with the Heinrich Hertz
Submillimeter Telescope (SMT) on Mt. Graham, Arizona, at an elevation of 3200
m. The facility is operated by the Arizona Radio Observatories (ARO), a
division of Steward Observatory at the University of Arizona. The basic flow
of observation and data reduction follows that described in Bieging et al.
(2010, 2014, 2016). In the following, we comment on several key aspects.

We map the southern part of the Cygnus-X complex from November 2019 to June
2021 with a beam size of 32″. Figure 1 shows the mapped region in the green
polygon. The region is divided into multiple 10′$\times$10′ rectangles, each
being mapped in an on-the-fly (OTF) session for about 1.5 hr. We use the dual-
polarization ALMA Band 6 prototype receiver with the spectrometer that
simultaneously observes CO(2-1) in the upper side-band and 13CO(2-1) in the
lower side-band. The data were calibrated and processed as described in
Bieging et al. (2014). Telescope pointing was checked and corrected as
necessary approximately every few hours. We use the quarter-MHz filter bank
for 13CO(2-1) which gives a spectral resolution of 0.34 km s-1. Because of the
broad line-width of the CO(2-1) line, we use the MHz filter bank for the line
to get the baseline fitting. Typical rms noise values per pixel were 0.09 K in
main-beam brightness temperature per 1.3 km s-1 for CO(2-1) and 0.10 K per
0.34 km s-1 for 13CO(2-1).

###  2.2 Nobeyama 45m data

In this paper, we make use of the Nobeyama 45m Cygnus-X survey data from
Yamagishi et al. (2018); Takekoshi et al. (2019). This survey utilized the
four-beam receiver FOREST on the NRO45 telescope to carry out a multi-line
mapping of the Cygnus-X North and South regions (see the white polygon in
Figure 1). The survey data is publicly available here. For our study, we
acquired the cubes for the CO(1-0) and 13CO(1-0) transitions. Out of the
several data products, we selected cubes with a beamsize of $23\arcsec~{}$, a
spectral resolution of 0.25 km s-1 , and a velocity range of $-40$ to 40 km
s-1. The choice gives the best match with our SMT data and covers the velocity
range of Cygnus-X (Schneider et al., 2006). The median rms noises for the
CO(1-0) and 13CO(1-0) data are 2.37 K and 0.95 K, respectively. Using Gaussian
kernel methods from the spectral-cube and FITS_tools Python packages, we
smoothed the NRO45 cubes to an effective angular resolution of $32\arcsec~{}$
and a spectral resolution of 0.34 km s-1 , to match the SMT 13CO(2-1) map. We
then reprojected them to the same spatial and spectral coordinates as the SMT
map, with a final pixel size of $10\arcsec~{}$ and a velocity range of $-15.9$
to 25.9 km s-1 .

##  3 Methods

We use the radiative transfer program RADEX (van der Tak et al., 2007) to
model the observation and estimate the gas volume density $n_{\rm H_{2}}$.
RADEX takes four parameters relevant to our study in its input: the kinetic
temperature of the gas ($T_{\rm kin}$), the volume density of the collision
partner ($n_{\rm H_{2}}$), the molecular column density ($N_{\rm{}^{13}CO}$),
and the FWHM line width of the observed velocity component.

For this pilot study, our strategy is to focus on gas parcels that are in
Local Thermodynamic Equilibrium (LTE) so that the volume density $n_{\rm
H_{2}}$ is the only unknown parameter for RADEX. Typically, LTE gas tends to
be relatively dense ($n_{\rm H_{2}}\gtrsim 10^{4}$ cm-3) with optically thick
CO lines. We analyze the aligned SMT and NRO45 cubes pixel by pixel, fitting
13CO lines for each. We can approximate the gas kinetic temperature with the
CO line excitation temperature. Then, the 13CO column density can be estimated
with an LTE equation for the relatively optically thin 13CO lines. The 13CO
line width can be obtained through line fitting. With this information,
$n_{\rm H_{2}}$ becomes the only free parameter, and by varying it we can find
a best-fit volume density for the observed emission. The complete procedure is
as the following:

  1. 1.

We use the astrodendro package (Rosolowsky et al., 2008) to find all velocity
components in CO(1-0), CO(2-1), 13CO(1-0), and 13CO(2-1).

  2. 2.

Using the astrodendro results as an initial guess, fit all four lines with
Gaussian profiles and get the peak intensity and the velocity dispersion for
each velocity component.

  3. 3.

Cross-match the velocity components between the four lines and keep the
components that exist in all of them. Such overlaps between components that
emit in all four lines are hereafter referred to as sets of components.

  4. 4.

Assuming an optically thick condition, compute the excitation temperature
$T_{\rm ex}$ for both CO(1-0) and CO(2-1) and keep the sets whose $T_{\rm ex}$
agrees within 30%.

  5. 5.

Compute 13CO column density N(13CO) with the CO $T_{\rm ex}$ assuming LTE. We
keep the sets whose N(13CO) are greater than $7\times 10^{15}\textrm{
cm}^{-2}$ to make sure we are in the power-law regime based on the N-PDF in
Schneider et al. (2016). We chose to focus on this power-law regime because
our method is less sensitive to low-density gas.

  6. 6.

With the derived $T_{\rm ex}$, N(13CO), and 13CO line dispersion, carry out
RADEX modeling (van der Tak et al., 2007) of the 13CO line emission with an
initial guess of $n_{\rm H_{2}}$.

  7. 7.

Iterate $n_{\rm H_{2}}$ to best match the 13CO(1-0) and 13CO(2-1) line
emission between RADEX and observations.

  8. 8.

Confirm the validity of LTE by checking the $T_{\rm ex}$ between RADEX and
observations.

Below, we expand on these steps in more details.

###  3.1 Fitting and selecting components

Figure 2: Plots illustrating Gaussian fit results for the pixel at (RA, DEC) =
(20h35m26s, 38°27′31″). The line data is shown as colored curves while the
black dashed line shows the fitted Gaussians. Our program assigned the peaks
centered at 2 and 14 km s-1 into two sets.

In each pixel, spectra in any or all of the four transitions could have
multiple peaks centered at different LSR velocities, assumed to correspond to
different dense regions along the line of sight. We isolated these peaks using
the astrodendro package (Step 1), using an absolute threshold of $3\sigma$
where $\sigma$ is the pixel-based noise. Each molecular line in each pixel had
its own value of $\sigma$, calculated from velocity ranges where the
corresponding cubes did not have any apparent bright regions. These ranges
were (21, 26) km s-1 for CO(1-0); (36, 100) km s-1 for CO(2-1); (-16, -11) and
(18, 26) km s-1 for 13CO(1-0); and (-16, -11) and (18, 26) km s-1 for
13CO(2-1). Pixels with any lines where $\sigma>2.5$ K, which appeared only
around the edges of the data cubes, were discarded. We use a relative
threshold of $1.5\sigma$ above the surroundings in velocity space, and a
minimum peak width of 3 spectral pixels (2 for the CO(2-1) map, which had a
factor of 4 lower spectral resolution) for astrodendro.

We then used the curve_fit method from the scipy.optimize package to fit a sum
of Gaussian components to each spectrum (Step 2), using the peaks identified
by astrodendro to provide initial values for the fit and determine the number
of peaks to fit. Figure 2 shows an example of the fitting. If any spectrum out
of the four needed more than five components to be adequately fit, we
discarded that pixel, under the assumption that the data there was too noisy
for a reliable analysis. The peak-fitting approach sometimes fits a low, wide
component to the data rather than a small but strong component to match a
peak. We excluded these components by discarding any Gaussians that did not
reach $3\sigma$ above the pixel’s noise level, for consistency with the
results from astrodendro.

Subsequent steps in our data analysis required peak temperatures from all four
transitions, so we removed Gaussian components that appeared in only some
lines but not in others at the same velocities (Step 3). To do this, we
constructed a simple hierarchy in which every Gaussian component in the
13CO(2-1) line was matched to the closest component in 13CO(1-0), then to
CO(2-1) and then CO(1-0). Using the standard velocity deviation $\sigma_{v}$
of each Gaussian in sequence, we counted pairs of components as aligned if the
”later” peak was within 2$\sigma_{v}$ in velocity of the ”earlier”. We also
excluded all components that shared matches with each other, to avoid
situations where a pair of components in one line appeared fully inside a
single component in another. After all of these exclusions, we were left with
only sets of nested, solitary components at similar velocities, assumed to
correspond to a single dense region each. We counted 195,344 sets across
175,592 pixels before further exclusions.

###  3.2 Physical conditions

####  3.2.1 Excitation temperature

Each set had distinct components in CO(2-1) and CO(1-0). Using their peak
temperatures, we estimated two CO excitation temperatures $T_{\rm ex}$ (Step
4), assuming $\tau\gg 1$, with the following equation (e.g., Eq. 3 in Kong et
al., 2015):

$T_{\rm ex}=\frac{h\nu/k}{\ln\left(1+\frac{h\nu/k}{T_{\rm mb,p}+J_{\nu}(T_{\rm
bg})}\right)},$ (1)

where $T_{\rm mb,p}$ is the peak main-beam temperature, $T_{\rm bg}$ is the
cosmic microwave background temperature of 2.73 K, and $J_{\nu}(T_{\rm
bg})=\frac{h\nu/k}{\left(e^{h\nu/kT_{\rm bg}}-1\right)}$ is the Rayleigh-Jeans
equivalent temperature of a blackbody of temperature $T_{\rm bg}$. To limit
our modeling to sets that were assumed close to thermal equilibrium and thus
amenable to an LTE approximation, we included only sets for which the $T_{\rm
ex}$ ratio of CO(2-1) and CO(1-0) was in the range 0.7-1.3. We counted 173,188
sets after this restriction.

####  3.2.2 13CO column density

We used the 13CO(2-1) line to compute the 13CO column density (Step 5),
assuming the LTE condition. The J=2-1 transition has the advantage of not
being sensitive to the excitation temperature between 10-40 K (see discussions
in Ginsburg et al., 2011; Bieging et al., 2016). We used an LTE equation to
estimate $N_{\rm{}^{13}CO}$ for each set (e.g., Eqs. 4, 5 in Kong et al.,
2015):

$\tau=-\ln\left(1-\frac{\frac{T_{\rm
mb}}{h\nu/k}}{\frac{1}{e^{\frac{h\nu}{kT_{\rm
ex}}}-1}-\frac{1}{e^{\frac{h\nu}{kT_{\rm bg}}}-1}}\right)$ (2)
$\begin{split}&N=\frac{8\pi k\nu^{2}}{hc^{3}A_{2\rightarrow
1}}\frac{e^{\frac{h\nu}{kT_{\rm bg}}}-1}{e^{\frac{h\nu}{kT_{\rm
bg}}}-e^{\frac{h\nu}{kT_{\rm
ex}}}}\frac{Q}{(2J+1)e^{\frac{-hB_{0}J(J+1)}{kT_{\rm ex}}}}\\\
&\times\frac{\int\tau(v)dv}{\int\left[1-e^{-\tau(v)}\right]dv}\int T_{\rm
mb}dv~{}({\rm cm}^{-2}),\end{split}$ (3)

where $A_{2\rightarrow 1}$ is the Einstein coefficient for the J=2-1
transition, and $B_{0}$ is the rotation constant of 13CO. The excitation
temperature $T_{\rm ex}$ used in this calculation was taken from the
calculated value for the CO(2-1) transition (if it was close enough to that
for the J=1-0 transition in Step 4). The integrals over the main beam
temperature $T_{\rm mb}$ and optical depth $\tau(v)$ were solved using the
Gaussian fits to the 13CO(2-1) spectrum, as integrating over the spectrum
itself was problematic when peaks overlapped. The partition function Q is
estimated by (e.g., Eq. 6 in Kong et al., 2015)

$Q\equiv\Sigma_{J}(2J+1)e^{\frac{-hB_{0}J(J+1)}{kT_{\rm
ex}}}\simeq\frac{kT_{\rm ex}}{hB_{0}}+\frac{1}{3}$ (4)

which is accurate to 10 percent for $T_{\rm ex}$ $>$ 5 K. See Mangum & Shirley
(2015) for a comprehensive derivation of the column density calculation.

We restricted our final analysis to sets with a column density
$N_{\rm{}^{13}CO}\geq 7\times 10^{15}\textrm{ cm}^{-2}$, equivalent to $N_{\rm
H_{2}}\gtrsim 5\times 10^{21}\textrm{ cm}^{-2}$ (assuming a H2-to-13CO ratio
of $\sim 7\times 10^{5}$, Wilson, 1999). These column densities are in
accordance with the regime Schneider et al. (2016) and S22 identified for a
power-law tail in the N-PDF.

####  3.2.3 FWHM line width

Each set had a single Gaussian component in the 13CO(2-1) transition. We used
this component to calculate a FWHM line width (from Step 2) to enter into
RADEX.

####  3.2.4 Kinetic temperature

Under the assumption that the gas in the sets of interest was in LTE, we set
$T_{\rm kin}$ for the RADEX model equal to the CO(2-1) excitation temperature
(from Step 4).

###  3.3 RADEX fitting

We aimed to constrain the $\text{H}_{2}$ volume density for a set by
reproducing the 13CO line intensities (Step 6). In RADEX, the background
radiation temperature was set to 2.73 K, and only $\text{H}_{2}$ was entered
as a collision partner.

We configured RADEX with the default escape probabilities for a uniform
sphere. As a result, we are essentially estimating the line-of-sight weighted
average of the volume density of the gas parcel that contributed to an
observed set of peaks. For a parcel with a monotonic density profile, our
estimation samples the volume density in a shell of the parcel. Since we are
sampling the volume density at different locations across Cygnus-X, the final
$\rho$-PDF is considered a reasonable representation of the underlying
$\rho$-PDF of the entire cloud at the scale of the beam size (32″,
corresponding to 0.22 pc at a distance of 1.4 kpc).

Figure 3: Plots illustrating the $n_{\rm H_{2}}$ fitting process. The left
panel is for the set at the pixel at (RA, DEC) = (20h37m59s, 39°55′01″). The
solid red and blue lines indicate the peak $T_{\rm mb}$ that RADEX outputs for
each value of $n_{\rm H_{2}}$, while the dashed lines are at the corresponding
peak $T_{\rm mb}$ observed along the line of sight. The dotted green line, at
$\textrm{n}=10^{4.5}~{}$cm-3, marks the best fit density for this set, where
the $\chi^{2}$, shown for each density in solid green and labeled on the
right, is minimized. The right panel is for the set at the pixel at (RA, DEC)
= (20h34m37s, 40°30′21″). The best fit density from RADEX fitting was
$\textrm{n}=10^{6}~{}$cm-3, which was the highest input in our grid.

Because the peak radiation temperatures of 13CO(1-0) and 13CO(2-1) were known
for each set, while $n_{\rm H_{2}}$ was unknown, we parameterized over H2
volume density and created a grid of identical inputs with only the volume
density changed. We tested a range for $n_{\rm H_{2}}$ of $10^{2}$-$10^{6}$
cm-3, with a logarithmic step size of 0.1 dex. We attributed a best-fit volume
density (Figure 3) to a set by using the value of $n_{\rm H_{2}}$ which
produced the smallest $\chi^{2}$ value between the observed and modeled peak
radiation temperatures in both 13CO lines (Step 7). The $\chi^{2}$ is computed
as

$\chi^{2}=\Sigma\frac{(O-E)^{2}}{\sigma^{2}},$ (5)

where O stands for observation and E stands for model. $\sigma$ is the
spectral noise. Currently, we are limited by just two 13CO line transitions.
Adding more transitions in the future will give stronger constraints. While
the vast majority ${\rm(\sim 88\%)}$ of sets had a local minimum $\chi^{2}$
within our test range, the remainder showed $\chi^{2}$ approaching an
asymptote as we approached the maximum test volume density of $10^{6}$ cm-3.
In the final output, these sets were placed in the $10^{6}$ cm-3 group. After
fitting, we counted 171,354 sets across 156,371 pixels (with the discrepancy
amounting to 1,834 sets that failed the integration for column density). Most
of the sets were poorly fit: only 5,480 of them reached $\chi^{2}<10$ at any
point (including 1,079 which did not have a local minimum), across 5,467
pixels. However, this subset covers the main spatial extent of the region we
have observed, sampling primarily from filaments and dense clouds. This
provides reason to believe that the H2 density distribution we measure
reflects the true distribution throughout the denser parts of Cygnus-X South.

Our choice of cutoff at $\chi^{2}<10$ was arbitrary; it represents the best
${\rm\sim 17\%}$ of fits meeting the column density requirement, but there is
no reason to suspect this is the best requirement for closeness of fit to use.
We repeated the final step of the procedure for different $\chi^{2}$ cutoff
values, spaced by powers of two from 2.5 to 1,280. Further increasing the
$\chi^{2}$ cutoff value did not increase the number of sets with
$N_{\rm{}^{13}CO}\geq 7\times 10^{15}\textrm{ cm}^{-2}$ (required by criterion
5 in §3). The results are shown in Figure 8 and discussed in Section 4.3.

##  4 Results

###  4.1 Global spectra

Figure 4: Averaged spectra for all four molecular lines across the entire
region of mutual overlap between cubes. Pixels with $\sigma>2.5$ K, which
appeared only around the edges of the data cubes, were discarded.

Figure 4 shows spectra of the four lines we studied, averaged across the area
where all four of these lines were present. All four lines have a strong
component at or near 0 km s-1 and a weaker component near 7 km s-1 , but with
significant emission across a much wider band of velocities, generally
spanning the range $-10$ to 15 km s-1 . Each of the 13CO lines peaks at about
one-fifth the main beam temperature of its counterpart in CO at the same
transition.

###  4.2 Integrated intensity and first-moment maps

Figure 5: Integrated intensity maps for the CO(2-1) and 13CO(2-1) data cubes.
The integrated velocity range is $-16$ to 26 km s-1 . The color bar has units
of K km s-1 . Figure 6: First-moment maps for the CO(2-1) and 13CO(2-1) data
cubes. The integrated velocity range is $-16$ to 26 km s-1 . Velocities in
pixels with a peak intensity less than $3\sigma$ above the noise level were
masked out. The color bar has units of km s-1 .

In Figure 5, we show maps of the integrated line intensity for CO(2-1) and
13CO(2-1), both from SMT observations. In general, the integrated intensity
traces the same features in both maps, and it is about five times higher in CO
than in 13CO everywhere.

Figure 6 shows maps of the first moment, the intensity-weighted mean velocity
of the gas along the line of sight, for CO(2-1) and 13CO(2-1). Both maps show
large regions where the mean velocity is zero or slightly negative, and other
regions where it is in the vicinity of 10 km s-1 , which is consistent with
the averaged spectra in Figure 4. Some features are visible in both integrated
intensity and first-moment maps: in particular, the filament at the bottom of
the images is visible as a bright spot in both integrated intensity maps and
as a negative-velocity region in the first moment maps. The map in 13CO is
missing many pixels because the emission intensity in those pixels was never
more than $3\sigma$ above the noise level. This is a consequence of emission
in 13CO being weaker overall.

###  4.3 $\rho$-PDF

Figure 7: $\rho$-PDF derived from best-fit volume densities for each set. Only
those with a best-fit $\chi^{2}<10$ were included. The blue line displays a
linear-regression power law fit to the data between log($n_{\rm H_{2}}$) = 4.5
and 5.5.

In Figure 7, we show the $\rho$-PDF of the best-fit volume densities for the
5,480 sets that reached $\chi^{2}<10$. The plot, in logarithmic scales, covers
a range in log($n_{\rm H_{2}}$) from 3.6–6.0, in logarithmic increments of 0.1
dex; no well-fit sets had best-fit volume densities less than $10^{3.6}$ cm-3.
A large number of sets whose $\chi^{2}$ values did not reach a local minimum
within our test range, and were therefore placed in the $10^{6}~{}$cm-3 group,
can be attributed to saturation while performing the fitting with RADEX. For
example, in Figure 3 right panel, the calculated peak $T_{\rm mb}$ for each
transition gradually approached the observed peak $T_{\rm mb}$, but stopped
rising before they could match. This could indicate that the density is
actually over $10^{6}~{}$cm-3 in that area. For instance, Schneider et al.
(2016) detected high-density tracers in the cloud. In our case, a higher
transition line would break the degeneracy.

The value of the probability distribution function P at each bin was
calculated from the ratio of the number of sets in it divided by the total (in
this case 5,480), and then divided by the bin width of 0.1 accordingly to make
$\int_{3.6}^{6.0}P~{}\mathrm{d}\log(n_{\rm H_{2}})=1$. We calculated a line of
best fit with linear regression, using the binned volume densities as the
x-coordinate in $\log(n_{\rm H_{2}})$-space, and the bin heights as the
y-coordinate in $\log(P)$-space. This is consistent with the corresponding
column-density log-log fit between $A_{V}\propto N_{\rm H_{2}}$ and $P$ for
power-law tails seen in e.g. S22.

At low densities, the distribution appears lognormal, but the shape is more
likely a result of a lack of detected sets of low enough intensity to
correspond to these densities. At high densities, the distribution appears to
be a straight line in the log-log plot, indicating a power-law relation in the
$\rho$-PDF ($P\propto\rho^{\alpha}$). A linear-regression fit between
log($n_{\rm H_{2}}$) = 4.5 and 5.5 gives a power-law slope of
$\alpha=-1.122\pm 0.048$. Above these densities, the frequency rises again
until our grid’s maximum volume density of $10^{6.0}$ cm-3, where the greatest
number of sets matched.

We produced additional plots with the same procedure but different values for
the largest permitted local minimum $\chi^{2}$, from 2.5 to 1,280, with the
latter limit just sufficient to encompass all of the sets that met the column
density requirement $N_{\rm{}^{13}CO}\geq 7\times 10^{15}\textrm{ cm}^{-2}$.
All of these plots showed some power-law slope between log($n_{\rm H_{2}}$) =
4.5 and 5.5, with a value between $\alpha$ = $-1.0$ and $-1.25$ (Figure 8).
Uncertainty in the slope is much greater for cutoffs below our main value of
$\chi^{2}=10$ but somewhat less for cutoffs above it, with the highest cutoff,
including all sets of sufficient column density, having a slope of
$\alpha=-1.155\pm 0.028$.

If there were a strong power-law in the true $\rho$-PDF of Cygnus-X, we might
expect the slope of the power-law fit to approach this value more closely for
more restrictive cutoffs of $\chi^{2}$, as we would be sampling from the sets
for which our modeling of volume density was most accurate. While there is no
clear trend towards a specific value in Figure 8, our power-law fit slopes
stay within the range between $\alpha$ = $-1.0$ and $-1.25$, which suggests
that the approximate value is still relevant to our conclusions. Running our
procedure on synthetic data with less natural noise could help us determine
whether there is a systematic bias towards this range of slopes.

###  4.4 $T_{\rm ex}$ Consistency

To confirm the accuracy of our $T_{\rm ex}$ estimates for CO(2-1) and CO(1-0)
for each set that we made in Section 3.2, we checked them against the $T_{\rm
ex}$ values for 13CO(2-1) and 13CO(1-0) that RADEX produced alongside its
calculated peak radiation temperatures for the best-fit volume density (Step
8). Under LTE conditions, all four excitation temperatures should be equal to
each other and to the gas kinetic temperature in the cloud.

The results are shown in Figure 9 for all 4,401 sets with minimum
$\chi^{2}<10$ but that did not reach our maximum density of $10^{6}~{}$cm-3.
(The plot excludes those that reached the highest density since the degeneracy
caused an anomalous fit; those sets tend to have error close to zero.) Almost
all ${\rm(\sim 87\%)}$ of RADEX’s 13CO excitation temperatures for these sets
were higher in the J=1-0 line, and all were lower in the J=2-1 line, than our
LTE estimates. About 92% of the sets had errors below 0.3 in both lines, being
consistent with our 30% $T_{\rm ex}$ error threshold for including them in the
first place.

Figure 8: Slopes of the linear regression power-law fits with varying cutoffs
for the largest permitted minimum $\chi^{2}$. In all cases, the power-law fit
was between log($n_{\rm H_{2}}$) = 4.5 and 5.5, and only sets with column
density $N_{\rm{}^{13}CO}\geq 7\times 10^{15}\textrm{ cm}^{-2}$ were included.
Numbers on the plot indicate the number of sets falling within that cutoff,
with 31,343 being the total number with the required column density. Figure 9:
Histogram of relative errors between $T_{\rm ex}$ calculated for CO transition
lines and those estimated for the corresponding 13CO lines by RADEX. The
horizontal scale goes up to 0.85 to include the highest outlier in the J=1-0
plot (not visible).

Because we assumed LTE when calculating $T_{\rm ex}$ in the first place for
the CO J=1-0 and J=2-1 transitions, the errors between our computations for CO
and RADEX’s for 13CO are understandable. RADEX systematically returns lower
$T_{\rm ex}$ for the J=2-1 transition than our calculations because under the
non-LTE conditions it simulated, that line is not excited as much, while it
usually gives higher $T_{\rm ex}$ for the J=1-0 line because stimulated
emission at high densities contributes to $T_{\rm ex}$ above what LTE would
predict from gas kinetic temperature alone.

##  5 Discussion

The shallowest $\rho$-PDF slope in Kritsuk et al. (2011, hereafter K11) was
$-1$, which was caused by rotationally supported cores111Note that the
simulations in K11 were also limited by resolutions. at a very high density of
order $n_{\rm H_{2}}\sim 10^{9}$ cm-3. Such a density is not observed in our
$\rho$-PDF in Figure 7 and could not be identified with these CO lines. In
addition, the SMT angular resolution of $\sim$32″ corresponds to a spatial
scale of $\sim$0.22 pc (assuming a distance of 1.4 kpc). Such a physical scale
is not able to resolve the sub-structure of a dense core. Thus, the shallow
slope in our measurement ($-1.12\pm 0.05$) is not explained by core rotation.

Stutz & Kainulainen (2015), when discussing flattened N-PDFs in the Orion A
molecular cloud, suggested that feedback from O-stars in the cloud could be
responsible (also see Federrath & Klessen, 2013). Similarly, Russeil et al.
(2013); Tremblin et al. (2014); Schneider et al. (2015); Pokhrel et al. (2016)
saw hints of shallow PDFs that were caused by stellar feedback. Cygnus X is a
region with its own O-star feedback mainly from the Cygnus OB2 Association,
which is to the north of the area we studied. So we tentatively attribute the
shallow power-law to the strong feedback environment.

Recently, S22 measured the N-PDF in Cygnus North and Cygnus South regions.
Both regions were shown to have two power-law tails in the N-PDF. Cygnus North
has a relatively shallow first power-law tail (slope $-1.83$), and its second
power-law is even shallower (slope $-1.40$), with the transition being at
$A_{V}\sim 80$ mag. The first N-PDF power-law tail for Cygnus South has a
slope of $-2.37$ which is the median value for their high-mass star-forming
sample, and the second power-law tail is steeper with a slope of $-2.66$, with
the transition being at $A_{V}\sim 37$ mag. Our mapped region has more overlap
with their Cygnus South region. Due to the limited dynamic range, our fitted
$\rho$-PDF (Figure 7) only shows one power-law. So it is not clear if there is
a second power-law at higher densities as shown by the N-PDF in S22. If we
still follow the K11 scheme, the slope of $-1.12$ for the $\rho$-PDF should
correspond to a slope of $-1.19$ for the N-PDF, which is shallower than the
measurement from S22. While the comparison between $\rho$-PDF and N-PDF is not
straightforward, it seems our measurement indicates a shallower PDF power-law
slope than that from S22.

The apparent difference could be due to the different sampling methods (§3).
We are distinguishing different ”sets” of aligned velocity components along
the same line-of-sight while S22 included the total column in the Herschel
map. We are only including sets in LTE ($T_{\rm ex}$ agrees within 30% for two
12CO transitions) which intrinsically prefers denser gas. We are keeping
preferentially higher column densities to ensure we are in the power-law
regime. Considering these limits, we are effectively sampling deeply embedded
denser gas that is temporarily shielded from the external radiation. But this
gas probably feels the compression due to the high pressure from hot gas
outside. This special sampling probably leads to the relatively shallower
$\rho$-PDF. In the future, a comparison with synthetic observations of gas
that is strongly impacted by feedback could help distinguish between such
processes.

Our finding of a power-law $\rho$-PDF generally agrees with S22. However,
Murase et al. (2023) did not find power-law profiles in 13CO N-PDF in the
majority of their sample. Note that they broke the Cygnus-X complex into
multiple individual clouds in the position-position-velocity space. While
Murase et al. (2023) considered the optical depth effect for the J=1-0
transition and argued that it was a minor effect, they did not clearly show
how they address the depletion effect of CO. Since they assumed LTE with
median $T_{\rm ex}\sim$10-16K, the depletion of 13CO could impact the column
density derivation and thus the N-PDF, assuming the dust temperature is
similar to $T_{\rm ex}$ under LTE. For instance, Lewis et al. (2021) showed
noticeable depletion when dust temperature was below 20 K, with the depletion
factor increasing from 2 to 10 at a temperature of 12 K. If higher column
density regions have lower temperature, the correction for depletion should be
larger at high column density. Indeed, Schneider et al. (2016) found that CO
depletion probably caused the N-PDF falloff at $A_{V}\gtrsim 40$ mag. They
concluded that N-PDFs based on molecular lines were not well constrained.

##  6 Summary and Conclusion

We have presented OTF maps of the Cygnus-X molecular cloud complex in the
J=2-1 transition of both CO and 13CO using the ARO SMT telescope. The maps
extend from 20h24m to 20h39m in RA and from 38° to 40°43′ in DEC, covering an
area of $\sim$10 square degrees. The mapped area overlaps with the Nobeyama
45m Cygnus-X CO Survey in the southern portion of the complex in which we
carry out a multi-transition analysis for CO(1-0), 13CO(1-0), CO(2-1), and
13CO(2-1). Through the analysis, we have derived physical properties of the
molecular gas, especially the gas volume density.

We use astrodendro to identify velocity-coherent gas components (sets in each
pixel) in all four lines. For each set we fit a Gaussian profile to each line
spectrum and obtain the peak intensity and line dispersion. We derive the
excitation temperature $T_{\rm ex}$ for each CO line assuming optically thick
conditions. We select sets whose $T_{\rm ex}$ agree within 30% between the two
CO lines so that LTE is a reasonable assumption for the sets. We then derive
the 13CO column density for each set. These physical properties are fed to the
RADEX program with which we iteratively find the best H2 volume density by
matching the intensities between the model and the observation for the two
13CO lines. We construct the volume density probability distribution function
($\rho$-PDF) which shows a single power-law tail between log($n_{\rm H_{2}}$)=
4.5 and 5.5 with a slope of $\sim-1.12\pm 0.05$. Higher transition lines are
needed to trace higher densities.

The $\rho$-PDF power-law slope is shallower than the theoretical values
between $-1.5$ and $-2.0$ in the context of self-similar gravitational
collapse (e.g., Kritsuk et al., 2011). Within the same context, the
corresponding column density PDF power-law tail is also shallower than the
measurements from Schneider et al. (2022) using Herschel maps. Possible causes
include the fact that Cygnus-X is strongly impacted by feedback from the OB2
association (similar to what was found in Stutz & Kainulainen, 2015) and we
are specifically sampling deeply embedded dense gas. Future comparison with
synthetic observations with feedback may help clarify these causes. Our
observational study of $\rho$-PDF provides another useful test for theories.

We thank the anonymous referee for a constructive report. We acknowledge
fruitful discussions with Yancy Shirley and Samantha Scibelli. The Heinrich
Hertz Submillimeter Telescope is operated by the Arizona Radio Observatory,
which is part of Steward Observatory at The University of Arizona. We thank
Patrick Fimbres, Christian Holmstedt, Robert Moulton, Blythe Guvenen, and all
the observatory staff members for helping with the data acquisition and
maintaining a smooth telescope operation.
SMT,NRO45

## References

  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Ballesteros-Paredes et al. (2011) Ballesteros-Paredes, J., Vázquez-Semadeni, E., Gazol, A., et al. 2011, MNRAS, 416, 1436, doi: 10.1111/j.1365-2966.2011.19141.x
  * Beerer et al. (2010) Beerer, I. M., Koenig, X. P., Hora, J. L., et al. 2010, ApJ, 720, 679, doi: 10.1088/0004-637X/720/1/679
  * Bieging et al. (2016) Bieging, J. H., Patel, S., Peters, W. L., et al. 2016, ApJS, 226, 13, doi: 10.3847/0067-0049/226/1/13
  * Bieging et al. (2010) Bieging, J. H., Peters, W. L., & Kang, M. 2010, ApJS, 191, 232, doi: 10.1088/0067-0049/191/2/232
  * Bieging et al. (2014) Bieging, J. H., Revelle, M., & Peters, W. L. 2014, ApJS, 214, 7, doi: 10.1088/0067-0049/214/1/7
  * Cho & Kim (2011) Cho, W., & Kim, J. 2011, MNRAS, 410, L8, doi: 10.1111/j.1745-3933.2010.00968.x
  * Collins et al. (2011) Collins, D. C., Padoan, P., Norman, M. L., & Xu, H. 2011, ApJ, 731, 59, doi: 10.1088/0004-637X/731/1/59
  * Federrath (2013) Federrath, C. 2013, MNRAS, 436, 1245, doi: 10.1093/mnras/stt1644
  * Federrath & Klessen (2013) Federrath, C., & Klessen, R. S. 2013, ApJ, 763, 51, doi: 10.1088/0004-637X/763/1/51
  * Federrath et al. (2008) Federrath, C., Klessen, R. S., & Schmidt, W. 2008, ApJ, 688, L79, doi: 10.1086/595280
  * Froebrich & Rowles (2010) Froebrich, D., & Rowles, J. 2010, MNRAS, 406, 1350, doi: 10.1111/j.1365-2966.2010.16769.x
  * Ginsburg et al. (2011) Ginsburg, A., Bally, J., & Williams, J. P. 2011, MNRAS, 418, 2121, doi: 10.1111/j.1365-2966.2011.19279.x
  * Ginsburg et al. (2014) Ginsburg, A., Robitaille, T., Beaumont, C., & ZuHone, J. 2014, Release Candidate 2 - includes yt interop, v0.2-rc2, Zenodo, doi: 10.5281/zenodo.11485
  * Girichidis et al. (2014) Girichidis, P., Konstandin, L., Whitworth, A. P., & Klessen, R. S. 2014, ApJ, 781, 91, doi: 10.1088/0004-637X/781/2/91
  * Glover et al. (2010) Glover, S. C. O., Federrath, C., Mac Low, M. M., & Klessen, R. S. 2010, MNRAS, 404, 2, doi: 10.1111/j.1365-2966.2009.15718.x
  * Goodman et al. (2009) Goodman, A. A., Pineda, J. E., & Schnee, S. L. 2009, ApJ, 692, 91, doi: 10.1088/0004-637X/692/1/91
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jaupart & Chabrier (2020) Jaupart, E., & Chabrier, G. 2020, ApJ, 903, L2, doi: 10.3847/2041-8213/abbda8
  * Jones et al. (2001–) Jones, E., Oliphant, T., Peterson, P., et al. 2001–, SciPy: Open source scientific tools for Python.  http://www.scipy.org/
  * Joye & Mandel (2003) Joye, W. A., & Mandel, E. 2003, in Astronomical Society of the Pacific Conference Series, Vol. 295, Astronomical Data Analysis Software and Systems XII, ed. H. E. Payne, R. I. Jedrzejewski, & R. N. Hook, 489
  * Kainulainen et al. (2009) Kainulainen, J., Beuther, H., Henning, T., & Plume, R. 2009, A&A, 508, L35, doi: 10.1051/0004-6361/200913605
  * Kainulainen et al. (2014) Kainulainen, J., Federrath, C., & Henning, T. 2014, Science, 344, 183, doi: 10.1126/science.1248724
  * Khullar et al. (2021) Khullar, S., Federrath, C., Krumholz, M. R., & Matzner, C. D. 2021, MNRAS, 507, 4335, doi: 10.1093/mnras/stab1914
  * Kong et al. (2015) Kong, S., Lada, C. J., Lada, E. A., et al. 2015, ApJ, 805, 58, doi: 10.1088/0004-637X/805/1/58
  * Kritsuk et al. (2011) Kritsuk, A. G., Norman, M. L., & Wagner, R. 2011, ApJ, 727, L20, doi: 10.1088/2041-8205/727/1/L20
  * Kryukova et al. (2014) Kryukova, E., Megeath, S. T., Hora, J. L., et al. 2014, AJ, 148, 11, doi: 10.1088/0004-6256/148/1/11
  * Lemaster & Stone (2008) Lemaster, M. N., & Stone, J. M. 2008, ApJ, 682, L97, doi: 10.1086/590929
  * Lewis et al. (2021) Lewis, J. A., Lada, C. J., Bieging, J., et al. 2021, ApJ, 908, 76, doi: 10.3847/1538-4357/abc41f
  * Mangum & Shirley (2015) Mangum, J. G., & Shirley, Y. L. 2015, PASP, 127, 266, doi: 10.1086/680323
  * Murase et al. (2023) Murase, T., Handa, T., Matsusaka, R., et al. 2023, MNRAS, 523, 1373, doi: 10.1093/mnras/stad1451
  * Murray & Chang (2015) Murray, N., & Chang, P. 2015, ApJ, 804, 44, doi: 10.1088/0004-637X/804/1/44
  * Oliphant (2007) Oliphant, T. E. 2007, Computing in Science & Engineering, 9, 10, doi: 10.1109/MCSE.2007.58
  * Padoan & Nordlund (2011) Padoan, P., & Nordlund, Å. 2011, ApJ, 730, 40, doi: 10.1088/0004-637X/730/1/40
  * Padoan et al. (1997) Padoan, P., Nordlund, A., & Jones, B. J. T. 1997, MNRAS, 288, 145, doi: 10.1093/mnras/288.1.145
  * Passot & Vázquez-Semadeni (1998) Passot, T., & Vázquez-Semadeni, E. 1998, Phys. Rev. E, 58, 4501, doi: 10.1103/PhysRevE.58.4501
  * Pokhrel et al. (2016) Pokhrel, R., Gutermuth, R., Ali, B., et al. 2016, MNRAS, 461, 22, doi: 10.1093/mnras/stw1303
  * Reipurth & Schneider (2008) Reipurth, B., & Schneider, N. 2008, in Handbook of Star Forming Regions, Volume I, ed. B. Reipurth, Vol. 4, 36
  * Rosolowsky et al. (2008) Rosolowsky, E. W., Pineda, J. E., Kauffmann, J., & Goodman, A. A. 2008, ApJ, 679, 1338, doi: 10.1086/587685
  * Russeil et al. (2013) Russeil, D., Schneider, N., Anderson, L. D., et al. 2013, A&A, 554, A42, doi: 10.1051/0004-6361/201219971
  * Rygl et al. (2012) Rygl, K. L. J., Brunthaler, A., Sanna, A., et al. 2012, A&A, 539, A79, doi: 10.1051/0004-6361/201118211
  * Schneider et al. (2006) Schneider, N., Bontemps, S., Simon, R., et al. 2006, A&A, 458, 855, doi: 10.1051/0004-6361:20065088
  * Schneider et al. (2013) Schneider, N., André, P., Könyves, V., et al. 2013, ApJ, 766, L17, doi: 10.1088/2041-8205/766/2/L17
  * Schneider et al. (2015) Schneider, N., Bontemps, S., Girichidis, P., et al. 2015, MNRAS, 453, L41, doi: 10.1093/mnrasl/slv101
  * Schneider et al. (2016) Schneider, N., Bontemps, S., Motte, F., et al. 2016, A&A, 587, A74, doi: 10.1051/0004-6361/201527144
  * Schneider et al. (2022) Schneider, N., Ossenkopf-Okada, V., Clarke, S., et al. 2022, A&A, 666, A165, doi: 10.1051/0004-6361/202039610
  * Seligman & Laughlin (2020) Seligman, D., & Laughlin, G. 2020, ApJ, 896, L8, doi: 10.3847/2041-8213/ab963f
  * Shu (1977) Shu, F. H. 1977, ApJ, 214, 488, doi: 10.1086/155274
  * Stutz & Kainulainen (2015) Stutz, A. M., & Kainulainen, J. 2015, A&A, 577, L6, doi: 10.1051/0004-6361/201526243
  * Takekoshi et al. (2019) Takekoshi, T., Fujita, S., Nishimura, A., et al. 2019, ApJ, 883, 156, doi: 10.3847/1538-4357/ab3a55
  * Tremblin et al. (2014) Tremblin, P., Schneider, N., Minier, V., et al. 2014, A&A, 564, A106, doi: 10.1051/0004-6361/201322700
  * van der Tak et al. (2007) van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627, doi: 10.1051/0004-6361:20066820
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Vazquez-Semadeni (1994) Vazquez-Semadeni, E. 1994, ApJ, 423, 681, doi: 10.1086/173847
  * Wakelam et al. (2017) Wakelam, V., Bron, E., Cazaux, S., et al. 2017, Molecular Astrophysics, 9, 1, doi: 10.1016/j.molap.2017.11.001
  * Wang et al. (2020) Wang, Y., Beuther, H., Schneider, N., et al. 2020, A&A, 641, A53, doi: 10.1051/0004-6361/202037928
  * Ward et al. (2014) Ward, R. L., Wadsley, J., & Sills, A. 2014, MNRAS, 445, 1575, doi: 10.1093/mnras/stu1868
  * Wilson (1999) Wilson, T. L. 1999, Reports on Progress in Physics, 62, 143, doi: 10.1088/0034-4885/62/2/002
  * Wright et al. (2015) Wright, N. J., Drew, J. E., & Mohr-Smith, M. 2015, MNRAS, 449, 741, doi: 10.1093/mnras/stv323
  * Yamagishi et al. (2018) Yamagishi, M., Nishimura, A., Fujita, S., et al. 2018, ApJS, 235, 9, doi: 10.3847/1538-4365/aaab4b