MNRAS 000 ,

1-14 (2018)

## On the extraction of power-law parts of the probability density functions in star-forming clouds

Todor V. Veltchev 1 2 , glyph[star] , Philipp Girichidis , Sava Donkov , Nicola Schneider , Orlin 3 4 5 Stanchev , Lyubov Marinkova , Daniel Seifried 1 1 5 and Ralf S. Klessen 2

- 1 University of Sofia, Faculty of Physics, 5 James Bourchier Blvd., 1164 Sofia, Bulgaria
- 2 Universit¨t Heidelberg, Zentrum f¨r Astronomie, Institut f¨r Theoretische Astrophysik, Albert-Ueberle-Str. 2, 69120 Heidelberg, Gera u u

many

3 Leibniz-Institut f¨r Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany u

4 Department of Applied Physics, Technical University, 8 Kliment Ohridski Blvd., 1000 Sofia, Bulgaria 5

- I. Physik. Institut, University of Cologne, Z¨lpicher Str. 77, 50937 Cologne, Germany u
- 6 OASU/LAB, Universit´ e de Bordeaux, 33615 Pessac, France

Submitted 2019 March 18

## ABSTRACT

We present a new approach to extract the power-law part of a density/column-density probability density function ( ρ -pdf/ N -pdf) in star-forming clouds. It is based on the mathematical method bPlfit of Virkar &amp; Clauset (2014) and assesses the power-law part of an arbitrary distribution, without any assumptions about the other part of this distribution. The slope and deviation point are derived as averaged values as the number of bins is varied. Neither parameter is sensitive to spikes and other local features of the tail. This adapted bPlfit method is applied to two different sets of data from numerical simulations of star-forming clouds at scales 0.5 and 500 pc and displays ρ -pdf and N -pdf evolution in agreement with a number of numerical and theoretical studies. Applied to Herschel data on the regions Aquila and Rosette, the method extracts pronounced power-law tails, consistent with those seen in simulations of evolved clouds.

Key words: ISM: clouds - Physical data and processes: gravity - Physical data and processes: turbulence - - methods: data analysis - methods: statistical

## 1 INTRODUCTION

The analysis of the probability density functions (pdf) of mass density ( ρ -pdf) and of column density ( N -pdf) has been increasingly recognized as a key to understand the physics of star-forming regions. About two decades ago, it was shown from theoretical modelling (V´zquez-Semadeni a 1994; Passot &amp; V´ azquez-Semadeni 1998) that supersonic turbulence in an isothermal and non-gravitating medium produces a lognormal ρ -pdf. This was confirmed by many numerical studies with varying characteristics of the environment like magnetic field, equation of state, Mach number, turbulence driving (e.g., Padoan, Nordlund &amp; Jones 1997; Klessen 2000; Li, Klessen &amp; Mac Low 2003; Kritsuk et al. 2007; Federrath et al. 2010; Konstandin et al. 2012; Molina et al. 2012). When the importance of gravity gradually increases in dense zones of the considered region, the ρ -pdf develops a power-law (PL) tail at its high-density end (Klessen 2000; Slyz et al. 2005; V´ azquez-Semadeni et al. 2008; Kritsuk, Norman &amp; Wagner 2011; Collins et al. 2012; Federrath &amp; Klessen 2013). In the course of further evolution of self-gravitating clouds, the point of deviation from the lognormal shape shifts to lower densities and the slope of the PL tail (hereafter, PLT) flattens towards a constant value (Girichidis et al. 2014).

Could one discern similar signatures in the corresponding column-density distributions, derived from plane projections of the numerical cubes? Recent simulations of formation and contraction of molecular clouds (MCs) clearly hint at an evolution of the N -pdf, morphologically analogous to that of the ρ -pdf. The obtained N -pdfs were found to have a lognormal shape in simulated turbulent clouds (Federrath &amp; Klessen 2013; Ward, Wadsley &amp; Sills 2014; Matsumoto, Dobashi &amp; Shimoikura 2015) while in contracting/starforming clouds their main part is still lognormal, but with a pronounced PLT (Ballesteros-Paredes et al. 2011; Kritsuk, Norman &amp; Wagner 2011; Federrath &amp; Klessen 2013; Auddy, Basu &amp; Kudoh 2018; K¨ ortgen, Federrath &amp; Banerjee 2019). In the latter case, there are asymptotic thresholds for slope and deviation point of the PLT at later evolutionary stages (Ward, Wadsley &amp; Sills 2014). Their values are not universal but depend on the sonic Mach numbers and the type of turbulence forcing (see Fig. 4 in Federrath &amp; Klessen 2013).

This general picture appears to be in good agreement with observations of star-forming clouds whose N -pdfs could be decomposed to a lognormal part and a PLT at their highdensity ends (e.g., Schneider et al. 2013, 2015a; Pokhrel et al. 2016). Flattening of the slope with the cloud evolution seems also to be supported by the work of Abreu-Vicente et al. (2015) who found shallow slopes for H ii regions and much steeper - for MCs at earlier stages of star formation. On the other hand, clear observational examples of a (quasi)lognormal N -pdf for a quiescent cloud are only a few; for instance, in Polaris (Schneider et al. 2013) and in Draco and Spider (Schneider et al. 2019).

In studies of a thermally bistable/multiphase interstellar medium, it has been found that the ρ -pdf and the N -pdf have more complex shapes accounted for by various physical factors. Thermal instabilities in purely hydrodynamic flows may generate bimodal ρ -pdfs while turbulent energy injections at small scales and/or including magnetic fields, Coriolis force and stellar feedback lead to evolutionary transition to a single-peaked distribution, possibly with a PLT (V´ azquez-Semadeni, Gazol &amp; Scalo 2000). The type of turbulent forcing modifies additionally the high-density part of the ρ -pdf: a PLT in the case of compressive forcing and a single lognormal shape, i.e. no bimodality, in the case of purely solenoidal forcing (Seifried, Schmidt &amp; Niemeyer 2011). Ballesteros-Paredes et al. (2011) demonstrated from numerical simulations of a thermally unstable gas that the N -pdf is also bimodal - this is clearly visible at early evolutionary stages and still detectable at later stages, as a highdensity PLT develops when self-gravity becomes important. These two cases are recently testified from Herschel observations of the diffuse cloud Draco (Schneider et al. 2019) and of MCs associated with H ii regions (Tremblin et al. 2014). In general, the emergence of a PLT in the (column-)density distribution is an established phenomenon in evolved thermally bistable/multiphase medium.

The decomposition of N -pdfs in star-forming clouds in to a (quasi-)lognormal part and a PLT is also a challenge from the perspective of observations. It has been currently debated whether a lognormal part can be extracted reliably. Indeed, the very low column-density range is affected by noise biases, uncertainties due to map area and small dynamic range of dust observations while the lognormal fit remains the best model for this part of the N -pdf (OssenkopfOkada et al. 2016; Chen et al. 2018). Lombardi, Alves &amp; Lada (2015) and Alves, Lombardi &amp; Lada (2017) put to question the detectability of a lognormal part of observational N -pdfs at all, attributing seemingly lognormal shapes at low column densities to incompleteness of data and sampling effects ('last closed contour'). However, K¨rtgen, Fedo errath &amp; Banerjee (2019) demonstrate from magnetohydrodynamical simulations, that in a turbulently mixed medium there is no 'last closed contour' for a sufficiently low columndensity threshold.

On the other hand, Brunt (2015) question whether a PL fit is indeed an appropriate description of the high columndensity part. He argues that N -pdfs, derived on the base of dust-extinction maps of MCs, could be reproduced combining two lognormals that account respectively for the warm and the cold gas included the selected region.

But even taking for granted that the N -pdf in a starforming cloud consists of two clearly recognizable parts, there is another, methodological issue. The characteristics of the two fitting functions are obviously interdependent. For instance, the usual procedure to extract the PLT is as follows:

- (i) Find the best lognormal fit of the main (low-density) N -pdf part (e.g. by χ 2 goodness).
- (ii) Estimate the 'deviation point' (DP) of the distribution from the lognormal fit (e.g. by 3 σ criterion, with σ being the Poissonian data uncertainty in the considered bin).
- (iii) Fit the rest of the distribution with a PL function.

