CO Mapping of Cygnus-X—Volume Density Distribution

We present CO(2-1) and 13CO(2-1) maps of the Cygnus-X molecular cloud complex using the 10 m Heinrich Hertz Submillimeter Telescope. 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 45 m Cygnus-X CO Survey, we carry out a multitransition 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 (ρ-PDF). The ρ-PDF exhibits a power-law shape in the range of 104.5 cm−3 ≲nH2≲ 105.5 cm−3 with a fitted slope of α = −1.12 ± 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 ρ-PDF provides another useful benchmark for testing models of molecular cloud formation and evolution.


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 ρ-PDF) is particularly helpful for the census.Due to the multiplicative nature of the density fluctuation, the ρ-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 ρ-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 ρ ∝ r −n , the ρ-PDF has the form of P (ρ) ∝ ρ −3/n , i.e., a Baade steeper collapsing density profile gives rise to a shallower power-law ρ-PDF.For instance, a collapsing isothermal sphere with a profile of ρ ∝ r −2 (e.g., Shu 1977) would show a powerlaw ρ-PDF of P (ρ) ∝ ρ −1.5 .For a collapsing turbulent core with ρ ∼ r −1.5 (Murray & Chang 2015), the ρ-PDF is P (ρ) ∝ ρ −2 .K11 also observed a density profile of ρ ∝ r −3 for rotationally supported cores at higher densities (n H 2 > 5 × 10 9 cm −3 ) which indicated a ρ-PDF of P (ρ) ∝ ρ −1 .The PDF at the end of the simulation in K11 showed two power-laws, one being P (ρ) ∝ ρ −1.67 at intermediate densities, the other P (ρ) ∝ ρ −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 ρ-PDF power-law slope of -1.64, which is consistent with the power-law at intermediate densities in K11.
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 ρ-PDF, instead of the N-PDF, of the massive starforming cloud Cygnus-X.While N-PDF is more straightforward to measure in observations, ρ-PDF is more directly linked to the underlying physics of the density structure.The ρ-PDF also provides a pathway to estimating the star formation efficiency, especially in active starforming clouds that are dominated by gravity (Kainulainen et al. 2014).
To measure the ρ-PDF, we simply estimate the gas density (n 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 13 CO lines for the modeling.The assumption is that the molecular line emission is solely excited through the collisional interaction with the partner H 2 .The advantage of the method is that the derived H 2 volume density n 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 H 2 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 ∼1.4 kpc (Rygl et al. 2012), is an active star-forming region with ample dense gas (total mass ∼ 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 ρ-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 ρ-PDF.Qualitatively, for instance, an expanding H II 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 ρ-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 ρ-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 13 CO(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 13 CO(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.

OBSERVATIONS AND DATA REDUCTION
2.1.SMT The CO(2-1) and 13 CO(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. (2010Bieging et al. ( , 2014Bieging et al. ( , 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 ′ ×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 13 CO(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 13 CO(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 13 CO(2-1).

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 13 CO(1-0) transitions.Out of the several data products, we selected cubes with a beamsize of 23 ′′ , 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 13 CO(1-0) data are 2.37 K and 0.95 K, respectively.Using Gaussian  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.
kernel methods from the spectral-cube and FITS tools Python packages, we smoothed the NRO45 cubes to an effective angular resolution of 32 ′′ and a spectral resolution of 0.34 km s −1 , to match the SMT 13 CO(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 ′′ and a velocity range of −15.9 to 25.9 km s −1 .

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 H 2 .RADEX takes four parameters relevant to our study in its input: the kinetic temperature of the gas (T kin ), the volume density of the collision partner (n H 2 ), the molecular column density (N13 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 H 2 is the only unknown parameter for RADEX.Typically, LTE gas tends to be relatively dense (n H 2 ≳ 10 4 cm −3 ) with optically thick CO lines.We analyze the aligned SMT and NRO45 cubes pixel by pixel, fitting 13 CO lines for each.We can approximate the gas kinetic temperature with the CO line excitation temperature.Then, the 13 CO column density can be estimated with an LTE equation for the relatively optically thin 13 CO lines.The 13 CO line width can be obtained through line fitting.With this information, n 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. We use the astrodendro package (Rosolowsky et al. 2008) to find all velocity components in CO(1-0), CO(2-1), 13 CO(1-0), and 13 CO(2-1).
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. 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. Assuming an optically thick condition, compute the excitation temperature T ex for both CO(1-0) and CO(2-1) and keep the sets whose T ex agrees within 30%.
5. Compute 13 CO column density N( 13 CO) with the CO T ex assuming LTE.We keep the sets whose N( 13 CO) are greater than 7 × 10 15 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.With the derived T ex , N( 13 CO), and Below, we expand on these steps in more details.

Fitting and selecting components
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σ where σ is the pixel-based noise.Each molecular line in each pixel had its own value of σ, 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 13 CO(1-0); and (-16, -11) and (18, 26) km s −1 for 13 CO(2-1).Pixels with any lines where σ > 2.5 K, which appeared only around the edges of the data cubes, were discarded.We use a relative threshold of 1.5σ 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.optimizepackage 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σ 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 13 CO(2-1) line was matched to the closest component in 13 CO(1-0), then to CO(2-1) and then CO(1-0).Using the standard velocity deviation σ v of each Gaussian in sequence, we counted pairs of components as aligned if the "later" peak was within 2σ 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.

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 ex (Step 4), assuming τ ≫ 1, with the following equation (e.g., Eq. 3 in Kong et al. 2015): where T mb,p is the peak main-beam temperature, T bg is the cosmic microwave background temperature of 2.73 K, and J ν (T bg ) = hν/k (e hν/kT bg −1) is the Rayleigh-Jeans equivalent temperature of a blackbody of temperature T 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 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.

13 CO column density
We used the 13 CO(2-1) line to compute the 13 CO 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 N13 CO for each set (e.g., Eqs. 4, 5 in Kong et al. 2015): where A 2→1 is the Einstein coefficient for the J=2-1 transition, and B 0 is the rotation constant of 13 CO.The excitation temperature T 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 mb and optical depth τ (v) were solved using the Gaussian fits to the 13 CO(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) which is accurate to 10 percent for T 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 N13 CO ≥ 7 × 10 15 cm −2 , equivalent to N H 2 ≳ 5 × 10 21 cm −2 (assuming a H 2 -to-13 CO ratio of ∼ 7 × 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.

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

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

RADEX fitting
We aimed to constrain the H 2 volume density for a set by reproducing the 13 CO line intensities (Step 6).In RADEX, the background radiation temperature was set to 2.73 K, and only 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-ofsight 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 ρ-PDF is considered a reasonable representation of the underlying ρ-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).= (20 h 37 m 59 s , 39 • 55 ′ 01 ′′ ).The solid red and blue lines indicate the peak T mb that RADEX outputs for each value of n H 2 , while the dashed lines are at the corresponding peak T mb observed along the line of sight.The dotted green line, at n = 104.5 cm −3 , marks the best fit density for this set, where the χ 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) = (20 h 34 m 37 s , 40 • 30 ′ 21 ′′ ).The best fit density from RADEX fitting was n = 10 6 cm −3 , which was the highest input in our grid.
Because the peak radiation temperatures of 13 CO(1-0) and 13 CO(2-1) were known for each set, while n H 2 was unknown, we parameterized over H 2 volume density and created a grid of identical inputs with only the volume density changed.We tested a range for n 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 H 2 which produced the smallest χ 2 value between the observed and modeled peak radiation temperatures in both 13 CO lines (Step 7).The χ 2 is computed as where O stands for observation and E stands for model.σ is the spectral noise.Currently, we are limited by just two 13 CO line transitions.Adding more transitions in the future will give stronger constraints.While the vast majority (∼ 88%) of sets had a local minimum χ 2 within our test range, the remainder showed χ 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 χ 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 H 2 density distribution we measure reflects the true distribution throughout the denser parts of Cygnus-X South.
Our choice of cutoff at χ 2 < 10 was arbitrary; it represents the best ∼ 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 χ 2 cutoff values, spaced by powers of two from 2.5 to 1,280.Further increasing the χ 2 cutoff value did not increase the number of sets with N13 CO ≥ 7 × 10 15 cm −2 (required by criterion 5 in §3).The results are shown in Figure 8 and discussed in Section 4.3.

Global spectra
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 13 CO lines peaks at about one-fifth the main beam temperature of its counterpart in CO at the same transition.

Integrated intensity and first-moment maps
In Figure 5, we show maps of the integrated line intensity for CO(2-1) and 13 CO(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 13 CO 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 13 CO(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 negativevelocity region in the first moment maps.The map in 13 CO is missing many pixels because the emission intensity in those pixels was never more than 3σ above the noise level.This is a consequence of emission in 13 CO being weaker overall.

ρ-PDF
In Figure 7, we show the ρ-PDF of the best-fit volume densities for the 5,480 sets that reached χ 2 < 10.The plot, in logarithmic scales, covers a range in log(n H 2 ) from 3.6-6.0, in logarithmic increments of 0.1 dex; no well-fit sets had bestfit volume densities less than 10 3.6 cm −3 .A large number of sets χ 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 mb for each transition gradually approached the observed peak T 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 highdensity 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 6.0 3.6 P d log(n 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 H 2 )-space, and the bin heights as the ycoordinate in log(P )-space.This is consistent with the corresponding column-density log-log fit between A V ∝ N 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 ρ-PDF (P ∝ ρ α ).A linear-regression fit between log(n H 2 ) = 4.5 and 5.5 gives a powerlaw slope of α = −1.122± 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 Averaged spectra for all four molecular lines across the entire region of mutual overlap between cubes.Pixels with σ > 2.5 K, which appeared only around the edges of the data cubes, were discarded.
permitted local minimum χ 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 N13 CO ≥ 7 × 10 15 cm −2 .All of these plots showed some power-law slope between log(n H 2 ) = 4.5 and 5.5, with a value between α = −1.0 and −1.25 (Figure 8).Uncertainty in the slope is much greater for cutoffs below our main value of χ 2 = 10 but somewhat less for cutoffs above it, with the highest cutoff, including all sets of sufficient column density, having a slope of α = −1.155± 0.028.If there were a strong power-law in the true ρ-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 χ 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 α = −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.

T ex Consistency
To confirm the accuracy of our T 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 ex values for 13 CO(2-1) and 13 CO(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 χ 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 (∼ 87%) of RADEX's 13 CO 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 ex error threshold for including them in the first place.
Because we assumed LTE when calculating T 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 13 CO are understandable.RADEX systematically returns lower T 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 ex for the J=1-0 line because stimulated emission at high densities contributes to T ex above what LTE would predict from gas kinetic temperature alone.

DISCUSSION
The shallowest ρ-PDF slope in Kritsuk et al. (2011, hereafter K11) was −1, which was caused by rotationally supported cores1 at a very high density of order n H 2 ∼ 10 9 cm −3 .Such a density is not observed in our ρ-PDF in Figure 7 and   .Slopes of the linear regression powerlaw fits with varying cutoffs for the largest permitted minimum χ 2 .In all cases, the power-law fit was between log(n H 2 ) = 4.5 and 5.5, and only sets with column density N13 CO ≥ 7 × 10 15 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.
could not be identified with these CO lines.In addition, the SMT angular resolution of ∼32 ′′ corresponds to a spatial scale of ∼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 ± 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 ∼ 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 ∼ 37 mag.Our mapped region has more overlap with their Cygnus South region.Due to the limited dynamic range, our fitted ρ-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 ρ-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 ρ-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 ex agrees within 30% for two 12 CO 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 ρ-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 ρ-PDF generally agrees with S22.However, Murase et al. (2023) did not find power-law profiles in 13 CO N-PDF in the majority of their sample.Note that they broke the Cygnus-X complex into multiple in-dividual 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 ex ∼10-16K, the depletion of 13 CO could impact the column density derivation and thus the N-PDF, assuming the dust temperature is similar to T 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 ≳ 40 mag.They concluded that N-PDFs based on molecular lines were not well constrained.

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 13 CO using the ARO SMT telescope.
The maps extend from 20 h 24 m to 20 h 39 m in RA and from 38 • to 40 • 43 ′ in DEC, covering an area of ∼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), 13 CO(1-0), CO(2-1), and 13 CO(2-1).Through the analysis, we have derived physical properties of the molecular gas, especially the gas volume density.
We use astrodendro to identify velocitycoherent 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 ex for each CO line assuming optically thick conditions.We se-Baade lect sets whose T ex agree within 30% between the two CO lines so that LTE is a reasonable assumption for the sets.We then derive the 13 CO column density for each set.These physical properties are fed to the RADEX program with which we iteratively find the best H 2 volume density by matching the intensities between the model and the observation for the two 13 CO lines.We construct the volume density probability distribution function (ρ-PDF) which shows a single power-law tail between log(n H 2 )= 4.5 and 5.5 with a slope of ∼ −1.12 ± 0.05.Higher transition lines are needed to trace higher densities.
The ρ-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 syn-thetic observations with feedback may help clarify these causes.Our observational study of ρ-PDF provides another useful test for theories.

Figure 1 .
Figure1.Spitzer 8 µ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.

Figure 2 .
Figure 2. Plots illustrating Gaussian fit results for the pixel at (RA, DEC) = (20 h 35 m 26 s , 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.

Figure 3 .
Figure3.Plots illustrating the n H 2 fitting process.The left panel is for the set at the pixel at (RA, DEC) = (20 h 37 m 59 s , 39 • 55 ′ 01 ′′ ).The solid red and blue lines indicate the peak T mb that RADEX outputs for each value of n H 2 , while the dashed lines are at the corresponding peak T mb observed along the line of sight.The dotted green line, at n = 10 4.5 cm −3 , marks the best fit density for this set, where the χ 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) = (20 h 34 m 37 s , 40 • 30 ′ 21 ′′ ).The best fit density from RADEX fitting was n = 10 6 cm −3 , which was the highest input in our grid.
Figure4.Averaged spectra for all four molecular lines across the entire region of mutual overlap between cubes.Pixels with σ > 2.5 K, which appeared only around the edges of the data cubes, were discarded.

Figure 5 .
Figure 5. Integrated intensity maps for the CO(2-1) and 13 CO(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 .
Figure6.First-moment maps for the CO(2-1) and 13 CO(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σ above the noise level were masked out.The color bar has units of km s −1 .

Figure 7
Figure7.ρ-PDF derived from best-fit volume densities for each set.Only those with a best-fit χ 2 < 10 were included.The blue line displays a linear-regression power law fit to the data between log(n H 2 ) = 4.5 and 5.5.
Figure8.Slopes of the linear regression powerlaw fits with varying cutoffs for the largest permitted minimum χ 2 .In all cases, the power-law fit was between log(n H 2 ) = 4.5 and 5.5, and only sets with column density N13 CO ≥ 7 × 10 15 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 .
Figure9.Histogram of relative errors between T ex calculated for CO transition lines and those estimated for the corresponding 13 CO lines by RADEX.The horizontal scale goes up to 0.85 to include the highest outlier in the J=1-0 plot (not visible).
H 2 to best match the 13 CO(1-0) and 13 CO(2-1) line emission between RADEX and observations.8.Confirm the validity of LTE by checking theT ex between RADEX and observations.
13 CO line dispersion, carry out RADEX modeling (van der Tak et al. 2007) of the 13 CO line emission with an initial guess of n H 2 .7. Iterate n