QUANTIFYING OBSERVATIONAL PROJECTION EFFECTS USING MOLECULAR CLOUD SIMULATIONS

We begin with a broad overview of how cloud information can be distorted during the observation process. A spectral line observation of a molecular cloud can be thought of as a transformation from a set of intrinsic quantities—density, velocity, temperature, chemical abundance—to a map of intensity in PPV. Information about the original cloud structure is lost during several steps of this transformation. The first problem is the projection from PPP space to the PPV space of the observation. This step is described by the following equation:

$$\begin{equation} \rho _{\rm PPV}(x, y, v) = \sum _{v_z(x, y, z) = v} \rho _{\rm PPP} (x, y, z) \left| \frac{\partial z}{\partial v_z(x, y, z)} \right|, \end{equation} \tag{ 1 }$$

where ρ_PPV is the density of material in PPV space (g cm⁻² km⁻¹ s), ρ_PPP is the density in PPP (g cm⁻³), and the derivative is the standard Jacobian used when transforming densities between coordinate systems.

From the perspective of feature identification, two aspects of this transformation break the correspondence between PPP structures and PPV features. First, distinct positions along the same line of sight that move at similar velocities will project to the same region in PPV. Thus, a feature in PPV may sample two or more density structures. This is the problem of superposition, and is illustrated in Figure 1(a). Second, spatial variations in v_z affect the gradient term in Equation (1), and can modulate the ρ_PPV field independently of ρ_PPP. In other words, a single density structure can map to multiple velocity-induced PPV features. This is shown schematically in Figure 1(b).

Zoom In Zoom Out Reset image size

In addition to projection, observations are also subject to chemical and radiative transfer effects, which further distort the intensity field from the density field via spatially variable excitation, abundance, ionization, and opacity conditions (Bell et al. 2006; Lee et al. 2013; Pineda et al. 2008). The ρ_PPV field determines the column density and, along with the temperature, the collision rate of the gas. Both of these affect the intensity field—the column density sets the opacity and can obscure background features, while the collision rate affects the excitation state and emissivity of the gas. Finally, observations are subject to noise and spatial filtering, which further degrade the data, increasing opportunities for confusion.

The net effect of these phenomena is too complicated to characterize analytically. Instead, we turn to numerical simulations, where cloud properties can be compared before and after synthetic "observations." With simulations we also have the freedom to disable individual aspects of the observation process, to better isolate the influence of each factor.

Several authors have investigated the relationship between PPP structures and PPV structures, using different simulations and analysis techniques. An early study by Adler et al. (1992) investigated structures identified in longitude–velocity diagrams of a synthetic model of the Galaxy. They found that many of these identifications were superpositions of separate PPP regions. On a smaller scale, Ballesteros-Paredes & Mac Low (2002) simulated observations of molecular cloud clumps in local thermodynamic equilibrium. They measured a number of canonical relationships in both PPP and PPV, including the clump mass spectrum, size–linewidth relationship, and size-density relationship. They found that the amount of confusion due to superposition depends on both the strength and spatial scale of turbulent driving—if turbulence induces stronger or smaller-scale density perturbations, confusion worsens. Despite confusion problems, both of these papers recovered mass–size and size–linewidth scaling relationships similar to real clouds. Issa et al. (1990) also considered how cloud superposition affects the size–linewidth relationship in Galactic CO surveys, and demonstrated that the slope of the relationship is quite robust to crowding.

Gammie et al. (2003) analyzed the aspect ratios of molecular cloud clumps projected onto the sky. The observed distribution of cloud aspect ratios is affected both by projection into 2D, as well as superposition of distinct cloud features. Offner & Krumholz (2009) also studied the intrinsic and projected distribution of simulated core shapes, noting that the intrinsic triaxial shape of cores is lost during projection. Using a similar analysis, Jones & Basu (2002) attempted to invert the distribution of apparent axis ratios to recover the intrinsic shape distribution of clouds. This inversion assumed that observed cloud shapes are unaffected by superposition, however.

Pichardo et al. (2000) correlated PPV structures in MHD simulations with the original density field and velocity field. The features in their PPV maps resemble the patterns in the velocity field more than the density field. There were also more small-scale PPV structures than there were PPP structures. The authors attributed these small PPV structures to the velocity-induced structures shown in Figure 1(b).

Ostriker et al. (2001) cataloged "observed" structures in their 3D MHD simulations by identifying regions of contrast in 2D projections. They, too, noted that features identified in this way often consist of several superposed density structures. Even though their feature extraction process ignored line-of-sight velocity information, they reported that velocity information is often unable to disambiguate superpositions in 2D.

In a series of papers (Lazarian & Pogosyan 2000, 2004, 2006), Lazarian et al. developed a mathematical formalism to describe how to recover statistical properties of turbulence (namely, the turbulent velocity and density power spectrum) from observations. This approach differs from the previous references in that it does not focus on the reality of observed structures. Instead, it considers the structure functions of spectra and slices or slabs of PPV cubes. These papers derive the expected shape of observable spatial and spectral structure functions, for idealized turbulence. The advantage of this analysis is that it explicitly treats PPV superposition, though other factors like spatially varying excitation and abundance conditions are not treated.

