The Anatomy of the Column Density Probability Distribution Function (N-PDF)

The column density probability distribution function (N-PDF) of GMCs has been used as a diagnostic of star formation. Simulations and analytic predictions have suggested the N-PDF is composed of a low density lognormal component and a high density power-law component, tracing turbulence and gravitational collapse, respectively. In this paper, we study how various properties of the true 2D column density distribution create the shape, or"anatomy"of the PDF. We test our ideas and analytic approaches using both a real, observed, PDF based on Herschel observations of dust emission as well as a simulation that uses the ENZO code. Using a dendrogram analysis, we examine the three main components of the N-PDF: the lognormal component, the power-law component, and the transition point between these two components. We find that the power-law component of an N-PDF is the summation of N-PDFs of power-law substructures identified by the dendrogram algorithm. We also find that the analytic solution to the transition point between lognormal and power-law components proposed by Burkhart, Stalpes&Collins (2017) is applicable when tested on observations and simulations, within the uncertainties. We reconfirm and extend the results of Lombardi, Alves&Lada (2015), which stated that the lognormal component of the N-PDF is difficult to constrain due to the artificial choice of the map area. Based on the resulting anatomy of the N-PDF, we suggest avoiding analyzing the column density structures of a star forming region based solely on fits to the lognormal component of an N-PDF. We also suggest applying the N-PDF analysis in combination with the dendrogram algorithm, to obtain a more complete picture of the global and local environments and their effects on the density structures.


INTRODUCTION
Star formation occurs in dense filamentary structures within molecular environments that are governed by the complex interaction of gravity, magnetic fields, and turbulence (McKee & Ostriker 2007). The initial distribution of the gas density at parsec scales, which is affected by the average density, level of turbulence and magnetic field strength, may determine the AU scale properties of star formation such as the initial mass function (IMF) and the overall star formation rate (Krumholz & McKee 2005;Hennebelle & Chabrier 2011;Padoan & Nordlund 2011a,b;Federrath & Klessen 2012;Mocz et al. 2017).

