Theory of Turbulent Equilibrium Spheres with Power-Law Linewidth-Size Relation
Sanghyuk Moon Eve C. Ostriker




0000-0002-6302-0485]Sanghyuk Moon
Department of Astrophysical Sciences, Princeton University,
  Princeton, NJ 08544, USA
0000-0002-0509-9113]Eve C. Ostriker
Department of Astrophysical Sciences, Princeton University,
  Princeton, NJ 08544, USA
  Institute for Advanced Study, 1 Einstein Drive, Princeton, NJ 08540, USA
sanghyuk.moon@princeton.edu, eco@astro.princeton.edu


    BEBonnor-Ebert
    TESturbulent equilibrium sphere
    SPSsingular polytropic sphere
    ISMinterstellar medium
    GMCgiant molecular cloud



§ ABSTRACT 
  Dense cores inherit turbulent motions from the interstellar medium in which they form.
  As a tool for comparison to both simulations and observations, it is valuable to construct theoretical core models that can relate their internal density and velocity structure while predicting their stability to gravitational collapse.
  To this end, we solve the angle-averaged equations of hydrodynamics under two assumptions: 1) the system is in a quasi-steady equilibrium; 2) the velocity field consists of radial bulk motion plus isotropic turbulence, with turbulent dispersion increasing as a power-law in the radius.
  The resulting turbulent equilibrium sphere (TES) solutions form a two-parameter family, characterized by the sonic radius r_s and the power-law index p.
  The TES is equivalent to the Bonnor-Ebert (BE) sphere when r_s→∞.
  The density profile in outer regions of the TES is slightly shallower than the BE sphere, but is steeper than the logotropic model.
  Stability analysis shows that the TESs with size exceeding a certain critical radius are unstable to radial perturbations. 
  The center-to-edge density contrast, mass, and radius of the marginally stable TES all increase with increasing average velocity dispersion.
  The FWHM of the column density profile is always smaller than the critical radius, by a larger factor at higher velocity dispersion, suggesting that observations need to probe beyond the FWHM to capture the full extent of turbulent cores.
  When applied to the highly turbulent regime typical of cluster-forming clumps, the critical mass and radius of the TES intriguingly resembles the typical mass and radius of observed star clusters.




§ INTRODUCTION


Dense prestellar cores are roughly spherical, compact (≲ 0.1 pc), centrally-concentrated objects found at the bottom of the hierarchy of ISM structure.
Due to gravity, cores are internally stratified, with 
radial profiles of density ρ characterized by a central plateau surrounded by an outer envelope approximately following ρ∝ r^-2 <cit.>.
Some cores exhibit molecular line profiles that cannot be explained by thermal motions alone, but could be interpreted as signatures of infall or turbulence <cit.>.
While the true dynamical nature of dense cores – individually and as a class – is still under debate, it is clear that a fraction of them undergo gravitational collapse, leading to star formation.


The radial column density profile of some observed cores closely matches the theoretical prediction of the isothermal equilibrium known as the BE sphere <cit.>, suggesting that those particular cores obey hydrostatic equilibrium, in which self-gravity is balanced entirely by a radial thermal pressure gradient <cit.>.
Stability analyses have shown that such equilibria are unstable unless they are truncated at a small enough radius <cit.>.

The mass and radius of the marginally stable BE sphere are called the BE mass, M_BE, and the BE radius, R_BE, respectively, given by 

    M_BE   = 4.43 M_G(ρ_c),
       = 1.18 M_G(ρ_e),
       = 1.86 M_G(ρ),


    R_BE   = 1.82 R_G(ρ_c),
       = 0.486 R_G(ρ_e),
       = 0.762 R_G(ρ),

where ρ_c, ρ_e, and ρ≡ 3M_BE/(4π R_BE^3) are the center, edge, and average density of a critical BE sphere, respectively, with ρ_c/ρ_e=14.0 and ρ̅/ρ_e=2.45.
Here, M_G and R_G are the characteristic gravitational mass and radius for an isothermal system obtained from dimensional analysis as a function of density ρ, 

    M_G    ≡c_s^3/G^3/2ρ^1/2
       = 1.27 M_⊙( T/10 K)^3/2( n_H/10^4 cm^-3)^-1/2,
      
    
        R_G    ≡c_s/G^1/2ρ^1/2
       = 0.15 pc( T/10 K)^1/2( n_H/10^4 cm^-3)^-1/2,

where c_s=(kT/μ)^1/2 is the sound speed, G is the gravitational constant, and n_H = ρ/μ_H is the hydrogen number density, with μ = 2.3m_H and μ_H =1.4m_H assuming 10% helium abundance by number.


Although the BE sphere provides essential physical insight regarding the onset of the core collapse, its applicability is limited because it neglects non-thermal pressure arising from internal turbulent motions, which are known to be present in real dense cores.
In particular, the observed turbulent velocity dispersion within dense cores is known to increase with distance from the core center approximately as a power-law <cit.>[Although some studies indicate that the turbulent velocity dispersion flattens out near the core center, forming a “coherent” region <cit.>, this could be due to projection effects (see <ref> for related discussion).].
In order to take into account the effect of turbulence, <cit.> constructed a theoretical model of isothermal spheres supported by non-thermal pressure using a phenomenological, “logotropic” equation of state P/P_c = 1 + A log (ρ/ρ_c), where P_c and ρ_c are the central pressure and density, respectively, and A is a free dimensionless parameter.
However, such an equation of state is not derived from first principles, and the linewidth-size relation that it predicts is inconsistent with observations except for a limited range of radius, as pointed out by <cit.>.
Other models simply replaced the pressure term in the hydrostatic equation with assumed nonthermal pressure <cit.> or employed a composite polytrope <cit.>.


