Density profile of a self-gravitating polytropic turbulent fluid in a rotating disk near to the cloud core

We obtain two equations (following from two different approaches) for the density profile in a self-gravitating polytropic cylindrically symmetric and rotating turbulent gas disk. The adopted physical picture is appropriate to describe the conditions near to the cloud core where the equation of state of the gas changes from isothermal (in the outer cloud layers) to one of"hard polytrope", and the symmetry changes from spherical to cylindrical. On the assumption of steady state, as the accreting matter passes through all spatial scales, we show that the total energy per unit mass is an invariant with respect to the fluid flow. The obtained equation describes the balance of the kinetic, thermal and gravitational energy of a fluid element. We also introduce a method for approximating density profile solutions (in a power-law form), leading to the emergence of three different regimes. We apply, as well, dynamical analysis of the motion of a fluid element. Only one of the regimes is in accordance with the two approaches (energy and force balance). It corresponds to a density profile of a slope -2, polytropic exponent 3/2, and sub-Keplerian rotation of the disk, when the gravity is balanced by the thermal pressure. It also matches with some observations and numerical works and, in particular, leads to a second power-law tail (of a slope approx. -1) of the density distribution function in dense, self-gravitating cloud regions.


INTRODUCTION
Star formation in galaxies goes in molecular clouds (MCs): very cold (T ∼ 10 − 30 K) objects of irregular shape consisting of molecular gas well mixed with small amounts of dust (see Ballesteros-Paredes et al. 2020, for a review).Their internal structure is fractal within a large range of spatial scales 0.001 pc L 100 pc (Elmegreen 1997;Hennebelle & Falgarone 2012), with mean number densities varying from about 10 2 cm −3 at L ∼100 pc up to > 10 5 cm −3 in the so-called prestellar cores (L 0.1 pc).MC physics is complex due to the interplay of gravity, supersonic turbulence and magnetic fields; at late stages of cloud's evolution, accretion from the surrounding medium and feedback from newborn stars and supernovae plays also an important role (Mac Low & Klessen 2004;McKee & Ostriker 2007;Klessen & Glover 2016).Gas thermodynamics in MC is generally characterized by an isothermal equation of state (EOS) Pgas ∝ ρ Γ where ρ is the density and Γ ≈ 1.However, in denser substructures like protostellar cores the effective EOS is rather one of a 'hard polytrope' with Γ > 1 (Federrath & Banerjee 2015;Kritsuk, Norman & Wagner 2011).The observable features at such small scales may also bear imprints of strong magnetic fields and/or disk rotation.
An increasingly popular approach to study MC physics is the analysis of the probability distribution function (PDF) of density ρ-PDF (which could be obtained from numerical simulations) and of column density N -PDF (derived from observational data and simulations).A correspondence between these statistical characteristics can be obtained assuming that the cloud is roughly spherically symmetric, with a radial density profile ρ(l) ∝ l −p where l is a given radius.The shape of the ρ-PDF/N -PDF reflects the general structure and the evolutionary stage of the cloud.Isothermal supersonic turbulence usually dominates MC physics at scales above 1 pc -in that case, a nearly lognormal1 ρ-PDF (hereafter, simply PDF) is predicted from theoretical considerations (Vazquez-Semadeni 1994;Passot & Vazquez-Semadeni 1998) which is also confirmed from numerous simulations (e.g., Klessen 2000;Li, Klessen & MacLow 2003;Kritsuk et al. 2007;Federrath et al. 2010).
Most authors attribute the emergence of a PLT to the increasing role of self-gravity in the energy balance at the corresponding small spatial scales (Klessen 2000;Dib & Burkert 2005;Slyz et al. 2005;Vazquez-Semadeni et al. 2008;Kritsuk, Norman & Wagner 2011;Collins et al. 2012).Examples of the PLT evolution from numerical simulations at scales of giant MCs and of collapsing clumps are discussed in Veltchev et al. (2019).In general, the PLT slope q becomes shallower in evolving self-gravitating clouds due to formation of denser substructures (clumps, cores), with an upper limit q ≃ −1.5.The latter can be explained within various models: collapse of the so called singular isothermal spheres (Penston 1969a,b;Larson 1969;Shu 1977;Hunter 1977;Whitworth & Summers 1985), pressure-less gravitational collapse in strongly self-gravitating systems (Girichidis et al. 2014), scale-free gravitational collapse (Li 2018), collapse and dynamics of isolated gravo-turbulent cloud (Jaupart & Chabrier 2020) or dynamical equilibrium between gravity and accretion (Donkov &Stefanov 2018, 2019, hereafter, Paper I andPaper II).Some numerical (Kritsuk, Norman & Wagner 2011) and observational studies (Schneider et al. 2015c(Schneider et al. , 2022) ) indicate existence of a second PLT at the high-density end of the PDF associated with scales at which local collapses take place (Fig. 1, bottom).This PLT is usually shallower 2 , with slope q ′ ≃ −1 (corresponding to a density-profile exponent p = 3, in the case of spherical symmetry) and hints at a slower gas accretion on the small-scale substructures in the cloud.This phenomenon still lacks a clear physical explanation; possible reasons might be rotation of contracting core(s), strong magnetic fields or change in the thermodynamic EOS (Kritsuk, Norman & Wagner 2011;Schneider et al. 2015c;Donkov et al. 2021).Donkov et al. (2021) proposed a model in which a flatter second PLT is associated with the densest cloud parts (e.g., protostellar cores).The authors obtain an equation for the conservation of energy of a fluid element per unit mass.It is the analogue of Bernoulli equation in this context.In this study spherical symmetry (without rotation), steady state of both macroscopic and microscopic motions, and a polytropic EOS with Γ > 1, are assumed.The leading order solutions correspond to four cases which have different density profiles, polytropic exponents, and energy balance equations for the  fluid element.The solution with p = 3, (q = −1) and Γ = 4/3 can provide an explanation of a flatter second PLT.
In this paper we continue the efforts to build an hydrodynamical model, which explain the second PL-tail in the vicinity of emerging protostar/protocluster.We consider a rotating disk (instead of static sphere) around central core (protostar), accreting material from its surrounding medium, and the fluid in the disk is in steady state in regard to macroscopic and microscopic motions, obeying also a polytropic EOS with exponent Γ > 1.The analogue of Bernoulli equation for the current model is derived from Euler equation with the application of the assumptions mentioned above.It describes, as in the previous study, the conservation of energy of a fluid element per unit mass.We find an approximate solution keeping only the leading order terms.It describes three regimes characterized by different constraints for the density profiles, the polytropic exponents, the angular velocity scaling exponents, and energy balance equations for the fluid element.Performing, in addition, dynamical analysis through the force balance equation, we conclude that only one regime is in agreement with two approaches (energy and force balance equations).Namely, the second regime, where the gas pressure force balances the gravity and the rotation of the disk is sub-Keplerian.In particular it gives, as a possible so-lution, density profile of a slope p = 2, which in the regarded cylindrical symmetry of the disk raises a PL-tail of a slope q = −1, and polytropic exponent Γ = 3/2.We consider this as a possible explanation of the second PLT.
The paper has the following structure.The model of the gas medium in the vicinity of the core of the cloud is presented in Section 2. The equation for the density profile is derived in Section 3. In more details: Section 3.1 contains comments on the equations of the medium; the derivation of the general form of the equation of energy conservation for a fluid element (per unit mass) is done in Section 3.2; each term of the latter equation is given in explicit form in Section 3.3; the equation for the density profile near to the cloud core is finally obtained in 3.4.The solution is presented in Section 4, as the energy balance equation is analysed and its solutions in three different physical regimes are obtained in Section 4.1, while the dynamical analysis is presented in Section 4.2.A discussion on the obtained solutions and on the model as a whole is given in Section 5. Our conclusions are presented in Section 6.

