Density distribution function of a self-gravitating isothermal turbulent fluid in the context of molecular cloud ensembles – III. Virial analysis

ABSTRACT

In this work, we apply virial analysis to the model of self-gravitating turbulent cloud ensembles introduced by Donkov & Stefanov in two previous papers, clarifying some aspects of turbulence and extending the model to account not only for supersonic flows but for trans- and subsonic ones as well. Making use of the Eulerian virial theorem at an arbitrary scale, far from the cloud core, we derive an equation for the density profile and solve it in approximate way. The result confirms the solution ϱ(ℓ) = ℓ⁻² found in the previous papers. This solution corresponds to three possible configurations for the energy balance. For trans- or subsonic flows, we obtain a balance between the gravitational and thermal energy (Case 1) or between the gravitational, turbulent, and thermal energies (Case 2) while for supersonic flows, the possible balance is between the gravitational and turbulent energy (Case 3). In Cases 1 and 2, the energy of the fluid element can be negative or zero; thus the solution is dynamically stable and shall be long lived. In Case 3, the energy of the fluid element is positive or zero, i.e. the solution is unstable or at best marginally bound. At scales near the core, one cannot neglect the second derivative of the moment of inertia of the gas, which prevents derivation of an analytic equation for the density profile. However, we obtain that gas near the core is not virialized and its state is marginally bound since the energy of the fluid element vanishes.

methods: analytical, methods: statistical, ISM: clouds, ISM: structure

1 INTRODUCTION

Molecular clouds (MCs), the birthplaces of stars, are described by complex physics. Turbulence, gravity, thermodynamics, and magnetic fields are regarded to be the main physical agents that determine the structure and the evolution of MCs from their formation in the interstellar medium to their destruction through feedback by the newborn stars (Elmegreen & Scalo 2004; Hennebelle & Falgarone 2012; Klessen & Glover 2016). This complexity can be investigated by various statistical tools that contain various levels of information about the cloud structure and the physics behind it. One appropriate tool is the probability density function (PDF) of column (N-PDF) or mass (ρ-PDF) density. The physics and structure of an MC are imprinted, although with some limitations and degeneracies in its PDF.

Also, the PDF shape may indicate the evolutionary stage of the cloud (Elmegreen & Scalo 2004; Vázquez-Semadeni 2010; Hennebelle & Falgarone 2012; Klessen & Glover 2016). At the largest scales and at early evolutionary stages of MCs, supersonic isothermal turbulence dominates their physics that produces a lognormal PDF (Vázquez-Semadeni 1994; Federrath et al. 2010; Molina et al. 2012). If the cloud mass exceeds its Jeans mass, gravitational collapse starts at larger scales, accompanied by ongoing accretion from the cloud surroundings. This process creates a mass and energy cascade from larger to smaller scales and matter is transferred to small substructures originally formed by supersonic turbulence: clumps and cores. The latter contract in turn as their masses and densities grow (Vázquez-Semadeni et al. 2019, and the references therein). At this stage power-law tail(s) (PLTs) emerge and develop at the high-density end of the initially lognormal ρ-PDF (Kritsuk, Norman & Wagner 2011; Girichidis et al. 2014; Jaupart & Chabrier 2020). The PLT regime corresponds to the range of densities of the clumps/cores where star formation takes place.

The goal of this paper is to investigate the relation between cloud physics and its structure imprinted in the ρ-PDF. We focus on advanced stages of MC evolution that correspond to a pronounced PLT of the ρ-PDF. In Donkov & Stefanov (2018, hereafter, Paper I) and Donkov & Stefanov (2019, hereafter, Paper II), we elaborated a model based on the notion of average representative object. It is a model of spherically symmetric and isotropic cloud that is statistically averaged over the so-called MCs ensemble – the set of all isothermal MCs that have the same ρ-PDF, the same temperature, the same size, and the same densities at the outer edge and in the centre. The average representative object is the simplest abstract member of the ensemble, with effective size derived from the volume weighted ρ-PDF in a natural way (see equation 1 in Section 2). This object retains the main structural and physical characteristics of all members of the ensemble (in terms of abstract scales) while the specific morphology and dynamics of each individual are lost. In regard to the dynamics of the average representative object, we assume that the accretion flow is radially symmetric and stationary, i.e. a steady state for the motion of the fluid elements. It has been shown that the steady state is indeed the stable solution for isothermal gas flows accreting radially on to a central core (Keto 2020).

These simplifications enable us to derive an equation for the ρ-PDF from the governing equations of the medium: the Euler system (the equation of motion and the continuity equation), the equation of state of the gas and the Poisson equation for the gravitational potential. In addition, the assumption of steady state is consistent with the approximately constant slopes of the PLTs observed at advanced stages of cloud evolution (Kritsuk et al. 2011; Girichidis et al. 2014; Veltchev et al. 2019; Khullar et al. 2021).

The equation of conservation of the total energy, as the sum of the kinetic, thermal, and gravitational energy, of the fluid element per unit mass was obtained in Papers I and II from the equations of the medium. Furthermore, through substitution with the explicit form of the energies, one gets a non-linear integral equation for the dimensionless density profile ϱ(ℓ) (Paper II). We obtained approximate solutions of this equation in two physical regimes: far and near to the cloud centre. The first solution, characterized by density profile ϱ(ℓ) = ℓ⁻², corresponds to balance between the accretion (or, accretion plus turbulent) kinetic energy and the gravitational energy of the cloud’s shells located between the fluid element and the central dense region of the cloud, which we labelled ‘core’. This corresponds to a steady-state solution, with the fluid element moving inwards with a constant velocity. The second solution corresponds to a free-fall of the fluid element for which the gravitational energy of the core and the accretion kinetic energy balance each other. In this case, one gets for the density profile ϱ(ℓ) = ℓ^−3/2; this has similarities to the (singular) isothermal sphere solutions (see e.g. Whitworth & Summers 1985).

This work is focused on the virial analysis of the suggested model as the latter is extended to account not only for supersonic (turbulent and/or accretion) flows (like in Papers I and II) but for trans- and subsonic ones as well. Applying the Eulerian virial theorem (EVT; McKee & Zweibel 1992; Ballesteros-Paredes 2006; Krumholz 2015) to arbitrary radius (scale) in the regime far from the core, under the model assumptions for spherically symmetric cloud, we derive an equation for energy balance per unit mass at this scale. The obtained solution of this equation (equation 22) confirms the profile form ϱ(ℓ) = ℓ⁻² in the considered regime and enables us to assess the energy per unit mass of the fluid element, which remained indeterminate in Papers I and II. In regard to the energy balance, we find three possible cases depending on the scaling exponent of the turbulent velocity and on the range of flow velocities in regard to the sonic speed. In two cases, the energy of the fluid element can be negative which means that the solution ϱ(ℓ) = ℓ⁻² corresponds to a stable state in the outer cloud shells and suggests that these cases must be confirmed from observations. On the other hand, one cannot arrive at an equation for the density profile in the near-to-the-core regime, at least because the second derivative of the moment of inertia of the fluid is not negligible. Qualitative considerations show that the gas around the core is not virialized and it is marginally bound.

The structure of the paper is as follows. Our model is briefly recalled in Section 2. In Section 3, we derive the equation for the density profile and obtain an approximate solution. This is the core of our analysis and is presented in four steps: formulation of the EVT (Section 3.1), discussion on its applicability to the suggested model and derivation of the equation for density profile in the far-from-the-core regime in general (Section 3.2) and explicit form (Section 3.3), and obtaining its approximate solution (Section 3.4). Section 4 contains a discussion on our model and the main results. Finally, in Section 5 we present the conclusions and summarize.

In Table 1, we provide a list of frequently used symbols and notations in the presented model (Papers I and II; this work).

Table 1.

Frequently used notations in the model.

Quantity	Description
lc	Total size of the modelled cloud
ρc	Density at the outer cloud boundary
Mc	Mass of the cloud
|$M_{\rm c}^{\ast }$|	Normalizing coefficient (≲ Mc)
l0	Size of the cloud core
ρ0, s0	Density and logdensity of the cloud core
ℓ	Dimensionless scale (normalized to lc)
ϱ	Dimensionless density (normalized to ρc)
p	Exponent of the power-law density profile
q	Slope of the ρ-PDF
cs	Sound speed in the cloud
ut	Turbulent velocity
ua	Accretion velocity
u0	Normalizing factor of the turbulent velocity fluctuations in the standard scaling law
β	Scaling exponent of the turbulent velocity fluctuations
vt	Dimensionless turbulent velocity (normalized to cs)
va	Dimensionless accretion velocity (normalized to cs)
τ	Total kinetic plus thermal energy of the fluid
w	Gravitational energy of the cloud
Ratios of energies (per unit mass) of the fluid element
E0	Total energy of a fluid element moving through the cloud scales, versus thermal energy
G0	Gravitational energy of the entire cloud (excluding the core) versus thermal energy
G1	Gravitational energy of the core versus thermal energy, calculated at the outer cloud boundary
T0	Turbulent kinetic energy at the outer cloud boundary versus thermal energy
A0	Accretion kinetic energy at the outer cloud boundary versus thermal energy
〈G〉, 〈T〉, 〈A〉	Values of G0, T0, and A0, respectively, averaged over the PDF of the abstract object.