For application to cores within molecular clouds, the
BE sphere is also problematic in that it assumes a definite outer boundary at which the density profile is sharply truncated, which is in stark contrast to real cores that continuously merge into the ambient surrounding cloud
[Cores embedded in much warmer gas can in principle be truncated by the thermal pressure of hotter phase gas, but this is not a realistic star formation environment.].
For example, it has been pointed out that the critical center-to-edge density contrast of 14 is too small compared to the large dynamic range of volume density in GMC, such that it would be extremely difficult for a core to remain stable <cit.>.
However, this objection does not account for the fact that a core never exists in isolation: the gravitational field away from the center of the core is increasingly governed by surrounding structures (dense filaments, other cores, and stars).
Thus, it is entirely possible that there exists an effective outer “tidal boundary” determined by the structure of gravitational potential around a core.




The primary goal of this paper is to construct a physical model of a dense core supported by both thermal and turbulent pressure.
Instead of assuming a relationship between pressure and density, we directly solve the angle-averaged equations of hydrodynamics assuming a statistical quasi-equilibrium and a power-law linewidth-size relationship.


<ref> visually illustrates the basic concept of our TES model, juxtaposed with the traditional BE sphere with no turbulent motion.
The results presented here quantify how the critical density contrast, radius, and mass change with the properties of underlying turbulent velocity field.
The BE sphere naturally appears as a limiting solution with vanishing turbulent velocities.

In a companion paper (Paper II), we will show that quasi-equilibrium structures resembling the TES emerge in three-dimensional simulations of self-gravitating, isothermal turbulence with conditions similar to star-forming GMCs.
The simulations presented in Paper II demonstrate that the size of the local potential well associated with a core is limited by the gravity from neighboring structures, imposing an effective maximum core radius.
This allows stable equilibria to exist in an isothermal medium with continuous density distribution.
Paper II will analyze evolution of individual cores and gravitational potential structure around them, to identify critical conditions that determine the onset of the collapse.


The remainder of this paper is organized as follows.
In <ref>, we derive the equation of motion governing the radial dynamics of a turbulent core.
In <ref>, we introduce dimensionless variables and derive equilibrium solutions.
In <ref>, we perform a stability analysis to obtain the critical radius beyond which the equilibrium becomes unstable. Equilibrium solutions truncated at this radius are termed “critical cores.”
<ref> presents the physical properties of critical cores as a function of the average velocity dispersion.
Finally, we summarize our work and discuss its implications in <ref>.




§ ANGLE-AVERAGED EQUATION OF MOTION


In this section, we derive the angle-averaged Lagrangian equation of motion (<ref>) satisfied by a spherical, isothermal region pervaded by turbulence.
We start by writing the continuity equation and the radial component of the momentum equation in  spherical coordinates as follows:

    ∂ρ/∂ t
      + 1/r^2∂/∂ r( r^2ρ v_r )
      + 1/rsinθ∂/∂θ( sinθρ v_θ)
    
      + 1/rsinθ∂/∂ϕ( ρ v_ϕ)
      = 0,


    ∂( ρ v_r)/∂ t
      + 1/r^2∂/∂ r( r^2 ρ v_r^2 + r^2 ρ c_s^2 )
    
      + 1/r sinθ∂/∂θ(ρ v_r v_θsinθ)
      + 1/rsinθ∂/∂ϕ(ρ v_r v_ϕ)
    
      - ρ v_θ^2 + ρ v_ϕ^2 + 2ρ c_s^2/r = ρ g_r.

Here, ρ is the gas density, v_r, v_θ, and v_ϕ are the radial, meridional, and azimuthal components of the gas velocity, and g_r is the radial component of the gravitational acceleration.
Integrating <ref> over a full solid angle and dividing by 4π leads to the angle-averaged continuity and momentum equation,

    ∂<ρ>/∂ t + 1/r^2∂/∂ r( r^2 <ρ v_r > ) = 0,


    ∂<ρ v_r >/∂ t
      + 1/r^2∂/∂ r( r^2 <ρ v_r^2 > + r^2 <ρ c_s^2 > )
    
      - <ρ v_θ^2 > + <ρ v_ϕ^2 > + 2 <ρ c_s^2 >/r
      = <ρ g_r >,

in which angle brackets denote the averaging operation

    <Q >≡1/4π∫_0^2π∫_0^π Q sinθ dθ dϕ

for any physical quantity Q.
It is useful to define a related operation

    <Q >_ρ≡∫_0^2π∫_0^πρ Q sinθ dθ dϕ/∫_0^2π∫_0^πρsinθ dθ dϕ = <ρ Q>/<ρ>

which is the mass-weighted angle-average of Q.

Without loss of generality, we decompose the velocity fields into mean and turbulent components,

    v_r    = <v_r >_ρ + δ v_r,
    
      v_θ   = <v_θ>_ρ + δ v_θ,
    
      v_ϕ   = <v_ϕ>_ρ + δ v_ϕ.,

such that <δ v_r >_ρ = <δ v_θ>_ρ = <δ v_ϕ>_ρ = 0 by definition.
Using <ref>, one can recast <ref> into a Lagrangian equation of motion 

    <D <v_r >_ρ/Dt>_ρ
      = f_thm + f_trb + f_grv + f_cen + f_ani

where

    <D <v_r >_ρ/Dt>_ρ
        = ∂<v_r >_ρ/∂ t + <v_r >_ρ∂<v_r >_ρ/∂ r

is the radial acceleration.
The individual force components (thermal, turbulent, gravitational, centrifugal, and anisotropic terms) are given by

    f_thm   = -1/<ρ>∂ P_thm/∂ r,
    
      f_trb   = -1/<ρ>∂ P_trb/∂ r,
    
      f_grv   = < g_r >_ρ,
    
      f_cen   = <v_θ>_ρ^2 + <v_ϕ>_ρ^2/r,
    
      f_ani   = <δ v_θ^2 >_ρ + <δ v_ϕ^2 >_ρ - 2 <δ v_r^2 >_ρ/r,

in which the thermal and turbulent pressures are defined by

    P_thm   ≡<ρ>c_s^2,
    
      P_trb   ≡<ρ><δ v_r^2 >_ρ,