SETUP OF THE MODEL
We consider a system of molecular hydrogen and dust, which consists of a rotating spherical gaseous core and an accretion gaseous disk rotating around the core.Feedback from the new-born stars in the core are neglected.We consider the system at the earliest period of its life.The outer size of the disk we label with lc and assume that lc ≫ h and also lc ≫ l0, where h is the thickness of disk and l0 is the inner size of disk, which is of the order of core's size.The core is assumed to be homogeneous, for simplicity.We, also, presume the disk is rotating differentially, and the gas in it is turbulent, self-gravitating, and accreting onto the core.We presume that there is a weak friction between the neighbouring slices (rings) of the disk, which causes the redistribution of angular momentum of the fluid elements while they are falling onto the core.Although this friction plays a role for the fluid dynamics in the disk we assume that in terms of the energy balance of each fluid element its contribution is negligible.So one can expect that the rotational kinetic energy is much larger than the kinetic energy due to the radial accretion.This is also confirmed from the results (see Section 4) and so the presented here model is self-consistent.We also suggest there exists a radial density profile ρ(l) for the disk, where lc l l0 is an arbitrary radius (or scale).
We normalize the gas density and scale to ρc (the density at disk's edge) and lc, respectively.From now on we will work with dimensionless density ̺ ≡ ρ/ρc and radius ℓ ≡ l/lc, where it is more appropriate.Accordingly the dimensionless density profile in the disk is ̺(ℓ).
We assume the fluid in the disk to be in a steady state regarding to both: its macro-motion (the motion of fluid elements), and its micro-motion (the thermal motion of gaseous molecules), and so the temperature is locally constant in time.An arbitrary fluid element, residing at scale (radius) ℓ, takes part in three independent motions -accreting radially, rotating around the central core, and moving turbulently.So its velocity field, respectively, reads: u = ua + urot + ut, where ua is the accretion (or radial in-fall) velocity, urot is the velocity of rotation, and ut is the turbulent velocity, accord-ingly.The turbulence, we presume, to be locally homogeneous and isotropic for scales lc l l0.Combining the latter presumptions with the steady state assumption, one can conclude that averaging the model (and hence the velocity field) in regard to time can be replaced by ensemble averaging in regard to an ensemble of copies of the system raised by possible different turbulent velocities of fluid elements at a fixed moment.So, if one averages turbulent velocity of a given fluid element, one obtains: ut = 0 (due to the local isotropy), and hence the averaged total velocity field reads: u = ua + urot = ua + urot.We note also that averaged scalar products between every two of the three velocity fields are zero: hence the three averaged velocity fields are orthogonal to each other.From now on we will use, where it is appropriate, dimensionless velocity fields substituting u for v ≡ u/cs, where cs is the isothermal sound speed explained below.
The entire system (core plus disk) is a central part of a larger cloud which material obeys radial symmetry and has its own density profile (probably different from the profile in disk).This cloud consists of molecular hydrogen, and is isothermal, self-gravitating and turbulent.It accretes matter from its surroundings and this material moves through the cloud shells to reach the central part where our system is located.We presume that the cloud shells are in the steady state, as well, regarding accretion, turbulence, and thermodynamics.These shells cause a constant gravitational potential Φ ext in the volume where our system resides, and therefore it does not affect the motion of fluid elements in the disk (for further details see Paper II).
Regarding thermodynamics of the gas in the disk we assume a polytropic equation of state (EOS): (1) It reflects the typical assumption for the medium in the densest part of the cloud, near to the core (Federrath & Banerjee 2015).When Γ = 1 this EOS reduces to the isothermal one Pgas = c 2 s ρ, with cs being the speed of sound (Federrath & Banerjee 2015) at some fixed temperature (10K T 30K), which characterizes the outer cloud matter.From the comparison of the two equations we see that P0 = c 2 s ρc.With this choice we can connect continuously the two EOS of the molecular gas: the isothermal one (for the outer shells) and the 'hard polytropic' one (for the disk).
We neglect magnetic fields and stellar feedback (in spite of they proved to be of great importance; see for instance: Lebreuilly et al. (2021Lebreuilly et al. ( , 2023)); Hennebelle et al. (2022)), and we note that this is a considerable simplification of the reality.Also the disk is presumed to be very thin, although, from simulations (see for example Bate (2018); Lebreuilly et al.Hence we deal with a very simple hydrodynamic model of protostellar system, but we consider this is admissible in regard to our efforts to obtain the density profile ̺(ℓ) (and hence the PDF) near the core, in first approximation.

DERIVATION OF THE EQUATION FOR
DENSITY PROFILE NEAR TO THE CORE

The fluid equations
We start with the equations of the medium, as a strong base in our aim to obtain the density profile ̺(ℓ) near to the core.They read as follow: -The Euler equations where the first one is the mass conservation law, while the second one is the equation of motion for a fluid element (here Φ is the total gravitational potential).The terms of equation ( 3) will be specified in Section 3.3.The dissipative terms and the friction between neighbouring rings have been neglected.The latter was also mentioned in Section 2. The former is justified by our assumption of statistical equilibrium in regard to turbulence.
-The polytropic EOS (1) of the gas introduced in the previous Section 2.
-The equation which determines the gravitational potential (the Poisson equation), according to the given density distribution ∆Φ = 4πGρ . (4)

General form of the equation of energy conservation for a moving fluid element
The equation for ̺(ℓ) must take into account the assumptions for the symmetry and the physics of the idealized object that we describe.First we remind the reader, according to Donkov et al. (2021), the derivation of the equation in its general form (here equation 5).The derivation is independent of the symmetry of the system in regard.This repetition we considered as very useful for readability of the text.It is as follows: i) Starting with some algebra of polytropic EOS: From ii) Then we multiply the above equation by an infinitesimal displacement d r = udt directed to the vector field u: and obtain for the left hand side the following expression iii) At that point we introduce dimensionless variables v 2 ≡ u 2 /c 2 s and φ ≡ Φ/c 2 s and hence arrive at: Up to this point the considerations were rather general.From now on, we will focus on our idealized object and on the motion of the fluid elements through the disk, in particular.We use brackets to designate the quantities which belong to it and are received after ensemble averaging with respect to the chaotic turbulent motion of the fluid element in a given space volume.By their use, the equation obtained at step iv) above can be rewritten: Here we apply the assumption for steady state, from which stems that ∂ ̺ Γ−1 /∂t = 0, and ∂ φ /∂t = 0.This leads to: We assumed that turbulence is locally homogeneous and isotropic.Hence, it does not affect the ensemble averaged motion of the fluid elements.As a result, it is a superposition of radial in-fall and rotation about the centre.Turbulence, however, must be taken into account when the kinetic energy term v 2 /2 is calculated, because the velocity v is first squared and then averaged.So one has: Moreover, the system in regard is described by an averaged density profile, which is closely connected with the PDF, which, in its turn, is averaged by assumption (see Paper I).Then, ̺ Γ−1 = ̺ Γ−1 .Following the above reasoning one gets an equivalent form of equation ( 5): 3.3 Explicit form of the terms in equation ( 6) The derivation of the explicit form of the terms in equation ( 6) is done as follows.