2 ON THE MODEL

The model presented in Papers I and II is a simplified description of real MCs but, on the other hand, it is also a generalization, as an average representative of a class of clouds with the same PDFs, sizes, boundary densities, and temperatures (see Donkov, Veltchev & Klessen 2017, and Paper I for more details). It is not aimed at a detailed description of MC morphology and dynamics but, rather, at reproduction of the main physical properties typical for an entire cloud ensemble. We briefly recall it focusing on some aspects of turbulence that need further clarification.

The average representative of the class is modelled by a spherically symmetric, isotropic ball of isothermal (Ferriere 2001) and self-gravitating gas with super-, trans-, or subsonic turbulent flows therein.¹ By use of the PDF, one can introduce dimensionless radii of the cloud, which we consider also as spatial scales:

$$\begin{eqnarray} \ell (s)=\Bigg (\int \limits _{s}^{s_{0}}p_v(s^{\prime })\mathrm{d}s^{\prime }\Bigg)^{1/3}~. \end{eqnarray}$$

(1)

Here, p_v(s) is the volume-weighted PDF, s ≡ ln (ρ(ℓ)/ρ_c) ≡ ln (ϱ(ℓ)) is the log-density, where ρ_c is the density at the outer boundary of the cloud, ϱ ≡ ρ/ρ_c is the dimensionless density. The relation between the mass density profile ρ(ℓ) and the PDF² in differential form is

$$\begin{eqnarray} p_v(s)\mathrm{d}s=-3 \ell ^2 \mathrm{d}\ell . \end{eqnarray}$$

(2)

The dimensionless radius (scale) ℓ = l/l_c varies in the range ℓ₀ < ℓ ≤ 1, where l is an arbitrary physical scale and l_c is the physical scale of the whole cloud. The lower boundary, ℓ₀ ≪ 1, corresponds to the radius of a homogeneous small dense core of density ρ₀ and log-density s₀ = ln (ϱ₀) = ln (ρ₀/ρ_c) located at the centre of the modelled cloud.

The model includes also (super-, trans-, or subsonic) accretion of matter, which flows from the outer cloud boundary down through all scales to the central core. In addition to the inflow velocity, we also consider turbulent velocity fluctuations that are homogeneous and isotropic at each scale ℓ in the interval ℓ₀ < ℓ ≤ 1 presumed to lie within the inertial range of the turbulent cascade where dissipation can be neglected. Therefore, it is convenient to use velocity decomposition |$\vec{u}=\vec{u}_{\rm a}+\vec{u}_{\rm t}$| with |$\vec{u}_{\rm a}=-u_{\rm a}\hat{e}_r$| (only infall) and |$\vec{u}_{\rm t}$| being isotropic. The entire cloud is at steady state at each scale both in regard to the macro-states (the motion of the fluid elements) and the micro-states (the thermal motion of the molecules). We consider it as averaged over an ensemble of spherical objects (clouds) that have the same PDFs (also sizes, boundary densities, and temperatures) and dynamics as described above. Due to the assumption of steady state, one can apply the ergodic hypothesis and hence substitute the averaging with regard to time for averaging over the ensemble of copies of our system that describes all dynamical states of the fluid elements. Let us then consider, at a fixed moment in time, a fluid element that resides in a given fluid volume of the cloud, at scale ℓ. The averaged velocity field is equal to the infall velocity |$\vec{u}_{\rm a}=-u_{\rm a}\hat{e}_r$|⁠, because the averaged turbulent velocity vanishes due to the assumption of isotropy of the turbulent motions.³ The density field must be averaged too, which yields a corresponding density profile – if the PDF is of power-law type (i.e. monotonically decreasing function of density).

The model was further developed in Paper II to account for the gravity of matter external in regard to the cloud. In the spirit of our model, we assume that the material outside the cloud boundary obeys a radial symmetry and gives rise to a potential φ^ext in the cloud volume, where it is constant and does not contribute to the gravitational force acting on the fluid elements.

3 DERIVATION OF THE EQUATIONS FOR THE MASS DENSITY PROFILE FROM EULERIAN VIRIAL THEOREM

3.1 Eulerian virial theorem

Let us briefly consider the EVT, which is the virial theorem for a fixed volume V, surrounded by a surface S. We focus on the simplified scalar form of the tensorial virial theorem that follows from the Boltzmann equation. For the derivation of this form and comments on the physical nature of its terms, we refer the reader to McKee & Zweibel (1992), Ballesteros-Paredes (2006), and Krumholz (2015).

The use of EVT, instead of Lagrangian virial theorem (LVT), is essential for the treatment in this work because only the EVT allows for the description of a non-equilibrium but stationary system. Since we consider a steady-state gas flow through mass shells at a given radius (and not the inward motion of the shells themselves), we naturally operate from Eulerian perspective. The EVT refers to fixed volume elements in space, through which the gas can flow, rather than to fixed mass elements that necessarily change their location as the gas flows. Moreover, in case the second time derivative of the moment of inertia is zero, the system in the LVT treatment is necessarily in equilibrium, or in a ballistic trajectory at most, with zero net forces (see the discussion in Ballesteros-Paredes 2006), whereas in the EVT treatment it is either in static equilibrium without any inside flows (e.g. gravity is balanced by gas pressure), or in steady state. Thus, the application of the EVT to a cloud in steady state allows for obtaining profiles of density and infall speed that depend only on radius, but not on time.

The mathematical expression of the EVT is

$$\begin{eqnarray} \frac{1}{2} \ddot{I} = 2\big (\tau -\tau _{S} \big) + w - \frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\int \limits _{S}\rho r^2 \vec{u} \cdot \vec{\mathrm{d}S}~, \end{eqnarray}$$

(3)

where |$\vec{r}$| is a radius vector (with absolute value |$r=\vert \vec{r}\vert$|⁠) from the origin of the Eulerian coordinate system, |$\vec{u}(\vec{r},t)$| is the velocity field (with absolute value |$u=\vert \vec{u}\vert$|⁠) in the fluid, and |$\rho (\vec{r},t)$| is the density field. The terms in equation (3) have the explicit form and physical meaning as follows:

$$\begin{eqnarray} I= \int \limits _{V} \rho r^2 \mathrm{d}V~, \end{eqnarray}$$

(4)

is the moment of inertia of the fluid in the volume V, calculated about the origin, and |$\ddot{I}$| is its second time derivative,

$$\begin{eqnarray} \tau = \int \limits _{V} \bigg (\frac{1}{2} \rho u^2 + \frac{3}{2} P_{\rm th} \bigg) \mathrm{d}V~, \end{eqnarray}$$

(5)

is the total kinetic plus thermal energy of the fluid (P_th is the thermal pressure),

$$\begin{eqnarray} \tau _{S}= \frac{1}{2} \int \limits _{S} \vec{r}\cdot \Pi \cdot \vec{\mathrm{d}S}~, \end{eqnarray}$$

(6)

is the confining pressure on the volume surface, including both the thermal pressure and the ram pressure of any gas flowing across the surface, where Π_ij = ρu_iu_j + P_thδ_ij is the general form of the fluid pressure tensor,

$$\begin{eqnarray} w= - \int \limits _{V} \rho \vec{r} \cdot \vec{\nabla } \phi \, \mathrm{d}V~, \end{eqnarray}$$

(7)

is the gravitational energy of the cloud, where the potential ϕ accounts both for self-gravity of the cloud and gravity of the external matter. The last term on the r.h.s. in equation (3) represents the rate of change of the momentum flux across the cloud surface.

We note that we neglect the magnetic fields in our model (see Paper I, section 3.1), in order to simplify the first step to obtain the density profile of the fluid from first principles.

3.2 Applicability of the EVT to our model: General consideration

