Multiple power-law tails in the density and column-density distribution in contracting star-forming clumps

We present a numerical study of the evolution of power-law tails (PLTs) in the (column-)density distributions ($N$-PDF, $\rho$-PDF) in contracting star-forming clumps in primordial gas, without and with some initial rotational and/or turbulent support. In all considered runs multiple PLTs emerge shortly after the formation of the first protostar. The first PLT (PLT 1) in the $\rho$-PDF is a stable feature with slope $q_1\simeq -1.3$ which corresponds -- under the condition of preserved spherical symmetry -- to the outer envelope of the protostellar object with density profile $\rho\propto l^{-2}$ in the classical Larson-Penston collapse model, where $l$ is the radius. The second PLT (PLT 2) in the $\rho$-PDF is stable in the pure-infall runs but fluctuates significantly in the runs with initial support against gravity as dozens of protostars form and their mutual tidal forces change the density structure. Its mean slope, $\langle q_2\rangle\simeq -2$, corresponds to a density profile of $\rho\propto l^{-3/2}$ which describes a core in free fall in the classical Larson-Penston collapse model or an attractor solution at scales with dominating protostellar gravity. PLT 1 and PLT 2 in the $N$-PDFs are generally consistent with the observational data of Galactic low-mass star-forming regions from {\it Herschel} data. In the runs with initial support against gravity a third PLT (PLT~3) in the $\rho$-PDFs appears simultaneously with or after the emergence of PLT 2. It is very shallow, with mean slope of $\langle q_3\rangle\simeq -1$, and is associated with the formation of thin protostellar accretion disks.