Kinetic term
As shown in the previous Section 3.2, the kinetic energy term reads: where v 2 a is the accretion kinetic energy per unit mass, v 2 rot is the rotational kinetic energy per unit mass, and v 2 t is the turbulent kinetic energy per unit mass.
In general, turbulent velocity fluctuations in molecular clouds obey a power-law scaling relation u = u0(L/1 pc) β (Larson 1981;Padoan et al. 2006;Kritsuk et al. 2007;Federrath et al. 2010), with some normalizing factor u0 and scaling index 0 β 1.We adopt this relation and assume v 2 t depends on the scale ℓ as: where the constant T0 = (u 2 0 /c 2 s )(lc/pc) 2β is the ratio of the turbulent kinetic energy per unit mass of a fluid element to the thermal energy per unit mass, at the largest scale of the disk; and 0 β 1 is a scaling exponent of the turbulent flow.
The accretion kinetic term v 2 a in its explicit form can be got from the continuity equation, as it is shown in Paper I (see Section 3.3 there).Here we apply a similar treatment but adopting cylindrical symmetry, which is more appropriate for the disk.The main steps, in cylindrical coordinates (l, ϕ, z), are as follows.The continuity equation (Eq.2) is averaged with respect to the ensemble of the micro-states of the system, caused by the turbulence: Applying the steady state assumption one obtains ∂ρ/∂t = 0, and hence: If the motion of the fluid element is averaged the turbulent velocity vanishes (i.e: u = ua + urot = ua + urot = (ua, urot, 0)).And also taking into account that the density profile is averaged by assumption, we further obtain: where we implicitly assume that ρ, ua and urot depend only on the scale l, but not on the polar angle ϕ and on the z−coordinate.Finally, one obtains the equation: and, hence, a formula for v 2 a : Here the dimensionless coefficient A0 is obtained as the ratio of the accretion kinetic energy term at the disk edge to the thermal kinetic energy, at isothermal EOS, per unit mass.
The rotation kinetic energy v 2 rot can be written in the form: v 2 rot = Arotω(ℓ) 2 ℓ 2 , where ω(ℓ) is the dimensionless angle velocity at scale ℓ, and Arot is the ratio of the rotational kinetic energy at the disk's edge to the thermal energy (at isothermal EOS), per unit mass.Here we presume an explicit formula for ω(ℓ) in the form of power-law: ω(ℓ) = ℓ −c , where c > 0 is a positive exponent.The latter means that the closer is a ring with radius ℓ to the centre, the larger is its angular velocity ω(ℓ) (remember that ℓ < 1).Hence the rotational kinetic energy reads: Therefore, finally, the full kinetic energy term is as follows: (11)