The EVT is a powerful tool as it describes the global energy balance in the objects/systems to which it is applied. We use the EVT for obtaining information about the energy balance for the entire cloud as well as for an arbitrary scale ℓ in its interior. Moreover, the modelled cloud is spherically symmetric and, with the origin of the Eulerian coordinate system set in the centre, the terms of the EVT are easy to calculate – in view of the radial accretion flow and the assumed steady state for macro- and micro-states. We stress that all terms in equation (3) are assumed to be ensemble averaged as described in Paper I (section 3.1) and in Section 2.⁴

There are two points to comment on the application of the EVT to the cloud considered here, i.e. to scales ℓ₀ < ℓ ≤ 1.

First, one can see that the last term in equation (3) becomes zero (for ℓ = 1 as well) since |$\langle \rho r^2 \vec{u} \rangle = {\rm const}(t)$| due to the steady-state condition.

Second, let us consider the term |$\langle \ddot{I} \rangle$|⁠. The model assumes that – despite the accretion flow through the scales, – no material is being accumulated in the cloud shells; all the gas eventually flows into the cloud core. The latter is presupposed to be very small: of nearly constant size and density within the characteristic accretion time. Therefore, considering a scale ℓ far from the core (ℓ ≫ ℓ₀), one can suppose |$\langle \ddot{I} \rangle \simeq 0$| and equation (3) takes the form

$$\begin{eqnarray} 0 = 2\big (\langle \tau \rangle - \langle \tau _{S} \rangle \big) + \langle w \rangle ~. \end{eqnarray}$$

(8)

On the other hand, in the regime near to the core (ℓ ≳ ℓ₀) one should be cautious and not adopt |$\langle \ddot{I} \rangle \simeq 0$|⁠, regardless of the presuppositions in the model.

3.3 The explicit form of the EVT equation

In this section, we derive the explicit form of equation (8) in the framework of the model (Section 2) and using the explicit form of the terms 〈τ〉 (formula 5), 〈τ_S〉 (formula 6), and 〈w〉 (formula 7).

3.3.1 The explicit form of 〈τ〉

It is convenient to normalize all terms in equation (8) to the thermal energy per unit mass: |$c_{\rm s}^2$|⁠. Introducing a dimensionless velocity field v ≡ u/c_s, the kinetic energy term can be written (see Paper I, sections 3.2 and 3.3) as

$$\begin{eqnarray} \left\langle \varrho v^2 \right\rangle =\varrho \left\langle v^2 \right\rangle =\varrho \left\langle v_{\rm t}^2 \right\rangle + \varrho \left\langle v_{\rm a}^2 \right\rangle ~, \end{eqnarray}$$

(9)

where |$\varrho \langle v_{\rm t}^2 \rangle$| is the kinetic energy per unit volume due to turbulent velocity fluctuations and |$\varrho \langle v_{\rm a}^2 \rangle$| is the kinetic energy per unit volume due to accretion of matter towards the cloud centre.

Since the considered spherically symmetric cloud is ensemble averaged we choose to apply a standard scaling relation for |$\langle v_{\rm t}^2 \rangle$| (see Paper II, section 2.1):

$$\begin{eqnarray} \left\langle v_{\rm t}^2 \right\rangle = \frac{u_0^2}{c_{\rm s}^2}\bigg (\frac{l_{\rm c}}{1~{\rm pc}}\bigg)^{2\beta } \ell ^{2\beta } = T_0 \ell ^{2\beta }~, \end{eqnarray}$$

(10)

where u₀ and 0 ≲ β ≲ 1 are the normalizing factor and the scaling exponent of the turbulent velocity fluctuations in the standard scaling law:

$$\begin{eqnarray} u=u_0 L^{\beta } ~, \end{eqnarray}$$

(11)

(Kolmogorov 1941; Larson 1981; Padoan et al. 2006; Kritsuk et al. 2007; Federrath et al. 2010), and |$T_0\equiv (u_0^2/c_{\rm s}^2)(l_{\rm c}/{1~{\rm pc}})^{2\beta }$| is the ratio of the turbulent kinetic energy per unit mass of the fluid element at the cloud boundary to the thermal energy per unit mass. The explicit form of |$\langle v_{\rm t}^2 \rangle$| presupposes that the turbulence is formally independent of the accretion.

The explicit form of the accretion term can be obtained from the continuity equation:

$$\begin{eqnarray} \left\langle v_{\rm a}^2 \right\rangle = A_0 \varrho (\ell)^{-2} \ell ^{-4}~, \end{eqnarray}$$

(12)

where A₀ = const(ℓ) (see Paper I, section 3.3). Since the quantities are dimensionless and normalized, one can obtain A₀ choosing ℓ = 1, i.e. at the cloud boundary. Then ϱ = 1 and one gets |$\langle v_{\rm a}^2\rangle = u_{\rm a,c}^{2}/c_{\rm s}^{2}$|⁠: the ratio of the accretion kinetic energy term at the cloud boundary to the thermal kinetic energy per unit mass.

Now, taking into account the presuppositions in Section 2 and assuming a density profile for the cloud in the form ϱ(ℓ) = ℓ^−p and the normalized thermal energy per unit volume P_th = ϱ(ℓ) = ℓ^−p, the term 〈τ〉 can be calculated as follows:

$$\begin{eqnarray} &&\langle \tau \rangle = 3M_{\rm c}^{*} c_{\rm s}^2 \int \limits _V \left\langle \frac{1}{2}\varrho v^2 + \frac{3}{2} P_{\rm th} \right\rangle \frac{\mathrm{d}V}{4\pi l_{\rm c}^3} = \\ &&3M_{\rm c}^{*} c_{\rm s}^2 \int \limits _{\ell _0}^{\ell } \bigg (\frac{1}{2} T_0 \ell ^{\prime (-p+2\beta)} + \frac{1}{2} A_0 \ell ^{\prime (p-4)} + \frac{3}{2}\ell ^{\prime (-p)} \bigg) \ell ^{\prime 2} \mathrm{d}\ell ^{\prime } \\ &&+ \langle \tau \rangle _{0} = \\ &&= 3M_{\rm c}^{*} c_{\rm s}^2 \bigg [\frac{1}{2}T_0 \frac{\ell ^{2\beta +3-p} - \ell _0^{2\beta +3-p}}{2\beta +3-p} + \frac{1}{2}A_0 \frac{\ell ^{p-1} - \ell _0^{p-1}}{p-1} + \\ &&+ \frac{3}{2}\frac{\ell ^{3-p} - \ell _0^{3-p}}{3-p}\bigg ] + \langle \tau \rangle _{0}~. \end{eqnarray}$$

(13)

The normalizing coefficient |$M_{\rm c}^{*}=(4/3)\pi \rho _{\rm c} l_{\rm c}^3$| is lower than the cloud mass by a factor of few; it is the mass of a ball with radius equal to the radius of the cloud boundary and with homogeneous density equal to the density at this boundary. The term 〈τ〉₀ denotes the kinetic and thermal energy of the core. Its exact expression requires additional assumptions about the physics of the core – at least, its thermodynamic state and the velocity field in it. We avoid such complication of the model and note only that one should expect that |$\langle \tau \rangle _{0} \sim \ell _0^{\alpha }$|⁠, where the exponent α is positive and of order unity.⁵ Considering the terms in equation (13) that contain explicitly ℓ₀ to some power δ, one can notice that the exponent δ ∼ 1 since observational and numerical studies show that 2 ≳ p ≳ 1 (Vázquez-Semadeni et al. 2019; Veltchev et al. 2019) and 1 ≳ β ≳ 0 (Goodman et al. 1998; Elmegreen & Scalo 2004; Shetty et al. 2012; Klessen & Glover 2016).

In view of the relations ℓ ≫ ℓ₀ and α ∼ δ ∼ 1, all terms in equation (13) containing |$\ell _0^{\delta }$| (0 < δ ∼ 1) and |$\ell _0^{\alpha }$| (0 < α ∼ 1) can be neglected in comparison to the ones containing ℓ^δ (0 < δ ∼ 1). From these considerations, we arrive at a more simple form for 〈τ〉:

$$\begin{eqnarray} \left\langle \tau \right\rangle \simeq 3M_{\rm c}^{*} c_{\rm s}^2 \left [ \frac{1}{2}T_0 \frac{\ell ^{2\beta +3-p}}{2\beta +3-p} + \frac{1}{2}A_0 \frac{\ell ^{p-1}}{p-1} + \frac{3}{2}\frac{\ell ^{3-p}}{3-p} \right ] ~. \end{eqnarray}$$

(14)

3.3.2 The explicit form of 〈τ_S〉