INTRODUCTION
A large variety of investigation tools have been proposed to study the physics and the evolution of molecular clouds in their variety of star-forming activity.Widely used in the last two decades is the analysis of the probability density function (PDF) of mass density (ρ-PDF; easily constructed from numerical data) and of column density (N -PDF).The popularity of this approach rests on many studies which have shown that the PDF shape and its distinctive parts bear imprints of the physical conditions and processes in the clouds.Supersonic turbulent flows in isothermal, nongravitating medium lead to a lognormal ρ-PDF (Vázquez-Semadeni 1994; Padoan, Nordlund & Jones 1997;Passot & Vázquez-Semadeni 1998;Klessen 2000;Li, Klessen & Mac Low 2003;Kritsuk et al. 2007;Federrath, Klessen & Schmidt 2008;Molina et al. 2012), with some possible deviations from ⋆ E-mail: eirene@phys.uni-sofia.bglognormality in the low-density range (Federrath et al. 2010;Konstandin et al. 2012;Pan, Padoan & Nordlund 2019).A more complex, bimodal ρ-PDF with two peaks might result from thermal instabilities in purely hydrodynamic flows (Vázquez-Semadeni, Gazol & Scalo 2000).In the high-density range, a power-law tail (PLT) of the ρ-PDF can appear in case of non-isothermal compressive turbulence, with polytropic index less than unity (Scalo et al. 1998;Passot & Vázquez-Semadeni 1998).But most typically, the development of high-density PLT out of a quasi-lognormal ρ-PDF reflects the increasing role of self-gravity in evolving starforming clouds (Klessen 2000;Vázquez-Semadeni et al. 2008;Kritsuk, Norman & Wagner 2011;Collins et al. 2012;Federrath & Klessen 2012, 2013;Burkhart et al. 2017).The emerging PLT is steep and hardly distinguishable from a lognormal wing; later, in the course of local collapses, its slope becomes shallower and tends toward a constant value whereas its span increases toward lower densities (Veltchev et al. 2019;Khullar et al. 2021).This evolution was studied and substan-tiated theoretically by Girichidis et al. (2014), from analytical model of a collapsing homogeneous sphere, and by Jaupart & Chabrier (2020) who showed that the PLT density range corresponds to regions in which gravity is the dominant agent in the statistics of turbulence.
The derivation of N -PDF from observations of star-forming regions is not straightforward.Often, the column density is calculated from the measured visual extinction or thermal dust emission; this involves some assumptions about dust properties and constraints on the column-density range and introduces uncertainties which are difficult to quantify.Nevertheless, observational N -PDFs can be appropriately corrected for effects of foreground contamination, observational noise and incompleteness (Ossenkopf-Okada et al. 2016).They display a similar functional form like the ρ-PDFs predicted from numerical and theoretical studies: a largely lognormal main part and a single PLT at the high-density end (e.g., Kainulainen et al. 2009;Froebrich & Rowles 2010;Pokhrel et al. 2016;Schneider et al. 2013Schneider et al. , 2015b) whose slope appears to correlate with the current stage of the cloud in the star formation process (Abreu-Vicente et al. 2015).This overall picture is in agreement with the interpretation that, in general, the PLT phenomenon indicates the dominance of self-gravity during collapse.Indeed, the N -PDF evolution as investigated in simulations of contracting clouds is characterized by the emergence and further growth of a PLT (Ballesteros-Paredes et al. 2011;Kritsuk, Norman & Wagner 2011;Federrath & Klessen 2013;Auddy, Basu & Kudoh 2018;Körtgen, Federrath & Banerjee 2019;Veltchev et al. 2019).If spherical symmetry and isothermality are acceptable approximations for the collapsing region which accounts for the PLT, the slope of the latter is interrelated to the density profile (Donkov, Veltchev & Klessen 2017).
Recently, there is a growing interest to the existence of double PLTs of (column-)density distributions in evolved starforming clouds.First indications of a second PLT of ρ-PDFs emerging at very high densities have been found by Kritsuk, Norman & Wagner (2011) from purely hydrodynamic simulations of supersonic, isothermal and self-gravitating turbulent medium with grid refinement reaching AU scales in the dense cores.A small single PLT (PLT 1) was detected at the formation time of first collapsing objects; after ∼ 0.4 free-fall times it spans already 6 orders of magnitude and a second, shallower PLT (PLT 2) grows out of its high-density end.The authors attributed this PLT 2 to support through centrifugal forces within rotating contracting cores.Further numerical magneto-hydrodynamic research of supersonic and selfgravitating turbulence confirmed the development of a PLT 2 of ρ-PDF after the first PLT had been established (Collins et al. 2011(Collins et al. , 2012)).High-resolution simulations of turbulent collapse by Murray et al. (2017) yield a first PLT with slope corresponding to an attractor solution for the density profile at small scales and a much shallower PLT 2 which is always associated with close vicinities of the forming protostars, within the typical radius of their disks.Thorough analysis of the double-PLT phenomenon in ρ-PDFs of simulated isothermal gravoturbulent fluids was performed by Khullar et al. (2021).Those authors conclude that the density at which the first PLT deviates from the lognormal main part accounts for the transition from unbound to bound gas while the deviation point of the PLT 2 delineates structures in rotationally flattened disks.
In addition, significant theoretical efforts were made to explain the PLT 2. Donkov & Stefanov (2019) model an ensemble of self-gravitating, isothermal turbulent clouds with small homogeneous core, power-law ρ-PDF and steady-state accretion.They obtain two different slopes in the regimes far from the core and near to the core; the latter is steeper and indicative of free fall.In a follow-up study (Donkov et al. 2021), their model was modified through a different treatment of gas thermodynamics in the regime far (isothermality) and near (an equation of state of a "hard polytrope") to the cloud core; in the second case, a polytropic exponent of 4/3 can account for a shallower PLT 2. In contrast to these two works, Jaupart & Chabrier (2020) consider the ρ-PDF as a dynamical entity and develop an analytical theory to describe its evolution in supersonic star-forming cloud as the balance between gravity and turbulence changes.In their framework, the PLT emerges at the time of starting collapse in some regions while a shallower PLT 2 forming at later stages is associated with dense clumps in free-fall collapse.
On the part of observations, there is a growing evidence for double PLTs in N -PDFs in star-forming regions.The first reported cases were three high-mass regions whose N -PDFs display combination of a PLT and a shallower PLT 2 (Schneider et al. 2015a).A more recent study from Herschel imaging of dozens Galactic regions with various mass and star-forming activity shows that the column-density distribution in most of them has a double PLT, with a shallower or a steeper PLT 2, while the main part is best fitted with two lognormals (Schneider et al. 2022).The double PLT phenomenon at the typical spatial scales of Herschel observations is still poorly understood.In clouds with a low mass fraction of dense cores, PLT 1 cannot be attributed (only) to local collapses but rather to gravitational collapse and accretion at larger scales as indicated by a number of observational studies (see references in Schneider et al. 2022, Sect. 4.2).Some suggested explanations of a flatter PLT 2 are feedback effects from young stars, such as ionization compression (Tremblin et al. 2014) or strong magnetic fields in initially subcritical clouds (Auddy, Basu & Kudoh 2018) while a steeper PLT 2 might result from the orientation of magnetic field with respect to the (column-) density distribution (Soler 2019).
The presented paper aims to contribute to the current research of the double PLT phenomenon and its suggested interpretation as signature of rotational support at the scales of protostellar disks.We study the evolution of PLTs both in ρ-PDF and N -PDF from high-resolution simulations of contracting star-forming clumps with varied initial balance of energies which allows to distinguish the effect of rotation on the high-density end of the PDF.The applied technique for multiple PLT extractions does not depend on fitting of the main PDF part, in view of the uncertainty regarding its shape.Details on the used numerical data are provided in Section 2. Section 3 informs the reader of the method to assess the parameters of PLTs.The results on the PLT evolution are presented in Section 4; their implications, in comparison to other works, are discussed in Section 5. Section 6 gives a brief summary of this work.
We use data from a large set of numerical simulations of protostar formation in primordial gas by Wollenberg et al. (2020, hereafter, W20).They have been performed with the Voronoi moving-mesh code Arepo (Springel 2010), adapted with a number of modifications (for details, see W20 and the references therein).The code enables us to adjust the grid resolution according to the growth of density fluctuations and thus is very suitable to investigate collapse in star-forming regions (e.g., Greif et al. 2011Greif et al. , 2012) ) and, in particular, the evolution of high-density parts of the (column-)density distribution with a good enough resolution.
The W20 simulations follow the gas collapse in a box with a sidelength of 13 pc which initially contains a single Bonnor-Ebert (BE) sphere at its centre, with the respective density profile (Bonnor 1956;Ebert 1955).For simplicity, we label this object hereafter as contracting star-forming clump.The linear radius of the sphere is RBE = 1.87 pc, corresponding to the critical dimensionless value for stability, and the mass within it is MBE = 2671 M⊙.Since W20 have explored the formation of Population III stars, the gas metallicity is zero and the initial temperature inside the sphere was set to 200 K which yields its Jeans content: MBE/MJ ≃ 3, where MJ is the critical mass for the stability of self-gravitating isothermal spheres (Jeans 1902).A network of 45 chemical reactions between different species of H and He and free electrons is used to compute cooling and the polytropic index γ (see W20, Sect.2.1).The numerical box outside the BE sphere is filled with homogeneously distributed gas of density ρ(RBE) = 2.7 × 10 −21 g cm −3 .With this setup, the impact of multi-scale flows and externally driven accretion onto the central clump are neglected.
We stress that we are primarily interested in the global collapse problem of a gravitationally unstable structure rather than in the details of the primordial conditions in terms of metallicity, chemical reactions and the resulting cooling1 .We further note that the clouds in W20 are strongly dominated by gravitational contraction.Consequently, it might be expected that the main features of this process do not differ substantially from those in collapsing clumps in the presentday star-forming regions.Hence, a comparison of the W20 simulations with observational data in the local Universe is instructive and appropriate.
Different physical setups for the BE sphere were chosen through varying the initial amounts of rotational and/or turbulent energy: • Pure infall (PI) motions: no turbulence and no rotation.
• Rotation only (RO), no turbulence: we choose between a low (0.01) or high (0.10) value of the ratio β between rotational kinetic and gravitational potential energy.In these cases, the rotation around the z-axis is accounted for by setting a uniform angular velocity.
• Turbulence only (TO), no rotation: here, we select a low (0.05) or high (0.25) value of the ratio α between turbulent and gravitational potential energy.
• Rotation and turbulence (RT): these are simulations Table 1.Runs from the W20 simulations selected to study the PLT evolution.Notation: PI -pure infall, α -the ratio between turbulent and gravitational potential energy, β -the ratio between rotational and gravitational potential energy, t SF -the time of formation of the first sink particle, τ ff ≃ 0.34 Myr -free-fall time of the simulated cloud, t end -run duration after t = t SF , N sinknumber of sink particles at t end and their total mass, M sink .with all possible combinations of the values of α and β specified above.
Five runs (realizations) have been performed for each setup; those without turbulence differed only in the random mesh positions for initialization of the BE density distribution, whereas those with turbulence differed in their random seeds.To avoid artificial fragmentation, a Jeans refinement criterion is applied so that the local Jeans length is always resolved with at least 16 cells.The formation of protostars in the collapsing gas was traced through sink particles.The latter were defined on several standard conditions (see W20, Sect.2.2); for the goals of the present study it is important to highlight the threshold (minimal) density for sink formation n th = 2.4 × 10 15 cm −3 (ρ th ∼ 4 × 10 −9 g cm −3 ) and the accretion radius racc = 2 AU which controls the collapse and the gravitational boundedness of gas in the close vicinity of a sink particle.In fact, racc defines the effective spatial resolution achieved in the densest parts of the collapsing region, and ρ th constitutes the upper density limit for our investigation of PLTs in the ρ-PDF.
For our study, we select a sample of eight W20 runs: two simulations of type PI, four of type RO, as well as one TO and one RT realization.Their main parameters are given in Table 1.The point in time, at which the first sink particle forms is labeled 'star formation time', tSF (Column 3).A single sink particle forms and grows in mass by accretion from the surrounding envelope material in the PI runs, whereas in all other ones there is significant fragmentation, leading to the formation of small clusters of accreting objects (Column 5).The simulation runs up to t = t end ∼ 10 3 yr after tSF (Column 4), when the ionizing feedback from the formed stars can no longer be neglected.At that time the accretion rate onto sink particles effectively drops to zero (cf.Fig. 2 in W20) and the star-formation efficiency (Column 7) approaches a constant value of order of a few percent.