An approach based on Principal Components Analysis has similarly been used to measure cloud statistics without identifying specific structures with clear boundaries (Heyer & Schloerb 1997; Brunt & Heyer 2002, 2013). This method decomposes PPV datacubes into linear superpositions of "eigenimages" with different spatial and spectral extents. These extents, derived from the spatial and spectral autocorrelation of the decomposed data, are used to reconstruct scaling relationships like the velocity power spectrum.

Recently, Shetty et al. (2010) carried out an analysis similar to the work by Ballesteros-Paredes & Mac Low (2002), and measured how projection affects the measurement of size–linewidth relationships in molecular cloud substructures. They identified structures using the dendrogram algorithm, which is explicitly designed to characterize hierarchical structures like molecular clouds. They concluded that superposition can change the power law scaling coefficient for the mass–size and size–virial relationships by ∼ ± 0.5, and the linewidth-size relationship by ∼ ± 0.05 (see Table 1 of that paper).

In summary, a broad range of work using a variety of numerical simulations has found that projection effects impact the study of cloud structures. Even though we can't explicitly measure this effect in real clouds, this work suggests it is important to quantify projection effects in the context of typical observed quantities. The analysis presented in this paper builds upon these previous studies in a few key aspects. First, we develop a method to systematically cross-match observed and real cloud structures, for all structures in a cloud simulation. This provides a more detailed view into how structure analysis is affected by factors like superposition. We also carry out a more detailed, non-LTE radiative transfer to better model the important effects of excitation and opacity.

We have applied the similarity analysis described above to two cloud simulations. Each of these simulations is meant to broadly represent the conditions in a molecular cloud like Perseus (Ridge et al. 2006; Bally et al. 2008). However, the mean simulation temperature, density, and line-of-sight dimension may differ from the true values in Perseus by factors of two.

The first simulation, henceforth O1, is performed with the orion adaptive mesh refinement (AMR) code (Truelove et al. 1998; Klein 1999). The simulation assumes a simple isothermal equation of state, which means that it is scale-free for density and temperature (e.g., Offner et al. 2008). The simulation is produced following the same procedure in Offner et al. (2013), which we briefly summarize below.

The simulation domain begins with a uniform density, which we perturb with a random velocity field for two crossing times. The input field has a flat power spectrum for large wavenumbers, 1 < k < 2, and we normalize the perturbations to maintain a constant 3D Mach number, $\mathcal {M}=22$. This Mach number was chosen to reproduce the observed velocity dispersion in Perseus. After the gas achieves a well-mixed turbulent state, we turn on self-gravity and allow collapse to proceed. The simulation has a 256³ base grid, four levels of AMR refinement, and employs periodic boundary conditions. New grids are added automatically to satisfy the Jeans criterion for a Jeans number (the ratio of cell size to the local jeans length) of 0.25 (Truelove et al. 1997). When the Jeans criterion is violated on level four within a collapsing region, a sink particle is introduced (Krumholz et al. 2004). Our similarity analysis is performed at half a freefall time, when ∼700 M_☉ is contained in sink particles (2.3% of the gas).

The second simulation (hereafter S11) is an updated version of a model originally presented in Shetty et al. (2011b) (specifically, the n100 simulation in that paper). It was generated using a modified version of the zeus-mp MHD code (Stone & Norman 1992a, 1992b; Norman 2000; Hayes et al. 2006). S11 differs from O1 in that it ignores gravity, includes a 5.85 μG magnetic field and includes treatments of the non-equilibrium heating and cooling of the gas, the penetration of UV radiation into the cloud, and also a simplified treatment of the formation and destruction of H₂ and CO. The original Shetty et al. (2011b) simulation used the chemical model presented in Glover et al. (2010), but the updated version presented here uses instead a treatment based on Nelson & Langer (1999), as described in Glover & Clark (2012). However, as explored in some detail in Glover & Clark (2012), this change in chemical networks does not significantly affect the CO distribution in the gas. Our updated version of the Shetty et al. (2011b) simulation also includes a number of improvements in the way in which the thermal evolution of the gas is modeled, as described in Appendix A of Glover & Clark (2012).

The S11 simulation begins with uniform density which is perturbed with a random turbulent velocity field for three turbulent crossing times. The input field is similar to that in the O1 simulation, and is normalized to maintain a constant 3D rms velocity dispersion of 5 km s⁻¹. Converting this value to a Mach number is complicated by the fact that the gas in the S11 simulation is not isothermal and hence has a spatially varying sound speed. The volume-weighted mean Mach number is relatively low, ${\cal M} \simeq 6$, because much of the cloud volume is filled by warm, CO-poor gas with T ∼ 60–70 K. If, however, we compute ${\cal M}$ only for gas with more than 10% of its carbon in CO, we find a much higher value, ${\cal M} \simeq 14$, as this gas is much colder, with T ∼ 10–20 K.