History
The column density probability distribution function (N-PDF) is a commonly used tool for quantifying the distribution of gas. Simulations and observations have shown N-PDFs to be an important diagnostic of turbulence and star formation efficiency in local star forming clouds (Federrath & Klessen 2012;Collins et al. 2010;Burkhart et al. 2015a;Myers 2015). This is because N-PDFs can constrain the fraction of dense gas within molecular clouds and provide a means of comparison with analytic models as well as numerical simulations (via synthetic observations) of star formation.
N-PDFs, as available from observations, have been utilized extensively for many different tracers of the ISM. This includes molecular line tracers such as CO (Lee et Burkhart et al. 2013b) and column density tracers such as dust (Kainulainen et al. 2009;Froebrich & Rowles 2010;Schneider et al. 2013Schneider et al. , 2014Schneider et al. , 2015bLombardi et al. 2015). Tracing the N-PDF using dust emission and absorption provides the largest dynamic range of densities, in contrast to molecular line tracers such as CO, which do not trace the true column density distribution due to depletion and opacity effects (Goodman et al. 2009a;Burkhart et al. 2013a,b).
Both the true 3D density (volume density) PDF and N-PDFs have been used to understand the properties of galactic gas dynamics, from the diffuse ionized medium to dense star-forming clouds (Hill et al. 2008;Burkhart et al. 2010;Maier et al. 2017). This is because the shape of the density/column density PDF is expected to be related to the underlying physics of the cloud and linked to the kinematics and the chemistry of the gas (Vazquez-Semadeni 1994;Padoan et al. 1997;Kritsuk et al. 2007;Burkhart et al. 2013bBurkhart et al. , 2015a. The low column density gas in molecular clouds, as well as in the diffuse neutral and ionized ISM, takes on the form of a lognormal (Vazquez-Semadeni 1994;Hill et al. 2008;Kainulainen & Tan 2013). This is primarily attributed to the application of the central limit theorem to a hierarchical (e.g. turbulent) density field generated by a multiplicative process, such as shocks. If the width of the lognormal portion of the N-PDF can be measured, it may be related to the sonic Mach number of the gas in a nearly isothermal cloud (Federrath et al. 2008;Burkhart et al. 2009;Kainulainen & Tan 2013;Burkhart & Lazarian 2012;Burkhart et al. 2015a).
However, there are known caveats to constraining the exact shape of the lognormal N-PDF. Below the sensitivity limit, observations are incomplete and do not constitute a statistically meaningful representation of the column density distribution. The sampled N-PDF is then subject to the uncertainty due to the choice of the map area and the foreground/background contamination (see Figure 1; Lombardi, Alves, & Lada 2015). Lombardi, Alves, & Lada (2015) found that both changing the map area used to derive the PDF and subtracting the foreground/background uniformly over the map affect the PDF shape at the low column density end (around A K ∼ 0.1, or A V ∼ 1, for nearby molecular clouds). Lombardi, Alves, & Lada (2015) further suggests that the sampling is statistically unbiased only when the region used to sample the PDF is defined by a "closed contour" (Alves et al. 2017).
In evolved molecular clouds with active star formation, the PDF shape of the densest gas/dust has been observed to develop a power-law form (Kainulainen & Tan 2013;Schneider et al. 2015a;Lombardi et al. 2015). This is expected from both numerical simulations and the analytic theory which suggest that the PDF of self-gravitating gas should develop a power-law tail (Shu 1977;Federrath & Klessen 2012;Burkhart et al. 2015a;Meisner & Finkbeiner 2015). The slope of the PDF power law tail depends on self-gravity and the magnetic pressure (Kritsuk et al. 2011;Ballesteros-Paredes et al. 2011;Collins et al. 2012;Federrath & Klessen 2013;Burkhart et al. 2015a), and can be analytically related to the powerlaw index of a collapsing isothermal sphere (Shu 1977;Girichidis et al. 2014).
Even in the most evolved molecular clouds which exhibit a power law tail towards the dense gas, a lognormal portion of the PDF can still be observed at low total gas and dust column densities. Recently, the HI PDF in and around GMCs has been measured and shown to be the primary component of the lognormal portion of the total gas+dust PDF (Burkhart et al. 2015c;Imara & Burkhart 2016). Burkhart et al. (2015c); Imara & Burkhart (2016) have shown that the lognormal portion of the column density PDF in a sample of Milky Way GMCs is comprised of mostly atomic HI gas while the power-law tail is built up by the molecular H 2 , with no contribution from HI. These studies, including an analytic study by Burkhart, Stalpes, & Collins (2017), suggest that the transition point in the column density PDF between the lognormal and power-law portions of the column density PDF traces important physical processes. These include the HI-H 2 transition and the so-called "postshock density" regime where the background turbulent pressure equals the thermal pressure and therefore selfgravity becomes dynamically important in the molecular gas (Burkhart, Stalpes, & Collins 2017;Li, McKee, & Klein 2015;Kritsuk, Norman, & Wagner 2011). We illustrate the idealized PDF and the known physics associated with different components in Figure 1.

Understanding Anatomy
What sets the shape of the N-PDF of a star forming molecular cloud? In this paper we seek to address this question by comparing dendrogram-based N-PDFs of observations with simulations that include turbulence, magnetic fields, and self-gravity.
Dendrograms are hierarchical tree-diagrams composed of branches, which are split into multiple substructures, and leaves, which have no measurable substructure (Rosolowsky et al. 2008). Dendrograms have been applied to star forming turbulent clouds in the past towards understanding the hierarchical properties of star forming regions (Rosolowsky et al. 2008;Goodman et al. 2009b;Beaumont et al. 2012) and supersonic turbulence (Burkhart et al. 2013a). We can use dendrograms to break up the column density PDF into its hierarchical constituent parts and relate these individual components to the underlying physics.
In this work, we use the column density PDF derived from dust tracers based on the Herschel observations of the dust thermal emission (the Gould Belt Survey, André et al. 2010), as well as the column density PDF from simulations that include turbulence and self-gravity (Collins et al. 2012;Burkhart et al. 2015a). The paper is organized as follows: In §2 we describe the data we use in this paper, namely Herschel observations of L1689 in Ophiuchus and Enzo MHD simulations, as well as the methods used in our analysis of the N-PDF, including the dendrogram algorithm and the fitting to the lognormal and the power-law models. We present our analysis of the lognormal component in §3 for the L1689 region in Ophiuchus and for the Enzo simulations. We then present our dendrogram lognormal + power-law PDF analysis for both simulations and observations in §4. We consider the relevance of the analytic solution to the column density at the transition point, as proposed by Burkhart, Stalpes, & Collins (2017), in §5.1 and in §5.2, respectively for observations and simulations. Finally, we discuss our results in §6 followed by our conclusions in §7.

Observation
The observed column density for the L1689 region in Ophiuchus (Oph L1689) is derived from data taken by the Herschel Space Observatory. Herschel was a satellite operated by the European Space Agency. Its Photodetecting Array Camera and Spectrometer (PACS) and Spectral and Photometric Imaging Receiver (SPIRE) covered wavelengths from 55 to 670 µm, with six broad spectral bands identified by their nominal central wavelengths: PACS 70 µm, 100 µm, and 160 µm, and SPIRE 250 µm, 350 µm, and 500 µm. The data used to derive the column density presented in this paper were obtained as part of the Herschel Gould Belt Survey (André et al. 2010).
To derive the column density, we make use of the maps at 160 µm, 250 µm, 350 µm, and 500 µm, produced by the Herschel Interactive Processing Environment (HIPE; Version 11.1.0). The maps produced by HIPE were not absolutely calibrated, in the sense that there may be a per-band additive offset needed to correct the Herschel zero level. This issue has traditionally been addressed for small 1 square degree regions by adding a scalar offset to the Herschel mosaic under consideration at each wavelength. However, in this study we seek to calibrate Herschel data for a large region of the Ophiuchus cloud dozens of square degrees in size. Over this sizeable footprint, we found that a single per-band scalar offset could not satisfactorily correct the Herschel zero level.
Instead, we chose to allow for a spatially varying zero- ) and sampled by a histogram. The red curve shows a perfect lognormal + power-law distribution, denoting a fitted model. The light gray areas at the lowest column density end indicate the uncertainty in the N-PDF below the sensitivity limit. In the intermediate column density range, the light gray areas indicate the uncertainty in the N-PDF due to the choice of the map area. Toward the high column density end, the colored power-law distributions in the power-law regime of the schematic N-PDF are the N-PDFs expected from self-gravitating substructures. The summation of N-PDFs of all self-gravitating substructures is shown as the gray area in the power-law regime, making up the power-law component of the full N-PDF. Lastly, the bluefilled lognormal curve denotes the purely lognormal distribution observable with atomic gas tracers such as HI 21cm emission (Burkhart et al. 2015c). level offset in each Herschel band. Specifically, we used the Meisner & Finkbeiner (2015) Planck -based thermal dust emission model to predict the Herschel emission at 10 FWHM over our entire Ophiuchus footprint. These low-resolution predictions incorporated color corrections to account for the Herschel bandpasses. To correct the Herschel zero level, we then high-pass filtered each Herschel mosaic at 10 , and replaced the low-order spatial modes (≥10 FWHM) with the corresponding Planckbased predictions. We thereby achieved dust emission maps at 160µm, 250µm, 350µm, 500µm which retain the high angular resolution of Herschel, but inherit the reliable zero level of Planck.
To derive reliable column density maps, we assume that the Herschel maps from 160µm to 500µm, following the zero level corrections, are dominated by "big grain (BG; Stepnik et al. 2003)" thermal dust emission. We also assume that the thermal dust emission can be described by a single modified blackbody (MBB). This assumption is only valid under certain conditions and becomes inaccurate at wavelengths shorter than 100µm due to contamination from "very small grain" (VSG) emission. Adopting a single-component modified blackbody model also presumes that the material along the line of sight can be characterized well by a single temperature. We incorporate a spatially varying value of dust emissivity power-law index, β in the following SED fitting.
We smooth the Herschel maps at 160µm, 250µm and 350µm to 36.1 FWHM, to match the angular resolution of the SPIRE 500µm map. In the SED fitting, we assume a 10% fractional uncertainty on each Herschel intensity measurement. For each pixel, we then derive the optimal temperature and 350µm intensity via simple χ 2 minimization. Using the equation I ν =τ ν B ν (T ), we can then derive the optical depth at any frequency based on our two fitted parameters. Then, a conversion from the 350µm optical depth to column density units are obtained by convolving the map of the 350µm optical depth (γ 350 ) to the nominal resolution of the NICEST map (Lombardi 2009) and comparing to the NICEST extinction map (in the unit of K-band extinction magnitude, A K ). A simple power law, A K = γτ c 350 , is fitted for data points around the median value of 350µm optical depth, ∼ 2.5×10 −4 . The resulting parameters, γ ∼ 2520 and c ∼ 1.11, are consistent with the solution suggested by Lombardi et al. (γ = 2500 with an almost perfect linearity, see also Eq.11-14;. Figure 2 is the final column density map of the L1689 region. The map covers a region of ∼ 2.5 pc by 2.5 pc, and includes all material with column density larger than A K = 0.8 mag (used to define the "dense molecular cloud" where the star formation almost exclusively occurs Lada 1992;Lada et al. 2010).
The resulting column density measurements can then be converted to the unit of equivalent V-band extinction magnitudes, using A V = A K /0.112 (Rieke & Lebofsky 1985), and to the number column density, using N(H 2 )/A V = 9.4×10 20 cm −2 (Bohlin et al. 1978) or N(H 2 )/A V = 6.9×10 20 cm −2 (Draine 2003;Evans et al. 2009). In this paper, for easier comparison with other works in various column density units and between simulation and observation, the column density is expressed in a dimensionless, normalized unit where the column density is divided by the median value above the detection limit. For L1689, the median column density is ∼ 7.3 mag in the unit of equivalent V-band extinction (A V ). Figure 2 shows the map used in the following analysis, and Figure 5 shows the N-PDF of the entire L1689 region in the normalized units. Note that the conversion from physical to dimensionless units does not change the shape of the N-PDF on the logarithmic scale.
An overview of physical properties of Oph L1689 is listed in Table 1, in comparison to the Enzo simulation.

Turbulence in Oph L1689
To estimate the magnitude of turbulent motions in Oph L1689, we fit Gaussian profiles for spectra from the FCRAO observations of 13 CO (1-0) molecular line emission (using data from the COMPLETE Survey, Ridge et al. 2006). Assuming that the dust temperature derived from Herschel observations of thermal dust emission is representative of the gas temperature (i.e. the gas and the dust are in thermal equilibrium; see §2.1 for details on the fitting of dust properties), we calculate the sonic Mach number for Oph L1689, M s = 4.50 +1.43 −1.09 . Since the 13 CO (1-0) line emission traces a density range similar to the density range (as traced by thermal dust emission) in question in this paper, we use the Gaussian line widths of the 13 CO (1-0) transition and the estimated sonic Mach number to assess the dynamics of Oph L1689 in the following analyses (see §4.1). The estimated sonic Mach number is also used to calculate the transitional column density ( §5.1), between the lognormal and the power-law components, in the analytic model proposed by Burkhart, Stalpes, & Collins (2017). See §2.2.2 and Table 1 for a comparison of physical properties to the Enzo simulation used in this paper.

Simulation
In order to investigate the physics of the N-PDF without the uncertainties that plague observations, we use MHD simulations of a collapsing turbulent cloud produced by the Enzo code (Collins et al. 2010(Collins et al. , 2012. The Enzo simulations used in this paper are generated by solving the ideal MHD equations with large-scale solenoidal forcing (b ≈ 1/3, where b is the forcing parameters; see §5 in this paper, and Collins et al. 2012;Burkhart et al. 2017). The simulations have a sonic Mach number of M s = 9 and an Alfvénic Mach number of M A = 12, which scales to ∼ 4.4 µG assuming a sound speed of c s = 0.2 km s −1 and an average volume density of n H = 1000 cm −3 (the "mid" case in Burkhart, Collins, & Lazarian 2015a). To investigate properties of the column density profile under the influence of gravitational collapse, the snapshot at the simulation time t = 0.6 t free-fall is taken. Each side of the simulation cube is ∼ 4.6 pc in length, and the density and the velocities are sampled on a 512 3 grid (resulting in a coarser grid cell size of ∼ 9×10 −3 pc, or ∼ 1.8×10 3 AU; compared to the smallest physical scale the simulation resolves with adaptive mesh refinement of ∼ 500 AU).
To obtain the column density, the density cube is integrated along one of the three axes of a density cube, through the full 4.6 pc. The resulting column density map is then convolved with a Gaussian beam with the same size as that of the Herschel 500-µm beam. For  Figure 7.) The gray contour corresponds to N/N 0 = 0.95, which is used as the minimum value (min_value in astrodendro; see §2.3) in the dendrogram analysis. N 0 is ∼ 7.3 mag in the unit of equivalent V-band extinction (A V ) for L1689. subsequent analyses involving the dendrogram algorithm, we take a 2D map region of 2.3pc by 2.3pc, in order to compare to Oph L1689. Similarly, the column density is expressed in the normalized units, where the column density is divided by the median value. (See §2.1 for details about the "normalization".) The median value of the simulated column density map is ∼ 1.3×10 22 cm −2 , or ∼ 13.5 mag in the unit of equivalent V-band extinction (A V ). See Collins et al. (2012) and Burkhart, Collins, & Lazarian (2015a) for details on scalings to physical units.
An overview of physical properties of the Enzo simulation is listed in Table 1, in comparison to Oph L1689.

Turbulence in the Enzo simulation
To estimate the magnitude of turbulence in the 2.3pc by 2.3pc by 4.6pc cube (a 2.3pc by 2.3pc map with a 4.6pc line of sight) from which the column density map ( Figure 3) was derived, we calculate the sonic Mach number, M s , based on the 3D velocity dispersion in the 2.3pc by 2.3pc by 4.6pc cube. Since we chose a region where gravity is particularly dominant (see Figure 3, in which multiple high-density structures are identifiable), we expect a smaller sonic Mach number (less turbulent material) than expected of the entire 4.6pc by 4.6pc by 4.6pc cube (M s = 9; Collins et al. 2012). For the 2.3pc by 2.3pc column density map, we find M s = 8.14 +2.95 −1.81 . The estimated sonic Mach number is used to calculate the transitional column density in the analytic model proposed by Burkhart, Stalpes, & Collins (2017). See §5 for the analysis of the analytic model. Table 1 shows that Oph L1689 and the Enzo simulation used in this paper have different median column densities and sonic Mach numbers, indicating that Oph L1689 and the Enzo simulation are at different stages of star formation and/or under the effects of different levels of gravity, turbulence, and likely also magnetic field. But, since the main goal of this paper is to examine the uncertainties in the N-PDF analysis and the origins of the different components of the N-PDF, we are more concerned with the relative column density structures than the absolute dynamics of Oph L1689 and the Enzo simulation. The dendrogram analyses of Oph L1689 and the Enzo simulation show that the two do have similar hierarchical density structures (see Table 1, and compare Figure 7 to Figure 8). As a side note, Beaumont et al. (2013) have also demonstrated that it remains difficult to make simulations "look like" a real molecular cloud, even in terms of the simplest observable diagnostics including the column density distribution and the distribution of the CO line widths.

Dendrogram
We use dendrograms (Rosolowsky et al. 2008) to identify substructures inside molecular clouds. The dendrogram is a structure finding algorithm that decomposes an N-dimensional (N ≥ 1) entity into smaller substructures based on the positions and the values of the pixels. The dendrogram algorithm can build a tree of substructures representative of the hierarchy of nested column density features inside a cloud (for previous examples of application in astrophysics and a detailed description of the algorithm, see Rosolowsky et al. 2008;Goodman et al. 2009b;Burkhart et al. 2013a). This hierarchical representation is ideal for the following analysis, where we hope to find the composition of various components of the N-PDF. Figure 4 shows a cartoon demonstrating how a dendrogram is computed from a two-dimensional map. The dendrogram (shown on the left-hand side of Figure 4) is a tree diagram where each structure can harbor exactly two substructures. The structures can then be categorized into leaf, branch, or trunk structures. The trunk structure is a special case, as it is the bottommost structure in the dendrogram.
We define the height of any structure (leaf, branch, or trunk) to be its vertical extent as shown in Figure 4. Three input parameters effect the outcome of the dendrogram algorithm. First, the minimum value in a trunk structure has to be larger than min_value 3 (the light blue dashed line in the dendrogram on the left and the light blue dashed contour in the 2D map on the top right of Figure 4). Second, the height of any substructure in the dendrogram has to be larger than min_delta (the left inset diagram on the bottom right of Figure 4). Lastly, the area of a leaf structure must be larger than min_npix (the right inset diagram on the bottom right of Figure  4).
In this paper, we use the Python-based astrodendro package 3 to compute and analyze the dendrograms of the column density maps. (See §2.1 and §2.2 for details on how we derive the column density maps from observations and simulations.) The astrodendro package offers a simple control of the three essential parameters-min_value, min_delta, and min_npix-in the dendrogram analysis, each defining one of the three criteria described in the above paragraph (see Figure 4). We use the same set of parameters in normalized units for observations and simulations in order to make sure that the analyses are consistent. We choose a minimum value (min_value) of N/N 0 = 0.95, which corresponds to A K ∼ 0.8 mag in observations of Oph L1689. Note that A K = 0.8 mag is also found to define the "dense cloud" where the star formation activity occurs (Lada, Lombardi, & Alves 2010). We then choose min_npix that corresponds to an area equivalent to a 0.05pc by 0.05pc square. And to guarantee a representative sampling of the N-PDF, we also make sure that each substructure has more than 200 pixels (for discussions on the sampling, see Clauset, Shalizi, & Newman 2009). Lastly, min_delta is selected to be N/N 0 = 0.475, roughly corresponding to the 3-σ level in observations. The final dendrograms are shown as connections between panels in Figure 7 and Figure 8. 2.4. Fitting the N-PDF In this paper, the term "N-PDF" is used to indicate the frequency distribution of the column density in a map or a binned histogram representation of this frequency distribution. This definition is mathematically different from a well-defined probability distribution, where the 8 10 a. The dendrograms are computed in normalized units, N/N 0 . b. The Mach number is derived from the average Gaussian line widths for the 13 CO (1-0) molecular line emissions in Oph L1689. c. The Mach number is calculated for the 2.3pc by 2.3pc by 4.6pc cube, from which the 2.3pc by 2.3pc column density map shown in Figure 3 is derived.
However, this discrepancy in definition does not affect the validity of results presented in this paper. In practice, since one can only sample the mathematical probability distribution function with a finite number of independent pixels (in both observations and simulations), the results in this paper can be directly compared to other work where the distributions of values in column density maps are examined.
Following results presented by Schneider et al. (2013); Myers (2015) and Burkhart, Stalpes, & Collins (2017) and the anatomical diagram shown in Figure 1, we assume that the N-PDF has a lognormal component in the low column density regime and a power-law component toward the high column density end. The two components are continuous at the transitional column density (Myers 2015;Burkhart, Stalpes, & Collins 2017). By defining the normalized column density on the logarithmic scale: the distribution can be written as where s t is the transitional column density in the logarithmic normalized units, p 0 is the amplitude of the N-PDF at the transition point, and M is the normalization/scaling parameter. To carry out the least-squares fitting, the N-PDF is first sampled using a binned histogram. The χ 2 -residual between the model and the histogram is minimized over a range of s t , with a step size in s t equal to one third of the histogram bin size. For convenience, the modeled N-PDF based on this equation is called a "lognormal + power-law" model in discussions below. We note that the least-squares fitting cannot be used to determine whether a lognormal + power-law model is a better fit to the observed distribution than a simple lognormal model (Clauset, Shalizi, & Newman 2009) and is subject to various fitting/binning choices. Thus, we limit our analyses involving χ 2 -fitting to the N-PDFs of the cloud-scale regions (the entire L1689 region or the integrated column density maps derived from the simulation cubes; see §2.1 and §2.2 for how we define the regions in observations and simulations). At the cloud scale, the large number of independently sampled pixels (on the order of 10 4 to 10 5 ) makes the χ 2 -fitting results less dependent on fitting/binning choices (Stutz & Kainulainen 2015), and the existence of a power-law component is evident by eye (see Figure 5 and Figure 6).
For smaller substructures within the cloud, we base our analysis on comparing (without fitting) the individual N-PDFs to the N-PDFs of entire regions. Thus, the results presented in this paper are essentially independent of the uncertainties in fitting. (See Stutz & Kainulainen (2015) for a comparison of various fitting/binning schemes, and Burkhart, Collins, & Lazarian (2015a), for alternative N-PDF diagnostics that do not require fitting.)

Observation
We investigate the lognormal component of the N-PDF of Oph L1689 and plot the total PDF in Figure 5.  Figure 5 also demonstrates that sampling of the lognormal component of the N-PDF is subject to the uncertainty due to the map area between N/N 0 ∼ 0.1 and ∼ 0.6 (A V ∼ 0.7 mag and ∼ 4.4 mag, consistent with Lombardi, Alves, & Lada 2015), and that the lognormal component is unobservable below the detection limit (N/N 0 ∼ 0.1 in this case). The uncertainty due to the map area and the detection limit leaves us with a very narrow range of well-sampled column density that can be used in an attempt to fit for the lognormal component. In the case of Oph L1689, the range of column density that is not affected by the uncertainties is 50% of the full range the fitted lognormal component spans (which includes the pixels outside the last closed contour; Alves et al. 2017). Thus, a direct comparison of the width of the lognormal component to the cloud dynamics such as the sonic Mach number is difficult and unreliable using the column density maps based on dust tracers. Figure 6 shows that, at t = 0.6 t free-fall , the 2.3pc by 2.3pc integrated column density map derived from the Enzo simulation has a lognormal component toward the low column density. The shape of the N-PDF is consistent with the N-PDF of the full 4.6pc by 4.6pc column density map derived from the entire cube (Burkhart et al. 2015a;Collins et al. 2012). Figure 6 also shows that, when we mimic the observational constraint on the map area by changing map size, the change in the N-PDF occurs at the low column density side of the peak column density (consistent with observations presented in §3.1 and Figure 5 in this paper and in Lombardi et al. 2015), even though in simulations a complete sampling is possible by including the entire cube. There is no detection limit in simulations. Fitting of the lognormal component is thus stable and can be shown to correlate with physical properties (Burkhart, Collins, & Lazarian 2015a).  Figure 5 shows that the N-PDF of Oph L1689 has a power-law component above N/N 0 ∼ 1.4 (the transition point). In an attempt to understand the composition of the power-law component, we apply the dendrogram algorithm to find structures with physical sizes larger than 0.05 pc above a minimum value of N/N 0 = 0.95 (A K ∼ 0.8, consistent with the "dense cloud" definition suggested by Lada, Lombardi, & Alves 2010), with min_delta = 0.475 in the normalized unit. (See §2.3 for the specifics of the dendrogram algorithm and the setup parameters.) The resulting Ophiuchus dendrogram is shown as the connections between panels in Figure 7.

Observation
When we examine the N-PDFs of individual "leaf" structures (color coded according to their levels in the dendrogram in Figure 7), we find that the N-PDFs of the individual "leaf" structures can be roughly categorized into three categories. First, there are "leaf" structures such as Structure 9, 12, and 18 (marked with "LN" in Figure 7), which sit at the lower levels in the dendrogram and have N-PDFs mostly in the lognormal regime of the entire cloud. These are likely transient (unbound) overdensities, and their N-PDFs do not always look like a lognormal function because that the sample size is small. Secondly, Structure 14 (marked with "LN/PL" in Figure  7) could be at an early stage of gravitational collapse. It sits at a low level in the dendrogram tree, and has an N-PDF consisting of both a lognormal component and a power-law component. Lastly, there are several structures sitting toward the top of the dendrogram tree with N-PDFs almost entirely in the power-law regime of the entire Oph L1689. Structure 6, 8, 10, 11, 16 and 17 belong to the last category (marked with "PL" in Figure  7). The N-PDF of each of these structures has a shape roughly resembling the power-law distribution. Some of these structures (e.g., Structure 6, 8, and 10 in Figure  7) span a smaller range of column density, while others (e.g., Structure 11 and 17 in Figure 7) are much denser with narrow 13 CO (1-0) line widths suggesting that they are gravitationally bound (see a full virial analysis of the dynamics to be presented in Chen et al. 2017, in preparation). Notice that the slopes of the N-PDFs of the dense substructures are not necessarily the same as that of the entire Oph L1689. The different slopes may indicate that the structures are undergoing gravitational collapse at different evolutionary stages (Stutz & Kainulainen 2015;Burkhart et al. 2015c). The sum of all the leaf N-PDFs, given by the grey shaded area in Figure 5 (and in Figure 1 and Figure 6), is a very good approximation to the observed PDF.
Following the nested structures of the dendrogram from the top-level "leaf" structures down to the "branch" structures, we find that the N-PDF of a "branch" structure containing several or more "leaf" structures with power-law N-PDFs has a power-law component similar to the power-law component of the N-PDF of the entire region. (For example, see the N-PDF of the branch that contains Structure 12, 11, 6, 8, and 10 in Figure 7.) And, Figure 5 shows that most of the pixels in the powerlaw component of the entire region are in the independent "leaf" structures in the dendrogram. Since a power-law N-PDF is analytically expected from the self-similar gravitational collapse of a cloud (Shu 1977), these results suggest that the power-law component of an extended region is the summation of N-PDFs of dense, probably self-gravitating, substructures within the region.

Simulation
Similar to observations (Figure 5), Figure 6 shows that the N-PDF of the Enzo simulation at t = 0.6 t free-fall has a power-law component above N/N 0 ∼ 1.2. We then apply the dendrogram algorithm on the simulated column density map with the same set of setup parameters as that applied on observations. (See §2.3 for details on the dendrogram setup parameters.) The dendrogram is shown as connections between panels in Figure 8. Notice that the dendrogram of the simulated column density map has a similar complexity as that of Oph L1689 (Figure 7 and Figure 8). Figure 8 shows the N-PDFs of substructures in the dendrogram. The N-PDFs of all "leaf" structures seem to sit in the power-law regime of the lognormal + power-law model fitted to the N-PDF of the entire region (Equation 2; see §2.4 for details on the lognormal + power-law model). However, we can still identify "leaf" structures that have N-PDFs with shapes of a lognormal distribution (Structure 2), a lognormal + power-law distribution (Structure 3), or a power-law distribution (Structure 6,12,13,15,16,17 and 18). The slopes of the power-law component of the denser substructures are different from each other and from the slope of the power-law component of the entire region. And similar to observations, Figure 6 shows that most of the pixels in the powerlaw component of the entire region are in the independent "leaf" structures. Most of these "leaf" structures also have power-law N-PDFs, albeit with different slopes. And, similar to Oph L1689, the power-law component of the N-PDF of the selected region in the Enzo simulation is a summation of the N-PDFs of likely self-gravitating substructures within the cloud (Shu 1977). Burkhart, Stalpes, & Collins (2017) interpret the transition point as the density at which the turbulent energy density is equal to the thermal pressure and provide an analytic model of the transitional column density value (see also Myers 2015). In this section, we test this analytic model and investigate the physics behind the transitional column density value between lognormal and power-law distributions. In the following paragraphs, we present the analytic model derived by Burkhart, Stalpes, & Collins (2017) assuming an equilibrium between the turbulent energy density and the thermal pressure at the transitional column density. See Myers (2015) for a more general mathematical description of the transition point.

THE TRANSITION POINT
We consider a piecewise form of the N-PDF as in Equation 2. (See §2.4 for more details on the lognormal + power-law model.) By assuming that the N-PDF at the transition point is continuous and differentiable, the transition point was derived in Burkhart, Stalpes, & Collins (2017) as: where s t is the logarithmic normalized column density at Fig. 4.-This cartoon shows how the three criteria used in the computation of a dendrogram, which are described in §2.3, are defined. On the left-hand side of the cartoon, a dendrogram is shown using a vertical axis of column density. The branch, leaf, and trunk structures and their heights are labeled according to their definitions (see §2.3). On the top right is a cartoon 2D map from which the mock dendrogram on the left-hand side is derived. Contours on the 2D map define the boundaries of dendrogram substructures of corresponding colors. The three criteria-1) the minimum value of a trunk structure has to be larger than min_value, 2) the height of any structure in the dendrogram has to be larger than min_delta, and 3) the area of a leaf structure has to be larger than min_npix-are shown in light blue with bold-faced annotations. See §2.3 for details on the definitions. the transition point, α is the slope of the power law tail, and σ s is the width of the lognormal component of the N-PDF.
In the strong collapse limit, where a well defined powerlaw tail is formed, the slope of the power-law component, |α|, assumes a value of ∼ 1.5. In this limit, Burkhart, Stalpes, & Collins (2017) showed that the normalized column density at the transition point, N t /N 0 , can be expressed as a function of the sonic Mach number and the forcing parameter: where b is the forcing parameter, varying from ≈ 1/3 (purely solenoidal forcing) to 1 (purely compressive forcing), and A = 0.11 is the scaling constant from volume density to column density (Federrath et al. 2008;Burkhart et al. 2017). We can then test Equation 4 by comparing this modeled transitional column density with the fitted transitional column density. (See §2.4 for details on how we fit for the transition point.) In the following subsections, we present the results of the tests in observations and simulations.

