Star formation in self-gravitating turbulent fluids

# Star formation in self-gravitating turbulent fluids

Norman Murray11affiliation: Canadian Institute for Theoretical Astrophysics, 60 St.George Street, University of Toronto, Toronto ON M5S 3H8, Canada; murray@cita.utoronto.ca 22affiliation: Canada Research Chair in Astrophysics & Philip Chang33affiliation: Department of Physics, University of Wisconsin-Milwaukee, 1900 E. Kenwood Blvd., Milwauke, WI 53211, USA; chang65@uwm.edu ,
###### Abstract

We present a model of star formation in self-gravitating turbulent gas. We treat the turbulent velocity as a dynamical variable, and assume that it is adiabatically heated by the collapse. The theory predicts the run of density, infall velocity, and turbulent velocity, and the rate of star formation in compact massive gas clouds. The turbulent pressure is dynamically important at all radii, a result of the adiabatic heating. The system evolves toward a coherent spatial structure with a fixed run of density, ; mass flows through this structure onto the central star or star cluster. We define the sphere of influence of the accreted matter by , where is the stellar plus disk mass in the nascent star cluster and is the gas mass inside radius . The density is given by a broken power law with a slope inside and to outside . Both and the infall velocity decrease with decreasing for ; , the size-linewidth relation, with , explaining the observation that Larson’s Law is altered in massive star forming regions. The infall velocity is generally smaller than the turbulent velocity at . For , the infall and turbulent velocities are again similar, and both increase with decreasing as , with a magnitude about half of the free-fall velocity. The accreted (stellar) mass grows super-linearly with time, , with a dimensionless number somewhat less than unity, the clump mass and the free-fall time of the clump. We suggest that small values of can be used as a tracer of convergent collapsing flows.

## 1. Introduction

The importance of massive stars in a human context became apparent when Burbidge et al. (1957) showed that such stars were responsible for the production and distribution of most of the heavy elements that composed the Earth, and which form the building blocks of life.

More recently it has emerged that massive stars are important on a less parochial scale. Star formation on galactic scales is observed to be slow, in the sense that the time to deplete the supply of gas from which stars are made is some fifty times the dynamical time. Quantitatively,

 ˙M∗=ηMgτDYN, (1)

where is a dimensionless constant, and is the dynamical time of a galactic disk of radius with circular velocity (Kennicutt, 1998; Leroy et al., 2008). Massive stars provide the feedback that limits the rate of star formation in galaxies. The feedback comes in the form of radiation pressure, gas pressures in ionized regions, blast waves from supernovae, and cosmic rays. The feedback also drives gas out of galaxies, setting the stellar mass-halo mass relationship for halo masses , e.g., Hopkins et al. (2013); Trujillo-Gomez et al. (2013); Agertz & Kravtsov (2014); Ceverino et al. (2014). Understanding how massive stars form is clearly important to our comprehension of the evolution of the universe, and many of its components.

Stars more massive than sixty solar masses () live only . Since we see stars with masses in excess of outside their natal dense gas clumps, the early accretion rate onto such a star must exceed , where we have assumed that the most massive stars are visible for , hence they can accrete for no more than . In the Milky Way, these massive stars are usually found in compact radius clusters containing up to of stars. The accretion rates for clusters must exceed a few hundredths of a solar mass per year, . For example, Cygnus OB2 contains inside a half light radius (Knödlseder, 2000), yielding . The free-fall time of the natal cluster is

 τff ≡ √3π32Gρ = 1(3×104M⊙M∗)−1/2(R1/26.4pc)3/2(ϵ0.5)Myr,

where all quantities are estimated at the half light radius, and is the fraction of natal gas in a dense cluster-forming region that is converted into stars. The star formation rate per free fall time, the small-scale analog of in eqn (1), is

 ϵff ≡ ˙M∗τffMg=τff3Myrϵ (2) ≈ 0.16(R1/26.4pc)3/2(M∗,1/23.5×104M⊙)1/2(ϵ0.5)1/2.

This estimate of is consistent with many estimates in the literature, for example, Wu et al. (2005), Lada et al. (2010), Heiderman et al. (2010), and Murray (2011), although there is some disagreement, e.g., (Krumholz et al., 2012a).