Gravity potential term
In order to calculate the gravitational term φ in equation ( 6) we consider the disk as built of narrow rings with masses ∆M and radii l.Then an arbitrary ring causes a potential φ(l, a) at a distance a from the centre.There are two possible approximations3 : i) If l > a, then: ii) If l < a, then: In both cases we confine the considerations only to the leading order term in regard to simplicity, and due to the series converge fast.Moreover, the higher order terms will raise subdominant terms in the explicit form of equation ( 6).Hence, if the fluid element resides at scale a, the first case will help us to calculate the gravitational potential caused by the rings located in disk plane at larger radii, and the second case to account for those located below it.Making use of the above treatment we are able to write the explicit formula for the gravitational term: where the first term is caused by the rings above the fluid element, the second term is due to the rings below the fluid element, and the last term presents the gravity of the spherical core.The dimensionless coefficients G0 and G1 are defined as follows: G0 ≡ 2πGlchρc/c 2 s and it means (with a precision to a factor of order unity) the ratio of the gravitational energy of homogeneous disk (in the first order approximation) to the thermal energy, for an isothermal EOS, per unit mass; G1 ≡ GM0/lcc 2 s and this coefficient has the same physical meaning but for the central spherical homogeneous core.
To calculate the gravitational term in an explicit form one needs of formula for ̺(ℓ).In Section 2 we assume that there exists a density profile for the disk, and now we extend it to a power-law form: ̺(ℓ) = ℓ −p , where p > 0 is a positive exponent 4 .The latter raises three possible cases: 1) If p = 1, then one has: 2) If p = 2, then one has: 3) If p = 1 or 2, then one has: Although there exist three possible cases for p > 0, the most probable range for the density exponent is: 2 p 1 (see simulation by Bate (2018) and Kritsuk, Norman & Wagner (2011); and observations by Schneider et al. (2015cSchneider et al. ( , 2022))).And if we regard the scales near to the core: ℓ ℓ0, then in all the three cases for φ the leading order term(s) scales with ℓ −1 .One can see this taking into account that 0 < ℓ < 1 and the leading order term(s) must have the smallest exponent.So making use of the latter considerations, in these three possible cases we have: the first term is subdominant and the second vanishes; 2) If p = 2, then: the second term vanishes; 3) If p = 1 or 2, then: the first term is subdominant and the second is too (because the second addend in the parentheses ℓ 2−p 0 ℓ −1 ∼ ℓ 1−p < ℓ −1 for ℓ ℓ0 and 2 > p > 1 ).
To simplify further considerations we can present all the three cases with the term: −G (i) ℓ −1 , where G (1) = G1 in the first case; G (2) = G0 + G1 in the second case; and G (3) = G1 in the third case, accordingly. 4This immediately leads to an explicit formula for the accretion kinetic term: 3.4 Derivation of the equation for ̺(ℓ) At last Eq.( 6) can be written in a appropriate form, allowing us to derive an equation for the density profile: where we have substituted the total differential 'd' for the derivative 'd/dℓ', because all the terms depend only on ℓ.The expression in the brackets in Eq. ( 13) is the total energy per unit mass of a fluid element.Labelling it by E0, one receives: This is the equation for the dimensionless density profile ̺(ℓ) = ℓ −p .As one can easily see this equation depends, also, on two more unknown parameters (except for p): the hard polytropic exponent Γ > 1, and exponent c > 0, determining the power-law of differential rotation.The analysis and possible solutions of Eq. ( 14) we present in the next Section 4.
To analyse the equation ( 14) easier let us replace all the three gravitational terms with the term: −2G (i) ℓ −1 , which implies that we confine our considerations to the range: 2 p 1 for the density profile exponent, and to the scales near the core: ℓ ℓ0.The obtained equation reads: 4 STUDY OF THE DERIVED EQUATION FOR DENSITY PROFILE