Observation
We compare the fitted transitional column density value to the modeled value (see Table 2 in this paper; Burkhart, Stalpes, & Collins 2017). Equation 4 shows that the normalized transition column density, N t /N 0 , is dependent on the sonic Mach number, M. We fit Gaussian line profiles for spectra from FCRAO observations of 13 CO (1-0) line emission (the COMPLETE Survey, Ridge et al. 2006) and find an average Mach number of Oph L1689, M s = 4.50 +1.43 −1.09 (Table 2; see §2.1.1 for details on the estimation of the sonic Mach number in observation). Using Equation 4, we then derive the modeled transitional column density, (N/N 0 ) trans, model = 1.14 +0.12 −0.08 with the forcing parameter b = 1/3 (purely solenoidal), and (N/N 0 ) trans, model = 1.40 +0.18 −0.15 with b = 1 (purely compressive). Since, in observations, we do not know the forcing parameter, b, which varies between 1/3 and 1, the modeled transitional column density, (N/N 0 ) trans, model , has a large uncertainty and ranges between 1.14 and 1.40.
To obtain the fitted transitional column density from the N-PDF, we first derive a binned histogram representation of the N-PDF and fit the histogram to the lognormal + power-law model. (See §2.4 for details on the χ 2 -fitting.) For the L1689 data, we find a fitted transitional column density, (N/N 0 ) trans, fit = 1.43±0.14 (Table 2). Compared to the modeled values presented in the above paragraph, we can say, with some uncertainties, that the observed N-PDF has a transition point mildly more consistent with the analytic model with a purely compressive forcing (b = 1). Burkhart, Stalpes, & Collins (2017) have verified that the transitional column density can be described by the above analytic expression ( §5), using a set of Enzo simulations with various setups of the sonic Mach number and the Alfvénic Mach number (see Figure 3 in Burkhart, Stalpes, & Collins 2017). Since it is possible to completely sample the lognormal + power-law distribution in simulations, Burkhart, Stalpes, & Collins (2017) followed Equation 3 (that is Equation 6 in Burkhart, Stalpes, & Collins 2017) and demonstrated that the fitted transitional column density matches the modeled value as a function of the fitted width of the lognormal component, σ s , and the fitted slope of the power-law component, α (Equation 3).