The fact that is one of the reasons that astronomers believe that stellar feedback is necessary to explain the slow global rate of star formation; if , no feedback would be required. A second motivation for feedback is that simulations on the scale of galaxies, and on the scale of star forming clumps () both find very rapid star formation if no feedback is included.

In this paper we are interested in the physics behind rapid star formation rather than the physics of feedback. The problem arises from the high accretion rates, discussed above, required to make massive stars. The early model of Shu (1977) gave an estimate of the star formation rate. His model assumed that stars formed from hydrostatic cores supported by thermal gas pressure. The accretion rate in his model was independent of time, given by , where is the sound speed in molecular gas, and . Thus the maximum accretion rate was , too small by a factor of ten to explain the origin of massive stars.

The difficulty with the accretion rate was overcome by Myers & Fuller (1992), who noted that cores in massive star forming regions had linewidths that exceeded the thermal line width by factors of several or more. They pointed out that the accretion rate in Shu’s model was limited by the signal speed of the expanding collapse wave, then took advantage of the larger signal speed available if they assumed that the pressure in the fluid was given by

 P(r,t)=ρ(r,t)[c2s+v2T(r)], (3)

where is the local density, is the sound speed, and is the turbulent velocity a distance from the center of the collapse.

The model was further developed by McLaughlin & Pudritz (1997) and McKee & Tan (2003). These models kept some of the features of Shu’s model, including the assumption of a hydrostatic initial core, and the static nature of the turbulent velocity in the equation of state (3). In particular, they assumed that was not affected by the collapse, and that it was independent of time. They did not address the question of how a hydrostatic core was assembled.

Observations in large GMCs, and low mass cores, find

 vT(r)=vT(R)(rR)p (4)

with (Larson, 1981; Myers & Goodman, 1988); as we will discuss further below, in massive star forming regions behaves differently, with and much larger turbulent velocities at a given separation than in low mass star forming regions. The turbulent core models either explicitly assume at turbulent velocity of the form (4), or assumed implicitly that it follows a power law in radius so that .

In a companion paper (Lee et al., 2014) (Paper I) we present three dimensional hydro and MHD simulations which recover the result in the bulk of the gas. However, we also show that in small regions around density peaks, which we identify as cores, .

In this paper we assume that massive stars form in regions of turbulent flows in which the gas is converging, so that there is never any time at which the region is in hydrostatic equilibrium. The infall therefore extends over a much larger region than that envisioned in the hydrostatic turbulent core models.

A more profound difference arises from the fact that we model the turbulent velocity as a dynamical variable, changing in response to the evolution of the system. This requires that we introduce a third dynamical equation, the energy equation. However, we are unable to solve the equations that result. This is not surprising, because numerical solutions of the full set, carried out in three dimensions, show that the motion becomes turbulent.

Suitably chastened, we then make two simplifications. First, we assume spherical symmetry, reducing the problem from three spatial dimensions to one. Second, we introduce a closure for the energy equation, that proposed by Robertson & Goldreich (2012).

We import the notion of the sphere of influence of the star from galactic dynamics, and define the radius at which the mass of gas just equals the (time dependent) mass of the star. With our two assumptions, and the concept of the sphere of influence, we are able to find self-similar solutions that describe the collapse of a self-gravitating turbulent gas. We calculate the run of density , of infall velocity , and of the turbulent velocity . We find the striking result that for , the density approaches a fixed function of radius, independent of time, i.e., . For , , with .

We show that the acceleration due to the pressure gradient always tracks the acceleration of gravity. For , the result is that both the infall velocity and the turbulent velocity increase with decreasing radius, as . But for , the pressure force exceeds that due to gravity. As a result, both and increase slowly with increasing radius, , with . This appears to explain the deviations from Larson’s size-linewidth relation seen in massive star forming regions.

It follows from the fact that the density is independent of time for , and from , that the mass accretion rate increases linearly with time, so that

 M∗=ϕMcl(t−t∗τff)2, (5)

where is the mass of the star forming clump, and is the free fall time of the clump.