Table 2 summarizes the properties of each simulation.

We used the radiative transfer program RADMC-3D (Dullemond 2012) to generate synthetic observations of each simulation in ¹²CO (J = 1–0), ¹²CO (J = 3–2), and ¹³CO (J = 1–0), using the large-velocity-gradient (LVG) approximation (Sobolev 1957; Shetty et al. 2011a). The observations were gridded to a spatial resolution of 0.1 pc pixel⁻¹, and velocity resolution of 0.05 km s⁻¹. Finally, we added noise to each cube (0.6K, 0.15 K and 0.25 K for the ¹²CO (J = 1–0), ¹²CO (J = 3–2), and ¹³CO (J = 1–0)  transitions, respectively). These are representative of the noise values of present-day cloud surveys in these transitions (Ridge et al. 2006; Davis et al. 2010).

We used the dendrogram algorithm to catalog intensity structures in each synthetic observation, as well as the density structures in the PPP density fields. The dendrogram algorithm is described in detail in Rosolowsky et al. (2008). Briefly, each structure in a dendrogram corresponds to a surface of constant intensity (in the observation) or density (in the PPP cube). The name dendrogram refers to the fact that these surfaces are hierarchically nested inside each other and representable via tree diagrams (Figure 3). Dendrograms capture the hierarchical structure of molecular clouds—a clear advantage over non-hierarchical clump-finding algorithms—and dendrogram decompositions are not sensitively dependent on tuning parameters of the algorithm (Pineda et al. 2009).

Zoom In Zoom Out Reset image size

When constructing a dendrogram, the main freedom one has is the degree to which structures are further decomposed into nested substructures. This process is called "pruning" the dendrogram, since it amounts to controlling how many "branches" (structures) are in the decomposition. The main purpose of pruning is to suppress the extraction of insignificant structures that are poorly resolved or are possible noise fluctuations. A dendrogram constructed with no pruning assigns every local intensity or density maximum (including every noise spike) to a unique structure. The effect that pruning has on the statistical properties of a dendrogram has been studied in detail by Burkhart et al. (2013).

Each dendrogram in this work is pruned such that every leaf contains a local intensity maximum that is brighter than the neighboring 7 voxels in any direction. Each leaf (the brightest, most-compact structure in a hierarchy) also contains at least 800 voxels (for spectral line cubes in PPV) or 100 voxels (for density cubes in PPP) and contains a voxel that is at least 7σ brighter than the contour at which the leaf merges with its neighboring structure. The PPP dendrogram is pruned less heavily than the PPV dendrograms, yielding a catalog with more structures. This prevents PPV structures from being poorly matched to PPP structures simply because the density structure decomposition is too coarsely grained. One of the convenient aspects of dendrograms is that the boundaries of non-pruned structures are independent of the pruning; in other words, while pruning can add or remove structures from a catalog, it does not affect how the included structures are defined. Thus our choice of pruning has little effect on subsequent analysis, other than to exclude from consideration the smallest cloud substructures. We explore how sensitive our analysis is to our pruning choice in Section 3.4.

We compare synthetic CO observations of the simulations to the COMPLETE Perseus data (Ridge et al. 2006) using several diagnostics: namely, the distribution of column density, velocity dispersion, and line intensity in various CO transitions. Both the O1 and S11 simulations represent the general physical properties of Perseus. Appendix A provides details about how each quantity was extracted from the data. Neither simulation agrees with Perseus when these diagnostics are examined in detail, but they are the closest available approximations. We discuss the limitations imposed by the suitability of the simulations in Section 3.10.

Figure 4 shows isosurface renderings for Perseus, O1, and S11, in the ¹²CO (J = 1–0)  transition. Isosurfaces are drawn at 3, 8, and 15 K. The O1 simulation (panel b) stands out from the other two panels in this figure. Compared to Perseus and S11, it has more space-filling emission at 3K and a lack of emission at 15 K.

Zoom In Zoom Out Reset image size

To make the differences between the simulations and Perseus more precise, Figure 5 shows three statistical comparisons between the S11 simulation and Perseus: the distribution of column density, ¹²CO (J = 1–0)integrated intensity, and line-of-sight velocity dispersion. The simulation has a higher average column density than Perseus and fainter CO lines. To rough approximation, integrated line intensity increases with gas density, temperature, and velocity dispersion. Since the S11 simulation is at a higher column density than Perseus (panel a), the stronger lines in Perseus are probably due to hotter gas, higher turbulence, and/or poor modeling of CO abundance.

Zoom In Zoom Out Reset image size