Let us first consider the integrand in equation (6) in dimensionless form, according to the model assumptions:

$$\begin{eqnarray} &&\left\langle \vec{\ell }\cdot \Pi \cdot \vec{\mathrm{d}S} \right\rangle = \\ &&\big [ \varrho \left\langle v_{\rm a}^2 \right\rangle + \varrho \left\langle v_{\rm t}^2 \right\rangle \left\langle \cos \alpha ^2 \right\rangle + 2 \varrho \left\langle v_{\rm a} v_{\rm t} \right\rangle \left\langle \cos \alpha \right\rangle + \varrho \big ] \\ &&\ell ^3 \sin \theta \mathrm{d}\theta \mathrm{d}\varphi ~, \end{eqnarray}$$

(15)

where |$\vec{\mathrm{d}S} = \hat{\ell } \ell ^2 \sin \theta \mathrm{d}\theta \mathrm{d}\varphi$|⁠, |$\hat{\ell }$| is unity vector pointing outwards, θ is the polar angle in the range [0, π], and φ is the azimutal angle in the range [0, 2π]. The last addend in equation (15) ϱ originates from the thermal pressure term. According to the model, the tensor product in Π can be presented as follows: |$\vec{v}\otimes \vec{v} = \vec{v_{\rm a}}\otimes \vec{v_{\rm a}} + \vec{v_{\rm t}}\otimes \vec{v_{\rm t}} + \vec{v_{\rm a}}\otimes \vec{v_{\rm t}} + \vec{v_{\rm t}}\otimes \vec{v_{\rm a}}$|⁠. Hence

$$\begin{eqnarray} \vec{\ell }\cdot \vec{v}\otimes \vec{v}\cdot \hat{\ell } = \ell v_{\rm a}^2 + \ell v_{\rm t}^2 \cos \alpha ^2 + 2 \ell v_{\rm a} v_{\rm t} \cos \alpha ~, \end{eqnarray}$$

where |$\vec{\ell }\cdot \vec{v_{\rm t}} = \ell v_{\rm t} \cos \alpha$|⁠, |$\hat{\ell }\cdot \vec{v_{\rm t}} = v_{\rm t} \cos \alpha$|⁠, and α is an accidental angle that takes values in the range [0, π] with equal probability, due to the assumption of isotropic turbulence. Also, one obtains for the accretion velocity |$\vec{\ell }\cdot \vec{v_{\rm a}} = - \ell v_{\rm a}$| and |$\hat{\ell }\cdot \vec{v_{\rm a}} = - v_{\rm a}$| because of the opposite directions of the vectors |$\vec{\ell }$| and |$\hat{\ell }$|⁠, on the one hand, and of |$\vec{v_{\rm a}}$|⁠, on the other hand.

In view of the statistical nature of α, one sees easily that:

$$\begin{eqnarray} \left\langle \cos \alpha \right\rangle = \frac{1}{\pi } \int _{0}^{\pi } \cos \alpha \mathrm{d}\alpha = 0 ~, \end{eqnarray}$$

$$\begin{eqnarray} \left\langle \cos \alpha ^2 \right\rangle = \frac{1}{\pi } \int _{0}^{\pi } \cos \alpha ^2 \mathrm{d}\alpha = \frac{1}{2} ~. \end{eqnarray}$$

The integration over the sphere with radius ℓ in equation (6) yields simply a factor 4π, because there is no any addend in the parentheses in equation (15) that depends on θ or φ. Thus one arrives at

$$\begin{eqnarray} \left\langle \tau _S \right\rangle \sim 4\pi \bigg (\ell ^3 \varrho \left\langle v_{\rm a}^2 \right\rangle + \frac{1}{2} \ell ^3 \varrho \left\langle v_{\rm t}^2 \right\rangle + \ell ^3\varrho \bigg) ~. \end{eqnarray}$$

Finally, plugging in the density profile of the presumed type ϱ(ℓ) = ℓ^−p and the expressions for |$\langle v_{\rm t}^2 \rangle$| and |$\langle v_{\rm a}^2 \rangle$| from equations (10) and (12), and normalizing in respect to density, velocity, and scale (like in Section 3.3.1), we obtain the explicit form of the surface term in the EVT as follows:

$$\begin{eqnarray} \left\langle \tau _S \right\rangle = 3M_{\rm c}^{*} c_{\rm s}^2 \bigg [ \frac{1}{2} A_0 \ell ^{p-1} + \frac{1}{4} T_0 \ell ^{2\beta +3-p} + \frac{1}{2} \ell ^{3-p} \bigg ] ~. \end{eqnarray}$$

(16)

3.3.3 The explicit form of 〈w〉

The gravitational term 〈w〉 is calculated in the EVT equation as follows. Note that |$-\vec{\nabla } \phi$| in equation (7) is the gravitational force per unit mass acting on a fluid element. Due to the Gauß theorem, in spherical systems only mass inside a sphere with radius ℓ contributes to the force, and we can write

$$\begin{eqnarray} -\vec{\nabla }\phi = c_{\rm s}^2 \left [-\frac{G}{c_{\rm s}^2}\frac{M(\ell)}{l_{\rm c}^2}\frac{\vec{\ell }}{\ell ^3} - \frac{G}{c_{\rm s}^2}\frac{M_0}{l_{\rm c}^2}\frac{\vec{\ell }}{\ell ^3}\right ]~, \end{eqnarray}$$

(17)

where |$\vec{\ell }$| is a dimensionless vector with absolute value ℓ and directed radially outwards, G is the gravitational constant, M₀ is the core mass,

$$\begin{eqnarray} M(\ell)=3M_{\rm c}^{*}\int \limits _{\ell _0}^{\ell } \ell ^{\prime 2}\varrho (\ell ^{\prime })d\ell ^{\prime } = \frac{3}{3-p}M_{\rm c}^{*}\left(\ell ^{3-p}-\ell _0^{3-p}\right) \end{eqnarray}$$

is the mass of the inner shells in respect to ℓ, and |$M_{\rm c}^{*}= (4/3)\pi l_{\rm c}^3 \rho _{\rm c}$| is the normalizing coefficient (see Section 3.3.1). Formula (17) can be rewritten in a more simple form recalling the dimensionless coefficients |$G_0=(2G/c_{\rm s}^2)(M_{\rm c}^{*}/l_{\rm c})$| and |$G_1=(2G/c_{\rm s}^2)(M_0/l_{\rm c})$| introduced in Paper I (section 3.4) and evaluated for the outer cloud boundary l_c. They are interpreted, respectively, as the ratio of the gravitational energy of the entire cloud (excluding the core) to its thermal energy per unit mass (G₀) and the ratio of the gravitational energy per unit mass of the core to the thermal energy per unit mass (G₁).

Now the equation (17) is transformed to

$$\begin{eqnarray} -\vec{\nabla }\phi = \left(c_{\rm s}^2/l_{\rm c} \right) \left [-\frac{1}{2} G_0 \frac{3}{3-p} \left(\ell ^{3-p}-\ell _0^{3-p} \right) \frac{\vec{\ell }}{\ell ^3} - \frac{1}{2} G_1 \frac{\vec{\ell }}{\ell ^3}\right ]~, \end{eqnarray}$$

(18)

and, after the corresponding substitutions in equation (7), one gets

$$\begin{eqnarray} &&\left\langle w \right\rangle = - 3M_{\rm c}^{*} \int \limits _{V} \frac{\rho }{\rho _{\rm c}} \vec{l} \cdot \vec{\nabla }\phi \frac{\mathrm{d}V}{4\pi l_{\rm c}^3} = \\ &&= 3M_{\rm c}^{*} c_{\rm s}^2 \int \limits _{\ell _0}^{\ell } \ell ^{\prime (-p)} \vec{\ell ^{\prime }} \cdot \bigg (-\frac{1}{2} G_0 \frac{3}{3-p} \left(\ell ^{\prime (3-p)}-\ell _0^{3-p} \right) \frac{\vec{\ell ^{\prime }}}{\ell ^{\prime 3}} - \\ &&- \frac{1}{2} G_1 \frac{\vec{\ell ^{\prime }}}{\ell ^{\prime 3}} \bigg) \ell ^{\prime 2} \mathrm{d}\ell ^{\prime } + \\ &&+ 3M_{\rm c}^{*} c_{\rm s}^2 \int \limits _{0}^{\ell _0} \ell _0^{(-p)} \vec{\ell ^{\prime }} \cdot \bigg (-\frac{1}{2} G_0 \frac{3}{3-p} \left(\ell ^{\prime (3-p)}-\ell _0^{3-p}\right) \frac{\vec{\ell ^{\prime }}}{\ell ^{\prime 3}} - \\ &&- \frac{1}{2} G_1 \frac{\vec{\ell ^{\prime }}}{\ell ^{\prime 3}} \bigg) \ell ^{\prime 2} \mathrm{d}\ell ^{\prime }~. \end{eqnarray}$$