Such approach rests on the assumption that the main N -pdf part is lognormal and thus the resulting DP and the PL slope depend on the lognormal parameters. This is problematic not only in view of the abovementioned debates on lognormality of the N -pdf at lower densities. If the PL regime is to be interpreted as a signature of the impact of self-gravity, then the slope value is an indicator of their evolutionary stage. It has been demonstrated that the PL slope is interrelated to the density profile of spherically symmetric, collapsing, isothermal clouds (Kritsuk, Norman &amp; Wagner 2011) and can be representative for the general structure of real MCs (Donkov, Veltchev &amp; Klessen 2017). Therefore one needs a method to extract the PLT on minimal assumptions about the rest of the column-density distribution.

In this paper we propose an objective method to extract the PLT of the ρ -pdf and N -pdf in star-forming MCs. The method is based on the statistical approach to test the power-law features of an arbitrary distribution, developed by Virkar &amp; Clauset (2014), as the PL slope and the DP are estimated simultaneously. In Section 2 we discuss the original approach and comment on its applicability on large datasets. Section 3 reviews the numerical and observational data used in this work. The proposed method itself is described in Section 4. The results from its application to (column-)density distributions from simulations and in two Galactic star-forming regions are presented in Sections 5 and 6, respectively. We conclude this work with discussion (Sect. 7) and summary (Sect. 8).

## 2 MATHEMATICAL METHOD

## 2.1 Application to unbinned data

The original statistical approach 'to discern and quantify power-law behaviour in empirical data' was proposed by Clauset, Shalizi &amp; Newman (2009). It deals with unbinned observational data and thus avoids possible subjectivity of the PL fit introduced through the choice of binning size. Assuming that the data above a chosen lower limit x min follow a power law p x ( ) ∝ x -α , the method derives an estimate of this limit ˆ x min through minimization of the Kolmogorov-Smirnov (KS) goodness-of-fit statistics and of the slope ˆ = ˆ(ˆ α α x min ) through maximum-likelihood estimators (see Clauset, Shalizi &amp; Newman 2009 for details). The procedure does not rule out that other (non-PL) functions might better fit the observed distribution - it derives the range and the slope of the best possible PL fit.

- A software implementation of the method, called

Plfit 1 , has been applied for PLT extraction from N -pdfs in zones of the Perseus star-forming region (Stanchev et al. 2015). Unfortunately, the use of this program is technically limited to smaller datasets.

## 2.2 Application to binned data

To deal with large datasets of size glyph[greaterorsimilar] 10 5 points from numerical simulations and high-resolution imaging of MCs, one could use a version of the approach, described in the previous subsection, adapted for binned distributions (Virkar &amp; Clauset 2014; software implementation is also available ). 2 The technique does not depend on assumptions about the binning scheme and is applicable to linear, logarithmic and arbitrary bins. Hereafter, we call this method bPlfit .

Like its unbinned counterpart, bPlfit uses maximumlikelihood estimators to assess the PLT. The required input is a binned distribution ( { b i } , 1 glyph[lessorequalslant] i glyph[lessorequalslant] k +1; { h i } , 1 glyph[lessorequalslant] i glyph[lessorequalslant] k ) where b i are the bin boundaries and h i are the corresponding counts for k bins. The approach is realized in several basic steps:

- 1. Choose consecutively test values b min , i for the DP, from the set ( b 1 , ..., b k ).
- 2. Calculate the slope α i ( b min , i ) for each test DP value.
- 3. Compute the KS goodness-of-fit statistic for each ( α , b i min , i ).
- 4. Find the PL parameters (ˆ α, b ˆ min ) that minimize the KS statistic.

In case b min is set by the user to some constant value, the method yields a single corresponding slope α .

Using synthetic data, the authors of the method performed a number of tests of the PL hypothesis against alternative distributions like lognormal, exponential and a few others. Using the log-likelihood ratio as criterion, they found that the PL distribution can be correctly distinguished from lognormal for dataset sizes as large as ∑ k i =1 h i glyph[greaterorsimilar] 2 × 10 4 and bin-sizes b i +1 /b i glyph[lessorequalslant] 2 (Virkar &amp; Clauset 2014; Sect. 5).

Therefore bPlfit can be adopted as a reliable approach to extract PLTs from ρ -pdfs or N -pdfs given that the datasets are large enough. The critical issue is how to minimize the possible subjectivity introduced through the choice of total number of bins k . We develop a procedure to derive average PLT parameters, using logarithmic binning. The procedure will be presented in Section 4 and then applied to simulated data on evolving, self-gravitating MCs (Section 5) and on Herschel data for two Galactic regions of active star formation (Section 6).

## 3 USED NUMERICAL AND OBSERVATIONAL DATA

We make use of data from two sets of simulations (Sec. 3.1 and 3.2) and of dust-continuum observations by Herschel (Sec. 3.3). The simulations were performed using the hydrodynamical code Flash (Fryxell et al. 2000; Dubey et al. 2008). The first numerical data set consists of

1 http://tuvalu.santafe.edu/ ∼ aaronc/powerlaws/

2 http://tuvalu.santafe.edu/ ∼ aaronc/powerlaws/bins

a magneto-hydrodynamical (MHD) simulations from the SILCC project (Section 3.1), including feedback effects and tracing the evolution of giant MCs in a galactic environment. The scales cover range from 500 pc down to 0.12 pc. The second set are isothermal hydrodynamical (HD) simulations of self-gravitating clouds with supersonic turbulence at scales of typical large clumps (0 5 pc) in MCs, with adap-. tive resolution going as high as ∼ 3 au in high-density zones. The very high resolution allows to resolve the high-density end of the PLT; see Section 3.2 for details. Hereafter, this set of simulations is called HRIGT (High-Resolution Isothermal Gravo-Turbulent).

The observational data stem from Herschel imaging that provides dust column-density maps at angular resolutions lower than 1 . ′ N -pdfs obtained from these maps were presented for a large number of molecular cloud types, ranging from quiescent regions (Schneider et al. 2013) to lowand high-mass star-forming regions (e.g. Russeil et al. 2013; K¨ onyves et al. 2015; Schneider et al. 2015a,b; Stutz &amp; Kainulainen 2015; Rayner et al. 2017).

## 3.1 Galactic-scale (SILCC) simulations

The simulation framework including the physical processes is based on the SILCC project (Walch et al. 2015; Girichidis et al. 2016). The details of the particular setup as well as the dynamical evolution are described in Seifried et al. (2017) and Girichidis et al. (2018), so we only highlight the most important aspects of the simulations here. We set up a stratified gas distribution in a simulation box of 0 5 . × 0 5 . ×± 5 kpc 3 . The total gas surface density is Σ gas = 10M glyph[circledot] pc -2 , which results in a total mass of 2 5 . × 10 6 M in the box. The gas iniglyph[circledot] tially follows a Gaussian distribution with a scale height of 30 pc. Gravitational effects are included in two ways. We account for the stellar component of the galaxy using an external potential, which is implemented as an isothermal sheet with a stellar surface density of 30 M glyph[circledot] pc -2 . The second component of the gravitational force is self-gravity, which is computed with a tree-based solver (W¨nsch et al. 2018). u Radiative cooling is computed via a non-equilibrium chemical network (see e.g. Glover et al. 2010, for reference) that actively follows the abundances of ionized (H + ), atomic (H) and molecular (H ) hydrogen as well as carbon monoxide 2 (CO) and singly ionized carbon (C + ). The gas is initially purely atomic and at rest.

We let supernovae (SNe) explode at a constant rate of 15 Myr -1 . The latter is derived using the Kennicutt-Schmidt relation (Schmidt 1959; Kennicutt 1998) and assuming the initial mass function by Chabrier (2003). We follow the evolution of the gas for about 50 Myr at a maximum resolution of ∆ x = 3 9pc in order to identify regions in . which molecular clouds form. We then restart the simulations at t = 11 9 Myr and progressively increase the adaptive refine-. ment up to a maximum resolution of ∆ x = 0 12 pc while en-. suring that the Jeans length is always refined with at least 16 cells. The resolution needs to increase gradually as described in detail in Seifried et al. (2017). The simulations do not include sink particles but focus on the evolution of the gas phase.

## 3.2 Clump-scale (HRIGT) simulations