The velocity dispersion of the spatially-averaged spectrum is shown as a vertical line in Figure 5(c). This number is larger than the typical line-of-sight velocity dispersion (the histograms in panel c), due to spatial velocity gradients across the region. In other words, the linewidth of the spatially-averaged spectrum is different from—and larger than—the mean of the line-of-sight linewidth distribution. While the spatially-averaged velocity dispersion in S11 is comparable to Perseus, the individual line-of-sight velocity dispersions in Perseus are skewed towards higher values. We speculate that this is related to the characteristic depth of each line of sight; lines of sight in S11 often intersect a single, ∼1 pc-thick filament of material, which moves more coherently—i.e., with a smaller velocity dispersion—than the cloud as a whole. It may be that typical lines of sight in Perseus pass through material spread out over a longer column, and hence have larger dispersions. This may explain the higher integrated intensities in Perseus, since W = ∫Tdv.

The O1 simulation exhibits similar discrepancies (Figure 6). Its mean density and Mach number were chosen to match Figures 6(a) and (c). It better reproduces the column density distribution in Perseus by construction, but like S11, the integrated emission is too faint. The mode of the line-of-sight velocity dispersion matches Perseus moderately well (panel 6c), but Perseus has a longer tail of high-velocity dispersion material. The velocity dispersion of the cloud-averaged spectrum (dashed lines in 6c) is 50% larger in Perseus.

Zoom In Zoom Out Reset image size