respectively.
In what follows, we assume the centrifugal force is negligible compared to other forces and the turbulence is statistically isotropic, such that f_cen = f_ani = 0.







§ FAMILY OF EQUILIBRIA



We now consider a roughly spherical, radially stratified region, and seek a solution in which thermal and turbulent pressure gradient forces balance self-gravity at every radius, i.e.

    f_thm + f_trb = -f_grv.

Within a region relatively far from surrounding gravitating masses, the gravitational force can be approximated by

    f_grv≈ -GM_enc(r)/r^2,

where

    M_enc(r) ≡ 4π∫_0^r  r'^2⟨ρ⟩ dr'

is the enclosed mass within the radius r.
Using <ref>, <ref> can be written as

    1/r^2∂/∂ r( r^2/<ρ>∂ P_eff/∂ r) = -4 π G <ρ>

where 

    P_eff≡ P_thm + P_trb

is the total effective pressure.
Motivated by observations, we further assume that the turbulent velocity dispersion increases with radius as a power law,

    <δ v_r^2 >_ρ^1/2 = c_s ( r/r_s)^p,

in which r_s is the sonic radius and p is the power law index, such that P_trb = P_thm(r/r_s)^2p.


To recast <ref> in a dimensionless form, we define a dimensionless radial coordinate ξ and logarithmic density contrast u by[See Appendix <ref> for an alternative dimensionless formulation based on given external pressure.]

    ξ   ≡(4π Gρ_c)^1/2/c_sr,
    
      u(ξ;ξ_s,p)    ≡lnρ_c/<ρ>,

where ρ_c is the density at the center (i.e., r=0).
In <ref>, it is made explicit that u is a function of ξ and involves two dimensionless parameters ξ_s = r_s(4π G ρ_c)^1/2/c_s and p.
Here, ξ_s is the dimensionless sonic radius defined by r=r_s in <ref>.
We also define the dimensionless enclosed mass

    m(ξ;ξ_s,p) = (4π G^3ρ_c)^1/2/c_s^3M_enc

and the function

    χ(ξ;ξ_s,p) ≡ 1 + (ξ/ξ_s)^2p

that relates thermal pressure to total pressure through P_eff = χ P_thm.
In terms of these dimensionless variables, <ref> becomes

    1/ξ^2∂/∂ξ[ ξ^2 (χ∂ u/∂ξ - ∂χ/∂ξ) ] = e^-u.

It is useful to note that the term in the square bracket is identical to the dimensionless enclosed mass, i.e.,

    m(ξ;ξ_s,p) ≡ξ^2 (χ∂ u/∂ξ - ∂χ/∂ξ).

We remind the reader that <ref> reduces to the usual isothermal Lane-Emden equation <cit.> in the limit of vanishing turbulent pressure.




For a solution to have a finite central density (i.e., a plateau), it must satisfy the boundary conditions

    u|_ξ=0   = 0,
    ∂ u/∂ξ|_ξ=0   = u'_0,

with a finite value of u'_0 which does not necessarily equal zero.
For example, for p=0.5, one can find a series solution of <ref> near ξ=0,

    u = 1/ξ_sξ + ( 1/6 - 1/2ξ_s^2) ξ^2 + ( 1/3ξ_s^3 - 7/36ξ_s) ξ^3 + ⋯,

which yields u'_0 = ξ_s^-1.