The second set of simulations focuses on smaller scales of star-forming regions. For this purpose we set up cubic periodic boxes with a length of 0 5 pc. . We reduce the complexity in physical modeling in favour of higher resolution and significantly higher adaptive refinement. We therefore omit the details of chemical phases and local shielding effects. The gas is isothermal at a temperature of 10 K. We assume a mean molecular weight of 2.3, which results in an isothermal sound speed of 0 19 km s . -1 . We set up a homogeneous gas distribution with two mean densities 〈 ρ 〉 (4 6 . × 10 -20 g cm -3 and 2 3 . × 10 -19 g cm -3 ) corresponding to a total mass of 85 and 426 M glyph[circledot] , respectively. Due to the initially uniform density and the isothermal temperature we can express the total mass in units of the initial average Jeans mass M J = 〈 ρ 〉 ( π/ 6) λ 3 J , where λ J = ( πc /G ρ s 〈 〉 ) with the isothermal sound speed c s . For our setup the Jeans masses are M J 0 , = 2 7 M . glyph[circledot] and M J 0 , = 1 2 M . glyph[circledot] for the two different densities. The total mass in the box is then 32 and 354 Jeans masses. The initial turbulent velocities are constructed in Fourier space with a peak of the power spectrum at k = 2, i.e. half of the box size. The spectral index for the velocity is set to -2 representing Burgers turbulence for strong shocks. We distinguish between purely compressive, purely solenoidal and naturally mixed velocities (Federrath, Klessen &amp; Schmidt 2008). The first case corresponds to radial projections in Fourier space, the second to tangential projections, and for the mixed case we do not apply any projection, which statistically results in approximately one third of the power in compressive and the remaining two thirds in solenoidal motions (Federrath, Klessen &amp; Schmidt 2008). We would like to emphasize that the nature of turbulent modes (compressive vs. solenoidal) has a strong impact on the time scales needed to trigger the formation of clouds and the onset of gravitational collapse in over-densities, see table 1. We only set the initial motions and let the turbulence decay. The root-mean-square velocity of the gas is set to 0 38 km s . -1 , i.e. the initial gas motions are supersonic with a Mach number of M = 2. This corresponds to a crossing time of approximately 1 3 Myr. We start with a uniform . grid of 256 3 cells and progressively refine up to an effective resolution of 32768 3 cells, such that the Jeans length is always refined with at least 16 cells. At the highest densities the cells thus have a size of 3 15 au. Self-gravity is included . from the beginning of the simulations, which causes the the gas to become locally unstable and collapse. At a threshold density of ρ = 2 12 . × 10 -13 g cm -3 we introduce Lagrangian sink particles based on Federrath et al. (2010).

## 3.3 Herschel observations of star-forming regions

We produced dust column-density maps of the Aquila and Rosette regions from publicly available level-3 data products produced with HIPE 13 from the Herschel archive, using the observations at 160, 250, 350 and 500 µ m. Aquila was observed within the Gould Belt keyprogram (Andr´ e et al. 2010) and presented in Bontemps et al. (2010), K¨nyves et o al. (2010, 2015) and Schneider et al. (2013). Rosette was part of the HOBYS keyprogram (Motte et al. 2010) and is discussed in Schneider et al. (2010), Schneider et al. (2012), Hennemann et al. (2010) and Di Francesco et al. (2010).

Table 1. HRIGT runs selected to elaborate the method for PLT extraction. Notation: v1, v2, v3 - chosen velocity field; vfmode - velocity field mode; m(ixed), c(ompressive), s(olenoidal) - acronyms for the vfmode (see text); M J 0 , - initial Jeans mass in the box; τ ff ( 〈 ρ 〉 ) - free-fall time of the mean density; t sim - run duration.

| Name     |   Total mass [ M glyph[circledot] ] [ M J , 0 ] |   Total mass [ M glyph[circledot] ] [ M J , 0 ] | vfmode      |   τ ff ( 〈 ρ 〉 ) [ kyr ] | t sim [ τ ff ( 〈 ρ 〉 ) ]   |
|----------|-------------------------------------------------|-------------------------------------------------|-------------|----------------------------|------------------------------|
| v1m-M085 |                                              85 |                                              32 | mixed       |                        310 | 1 . 05                       |
| v1c-M426 |                                             426 |                                             354 | compressive |                        139 | 0 . 41                       |
| v1s-M426 |                                             426 |                                             354 | solenoidal  |                        139 | 1 . 57                       |
| v2m-M085 |                                              85 |                                              32 | mixed       |                        310 | 2 . 50                       |
| v2m-M426 |                                             426 |                                             354 | mixed       |                        139 | 0 . 80                       |
| v3m-M085 |                                              85 |                                              32 | mixed       |                        310 | 1 . 38                       |

The angular resolution of the maps is 11 . ′′ 7, 18 . ′′ 2, 24 . ′′ 9, and 36 3 for 160, 250, 350, and 500 . ′′ µ m, respectively. For an absolute calibration of maps (included in the SPIRE level-3 data), the Planck -HFI observations were used for the HIPEinternal zeroPointCorrection task that calculates the absolute offsets, based on cross-calibration with HFI-545 and HFI-857 maps, including colour-correcting HFI to SPIRE wavebands, assuming a grey-body function with fixed β . The zero-point determination for the PACS 160 µ mmapwas performed using IRAS data, following the procedure described in Bernard et al. (2010). Column density and temperature maps were then produced at an angular resolution of 18 ′′ , following the procedure outlined in Palmeirim et al. (2013) which employs a multi-scale decomposition of the imaging data and assuming a constant line-of-sight temperature. We performed a pixel-by-pixel SED fit from 160 to 250 µ m, using a dust opacity law κ 0 = 0 1 . × ( ν/ 1000 GHz) β cm 2 g -1 which is similar to that of Hildebrand (1983) adopting β = 2 and assuming a gas-to-dust ratio of 100. (This dust opacity law is commonly adopted in other papers on the Gould Belt survey and HOBYS.) We estimate that the final uncertainties of the column-density map are around 20 - 30 %.

We selected the Aquila and Rosette molecular clouds because they represent examples of molecular clouds forming mostly stars of low and high mass, respectively. The area of the Aquila cloud that was covered with Herschel imaging is the western high extinction region of the Aquila Rift that extends over ∼ 20 ◦ located above the Galactic plane (b ∼ 4 ). It contains the young star cluster W40, associated ◦ with the H ii region Sharpless 2-64. The molecular mass of the Aquila cloud determined from the Herschel column density map within a radius of 9.8 pc is ∼ 2 × 10 4 M . The glyph[circledot] mass value assumes a distance of 260 pc (Strai˘ys z et al. 2003) for the entire Aquila complex, though a larger distance of around 415 pc (Dzib et al. 2010) is also discussed in the literature and recently confirmed by Gaia measurements to be 439 pc (Ortiz-Le´n et al. 2018). The Rosette o molecular cloud is an evolved region that already formed a cluster of 17 OB stars (NGC 2244), located in the center of the cloud. Nevertheless, there is continuous star formation along the central molecular ridge of the Rosette cloud (Hennemann et al. 2010; Schneider et al. 2012). The total mass of the complex from the Herschel column density map is ∼ 10 5 M , assuming a distance of 1.6 kpc (e.g. Turner glyph[circledot]

Figure 1. On the dependence of bPlfit output from the chosen number of logarithmic bins k . The HRIGT run v3m-M085 is selected for illustration. DPs (top) and slopes (bottom) of the ρ -pdf PLT as functions of k are plotted in the left column. The extracted short ( glyph[lessorequalslant] 5 bins) and long tails are indicated with open and filled circles, respectively. The open triangle denotes an example of ρ -pdf with extracted short tail (top right) while open square stands for a ρ -pdf with long tail (bottom right). In the right panels, short solid lines indicate the DP estimates and the derived slopes are plotted with dotted lines.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00489v1-figure-dd10143d-deca-473d-bc93-f44014ad69b7.jpg)

1976). With the adopted distance estimates, the linear resolutions of the maps are 0.04 pc for Aquila (1 pc=7 . ′′ 9) and 0.14 pc for Rosette (1 pc=2 . ′′ 15). The column-density maps are provided in Appendix A for reference.

## 4 ADAPTED METHOD TO EXTRACT PLTS FROM PDFS

## 4.1 Criterion for plausible PLT extraction