METHOD FOR PLT EXTRACTIONS
For the analysis of the high-density parts of ρ-/N -PDFs from the selected W20 runs, we apply a technique which allows for extraction of multiple PLTs.It is an extension (Marinkova et al. 2021) of the so-called adapted bPlfit method for PLT extraction (Veltchev et al. 2019).Below we summarize the most important features of this technique.Items (i) to (iv) concern the adapted bPlfit method in general and item (v) describes the used technique to extract multiple PLTs.
(i) The adapted bPlfit method is based on the statistical approach bPlfit (Virkar & Clauset 2014) aimed to assess the assumed power-law part2 of an arbitrary binned distribution via Kolmogorov-Smirnov goodness-of-fit statistics.In contrast to other commonly used techniques, the PLT is extracted without any assumptions about the possible functional shape (e.g., lognormal) of the rest of the distribution.This is an advantage in comparison with other approaches, especially when the main part of the PDF suffers from incompleteness and/or has a complex shape, probably affected by multiple competing physical agents.
(ii) The input for the method consists of three parameters: (a) the lower density cutoff controls the range to search for a PLT; (b) the upper density cutoff excludes the poorly resolved high-density tail of the PDF (e.g., the threshold for sink formation ρ th in analysis of ρ-PDFs); (c) range of variation of the total number of bins N bins within the chosen cutoffs, depending on the data volume and on the span of the expected PLT: e.g., 30 ⩽ N bins ⩽ 100.
(iii) The output of the method for fixed N bins are the two parameters of the suggested PLT: deviation point (DP) and slope.This output is not sensitive to the chosen binning scheme: linear, logarithmic, etc.In general, it is also weakly affected by the choice of N bins / bin size -as demonstrated via applications of bPlfit to analytical PDFs composed of a lognormal function and a PLT (Veltchev et al. 2019, Appendix B).Nevertheless, working with PDFs from discrete data sets, some choices of N bins may lead to extraction of unreliable PLTs as the method fits a few bins at the end of the distribution.Such false PLTs are excluded from further consideration by setting some minimal span of the extracted PLT in units of bin size (see Fig. 2 in Veltchev et al. 2019).
(iv) The obtained set of reliable PLTs (after exclusion of false PLTs) is used for averaging the PLT parameters over N bins within the range of variation.
(v) The extraction of multiple PLTs proceeds in the following main steps (Sect.4.2 in Marinkova et al. 2021): (a) The lower cutoff is gradually increased, starting from the lower data limit, and averaged PLT parameters are obtained for each lower cutoff choice (output sample).The upper density cutoff is set at ρ th ; (b) The number of PLTs and their approximate parameters are made by eye; points with high standard deviation from these initial estimates are excluded from the output sample (see Fig. 3 in Marinkova et al. 2021); (c) The DP of a PLT is set as the upper density cutoff of the previous PLT (if there are indications of its existence) and step A is re-done for the latter; (d) The data from the final output sample are used to average the PLT parameters over the groups associated with the presumed PLTs.
In principle, the technique can distinguish between PLTs with very similar slopes.Therefore additional criteria for plausible extractions of multiple PLTs might be imposed depending on the PDF shape, the data resolution and/or the expected values.For instance, one can require a minimal difference in slopes between two extracted PLTs and a minimal density span of each one, in order to consider both as 'true detections'.In Marinkova et al. (2021), these control parameters were set to 0.4 and one order of magnitude, respectively.
In this work, we do not impose any preliminary criteria for plausible PLT extractions, in view of the dynamics and complexity of the PDF shape and of our goal to trace the emergence and the development of a PLT in small time steps in regard to the dynamical time.Additional inspection of the PLT parameters was performed once the results for the whole considered run of W20 had been obtained.Then, at time points when both the differences in slope and DP of two successive PLTs turn out to be less than the uncertainties of these PLT parameters, the PLTs were merged into a single one.For examples, see the evolution of ρ-PDFs for the runs with β = 0.10 in Fig. 5, right.