Simulation
In this section, we present a separate verification of the analytic model of the transition point, using the 2.3pc by 2.3pc integrated column density map (Figure 3; see also §2.2 for details on how the map was made). Instead of modeling the transitional column density value based on the fitted width of the lognormal component and of the fitted slope of the power-law component (see Equation  3; Burkhart, Stalpes, & Collins 2017), we follow Equation 4 and model the transitional column density based on the Mach number, M s , and the forcing parameter of the simulations, b. Since the Mach number and the forcing parameter are independent from the fitting of the N-PDF, the test presented below following Equation 4 is stricter than (and independent from) the one presented in Burkhart, Stalpes, & Collins (2017), which follows Equation 3 and involves fitting the N-PDF on both sides of the equation.
We calculate the sonic Mach number, M s , based on the 3D velocity dispersion of the 2.3pc by 2.3pc by 4.6pc cube (a 2.3pc by 2.3pc map with a 4.6pc line of sight) from which the column density map was derived. Since the gravity in the 2.3pc by 2.3pc region is particularly dominant (see Figure 3 in which multiple highdensity structures are identifiable), we expect a smaller sonic Mach number (less turbulent materials) than expected of the entire cube (M s = 9; Collins et al. 2012). For the 2.3pc by 2.3pc region shown in Figure 3, we find M s = 8.14 +2.95 −1.81 (see §2.2.1 for details on the estimation of the sonic Mach number in simulation, and also §2.2.2 for a discussion on differences in the physical properties between Oph L1689 and the Enzo simulation used in this paper). Knowing that the simulation has purely solenoidal forcing (b ≈ 1/3), we then derive the transitional column density based on the analytic model (Equation 4), (N/N 0 ) trans, model = 1.26 +0.11 −0.08 . This is consistent with the fitted transitional column density, (N/N 0 ) trans, fit = 1.18±0.12 ( Table 2). The result again verifies that the transition point of the observed N-PDF is well described by the analytic model proposed by Burkhart, Stalpes, & Collins (2017), especially in simulations. Table 2 gives an overview of comparisons between the fitted and modeled transitional column densities in observation and simulation. We see that when the transition point between the lognormal and the power-law components of an N-PDF can be fitted, the analytic model could be potentially useful for deriving the dynamics of star-forming materials from the column density distribution. Unfortunately, the ability to estimate the transitional column density in the lognormal + powerlaw model is limited by the uncertainty due to changing the map area used to derive the N-PDF in observations. On top of the difficulty in fitting, the forcing parameter, b, is usually difficult to measure in observations, adding 1.18±0.12 Yes a. The forcing parameter ranges from 1/3 (purely solenoidal forcing) to 1 (purely compressive forcing). b. The uncertainty is estimated based on the inherent uncertainties of Herschel observations. c. The two values are said to agree when the ranges enclosed by the uncertainties overlap. d. The Mach number is derived from the average velocity dispersion of Gaussian fits to the 13 CO (1-0) molecular line emissions in the region.
a large uncertainty to the modeled transitional column density (see recent attempts by Orkisz et al. 2017;Herron et al. 2017;Otto et al. 2017).