(19)

Note that in the second integral, we take into account the assumption that the core is homogeneous and thus the density profile is simply |$\varrho =\ell _0^{-p}={\rm const}(\ell)$|⁠. After some algebra, we arrive at

$$\begin{eqnarray} \left\langle w \right\rangle &=& 3M_{\rm c}^{*} c_{\rm s}^2 \bigg [ - \frac{3}{3-p} \frac{G_0}{2} \bigg (\frac{\ell ^{5-2p}-\ell _0^{5-2p}}{5-2p} - \\ &&\frac{\ell ^{2-p}-\ell _0^{2-p}}{2-p} \ell _0^{3-p} + \frac{p-3}{2(5-p)}\ell _0^{5-2p} \bigg) - \\ &&- \frac{G_1}{2} \bigg (\frac{\ell ^{2-p}-\ell _0^{2-p}}{2-p} + \frac{\ell _0^{2-p}}{2} \bigg) \bigg ]~. \end{eqnarray}$$

(20)

The same reasoning about sub-terms of the kinetic energy term 〈τ〉 containing ℓ^δ or |$\ell _0^{\delta }$| (0 < δ ∼ 1) presented in Section 3.3.1 is applicable for the gravitational energy term 〈w〉 as well. One can neglect the terms containing |$\ell _0^{\delta }$| and retain the leading ones with ℓ^δ because ℓ ≫ ℓ₀. Then equation (20) is simplified to the form

$$\begin{eqnarray} \left\langle w \right\rangle \simeq 3M_{\rm c}^{*} c_{\rm s}^2 \bigg [- \frac{3}{3-p} \frac{G_0}{2} \frac{\ell ^{5-2p}}{5-2p} - \frac{G_1}{2} \frac{\ell ^{2-p}}{2-p} \bigg ] ~. \end{eqnarray}$$

(21)

Eventually, one obtains equation (8) in an explicit normalized form

$$\begin{eqnarray} &&\frac{[2\left\langle \tau \right\rangle - 2\left\langle \tau _S \right\rangle + \left\langle w \right\rangle ]}{3M_{\rm c}^{*} c_{\rm s}^2} \simeq \\ &&\bigg (\frac{1}{2\beta +3-p} - \frac{1}{2} \bigg) T_0 \ell ^{2\beta +3-p} + \bigg (\frac{1}{p-1} - 1 \bigg)A_0 \ell ^{p-1} + \\ &&\bigg (\frac{3}{3-p} - 1 \bigg) \ell ^{3-p} - \frac{3}{3-p} \frac{G_0}{2} \frac{\ell ^{5-2p}}{5-2p} - \frac{G_1}{2} \frac{\ell ^{2-p}}{2-p} \simeq 0 ~. \end{eqnarray}$$

(22)

3.4 Solutions of equation (22)

To obtain possible solutions of equation (22), we use the approach of comparing the exponents of its terms. The method was already used in Papers I and II and explained in details in section 3.1 of Paper II. Here, we briefly recall the main idea, for clarity.

After the approximations discussed in the previous section, the EVT equation (22) contains only terms with power-law dependence on ℓ. Since 0 < ℓ ≤ 1, the terms with lower powers dominate over the higher ones. The idea of the approach is to retain only the dominant terms (of the lowest power) and to neglect the remaining ones. The dominant terms must be at least two: a single dominant term would remain unbalanced and then the only solution would be the trivial one. The powers of ℓ factor out and only constants remain. At the end, one obtains an equation that gives the balance between the different energy terms.

To identify the leading order terms and to obtain an equation for p, one chooses a pair of terms with the hypothesis that they are dominant and have equal powers. The equality of the powers results in a simple equation for p. Adopting the obtained root for p, the values of all exponents in the equation are evaluated and the validity of the considered hypothesis is confirmed or rejected. The same recipe is applied to all possible pairs of terms (Zhivkov 1999; Riley, Hobson & Bence 2006).

In our model, the exponents of the terms in equation (22) are: p − 1, 3 − p, 2β + 3 − p, 5 − 2p, and 2 − p. The roots for p obtained by the algorithm described above are: 2, 3/2, 3, 2 + β, 2 − 2β, respectively (see Table 2), and the values of the exponents for each of them are given in Table 3. We remind the reader that 0 ≲ β ≲ 1 (Section 3.3.1).

Table 2.

Possible values of the density profile exponents p obtained from comparison between the exponents of the terms in the equation (22).

Exponents	p − 1	3 − p	2β + 3 − p	5 − 2p	2 − p
p − 1	–	2	2 + β	2	3/2
3 − p	2	–	–	2	–
2β + 3 − p	2 + β	–	–	2 − 2β	–
5 − 2p	2	2	2 − 2β	–	3
2 − p	3/2	–	–	3	–

Table 3.

Calculated exponents of the terms in equation (22), corresponding to the roots for the density profile p from Table 2.

p \ exp.	p − 1	3 − p	2β + 3 − p	5 − 2p	2 − p
2	1	1	2β + 1	1	0
3/2	1/2	3/2	2β + 3/2	2	1/2
3	2	0	2β	−1	−1
2 + β	β + 1	1 − β	β + 1	1 − 2β	−β
2 − 2β	1 − 2β	1 + 2β	1 + 4β	1 + 4β	2β

As discussed in Section 3.2, there are two different physical regimes in the framework of our cloud model: far from the core and near to the core. We explore the solutions for the density profile in the former case and comment on the difficulty to obtain such in the regime near to the core.

3.4.1 The regime far from the core

Here, one neglects the last term in equation (22), i.e. the core contribution in the gravitational term. This reduces the set of exponents to p − 1, 3 − p, 2β + 3 − p, and 5 − 2p and yields roots for p: 2, 2 + β and 2 − 2β. Now there are two possibilities for the energy balance, depending on the turbulent velocity scaling:

• If β > 0, there is a single solution for p = 2 and the energy balance is in the form

  $$\begin{eqnarray} 4 - 3G_0 \approx 0 ~. \end{eqnarray}$$

  (23)
• If β = 0, there are three equal roots, corresponding to the single solution p = 2, whereas the energy balance is different

  $$\begin{eqnarray} T_0 + 4 - 3G_0 \approx 0 ~. \end{eqnarray}$$

  (24)

In case of a positive turbulent velocity scaling exponent (i) gravity is balanced by thermal energy; the accretion and turbulent terms do not play any role. Apparently such a balance is possible only if the accretion and turbulent velocities are of the order of the sound speed or less. This subcase yields G₀ ≈ 4/3. The result for the energy balance of the fluid element per unit mass from Paper II is not applicable here, since it was obtained only for supersonic (accretion and/or turbulent) flows. Hence, one is not able to obtain the energy E₀ of fluid element.

On the contrary, if the turbulent velocity does not depend on the spatial scale (ii), there exist two further possibilities. If the accretion and turbulent velocities are of the order of the sound speed or less, the energy balance reads: T₀ + 4 − 3G₀ ≈ 0; then the results from Paper II are inapplicable again (the energy E₀ of the fluid element remains undetermined). But if the turbulent flow is supersonic, the balance is: T₀ − 3G₀ ≈ 0, and the result from Paper II is valid: A₀ + T₀ − 3G₀ ≈ E₀. Combining these two equations, one obtains that E₀ ≈ A₀ > 0. The instance E₀ ≈ 0 is also possible – if the turbulent flow is supersonic while the accretion one is trans- or subsonic; then the accretion term does not play a role in the equation (see Paper II, section 4, fourth paragraph). Hence, the energy of fluid element is positive (an unstable energy state), if both the turbulence and accretion are supersonic, and is roughly zero (a marginally bound energy state), if only the turbulence is supersonic.