EVOLUTION OF THE PLTS
In this Section we present the emergence and the evolution of multiple PLTs in the course of the selected W20 runs.Logarithmic density s = ln(ρ/ρ0) and logarithmic column density z = ln(N/N0) are defined as normalization units using the mean density of the initial Bonnor-Ebert sphere ρ0 ≡ ⟨ρ⟩BE = 3MBE/(4πR 3 BE ) = 6.6 × 10 −21 g cm −3 and its corresponding column density N0 ≡ ρ0×13 pc = 0.26 g cm −2 .Multiple (up to three) PLTs of the ρ-PDF and the N -PDF are given by the equations: As,2 exp(q2s), s2 ⩽ s < s3 As,3 exp(q3s), s ⩾ s3 (1) and where (qi, si) and (ni, zi) are the corresponding two PLT parameters (slope and DP) and As,i and Az,i are constants, for 1 ⩽ i ⩽ 3. We point out that the upper (column-)density ends of PLT 1 and PLT 2 parts do not necessarily coincide with the DPs of the next PLT (see the comment in Sect.4.1).
Since the development of PLTs is expected to take place in the course of local collapse(s) we apply our technique to the evolutionary times after the formation of the first sink particle to the end of the run: tSF ⩽ t ⩽ t end .At each time step, after the initial detection of all PLTs through the adapted bPlfit, their PLT parameters were re-evaluated, setting the DP of the high-order PLT as upper density cutoff smaxe.g., taking s3 ≡ smax to obtain new estimates of (q2, s2).This is a small modification in comparison to Marinkova et al. (2021), to avoid possible effects of higher-order PLTs on the parameters of the lower-order ones.

Density PDFs
Examples of ρ-PDFs (for some fixed N bins ) at the onset of star formation and the end of the simulation are plotted in .Density PDFs at t = t SF (dashed) and t = t end (thick grey) in case of pure infall (top) and initial weak rotation support (bottom).Only densities below the sink particle threshold ln(ρ th /ρ 0 ) = 27.17 are considered.The DPs are marked with the corresponding line notation; the slopes at t = t end are plotted.
Fig. 1.The main, non-power-law part of the density distribution deviates substantially from lognormal due to the setup of the simulated clumps, immersed into an initially homogeneous medium (Sect.2).A single PLT 1 spanning several orders of magnitude is already present at the time of first sink particle formation tSF, both in case of pure infall and initial weak rotational support (β = 0.001).At the end of the simulation, there is clear evidence of multiple PLTs.The DP s1 of PLT 1 is approximately the same (cf.the dashed and solid red marks).A steeper PLT 2 is formed and spans more than 3 decades in density in the pure-infall runs (Fig. 1, top).In addition, in the runs with β = 0.001, a third, much shallower PLT (PLT 3) has grown out of PLT 2 at t = t end , reducing significantly its span (Fig. 1, bottom).
It should be noted that our technique for multiple PLT detection tends to extract shorter PLTs in two cases: i) if the PDF contains large bump(s) in the considered density range 3 , 3 Such features appear at the late stages of runs with initial rota-  and/or ii) there is a density range of smooth transition between two PLTs.For instance, after visual inspection of both panels in Fig. 1, one might put by eye s2 at some lower value and estimate q2 accordingly.However, such a fitting function would produce a larger Kolmogorov-Smirnov statistics than the PLT estimate from bPlfit does and hence would be rejected by the method (see Sect. 3 in this work and Sect. 3 in Marinkova et al. 2021).The result would be the same as smin > s1 is being varied -the adapted bPlfit method tends to extract a PLT with parameters which depend weakly on the smin value.Fig. 2 displays the spatial scales which correspond to the extracted PLTs from the studied runs with pure-infall motions as well those with some initial rotational support.The effective sizes are calculated as l eff = (V (si)/Vtot) 1/3 L where L = 13 pc is the sidelength of the numerical box and V (si)/Vtot is the fraction of the total volume defined by the isocontours corresponding to the DPs si (1 ⩽ i ⩽ 3).

