Hierarchical Structure of Magnetohydrodynamic Turbulence In PositionPositionVelocity Space
Abstract
Magnetohydrodynamic turbulence is able to create hierarchical structures in the interstellar medium that are correlated on a wide range of scales via the energy cascade. We use hierarchical tree diagrams known as dendrograms to characterize structures in synthetic PositionPositionVelocity (PPV) emission cubes of isothermal magnetohydrodynamic turbulence. We show that the structures and degree of hierarchy observed in PPV space are related to the presence of selfgravity and the global sonic and Alfvénic Mach numbers. Simulations with higher Alfvénic Mach number, selfgravity and supersonic flows display enhanced hierarchical structure. We observe a strong dependency on the sonic and Alfvénic Mach numbers and selfgravity when we apply the statistical moments (i.e. mean, variance, skewness, kurtosis) to the leaf and node distribution of the dendrogram. Simulations with selfgravity, larger magnetic field and higher sonic Mach number have dendrogram distributions with higher statistical moments. Application of the dendrogram to 3D density cubes, also known as PositionPositionPosition cubes (PPP), reveals that the dominant emission contours in PPP and PPV are related for supersonic gas but not for subsonic. We also explore the effects of smoothing, thermal broadening and velocity resolution on the dendrograms in order to make our study more applicable to observational data. These results all point to hierarchical tree diagrams as being a promising additional tool for studying ISM turbulence and star forming regions for obtaining information on the degree of selfgravity, the Mach numbers and the complicated relationship between PPV and PPP data.
1 Introduction
The current understanding of the interstellar medium (ISM) is that it is a multiphase environment, consisting of gas and dust, which is both magnetized and highly turbulent (Ferriere 2001; McKee & Ostriker 2007). In particular, magnetohydrodynamic (MHD) turbulence is essential to many astrophysical phenomena such as star formation, cosmic ray dispersion, and many transport processes (see Elmegreen & Scalo 2004; BallesterosParedes et al. 2007 and references therein). Additionally, turbulence has the unique ability to transfer energy over scales ranging from kiloparsecs down to the proton gyroradius. This is critical for the ISM, as it explains how energy is distributed from large to small spatial scales in the Galaxy.
Observationally, several techniques exist to study MHD turbulence in different ISM phases. Many of these techniques focus on emission measure fluctuations or rotation measure fluctuations (i.e gradients of linear polarization maps) for the warm ionized media (Armstrong et al.1995; Chepurnov & Lazarian 2010; Gaensler et al. 2011; Burkhart, Lazarian & Gaensler 2012) as well as spectroscopic data and column density maps for the neutral warm and cold media (Spangler & Gwinn 1990; Padoan et al. 2003), For studies of turbulence, spectroscopic data has a clear advantage in that it contains information about the turbulent velocity field as well as the density fluctuations. However, density and velocity are entangled in PPV space, making the interpretation of this type of data difficult. For the separation of the density and velocity fluctuations, special techniques such as the Velocity Coordinate Spectrum (VCS) and the Velocity Channel Analysis (VCA) have been developed (Lazarian & Pogosyan 2000, 2004, 2006, 2008).
Most of the efforts to relate observations and simulations of magnetized turbulence are based on obtaining the spectral index (i.e. the loglog slope of the power spectrum) of either the density and/or velocity (Lazarian & Esquivel 2003; Esquivel & Lazarian 2005; Ossenkopf et al. 2006). However, the power spectrum alone does not provide a full description of turbulence, as it only contains the Fourier amplitudes and neglects information on phases. This fact, combined with complexity of astrophysical turbulence, with multiple injection scales occurring in a multiphase medium, suggests researchers need additional ways of analyzing observational and numerical data in the context of turbulence. In particular, these technique studies are currently focused into two categories:

Development: Test and develop techniques that will complement and build off the theoretical and practical picture of a turbulent ISM that the power spectrum presents.