Study of the equation (15)
A general approach consists in obtaining an approximate solution of the equation ( 15), where only the leading order terms does matter.Remember that term(s) with smallest exponent is/are dominant, because ℓ: 0 < ℓ0 ℓ < 1.When one rejects the subdominant term(s), the rest must consist at least of two leading terms with opposite signs, in order to balance the equation.If the dominant term is only one, or the dominant terms have one and the same sign, then the equation does not have a solution.
Let us assume that the gravitational term cannot be neglected, i.e. it is of leading order 5 .Then the dominant power of ℓ is "−1".Obviously the turbulent kinetic term T0ℓ 2β is not important for the solution.Since 2 p 1, the accretion kinetic term can not be important, also.The gravitational pull must be balanced by gas thermal pressure and/or centrifugal force.Here one can see three possible regimes, for which the equation ( 15) has an approximate solution: I) Only the centrifugal force balances the gravity.Then 2(1 − c) = −1 and −p(Γ − 1) > −1.Hence c = 3/2 and 1 < Γ < 1 + 1/p.So, if 2 p 1, then 1 < Γ < 2. The corresponding energy balance equation is: It is worth to note that the possibilities: 2(1−c) < −1 and/or −p(Γ − 1) < −1 are forbidden, since if one of them or both are realized then the gravity term does not matter, and hence the equation does not have a solution.Therefore this analysis restricts the possible values for c and Γ, as it follows: 0 < c 3/2 and 1 < Γ 2.