DISCUSSION
Our study highlights the importance of comparing observations and simulations. For example, the lognormal portion of the PDF in observations suffers biases, such as boundary effects and unresolved foreground/background contributions (Schneider et al. 2015a;Lombardi et al. 2015). Simulations provide an avenue to study these effects as the density in the simulations is completely sampled due to mass conservation and periodic boundary conditions. However, the simulations are missing important physical effects such as feedback from stars (Offner & Arce 2015), non-isothermal effects (Nolan et al. 2015) and non-ideal MHD effects (Meyer et al. 2014;Burkhart et al. 2015b). No single simulation is able to capture all the physical processes and scales involved in star formation. Simulations, therefore, can serve only as a general reference to interpret the observations.
One aspect not addressed in our study of the anatomy of the N-PDF is the influence of the magnetic field. The effect of the magnetic field on the shape and behavior of the PDF of gravitating turbulent clouds has been studied in the past (Burkhart & Lazarian 2012;Mocz et al. 2017;Kritsuk et al. 2011;Collins et al. 2012;Federrath & Klessen 2012). In general, the lognormal portion of the N-PDF is not strongly affected by changing the magnetic field strength Burkhart & Lazarian 2012;Collins et al. 2012;Federrath & Klessen 2012;Burkhart et al. 2015a). However, studies have shown that the strength of the magnetic field can alter the slope of the power-law tail portion of the N-PDF (Burkhart et al. 2015a;Mocz et al. 2017). These studies found that a higher magnetic field strength produces an N-PDF with a steeper power-law tail slope since the magnetic field inhibits the collapse.
6.1. The dendrogram analysis and the N-PDF of an independent structure The analysis presented in this paper shows that the power-law component of the N-PDF is the sum of individual substructures with power-law PDFs. Past a transition point where the shape of the PDF changes from lognormal to power-law, individual dendrogram structures show clear power-law forms, while at or below the transition point the column density PDFs of substructures primarily take on a lognormal form (e.g. structures 11, 14, 17 in Figure 7).
One notable deviation from the above picture is Structure 6 (see Figure 2 and Figure 7) in the dendrogram analysis of the N-PDF of Oph L1689. Structure 6 sits above the transition point, so we would expect that its N-PDF should take on a power-law shape. However, the N-PDF of Structure 6 is visually indistinguishable from a lognormal distribution. Despite the fact that there are only ∼ 2×10 2 pixels in Structure 6 and that the sampling of the N-PDF is thus probably incomplete, we want to ask the question: Does the lognormal shape of the N-PDF of Structure 6 mean that the substructure is, in fact, not gravitationally bound? Ongoing star formation is evident in this region within a time scale ∼ 5×10 5 years, since there are multiple young stellar objects nearby (Gutermuth et al. 2009). Feedback from active star formation could have shredded the star forming materials in the region into smaller, unbound pieces like Structure 6. The shredded pieces would have N-PDFs that become increasingly steep, potentially mimicking a power-law distribution at high column density. (Similar effects of feedback on the N-PDF shape have been observed at cloud scales by Schneider et al. (2015a).) We also note that the Enzo simulation we compare to (Figure 8) lacks a similar lognormal substructure at these high densities, and that there is no feedback in the Enzo simulation. Local environmental effects such as feedback therefore most likely play a role in the unusual shape of structure 6.
If an independent structure (a leaf structure in the dendrogram) is axially symmetric (e.g., a filamentary structure) or radially symmetric (e.g., a spherical structure), Myers (2015) pointed out that the N-PDF is the inverse function of the column density profile along a cut perpendicular to the axis of symmetry or in the radial direction, respectively. Myers (2015) observed that, for a filamentary structure with a Plummer-like profile perpendicular to the axis of the filament (typical of star forming filamentary structures; Peretto et al. 2011), the inverse function of the Plummer-like column density profile deviates from the lognormal+power-law model of the N-PDF. Myers (2015) suggested that this is likely due to the nearly constant column density within the first thermal scale length. Similarly, the inverse function of the narrow N-PDF of structure 6 in Figure 7 would have a flat shape, consistent with a nearly constant column density profile. Is the thermal length scale then the scale limit for the N-PDF to be representative of the gravitational effect? A more detailed analysis of the velocity structure of dendrogram features such as Structure 6 in Oph L1689 (see Figure 7) is needed to answer this question.  Figure 1, but using the real Herschel-based column density map of Oph L1689 instead. The solid black histogram represents the N-PDF of the Hershcle-based map. The red curve is a lognormal + power-law fit to the histogram above the range of column density affected by the uncertainty due to the choice of the map area and the detection limit. The light gray area toward the low column density end shows the change in N-PDF following the change in the map area by ±20%. Toward the high column density end, the colored power-law distributions in the power-law regime of the full N-PDF are the N-PDFs of independent "leaf" structures in the dendrogram. The summation of N-PDFs of all "leaf" substructures is shown as the gray area in the power-law regime, making up most of the power-law component of the full N-PDF.