Turbulence in both simulations is arbitrarily driven, and this non-physical prescription does not reproduce statistics of the 1D line-of-sight velocity field in Perseus. One possibility for this discrepancy is that the simulations are not fully resolving the velocities on the smallest scales. Resolution studies of grid based codes indicate that on small scales the amplitude of the velocity power spectrum is larger in simulations with higher resolution (e.g., Vestuto et al. 2003). For the intermediate resolutions considered here, therefore, the velocity dispersions on the smallest scales may be underestimated. Additionally, some of the high-dispersion emission from Perseus coincides with the parsec-scale, stellar-driven shells studied in Arce et al. 2011 (in particular, Arce's CPS 4 and CPS 5). Neither simulation considers the effect of stellar feedback.

Lee et al. (2013) have also recently compared a spatially-truncated version of the S11 simulation to Perseus. They too note that the S11 simulation is under-luminous. They report a more extreme under-luminosity in W(¹²CO (J = 1–0)) of eight times relative to Perseus, though their truncation acts to further diminish the line intensity. They posit that velocity crowding is responsible for this discrepancy–the lower velocity dispersion in S11 makes superposition more likely, and optically thick ¹²CO features are more likely to obscure each other. Indeed, we find that the peak of the W(¹³CO (J = 1–0)) distributions in both O1 and S11 better reproduce the values in Perseus (Appendix A). The optical depth in W(¹³CO (J = 1–0)) is lower and should be less affected by velocity crowding. Nevertheless, Perseus still shows an excess of large ¹³CO (J = 1–0)  intensities that neither simulation reproduces.

We hope that these discrepancies will motivate further efforts to generate simulated clouds which better agree with the statistical properties of Perseus. To that end, we discuss a suite of diagnostics in Appendix A, to facilitate standardized comparisons in the future. For the purposes of this paper, it is important to bear in mind that the "observed" properties in both simulations are under-luminous and under-dispersed compared to Perseus.

As a first comparison between the properties of the O1 and S11 simulations, Figure 9 compares the ¹²CO (J = 1–0)  transitions for each simulation, and Figure 10 shows a color-coded PP slice of the S11 simulation. The O1 and S11 simulations show marked differences. The S11 simulation has overall better match quality. The S11 simulation also has a higher dynamic range of structure brightnesses; this is probably due to the fact that that simulation explicitly treated gas heating and CO dissociation, whereas the O1 simulation is isothermal and assumes a constant CO abundance. Heating and dissociation give the S11 simulation more freedom to affect line intensity (by raising the excitation temperature as well as the abundance of emitting material). CO is dissociated in low column density regions of S11, and the simulation has large regions devoid of emission.

Noise has been added to each synthetic observation, to match the noise levels in present-day cloud observations in these transitions (Section 2.1). This raises the following question: to what extent does noise make it more difficult to extract cloud features, and hence lower the match quality? To address this, we also show in Figure 8 the same PP slices without noise. Note that the presence or absence of noise has little bearing on the match quality of most structures. Instead, the lower match quality in the ¹²CO (J = 1–0)  transition seems dominated by the high filling factor and opacity of emission, which crowds features in PPV and leads to superposition.

The ability to recover structures varies with tracer; Figure 7 shows that structures in ¹²CO (J = 1–0)  are most affected by projection problems, followed by ¹²CO (J = 3–2)  and ¹³CO (J = 1–0). The latter two transitions trace higher densities and are optically thinner than ¹²CO (J = 1–0). It is unclear from Figure 7 alone whether the poor match quality in ¹²CO (J = 1–0)  is the result of the higher filling factor or higher opacity in that line.

We can partially decouple the effects of filling factor and opacity by disabling absorption in the radiative transfer. This doesn't prevent crowding or superposition in PPV, but it does prevent structures from blocking radiation. If opacity is the primary problem with ¹²CO (J = 1–0)  observations, then disabling absorption should increase the overall match quality.

To test the effects of filling factor alone, we perform a modified radiative transfer calculation on the O1 simulation, where opacity is disabled. The full equation of radiative transfer is given by

$$\begin{equation} I_\nu = \int e^{-\tau (z)} B_\nu \left(T_{\rm ex}\left(z\right)\right) \alpha (z)\, dz, \end{equation} \tag{ 5 }$$

where τ(z) is the optical depth to depth z, B_ν is the Planck function, T_ex(z) is the excitation temperature of the gas, and α the absorption coefficient. T_ex and α are functions of the gas density at each energy level.

For our modified radiative transfer calculation, we run RADMC-3D as normal to compute the level populations and hence T_ex and α throughout the simulation volume. However, we then integrate a modified equation of radiative transfer with no absorption:

$$\begin{equation} \tilde{I}_\nu = \int B_\nu \left(T_{\rm ex}\left(z\right)\right) \alpha (z)\, dz. \end{equation} \tag{ 6 }$$

This produces a modified version of the O1 synthetic observations, which we label as O2. The only difference between O1 and O2 is that O2 includes no absorption. Figure 11 compares the match quality for the ¹²CO (J = 1–0)  transition in simulations O1 and O2. The O2 structures are markedly higher quality. Figure 11 indicts the e^−τ(z) term in Equation (5) as a primary reason for low match qualities in the ¹²CO (J = 1–0)  transition.

Note that this experiment does not fully disable the effects of opacity. In addition to the e^−τ(z) term, opacity acts to increase the excitation temperature of the gas by absorbing radiation emitted by other parts of the cloud. Disabling this absorption could lead to de-excitation of parts of the cloud, and this spatially-varying excitation would partially decouple the intensity field from the density field and decrease the filling factor of emission.

In Section 2.2, we described our pruning strategy—namely, we require that each structure contains a voxel 7σ above the ambient intensity and contains at least 800 voxels altogether. Figure 12 shows a less aggressive pruning strategy, where we relax the N > 800 voxel criterion to N > 400. There are more structures in this dendrogram (283 structures in the ¹³CO (J = 1–0)  transition, compared to 191 in the original pruning). Compared to Figure 7(c), these additional points are concentrated at Areas <0.5 pc². We reiterate that the points in Figure 12 are a superset of Figure 7(c), and the extra points are substructures nested inside the structures from the original pruning.

On average, the new structures have modestly lower match qualities: 38% of the new structures have q < 0.5, compared to 27% of structures in the original pruning. However, we interpret Figure 12 as evidence that the reality or quality of dendrogram-identified structures is fairly insensitive to the details of pruning, provided that statistically-insignificant noise-spikes are not identified as structures.

What impact does confusion caused by projection and radiative transfer have on subsequent analyses? This is a problem-dependent question, and we focus here on the fairly common virial analysis. The virial parameter, often defined as $\alpha = 5 \sigma _v^2 R / {G M}$, gives the approximate ratio of kinetic to gravitational potential energy (McKee & Zweibel 1992). The value α ∼ 2 denotes the approximate equipartition between these two energy terms and is often used to assess the boundedness of a given structure. However, the true virial state of an object is affected by several additional unobservable terms (for example, surface terms and magnetic energy; Ballesteros-Paredes 2006; Dib et al. 2007; Bertoldi & McKee 1992), and we emphasize that the α < 2 threshold is a crude proxy for boundedness. Furthermore, the virial analysis implicitly assumes that structures are roughly spherically-symmetric, which does not well-describe the larger features in a dendrogram of a filamentary cloud.

Figure 13 shows, for the ¹³CO (J = 1–0)  transition of simulation O1, the relationship between size, velocity dispersion, mass, and virial parameter for each structure in the dendrogram. Each point is color-coded by match quality as before. Appendix B describes how each quantity is measured.

We also show a simple power law fit to the size–linewidth, mass–size, and virial–size relationship, which are frequently measured in molecular cloud studies. The black line shows the fit to all points, while the blue line shows structures with q > 0.5. The slopes of these lines essentially reproduce Larson's classical relationships (M ∼ R², V_rms ∼ R^0.5, α ∼ R⁰; Larson 1981). Ignoring the low-quality match structures affects the slope by ≲ 0.05.

Assessing the uncertainty in these scaling relationships is subtle. A naive least-squares error analysis suggests a small uncertainty for the scaling exponents (∼.025). However, these data points correspond to nested structures and are not independent of each other. Consequently, the least-squares error estimate is overly optimistic. We have experimented with different sensible strategies for pruning the dendrogram, as well as different definitions for the size of an irregular structure (see Appendix B). Varying these options can change the slope of the scaling relationships by ∼ ± 0.2 (see, for example, Table 3), and we feel this is a more appropriate estimate for how precisely the scaling relationships are constrained.

For the O1 simulation, then, filtering based on q does not significantly affect the scaling relationships one obtains from PPV data.

PPV-derived properties are most often used as approximations for properties of the (partially un-measurable) 6D spatial-kinematic state of the cloud. How accurate are these approximations?

Figure 14 compares, for the ¹³CO (J = 1–0)  transition of O1, quantities measured in PPV with the equivalent measurement of each PPV structure's nearest PPP match (measured in PPP). The point sizes indicate structure size, and the dashed lines are a factor of two above and below the 1:1 line. In the lower right corner of each panel, we also report a few summary statistics: the geometric mean of the ratio of PPV/PPP measurements and the geometric standard deviation of this ratio. The geometric mean is defined as μ_g = (∏x_i)^1/N and the geometric standard deviation as $\sigma _g = \exp {(\sqrt{\vphantom{A^A}\smash{{({{1}/{N}) \sum {\ln (x_i / \mu _g)}}}}})}$. The base-10 log of σ_g is the scatter about the ratio μ_g in dex. These numbers measure the fractional bias from and scatter about the 1:1 line, respectively. We report the geometric mean and standard deviation for all points, as well as the subset of points with q > 0.5. In this analysis, we exclude structures which touch the edge of the cube, since their full extent is not measured.

There are several features of Figure 14 worth commenting on. First, the strongest outliers are the reddest, lowest-quality PPV structures, shown as the red points in the upper-left corner of panels (a)–(c). These structures have no correspondence to any PPP structure. Because these artificial features cannot be sensibly matched to anything in the PPP cube, they are arbitrarily matched to the largest density structures and occupy the upper left corners of panels (a)–(c).

Second, structures with q > 0.5 are dispersed about the 1:1 line by σ_g = 1.4–1.7 (∼0.1–0.2 dex) in panels a–c. The act of filtering on q reduces scatter by 15–30% for mass and velocity. The reduction is larger for masses in panel a, but the dispersion in these points is dominated by the handful of outliers in the upper left corner.

Finally, the virial parameter plot (panel d) shows higher scatter for q > 0.5 structures—0.34 dex, or a factor of 2.2—than do the mass, size, or linewidth plots. The individual errors in the mass, size, and velocity dispersion measurements compound when measuring the virial parameter, and this property is the least-faithfully recovered.

To summarize Figure 14, radiative transfer and projection effects produce a factor of 1.4–2 uncertainty on kinematic properties derived from the O1 simulation using ¹³CO (J = 1–0)  emission.

The O1 simulation includes gravity. Gravity acts to gather and collapse gas, which may create locally crowded regions of high confusion. On the other hand, gravitational collapse should also gather diffuse material on large scales. This may decrease the filling factor on large scales, in mitigate confusion. To probe how important each of these effects are, Figure 15 shows the equivalent set of comparisons of O1 at an earlier epoch, at the instant gravity is enabled. This timestamp captures the steady-state turbulent structure of the O1 simulation, before gravity has had an influence. The column density of the simulation at the original and earlier timestamp is shown in Figure 16.

Gravity causes structure collapse at small scales, producing more dense, small regions (Figure 16(a)). Before gravity is enabled, structures on average are larger, more diffuse, and overlap more. This effects the kinematic properties in Figure 15 in the following ways: without gravity, there are fewer structures overall (147 versus 191 at the original simulation time). Because structures overlap more without gravity, there is a greater fraction of q < 0.5 structures (37% versus 28%). Finally, the virial parameter averaged over all structures is 5.8 without gravity, compared to 3.0 for the original O1 simulation.

Despite the moderately worse confusion, the scatter and bias in Figures 14 and 15 are remarkably similar. In other words, the influence of gravity in the O1 simulation does not greatly impact the ability to recover physical properties.

A main difference between the O1 and S11 simulations is the inclusion of limited CO chemistry in S11. Spatial abundance variations in CO can decouple the CO density from the H₂ density. This, in turn, can break the correspondence between CO intensity and PPP density.

Figure 17 shows the same comparisons for the S11 simulation. Panels b and c have a comparable amount of bias and scatter as the previous figure for O1, albeit with fewer structures overall. The most dramatic difference between this plot and Figure 14 is the mass comparison in panel a. The points in the S11 simulation are shallower than the 1:1 line—while PPV-structures recover sensible masses below ∼100 M_☉, M_PPV overestimates M_PPP for larger structures. This directly affects the virial plot in panel d, which exhibits a bias towards PPV-underestimates of α.

One may also wonder if the approximation to estimate mass from intensity (given by Equation (B6)) contributes to the mass discrepancy. This is unlikely. Masses are estimated by measuring the mass-to-light ratio of the brightest 5% of the pixels and using this as a conversion factor. This conversion factor underestimates the mass-to-light ratio for low-column density lines of sight, where gas is sub-thermally excited and emits inefficiently. Adopting a more appropriate mass-to-light ratio for these lines of sight would actually lead to even larger M_PPV measurements, which exacerbates the discrepancy.

The mass discrepancy for large structures is caused by CO dissociation at low column densities, which S11 includes but O1 does not. Figure 18(b) shows that the S11 simulation contains large regions devoid of CO. There is H₂ in these regions (panel a), but little CO due to dissociation. In other words, the topologies of CO and H₂ gas partially diverge on large scales. It is thus more difficult to find exact PPP-equivalents of large-scale PPV structures in S11. This leads to discrepancies which, evidently, are most pronounced for mass measurements. This doesn't explain why the bias is towards PPV mass overestimates, as opposed to underestimates. We do not have a simple explanation for the direction of the bias. For whatever reason, PPV structures are better matched (as quantified by Equation (3)) to slightly too-small PPP structures than they are to slightly too-large PPP structures. We speculate this has to do with the detailed topology of gas in PPP versus PPV.

We can verify that chemistry in S11 causes the mass discrepancy by repeating the analysis on a version of S11 without chemistry. To do this, we re-compute the radiative transfer on S11, assuming a constant temperature of 15 K, and a constant CO/H₂ abundance of 10⁻⁴. The resulting integrated intensity map is shown in Figure 18(c), and the kinematic comparisons in Figure 19. There are more structures in this simulation due to extra CO emission. Furthermore, the mass-mass plot follows the 1:1 line much better.

At the largest scales, the S11 simulation (which does not include gravity) has a virial parameter of α_PPV ∼ 9. This is higher than the O1 simulation, for which the largest-scale structures have a virial parameter of α_PPV ∼ 1. The S11 structures with the smallest values of α_PPV ∼ 1 are in fact cause for concern, since they indicate that gravitational and kinetic energies are comparable in magnitude. Because S11 does not include the effects of gravity, the dynamical nature of these structures is less faithfully modeled.

The virial parameter is most often used to estimate the gravitational boundedness of a structure, with values of α < 2 interpreted to indicate that a structure is bound. This interpretation is problematic, as it ignores other forces and oversimplifies the role of turbulence as a support against collapse (Bertoldi & McKee 1992; Ballesteros-Paredes 2006). Our kinematic analysis of S11 and O1 add another cause for concern: measurements of the virial parameter based on PPP and PPV data differ by a factor of two to three, and cluster around α_PPV = 1–5. This implies that, using CO data alone, it is unclear on which side of the α_PPP = 2 boundary many structures fall.

To illustrate this, we repeat an analysis similar to Goodman et al. (2009), who measured the fraction of low-virial parameter structures for the L1448 subregion of Perseus as a function of size.

For each simulation, we assign a virial parameter to each voxel according to the smallest structure to which that voxel belongs. Next, we bin the structures by size, and in each bin, make a mask of all the pixels associated with these structures; we refer to this set as $\mathcal {S}$. Finally, we compute the fraction of emission in these pixels with α_PPV < 2:

$$\begin{equation} f = \frac{\sum {\lbrace L(\boldsymbol {r}) | \boldsymbol {r} \in \mathcal {S}, \alpha (\boldsymbol {r}) < 2}\rbrace }{\sum {\lbrace L(\boldsymbol {r}) | \boldsymbol {r} \in \mathcal {S}}\rbrace } \end{equation} \tag{ 7 }$$

We plot this fraction as a function of size for O1, S11, Perseus, and L1448 in Figure 20. Note that the L1448 plot is slightly different from Figure 4 in Goodman et al. (2009) because we use a different scheme for measuring the fraction of α_PPV < 2 emission.⁸ Also, remember that α_PPV underestimates α_PPP in the S11 simulation, pushing the line higher than it would otherwise be.

Zoom In Zoom Out Reset image size

For the O1 and S11 simulations in Figures 20(a) and (b), we plot the relationship for all structures (black), as well as those with high match qualities (q > 0.5, blue). Because estimates of α in PPV are scattered by a factor of ∼2 about the corresponding PPP measurements, the grey bands show the range of possible values the black line can take if each value of α is mis-estimated by a factor of two. Because so many structures fall within a factor of two of α_PPV = 2, the grey band covers a large swath of the plot. Thus, in addition to the conceptual problems associated with inferring boundedness from the value of α_PPV, there is an intrinsic observational ambiguity associated with reliably determining structures to be above or below α = 2.

As discussed in Section 2.3, the O1 and S11 simulations do not reproduce several statistical properties of Perseus. In particular, both simulations tend to have lower line-of-sight velocity dispersions than Perseus. This may act to increase superposition effects in the simulations since the material is more crowded in PPV space. Similarly, the lack of external radiation fields in O1 (with no chemistry) suppresses dissociation of low-density CO, producing an artificially high filling factor of emission and greater likelihood for superposition. Clouds with significant excitation variation due to external heating may have less PPV superposition of features.

As the field of molecular cloud simulation and synthetic observation advances, we should be able to make new quantitative estimates of the degree to which various tracers are confused under different physical and observing conditions. The suite of diagnostics presented in Appendix A enables a standardized comparison between a simulated and observed dataset. These comparisons address the main statistical cloud properties most relevant for superposition analysis—column density, intensity, and linewidth—and can help to assess how well a given simulation acts as a surrogate for studying unobservable projection effects in a real dataset.

The metrics that follow measure properties along lines-of-sight. As such, we try to focus only on lines-of sight with substantial cloud emission and mask out regions with little cloud material. We base our masking on the process described in Pineda et al. (2008) and require each line of sight to satisfy the following inequalities:

• 1.  
  The peak ¹²CO (J = 1–0)  line intensity is at least 10σ above the T = 0.
• 2.  
  The peak ¹³CO (J = 1–0)  line intensity is at least 5σ above T = 0.
• 3.  
  The velocity dispersion of ¹²CO (J = 1–0)  is at least 0.8 × the velocity dispersion of ¹³CO (J = 1–0).

These cuts are applied both to the Perseus data and to each simulation. The first two cuts are self explanatory. Pineda et al. proposed the last cut to further filter noisy or pathological lines of sight; the rationale is that, since ¹²CO is more abundant and opaque than ¹³CO, it should always have a larger spatial and kinematic extent.

The distribution of column densities is reasonably well-characterized for nearby clouds: near-infrared extinction measurements trace column densities across the range $10^{20} \lesssim N_{\rm H_2} [{\rm cm}^{-2}] \lesssim 10^{23}$, and far infrared dust emission probes higher column densities (Goodman et al. 2009; Kelly et al. 2012; Lombardi 2009; Kainulainen et al. 2009).

Most cloud column density distributions are approximately log-normal, with mean column densities of $N_{\rm H_2} \sim 10^{21}$ cm⁻² and width parameters σ ∼ 0.3–0.5. In addition, some clouds (especially those currently undergoing star formation) display excess power-law tails at column densities above ∼3 × 10²¹ cm⁻² (Goodman et al. 2009; Kainulainen et al. 2009).

The column density distribution is shown in Figure 21(a) for O1, and 22(a) for S11.

The distribution of line-of-sight integrated intensity W = ∫I(x, y, v) dv is also straightforward to compute from spectral line observations. Its distribution is shown in Figures 21(b) and 22(b) for ¹²CO (J = 1–0), and Figures 21(c) and 22(c) for ¹³CO (J = 1–0).

The distribution of peak line-of-sight intensity is a crude measure of the excitation state of the gas; it breaks the degeneracy between excitation state and column density in the integrated line intensity. Its distribution is shown in Figures 21(d) and 22(d) for ¹²CO (J = 1–0), and Figures 21(e) and 22(e) for ¹³CO (J = 1–0).

Likewise, we can compute the distribution of line-of-sight velocity dispersions. We compute the velocity dispersion by computing the intensity-weighted second moment of velocity along each line of sight. Since the second moment is sensitive to faint emission at large velocity offsets, we only consider pixels 3σ above the background. Its distribution is shown in Figures 21(f) and 22(f) for ¹²CO (J = 1–0), and Figures 21(g) and 22(g) for ¹³CO (J = 1–0).

The previous diagnostics are 1D distributions, and say nothing about the correlation among different quantities. The joint distribution of line intensity and column density is particularly interesting, since the ratio of these quantities defines the much-studied X-factor. Higher line intensities at a given column density indicate higher excitation levels, lower opacity, greater abundance of the exciting molecule, and/or greater linewidth (if the line is opaque). The joint distributions are shown in Figures 21(h) and 22(h) for ¹²CO (J = 1–0), and Figures 21(i) and 22(i) for ¹³CO (J = 1–0).

Each of the above diagnostics can be converted into a numerical score, to quickly summarize how well a given simulation reproduces a given diagnostic. We base our score on the Kuiper statistic for two cumulative distribution functions:

$$\begin{equation} K = \max {\rm \left(CDF_A - CDF_B \right) + \max \left(CDF_B - CDF_A \right). } \end{equation} \tag{ A1 }$$

The Kuiper statistic is a modification of the well-known Kolmogorov–Smirnov statistic, and is more sensitive to discrepancies in the tails of distributions (see the discussion in Section 14.3.4 of Press et al. 2007).

For every comparison of 1D distributions, we define the score as 1 − K, where smaller scores indicate less similarity between the simulation and Perseus.

There are a few ways to generalize the Kuiper statistic to the 2D joint distribution of column density and line intensity (see Section 14.8 of Press et al. 2007). For each of these 2D distributions, we compute 4 cumulative distribution functions

$$\begin{eqnarray} {\rm CDF}_1(X, Y) &=& P(x < X, y < Y), \end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray} {\rm CDF}_2(X, Y) &=& P(x > X, y < Y), \end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray} {\rm CDF}_3(X, Y) &= &P(x < X, y > Y), \end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray} {\rm CDF}_4(X, Y) &=& P(x > X, y > Y), \end{eqnarray} \tag{ A5 }$$

compute the K statistic for each CDF and save the largest statistic. Our final score is defined as 1 − K_max/2. The factor of 2 correction is included because, in 2D, the Kuiper statistic varies between 0–2 whereas, in 1D, it varies between 0–1.

Table 4 summarizes these scores for the O1 and S11 simulation, using Perseus as the benchmark. We encourage other researchers to generate molecular cloud simulations that better reproduce these observational diagnostics.