Synergy: Use several techniques simultaneously to obtain an accurate picture of the parameters of turbulence in the observations.
In regards to the first point, in the last decade there has been substantial progress in the development of techniques to study turbulence. Techniques for the study of turbulence can be tested empirically,
using parameter studies of numerical simulations,
or with the aid of analytical predictions (as was done in the case of the VCA). In the former, the parameters to be varied (see Burkhart & Lazarian
2011) include the Reynolds number, sonic and Alfvénic Mach number, injection scale, equation of state, and, for studies of molecular clouds,
should include radiative transfer and selfgravity (see Ossenkopf 2002; Padoan et al. 2003; Goodman et al. 2009).
Some recently developed techniques include the application of
probability distribution functions (PDFs), wavelets, spectral correlation function (SCF),
The latter point in regard to the synergetic use of tools for ISM turbulence is only recently being attempted as many techniques are still in developmental stages. However, this approach was used in Burkhart et al. (2010), which applied spectrum, bispectrum and higher order moments to HI column density of the SMC. The consistency of results obtained with a variety of statistics, compared with more traditional observational methods, made this study of turbulence in the SMC a promising first step.
The current paper falls under the category of “technique development.” In particular, we investigate the utility of dendrograms in studying the hierarchical structure of ISM clouds. It has long been known that turbulence is able to create hierarchical structures in the ISM (Scalo 1985, 1990; VazquezSemadeni 1994; Stutzki 1998), however many questions remain, such as what type of turbulence is behind the creation of this hierarchy and what is the role of selfgravity and magnetic fields? Hierarchical structure in relation to this questions is particularly important for the star formation problem (Larson 1981; Elmegreen & Elmegreen 1983; Feitzinger & Galinski 1987; Elmegreen 1999; Elmegreen 2011).
Early attempts to characterize ISM hierarchy utilized tree diagrams as a mechanism for reducing the data to hierarchical “skeleton images” (see Houlahan & Scalo 1992). More recently, dendrograms have been used on ISM data in order to characterize selfgravitating structures in star forming molecular clouds (Rosolowsky et al. 2008 and Goodman et al. 2009). A dendrogram (from the Greek dendron “tree”, gramma “drawing”) is a hierarchical tree diagram that has been used extensively in other fields, particularly in computational biology, and occasionally in galaxy evolution (see Sawlaw & HaqueCopilah 1998 and Podani, Engloner, & Major 2009, for examples). Rosolowsky et al. (2008) and Goodman et al. (2009) used the dendrogram on spectral line data of L1448 to estimate key physical properties associated with isosurfaces of local emission maxima such as radius, velocity dispersion, and luminosity. These works provided a new and promising way of characterizing selfgravitating structures and properties of molecular clouds through the application of dendrogram to CO(J=10) PPV data.
In the current paper we apply the dendrogram to synthetic observations (specifically PPV cubes) of isothermal MHD turbulence in order to investigate the physical mechanisms behind the gas hierarchy. Additionally, we are interested in the nature of the structures that are found in PPV data and how these structures are related to both the physics of the gas and the underlying density and velocity fluctuations generated by turbulence. Simulations provide an excellent testing ground for this problem, as one can identify which features in PPV space are density features and which are caused by velocity crowding. Furthermore, one can answer the question: under what conditions do the features in PPV relate back to the 3D density or PPP cube?
In order to address these questions we perform a parameter study using the dendrogram. We focus on how changing the global parameters of the turbulence, such as the sonic Mach number, Alfvénic Mach number, affects the amount of hierarchy observed, the relationship between the density and velocity structures in PPV, and the number and statistical distribution of dominant emission structures.
The paper is organized as follows. In § 2 we describe the dendrogram algorithm, in § 3 we discuss the simulations and provide a description of the MHD models. We investigate the physical mechanisms that create hierarchical structure in the dendrogram tree and characterize the tree diagrams via statistical moments in § 4. In § 5 we compare the dendrograms of PPP and PPV. In § 6 we discuss applications and investigate issues of resolution. Finally, in § 7 we discuss our results followed by the conclusions in § 8.
2 Dendrogram Algorithm
The dendrogram is a tree diagram that can be used in 1D, 2D or 3D spaces to characterize how and where local maxima merge as a function of a threshold parameter. Although the current paper uses the dendrogram in 3D PPV space to characterize the merger of local maxima of emission, it is more intuitive to understand the 1D and 2D applications. A 1D example of the dendrogram algorithm for an emission profile is shown in Figure 1. In this case, the threshold value is called , and is the minimum amplitude above a merger point that a local maximum must have before it is considered distinct. That is, if a merger point (or node) is given by and a local maximum is given by then in order for a given local max to be considered significant, . If then would merge into and no longer be considered distinct.
For 2D data, a common analogy (see Houlahan & Scalo 1992; Rosolowsky et al. 2008) is to think of the dendrogram technique as a descriptor of an underwater mountain chain. As the water level is lowered, first one would see the peaks of the mountain, then mountain valleys (saddle points) and as more water is drained, the peaks may merge together into larger objects. The dendrogram stores information about the peaks and merger levels of the mountain chain.
For our purposes, we examine the dendrogram in 3D PPV space (see Rosolowsky et al. 2008; Goodman et al. 2009 for more information on the dendrogram algorithm applied in PPV). In the 3D case, it is useful to think of each point in the dendrogram as representing a 3D contour (isosurface) in the data cube at a given level.
Our implementation of the dendrogram is similar to many other statistics that employ a user defined threshold value in order to classify structure. By varying the threshold parameter for the definition of “local maximum” (which we call , see Figure 1), different dendrogram tree diagrams and distributions of local maximum and merger points are created. An example of another statistic that utilizes a density/emission threshold value is the Genus statistic, which has proven useful for studying ISM topology (Lazarian, Pogosyan & Esquivel 2002; Lazarian 2004; Kim & Park 2007; Kowal, Lazarian & Beresnyak 2007; Chepurnov et al. 2008). For the Genus technique, the variation of the threshold value is a critical point in understanding the topology of the data in question.
As sets the definition for “local maximum,” setting it too high will produce a dendrogram that may miss important substructures, while setting it very low may produce a dendrogram that is dominated by noise. The issues of noise and the dendrogram were discussed extensively in Rosolowsky et al. 2008. While the dendrogram is designed to present only the essential features of the data, noise will mask the lowamplitude or high spatial frequency variation in the emission structures. In extreme cases, where the threshold value is not set high enough or the signaltonoise is very low, noise can result in local maxima that do not correspond to real structure. As a result, the algorithm has a built in noise suppression criteria which only recognizes structures that have 4 significance above . Such a criterion has been previously used in data cube analysis, as noise fluctuations will typically produce 1 variations (Brunt et al. 2003; Rosolowsky & Blitz 2005; Rosolowsky et al. 2008).
The algorithm we use is extensively described in Goodman et al. (2009) in the Supplementary Methods section and in Rosolowsky et al. (2008), however we describe its main points here. To produce the dendrogram, we first identify a population of local maxima as the points which are larger than all surrounding voxels touching along the face (not along edges or corners). This large set of local maxima is then reduced by examining each maximum and searching for the lowest contour level that contains only that maximum. If this contour level is less than below the local maximum, that local maximum is removed from consideration in the leaf population. This difference in data values is the vertical length of the “leaves” of the dendrogram. Once the leaves (local maxima) of the dendrogram are established, we contour the data with a large number of levels (500 specifically, see Rosolowsky et al. 2008; Goodman et al. 2009). The dendrogram “branches” are graphically constructed by connecting the various sets of maxima at the contour levels where they are joined (see Figure 1 for a 1D example). For graphical presentation, the leaves of the structure tree are shuffled until the branches do not cross when plotting. As a result, the axis of the dendrogram contains no information.
Once the dendrogram is created, there are multiple ways of viewing the information it provides such as:

A tree diagram (the dendrogram itself).

3D viewing of the isocontours and their connectivity in PPV space.