Pure infall
The evolution of the double PLT in the pure-infall runs is shown in more details in Fig. 3.It is very similar for the two studied realizations.PLT 1 has a slope q1 ≃ −1.3 at t = tSF which remains unchanged in the course of the run; its DP s1 varies as PLT 2 is emerging and then decreases slowly.A PLT 2 is detected early after the first sink formation.Its slope q2 rapidly steepens to −2 and then stabilizes at this value.In contrast, the DP s2 initially swings within 2 orders of magnitude and seems to tend toward a constant value at t ∼ t end .
In the case of pure infall, a single sink particle forms (cf.Table 1, column 5) in the center of the clump.The spherical symmetry is preserved thereafter and it is instructive to compare the evolution of the radial density profile ρ(l) (where l is a given radius) to that of the ρ-PDF (Appendix A).In case of a power-law density profile ρ ∝ l −p , its exponent is tional and/or turbulent support as dozens of protostars form and their mutual tidal forces cause redistribution of material in their vicinities.p = −3/q (see Donkov, Veltchev & Klessen 2017, and the references therein).Hence the obtained PLT slopes q2 ≃ −2 and q1 ≃ −1.3 in the pure-infall case correspond to p2 = 3/2 and p1 ≃ 2.3, respectively.For reference, in the classical Larson-Penston (LP) models of isothermal spherical collapse (Larson 1969;Penston 1969), the structure of a protostellar object is characterised by profiles p LP, 2 = 3/2 (a core in free fall) and p LP, 1 = 2 (its envelope where the infall motions are retarted due to thermal pressure support4 ).The profile p1 ≃ 2.3 in the outer zone around the sink particle, calculated directly (Fig. A1, right) or assessed from the slope q1 of PLT 1 (Fig. 3, bottom), is steeper than the theoretical one.However, unlike the LP models, the gas in the simulations of W20 is nonisothermal, as illustrated in Fig. 4, where we plot the density PDF together with the median gas temperature.Hence the support of thermal pressure against gravitational contraction is much stronger which leads to a flatter PLT 1.We point out that such an effect is not expected in local star-forming clumps wherein the temperatures are much lower due to the efficient cooling.
To assess the dynamics of the PDF evolution, it is appro- priate to calculate (t end − tSF) in units of the free-fall time of the collapsing clump core.The latter is associated with the volume delineated by isodensity contour ρ = ρ(⟨s2⟩) where the parameters of PLT 2 are averaged in time: ⟨s2⟩ = 8 ln(10) and ⟨q2⟩ = −2 (see Fig. 3).With the mean density within this contour the free-fall time is: Using the latter as an additional unit in Fig. 3 (top xaxis), one sees that both detected PLTs preserve nearly constant values of their parameters within 10 free-fall times of the collapsing clump core, i.e. the two PLTs reflect a quasistationary regime at this stage of the contracting clump evolution.

Local collapses with initial rotational support
In the runs with initial rotational support, dozens of sink particles form, due to emergence of an unstable accretion disk and its subsequent fragmentation (see also Clark et al. 2011;Prole et al. 2023).Their evolution in space, time and number is complex which leads to a permanent, dynamical redistribution of material in the simulated contracting clump.This is clearly noticeable in the evolution of the PLTs shown in Fig. 5.In all realizations with β = 0.01 and β = 0.10 we study, the slope of PLT 1 is the only PLT parameter whose value remains relatively stable in time.It varies within a narrow range around the mean value obtained in the pure-infall case (−1.25 ≳ q1 ≳ −1.5) and appears not to be significantly affected by the location of the DP s1 and/or the evolution of other detected PLTs.
Both PLT 2 and PLT 3 appear early after t = tSF in the realizations with β = 0.01 and are detectable up to t ≲ t end (Fig. 5, left).They are of similar length, with PLT 3 being  much shallower.Although s2 and s3 vary permanently, the minimal length of PLT 2 (about 1 decade) is being preserved.PLT 2 steepens slowly in the course of the run, slower than in the pure-infall case (cf.Fig. 3, bottom), but eventually reaches larger slopes, still comparable to −2 -i.e., the result from the pure infall is roughly reproduced.The slope of PLT 3 varies around −1, with deviations up to 0.3.In terms of scale, PLT 2 and PLT 3 respectively correspond to effective sizes of dozens to a few hundred AU (vicinity/envelope of the emerging protocluster) and ∼ 10 − 20 AU; see Fig. 2. The latter scales are typical for thin accretion disks which can form from 2D turbulence -allegedly one of the mechanisms for transfer of angular momentum outwards (e.g., Dubrulle & Valdettaro 1992;Klahr & Bodenheimer 2003).In agreement with such interpretation, Donkov et al. (2024) showed that a PLT with q ∼ −1 should emerge at the high-density end of the ρ-PDF in self-gravitating, polytropic turbulent clouds -this PLT corresponds to a thin, rotating accretion disk in the close vicinity of a protostellar core.
In the runs with stronger initial rotational support (β = 0.10), a PLT 3 emerges and behaves in a similar way as for β = 0.01.However, here PLT 2 is very unstable and merges with PLT 1 at some points in time (Fig. 5, right).The latter phenomenon can be explained from inspection of the (column-)density and temperature maps at late evolutionary stages of the collapsing clump (Fig. B1 in the Appendix).In the runs with β = 0.10 the emerging protocluster is more dispersed in space and consists of several subgroups with total effective size of the PLT 2 regime -their interaction causes complex dynamics of gas flows and perturbations of the density field.