As bPlfit is applied to a given density distribution, one asks how does the result depend on the chosen total number of bins k . Fig. 1 (left) shows the bPlfit outputs with k varied in a large range, taking an example from the HRIGT runs. One clearly distinguishes two subsets of extremely different slopes/DPs. Within them, both PLT parameters do not seem to correlate with k . Examples of pdfs from each subset are shown in Fig. 1, right (denoted by open square and open triangle in the left panel). Although differences between the plotted distributions are hardly discernable by eye, bPlfit extracts a long, shallow PLT from the pdf displayed in the bottom-right panel and a very short, steep PLT from the pdf in the top-right panel. In the former case, the derived slope is not affected by the drop of statistics close to the upper data limit. In the latter case, the PL fit is based exactly on the last few bins; thus it reflects the incompleteness of data and should be considered as false. Inspection of the bPlfit outputs from all sampled numerical and observational data (Sections 3.1-3.3) showed that such false PLT extractions occur in an unpredictable way, for various choices of k ; sometimes they alternate with extractions of long PLTs as k is being increased (see Fig. 1, left). This is a peculiarity of the bPlfit technique itself and can not be overcome, e.g., by introducing some cut-off of the pdf at the high-density end. The simplest and most straightforward solution of the problem is to exclude false PLT extractions from further consideration. Hence a necessary first step for plausible PLT extraction should be to select pdfs with large enough PLT spans in units of bins b k +1 /b min &gt; FALSE where FALSE is the maximal span of false PLTs. The latter typically is only two 3 or three bins; thus FALSE = 5 bins would be a conservative and secure choice. PLTs that pass the false-extraction test (filled symbols in Fig. 1, left) can be sampled to calculate averaged DP and PL slope in the considered density distribution.

Second, one needs to quantify the minimal PLT span

3 Since two points determine a straight line (a PL function on a log-log plot) and bPlfit 'fits' them in case a larger PLT is not recognized.

( b k +1 /b min ) lim (hereafter, lower plausibility limit ) which defines a plausible bPlfit output in terms of fit's goodness. The choice of ( b k +1 /b min ) lim should provide a sufficient number of plausible PLT extractions to calculate an average DP and slope. The total statistics within the supposed PLT density range is important here, as illustrated in Fig. 2, top. In the case of rich statistics (circles), i.e., long and shallow tails, the PLT span increases quasi-monotonically with the chosen total number of bins k . Adopting a high value of ( b k +1 /b min ) lim (dash-dotted line) is affordable - it corresponds to some minimal k above which all PLT extractions are plausible. However, in case of poor statistics (triangles), i.e., short, unevolved tails, the relation between k and the PLT span is more complicated. The latter could be underestimated due to emergence of 'wiggles' or 'spikes' as k is increasing (bin size is decreasing). Artificially shorter (but still plausible) tails can be detected for some choices of k (shown with arrows in Fig. 2, top) due to a relatively low pdf value in the bin, adjacent from left to the estimated DP. Then one needs to choose a lower ( b k +1 /b min ) lim (dotted line) to get a larger sample of detected plausible PLTs for derivation of averaged bPlfit outputs.

Fig. 2, bottom, shows the results from tests of the criterion for plausible PLT extraction. We selected two evolutionary stages of a HRIGT run. At t = 0 5 . 〈 τ ff 〉 the PLT is still emerging, with steep extracted slopes ( -3 &gt; α glyph[greaterorsimilar] -6). Its span is small and one needs to adopt a lower value of ( b k +1 /b min ) lim (red line) to obtain a richer sample of PLT extractions. At an evolutionary stage as late as 0 97 . 〈 τ ff 〉 the PLT span b k +1 /b min increases fast and almost monotonically with the bin size, except one false PLT extraction for k = 75 bins. In this case, one could adopt a lower plausibility limit up to three decades.

For comparison, we applied bPlfit on analytical pdfs that are combinations of a lognormal function and a PLT and are good approximations of the pdfs from the chosen HRIGT run. In contrast to distributions from numerical/observational data, analytical pdfs are smooth and possess a bin-independent statistics 4 . The results for two pairs of analytical pdfs are plotted in Fig. 2, bottom (grey dashed lines). Thin lines indicate steep PLTs ( α = -3 6), . corresponding to t = 0 50 . 〈 τ ff 〉 , and thick lines denote shallow PLT slopes ( α = -1 9), . corresponding to t = 0 97 . 〈 τ ff 〉 . As expected, the PLT span of the analytical pdfs behaves similarly to that of the numerical distributions but without detections of false tails. The bPlfit method extracts correctly the slope, with vanishing deviations from the real value, while the typical deviations of the DP are of order of the bin size (see Appendix B).

## 4.2 Average PLT parameters

Now we are able adopt a procedure for extraction of PLT with an average span and slope, by use of bPlfit . The steps are as follows:

- (i) Choose the range of total bin number k . (Normally, k could be varied from a few dozens to about two hundred, depending on the given dataset.)
- 4 The pdf value in given bin always reflects the real distribution, regardless of how small the bin size is.
- (ii) Construct the pdf for each choice of k .
- (iii) Run bPlfit on each pdf and get the PLT parameters (slope and DP).
- (iv) Set a maximal span of false PLTs FALSE in units of bin size and exclude false PLT detections from further consideration. (Recommended conservative choice: FALSE = 5 bins.)
- (v) Set the lower plausibility limit ( b k +1 /b min ) lim in units of bin size for the PLT span, depending on the statistics in the supposed PLT range.
- (vi) Select a sample of plausible PLTs with spans ( b k +1 /b min ) glyph[greaterorequalslant] ( b k +1 /b min ) lim .
- (vii) Calculate average slope and DP from the sample.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00489v1-figure-4d70fc40-ad83-45e2-8919-10af70f7d840.jpg)

Figure 2. Span of the detected PLT as a function of the total number of bins. Top: Illustration on the cases of rich (circles) and poor (triangles) statistics in the supposed PLT range. Open symbols denote false PLTs with spans glyph[lessorequalslant] FALSE while artificially shorter but plausible PLT extraction are marked with arrows. Two choices of a ( b k +1 /b min lim ) are shown (dash-dotted and dotted lines). The choice of the lower value provides larger samples of plausible PLT detections. Bottom: On the criterion for sampling of pdfs with reliable PLT parameters. Two examples of bPlfit outputs from the HRIGT run are selected for illustration. They are compared with analytically defined pdfs (grey dashed lines) with approximately the same PLT slopes: α = -1 9 (thick) . or α = -3 6 . (thin). (See main text, for further details.) The shaded area denotes the zone of binned distributions that fail the false-extraction test with FALSE ≡ 5 bins; the adopted choice of ( b k +1 /b min lim ) = 8 bins is plotted with a red line.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00489v1-figure-2429c3a4-706f-4fb5-8eed-41a8ab741251.jpg)

This adapted bPlfit method is appropriate for extrac-

tion of PLTs from ρ -or N -pdfs in star-forming clouds. In the next Section we apply it to the selected data from numerical simulations as well to column-density data on two Galactic regions.

## 5 APPLICATION OF THE METHOD TO NUMERICAL DATA

In this Section we present the results from application of the adapted bPlfit method to the selected simulation runs. We define logarithmic density s = log( ρ/ρ 0 ) and logarithmic column-density z = log( N/N 0 ), using some arbitrary normalization units. The PLTs of the ρ -pdf and the N -pdf are introduced as follows:

$$p ( s ) \ = \ A _ { s } ( \rho / \rho _ { 0 } ) ^ { q } \quad \ \ ( 1 ) \quad \ \ \ t i$$

$$p ( z ) \ = \ A _ { z } ( N / N _ { 0 } ) ^ { n } \, \quad \ \ ( 2 ) \quad m$$

where A s and A z are constants. The adopted lower plausibility limit ( b k +1 /b min ) lim for the PLT span is 8 bins, both for the SILCC and the HRIGT data.

## 5.1 Pdf evolution on galactic scales

The chosen time span of the SILCC run is between 12 and 16 Myr, i.e. from the stage of emergence of high-density, selfgravitating clouds to the stage of formation of the first stars in them. The density and column-density distributions at the end of the run are shown in Fig. 3. Here, we would like to highlight the two zoom region with higher resolution because they contribute mostly to the high density PLT. Therein one clearly identifies dense ( ρ glyph[greaterorsimilar] 10 -22 g cm -3 ), elongated and filamentary structures with size dozens of pc. They can be considered as giant star-forming clouds at an early stage of their evolution. High-density cores within them are visible in the zoomed-in frames in Fig. 3; they display sizes of several parsec and N glyph[greaterorsimilar] 10 -2 g cm -2 (3 × 10 21 cm -2 ), typical for molecular clouds.