This paper is organized as follows. In §2 we present a spherically symmetric one dimensional model for the collapse of a turbulently supported gas cloud consisting of the continuity and momentum equations, and a simple closure for the energy equation. We describe approximate analytic solutions to these equations, and comment on the effects of thermal and magnetic support. In §3 we solve the equations numerically. In §4 we briefly compare observations of infall in the massive star forming regions, with special attention given to G10.6-0.4. We also discuss previous analytic and numerical work. We wrap up with our conclusions in the final section.

## 2. Spherical Gravitational Collapse Models

The governing equations are the continuity equation

 ∂ρ∂t+∇⋅(ρu)=0, (6)

the momentum equation

 ∂ρu∂t+∇⋅ρuu=−∇P+ρg, (7)

and the energy equation

 ∂∂t [ρ(12u2+E)]+∇⋅[ρ(12u2+E)⋅u+uP] (8) =ρg⋅u.

In these equations is the density, is the fluid velocity, is the gas pressure, is the acceleration due to gravity, and is the internal energy of the gas.

In the companion paper (Paper I), we solve these equations via a three dimensional simulation. In this paper, we simplify our model to gain some physical insight into the collapse process. To do so we will simplify the energy equation. We start by noting that simulations of both supersonic and subsonic turbulence have generally found that the properties of the turbulent flow are universal, in the sense that certain relations are found under a wide variety of conditions. In particular, simulations find velocity power spectra that have approximate power law forms, e.g., for subsonic tubulence or for supersonic turbulence. Correspondingly, or .

It is tempting to conclude that the energy equation can be replaced by a simple closure relation specifying the turbulent velocity, namely , where is given by Larson’s law. This is what has been done in most previous analytic work on turbulent collapse.

This temptation should be resisted. It is well known that turbulence decays on a turnover time , and it is believed that in star forming galaxies the energy lost in this decay is replaced by stellar feedback, e.g., from supernovae. One could argue that this leads to a steady state, i.e., that . However, Robertson & Goldreich (2012) showed that compression of turbulent gas also drives turbulence. Turbulent gas in a converging flow will be compressed, which tends to increase , while will also decay due to dissipation. Therefore, we use a closure scheme suggested by Robertson & Goldreich (2012) that captures both the decay and the compressive driving.

In the absence of collapse, undriven turbulence decaying on an eddy turnover time is described by

 12dv2T(r)dt=−ηv3Tr, (9)

with a dimensionless constant of order unity.

Collapse alters the turbulence. Ignoring for the moment the cascade of energy from large scales to small that is responsible for (9), a reduction in the radius from to results in a scaling of the velocity , i.e.,

 (dvTdt)AH=−vTurr, (10)

where is negative for a collapse. Thus the scaling relation suggests that a collapse will tend to increase the turbulent velocity. Robertson & Goldreich refer to this as “adiabatic heating”, and show that such heating occurs in cosmological simulations in which the scale factor decreases with time, mimicking a gravitational collapse.

Combining this driving term with the turbulent damping,

 dvTdt=−(1+ηvTur)vTurr. (11)

We use this equation to replace the energy Equation (8).

### 2.1. Fundamental Equations

In Paper I we show that on parsec and smaller scales, the collapse is, very roughly, spherically symmetric. Motivated by this and our desire to simplify our model, we employ spherical coordinates. Both simulations and observations of star forming regions show the presence of star forming regions that have cylindrical geometry, as well as other regions with spherical geometry, or a combination of both. We have made some initial forays into cylindrical geometry, but for simplicity we stick to spherical geometry in this paper.

The continuity and momentum equations become

 ∂ρ∂t+1r2∂∂r(r2uρ)=0, (12)

and

 ∂ur∂t+ur∂ur∂r+1ρ∂ρ(v2T(r,t)+c2s)∂r+GM(r,t)r2=0, (13)

while the energy equation is

 ∂vT∂t+ur∂vT∂r+(1+ηvTur)vTurr=0. (14)

In these equations is the mean radial velocity of the fluid, excluding the turbulent velocity, which is denoted by . In the rest of this section only, we will neglect the sound speed to simplify the presentation. In any case, we find that its inclusion has little effect on the results.