Local collapses with initial turbulent support
Runs with α > 0 and β = 0 are also characterized by some non-zero initial angular momentum (see W20).For comparison with the results of the RO runs, we select one TO and one TR realization with strong turbulent support, i.e. with  ).In both cases, tSF and N sink are similar to those of the runs with initial rotational support, displaying a bit larger star-formation efficiency (Table 1), whereas the ρ-PDF evolution is more complex.In the TO run α025 − 1, PLT 2 emerges initially with slope |q2| ≲ 2 -like the one from the PI runs but significantly later in time (Fig. 6, left).Its further evolution is seemingly intertwined with the development of a PLT 3 which becomes detectable at times close to the end of the run.At t = t end , the ρ-PDF consists of 3 PLTs, with parameters quite similar to those from the RO runs. Figure 6, right, exemplifies the combined effect of strong initial rotational and turbulent support.The emergence of multiple PLTs occurs at about the same time as in the β010 runs (cf.Fig. 5, right).However, both PLT 2 and PLT 3 are unstable in time, merge and eventually disappear at t ∼ t end .Such a picture is expected in view of intermittent turbulence.The general, the conclusions are that the inclusion of turbulent support slows down the development of multiple PLTs and that the presence/absence of the latter cannot be used as strong indicators of physical processes in turbulent, collapsing regions.
Nevertheless, PLT 3 of slope q3 ∼ −1 appears in all runs with rotational and/or turbulent support.The spatial scale and the density range of this regime corresponds to the sizes of the formed dense accretion disks.Obviously the PLT 3 range is indicative of the role of centrifugal forces in the contracting clump, as suggested originally by Kritsuk, Norman & Wagner (2011).

Column-density PDFs
The N -PDFs in the simulations of W20 also display multiple PLTs.The evolution of the PLT parameters derived from column-density maps obtained from the projection along the rotational Z-axis is illustrated in Fig. 7.The essential difference with the ρ-PDF investigation is that a PLT 3 is not resolved at all.PLT 1 is apparently a permanent feature in the pure-infall case as well in those with initial rotational support.Its slope varies in the range −1.75 ≲ n1 ≲ −1.45.
PLT 2 is detectable shortly after t = tSF in the pure-infall case, with slope which increases towards −4 by the end of the run.In the runs with rotational support, PLT 2 is much flatter; for β = 0.10 it merges with PLT 1 -shortly after its emergence and analogously to its density counterpart (Fig. 7, top right).An interesting issue in the PDF analysis is the relationship between the PLT slopes q and n in the ρ-and N -PDF, respectively.It has been demonstrated in a number of works (see, e.g.Donkov, Veltchev & Klessen 2017, Sect.5.2, and the references therein) that, if the cloud structure can be described by a power-law density profile (i.e. in case of point symmetry in regard to the clump centre), those PLT paratemers should relate as: Since the clumps in this study are self-gravitating and contracting, a power-law density profile is expected, at least in the inner regions, while the presence of rotating disk-like structures should break the spherical symmetry.In Fig. 8 we juxtapose the derived slopes n1 (left) and n2 (right) plotted in Fig. 7 with their counterparts from the ρ-PDF analysis and plot the relation from equation (3) for comparison.In the runs with pure infall, the correlation between q1 and n1 is consistent with the theoretical expectation, with an apparent systematic shift with respect to the slope q LP, 1 = −3/2 (n LP, 1 = −2) corresponding to the asymptotic density profile p = 2 in the envelope of collapsing core in free fall (Larson 1969;Penston 1969).We attribute this phenomenon to the non-isothermal thermodynamics in the W20 simulations commented on in Sect.4.1.1 and leading to shallower slopes q1 and, hence, n1.The results from runs with low rotational support display similar behaviour with some larger scatter while those from the run with β = 0.10 deviate significantly from relation (3) due to the merger of PLT 2 with the PLT 1.
The diagram q2 vs. n2 (Fig. 8, right) displays an even clearer segregation between the pure-infall case and the ones with initial rotational support.The slopes of the second PLTs in the pure-infall case still relate in a good agreement with the expectation from the LP collapse model unlike those (if detected at all) in the runs β = 0.01 and β = 0.10.The latter finding indicates that in case of rotational support the spherical symmetry is broken at the spatial scales which correspond to PLT 2.