Our model allows to write down the energy-balance equations as the different forms of energy are averaged over the PDF of the abstract object. Let us recall the definition of G₀ and of the average density of the cloud which in the case of a long PL-tail can be evaluated from 〈ρ〉_c = ρ_cq/(1 + q) (Donkov et al. 2017), where q is the PDF slope and q = −3/p. Then in the case q = −3/2 (p = 2), the above formula gives 〈ρ〉_c = 3ρ_c, and, since |$G_0=(2G/c_{\rm s}^2)(M_{\rm c}^{*}/l_{\rm c})\sim \rho _{\rm c}$|⁠, we can obtain the average, over the entire cloud (our abstract object), gravitational energy per unit mass of the fluid element: 〈G〉 ∼ 〈ρ〉_c, hence 〈G〉 = 3G₀. For p = 2, the accretion energy per unit mass remains constant from one scale to the other: |$\langle v_{\rm a}^2\rangle \propto \varrho (\ell)^{-2} \ell ^{-4}\propto \ell ^{2p-4}\propto \ell ^0$|⁠. Then its average value is 〈A〉 = A₀. The turbulent therm only matters if β = 0, hence the turbulent energy does not scale, that is why 〈T〉 = T₀.

To sum up, in the physical regime far from the core (i.e. when its presence can be neglected) there exists a single solution for the density profile ϱ(ℓ) = ℓ⁻². The latter corresponds to three possible configurations for the energy balance as listed below.

• Balance between the gravitational and thermal energy (Case 1):

  $$\begin{eqnarray} 4 - \left\langle G \right\rangle \approx 0~~, \end{eqnarray}$$

  (25)

  which can be realized for trans- or subsonic flows, if the turbulent velocity scales with β > 0.
• Balance between the gravitational, turbulent, and thermal energies (Case 2):

  $$\begin{eqnarray} \left\langle T \right\rangle + 4 - \left\langle G \right\rangle \approx 0~~, \end{eqnarray}$$

  (26)

  which can be realized for trans- or subsonic flows, if the turbulent velocity does not scale (β = 0).

  The energy of the fluid element E₀ remains undetermined both in Cases 1 and 2.

• Balance between the gravitational and turbulent energy (Case 3):

  $$\begin{eqnarray} \left\langle T \right\rangle - \left\langle G \right\rangle \approx 0~~, \end{eqnarray}$$

  (27)

  which can be realized only for supersonic flows, if the turbulent velocity is independent of spatial scale (β = 0).

  The energy of the fluid element in Case 3 can be either positive (E₀ = 〈A〉 > 0; implying an unstable energy state), if the accretion flow is also supersonic, or approximately zero (implying a marginally bound energy state), if the accretion flow is trans- or subsonic.

3.4.2 The regime near to the core

A solution in this regime was obtained in Paper II only for supersonic flows – it is a free-fall case characterized by density profile ϱ(ℓ) = ℓ^−3/2 and an energy balance A₀ − G₁ ≈ 0. It was shown that the energy of the fluid element is zero and thus the state of the gas can be characterized as marginally bound.

This solution cannot be reproduced from virial analysis within our model. One cannot neglect all terms containing ℓ₀ in equations (13) and (20) to obtain their simplified forms (14) and (21), accordingly. These approximations are not valid in the whole range of scales near the core: ℓ ≳ ℓ₀. In regard to the upper limit of this range, one can speculatively look for a transitional regime described by the solution p = 3/2 according to the last row in Table 2. However, the assumption |$\langle \ddot{I} \rangle \simeq 0$| breaks in the regime near to the core and hence the equation (22) is not anymore valid.

4 DISCUSSION

4.1 Basis of the model and main results

In a series of three papers (Papers I and II and this work), we aim at modelling of the ρ-PDF of a self-gravitating isothermal turbulent fluid. This corresponds to late stages of the evolution of MCs when the dense part of the PDF transforms to PLT(s) (Kritsuk et al. 2011; Girichidis et al. 2014; Jaupart & Chabrier 2020). Our approach is statistical and analytical. We introduce an ensemble of MCs as an useful approach to describe a set of MCs with the same ρ-PDF, size, boundary densities, and temperature (Donkov et al. 2017). Furthermore, we introduce the concept of average representative object. This is the member of the ensemble with the simplest geometry and physics: spherically symmetric, radially isotropic, and in dynamical equilibrium (steady state) in regard to the motions of fluid elements and gas molecules. An equation for the density profile is derived from the equations of the medium written for the average representative object. Approximate solutions of this equation yield the density profile in two regimes (Papers I and II): far from the cloud core (ϱ(ℓ) = ℓ⁻²) and near to the cloud core (ϱ(ℓ) = ℓ^−3/2). In these two regimes, we presumed supersonic turbulent and accretion flows. In this study, we applied the EVT to derive the density profile from an extended model that allows for super-, trans-, and subsonic flows. A solution turns out to be possible only in the regime far from the core, where the profile ϱ(ℓ) = ℓ⁻² is reproduced by use of this independent approach.

The energy balance corresponding to this regime can be realized in three configurations (Cases 1, 2, and 3). Cases 1 and 2 are possible only for trans- or subsonic flows. In the former, gravity is balanced by the thermal energy, since the turbulent velocity depends on spacial scale (β > 0) and, hence, is subdominant in the balance equation (22). Thus one obtains: 〈G〉 ≈ 4. In Case 2, turbulence is important because it does not scale. One obtains the equation: 〈T〉 + 4 − 〈G〉 ≈ 0, i.e. gravity is balanced by the turbulent and thermal energies. Unfortunately, both in Cases 1 and 2, the energy of the fluid element E₀ cannot be determined since the energy balance obtained in Paper II applies only for supersonic flows.

Case 3 is possible for supersonic turbulent flows. The energy balance equation obtained in this work reads: 〈T〉 − 〈G〉 ≈ 0, i.e. gravity is balanced by turbulence. Here, one can use the energy balance equation from Paper II to compute the energy of the fluid element. There are two possibilities. If the accretion flow is also supersonic, then E₀ ≈ 〈A〉 > 0, i.e. the energy state is unstable. For trans- or subsonic accretion flows, E₀ ≈ 0 which corresponds to a marginally bound energy state.

In the regime near to the core, the second derivative of the moment of inertia in the EVT cannot be neglected. The only safe conclusion here is that the state of the gas is marginally bound and the core vicinity is not virialized.

An interesting question arises here. Can we extend the scope of validity of the solutions obtained in Papers I and II also to trans- and subsonic turbulent and accretion flows?

As demonstrated in Paper II, section 3.2, if the fluid element is near to the core, the leading term in equation (14) there must be of order ℓ⁻¹ and hence the thermal term 2(− p)ln (ℓ) is not important. Then, only accretion can balance the gravity of the core. Therefore, the solution p = −3/2 in the regime near to the core is valid not only for supersonic, but also for trans- or subsonic flows.

The solution in the regime far from the core is of more interest in the context of this work. Let us consider equation (11) in Paper II, section 3.1. The thermal term is again 2(− p)ln (ℓ) and is neglected due to the supersonic coefficients A₀ and T₀. However, in case of trans- or subsonic flows, A₀ and T₀ are of order unity or less, and then the term containing ln (ℓ) can be neglected only if the fluid element is very close to the outer cloud boundary, i.e. ℓ = 1 − ϵ, where ϵ ≪ 1. This is a very narrow range indeed, in contrast to the expectation that a ‘far from the core’ regime should span at least one order of magnitude or more. A possible solution of this problem can be to renormalize the model (and its equations), shifting the cloud boundary at l_c ≡ (1 − ϵ)l_c and then repeating all considerations. The removed shell with thickness ϵ is added accordingly to the matter located outside the cloud, which is assumed to obey radial symmetry (see Section 2). Thus this operation does not change anything in the model; it can be repeated n ≫ 1 times until the outer cloud boundary reaches values like (1 − nϵ) ∼ ℓ₀, i.e. until the fluid element in consideration approaches the regime ‘near to the core’ where the gravity of the core becomes important. Applying the above procedure, one can obtain the solution for the density profile ϱ(ℓ) = ℓ⁻² for trans- or subsonic flows and the energy balance equations (12) and (13) in Paper II will be also valid.

By the same reasoning, one can calculate the energy of the fluid element E₀ for trans- and subsonic flows in the context of this work. If the turbulent velocity scales (β > 0) and thus turbulence is not important, then E₀ ≈ 〈A〉 − 4, because from Paper II one has 〈A〉 − 〈G〉 ≈ E₀. If β = 0 and the turbulence matters, then again E₀ ≈ 〈A〉 − 4, because from Paper II one has 〈A〉 + 〈T〉 − 〈G〉 ≈ E₀, accordingly. In both cases, the dimensionless accretion energy 〈A〉 ≲ 1 and one is to expect that E₀ is negative or at least approximately zero. Therefore, the outer shells of the cloud should be stable or marginally bound.

4.2 The model framework in the context of previous studies