Dynamical analysis of the obtained results
We continue our study with an analysis of the obtained results in the previous Section 4.1.Here we change the point of view and regard the equations of motion for a fluid element in polar coordinates (equation trough the z−axis is not regarded, because through z only the turbulent motion is relevant and as it commented in Section 2 its averaged contribution is zero).We neglect the friction and dissipation forces according to the model assumptions.We account only for the gravitational and pressure-gradient forces.Then the gravity must be balanced by centrifugal and pressure forces, in radial direction.In the tangential to the circle orbit direction there are no forces.Therefore the normalized equations of motion are: and where the first equation is trough the disk radius and the second is in the tangential direction, accordingly.Make use from the scaling relation for the angular velocity ω(ℓ) = ℓ −c , one obtains that the equation ( 17) gives: (2 − c)ℓ 1−c l ≈ 0, hence l ≈ 0, and finally l ≈ 0. Using the latter we obtain from equation ( 16) the following force balance equality: One can easily see that the equation ( 18) reproduces the approximation obtained in the previous Section 4.1, but with different energy balance equations: 1) If only the centrifugal force balances the gravity, then 1−2c = −2, and hence c = 3/2; moreover −p(Γ−1)−1 > −2, and therefore 1 < Γ < 2; the corresponding energy balance is: 2) If only the gas pressure balances the gravity, then −p(Γ− 1) − 1 = −2 and 1 − 2c > −2, and hence 3/2 Γ 2 and 0 < c < 3/2; the corresponding energy balance is: 3) If both terms on the left hand side balance the gravity, then c = 3/2 and 3/2 Γ 2; and the corresponding energy balance is: The dynamical approach in this Section must be consistent with the energy equation which we solved in the previous Section 4.1, therefore the regimes I and 1, II and 2, III and 3, must be in agreement, respectively.Solving the energy balance equations in the regimes I and 1 together, we obtain that the only solution is: Arot = 0, G (i) = 0, which is trivial and not realistic.Solving simultaneously the energy balance equations in the regimes II and 2, we obtain G (i) = 1 + p together with Γ = 1 + 1/p.And, finally, in the regimes III and 3, the energy balance equations give: Arot = 0, G (i) = 1 + p, which is unrealistic like in the first regime.So one might conclude that the only reasonable approximate solution is the second regime, when the gas pressure force balances the gravity and the rotation of the disk is sub-Keplerian.
For better readability of the obtained results in this paragraph we summarized them in table 1.

Towards a fiducial model
In the presented work we suggest a model for a thin disk around protostar (the central core).Our goal is to describe the most inner part of a protostellar cloud (the larger cloud in which is embedded the regarded system, see Section 2) in a fiducial way, so that to obtain the PDF of density of this subsystem (namely the second PL-tail of the distribution) which is observed in simulations and observations.
Building the model we account for the self-gravity of the disk and core (due to the assumption for radial symmetry of the outstanding material its gravitational field does not affect the motion of fluid elements in the disk), for the thermodynamics, turbulence, and accretion.Our model is purely hydrodynamical (no magnetic fields are regarded), the disk is thin (its perpendicular size h ≪ l is much less than its radius at all regarded scales), and the feedback from new-born stars is also neglected.Due to the latter suppositions our treatment is highly idealized, comparing to the results from simulations Energy balance (from Eq. ( 15)) Energy balance (from Eq. ( 18)) and observations, where it is clear that magnetic fields essentially concern the physics and morphology of protostellar disks.Also the observed disks are rarely thin, usually they are fluffy.Finally the new-born stars in the centre produce, in addition, outflows and radiation, which make the physical picture more complicated.Furthermore we neglect the viscosity of the fluid and friction between the neighbouring rings of the disc in the considered equations.So the suggested model is very simple and reflects only a part of the real state in which resides the very inner part of a protostellar cloud.Nevertheless, given the highly difficult physics in the real protostellar disks, our approach looks justifiable as a first approximation.
Applying together the equation ( 15) for conservation of the total energy of a fluid element and equation ( 18) for the force balance, we concluded in the previous Section 4.2, that there exists only one approximate solution of these equations, which corresponds to the regime where the gravitational force is balanced by the gas pressure force and the rotation of the disk is sub-Keplerian (so the scaling index of the angular velocity is in the range: 0 < c < 3/2).Accordingly, the energy balance equation for a fluid element reads: pΓ − G (i) ≈ 0, and the polytropic exponent is given by the formula: Γ = 1 + 1/p, hence G (i) = 1 + p. Taking into account that we suppose the density profile exponent p to be in the range 1 p 2, this gives the ranges for Γ and G (i) .They are 2 Γ 3/2 and 2 G (i) 3, respectively.
We continue our analysis of the solutions comparing them with realistic power-law PDFs, obtained from numerical experiments.Our aim is to restore the density PDF near to the cloud core, where the EOS changes from isothermal to 'hard polytrope' one, and protostellar disk forms.To compare the model with observations and simulations, one has to account for the 2D symmetry of the disk.Then the slope q of the density PL-tail is related to the density exponent p, trough the formula: q = −2/p (Donkov, Veltchev & Klessen (2017), and the references therein).Then, if 1 p 2, hence −2 q −1.The slope q = −1 corresponds to density exponent p = 2, and therefore Γ = 3/2, G (i) = 3.Some observations (Schneider et al. 2015c) and simulations (Kritsuk, Norman & Wagner 2011;Veltchev et al. 2019;Marinkova et al. 2021) indicate the emergence of a second PLT of the PDF, with slope q ∼ −1.In the hydrodynamical simulation by Bate (2018), where more than 100 protostellar disks are studied, the author claim that the radial surface density profiles of isolated disks typically scaling as Σ ∼ ℓ −1 .This can be interpreted in our case as ̺ = ℓ −1 , due to the constant vertical size h of the disk.Therefore the PDF slope must be q = −2 (accordingly Γ = 2 and G (i) = 2), which is observed in some clouds by Schneider et al. (2022).
We, also, have to mention that Bate (2018), in their fig.16, present a cumulative statistics of the studied disks in regard to the radial surface density profile exponent (in our case the density exponent p), corresponding to different ranges for the radius (in regard to the mass-fraction of the disk, containing in these ranges).The main that one can conclude, from their graphics, is that the scaling exponent is not constant for all the radii in regard.Also, for all the disks p = 2 is an upper limit (and about 80 percent have for upper limit the value p = 1), and the values 0 p < 1 are presented, too.The latter means that one may expect the slopes q < −2, for the second PL-tail of the density PDF, to be observed.
For our model, the above results, would mean that the range for p can be extended to 0 < p 2. Hence in the case of the only valid approximate solution the gas polytropic exponent can take values in the interval: 3/2 Γ < ∞, and the energy coefficient accounting for the self-gravity will be in the range: 1 G (i) 3.