Comparison with other works on PLTs
In general, our result on the ρ-PDF in the pure-infall case (Sect.4.1.1)reproduces the solution from the classical LP protostellar collapse model (Larson 1969;Penston 1969): a protostar formed within a core (q2 ≃ −2; p2 ≃ 3/2) with its outer envelope (q1 ≃ −1.3; p1 ≃ 2.3), where p1 > p LP, 1 = 2 due to the non-isothermal thermodynamics and, hence, stronger pressure support.This picture is not surprising in view of the setup of the used W20 simulations: collapsing BE spheres put in a large, homogeneous gas reservoir.The two protostellar 'core + envelope' regimes (corresponding to PLT 2 and PLT 1, respectively) are recognizable even in the runs with initial rotational and/or turbulent support in which dozens of protostars form, with a complex and dynamical distribution in space.
Accretion on the protostar(s) is far from steady state and its total rate varies significantly in all studied runs (see Sect. 4.3 in W20).This clearly affects the structure and evolution of the PDF.Another and more appropriate model for comparison seems to be the one of Murray & Chang (2015) which is aimed to describe self-gravitating clumps wherein massive star formation takes place.Taking into account the driving of turbulence in the compressed accretion flows, those authors find an 'attractor solution' at small scales where the gravity of the protostar dominates: a time-independent density profile with p = 3/2, irrespectively of the type of support.At larger scales, the profile exponent should depend on position: 1.6 ≲ p(ρ) ≲ 1.8.This picture is supported by simulations of Murray et al. (2017).In terms of ρ-PDF (see their Fig.10), they found a PLT with slope ∼ −1.7 which seems to be a combination of PLT 1 and a steeper PLT 2 and another, very shallow one at the high-density end, with slope > −1.This general ρ-PDF shape becomes clearer when only a sphere of 1 pc radius around the star particle is considered and is approximately consistent with our results.
The effects of turbulence on the ρ-PDF evolution in isothermal gravoturbulent fluids were recently studied by Khullar et al. (2021), based on simulations with systematic variation of the Mach number and the virial parameter.The main PDF part is initially lognormal due to the driven turbulence.Later on two PLTs develop which reflects the increasing role of gravity and rotation -the DP of the first one marks the transition from unbound to bound gas and that of the second one delineates disk structures in the collapsing region.In terms of slope, the two PLTs detected in Khullar et al. (2021) remain relatively stable in the course of runs with stronger self-gravity and generally correspond to PLT 2 and PLT 3 extracted in this study: q2 ∼ −2 and q3 ∼ −1.We note, however, that the density span of PLT 2 in their work is much larger which can be attributed to the strong turbulent support against collapse.

Comparison with observations
Although many recent works discuss the features of N -PDFs derived from observations of star-forming regions, only a few are dedicated to the multiple-PLT phenomenon.The study of Schneider et al. (2022) based on Herschel maps is the only one (to our knowledge) which provides a large statisticsthose authors found evidence for double PLTs in N -PDFs of 28 Galactic regions of various mass, size and star-forming activity: from giant molecular clouds with signatures of highmass formation to molecular clouds wherein low-mass stars form (LM regions) and diffuse regions without star formation.The objects from their sample which might be relevant to compare with the simulated contracting clumps from this paper in terms of mass and size are the LM regions (cf.Table 1 in Schneider et al. 2022).The latter are also well resolved at the scales of starless/prestellar cores (< 0.1 pc).On the other hand, in view of the very different chemistry, temperature regime and timescales of star formation, the comparison presented below should be considered with caution.
We apply the adapted bPlfit technique to the N -PDFs of 10 LM regions to extract PLTs.In many cases, the derived PLT parameters differ significantly from those of Schneider et al. (2022) who fitted the N -PDFs by use of sets of models which combine lognormal and power-law functions and determined the best fitting model through the Bayesian information criterion.This is not surprising in view of the essential difference between the methods applied; see Appendix C for an illustration and comments.
In Fig. 9 the derived slopes are plotted vs. the spatial scales defined by the given PLT range of column density.The effective sizes of the scales are calculated analogously to those shown in Fig. 2: l eff = (S(PLTi)/π) 1/2 where S(PLTi) are the areas5 enclosed by the isocontours corresponding to the DPs of PLT 1 and PLT 2 (1 ⩽ i ⩽ 2).Single PLTs of similar slope are detected in three regions (open circles).The PLT 1 and PLT 2 regimes seem to be clearly separated in terms of spatial scale: the former are comparable to the size of the whole star-forming region whereas the steeper6 PLT 2 is associated with clumps of typical size a few tenths of pc.The median slope values of PLT 1 |n 1, med | = −2.13(red dotted line) and PLT 2 |n 2, med | = −3.2(blue dotted line) are much closer than those found from our PI runs but the scatter is too large to allow for a clear conclusion.
To sum up, the shape of the N -PDFs from this work is similar to the well resolved high-density part of the ones in all but one LM star-forming regions studied by Schneider et al. (2022).