A histogram of the dendrogram leaf and node values (i.e. intensities), which can then be further statistically analyzed.
We note that this third point is a novel interpretation of the dendrogram that we develop in this work. Here, the histogram will be composed of intensity values important to the hierarchical structure of the image. We define a distribution which includes the intensity values of the leaves (i.e. local maximum), denoted by , and intensity values of the nodes, denoted with . This interpretation is visualized in Figure 3 and is further described in Section 4.
The purpose of this paper is to use dendrograms to characterize the observed hierarchy seen in the data. While turbulence has often been cited as the cause of the observed hierarchical structure in the ISM (Stutzki 1998), it is unclear to what extent magnetic fields, gas pressure, and gravity play roles in the creation of ISM hierarchy even though these parameters are known to drastically change the PDF and spectrum of both column density and PPV data (see Falgarone 1994; Kowal, Lazarian & Beresnyak 2007; Tofflemire et al. 2011).
3 Data
We generate a database of sixteen 3D numerical simulations of isothermal compressible (MHD) turbulence by using the MHD code of Cho & Lazarian 2003 and vary the input values for the sonic and Alfvénic Mach number. We briefly outline the major points of the numerical setup.
The code is a thirdorderaccurate hybrid essentially nonoscillatory (ENO) scheme which solves the ideal MHD equations in a periodic box:
(1)  
(2)  
(3) 
with zerodivergence condition ,
and an isothermal equation of state , where
is the gas pressure.
On the righthand side, the source term is a random
largescale driving force. Although our simulations are ideal MHD, diffusion is still present in the form of numerical resistivity
acting on small scales.
The scale at which the dissipation starts to act is defined by the numerical diffusivity
of the scheme
We drive turbulence solenoidally
We divide our models into two groups corresponding to subAlfvénic and superAlfvénic turbulence. For each group we computed several models with different values of gas pressure (see Figure 2) falling into regimes of subsonic and supersonic. We ran one compressible MHD turbulent model with selfgravity at 512 resolution and a corresponding case for the same initial values of pressure and mean magnetic field without selfgravity (models 15 and 16).
We solve for the gravitation potential () using a Fourier method similar to that described in Ostriker et al. 1999. In this case, Equation 2 now has a term on the right hand side. The gravitational kernal used to provide a discrete representation of the Poisson equation is:
(4) 
We can set the strength of selfgravity by changing the physical scaling of the simulations, i.e. by changing the size scaling, cloud mass, and crossing time, which effectively changes the relation between in physical units and code units, which we term in Figure 2. We set these values to give a global virial number of , which is in the upper limit of the observed values found in GMCs (see Kainulainen et al. 2011). We choose a high virial value to investigate the minimum effect gravity might have on the hierarchical structure of clouds. More information about the scalings and their relation to the virial parameter can be found in Section 6.
The models are listed and described in Figure 2. Here we list the rootmeansquared values of the sonic and Alfvénic Mach numbers calculated in ever cell and then averaged over the box.
We use density and velocity perpendicular to the mean magnetic field in order to create fully optically thin synthetic PPV data cubes, although we also investigate dendrogram for other LOS orientations. The PPP and synthetic PPV cubes are all normalized by the mean value, i.e. , unless otherwise stated. Varying the optical depth will be done in a later work. We create cubes with a given velocity resolution of 0.07, which is approximately ten times smaller than the rms velocity of the simulation ( unity). For reference, the sound speed of the simulations varies from for our most subsonic to most supersonic simulations. PPV cubes are created by reorganizing density cubes into channel bins based on given velocity intervals. Additional discussion on comparing the simulations to observations is found in Section 6.
4 Characterizing Hierarchy and Structures Created by Turbulence
We applied the dendrogram algorithm on synthetic PPV cubes with various sonic and Alfvénic Mach numbers. An example of how the the tree diagram output changes with threshold value is shown in Figure 3. The top row of Figure 3 shows the isosurfaces with the colors relating back to the colors in the corresponding dendrogram shown in the middle row. As the threshold intensity value (which, shown here with a black line, sets the definition of the local maximum or “leaves of the tree”) increases, structures in the dendrogram begin to merge with each other. The leaf and branch length and number of structures provide information on the hierarchical nature of the PPV cube. The bottom row of Figure 3 shows the histograms of the dendrogram distribution of intensities (leaves and nodes). The red line is a reference line at intensity level 25. This distribution also changes with changing threshold value, as leaves merge with one another and the hierarchy changes.
In the next subsections, we investigate the effects of the compressibility, magnetization, and selfgravity on the number of structures, amount of hierarchical structure, and moments of the dendrogram distribution. We define a hierarchical dendrogram as one which has many segments on its paths and hence many levels above the root.
4.1 Sonic and Alfvénic Mach Numbers
Leaf and Branch Counting
We computed the dendrogram for all synthetic nonself gravitating PPV cubes with varying threshold values. Figure 4 top shows how the total number of structures (i.e., dominant emission contours including dendrogram “leaves and branches”) changes as we change . We plot the total number of structures vs. on a logarithmic scale (i.e. Log N vs Log ) for simulations with three differing values of sonic Mach number (the M7, M3, M0.7 models) and two values of Alfvénic Mach number. The left panel shows subAlfvénic models and the right shows superAlfvénic models. Error bars are created by taking the standard deviation between different time snapshots. We note that power law tails can be seen at values of past the mean value (i.e. past log ). We over plot the values of the slopes with solid black lines for reference. The symbols with no lines through them in Figure 4 are for PPV cubes with LOS taken perpendicular to the mean magnetic field. We tested our results for LOS taken parallel to the mean magnetic field and found similar results. We show examples of the results for cubes with LOS taken parallel to the mean magnetic field in Figure 4 for the M7 models using symbols connected with solid lines. The M7 models are highly supersonic and therefore will show the most deviation along different sight lines.
When is at or slightly above the mean value of the data cubes, there is little difference in the number of structures between simulations of different sonic Mach number. This is surprising, since the structures seen in subsonic turbulence are very different from the supersonic case. In the regime where is at the mean value, we are sampling most of the PPV cube emission and therefore are not sensitive to the differences seen at larger threshold values, which will merge low intensity structures. Once we increase beyond the mean however, the number of structures between the subsonic (black plus signs) and supersonic simulations (red stars and green diamonds) rapidly diverges. The larger number of structures in the supersonic case is a result of shocks creating higher intensity values in the PPV cubes. Additionally, the slopes in the subsonic simulations are much steeper as compared with the supersonic simulations since the number of structures the dendrogram considers significant at a given threshold value rapidly falls off to zero. Subsonic models have fewer significant emission contours since they do not have density enhancements created by shocks and therefore the density/velocity contrast between subsonic and supersonic turbulence becomes clear at higher threshold values. The higher the Mach number, the more small scale intensity enhancements we expect to see.
As increases, differences between supersonic (M3) and very supersonic (M7) cases become more apparent, as the slopes for the M3 case are steeper. This is because interacting shocks in the M7 case are much stronger, and hence there is more contrast in the emission contours. Thus, as we increase , the structures merge more rapidly for lower values of the sonic Mach number.
Comparison between the left and right top panels of Figure 4 shows that the magnetic field also affects the number of structures and the trend with the threshold value. When is low, the superAlfvénic cases (right panel) show slightly more structures than the subAlfvénic ones (left panel). However, as increases, the number of structures decreases more rapidly in the case of superAlfvénic turbulence, i.e. the slopes are steeper. Interestingly, when we take a sightline parallel to the mean magnetic field (green symbols with the solid line drawn through) we still see the same trend with Alfvénic Mach number. We suspect that these results are generally independent of time evolution in driven turbulence, as clump lifetimes have been shown to be determined by turbulence motions on large scales, rather then limited by their crossing times (see FalcetaGonçalves & Lazarian 2011)
The bottom plot of Figure 4 shows the number of segments from root to leaf on the largest branch vs. the threshold parameter . A test of hierarchy is to count the number of segments along the largest branch, from leaf to root. Similar to what was shown in the top figure, the sonic Mach number has a strong relation to the amount of hierarchical structure created in the gas. Higher sonic Mach number yields more shocks, which in turn produce more high density clumps and more hierarchical structures in PPV space. Interestingly, the magnetic field seems to play an even stronger role in the hierarchical branching than shocks. Comparison between the axis values of the left and right plots reveals that a larger Alfvénic Mach number allows for more hierarchical structure in the PPV dendrogram. In the case of superAlfvénic turbulence, magnetization is low and hence the structures created are closer to that of hydrodynamic turbulence, which is well known to show fractal behavior and hierarchical eddies. As turbulence transitions to subAlfvénic, it becomes magnetically dominated with fewer degrees of freedom.
We illustrate the results of Figure 4 as a cartoon shown in Figure 5. The cartoon shows various configurations of two ISM clouds. We label the first cloud as case A, which is a cloud with a global Alfvénic Mach number . The second cloud, which is a cloud with a global Alfvénic Mach number , is labeled as cases B and C, indicating shock compression parallel and perpendicular to the mean magnetic field, respectively. Both clouds are assumed to have the same supersonic value of the sonic Mach number. Case A shows hierarchical structure forming in clumps that are not affected strongly by the magnetic field. The clumping and hierarchy is due to compression via shocks and the shredding effect of hydrodynamic turbulence. The turbulent eddies for cloud A can evolve with a full 3D range of motion and have more degrees of freedom as compared with turbulence in the presence of a strong magnetic field. In light of this, when we consider strong magnetization, we must now investigate the effects of shock compressions oriented parallel and perpendicular to the mean magnetic field (case B and C). For shock compression parallel to the field lines (case B), the clumps will be confined in the direction perpendicular to the field, and thus the compression will be able to squeeze the clumps, decrease the hierarchy in the gas, and create additional large density contrast. For shock compression perpendicular to the field lines, the magnetic pressure relative to the shock compression is higher, and the clumps will not feel as much of the compression. Furthermore, the strong field creates anisotropy in the eddies, which are stretched along the direction of the mean field line. This limits the range of motion of the eddies, which in turn limits their ability to interact. Thus, in comparison of cloud B/C with cloud A the contrast is higher while hierarchical structure is less.
Statistics of the Dendrogram Distribution
A dendrogram is a useful representation of PPV data in part because there are multiple ways of exploring the information on the data hierarchy. In this section we investigate how the statistical moments of the distribution of the dendrogram tree (see bottom panels of Figure 3 for example) changes as we change the threshold parameter and how these changes depend on the compressibility and magnetization of turbulence. We consider a distribution containing all the intensity values of the leaves and merging intensity contour values in a given dendrogram. The question that forms the basis of our investigation in this section is: Do the moments of the distribution have any dependencies on the sonic and Alfvénic Mach numbers?
The first and second order statistical moments (mean and variance) used here are defined as follows: and , respectively. The standard deviation is related to the variance as: . The third and fourth order moments (skewness and kurtosis) are defined as:
(5) 
(6) 
We calculate the moments of the dendrogram tree distribution while varying our simulation parameter space. In particular, we vary the sonic Mach number, the Alvénic Mach number, and the threshold value. We find the moments vs. the threshold parameter to show linear behavior. As increases, the number of the intermediate intensity values that make up the branches and the hierarchical nesting (i.e. the intensity values between the high intensity local maximum and the low intensity values near the trunk) merge with each other. This effect can be seen visually in Figure 3. Thus, as increases, the mean and variance of the distribution (example show in the bottom of Figure 3) will increase.
We plot the moments vs. with in Figure 6. This figure shows the full range of our simulations with sub and superAlfvénic combinations. Generally, as the sonic Mach number increases so do the moments. We found this trend to be consistent over a range of values, and hence only plot one case here. Error bars, created by taking the standard deviation of the value between different time snapshots of the simulation, generally increase with sonic number as the fluctuations become increasingly stochastic and shock dominated. The increase of the moments of is related to the compressibility of the model and more supersonic cases display more prominent clumpy features, which drive up both the average and the variation from average. The tails and peak of the distribution also become increasingly skewed and kurtotic towards higher values of intensity and the distribution becomes increasingly peaked around the mean value.
It is interesting to note that a strong dependency on the magnetization of the model exists, particularly as the sonic number goes up. The subAlfvénic simulations show increased moments, which implies that they exhibit more contrast (mean value is higher) and more skewed/kurtotic distributions in their gas densities.
In the above analysis, the distribution included all leaves and branches of the dendrogram tree. We could further cut the tree into its respective branches and leaves and analyze the distributions separately, which provides additional constraints on the parameters. We investigated the statistical moments on the histograms of the branch lengths, leaf lengths, and leaf intensities and found the trends discussed above to be consistent with the results of Figure 6, and hence do not include the plots.
4.2 SelfGravity
Leaf and Branch Counting
The issues of the importance of selfgravity in simulations for comparisons with the molecular medium have been raised by a number of authors (Padoan et al. 2001; Li et al. 2004; Goodman et al. 2009; Federrath et al. 2010). While selfgravity is known to be of great importance to accretion disk physics and protostellar collapse, its role in the outer regions of GMCs and in diffuse gasses is less obvious. The dendrogram can potentially be used to explore the relative role between turbulence and gravity regarding both the structure of the hierarchy and the distribution of dominant emission contours.
Figure 7 shows tree diagrams at constant threshold value =45 for subAlfvénic supersonic simulations with and without selfgravity (models 15 and 16 in Figure 2). A large value of is used in order to not over crowd the dendrogram with branches. It is clear that the case with selfgravity (the tree diagram in the right panel of Figure 7) has a dendrogram with more significant structure and hierarchical nesting. Both models are supersonic with approximately the same sonic Mach number. Visually, Figure 7 makes it clear that selfgravity, even with a virial number of 90, is an important factor for the creation of hierarchical structure.
In order to quantify the number of structures and hierarchical nesting seen in Figure 7 we repeat the analysis performed in Figure 4. We show the number of structures vs. (top panel) and the number of segments on the largest branch of the tree vs. (bottom panel) in Figure 8. It is clear that the case with no selfgravity (black crosses) shows less overall structure (leaves and nodes, top plot) and less hierarchical structure (bottom plot) compared with the cases with selfgravity (red asterisk). Even with a high virial parameter, the supersonic selfgravitation simulation has significantly more nested structures and more contours considered to be areas of significant emission than supersonic simulation without gravity. Interestingly, the powerlaw slopes seen in the top plot of Figure 8 are not significantly different between the selfgravitating and non selfgravitating cases. In fact, the simulations appear to be only shifted vertically and no substantial difference is seen regarding the slopes. Therefore it is possible that the change in the slope, which was observed in Figure 4, can be attributed to the sonic and Alfvénic nature of the turbulence directly, and is not substantially affected by the inclusion of weak selfgravity. More work should be done in the future to determine how the dendrogram is influenced by the presence of strong selfgravitating flows. We discuss further the use of the dendrogram for quantifying the relative importance of selfgravity and turbulence in the discussion section (section 7).
Statistics of the Dendrogram Distribution
We show how selfgravity affects the moments of the dendrogram distribution as we vary in Figure 9 for models 15 and 16. Higher levels of selfgravity show increases in all four moments over a range . The presence of gravity increases the mean value and variance of emission contours as well as skews the distribution towards higher values due to the presence of attracting regions.
The results of Section 4.2 can be applicable to molecular clouds where, given a relatively constant value of the sonic Mach number, regions under the influence of stronger selfgravity could possibly be identified with the use of the analysis presented in Figures 8 and 9. This could be done by comparing the dendrogram moments and number of structures with different between various sections of a cloud. Areas that show increased moments and increased hierarchical clumping with no changes in the slopes of structures vs. , may indicate changes in selfgravity. Knowing the Mach number of these regions will greatly increase the reliability of such analysis and fortunately, other techniques exist to find these (see for example Burkhart et al. 2010; Esquivel & Lazarian 2011; Kainulainen & Tan 2012). We further discuss applying the dendrogram on observations in section 7.
5 Dendrograms of PPP vs. PPV
The issue of interpreting structures seen in PPV space has vexed researchers for over a decade (see Pichardo et al. 2000). How the structures in PPV translate to PPP depends on many factors, most importantly the nature of the turbulent environment. The dendrogram presents a unique way of studying how the hierarchy of structures seen in density space (PPP) relate to PPV space via simulations.
For turbulent clouds, it is never the case that the structures in PPV have a onetoone correspondence with the density PPP, although this assumption may be more appropriate for some environments than others. We show a simple example illustrating this in Figure 10, which shows two synthetic subsonic superAlfvénic PPV data cubes (left), which share the same velocity distribution but have different density distributions. The bottom left PPV cube has constant density/column density, while the top left PPV cube’s corresponding turbulent density cube (PPP cube) is shown on the right.
Interestingly, the bottom left PPV cube has a very similar level of structure as compared with the top PPV cube, despite the fact that the column density of the bottom cube is constant. This points out the well known fact that there is not a onetoone correspondence with PPV and PPP space. In fact, in this example (a subsonic model) most of the structures seen in PPV are due to the velocity rather than the density. Figure 10 illustrates the dominance of velocity in the subsonic case in the bottom PPV cube. Fluctuations in PPV here are entirely driven by the turbulent velocity field.
To illuminate this point further, Figure 11 shows PPP and PPV dendrograms for supersonic turbulence with (model 15, middle) and subsonic turbulence with (model 1, bottom). We also show the corresponding isosurfaces for the supersonic case in the top row. Comparing PPV and PPP should be done with care as they are different spaces. We increased the value of until the PPP dendrogram becomes mostly leaves, i.e. they have little hierarchy. The leaves are reached at . We took the corresponding optically thin PPV cube and applied the dendrogram with the same threshold value. If the dominant emission is due to density than the leaves should be similar for both PPV and PPP. All PPV and PPP cubes are normalized to have a mean value of unity.
Interestingly, the supersonic turbulence dendrogram for density looks very similar to the corresponding PPV dendrogram for the same at the level of the leaves. For the subsonic case we see that the dendrogram of density and PPV look nothing alike (same ). In this case, the velocity field dominates PPV space. Hence, we do not show the isosurfaces for the subsonic case. In supersonic turbulence, the highest density peaks correspond to the highest intensity fluctuation in the PPV. This implies that if one knows that the turbulence in question is supersonic, the structures in PPV space at the level of the leaves can generally be interpreted as 3D density structures. However, if the turbulence is subsonic in nature this assumption is not appropriate.
6 Application
The dendrogram shows dependencies on the parameters of turbulence that are important for studies of star forming regions and the diffuse ISM. When analyzing a particular data set, one should keep in mind that comparisons between the observational and scaled numerical data, or comparisons between different clouds or objects within the same data set, are the most useful means of extracting these parameters.
Our simulations can be scaled to observations by specifying the physical size of the simulation volume, the isothermal sound speed of the gas, and mass density. In particular, appropriate scalings must be made in the case of the selfgravitating simulation. We can set the strength of selfgravity by changing the physical scaling of the simulations, i.e. by changing the box size, cloud mass, and crossing time. In this case, the relevant scaling between physical and code units are the size scaling (), the velocity scale factor (), the time scale factor (), and the mass scale factor (). These are given by:
(7) 
where is the physical length of the box, and is in code units and is equal to unity.
(8) 
where and are the observed and code sound speeds, respectively,
(9) 
and
(10) 
where is the observed cloud mass and is the average density of the simulation (unity).
Using these relations between the code units and physical units we can define a relation between in code units and physical units, and then relate this to the free fall time. In this case,
(11) 
Here we use a value of (see Figure 2). In this case, our free fall time (in code units and therefore using of unity) is . The dynamical time (in code units with simulations box size of unity and rms velocity) of our selfgravitating simulation is . The ratio of the free fall time to the dynamical time of the simulation with self gravity gives the global virial parameter . Additional information on scaling simulations to observations can be found in Hill et al. 2008.
We include the effects of changing the velocity resolution, thermal broadening, and smoothing in the next subsection.
6.1 Smoothing
We investigate how smoothing and data resolution affect the dendrogram. When dealing with observational data, one must always consider the effect that the telescope beam smoothing will have on the measurement. The observations are rarely done with pencil beams, and the measured statistics change as the data is averaged. We expect the effect of smoothing to depend on a dimensionless number, namely, the ratio of the size of the turbulence injection scale to the smoothing scale.
We apply the same technique that was applied in the previous sections, i.e. exploring the number of structures and moments of dendrogram tree statistics. However, we now include a boxcar smoothing kernel (truncating the edges). We expect that smoothing will affect supersonic turbulence and cases of high selfgravity the most. In this case, shocks and smallscale gravitational clumps become smoothed out and more difficult for the algorithm to identify. In the subsonic or low gravity cases, smoothing makes less of a difference, since the gas is already diffuse and less hierarchical.
We show how the moments and number of structures changes with smoothing size (in pixels) in Figure 12. One could also discuss smoothing beam size in terms of the injection scale of the turbulence. For instance, 7 pixel smoothing represents a beam scale that is 30 times smaller than our injection scale of turbulence.
We found that in general, subsonic and transonic turbulence are not as affected by smoothing compared to highly supersonic models. In light of this, we plot the moments and number of structures vs. for different smoothing degrees for a highly supersonic model with (M7 models) in Figure 12. Two panels show different Alfvénic regimes with the axis the same for both for ease of comparison. Black lines indicate no smoothing, while red and blue indicate three and seven pixel smoothing, respectively. Error bars are produced by taking the standard deviation between different time snapshots of the simulations with well developed turbulence.
As smoothing increases for this supersonic model, we see that the values of the moments as well as the total number of structures decreases. However, even out to seven pixel smoothing the differences between the Aflvénic cases is evident in the mean and variance, respective of the error bars. Furthermore, the trends with the threshold parameter do not change when we introduce smoothing, which gives us further confidence that this technique can be applied to the observational data. Other than the change in amplitude, the trends remain the same as what was seen in Section 4.
6.2 Velocity Resolution and Thermal Broadening
In addition to smoothing, we must also consider the effects of velocity resolution. As the velocity resolution changes in PPV space, so do the structures observed. We investigated how the number of structures in the dendrogram distribution change when we vary the velocity resolution. We find that the number of substructures drops as the velocity resolution decreases, from several hundreds to several dozen when changing the velocity resolution from to . This effect corresponds to the channel sampling dropping from 60 to 15 channels. This may provide too low a number of statistics in the dendrogram distribution to look at the moments, however the general trends with the physical parameters stay consistent with section 4.
An additional observational consideration that should be made regards thermal broadening. The bulk of the this paper focuses on the effects of turbulence and magnetic fields in the creation of hierarchical structure in ISM clouds. However, for warm subsonic or transonic gas, thermal broadening effects should also be considered. Convolution with a thermal broadening profile (i.e. a Gaussian) will smooth out the velocity profiles in these cases. However, we expect the structures seen in supersonic gas to not be affected since turbulence dominates the line broadening.
To demonstrate this we convolve the line profiles of four of our simulations (models M0.7, M3, M7, M8) with Gaussian profiles to mimic the effects of thermal broadening. The thermal Gaussian has FWHM given as the ratio of the turbulent line width to the sonic Mach number.
We repeat the analysis done in Figure 4, using the models which now include thermal broadening, in Figure 13. We show the fitted slopes of the number of structures vs. from Figure 4 (top panel) as solid black lines accompanied by the numerical value of the slope. These black lines serve as a reference for which to compare with simulations without thermal broadening with the simulations that now include thermal broadening (colored symbols).
As expected, the supersonic dendrograms are mostly unaffected by the inclusion of thermal broadening, since turbulent dominates the line profits. We find the number of structures vs. and corresponding slopes (top plot) for the supersonic simulations to be very similar to those shown in Figure 4, with only a slight shallowing of the slopes. Additionally the amount of hierarchical branching (bottom plot) is also similar. However the subsonic PPV cube intensities are more drastically lowered as a result of the convolution with a broader Gaussian and thus the values of must be lowered as well. Additionally, the slopes seen for the subsonic simulations in the top plot of Figure 13 are shallower then the slopes of the simulations with no thermal broadening included (comparison with solid lines referencing Figure 4). Therefore when examining the the dendrograms of warm gas (such as HI and H) thermal broadening effects must also be taken into account.
7 Discussion
Hierarchical tree diagrams are finding more applications in interstellar studies, not only to locate clumps and calculate their properties, but also for characterizing properties of the physics present in interstellar and molecular gas. We used dendrograms to analyze how turbulence, magnetic fields and selfgravity shape the amount of structure and gas hierarchy in isothermal simulations. We examined the changes in the distribution of the dendrogram as we vary the threshold parameter . This is analogous to changing the corresponding threshold parameter in other techniques that rely on contouring thresholds, e.g. in the Genus analysis (see Chepurnov et al. 2009). By varying we obtained a new outlook on the technique; in particular, we found that the dendrogram distribution and hierarchy have a strong dependency on the magnetization and compressibility of the gas and are sensitive to the amount of selfgravity.
7.1 The Hierarchical Nature of MHD Turbulence
The number of structures and the amount of hierarchy formed by MHD turbulence has interesting implications for the evolution of ISM clouds and for the star formation problem. In Section 4 we found that more hierarchical structure and more overall structure was created in the presence of supersonic turbulence. We also found that the inclusion of selfgravity enhanced these trends. The magnetic field also had a strong influence in the creation of hierarchical nesting in PPV space. The relationship between the magnetization and the cloud dynamics is still not well understood, especially in regards to star formation. Star forming clouds are known to be hierarchical in nature and magnetized, but their exact Alfvénic nature is less clear. The results from this work seem to suggest that very hierarchical clouds might tend towards being superAlfvénic. Several authors have suggested a variety of evidence for molecular clouds being superAlfvénic. This includes the agreement of simulations and observations of Zeemansplitting measurements, vs. relations, vs. relations, statistics of the extinction measurement etc. (Padoan & Nordlund 1999; Lunttila et al. 2008; Burkhart et al. 2009; Crutcher et al. 2009; Collins et al. 2012). Furthermore a study done by Burkhart et al. 2009 found that even in the presence of globally subAlfvénic turbulence, the highest density regions tend towards being locally superAlfvénic. This suggests that even in the case of globally subAlfvénic turbulence, the densest regions might be superAlfvénic. It is interesting that the dendrogram technique also points to superAlfvénic turbulence as an avenue for hierarchical structure creation. This provides motivation for the dendrogram technique to be applied to the observational data with varying threshold value in order to see how the nature of the hierarchical structure and total number of structures change in the observations.
7.2 Characterizing selfgravity and obtaining the Sonic and Alfvénic Mach Numbers from the Observations
We provided a systematic study of the variations of the dendrogram thresholding parameter with the sonic and magnetic Mach numbers. We also included a simulation with weak selfgravity (global virial number of 90) in order to investigate the influence gravity has on the observed hierarchy. While real molecular clouds generally have virial numbers much lower than this value, we showed that the dendrogram is highly influenced by the inclusion of even weak selfgravity. The sonic and Alfvénic Mach numbers, as well as the virial parameter, are critical for understanding most processes in diffuse and molecular gas, including the process of star formation. Thus, the dendrogram, with its sensitivity to these parameters, provides a possible avenue of obtaining the characteristics of turbulence and the relative importance of selfgravity to turbulence in the ISM. For example, the dendrogram, coupled with a virial analysis, was already used to compare the relative importance of selfgravity within the L1148 GMC cloud using CO data in Rosolowsky et al. 2008 and Goodman et al. 2009.
We view this work as a springboard for applying this technique to the observational data, which is why we addressed the issues of smoothing and thermal broadening in Section 6. It is clear that the relation between the dendrogram structures, their statistics and the thresholding value , as explored in this work (for examples see Figures 4 and 8), do not yield universal numbers, i.e. they depend on the observational characteristics such as the velocity resolution and beams smoothing etc. Therefore, in order to apply the dendrogram to the observations, we feel that it is necessary to define a fiducial dendrogram for the data where some information is known about turbulence, and in the case of selfgravitating clouds, the virial parameter. The fiducial dendrogram can then be compared to other regions within the same data, which all contain the same observational constraints. Similarly, for comparison of simulations and observations, the simulated observations must be tailored to the resolution of the observational data.
In order to obtain a fiducial region for further dendrogram analysis and to increase the reliability of the parameters found via the dendrogram, it is advantageous to combine different techniques designed to investigate ISM turbulence. For instance, by applying the VCA and VCS techniques to PPV data (see Lazarian 2009 for a review), one can obtain the velocity and density spectra of turbulence. While these measures are known to depend on and to a lesser degree on (see Beresnyak, Lazarian & Cho 2005, Kowal, Lazarian & Beresnyak 2007, Burkhart et al. 2009), the utility of the spectra is not limited to measuring these quantities. Spectra provide a unique way to investigate how the energy cascades between different scales, and show whether comparing observations with the simulations with a single scale of injection is reasonable.
The analysis of the anisotropies of correlations using velocity centroids provides an insight into media magnetization, i.e., it provides (Lazarian et al. 2002, Esquivel & Lazarian 2005), which is complementary to the dendrogram technique. Studies of the variance, skewness and kurtosis of the PDFs (see Kowal, Lazarian & Beresnyak 2007; Burkhart et al. 2009, 2010, 2012) provides measures of the sonic Mach number . Similarly, Tsallis statistics measures (Esquivel & Lazarian 2010, Tofflemire et al. 2011) provide additional ways of estimating both and . We feel the approach to obtaining these parameters should be conducted with synergetic use of multiple tools, as was done in Burkhart et al. 2010 on the SMC. The dendrogram is a unique tool as it can classify the hierarchical nature of the data and that it should be added to a standard set of statisticaltools for studies of ISM data.
All these techniques provide independent ways of evaluating parameters of turbulence and therefore their application to the same data set provides a more reliable estimate of key parameters such as compressibility, magnetization, and degree of selfgravity. Dendrograms have some advantages over other statistics designed to search for turbulence parameters, in that one can analyze the resulting tree diagram in many different ways, as highlighted in here and in previous works. These include finding local maxima, calculating physical properties of dominate emission, exploring how those clumps are connected in PPV, varying the threshold and calculating moments and level of hierarchy. Of course, one should keep in mind that the medium that we investigate observationally is far from simple. Multiple energy injection sources, for example, are not excluded. Thus, obtaining a similar answer with different techniques should provide us with additional confidence in our results.
Finally, we should stress that for studies of astrophysical objects, the dendrogram and other statistical measures can be applied locally to different parts of the media. For instance, Burkhart et al. (2010) did not characterize the entire SMC with one sonic Mach number. Instead, several measures were applied to parts of the SMC in order to obtain a distribution of the turbulence in the galaxy. A similar local scale selection was applied also to the SMC in Chepurnov et al. (2008) using the Genus technique. Correlating the variations of the turbulence properties with observed properties of the media, e.g. star formation rate, should provide insight into how turbulence regulates many key astrophysical processes.
8 Summary
We investigated dendrograms of isothermal MHD simulations with varying levels of gravity, compressibility and magnetization using multiple values of the threshold parameter . The dendrogram is a promising tool for studying both gas connectivity in the ISM as well as characterizing turbulence. In particular:

We propose using statistical descriptions of dendrograms as a means to quantify the degree of hierarchy present in a PPV data cube.

Shocks, selfgravity, and superAlfvénic turbulence create the most hierarchical structure in PPV space.

The number of dendrogram structures depends primarily on the sonic number and selfgravity and secondarily on the global magnetization.

The first four statistical moments of the distribution of dendrogram leaves and nodes have monotonic dependencies on the inclusion of selfgravity and the sonic and Alfvén Mach numbers over a range of .

The dendrogram provides a convenient way of comparing PPP to PPV in simulations. Density structures are dominant in supersonic PPV and not in subsonic. Thus, it is more justifiable to compare PPV to PPP when the gas is known to be supersonic.
Footnotes
 affiliation: Astronomy Department, University of Wisconsin, Madison, 475 N. Charter St., WI 53711, USA
 affiliation: Astronomy Department, University of Wisconsin, Madison, 475 N. Charter St., WI 53711, USA
 affiliation: HarvardSmithsonian Center for Astrophysics, 60 Garden Street, MS78, Cambridge, MA 02138
 affiliation: University of British Columbia, Okanagan Campus, 3333 University Way, Kelowna BC V1V 1V7, Canada
 The similarities between VCA and SCF are discussed in Lazarian (2009).
 ENO schemes, such as the one employed here, are considered to be generally low diffusion ones (see Liu & Osher 1998; Levy, Puppo & Russo 1999, e.g )
 The differences between solenoidal and compressive driving is discussed more in Federrath et al. (2008). One can expect driving in the ISM to be a combination of solenoidal and compressive, however both types of driving will produce shocks on a range of scales, which is what we study here.