Caveats
We should mention several caveats regarding our model.The first one is the hypothesis of a thin disk, that helps us to simplify the equations and the calculations.Real disks (from simulations (Bate 2018;Lebreuilly et al. 2021) and observations (Andrews et al. 2018a,b;Sheehan et al. 2022a,b;Tobin et al. 2020a,b)) look more fluffy and do not have a perfect shape.These differences in shape and sizes certainly leads to different physics and more complicated equations (in real case).But our intention is to regard only the mid-plane of a fluffy disk where the more material is concentrated, and hence this thin volume dominates the physics of the whole disk, at least in the sense of energies per unit mass, what is a key point for our approach.
The second one is the neglecting of the magnetic fields, which are found to have a huge influence on the size, shape and dynamics of the fluid in disk (Lebreuilly et al. 2021;Hennebelle et al. 2022).Also, as one can see in the work by Mocz et al. (2017) the magnetic fields play significant role even at larger scales.In the latter work the authors present simulations of the collapse of prestellar cores in supersonic, turbulent, isothermal, magnetized environments, exploring the effect of the mean magnetic field strength.The density profile is ̺ ∝ ℓ −2 , which means that the scales in regard correspond to the outer shells in respect to the protostellar disk (which is not formed yet).This can be seen in their fig.5, where they present the scaling of the gas pressure, magnetic pressure, kinetic energy density, and gravitational potential energy density.Although, in their fig.4 one can see that the density PDFs are getting shallower for stronger Bfield simulations.Hence we expect that magnetic fields will have even greater impact at protostellar disk's scales.Their consideration needs, however, much more complicated physical model and hence more and more involved mathematics.So we regard this as a next step in the building of more reliable model.For now we are satisfied with a simpler hydrodynamical model.
The third one is the neglecting of a feedback from newborn stars in the core.This is related to magnetic fields, outflows and radiation that new-born stars cause.Of course, the latter factors strongly affect the shape and size of the disk, also its temperature, and the accretion rate of material (see the simulation of Lebreuilly et al. (2023)).This considerable simplification we regard as a necessary step to simplify our first approximation model.
The fourth one is the disregard of the viscosity in the disk, which is proved to explain the disk evolution at its late stages (see Lynden-Bell & Pringle (1974)), when the protostar is always burning.The latter simplification we justify with the assumption that our model regards the system at early stages after its formation.Also, if the viscosity scales according to the power-law ℓ γ , where γ 1 (see Lynden-Bell & Pringle (1974) and Hartmann et al. (1998)), then we expect it will be vanishing for ℓ ≪ 1, where our solutions are valid.
The fifth is the assumption for steady state.As it shown in Hartmann et al. (1998) the accretion rate of the material at the disk boundary is not a constant with time, but rather decline during the disk evolution.So the assumption for steady state will be broken.But if one takes into account that the changes of the large-scale velocity field (the accretion) are slow, and the establishing of local pressure equilibrium (through the sound speed), near to the core, is much faster, then comparing the typical time-scales of these processes one may conclude that the steady state hypothesis can be substantiated.We suppose that there exists a period of time when the cloud reaches a nearly steady state (Burkert 2017) and its PDF is characterized by two PLTs in the high density range, with roughly constant slopes (Kritsuk, Norman & Wagner 2011;Girichidis et al. 2014;Schneider et al. 2015c).Our model reproduces those two slopes in approximate solutions (Paper I, Paper II, Donkov et al. (2021Donkov et al. ( , 2022)), and this work).
We are aware of the considered model is highly idealized and also not applicable to many real objects.But, indeed, our goal is not to restore the complete dynamics and morphology of the fluid near to the core, but rather to assess in first approximation its general characteristics in purely hydrodynamical case.