By applying the adapted bPlfit to the entire SILCC cube, we study the evolution of the ρ -pdf of the gas and of its atomic and molecular phases separately. Fig. 4 (right) presents the results in terms of the PLT parameters as the total number of bins is varied in the range 15 glyph[lessorequalslant] k glyph[lessorequalslant] 150. In the very early phase ( t glyph[lessorsimilar] 12 8 Myr) the atomic gas domi-. nates and its ρ -pdf displays a shallow PLT with DP in the diffuse regime ( ∼ 10 -24 cm -3 ), i.e. not attributable to gravitational contraction (see Fig. C1 in Appendix C and the comment therein). On the other hand, the ρ -pdf of H 2 is close to lognormal and no PLT is detected. Examples of the pdfs at two later evolutionary stages are shown in Fig. 4 (left). PLTs of the molecular gas well distinguishable from lognormal wings, with slopes -3 glyph[greaterorsimilar] q glyph[greaterorsimilar] -4, are detected at t ∼ 13 Myr (top). At that time both H i and H 2 contribute to the PLT of the total gas as read from the derived PLT parameters (cf. Fig. 4, right).

In the course of further evolution, the molecular gas dominates the density regime at ρ glyph[greaterorsimilar] 10 -20 g cm -3 (Fig. 4 , bottom left). The latter corresponds mostly to the filamentary clouds, visible in Fig. 3. The slope of the ρ -pdfs, both of H 2 and of the total gas, gets shallower and tends to stabilize about q glyph[greaterorsimilar] -1 5 . (bottom right). Its variations around this value are probably due to the complex interplay between gravity, turbulence and external pressure at the late evolutionary stage of the clouds and the formation of prestellar cores therein. We note that the variations of the average slope are more pronounced in the case with H 2 due to the smaller density span of the PLT. This behaviour is a characteristic of the adapted bPlfit method - the lower the DP value is (the larger the PLT span) the less the slope depends on local features (wiggles etc.) of the PLT. There is a hint of a PLT also for the atomic gas. If real, it is much steeper ( q ∼ -3) than for the molecular gas, but with a similar DP.

We also study the evolution of the N -pdf derived from column-density maps along the axis perpendicular to the galactic plane. Fig. 5 displays the N -pdfs of the total gas and its atomic and molecular phase at chosen evolutionary times. In contrast to the ρ -pdf, the N -pdf PLT span is almost entirely dominated by molecular gas. The N -pdf of the atomic gas remains lognormal during the whole studied period, while a PLT in the H 2 N -pdf emerges as early as t glyph[lessorsimilar] 13 Myr. The latter develops further as the slope gets shallower and the span extends to over two orders of magnitude. The small spikes which at the high-density end of the PLT are resolution effects and do not influence the slope substantially. The evolution of the N -pdf parameters for H 2 and for the total gas are shown in details in Fig. 6. Initially, the slope is steep: n ∼ -3 for H ; 2 n glyph[lessorsimilar] -3 for the total gas. At late evolutionary times it gets shallower and seems to approach constant values: around -2 for H , and 2 slightly lower, for the total gas. The DP shifts towards higher densities in the course of evolution of the H 2 pdf and varies around a constant value, similar to the stabilization of the slope. On the other hand, the DP of the total-gas N -pdf is remarkably constant with the evolution. This results from the practically unchanging central, quasi-lognormal part of the atomic-gas N -pdf.

## 5.2 Pdf evolution in self-gravitating clumps

The HRIGT simulations allow us to follow the evolution of ρ -pdf and N -pdf in self-gravitating media at scales glyph[lessorequalslant] 0 5 pc, . typical for clumps within molecular clouds. Figure 7 presents the column density maps at the end of the HRIGT runs. Clumps with moderate column density form in the runs with lower density (top panels), with total mass in the box 32 M J 0 , (cf. Table 1). Filamentary structures containing a number of dense cores appear in the runs with many Jeans masses (354 M J 0 , ; bottom panels). The runs with lower total number of Jeans masses (top) show the formation of one to a few overdensities. The majority of the volume is filled with diffuse gas. If the box contains numerous Jeans masses of gas (bottom) relatively weak turbulent velocities are needed to trigger the formation of independently collapsing clumps. The driving mode of the turbulent velocity field results in different times to trigger the collapse. Compressive motions (left-hand panel) push the gas to form filaments with strong density contrast. In contrast, solenoidal forcing (right-hand panel) leads to collapsing clumps embedded in diffuse gaseous structures, see also Girichidis et al. (2011). Despite the differences in the spatial gas distribution, the evolution of the PLT parameters in all runs turns out to be quite uniform. It is illustrated in Fig. 8. In each

Figure 3. Snapshot of density (left) and column-density (right) field at the end of the SILCC run. Two zones of higher resolution, containing evolved GMCs are zoomed in. The slice of the density field is taken at the galactic plane ( z = 0). We highlight the highresolution regions, which mainly contribute to the high-density part of the PDF.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00489v1-figure-b8499f37-d8b3-4ecf-a9fd-f0c8b7070d17.jpg)

Figure 4. Evolution of the ρ -pdf PLT from the SILCC data. Left: ρ -pdfs at the evolutionary stages with an emerging (top) and well developed (bottom) PLT; Right: Evolution of the average PLT parameters. The error bars indicate the minimal/maximal value in the sample.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00489v1-figure-1374e8e3-6f31-440a-bcfd-1fe1b8abe595.jpg)

Figure 5. Evolution of the N -pdf as derived from the SILCC simulation. Colors denote the pdf of the total gas (black), of H i (red) and of H 2 (blue).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00489v1-figure-8667a700-ae3e-4877-b825-2291e3526e91.jpg)

run, the initial velocity field introduced in the data cube leads to formation of an approximately lognormal ρ -pdf (not shown) within a few crossing times. A distinguishable PLT with slope q glyph[greaterorsimilar] -4 appears later in the evolution. At the final stage, as first stars/sink particles form, q seems to vary around some constant value: ∼ -2 in the runs v1c-M426 and v1s-M426 and ∼ -1 5, in the other cases (Fig. 8, bot-. tom left). On the other hand, the DPs also tend towards some constant value, though with variations up to about one order of magnitude (top left). The latter could be due mostly to foreground and background contamination.

The evolution of the PLT parameters of the N -pdfs is qualitatively similar to the one of the ρ -pdfs (Fig. 8, right panels). However, the relative duration of the stage characterised by distinguishable PLTs ( n glyph[greaterorsimilar] -4) - if detected at all - is much shorter. At this stage the slope shallows rapidly, with a noticeable trend towards some constant value ( n ∼ -2 5) in the runs v1m-M085 and v2m-M085. The range . of DP variations at the final evolutionary stage is tighter compared to the ρ -pdf analysis: about half an order of magnitude in a chosen run.

## 6 APPLICATION OF THE METHOD TO HERSCHEL DATA

The adapted bPlfit method was also applied to the N -pdfs in the star-forming regions Aquila and Rosette (Sec. 3.3). Aquila is a filamentary cloud complex of recent starforming activity where hundreds of starless cores and young protostellar objects have been identified (K¨ onyves et al. 2010). Star formation in the Rosette molecular cloud probably follows two distinct mechanisms. It can be triggered by the expanding H ii region in a narrow area in the H ii region/molecular cloud interface zone, but it also takes place deep inside the molecular cloud where filaments merge and provide mass input for the formation of OB clusters (Schneider et al. 2012). The cloud structure is overall highly filamentary and the clumps identified from molecular-line and

Figure 6. Evolution of the N -pdf PLT parameters from the SILCC data.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00489v1-figure-adb358e2-b7d3-436f-bcab-2a5b72a6a63c.jpg)

dust emission turn out to be mostly gravitationally bound (Veltchev et al. 2018). We can thus expect that the combination of large-scale self-gravity (filament merging) and small-scale gravitational collapse leads to a clear PL tail; see Sec. 7 and the references therein.