References
 Armstrong, J. W., Rickett, B. J., Spangler, S. R., 1995, ApJ, 443, 209
 BallesterosParedes, J., Klessen, R. S., Mac Low, M.M., & VazquezSemadeni, E. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil (Tucson, AZ: Univ. of Arizona Press), 63
 Beresnyak, A., Lazarian, A., Cho, J., 2005, ApJ, 624, 93
 Berkhuijsen E., Fletcher, A., 2008, MNRAS, 390, 19
 Brunt, C. M., Kerton, C. R., & Pomerleau, C., 2003, ApJS, 144, 47
 Brunt, C., & Heyer, M., 2002, ApJ, 566, 27
 Brunt, C., 2010, A&A, 513, 67
 Burkhart, B., FalcetaGoncalves, D., Kowal, G., Lazarian, A., 2009, ApJ, 693, 250
 Burkhart, B., Stanimirovic, S., Lazarian, A., Grzegorz, K., 2010, ApJ, 708, 1204
 Burkhart, B., Lazarian, A., Gaensler, B., 2012, ApJ, 708, 1204
 Burkhart, B., Lazarian, A., 2012, ApJ, 755, 19
 Burkhart, B., & Lazarian, A., 2011, IAUS, 274, 365
 Chepurnov, A., Gordon, J., Lazarian, A., & Stanimirovic, S.,2008, ApJ, 688, 1021
 Chepurnov, A., & Lazarian, A. 2009, ApJ, 693, 1074
 Chepurnov & Lazarian, 2010, ApJ, 710, 853
 Cho, J. & Lazarian, A. 2003, MNRAS, 345, 325
 Collins et al., 2012, ApJ, 750, 13
 Crutcher, R., Hakobian, N., Troland, T., 2009, ApJ, 692, 844
 Elmegreen, B. G., & Elmegreen, D. M. 1983, MNRAS, 203, 31
 Elmegreen, B. G., & Scalo, J. 2004, ARA&A, 42, 211
 Elmegreen, B. G., 1999, ApJ, 527, 266
 Elmegreen, B. G., 2011, Star Formation in the Local Universe, Eds. C. Charbonnel & T. Montmerle, EAS Publications Series, 51,31
 Esquivel, A, & Lazarian, A., 2005, ApJ, 295,479
 Esquivel, A., & Lazarian, A., 2010, ApJ, 710, 125
 FalcetaGonçalves, D., Lazarian, A., 2011, ApJ, 735, 99
 Falgarone, E., Lis, D. C., Phillips, T. G., Pouquet, A., Porter, D. H., & Woodward, P. R. 1994, ApJ, 436, 728
 Ferriere, K., 2001, RvMP, 73, 1031
 Federrath, C., et al., 2008, ApJ, 688, 79
 Federrath, C., et al., 2010, ApJ, 713, 269
 Feitzinger, J. V., & Galinski, T. 1987, A&A, 179, 249
 Frisch, U., 1995, Turbulence, Univ. of Cambridge Press
 Gaensler et al., 2011, Nature, 478, 214217
 Gill, A.G., & Henriksen, R.N., 1990, ApJ, 365, L27
 Goodman et al., 2009, Nature, 457, 63
 Goodman, A., Pineda, J., Schnee S., 2009, ApJ, 692, 91
 Haverkorn, M., & Heitsch, F., 2004, A&A,421, 1011
 Hill, A. et al., 2008, ApJ, 686, 363
 Houlahan P., & Scalo J., 1992, ApJ, 393, 172
 Kang, H., Ryu, D., & Jones, T. W., 2009, ApJ, 695, 1273
 Kainulainen, J., Beuther, H., Banerjee, R., Federrath, C., Henning, T., 2011, A&A, 530, 64
 Kim, S., Park, G., 2007, ApJ, 663, 244
 Kowal, G., Lazarian, A. & Beresnyak, A., 2007, ApJ, 658, 423
 Larson, R.B. 1981, MNRAS, 194, 809
 Lazarian, A., 2004, J. Korean Astron. Soc., 37, 563
 Lazarian, A., 2009, SSR, 143, 357
 Lazarian A., & Esquivel, A., 2003, ApJ, 592, 37
 Lazarian, A. & Pogosyan,D., 2000, ApJ, 537, 720
 Lazarian, A. & Pogosyan, D., 2004, ApJ, 616, 943
 Lazarian,A. & Pogosyan,D., 2006,ApJ, 652, 1348
 Lazarian, A. & Pogosyan,D., 2008,ApJ, 686, 350
 Lazarian, A., Pogosyan, D., & Esquivel, A. 2002, in ASP Conf. Proc. 276, Seeing Through the Dust: The Detection of H I and the Exploration of the ISM in Galaxies, ed. A. R. Taylor, T. L. Landecker, & A. G. Willis (San Francisco: ASP), 182
 Levy, D., Puppo, G. & Russo, G., 1999, Mathematical Modeling and Numerical Analysis, 33, 547
 Li, P. S., Normand, M., Mac Low, M., Heitsch, F., 2004, ApJ, 605, 800
 Liu, X.D. & Osher, S., 1998, Journal of Computational Physics, 141, 1
 McKee, C., Ostriker, E., 2007, , ARA&A, 45, 565
 Ossenkopf, V., & Mac Low, M.M., 2002, A&A, 390, 307
 Ossenkopf, V., Esquivel, A., Lazarian, A., & Stutzki, J. 2006, A&A, 452, 223
 Ostriker, E., Gammie, C., Stone, J., 1999, ApJ, 513, 259
 Padoan, P., & Nordlund, A., 1999, ApJ, 526, 27
 Padoan, P., Rosolowsky, E. W., & Goodman, A. A. 2001, ApJ, 547, 862
 Padoan, P., Nordlund, A., Rognvaldsson, O. E., & Goodman, A. A., 2001, in ASP Conf. Ser. 243, From Darkness to Light: Origin and Evolution of Young Stellar Clusters, ed. T. Montmerle, & P. Andre (San Francisco: ASP), 279
 Padoan, P., Goodman, A., Juvela, Mika, 2003, ApJ, 588, 881
 Padoan, P., Juvela, M., Kritsuk, A. G., & Norman, M. L. 2009, ApJ, 707, L153
 Pichardo et al. 2000, ApJ, 532, 353
 Podani, J., Engloner, A., & Major, A., 2009, Stat. Appl. Genet. Mol. Biol., 22.
 Price, Federrath, & Brunt, 2011, ApJ, 727
 Rosolowsky, E., Goodman, A., Wilner, D., Williams, J.,1999, ApJ, 524, 887
 Rosolowsky, E., & Blitz, L., 2005, ApJ, 623, 826
 Rosolowsky, E. W., Pineda, J. E., Kauffmann, J., & Goodman, A. A. 2008, ApJ, 679, 1338
 Sawlaw W. & HaqueCopilah, S., 1998, ApJ 509, 595
 Spangler, S. R., & Gwinn, C. R. 1990, ApJ, 353, L29
 Toffelmire, B., Burkhart, B., Lazarian, A., 2011, ApJ, 736, 60
 Scalo, J.S. 1985, in Protostars and Planets II, ed. D.C Black and M. S. Matthews, (Tucson: Univ. of Arizona Press), p. 201
 Scalo, J. 1990, in Physical Processes in Fragmentation and Star Formation, eds. R. CapuzzoDolcetta, C. Chiosi, & A. Di Fazio, Dordrecht: Kluwer, p. 151
 Stutzki, J., Bensch, F., Heithausen, A., Ossenkopf, V., Zielinsky, M., 1998, A&A, 336,697
 VazquezSemadeni, E., 1994, ApJ, 423, 681