The presented results are linked to those from several previous studies. The case of collapsing homogeneous gas ball without accretion was originally investigated by Larson (1969) and Penston (1969), with self-gravity and isothermal gas pressure being the main acting forces. The numerical solutions of the equations of the medium yield a density profile with a slope of p = 2 in the outer shells. Shu (1977) and Hunter (1977) also studied the problem analytically and obtained two density profiles: p = 2 for the outer shells (but in static equilibrium: pressure provides support against gravity) and p = 3/2 for the free-falling inner shells near to the singularity. The latter solution describes the so-called inside-out collapse. Later, Whitworth & Summers (1985) have shown that the solutions of Larson–Penston and of Shu–Hunter are just the opposing ends of a two-parameter family of solutions: for intrinsically unstable and initially quasi-hydrostatic clouds, respectively. This general picture was reproduced by the numerical study of Ogino, Tomisaka & Nakamura (1999) who constrained its applicability to an epoch of accretion rate which is constant in time. (For stability analysis of the solutions, see also Hanawa & Nakayama 1997.)

Using numerical simulations, Naranjo-Romero, Vázquez-Semadeni & Loughnane (2015) explore a collapsing core embedded in a larger medium (called ‘cloud’) and accreting material from it. They also obtain a density profile p = 2 in the outer shells of the core as the latter collapse with the cloud. Li (2018) derives a density profile with slope p = 2 taking into account the interaction of gravity, accretion, and turbulence in the course of the collapse. He asserts that the isothermal state is not a necessary condition for such profile – the slope is universal in case of a scale-free gravitational collapse. Recently, Jaupart & Chabrier (2020) develop the approach of Pan et al. (2018, 2019) and derive an analytical theory describing the PDF of density fluctuations in supersonic turbulent star-forming clouds in the presence of gravity. The ρ-PDF is considered as a dynamical – rather than stationary – entity. The authors demonstrate the appearance of two PLTs in the ρ-PDF, with two characteristic exponents, corresponding to different evolutionary stages in regard to the balance between turbulence and gravity. It is argued that a second PLT, with a characteristic slope of −1.5 (p = 2), is a signature of regions in free-fall collapse.

A number of observations also confirm our result for the regime far from the core. The radial density profile in star cluster-forming molecular clumps is found to be very close to p = 2 (Mueller et al. 2002; Evans 2003; Wyrowski et al. 2012, 2016; Palau et al. 2014; Csengeri et al. 2017; Zhang & Li 2017). Similar results, with some larger spread of p, have been obtained for the envelopes of prestellar/protostellar cores whereas the central region (of small relative volume) exhibits a flat profile (e.g. Chandler & Richer 2000; Shirley et al. 2000; Motte & André 2001; Smith et al. 2011), like in our model. This is generally consistent with the shape of stable or unstable Bonnor–Ebert spheres.

4.3 Plausibility of the results

At least two issues are to be discussed in regard to the plausibility of the obtained solutions.

First, Cases 2 and 3 for the equations (26) and (27) of energy balance of the fluid element require lack of scaling (β = 0) of the turbulent velocity fluctuations. An essential feature of both subsonic and supersonic turbulence is the positive power-law scaling relation with the size of considered volume. The lack of such scaling is to be interpreted rather like an additional sound speed but due to a macro-motion of fluid elements while thermal motion occurs at micro-scales. But although β = 0 leads to a caveat, several observational works testify for such phenomenon. Goodman et al. (1998) found lack of scaling of the non-thermal velocity dispersion in the so-called coherent dense cores (extracted by isodensity contours of the NH₃ emission): low-mass (1–10 M_⊙) star-forming clumps with typical sizes of ∼0.1 pc. Similarly, Caselli et al. (2002) obtained a constant non-thermal velocity dispersion in the interiors of dense low-mass cores extracted from N₂H⁺ maps. Pineda et al. (2010) studied a sample of dense NH₃ cores, in whose inner parts the non-thermal velocity dispersion does not depend on scale. Gibson et al. (2009) observed dense high-mass CS cores wherein no velocity scaling is present. Analysing a sample of dense high-mass (∼10²–10³ M_⊙) clumps delineated through HCN and CS isodensity contours, with sizes of several tenths of pc, Wu et al. (2010) found a very weak scaling of velocity. By use of a different clump-finding technique (gaussclumps), Veltchev et al. (2018) extracted high-mass (∼10–10³ M_⊙) CO condensations of low density and sizes between 0.1 and 1 pc in the star-forming region Rosette, which exhibit no (β = 0) or a very weak scaling (β = 0.2; see their fig. 7).

A second issue is the condition for trans- or subsonic flows in Cases 1 and 2 (equations for energy balance 25 and 26, respectively). The Larson–Penston solution for isothermal self-gravitating sphere (Larson 1969; Penston 1969) predicts that in case of established density profile ϱ = ℓ⁻² at later stages of the collapse the velocity (accretion) field tends to a supersonic value (∼3.3c_s). Guerrero-Gamboa & Vázquez-Semadeni (2020) find that the turbulent velocity appears to approach a constant fraction of the infall speed and thus, if the infall speed is sufficiently supersonic, so should be the turbulent velocity. Then the found solutions in Cases 1 and 2 do not matter. Indeed, the dense high-mass cores observed by Gibson et al. (2009) and Wu et al. (2010) exhibit supersonic non-thermal velocity dispersions. However, dense ‘coherent cores’ (Goodman et al. 1998; Caselli et al. 2002; Pineda et al. 2010) and the high-mass cores of low density extracted by Veltchev et al. (2018) are characterized by trans- or subsonic non-thermal velocity dispersions.

To sum up, the lack of scaling of turbulent velocity in Cases 2 and 3 and the condition for trans- or subsonic flows in Cases 1 and 2 may cause concerns how plausible are the obtained results. Nevertheless, there are several observational studies indicating that the solutions from our model are valid for several different types of cores. Case 2 may apply to dense ‘coherent cores’ (Goodman et al. 1998; Caselli et al. 2002; Pineda et al. 2010) and to massive CO clumps associated with the low-density envelopes of star-forming cores (Veltchev et al. 2018). On the other hand, Case 3 is similar to the dense high-mass objects studied by Gibson et al. (2009) and Wu et al. (2010).

It should be acknowledged, however, that two of those data sets (those from Gibson et al. 2009; Wu et al. 2010), along with several others, were discussed by Ballesteros-Paredes et al. (2011) as possible evidence that the line widths are actually reflecting infall motions, rather than turbulence, so that the velocity dispersion obeys a scaling: |$\sigma _{\rm v} \sim \sqrt{\Sigma \ell }$|⁠, where Σ is the column density, rather than a simple turbulent scaling of the form: σ_v ∼ ℓ^β (see figs 1 and 2 of Ballesteros-Paredes et al. 2011).

4.4 Caveats

We recognize that our model is idealized and emphasize this mostly in view of four possible caveats.

First, in regard to the cloud geometry, spherical symmetry is not the exact shape of clumps and cores in MCs. However, a 3D object that is characterized by a single spatial size can be reconstructed from the volume-weighted ρ-PDF, if and only if, the latter is a monotonic function. In such case, a given scale as defined in equation (1) can be linked to the density and thus the necessary one to one relation between density and radial position is generated.

The choice of a sphere is appropriate since spherical geometry is a first-order approximation for a system without magnetic field and rotation. Moreover, the pressure forces tend to isotropize the system and its lowest-energy state corresponds to spherical shape. Also, a small deviation from the spherical symmetry does not cause significant changes in the results from the model. The considered average representative object is averaged in time and thus turbulence, which is homogeneous and isotropic by assumption, is not inconsistent with the spherical symmetry. We remind the reader that our aim is not to reproduce the complete dynamics and morphology of the fluid in an individual cloud, but rather to assess the general characteristics of a whole MC ensemble.

A second caveat – directly linked to the assumption of spherical symmetry – might be the neglection of rotation and magnetic fields in our model. In the first case, angular momentum conservation should be taken into account and the state of minimum energy is a rotationally supported disc and not a sphere. However, the typical sizes of discs around protostars (tens to a few hundreds of au) are much smaller than the spatial scales at which our model operates. Such discs form inside the cloud core, while outside it the velocity component (and the corresponding energy term) due to radial inflow dominates by far. Inclusion of magnetic fields in the model would also break the spherical symmetry and change the energy balance equation. Such construction would be justified in regard to an earlier epoch of cloud evolution characterized by possible magnetic support against collapse (Shu, Adams & Lizano 1987) while the pronounced PLT of the ρ-PDF in our study implies advanced stages of self-gravitating clouds. To sum up, although we recognize the importance of rotation and magnetic fields for a realistic description of cloud dynamics, their neglection cannot affect significantly the density profile at the considered spatial scales and evolutionary time of the cloud.