In practice, one can integrate <ref> in terms of the logarithmic radius ϖ≡lnξ, with the initial conditions u=0 and ∂ u/∂ϖ=ξ (∂ u/∂ξ) =0 for very small ξ.
The resulting solutions form a two-parameter family characterized by ξ_s and p.
We term this two-parameter family of solutions TES.
A python package that calculates these solutions is available in <cit.>.
[<https://github.com/sanghyukmoon/turbulent_equilibrium_sphere>.]




<ref> plots the radial density profiles of the TES having p=0.5 and ξ_s=∞, 7, and 3.
The profile with ξ_s=∞ is identical to that of the BE sphere.
For smaller values of ξ_s, the density profile in the inner part is steeper than the BE solution, while it is shallower in the outer part. 
In the innermost regions of cores where gravity is weak and the  local density profile ρ∝ r^-q is shallow (small q), thermal pressure gradients for the TES must be steeper than for the BE sphere in order to compensate for turbulent pressure forces, which are inward when q < 2p.
Conversely, in the outer regions of cores where the BE solution requires a steep density gradient for thermal pressure to balance gravity, the density gradient can be shallower for the TES because turbulent pressure support assists thermal pressure support provided q > 2p.






§ STABILITY



We assess the stability of the TES using an analysis similar to <cit.>.
We consider a sphere of radius r and imagine that the volume enclosed is slightly compressed.
In reality, this compression can be driven by either inflows from larger radii or random turbulent motions.
As the gas distribution interior to this Lagrangian boundary surface adjusts to a new equilibrium, the effective pressure at the boundary will change according to the normalized bulk modulus defined by

    K ≡ -(∂ln P_eff/∂ln V_r)_M_enc,c_s,r_s,p

where

    V_r ≡4π r^3/3

is the spherical volume enclosed within the radius r.
Note that we keep the sonic radius r_s and index p constant in <ref>, assuming that the background turbulent flow is not affected by the perturbations applied.[This implies that the turbulent kinetic energy contained within the Lagrangian radius (i.e., the radius containing constant M_enc) decreases under the compression. Although turbulence can in principle be amplified when the compression timescale is shorter than the flow crossing time <cit.>, in practice we do not expect that the overall compression speed of turbulent cores will exceed their velocity dispersion, unless they are already undergoing gravitational collapse.]





To calculate the bulk modulus K, we first write the effective pressure, volume, and mass as functions of ξ and ξ_s as follows:

    P_eff   =ρ_c c_s^2 χ e^-u= c_s^4/4π G r_s^2ξ_s^2χ e^-u,
    
      V_r    = 4π/3r^3=4π r_s^3/3ξ_s^-3ξ^3 ,
    
      M_enc   =c_s^3/(4π G^3 ρ_c)^1/2m= c_s^2 r_s/Gξ_s^-1 m .

We remind the reader that χ, u, and m appearing in <ref> are dimensionless functions of the variable ξ and the parameters ξ_s and p.


The total derivatives of ln P_eff, ln V_r, and ln M_enc at a fixed c_s, r_s, and p are

    δln P_eff   = ( ∂lnχ/∂lnξ - ∂ u/∂lnξ) δlnξ
          + ( 2+ ∂lnχ/∂lnξ_s - ∂ u/∂lnξ_s) δlnξ_s
       = -m/χξδlnξ
          + (2+ ∂lnχ/∂lnξ_s  - ∂ u/∂lnξ_s) δlnξ_s,
      
    δln V_r    = 3 δlnξ - 3 δlnξ_s,
    δln M_enc   = ∂ln m/∂lnξδlnξ + ( ∂ln m/∂lnξ_s - 1 ) δlnξ_s,

where we have used <ref> in the second equality of <ref>.
Because we are interested in the bulk modulus of a region of constant M_enc, <ref> requires

    δlnξ/δlnξ_s = m/ξ^3 e^-u( 1 - ∂ln m/∂lnξ_s)

to be satisfied.
Using <ref>, after some algebra, one obtains

    K    = 2/31 - 1/2∂ u/∂lnξ_s + 1/2∂lnχ/∂lnξ_s - m^2/2χ e^-uξ^4( 1 - ∂ln m/∂lnξ_s) /1 - m/ξ^3e^-u( 1 - ∂ln m/∂lnξ_s) 
       = 2/31 - 1/2∂ u/∂lnξ_s + 1/2∂lnχ/∂lnξ_s - ( 4π/3)^1/3GM_enc^2/6P_effV_r^4/3( 1 - ∂ln m/∂lnξ_s) /1 - M_encc_s^2/3P_thmV_r( 1 - ∂ln m/∂lnξ_s).

In the limit of negligible turbulence (ξ_s→∞), each of the derivatives with respect to ξ_s in <ref> may be set to zero, and we recover Equation 2.16 of <cit.> for an isothermal sphere.
In the simultaneous limit ξ_s→∞ and G→ 0, the density is uniform so that M_encc_s^2/(P_thmV_r) = 1, and one recovers the ideal gas K=1.
The equilibrium is unstable when K < 0, because a slight compression leads to a further decrease in the interior pressure.


<ref> plots density profiles for selected TES with p=0.5 and ξ_s=∞, 7, and 3, together with the radial profile of K in each case.
In all cases, for ξ≪ 1 the solutions have K≈ 1 because there is not enough mass to be self-gravitating and the thermal pressure dominates the turbulent pressure.
As ξ increases, however, self-gravity becomes more and more important and K decreases.
We define the critical radius ξ_crit as the radius where K first becomes negative. 
For ξ>ξ_crit, the magnitude of K rapidly increases, undergoing another sign change at ξ = ξ_crit,2 back to positive.
A TES with outer radius ξ > ξ_crit is unstable because a part of its interior has K < 0.





For each member of the family of solutions with given p and ξ_s, the TES with outer radius ξ = ξ_crit is identified as the critical solution.  
We further define the critical logarithmic density contrast u_crit≡ u(ξ_crit)=ln[ρ_c/ρ(ξ_crit)] and the critical mass m_crit≡ m(ξ_crit) associated with the critical radius ξ_crit.
<ref>(a)–(d) shows the parametric dependence of ξ_crit, m_crit, u_crit, and the average turbulent velocity dispersion σ_1D on ξ_s, for p=0.3, 0.5, and 0.7.  
Here, we define σ_1D using the mass-weighted average over the whole core within r_crit, assuming the turbulence is statistically isotropic:


    σ_1D   ≡(∭_r<r_critρδ v_r^2 dV/∭_r<r_critρ dV)^1/2
       = c_s [∫_0^ξ_crit e^-u (ξ/ξ_s)^2pξ^2 dξ/∫_0^ξ_crit e^-uξ^2 dξ]^1/2.









In the limit of ξ_s→∞, we find ξ_crit=6.45, m_crit=15.7, and u_crit=2.64, identical to the well-known result for the critical BE sphere <cit.>.

The point ξ_s = ξ_crit might be of some interest, because the sonic radius is within the core for smaller ξ_s, while it is outside the core for larger ξ_s.
For p=0.3, 0.5, and 0.7, ξ_s=ξ_crit occurs at ξ_s = 10.2, 8.99, and 8.21.
It is interesting to note that σ_1D≈ c_s at ξ_s ≈ 6, quite insensitive to p.
It is worth mentioning, as well, that while this and other curves plotted in <ref> give an impression that they all cross at the same locus, a detailed examination reveals that it is only approximately true, and the crossing points differ for different curves (e.g. ξ_crit≈ 14 when ξ_s ≈ 5, and m_crit≈ 65 when ξ_s ≈ 4).
The similar crossing points arise because the curves with different p (for panels a-c) go to the same limit ξ_crit→ξ_BE at ξ_s→∞, and more generally these solutions must be quite similar in the regime where ξ_s>ξ_crit because χ in <ref> is close to unity.

For p=0.5, <ref> provides tabulated values for ξ_s, ξ_crit, m_crit, u_crit, and σ_1D.
ccccc
  Properties of critical TES with p=0.5.
  
  
ξ_s    
ξ_crit    
m_crit    
u_crit    
σ_1D
  
  ∞     6.45     15.7     2.64     0.0 

  18.2     7.50     19.9     2.81     0.5 

  6.42     10.6     34.5     3.25     1.0 

  4.21     15.8     64.3     3.76     1.5 

  3.42     22.7     115     4.24     2.0 

  2.55     103     1.4× 10^3     6.17     5.0 

  2.44     385     1.3× 10^4     7.79     10 

  2.42     ∞     ∞     ∞     ∞
  <ref> is published in its entirety in machine-readable format. A few representative rows are shown here for guidance regarding its form and content.


We note that ξ_crit, m_crit, and u_crit can be translated to the dimensional critical radius r_crit, mass M_crit, and critical density contrast (ρ_c / ρ_e)_crit,  through <ref>.
For radius and mass, conversion between dimensionless and dimensional variables is given by:

    r =  ξ× 0.014  ( T/10 K)^1/2( n_H,c/10^5 cm^-3)^-1/2


    M =  m× 0.11  M_⊙( T/10 K)^3/2( n_H,c/10^5 cm^-3)^-1/2

where n_H,c is the central density and T is the temperature.




<ref> shows that the dimensionless critical radius ξ_crit steeply increases with decreasing ξ_s.
This indicates that a core forming in a region of strong local turbulence would initially be stable: large ξ_crit would make the critical radius lie beyond the effective tidal radius that limits the core, as imposed by the gravity of nearby structures in the GMC
However, evolution generally occurs in the direction of increasing ξ_s ∝ρ_c^1/2 r_s as the initial converging flows further compress the core and turbulence dissipates, leading to a decrease in ξ_crit.
When ξ_crit becomes small enough, the core will become unstable, triggering collapse.
More details on the collapse scenario will be discussed in <ref>; the quantitative test of this scenario using numerical simulations will be presented in Paper II.


We empirically find that, for p ≳ 0.4, there exists a minimum dimensionless sonic radius ξ_s,min below which K>0 at every radius.
In <ref>, ξ_s,min for a given p corresponds to the value at which ξ_crit→∞. 
The existence of ξ_s,min for a TES family means that for a given central density, every member of the family is stable if the turbulence is sufficiently strong that r_s < ξ_s,min c_s/ (4π G ρ_c)^1/2, where <ref> can be used to translate this into physical terms.

Alternatively, because ξ_s ∝ρ_c^1/2 r_s (<ref>), the existence of ξ_s,min implies that for a given dimensional sonic radius r_s, a TES with the central density lower than ρ_c,min≡ξ_s,min^2c_s^2/(4π G r_s^2) is always stable.
Let us suppose that the larger-scale environment of the cores is a spherical cloud of radius R_cloud, mass M_ cloud, and one-dimensional Mach number ℳ_1D = (R_ cloud/r_s)^p, with virial parameter

    α_vir,cloud≡5 ℳ_1D^2c_s^2 R_ cloud/G M_ cloud.

The minimum central density for instability to be possible can then be written as

    ρ_c,min/ρ_0 = ξ_s,min^2/15α_vir,cloudℳ_1D^2/p - 2,

where ρ_0 = M / (4π R^3/3) is the mean density of the cloud.


For turbulent power law index p=0.5, we find ξ_s,min = 2.42, which would correspond to a minimum sonic scale of r_s≈ 0.03 using the fiducial central density and temperature in <ref>.
Alternatively, the minimum molecular hydrogen number density at the core center, n_H2,c^min=n_ H,c^ min/2 ≡ρ_c,min / (2μ_Hm_H), that would allow for instability is given in terms of the physical sonic scale r_s by

    n_H2,c^min = 2.2× 10^4 cm^-3 (ξ_s,min/2.42)^2( T/10 K) ( r_s/0.05 pc)^-2.

We emphasize that these are the local values of T and r_s, the latter of which can vary considerably within a GMC (see Paper II).
With typical molecular cloud virial parameter α_vir,cloud∼ 2-4 on ∼ 100 scale <cit.>,  <ref> implies that the central density in the core would have to exceed the ambient cloud density by at least a factor ρ_c,min/ρ_0∼ M^2 for instability to be possible, when p=0.5. 


We stress that the condition ξ_s> ξ_s, min is a necessary but not a sufficient condition for collapse.
That is, ρ_c > ρ_c,min  does not guarantee instability and hence should not be interpreted as a critical density for collapse.
Instead, ρ_c > ρ_c,min simply means that there exists a critical radius beyond which a quasi-equilibrium is unstable.
For a core to collapse, its central density must exceed the minimum value, and the core's outer radius and total mass must exceed the critical values for TES solutions, as shown in <ref>.






For p=0.5, <ref> plots the critical radius and mass as a function of ρ_c / ρ_c,min=(ξ_s/ξ_s,min)^2.
We note that cases with large (small) ρ_c/ρ_c,min correspond to weak (strong) turbulence. 
A highly turbulent core would therefore have to be larger and more massive to become unstable, compared to quiescent cores.
Notably, <ref> shows that when ρ_c is close to ρ_c,min, r_crit is at least a factor of a few larger than the BE radius at that ρ_c, R_BE=1.82 R_G(ρ_c).
<ref> also shows that the mass begins to exceed M_BE = 4.43 M_G(ρ_c) by a factor two or more when ρ_c/ρ_c,min < 9.
At the same time, <ref> shows that for a given r_s (and therefore given ρ_c,min), the central density must become quite large in order for the unstable core mass to approach M_BE: ρ_c/ρ_c,min has to be at least 21 in order to have a core with M_crit/M_BE < 1.5.
Since real cores cannot have either too large a radius (this would exceed the cloud size) or too high a central density (difficult to reach with realistic larger-scale dynamics), the majority of cores that become unstable may have a moderate range of r_crit/R_BE and M_crit/M_BE, or equivalently, of σ_1D/c_s (see <ref>).
A comparison of gravitational and flow crossing timescales also suggests a highly turbulent core with σ_1D≫ c_s would be difficult to collapse as a whole (see below).
Although ξ_s,min is a strict lower limit for allowing instability, we expect collapse would mostly occur at ξ_s ≳ 6, for which turbulent velocity dispersion is at most transonic regardless of the power law index p (see <ref>(d)).
We further discuss scenarios for collapse in <ref>.




§ PHYSICAL PROPERTIES OF CRITICAL TES



While the structure of the TES is fully characterized by ξ_s and p, it is often useful to reparameterize the results in terms of the core-averaged one-dimensional velocity dispersion σ_1D defined in <ref>.
Here we will consider just the case of critical cores.
For p=0.5, <ref>(a) plots the density profiles of critical TES having different values of σ_1D/c_s, normalizing relative to the density at the outer radius.
Increased turbulent support makes the core larger and more centrally concentrated; a trans-sonic core has a factor of two higher density contrast than a purely thermal core, while increasing the internal Mach number to 2 leads to a factor ∼ 5 increase in the density contrast. Here we show just cases with p=0.5, but other indices produce similar profiles for low σ_1D since thermal pressure dominates.





<ref>(a) suggests that a TES develops a power law envelope as σ_1D increases.
To illustrate this more clearly, in <ref>(b) we plot the density profiles of highly turbulent TES with σ_1D = 10 c_s for selected linewidth-size indices p=0.3, p=0.5, and p=0.7.
We then calculate the local density slope q ≡ -∂lnρ/∂ln r for all cases shown in <ref>(a)-(b) and plot the profiles of q in <ref>(c).
While no density profile is described by a single power law, the outer density slopes generally lie in the typical range q ∼ 1–2 found in observations <cit.>.


<ref>(c) suggests that for σ_1D/c_s≳ 2, q becomes flat in the outer part of the core.
In fact, if in the limit σ_1D≫ c_s one substitutes ρ∝ r^-q into <ref>, one can show that the outer envelope of the TES becomes a power law with q → 2 - 2p when 0< p < 0.5 and q → 2p when 0.5 ≤ p < 1.
The corresponding limits in the polytropic index γ_p ≡∂ln P_eff / ∂lnρ = 1 - 2p/q are γ_p → (1 - 2p) / (1 - p) when 0< p < 0.5 and γ_p → 0 when 0.5 ≤ p < 1.
The former limit is equivalent to the SPS ρ∝ r^-q and P_eff∝ρ^γ_p with q = 2/(2-γ_p) <cit.>: the outer envelope of a highly turbulent TES with 0 < p < 0.5 approaches the SPS solution.
For 0.5 ≤ p < 1, a singular polytropic solution to <ref> does not exist.
We note that the convergence of q → 1 when p=0.5 is extremely slow: the local slope of the TES solution manages to reach q = 1.05 at ξ_crit = 10^9.
Because of this, the outer envelope of the p=0.5 TES shown in <ref>(c) has not yet reached the power-law limit of q→ 1, while the solutions with p=0.3 and 0.7 almost reach the limiting value q → 1.4.


For all these solutions, equilibrium is achieved by balancing the three forces: thermal pressure gradient (f_thm), turbulent pressure gradient (f_trb), and gravity (f_grv) (e.g., <ref>).
While f_thm is always positive and f_grv negative, the direction of f_trb can be either outward or inward depending on the sign of 2p - q, because locally P_trb∝ r^2p - q (see the last paragraph of <ref>).
<ref>(d) shows the relative importance of gravity compared to the other two forces in determining the equilibrium, by plotting the quantity

    ℛ_grv≡|f_grv|/max(|f_thm|, |f_trb|)

as a function of radius.
It shows that for all solutions (except the BE sphere
where f_trb = 0), ℛ_grv≪ 1 in the inner part: there, the thermal pressure gradient is balancing the inward turbulent pressure gradient, with gravity playing only a minor role.
However, as q increases toward the outer part, the turbulent pressure gradient force decreases in magnitude, and gravity starts to take over.



As mentioned in <ref>, a previous theoretical model that was introduced in order to represent nonthermal support was the logotrope <cit.>.   
<ref> compares the internal structure of the TES with that of the logotrope having the identical velocity dispersion.
It shows that compared to the TES, the logotrope has a significantly flatter outer density profile.
<ref> also compares the linewidth-size relation for the TES with that of the logotrope, where the unphysical turnover of the logotrope is evident at large scales.  












Using <ref>, the critical mass and radius of a TES can be expressed in terms of its center, edge, or mean densities as follows:

    M_crit   = m_crit/(4π)^1/2M_G(ρ_c),
       = m_crite^-u_crit/2/(4π)^1/2M_G(ρ_e),
       = (3/4π)^1/2( m_crit/ξ_crit)^3/2M_G(ρ),


    r_crit   = ξ_crit/(4π)^1/2R_G(ρ_c),
       = ξ_crite^-u_crit/2/(4π)^1/2R_G(ρ_e),
       = (3/4π)^1/2( m_crit/ξ_crit)^1/2R_G(ρ),

where ρ_e ≡<ρ>(r = r_crit) is the edge density, ρ≡ 3M_crit/(4π r_crit^3)=3ρ_c m_crit/ξ_crit^3 is the mean density, and M_G and R_G are the gravitational mass and radius defined in <ref>.
We remind the reader that u_crit, ξ_crit and m_crit are dimensionless functions of p and ξ_s, or alternatively, of p and σ_1D (see <ref>(d) for the relationship between ξ_s and σ_1D).
We note that <ref> reduce to <ref> in the limit of σ_1D→ 0.


<ref>(a),(b) plots M_crit/M_G(ρ) and r_crit/R_G(ρ) as functions of σ_1D/c_s, showing that both quantities increase with σ_1D/c_s, rather slowly when σ_1D≲ c_s and rapidly for σ_1D≳ c_s.
For p=0.5 and σ_1D = c_s, we find r_crit = 0.88 R_G(ρ) and M_crit = 2.86 M_G(ρ), which are factors of 1.15 and 1.54 larger than the corresponding values for the BE sphere having the same mean density (see also the profile comparison in <ref>).[We note that the difference between the TES and BE sphere is more pronounced when compared at the same central density (see <ref>).]
<ref>(a),(b) also plot the “turbulent BE” mass and radius resulting from a naive substitution c_s→ (c_s^2 + σ_1D^2)^1/2 in <ref>, which significantly overestimate the true critical mass and radius depending on the Mach number.
For example, the “turbulent BE” mass and radius
of mildly supersonic cores (σ_1D = 2c_s) are larger than M_crit and R_crit by a factor of 3.7 and 1.5 (assuming p=0.5), respectively, indicating such simple substitution would lead to incorrect results.


When σ_1D≫ c_s, the critical radius and mass as well as the critical density contrast all increase as a power-law in σ_1D/c_s.
For example, we find that the critical quantities for p=0.5 are well approximated by

    M_crit   ≈ M_BE(ρ) (1 + σ_1D^2/2c_s^2) 
    
                        
                           ≈ 30 M_⊙( T/20 K)^3/2( n_H/500 cm^-3)^-1/2
         ×[ 1 + 64(σ_1D/3 km s^-1)^2(T/20 K)^-1],


    r_crit   ≈ R_BE(ρ)( 1 + σ_1D^2/2c_s^2)^1/3
    
                        
                           ≈ 0.74 pc( T/20 K)^1/2( n_H/500 cm^-3)^-1/2
         ×[ 1 + 64(σ_1D/3 km s^-1)^2(T/20 K)^-1]^1/3,


    (ρ_c/ρ_e)_crit≡ e^u_crit≈ 14( 1 + 0.7σ_1D^2/c_s^2)^1.2

which are plotted in <ref>(a)–(c) with blue dashed lines.
We note that <ref> are valid within relative error of 5% for σ_1D < 9.5 c_s, σ_1D<13c_s, σ_1D<31c_s respectively.
Intriguingly, by evaluating  <ref> at σ_1D≳ 3 km s^-1, it is readily evident that the critical mass at high Mach number σ_1D≫ c_s becomes comparable to typical masses of star clusters, as will be discussed in <ref>.


<ref>(c) plots the center-to-edge density contrast for marginally stable equilibrium solutions, as a function of σ_1D/c_s.
The critical center-to-edge density contrast converges to the well-known value of 14 in the limit σ_1D→ 0, but then monotonically increases with σ_1D/c_s with very weak dependence on p, reaching 26 and 70 at σ_1D = c_s and 2c_s, respectively.
We note that for the equilibrium with larger p, the outer density profile has to have a steeper slope to generate outward pressure gradient, leading to smaller r_crit and M_crit despite having similar center-to-edge density contrast.

<ref>(c) also shows the ratio of the average density to the edge density (right axis). 
The average density is a factor of ∼ 1.5 - 2.5 larger than the edge density, with a mild dependence on σ_1D/c_s and p, and the largest ratio for the limiting non-turbulent case.
The ratio ρ̅/ρ_e is only ∼ 1.5-2 for the highly turbulent cores that have extreme center-to-edge density contrast ρ_c/ρ_e > 10^2; while central densities are large, very little mass is involved.


<ref>(d) plots the ratio of the gravitational free-fall time t_ff,avg≡[3π/(32 Gρ)]^1/2 at the mean core density to the effective crossing time t_cr≡ r_crit/(c_s^2 + σ_1D^2)^1/2.
The large t_ff,avg/t_cr for σ_1D≫ c_s and small p may indicate that turbulence would quickly rearrange fluid elements in the outer part of the core before gravity can bring about collapse.
For p<0.5, however, t_ff,avg/t_cr∼ 0.7 - 2, so turbulence would be less likely to disrupt cores before the structure as a whole collapses.





Observations with limited spatial resolutions are able to constrain ρ and therefore t_ff,avg.
However, the characteristic timescale of the evolution leading to the formation of a singularity (i.e., protostar) at the center is t_ff,c≡ [3π / (32 G ρ_c)]^1/2, which is shorter than t_ff,avg.
Using the data shown in <ref>(c), we find an analytic approximation

    (ρ_c/ρ)_crit≈ 5.70 ( 1 + 0.8 σ_1D^2/c_s^2)^1.2

analogous to <ref>, which indicates that the ratio

    t_ff,c/t_ff,avg = (ρ_c/ρ)_crit^-1/2≈ 0.42(1 + 0.8σ_1D^2/c_s^2)^-0.6

is a decreasing function of σ_1D/c_s (see <ref>(a)).
We note that t_ff,c is a factor of ∼ 3 times smaller than t_ff,avg for a transonically turbulent core.
Since t_ff,avg/t_cr∼ 1 for transonic cores, this implies that the center would collapse on a timescale shorter than the outer part would evolve either due to gravity or to turbulence-driven restructuring.
This would be all the more true for highly turbulent cores that reach critical conditions, as ρ_c/ρ̅≫ 1 in this limit.




Observational sensitivity limits generally make it difficult to probe the outer envelopes of cores.
It is common to define the core radius using the FWHM of the column density profile, which is known to approximate R_BE if the underlying density distribution is similar to the BE sphere <cit.>.
We introduce the notation R_FWHM to denote this radius, i.e., R_FWHM≡ 2R_HM where the column density at R_HM is the half of the central value.
To relate R_FWHM to the critical radius for cores having nonzero velocity dispersion, we first 
integrate to obtain 2D maps of TES, thereby obtaining column density profiles.
<ref> compares examples of the column density profiles for the critical TES having zero, transonic, and mildly supersonic turbulent velocity dispersions (see <ref> for the corresponding volume density profiles normalized to the edge density).
For the latter two, we assume the turbulent power-law index of p=0.5.
Based on the column density profiles, we plot the ratio R_FWHM/r_crit in <ref>(b).
This shows that while R_FWHM/r_crit∼ 0.8 when σ_1D≲ 0.5 c_s, this ratio is smaller by a factor of two or more for cores having supersonic turbulent velocity dispersions.
That is, when the turbulence level is high, the radius of the core must extend well beyond R_FWHM for a core to be unstable to collapse.
The dotted curve in <ref> shows an example of this.





§ DISCUSSION




 §.§ Collapse Scenario


In the idealized theory, a quasi-equilibrium core supported by thermal and turbulent pressure is unstable if its maximum radius r_max exceeds r_crit.
However, the TES model describes a single isolated object, without regard to any gravitational forces other than those produced by the core itself.
In reality, hierarchical structure around the core contributes to the landscape of the gravitational potential, and r_max must also be effectively tidally limited.
For example, in the case of two cores with equal mass, the effective tidal limit is half of their separation.
More generally, if we take r_max to be the distance to the nearest saddle point of the gravitational potential, it scales linearly with the distance to an external structure of mass M_external, and inversely with 1 + (M_external/M_core)^1/2.
While making a quantitative prediction of r_max is beyond the scope of this work, it is clear that the spatial extent of any realistic core forming inside a turbulent cloud is limited by the scale associated with the velocity perturbation creating the core.
In observations, transformation of column density maps of molecular clouds into a dendrogram <cit.> can provide one way of characterizing r_max.


<ref> schematically illustrates two alternative evolutionary scenarios leading to collapse or dispersal.
An overall assumption is that a core is formed by turbulent flows that are locally converging, and the density of the core is determined by the strength of the converging motion.
Larger-amplitude and longer-duration converging velocity perturbations will create higher local density contrast relative to ambient values.
At the same time, turbulence on scales smaller than the overall converging flow will help to protect against collapse by providing support against gravity.


If the converging flow is strong (creating large ρ_c) and/or the smaller-scale turbulence is weak (i.e., large r_s), the resulting core would have large enough ξ_s such that r_crit < r_max.
The core would then subsequently undergo collapse leading to star formation.
In contrast, when the converging flow is too weak or turbulence is too strong such that r_crit > r_max, the entire region within r_max is stable and the core would disperse unless it is further compressed by, e.g., passing shock waves.


In Paper II, we will test the above scenario by directly comparing r_crit and r_max for cores forming in a large suite of numerical simulations.
By tracking cores over time, we show that the net radial force within each core becomes negative (i.e., accelerating the inward bulk motions) when r_crit from TES theory first falls below r_max computed from the full gravitational potential.
We also show that the internal density profiles of the simulated cores at the time they initiate gravitational collapse are consistent with the TES solutions obtained in this paper.
We will connect the timescales of different stages of evolution to observed core lifetimes in future work. 






 §.§ Projection Effect on Linewidth-Size Relation


A basic assumption of the TES model is that the turbulent velocity dispersion increases with radius as a power law (<ref>).
Although this assumption is motivated by several observational studies <cit.>, a direct comparison between <ref> and the observed linewidth-size relation is difficult because the latter involves the integration along the line-of-sight, weighted by the density and emissivity of the transition.
To explore how the intrinsic turbulent velocity structures within a core are related to the observable line-of-sight velocity dispersion, it is necessary to perform numerical simulations in which prestellar cores form within the context of a larger scale GMC, which is done in Paper II.
Our preliminary analyses from those simulations show that the line-of-sight non-thermal velocity dispersion traced by dense gas is approximately constant within a core, reminiscent of the so-called “coherent cores” <cit.>, even though the intrinsic velocity dispersion obeys a power-law.
Careful modeling including density selection effect of different molecular tracers would be required to relate the observed linewidth-size relation to the intrinsic turbulent structures in three dimensions.




 §.§ Implication for Star Cluster Formation


Massive star clusters such as those found within Carina or the Orion Nebula contain a few thousand solar masses within a few parsecs and have velocity dispersions of a few kilometers per second <cit.>.
Given that the thermal Jeans mass under  GMC conditions on large scales is ∼ 10-10^2 M_⊙, a longstanding question is how material can avoid early fragmentation prior to assembling the massive, compact structures that are cluster progenitors.
Our TES solutions shed some light on this question: provided that the central density remains lower than the value given in <ref> or <ref>, a turbulent overdense patch of a GMC encompassing significant mass could in principle become progressively more concentrated without undergoing runaway collapse. 
Assuming T = 20 K, n_H = 500 cm^-3, and σ_1D = 3 km s^-1, typical of dense clumps in GMCs, <ref> yield M_crit∼ 2× 10^3 M_⊙ and r_crit∼ 3 pc, which are intriguingly similar to typical mass and radius of ONC-like clusters.
The quantitative results from our TES model may therefore be key to realizing the equilibrium cluster formation scenario of <cit.>.



This work was supported in part by grant 510940 from the Simons Foundation to E. C. Ostriker.
Computational resources for this project were provided by Princeton Research Computing, a consortium including PICSciE and OIT at Princeton University.
We are grateful to the anonymous reviewer for constructive comments on the manuscript.



aasjournal





sectionappsec



§ ALTERNATIVE FORMULATION


In this appendix, we provide an alternative nondimensionalization of <ref> based on the external pressure P_ext rather than the central density.
We find this formulation useful when probing the highly turbulent regime where the dependence of the critical quantities on ξ_s becomes extremely sensitive (see <ref>).
In order to distinguish the alternative dimensionless variables from those in the main text, we attach a prime to each variable and define

    ξ'    ≡G^1/2P_ext^1/2/c_s^2r,
    
      u'(ξ';ξ'_s,p)    ≡lnP_eff/P_ext,
    
      m'(ξ';ξ'_s,p)    ≡G^3/2P_ext^1/2/c_s^4 M_enc,
    χ'(ξ';ξ'_s,p)    ≡ 1 + ( ξ'/ξ'_s)^2p.

In terms of these dimensionless variables, <ref> becomes

    1/ξ'^2∂/∂ξ'( χ' ξ'^2 ∂ u'/∂ξ') = - 4π e^u'/χ'.

Starting from very small ξ', one can integrate <ref> in terms of the logarithmic radius t' ≡lnξ' with the initial conditions u' = u_c and ∂ u' / ∂ t = 0, until P_eff matches the external pressure P_ext (i.e., u' = 0) at the outer radius ξ_ext.
With all other parameters fixed, the total mass of the sphere m'(ξ_ext) initially increases with u_c (starting from u_c = 0), and then reaches the maximum when u_c = u_c,crit.
It is not difficult to see that this stationary point is equivalent to the point where K = 0 and thus defines the marginally stable solution.
This also implies that, for a given external pressure, there is a maximum mass above which no equilibrium solution exists.