Indeed, the PLTs extracted through the adapted bPlfit from the N -pdfs in Aquila and Rosette are evident and span about one order of magnitude (Fig. 9). The slopes and the DPs are in general agreement with the values obtained from the SILCC data, for evolved giant molecular clouds (cf. Fig. 6). Similar slopes have been found, by use of other techniques, by K¨nyves et al. (2015) in Aquila o ( n = -2 91) and by Schneider et al. (in prep.) in Rosette . ( n = -2 73). . A steeper slope of n = -3 23 was found in . Rosette by Tremblin et al. (2014) using a lower angular resolution map (36 ′′ ) and fitting a narrower range. We note that the PLT parameters are not affected by local features in the tail. There are some resolution effects at its high-density end but they do not influence the slope, as commented in the next Section. The non-PL part of the N -pdf in Rosette has a quasi-lognormal shape while the case in Aquila seems to be more complex. How the output of the adapted bPlfit compares with PL tails of observational N -pdfs, extracted by use of other approaches (typically, by assuming a combination of lognormal part and a PLT), is worth of a study on its own. We are going to present such comparison in a forthcoming paper (Schneider et al., in preparation).

Figure 7. Snapshots of the column density at the end of each HRIGT run. The upper panel shows the boxes with lower masses and mixed turbulent driving. In all runs we note the formation of one to a few dense clumps. The bottom panel depicts the runs with higher mass. Here we find numerous clumps, which we attribute mainly to the larger number of Jeans masses in the box. The difference in turbulent driving mode is visible in the filamentary structure. Compressive modes (lower left-hand panel) lead to stronger filaments with larger density contrast. Solenoidal driving (lower right-hand panel) forms more diffuse over-densities. The differences in turbulent driving also lead to significant differences in the total simulation time, see also Girichidis et al. (2011). Note the different scales in the column density range.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00489v1-figure-a858d37a-f044-42bb-ae7f-c6c09b216069.jpg)

## 7 DISCUSSION

An obvious advantage of the proposed method is that it rests on a sole assumption about the ρ -pdf/ N -pdf: the existence of a PLT . No additional conjectures about the distribution's shape at lower densities are required. The PLT parameters are derived simultaneously which removes the uncertainty of the DP introduced by its assessment by eye. The procedure of averaging of the PLT parameters minimizes the uncertainties caused by some arbitrariness of the binning scheme. Moreover, the adapted bPlfit method is not sensitive to the low-density completeness limit of the data, as long the latter is lower than the N -pdf peak. Therefore the assessed PLT parameters can be considered as reliable.

able assessment of the PLT parameters requires that the tail spans at least one order of magnitude, with a statistically significant number of measurement points providing very low Poisson noise. Poor resolution in dense zones of starforming regions could lead to distortions or artificial shortening of the tail and hence impede the procedure of averaging of the PLT parameters, even if one sets ( b k +1 /b min ) lim as low as a few bins (cf. Section 4).

A caveat of the method is the unpredictable but - typically - rare extraction of false, short PLTs by bPlfit for some choices of the total number of bins. In fact, this results from fitting at the very high-density end of the pdf, where the statistic drops. Setting of appropriate maximal span of false PLTs FALSE to exclude them solves the problem. Then a reli-

In case the observed N -pdfs display two PLTs (e.g. Russeil et al. 2013; Schneider et al. 2015b; Pokhrel et al. 2016,) the adapted bPlfit method extracts the DP of the lowerdensity one while the slope estimate is not significantly affected by the existence of another PLT. In other words, the method can not extract two PLTs simultaneously - given the first PLT is detected, the procedure should be reapplied only to the PL part of the N -pdf. The choice of the lower density cut-off is critical here and it needs further investigation, in a follow-up work of this study.

The PLT parameters obtained by applying adapted

Figure 8. Evolution of the parameters of PLTs extracted from the ρ -pdfs (left) and the N -pdfs (right) from the HRIGT runs. The time is normalised to the run duration, specified in Table 1. The time axis for t glyph[greaterorequalslant] 0 9 . 〈 〉 t ff is stretched exponentially.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00489v1-figure-7dbf6460-362d-4f08-89af-9f5ab2825dc6.jpg)

bPlfit to the selected HRIGT and SILCC data are consistent with the results from other numerical and observational studies. Values of the ρ -pdf slope q ∼ -1 5, both on galac-. tic and clump scales, are derived at late evolutionary stages of self-gravitating star-forming clouds (Kritsuk, Norman &amp; Wagner 2011; Collins et al. 2012; Federrath &amp; Klessen 2013), as expected also from theoretical considerations (Girichidis et al. 2014; Donkov &amp; Stefanov 2018). On the other hand, N -pdfs derived from observations of regions with star-forming activity display pronounced PLTs of slopes -2 glyph[greaterorsimilar] n glyph[greaterorsimilar] -4 (Abreu-Vicente et al. 2015; Schneider et al. 2013, 2015a; Pokhrel et al. 2016), also in agreement with our results on the PLT evolution from the SILCC runs (Fig. 6).

As shown in some of the abovementioned studies, the slope n of the N -pdf should be related to q as:

$$n = 2 q / ( 3 + q ) \, \quad \ \ \ ( 3 ) \quad \ ^ { \iota },$$

that the PLTs are clearly distinguishable from lognormal wings, i.e. for slopes glyph[greaterorsimilar] -4. Note that the latter does not hold for the earliest and relatively short stage characterised by | q | glyph[greaterorsimilar] 2 5. . The slopes derived from the SILCC runs are also in good agreement with formula (3); again, with a few exceptions at the earliest studied evolutionary stage. A systematic shift of a few dex from the model prediction (dashed line) is evident, mostly considering the pdf evolution of the total gas (Fig. 10, right; black symbols). We interpret this shift in the upward direction - the N -pdf PLT slopes | n | are larger than expected while the obtained values of q at later evolutionary times are about the expected one of -1 5 . (Girichidis et al. 2014). A reasonable explanation might be that at scales about the typical sizes of giant MCs, the relation (3) changes due to the contribution of collapsing filaments to the PLT of N -pdf. In general, we point out that the adapted bPlfit method extracts PLTs of the ρ -pdfs and N -pdfs with slopes which are mutually consistent.

on the assumption that the general cloud structure can be described through a power-law density profile (see Donkov, Veltchev &amp; Klessen 2017, and the references therein). Figure 10 presents a test of this relation using the derived values of q and n at different evolutionary times from the HRIGT and SILCC runs. Most of the clump-scale simulations (left) fit well with the prediction of the theory given

## 8 SUMMARY

We present a new approach to extract the power-law tails (PLTs) of the density ( ρ -pdf) or column-density distribution ( N -pdf) in star-forming clouds. It combines the mathemat-

Figure 9. Extracted PL tails (dashed lines) from N -pdfs in two Galactic star-forming regions. The conversion to column density is done using the standard factor N (H ) 2 /A V = 0 94 . × 10 21 cm -2 mag -1 (Bohlin et al. 1978); mean particle mass of 1 37 atomic . mass units is adopted to calculate the DPs in [ g/cm 2 ].

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00489v1-figure-3b4d154d-002a-4531-b15a-fce4ab534f5b.jpg)

ical method bPlfit (Virkar &amp; Clauset 2014) which assesses the power-law part of an arbitrary binned distribution with a procedure of averaging of the output parameters: slope and deviation point (DP), as the number of bins is varied. Reliable results are obtained as long as the PLT spans at least one order of magnitude.

Compared with other techniques for PLT extraction, the proposed adapted bPlfit method has two considerable advantages. First, it is not based on assumptions about the rest of the pdf (e.g. lognormal) to calculate the slope and the DP - the sole assumption is that a PLT actually exists. Second, the output parameters are not sensitive to local features ('spikes' etc.) of the distribution tail.

The adapted bPlfit was tested on data from numerical simulations of star-forming clouds at clump scale (0 5 pc; . self-gravitating isothermal medium) and at galactic scales (500 pc; stratified gas distribution, chemistry and feedback included). In both cases, the emergence of a PLT takes place at late evolutionary times when protostars (sink particles) and/or dense filamentary clouds are formed. In the course of further cloud evolution, the slopes and the DPs tend to approach constant values - in agreement with a number of theoretical and numerical studies on the ρ -pdf and N -pdf in star-forming regions.