The mass where

 Mg(r,t)=4π∫r0+r2ρ(r,t)dr (15)

is the gas mass.

We will employ rather than , with the understanding that may refer to the mass of a single star or to a star cluster. For later use we define the radius or sphere of influence of the nascent star or star cluster by

 M∗(t)=4π∫r∗(t)0+r2ρ(r,t)dr. (16)

Thus is the radius at which the enclosed gas mass equals the mass in stars; this makes sense if the distribution of young stars is more centrally concentrated than the gas, which is true on large enough scales. For a pure power law density distribution ,

 r∗(t)=R[M∗(t)Mcl]1/(3−kρ), (17)

where is a fiducial radius, and is the mass of gas inside that radius. We have in mind the thickness of the ribbon-like features found in hydrodynamic simulations, e.g., those in paper I, for .

### 2.2. Solution

We will concentrate on a solution from the time when a star first forms. We define scaled variables

 x≡rr0 (18)

and

 y≡tτDYN, (19)

where is a fiducial radius, and .

We define

 p≡rvT∂vT∂r, p′≡rur∂ur∂r, and kρ≡−rρ∂ρ∂r, (20)

anticipating that each will vary only slowly with radius. In other words, we look for solutions in which the density, infall velocity, and turbulent velocity are all approximately power laws in both and ,

 ρ(r,t)=ρ0x−kρyα, (21)
 ur(r,t)=u0xp′yβ, (22)

and

 vT(r,t)=v0xpyγ. (23)

where we have defined , with similar expressions for and .

For an infall solution, we expect that the magnitude of the infall velocity increases with time; mass accumulates in the central object, which then exerts a gravitational force on the surrounding material that increases with time. Our energy closure then ensures that the turbulent velocity will increase (at a fixed radius) with increasing time, so that

 β>0,γ>0. (24)

We define , so that .

### 2.3. Solution at Small Radius, r<r∗(t)

The continuity equation becomes

 ∂~ρ∂tx−kρ+u0ρ0r0xp′−kρ−1yα+β[2+p′−kρ]=0. (25)