SUMMARY
This work presents a numerical study of the evolution of power-law tails (PLTs) in (column-)density probability distribution functions (ρ-PDF and N -PDF).We make use of highresolution simulations of contracting star-forming clumps in primordial gas, with varied initial (rotational, turbulent) support against gravity (Wollenberg et al. 2020), and of the adapted bPlfit method for PLT extraction (Veltchev et al. 2019;Marinkova et al. 2021) which assesses the assumed PLT without any assumptions about the possible functional shape of the rest of the PDF.The clumps are strongly dominated by self-gravity and we are interested in those density ranges where it plays the major role.Therefore the main features of the collapse phase should not differ substantially from those in the present-day star-forming regions.
The primary results can be summarized as follows: • In all considered runs (pure infall, initial rotational support and/or initial turbulent support) multiple PLTs emerge shortly after the formation of the first protostar.The first PLT (PLT 1) in ρ-PDF is a stable feature with slope q1 ≃ −1.3 which corresponds to the outer envelope of the protostellar object (with density profile ρ ∝ l −2 ) in the classical Larson-Penston (LP) protostellar collapse model (Larson 1969;Penston 1969) if the assumption of spherical symmetry is valid.The shallower slope is explained with the nonisothermal thermodynamics which provides additional support against collapse.
• The second PLT (PLT 2) in ρ-PDF is a stable feature in the pure-infall runs (a single protostar formed) but fluctuates significantly in the runs with initial motions (rotation/turbulence) as dozens of protostars form.Its mean slope ⟨q2⟩ ≃ −2 corresponds to a density profile ρ ∝ l −3/2 which describes a core in free fall, in the LP collapse model, or an attractor solution at scales with domination of protostar's gravity, in the model of self-gravitating clumps of Murray & Chang (2015).
• PLT 1 and PLT 2 in the N -PDFs from the pure-infall runs and those with initial rotational support are generally consistent with the estimates obtained through the adapted bPlfit method from Herschel data for a number of Galactic low-mass star-forming regions (Schneider et al. 2022).
• In the runs with initial rotation and/or turbulence a third PLT (PLT 3) in ρ-PDFs appears simultaneously with or after the emergence of PLT 2. It is very shallow and is associated with formation of protostellar accretion disks.Its mean slope ⟨q3⟩ ≃ −1 agrees with the results from numerical studies of supersonic, isothermal, self-gravitating turbulence (Kritsuk, Norman & Wagner 2011;Khullar et al. 2021) .Evolution of the density profile from a pure-infall run (left) and correspondence between its PL exponents p and the PLT slopes q of the ρ-PDF (right).Vertical marks in the left panel denote the deviation points of the deviation points of the power-law parts for the selected evolutionary times (in units t − t SF ).Dashed line in the right is the expected relation q = −3/p; see text.
result from this approach is illustrated in Fig. C1, bottom.Note that no gaps between the spans of PLT 1 and PLT 2 are allowed.In contrast, the adapted bPlfit method rests on the assumption of a single PLT only and is searching for a possible PLT 2 (see Sect. 3).The column-density ranges of PLT 1 and PLT 2 can be separated due to N -PDF features (e.g., large bumps) which are poorly fitted by a power-law (Fig. C1, top).Both methods have their advantages in analysing the N -PDFs in star-forming regions.The approach of Schneider et al. (2022) is based on the pertinent physical assumption that the main part of the PDF is shaped primarily by isothermal, non-gravitating, supersonic turbulent gas which is characterized by a lognormal density distribution (see references in Sect.1).In some regions, this main PDF part can be fitted by multiple lognormals which correspond to some characteristic spatial scales in zones of diffuse gas (Stanchev et al. 2015).
On the other hand, the adapted bPlfit method can capture PLTs which are separated in terms of density range due to -for instance -projection effects and/or the evolutionary stage of collapsing clouds.
Figure1.Density PDFs at t = t SF (dashed) and t = t end (thick grey) in case of pure infall (top) and initial weak rotation support (bottom).Only densities below the sink particle threshold ln(ρ th /ρ 0 ) = 27.17 are considered.The DPs are marked with the corresponding line notation; the slopes at t = t end are plotted.

Figure 4 .
Figure 4. PDF of the density and the median temperature as a function of density from a pure-infall run.The averaged DPs of PLT 1 and PLT 2 are shown.The change in slope in the density PDF coincides with a strong increase in gas temperature.

Figure 5 .
Figure5.The same as Fig.3, but for two rotation-only runs: β001 − 4 (left) and β010 − 2 (right).Vertical dashed lines denote time spans when the PLT 1 and the PLT 2 merge into a single PLT.Other notations are the same like in Fig.3.

Figure 7 .
Figure 7. Evolution of PLTs of the N -PDF (projection along the Z axis) in the pure-infall case (left) and in cases with initial rotational support (right): β = 0.01 (dashed, open symbols) and β = 0.10 (solid, filled symbols).A single realization is chosen for each case.

Figure 8 .
Figure8.Correspondence between the slopes of PLT 1 (left) and PLT 2 (right) from ρ-PDFs and their N -PDF counterparts from maps averaged over the three main directions.The expected relation from equation 3 is plotted for reference (dashed).The black square in the left panel shows the expectation from the LP model of collapsing core; see text.

Figure 9 .
Figure 9. Slopes of PLT 1 (red) and PLT 2 (blue) of N -PDFs, derived from Herschel maps of low-mass star-forming regions (Schneider et al. 2022), plotted against the spatial scale.Open circles denote detections of a single PLT; median values |n 1, med | and |n 2, med | are shown with dotted lines.The ranges of slope variation in the pure-infall case (solid line) and in those with rotational support (dashed line) from this work (Sect.4.2) are given for comparison. .