Figure 10. Theoretical relation between the PLT slopes of ρ -pdf and N -pdf (dashed line; see equation 3), juxtaposed with the results from our numerical runs. The color and symbol notation in the right panel is the same as in Fig. 6 and in the left one like in Fig. 8.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00489v1-figure-b8ef63bd-43fe-48fb-ace3-25a314906025.jpg)

The adapted bPlfit method was applied also to N -pdf of the regions Aquila and Rosette, observed with Herschel . We found pronounced PLTs with slopes as shallow as ∼ -2 7, consistent with those for evolved clouds from the used . galactic-scale simulation.

## ACKNOWLEDGEMENT

We are grateful to the anonymous referee whose critical comments and constructive suggestions helped us to improve this paper.

T.V. and S.D. acknowledge support by the Deutsche Forschungsgemeinschaft (DFG) under grant KL 1358/201 and additional funding from the Ministry of Education and Science of the Republic of Bulgaria, National RI Roadmap Project DO1-157/28.08.2018. S.D. acknowledges support by the Bulgarian National Science Fund under Grant N 12/11 (20.12.2017). P.G. acknowledges funding from the European Research Council under ERC-CoG grant CRAGSMAN-646955. T.V. and N.S. acknowledge support from the ANR/DFG collaboration project GENESIS (ANR16-CE92-0035-01/DFG1591/2-1). O.S. and L.M. acknowledge support by the Scientific Research Fund of the University of Sofia under Grant #80-10-68/19.04.2018. R.S.K. thanks funding from the DFG in the Collaborative Research Center (SFB 881) 'The Milky Way System' (subprojects B1, B2, and B8) and in the Priority Program SPP 1573 'Physics of the Interstellar Medium' (grant numbers KL 1358/18.1, KL 1358/19.2) D.S. acknowledges the support of the Bonn-Cologne Graduate School, which is funded through the German Excellence Initiative. D.S. also acknowledges funding by the DFG via the Collaborative Re-

search Center SFB 956 Conditions and Impact of Star Formation (subprojects C5 and C6).

The Flash code used in this work was partly developed by the Flash Center for Computational Science at the University of Chicago. The authors acknowledge the LeibnizRechenzentrum Garching for providing computing time on SuperMUC via the project pr94du as well as the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu).

## REFERENCES

Abreu-Vicente, J., Kainulainen, J., Stutz, A., Henning, Th., Beuther, H., 2015, A&amp;A, 581, 74

- Andr´ e, P., Men'shchikov, A., Bontemps, S., K¨nyves, V., Motte, o F., Schneider, N., Didelon, P., Minier, V., et al., 2010, A&amp;A, 518, L102

Alves, J., Lombardi, M., Lada, C., 2017, A&amp;A, 606, L2

Auddy, S., Basu, S., Kudoh, T., 2018, MNRAS, 474, 400

Ballesteros-Paredes, J., V´zquez-Semadeni, E., Gazol, A., Hart-

a

mann, L., Heitsch, F., Col´ ın, P., 2011, MNRAS, 416, 1436

- Bernard, J.-Ph., Paradis, D., Marshall, D. J., et al. 2010, A&amp;A, 518, L88

Bohlin, R., Savage, B., Drake, J., 1978, ApJ, 224, 132

Bontemps, S., Andr´, Ph., K¨nyves, V., Men'shchikov, A., Schneie o der, N., Maury, A., Peretto, N., Arzoumanian, D., et al., 2010, A&amp;A, 518, L85

Brunt, C. M., 2015, MNRAS, 449, 4465

Chabrier G., 2003, 115, 763

- Chen, H., Burkhart, B., Goodman, A., Collins, D., 2018, ApJ, 859, 162
- Clauset, A., Shalizi, C. R., Newman, M. E. J., 2009, SIAM Rev., 51, No. 4, 661
- Collins, D., Kritsuk, A., Padoan, P., Li, H., Xu, H., Ustyugov, S., Norman, M., 2012, ApJ, 750, 13
- Di Francesco, J., Sadavoy, S., Motte, F., Schneider, N., Hennemann, M., Csengeri, T., Bontemps, S., Balog, Z., et al., 2010, A&amp;A, 518, L91

Donkov, S., Stefanov, I., 2018, MNRAS, 474, 5588

Donkov, S., Veltchev, T., Klessen, R. S., 2017, MNRAS, 466, 914 Dzib, S., Loinard, L., Mioduszewski, A. J., et al. 2010, ApJ, 718, 610

Federrath C., Klessen R. S., 2013, ApJ, 763, 51

- Federrath C., Klessen R. S., Schmidt, W., 2008, ApJ, 688, L79
- Federrath C., Roman-Duval J., Klessen R. S., Schmidt W., Mac Low, M. 2010, A&amp;A, 512, A81

Fryxell B., Olson K., Ricker P. et al. 2000, ApJS, 131, 273

Dubey A. et al., 2008, in Pogorelov N.V. Audit E., Zank G.P., eds, ASP Conf. Ser. Vol. 385, Numerical Modeling of Space Plasma Flows: Astronum 2007. Astron. Soc. Pac., San Francisco. p. 145

- Girichidis P., Federrath C., Banerjee R., Klessen R. S., 2011, MNRAS, 413, 2741
- Girichidis P., Konstandin L., Whitworth A. P., Klessen, R. S., 2014, ApJ, 781, 91
- Girichidis P., Walch S., Naab T., Gatto A., W¨nsch R., Glover u S. C. O., Klessen R. S., Clark P. C., et al., 2016, MNRAS, 456, 3432

Girichidis P., Seifried D., Naab T., Peters T., Walch S., W¨nsch u R., Glover S. C. O., Klessen R. S., 2018, MNRAS, 480, 3511 Glover S., Federrath C., Mac Low M.-M., Klessen R. S., 2010, MNRAS, 404, 2

- Federrath, C., Banerjee, R., Clark, P. C., Klessen, R. S., 2010, ApJ, 713, 269
- Hennemann M., Motte F., Bontemps S., Schneider N., Csengeri T., Balog Z., di Francesco J., Zavagno A., et al., 2010, A&amp;A, 518, L84

Hildebrand R., 1983, QJRAS, 24, 267

Kainulainen J., Tan J. C., 2013, A&amp;A, 549, A53 Kainulainen J., Beuther H., Henning T., Plume R. 2009, A&amp;A, 508, L35 Kainulainen J., Beuther H., Banerjee R., Federrath C., Henning T., 2011, A&amp;A, 530, A64

Kennicutt Jr. R. C., 1998, ApJ, 498, 541

Klessen R. S., 2000, ApJ, 535, 869

Konstandin L., Federrath C., Klessen R. S., Schmidt W. 2012, JFM, 692, 183

- K¨ onyves V., Andr´ Ph., Men'shchikov A., Schneider N., Arzoue manian D., Bontemps S., Attard M., Motte F., 2010, A&amp;A, 518, L106
- K¨ onyves V., Andr´ Ph., Men'shchikov A., Palmeirim P., Arzoue mania D., Schneider N., Roy A., Didelon P., et al., 2015, A&amp;A, 584, A91
- K¨ ortgen B., Federrath C., Banerjee R., 2019, MNRAS, 482, 5233 Kritsuk A., Norman M., Padoan P., Wagner R., 2007, ApJ, 665, 416

Kritsuk A., Norman M., Wagner R., 2011, ApJL, 727, L20 Li Y., Klessen R., Mac Low M.-M., 2003, ApJ, 592, 975 Lombardi M., Alves J., Lada C. J., 2015, A&amp;A, 576, L1 Molina F., Glover S., Federrath C., Klessen R. S.,2012, MNRAS,

- Matsumoto T., Dobashi K., Shimoikura T., 2015, ApJ, 801, 77 423, 2680
- Motte F., Zavagno A., Bontemps S., Schneider N., Hennemann M., Di Francesco J., Andr´ e Ph., Saraceno P., et al., 2010, A&amp;A, 518, L77
- Ortiz-Le´n G., Loinard L., Dzib S., Kounkel M., Galli Ph., Tobin o J., Evans N. J., II, Hartmann L., et al., 2018, ApJ, 869, 33O Ossenkopf-Okada V., Csengeri T., Schneider N., Federrath C., Klessen R. S., 2016, A&amp;A, 590, A104

Padoan P., Nordlund A., Jones B., 1997, MNRAS, 288, 145