As , the second term grows more rapidly than the first term, since ; in fact we will show that . Thus at small radii ( satisfying this equation demands that both terms go to zero independently. For the larger, second term, this implies

 kρ=2+p′. (26)

Similarly, for the first term, we find

 ∂ρ(r,t)∂t=0. (27)

This is a very striking result: self-similar gravitational collapse (in which the density and infall velocity follow power laws in radius) evolves to an attractor, in which the density does not vary with time. Note that the argument has made no reference to the type of pressure support, i.e., thermal, turbulent, or magnetic.

#### 2.3.1 Using the energy closure to relate ur to vT

Using the power law ansatz, Equation (14) becomes

 γ v0τDYNyγ−1xp+u0v0r0yβ+γxp+p′−1× (28) [p+1+(vT(r,t)ur(r,t))η]=0.

We can then show

 vT(r,t)=−ur(r,t)1+pη−γηrt. (29)

At small and large , the result (29) implies that , so that (29) is the analog of the result in Robertson & Goldreich (2012) that , where is the scale factor in their cosmological simulations and is given by their Equation (9). We note that this result implies

 p=p′,β=γ,for r

#### 2.3.2 The Momentum Equation at Small Radii

We start by noting that

 M(r,t)=M∗(t)+4π3−kρr30ρ0x3−kρ, (31)

where the second term on the right hand side does not depend on time, a result that follows from equation (27).

The momentum equation is

 (u0τDYNyβ−1xp′+4πGρ0r03−kρx1−kρ)+ [u20r0p′y2βx2p′−1+v20r0(2p−kρ)y2γx2p−1+ GM∗(t)r20x−2]=0. (32)

We have written equation (2.3.2) so as to group together terms that scale with in the same manner. From equations (25) and (30), ; since we expect gravity to cause the infall to increase the magnitude of as decreases, . It follows that , i.e., the pressure gradient term is larger in magnitude than the gas self-gravity term, and also larger than . The pressure gradient term and the advective term must together balance the stellar gravitational term, so we group those three terms together in the square brackets. On setting the sum of the three terms to zero, we have

 ur=−Γ√GM∗(t)r, (33)

where we have used Equation (29), and defined

 Γ≡[(kρ−2p)(1+pη)2−p′]−1/2. (34)

Combining Equation (33) with Equation (26) shows that

 p′=p=−12,kρ=3/2,(r

Equation (33) also shows that .

The remaining two terms, the time rate of change of and the gas self-gravity, inside the parentheses in equation (2.3.2) must separately cancel, or

 ∂ur∂t=−4πGρ(r0)r03−kρx1−kρ. (36)

This shows that the infall velocity increases linearly with time () after the central object forms. Integrating,

 ur(r,t)=ur(r,0)−4πGρ(r0)r03−kρx−1/2t. (37)

The two expression (33) and (37) must agree at late times, so that:

 M∗(t)=Γ−2(4πρ(r0)3−kρ)2Gr30t2. (38)

We show in the next section that for . This allows us to connect the density at any to the density at ,

 ρ(r) = ρ(R)(rr∗)−3/2(r∗R)−kρ (39) = ρ(R)(rR)−3/2⋅(r∗R)3/2−kρ

The fact that suggests that we define

 ψ(kρ,r∗(t))≡(r∗R)3/2−kρ. (40)

Because , the dependence of on is very weak; if , .

Using this, the stellar mass at time can be written in terms of the clump mass as

 M∗(t)=ψ⋅(π2Γ)2Mcl(t−t∗τff)2, (41)

where we now explicitly display , the time at which the central density diverges, or that at which the central star forms.

It follows that

 r∗(t)=ψ2/3⋅[π2Γ]4/3(t−t∗τff)4/3R, (42)

where to a first approximation it is acceptable to treat as a constant, independent of both and time.

It is instructive to compare the acceleration terms in the momentum equation, normalizing to . The acceleration due to the pressure gradient is

 1g(1ρ∂P∂r)≈(2p−kρ)(vTur)2Γ2≈−55+4η2, (43)

while the net acceleration of the gas is approximately

 1gdurdt≈p′Γ2≈−4η25+4η2. (44)

For , the ratio between the acceleration due to the pressure gradient and gravity is (equation [43]), i.e., they nearly balance each other. Hence, the net acceleration of the gas is only roughly one quarter that of gravity. In other words, for , gravity dominates the dynamics, but it is nearly balanced by turbulent pressure, and the infall velocity is substantially smaller than the free fall velocity.

In the next section we show that for , both and are greater than zero; the turbulent pressure decelerates the infall, and the turbulent velocity tracks the infall velocity, as prescribed by Equation (11). As the flow crosses , the acceleration due to gravity, which we will show scales as for , transitions to scaling as for . The infall velocity responds to this increased gravity rapidly, i.e. in a moderately small fraction of the local dynamical time; the turbulent velocity follows, but with a substantial lag (of order the local dynamical time), eventually reaching , as enforced by our energy closure, Equation(11). The acceleration due to the pressure gradient tracks that due to gravity at all radii. However for the pressure term dominates the gravity term, while for , the gravity term is roughly twice the pressure term, as we have just seen.

#### 2.3.3 Summary at r<<R

Gathering the results, we have shown that the turbulent velocity tracks the infall velocity, that (from the continuity equation), and for ; it follows that for .

The mass accretion rate for is roughly independent of radius,

 dMg(r,t)dt = 4πr2ρ(r,t)ur(r,t) (45) ≈ 4πr20ρf(r0)ur(r0,t),

where we used the fact that for .

Although the density , i.e., it is independent of time, the accretion rate does vary with time since the infall velocity does. In other words, the accretion rate does not vary with radius (for ) at a fixed time, but it does increase linearly with time at a fixed radius.

### 2.4. Solution for r∗<r≲R

It is useful to define times scales for the density, radial velocity, and turbulent velocity to change by a factor of order unity:

 τρ≡(1ρ∂ρ∂t)−1, (46)
 τur≡(1ur∂ur∂t)−1, (47)

and

 τvT≡(1vT∂vT∂t)−1. (48)

Well away from and , we expect that these should be of order the local dynamical time

 τdyn≡rvr. (49)

The momentum equation can be written as

 τdynτur+p′+(vTur)2(2p−kρ)+(vTur)2GM/rv2T=0. (50)

Since the turbulent velocity tracks the infall velocity, varies slowly with ; for , we expect to vary slowly with as well. Then

 vT(r,t)=Γ′√GM(r,t)r, (51)

where

 Γ′≡∣∣∣vTur∣∣∣Γ=[kρ−2p−(p′+τdynτur)(urvT)2]−1/2 (52)

is a slowly varying function of for . Note that , the virial parameter.

Because is a slowly varying function of , it follows from Equation (51) that

 2p=2−kρ,r>r∗. (53)

The continuity equation can be written as

 τdynτρ+2+p′−kρ=0. (54)

Combining this with Equation (53), we find

 p′=∣∣∣τdynτρ∣∣∣−2p,r>r∗. (55)

In the appendix we argue that if the outer boundary condition on the density is , we can approximate

 ∣∣∣τdynτρ∣∣∣≈∣∣∣lnr/r∗lnR/r∗−1∣∣∣. (56)

In the Appendix we show that

 p(r)=12kρ−12(urvT)2[τdynτur+p′+Γ−2(0)], (57)

and estimate .

From Equation (53), we find

 kρ≈1.6−1.8, (58)

again depending on the initial conditions and on .

Unlike the case at , here the dynamics are controlled by the pressure gradient, and not by gravity. The ratio

 1g(1ρ∂P∂r)=−(kρ−2p)Γ′2≈−3to−4 (59)

for , the fiducial value used in the numerical simulations below, and , motivated by the results of the Appendix.

The force due to the steep pressure gradient (or more precisely the density gradient, since ) is larger than the force of gravity. This outward-directed force causes the infall velocity to decrease inward, for . This is quantified by combining Equation (55) with Equation (56). At , we expect , so that decreases with decreasing radius.

This result, that decreases inward for , is in contrast to self-similar solutions with an assumed fixed turbulent velocity. In such theories, the infall velocity is either zero outside the radius of the expanding collapse wave (Shu, 1977; McKee & Tan, 2003) or has a magnitude that increases inward even at large radii (Fatuzzo et al., 2004).

In the next section, we use direct numerical integration of the 1D equations to show that , depending on the outer boundary conditions and on , similar to the result of 3D simulations in Paper I.

This value of the turbulent exponent is much smaller than the exponent in non-self-gravitating supersonic turbulence, , and smaller than the Kolmogorov value for subsonic turbulence. Recall that the Kolmogorov value follows from conservation of energy in the cascade from large scales to small scales. In the case of supersonic collapse, the exponent is smaller than the value one would derive under the assumption that the kinetic energy in the cascade is conserved. Thus, the kinetic energy on small scales is larger than that produced by an energy conserving cascade.

The reason is simple: a fraction of the ordered inflow kinetic energy, and a fraction of the potential energy released in the collapse, is being converted into turbulent kinetic energy. In our simple theory, this fraction is set by the closure relation, Equation (11), which captures the behavior seen in the numerical simulations of Robertson & Goldreich (2012). This extra energy input boosts the turbulent velocity relative to that in an energy conserving cascade.

Inside the sphere of influence the conversion of gravitational potential energy into turbulent energy is even more dramatic, as we saw above (), so that the turbulent velocity increases with decreasing .

Simulations of non-self-gravitating turbulence consistently show , as do observations of GMCs (Larson, 1981; Heyer & Brunt, 2004) and low mass star formation regions (Myers, 1983; Falgarone et al., 1992) in the Milky Way. We suggest that measurements of the scaling of turbulent velocity can be used as a sign post for collapse, a point we return to below.

As in the case , we have found that the dynamics converges to a run of density that varies only on the large scale dynamical time. In this case (for ) is slightly larger than the value found for .

Unlike the small case, however, we also find that the run of velocity (either or ) evolves rapidly to a solution that varies only on the dynamical time evaluated at the outer scale , i.e., for the velocity increases only very slowly with time. This variation, on the global dynamical time, is simply a result of the change in the mass contained inside . For the turbulent velocity is well described by a power law (with ) while the infall velocity shows substantial deviations from a simple power law, as given by eqns. (55) and (56).

It follows that the mass accretion rate for varies fairly rapidly with radius,

 dMg(r,t)dt=4πr2ρ(R)ur(R,t)(rR)2−kρ+p′, (60)

where . We have allowed for the possibility that the infall velocity at increases with time, as the growing mass of the clump will tend to accelerate the infall, but this increase is on the (large scale) dynamical time, and hence proceeds rather leisurely; essentially, the mass accretion rate is independent of time for , but not for !

In a more realistic model, including the effects of thermal gas pressure, magnetic pressure, and magnetic tension, as well as three dimensional effects, this mass accretion rate will be reduced. We discuss these effects briefly in the next subsection.

### 2.5. The Effects of Thermal and Magnetic Pressure, Rotation, and Magnetic Tension

There are several factors that will reduce the mass accretion rate in a more realistic model. First, we do not expect that gas will fall onto the star forming region from every direction, as we have assumed so far. In the three dimensional simulations in Paper I, we find that the infall covers a solid angle steradians, which we account for by a factor . In our MHD simulations, the fraction of the sky (as seen from the central object) over which accretion occurs is yet smaller, as the magnetic field line tension inhibits accretion across field lines, described by a factor . Finally, the pressure associated with the magnetic field energy density provides an outward pressure gradient much like the turbulent pressure.

Rotational support will tend to slow the rate of accretion. However, observations suggest that the rotational velocity in star forming cores is much smaller than the turbulent velocity, e.g., Beltrán et al. (2011); Palau et al. (2014), who find at . This appears to hold down to scales of order or , below which there is evidence for rotating disks, e.g., van der Tak et al. (2000).

We can include the effects of thermal pressure, and, in an approximate way, some of the effects of magnetic pressure and tension. All three tend to reduce the mass accretion rate, although for the massive star formation regions we consider here, with large virial velocities and high Mach numbers, the effects of thermal pressure are small.

The effects of magnetic pressure are potentially very significant. We expect that the strong turbulence generated in the converging flows we consider will tend to magnify the ambient magnetic field, possibly reaching equipartition, i.e., a magnetic pressure equal to , or a substantial fraction thereof. This effect can be incorporated in a very rough way by assuming a pressure of the form

 P(r)=ρ(r)[c2s+ϕBv2T], (61)

where . This expression does not entirely capture the effect of the magnetic pressure, since it implies that the magnetic field decays when the turbulent velocity does, contrary to what we expect. This might be important if the magnetic pressure becomes large enough to reduce the infall velocity of the flow; reducing will reduce , and hence in (61), the pressure associated with . In fact, once the turbulence generates the magnetic field, the decay time of the magnetic field need not be related to the decay time of the turbulence.

We denote the fraction by which the magnetic pressure decreases the accretion rate by ; we show below, using our simple model, that for a plasma , .

Accounting for all these effects, the expression for the stellar mass becomes

 M∗(t) = fΩfΩBfB2ψ(π2Γ)2Mcl(t−t∗τff)2 (62) = ϕMcl(t−t∗τff)2. (63)

## 3. Numerical results

We numerically integrate the full set of equations (12), (13), and (14), supplemented by (15) using a simple first order upwind code.

Our initial conditions are that the density is constant, , while the radial and turbulent velocities follow a power law,

 ur(r,t=0)=−|ur(R)|(rR)p′ (64)

and

 vT(r,t=0)=vT(R)(rR)p, (65)

with fiducial values , , and . The default value for the eddy decay time-scale is . We have explored variations around these fiducial values, and report on some of the results below. The boundary radius for most of our runs. The initial mass is , and the corresponding Keplerian velocity is

 vk(R,0)=√GM(R,0)/R≈2.8kms−1. (66)

Our outer boundary conditions specify the turbulent velocity and density, and . The inflow velocity is set to be a power law extrapolation of , where is the grid scale. The inner boundary conditions are , , and , i.e., we insist that the mass flowing into the innermost shell , which is deposited onto the central star, be the same as the mass flow rate into the surrounding shell.

We have integrated the equations both with and without including the effects of a finite sound speed. The results are virtually identical to those when the sound speed is neglected.

Figure 1 shows the run of density and infall velocity for fifteen snapshots between and , where . The first four snapshots show a core in the density distribution, , with the constant increasing rapidly from in the initial snapshot to in the sixth snapshot. By the eighth snapshot, the density distribution has evolved into a (slightly broken) power law. The break in the power law occurs at ; for , , while for , , with .

The infall velocity also evolves into a broken power law. For , the exponent , while for , , before flattening at .

Figure 2 shows the ratio

 (∇ρv2T)/ρGM(r,t)/r2 (67)

of the pressure gradient term to the gravity term in the momentum equation. We distinguish between times previous to the formation of the central object (shown by dotted lines) and times after the central object forms (shown by solid lines). These curves correspond to the same times as the curves in the previous figure.

At early times and small radii, the density is nearly constant, so , and the pressure gradient provides an inward force, i.e,., it acts in the same sense as gravity. Thus the pressure gradient force and the force due to gravity both act to move gas toward the origin. However, at large radii, the dotted curves drop down, eventually becoming negative at large radii, i.e., the pressure gradient term is directed outward, opposing the force due to gravity. Near the outer boundary, this outward pressure gradient term is larger in magnitude than the acceleration of gravity.

After the first time step, the radius at which the pressure gradient changes sign is located at ; well inside this radius, the density is uniform, At each subsequent plotted time, the transition radius decreases by a factor of about three. Once again, at a fixed time, but well inside this transition radius, the density is uniform.

After the central object forms the situation is very different. To begin with, the pressure gradient is directed outward at all radii. However while the pressure gradient force is directed outward, it is smaller in magnitude than the force due to gravity at radii . Thus the infall velocity of a (Lagrangian) mass shell increases inward. The ratio of the two forces given by Equation (43) is constant, and is depicted in Figure 2 by the horizontal dashed line. The numerical result at late times (after the star forms) tracks this prediction for . As a result, the infall velocity , but it is smaller than the free-fall velocity by a factor . This is different than the prediction of the turbulent core model, in which the pressure term becomes negligible as , so that the infall velocity approaches the free-fall value.

In contrast, for the pressure gradient term is larger in magnitude than the gravity term, so that the infall velocity of a mass shell decreases inward. Note that the inertial term is much smaller than either the pressure gradient term or the gravity term, or their difference, so that the infall velocity decreases in an Eulerian sense as well, as seen in Figure 1.

Figure (3) shows and . Note that and are nearly constant for , while is roughly constant for . The prediction of Equation (34) using the values for the exponents at is shown in the figure as the horizontal dotted line; it agrees well with the numerical results.

Figure (4) shows the exponent in the size-linewidth relation, . The exponent is positive for , but much smaller than the value seen in simulations of non-self-gravitating turbulence. For the value of quickly plunges to , the value of a Keplerian profile. However, the normalization is below that of the Keplerian free-fall velocity as we discussed above. The predictions of the analytic theory, equations (A9) for and (A11) for are also shown, and agree well with the numerical results except in the immediate vicinity of .

Figure (5) shows the result for ; the prediction for and the prediction of Equation (55) are shown as the dashed and dotted lines, respectively. Here we note again that the gravity from the central star drives the radial velocity to a Keplerian profile inside ; while we do not show it here, the magnitude of the infall velocity does not reach free-fall values.

Figure (6) plots the density exponent (solid line) and the prediction of Equation (35) for (dashed line), and of Equation (53) for (dotted line). The prediction of is recovered nicely in the numerical results.

Figure (7) shows the accreted mass versus , with defined to be the time when the accreted mass exceeds some minimum mass ; in the figure, . Varying alters the appearance of the plot at early times but has no effect at later times. The mass increases in proportion to , and the numerical result is in good agreement with the prediction of Equation (41).

We have calculated the run of density and velocity when the effects of a strong magnetic field, corresponding to a plasma relative to