CONCLUSIONS
In the current paper an equation for the density profile (see equation ( 15)) of a self-gravitating polytropic turbulent disk, accreting material from its vicinity and immersed in a molecular cloud is obtained.The derived profile describes the conditions of the medium near to the cloud core, where a protostar/protocluster is forming appropriately.The model is based on the framework established in Donkov & Stefanov (2018, 2019).Unlike these papers, here the EOS of the gas near to the very small and dense core has been changed from isothermal (Γ = 1) to polytropic (Γ > 1), like in Donkov et al. (2021), but compared to the latter paper the symmetry in the core vicinity is not spherical but rather cylindrical, namely rotating disk.In order to prove the invariance of the total energy per unit mass with respect to the fluid flow we assume that the accretion of matter through the boundary of the cloud and through all spatial scales down to the core is quasi-static, i.e. the cloud is in steady state.The equation that we obtain for the density profile function, or the generalized Bernoulli equation of our model, which describes the balance of the different types of energies per unit mass of a fluid element -the kinetic, the thermal, and the gravitational, is a non-linear integral one (see equations ( 13), ( 14) and ( 15)).
The method that we apply for the obtaining of an approximate solution for the density profile of the form ρ(l) ∝ l −p allows us to evaluate the power-law PDF which describes the dense cloud regions, where self-gravity cannot be neglected.The PDF slope is consistent with estimates from numerical simulations (Kritsuk, Norman & Wagner 2011;Veltchev et al. 2019;Marinkova et al. 2021) and some observations (Schneider et al. 2015c(Schneider et al. , 2022)).
Our results show three regimes characterized by different density profiles, polytropic exponents, angular velocity scaling exponents and energy balance equations for the fluid elements.Which one of them occurs depends of the balance of the centrifugal kinetic term and the thermal term, on the one side, and the gravitational term, on the other.Then applying to the model dynamical analysis (see equation ( 18)) and comparing the result from it to the regimes, obtained through the energy equation, one might conclude that the only reasonable approximate solution is given by the second regime, where the gas pressure force balances the gravity and the rotation of the disk is sub-Keplerian.Accordingly, the energy balance equation for a fluid element reads: pΓ − G (i) ≈ 0, and the polytropic exponent is given by the formula: Γ = 1 + 1/p, hence G (i) = 1 + p. Taking into account that we suppose the density profile exponent p to be in the range 1 p 2, this gives the ranges for Γ and G (i) .They are 2 Γ 3/2 and 2 G (i) 3, respectively.The above-mentioned second physical regime applies to the detection of such power-law tail in (column-)density PDFs in some observational and numerical studies of dense cloud regions (Schneider et al. 2015c(Schneider et al. , 2022;;Kritsuk, Norman & Wagner 2011;Marinkova et al. 2021).The extracted first power-law tail in PDFs at lower densities can be assessed from the consideration of the isothermal outer shells (Donkov & Stefanov 2018, 2019;Donkov et al. 2022).
At the end we would like to state that we are aware of the constraints of our results for the second power-law tail in the density PDF, which determines the structure of the inner regions of molecular clouds.In our consideration we have neglected factors such as magnetic fields, viscosity and feedback from new-born stars, which may be significant.A substantiated explanation of the second power-law tail should come from a more complex model and our work can be considered as a contribution to such efforts.

Figure 1 .
Figure 1.The figure presents examples of analytically obtained evolved density PDFs with one (top) and two (bottom) PLTs.Please, see the text for details about the slope values.The main part of the PDF is fitted with log-normal function (dotted).

Table 1 .
Summary of the results