- Palmeirim P., Andr´ Ph., Kirk J., Ward-Thompson D., Arzoumae nian D., K¨nyves V., Didelon P., Schneider N., et al., 2013, o A&amp;A, 550, A38

Passot T., V´zquez-Semadeni E. 1998, PhRvE, 58, 4501 a

Pokhrel R., Gutermuth R., Ali B., Megeath T., Pipher J., Myers P., Fischer W. J., Henning T., et al., 2016, MNRAS, 461, 22 Rayner, T., Griffin, M., Schneider, N., Motte, F., K¨nyves, V., o Andr´ e, Ph., Di Francesco, J., Didelon, P., et al., 2017, A&amp;A, 607, A22

Russeil D., Schneider N., Anderson L., Zavagno A., Molinari S., Persi P., Bontemps S., Motte F., et al., 2013, A&amp;A, 554, 42 Schmidt M., 1959, ApJ, 129, 243

- Schneider N., Motte F., Bontemps S., Hennemann M., di Francesco J., Andr´ Ph., Zavagno A., Csengeri T., et al., 2010, e A&amp;A, 518, L83
- Schneider N., Csengeri T., Hennemann M., Motte F., Didelon P., Federrath C., Bontemps S., Di Francesco J., et al., 2012, A&amp;A, 540, L11

Schneider N., Andr´ e Ph., K¨ onyves V., Bontemps S., Motte F., Federrath C., Ward-Thompson D., Arzoumanian D., et al. 2013, ApJ, 766, L17

- Schneider N., Ossenkopf V., Csengeri T., Klessen R. S., Federrath C., Tremblin P., Girichidis P., Bontemps S., et al., 2015a, A&amp;A, 575, A79

Schneider N., Bontemps S., Girichidis P., Rayner T., Motte F., Andr´ e Ph., Russeil D., Abergel A., et al., 2015b, MNRAS, 453, L41

- Schneider N., Ossenkopf-Okada V., Roy A., Kritsuk A., R¨llig M., o Seifried D., Franeck A., Klessen R.S., et al., 2019, submitted to MNRAS

Seifried D., Schmidt W., Niemeyer J., 2011, A&amp;A, 526, A14 Seifried D., Walch S., Girichidis P., Naab T., W¨nsch R., Klessen u R. S., Glover S. C. O., Peters T., Clark P., 2017, MNRAS, 472, 4797

- Slyz A., Devriendt J., Bryan G., Silk J., 2005, MNRAS, 356, 737 Stanchev O., Veltchev T., Kauffmann J., Donkov S., Shetty R.,

## 14 Veltchev et al.

- K¨ ortgen B., Klessen R. S., 2015, MNRAS, 451, 5575 Strai˘ys V., Cernis K., Barta˘i¯te, S., 2003, A&amp;A, 405, 585 z ˘ s u Stutz A., Kainulainen J., 2015, A&amp;A, 577, L6 Tremblin P., Schneider N., Minier V., Didelon P., Hill T., Anderson L. D., Motte F., Zavagno A., 2014, A&amp;A, 564, A106 Turner D.G., 1976, ApJ, 210, 65 V´ azquez-Semadeni E., 1994, ApJ, 423, 681 V´ azquez-Semadeni E., Gazol A., Scalo J., 2000, ApJ, 540, 271 V´ azquez-Semadeni E., Gonz´lez R., Ballesteros-Paredes J., Gazol a A., Kim J. 2008, MNRAS, 390, 769 Virkar Y., Clauset A., 2014, Annals of Appl. Stat., vol. 8, No. 1, 89 (arXiv:1208.3524) Veltchev T., Ossenkopf-Okada V., Stanchev O., Schneider N., Donkov S., Klessen R. S., 2018, MNRAS, 475, 2215 Walch S., Girichidis P., Naab T., Gatto A., Glover S. C. O., W¨ unsch R., Klessen R. S., Clark P. C., et al., 2015, MNRAS, 454, 238 Ward R. L., Wadsley J., Sills A., 2014, MNRAS, 445, 1575 W¨ unsch R., Walch S., Dinnbier F., Whitworth A., 2018, MNRAS, 475, 3393

Figure A1. Herschel maps of Aquila (top) and Rosette (bottom) used for construction of the corresponding N -pdfs (cf. Section 3.3).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00489v1-figure-2800adc7-baf8-423a-b240-405923c8abcb.jpg)

## APPENDIX A: COLUMN-DENSITY MAPS

Column-density maps of the studied star-forming regions are shown in Fig. A1. Note that they are rotated so that no absolute coordinates are given. We express the column density in visual extinction, using the conversion N H ( ) = A V 0 94 . × 10 21 cm -2 mag -1 Bohlin et al. (1978).

Figure B1. Test of bPlfit on analytical pdfs with known PLT parameters. See text.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00489v1-figure-75d6a357-f054-419d-9b74-796655225284.jpg)

## APPENDIX B: TEST OF bPlfit ON ANALYTICAL PDFS

We test the accuracy of the bPlfit technique by applying it to two series of analytical pdfs which are combinations of a lognormal function and a PLT. Pdfs with a PLT slope -3 6 . (short, unevolved tail) are illustrative for early evolutionary state of the cloud while those with slope -1 9 are typical for . later evolutionary stages. The results for two pairs from both groups are shown in Fig. B1; each pair is characterized by the same slope while the PLT span is different. The deviations of the extracted PLT parameters from the real values do not seem to depend on the total number of bins k . The slopes are derived with a high average accuracy of 10 -3 (bottom panel) whereas the typical deviation of the DP from the real value (top panel) is between 1 and 4 bins though it may be larger in some exceptional cases.

## APPENDIX C: PLTS AT AN EARLY EVOLUTIONARY STAGE

As seen in Fig. 4, the ρ -pdfs of the atomic and the total gas at early stages of the SILCC run are characterized by very shallow PLTs ( q &gt; -1). An example is shown in Fig. C1, bottom. The DP coincides with a pronounced maximum of the distribution and corresponds to a typical density of the diffuse zones ( ρ ∼ 10 -24 g cm -3 ) which separates large 'voids' of very diluted, mainly atomic gas from the emerging

)

(

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00489v1-figure-fe2e8f6a-8e19-409c-8d50-536bbda88c92.jpg)

t

= 11.91 Myr

Figure C1. Slice of the density field ( z = 0) of the total gas (H i + H + H 2 + ) in the SILCC run (top) and the corresponding ρ -pdfs of different tracers (bottom) before we start the additional refinements in the zoom-in regions. The DP (short vertical mark) and the shallow slope (dashed line) are denoted. See text.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00489v1-figure-3511601b-9627-43be-9f2e-9ee27acf65a6.jpg)

giant MCs (cf. top panel). The ρ -pdf at this stage is still influenced by the initial conditions. The initial setup results in a flat density ρ -pdf with an initial peak at ∼ 10 -23 g cm -3 . During the initial phase with cooling and a vertical contraction of the disc a non-negligible fraction of the gas is compressed to higher densities. At the same time the SN driving creates over-densities via shock waves, which results in an accumulation of gas at ∼ 10 -24 g cm -3 . The emerging peak in the pdf at that density effectively leads to flat slopes for the higher densities. As this behaviour mainly reflects the initial conditions rather than any physical processes in equilibrium or steady state we should not overinterpret the value of the slope.
Todor V. Veltchev Philipp Girichidis Sava Donkov Nicola Schneider Orlin Stanchev Lyubov Marinkova Daniel Seifried Ralf S. Klessen
On the extraction of power-law parts of the probability density functions in star-forming clouds
We present a new approach to extract the power-law part of a
density/column-density probability density function (rho-pdf/N-pdf) in
star-forming clouds. It is based on the mathematical method bPLFIT of Virkar &
Clauset (2014) and assesses the power-law part of an arbitrary distribution,
without any assumptions about the other part of this distribution. The slope
and deviation point are derived as averaged values as the number of bins is
varied. Neither parameter is sensitive to spikes and other local features of
the tail. This adapted bPLFIT method is applied to two different sets of data
from numerical simulations of star-forming clouds at scales 0.5 and 500 pc and
displays rho-pdf and N-pdf evolution in agreement with a number of numerical
and theoretical studies. Applied to Herschel data on the regions Aquila and
Rosette, the method extracts pronounced power-law tails, consistent with those
seen in simulations of evolved clouds.