CONCLUSION
We present the anatomy of the column density PDF (N-PDF) of star forming molecular clouds in both observations and simulations. By assuming a lognormal + power-law model, we examine how the lognormal component, the power-law component, and the transition point could be useful for estimating the dynamical properties of a star forming region. We also examine the uncertainties that could affect the N-PDF analysis.
With the help of the dendrogram algorithm, we demonstrate that the power-law component of an N-PDF is primarily a summation of N-PDFs of substructures inside the star forming cloud. Most of these substructures show N-PDFs following power-law distributions, with powerlaw indices different from the N-PDF of the entire region. The power-law shapes and varying indices of N-PDFs suggest that these substructures could be going through different stages of gravitational collapse.
The analytic prediction of a transition point between lognormal and power-law components proposed by Burkhart, Stalpes, & Collins (2017) is verified independently (within uncertainties) in this work using data from both observations and simulations. We show that the transition point is generally well described by the analytic model (see Equation 4 in this paper; Burkhart, Stalpes, & Collins 2017). The result suggests that finding the transitional column density value in a N-PDF could potentially provide information on the dynamics of the cloud. In particular, Equation 4 could potentially be useful for getting a general estimate of the sonic Mach number, even though it is usually difficult to determine the transition point and the forcing parameter in observations (see a recent attempt by Orkisz et al. 2017).
Based on the results presented in this paper, we give the following suggestions for future studies involving the N-PDF analysis. First, we do not suggest analyzing the dynamics based solely on fits to the lognormal component of the N-PDF derived from observations of dust tracers. Second, measuring the column density at the transition point between the lognormal and the powerlaw components could potentially provide information on the dynamical properties of a star forming region, but only when the transitional column density and the forcing parameter can be estimated. Lastly but most importantly, we recognize that the power-law component of an N-PDF is a summation of N-PDFs of substructures, which are likely going through various stages of gravitational collapse. Based on the results presented in this paper, we suggest combining the N-PDF analysis with the dendrogram algorithm to obtain a more complete picture of the effects of global and local environments, including gravity, turbulence, the magnetic field, and the feedback from star formation, in future attempts to analyze the column density structures of a star forming region.
B.B. acknowledges support from the NASA Einstein Postdoctoral Fellowship.
Computer time was provided through NSF TRAC allocations TG-AST090110 and TG-MCA07S014. The computations were performed on Nautilus and Kraken at the National Institute for Computational Sciences (http://www.nics. tennessee.edu/).  Figure 1, but using the integrated column density map based on the Enzo simulation (see §2.2). The solid black histogram represents the N-PDF of the entire column density map. The red curve is a lognormal + power-law fit to the histogram above the range of column density affected by the uncertainty due to the choice of the map area. The light gray area around N/N 0 = 1 shows the change in N-PDF following the change in the map area by ±20%, even though it is possible to completely sample the N-PDF by including the entire simulation cube. Toward the high column density end, the colored power-law distributions in the power-law regime of the full N-PDF are the N-PDFs of independent "leaf" structures in the dendrogram. The summation of N-PDFs of all "leaf" substructures is shown as the gray area in the power-law regime, making up most of the power-law component of the full N-PDF. -The connection between dendrogram structures are shown as the thick lines connecting the panels, and the N-PDF of each structure is shown inside each panel. The red curve in each panel is a lognormal + power-law fit to the N-PDF of the entire Oph L1689. The "leaf" structures (i.e. independent structures that cannot be broken down into smaller structures) in the dendrogram are color coded according to their levels in the dendrogram tree. The numberings of the substructures are produced by astrodendro. LN, LN/PL and PL denote the categories discussed in §?? and stand for "lognormal," "lognormal + power-law" and "power-law." See Figure 3 for locations of the "leaf" structures on the integrated column density map.  Figure 7, but for the Enzo simulation. The connections between dendrogram structures are shown as thick lines connecting the panels, and the N-PDF of each structure is shown inside each panel. The red curve in each panel is a lognormal + power-law fit to the N-PDF of the entire 2.3pc by 2.3pc region. The "leaf" structures in the dendrogram are color coded according to their levels in the dendrogram tree. The numberings of the substructures are produced by astrodendro. LN, LN/PL and PL denote the categories discussed in §4.2 and stand for "lognormal," "lognormal + power-law" and "power-law." See Figure 3 for locations of the "leaf" structures on the integrated column density map.