The third possible caveat concerns the assumption of steady state. Indeed, a steady state or an approach to it has been suggested by various authors (e.g. Murray & Chang 2015; Murray et al. 2017; Guerrero-Gamboa & Vázquez-Semadeni 2020; Gomez et al. 2021). We note also that a quasi-steady state is not inconsistent with the collapse, as long as there is an external mass reservoir that supplies the accretion (implicitly assumed in our model). On the other hand, many authors consider MCs as dynamically evolving objects that are never in a steady state (see Vázquez-Semadeni et al. 2019, and the references there in). We assert that our assumption is justified as long as the final stages of MC evolution are considered, when the high-density part of the PDF consists of (a) well developed PLT(s). Numerical simulations show that the structures corresponding to this/these PLT(s) are in steady state (Kritsuk et al. 2011; Girichidis et al. 2014; Khullar et al. 2021).

The fourth caveat arises from the requirement for lack of turbulent-velocity scaling (β = 0) in Cases 2 and 3. It seems strange in view of the typical picture in (subsonic and supersonic) turbulent media. Nevertheless, as discussed in the Section 4.3, there exist several observational examples for such phenomenon which suggest that it is still not completely understood. Goodman et al. (1998) and Pineda et al. (2010) found a jump of the non-thermal velocity dispersion at the boundary of their ‘coherent cores’: it is significantly larger outside them and scales positively with size. Goodman et al. (1998) explain this behaviour with ‘the marked decline in the coupling of the magnetic field to gas motions due to a decreased ion/neutral ratio in dense, high filling factor gas’ and hence the appearance of a dissipation threshold of MHD turbulence at the cores’ edges. Pineda et al. (2010) discuss the absence of velocity scaling and attribute it to possible changes in density, pressure and/or magnetic field at the cores’ boundaries.

Last, we recognize that the presented model (Papers I and II, and this work) is applicable only to the protostellar stage of MC evolution. The obtained solutions for the density (∼ℓ^−3/2) and velocity (∼ℓ^−1/2) profiles in the regime near to the core are monotonically increasing functions. This behaviour of the flow characterizes the protostellar cores while in the prestellar case the density is constant while the velocity decreases, in the region near to the centre, where the singularity (protostar) does not form yet (Whitworth & Summers 1985). It should be noted though, that the approach used to derive the solution in the regime far from the core can be applied to the outer regions of prestellar cores. In Papers I and II, the conditions to obtain this solution were: (i) a steady state, and (ii) a negligible effect of the gravity of the central region on a moving fluid element in the outer regions. In the framework of this paper, one should add a requirement for negligible internal energy of the central regions and for negligible change of the moment of inertia (see Sections 3.2, 3.3.1, and 3.3.3). All these conditions are satisfied in the case of prestellar core.

4.5 On the correspondence between the profile ϱ(ℓ) = ℓ⁻² and the assumption of steady state

Finally, we note an important feature of the density profile obtained for the regime far from the core. Considering a spherically symmetric cloud that accretes matter radially from its surroundings, of a density profile of type ϱ(ℓ) = ℓ⁻² is equivalent to a steady state flow through the cloud shells. In Papers I and II and in this work, we derived a profile ϱ(ℓ) = ℓ⁻² in the regime far from the core if the steady state flow is assumed.

Let us demonstrate that the opposite is true as well: an accreting cloud characterized by a density profile ϱ = ℓ⁻² and approximate spherical symmetry should be in a steady state. Let v_nt be the non-thermal velocity dispersion in the cloud dominated by the infall motions. Then, as shown by Heyer et al. (2009) and Ballesteros-Paredes et al. (2011), |$v_{\rm nt}^2\propto G \Sigma R$|⁠, where Σ is the surface mass density and R is the radius of the considered shell (i.e. R corresponds to ≡ ℓ). Combining this relation with the one for the density profile and taking into account that Σ ∝ ϱℓ ∝ ℓ⁻¹, one obtains v_nt = const(ℓ). On the other hand, the continuity equation in spherical coordinates reads

$$\begin{eqnarray} \frac{\partial \varrho }{\partial t} + \vec{\triangledown } \cdot (\varrho \vec{v}) =\frac{\partial \varrho }{\partial t} + \frac{1}{\ell ^2} \frac{\partial }{\partial \ell } \left(\ell ^2 \varrho v_{\rm nt} \right) = 0~, \end{eqnarray}$$

where we assume for the averaged velocity: |$\vec{v} \approx (v_{\rm nt},0,0)$|⁠. Using the considerations above, it is obvious that ℓ²ϱv_nt = const(ℓ). Hence ∂ϱ/∂t = 0, i.e. the density field does not depend explicitly on the time. Thus the cloud is in approximate steady state.

The equivalence between the assumption of steady state and the density profile ϱ = ℓ⁻² has been also shown by Li (2018). Recently Gomez et al. (2021) demonstrated as well that such slope is an attractor – other slopes approach this value in the course of the cloud evolution. This is also consistent with the obtained possibility in the case of trans- or subsonic flow for a negative total energy of the fluid element in the regime far from the core in our model which implements stability of the outer shells.

5 CONCLUSIONS

In this study, we extend the model of ensembles of self-gravitating turbulent clouds developed in Papers I and II, applying the EVT and broadening the scope of consideration to include super-, trans-, and subsonic flows. Our analysis shows that the solution ϱ(ℓ) = ℓ⁻² is confirmed in the regime far from the core. Moreover, the obtained energy balance equations yield the value of the energy of the fluid element which remains indeterminate in our original approach, by use of the differential equations of the medium (Papers I and II). This energy can be negative in the cases of trans- or subsonic flows, i.e. the energy balance in the solution far from the core shows that the outer cloud shells can be stable. Thus this solution for the density profile might correspond to a stable dynamical state and must be observable. In the regime near to the core, one cannot derive an equation for the density profile, this is because the second derivative of the moment of inertia of the gas could not be neglected. In this regime, we arrive only at the qualitative conclusion that the gas near to the core is not virialized. Also, the energy of the fluid element in this regime is close to zero which suggests that the gas state can be characterized as marginally bound. In this way, the model is further developed by use of the EVT as an independent approach.

ACKNOWLEDGEMENTS

The authors are thankful to the anonymous referee whose remarks and suggestions helped us to substantially improve this paper. SD and TV acknowledge support by the DFG (Deutsche Forschungsgemeinschaft) under grant KL 1358/20-3. TV thanks for the additional funding from the Ministry of Education and Science of the Republic of Bulgaria, National RI Roadmap Project DO1-176/29.07.2022. RSK acknowledges support from the Deutsche Forschungsgemeinschaft (DFG) via the collaborative research centre (SFB 881, Project-ID 138713538) ‘The Milky Way System’ (sub-projects A1, B1, B2, and B8) and from the Heidelberg cluster of excellence (EXC 2181 - project-ID 390900948) ‘STRUCTURES: A unifying approach to emergent phenomena in the physical world, mathematics, and complex data’. He also thanks for funding from the European Research Council in the ERC synergy grant ‘ECOGAL - Understanding our Galactic ecosystem: From the disc of the Milky Way to the formation sites of stars and planets’ (project-ID 855130).

DATA AVAILABILITY

No new data were generated or analysed in support of this research.

Footnotes

1

In Papers I and II, we considered only the case of supersonic turbulence and accretion, but here the scope is extended, in view of the obtained possibilities for the energy balance which correspond to the solution for the density profile.

2

Strictly speaking, this relation makes sense only if the density PDF is monotonic – which is the case considered in this work.

3

After the averaging, the turbulent velocity does not contribute to the motion of the fluid element. However, its contribution to the kinetic energy term 〈u²〉 is non-vanishing due to the scalar nature of the latter. The velocity is first squared, then averaged.

4

All terms are averaged in respect to the micro-states determined by the turbulence at a given scale ℓ. The averaged terms are put in 〈...〉.

5

We remind the reader that the core is homogeneous with density |$\rho _0\propto \ell _0^{-p}$| (Section 2). Then, on the general assumption of polytropic equation of state, the thermal energy of the gas is |$\sim \rho _0^\Gamma l_0^3\propto \ell _0^{3-p\Gamma }$|⁠, with Γ ∼ 1. In addition, if turbulence is the sole contributor to the kinetic energy the latter is |$\sim \rho _0 l_0^{3+2\gamma }\propto \ell _0^{3+2\gamma -p}$|⁠, where γ ≤ 1 is the turbulent velocity scaling exponent in the core.