THE SPHERICALLY SYMMETRIC GRAVITATIONAL COLLAPSE OF A CLUMP OF SOLIDS IN A GAS

A layer of particles can settle toward the disk midplane if this region of the disk is sufficiently free of turbulent velocity fluctuations, as may happen (Gammie 1996) when there is insufficient ionization to sustain magneto-rotational instability. A rich set of mechanisms operate in this layer, which we now discuss.

Seeking an alternative to coagulation, Safronov (1969) and Goldreich & Ward (1973) showed that in the absence of gas, the midplane layer undergoes (axisymmetric) gravitational instability at a sufficiently short radial wavelength. Sekiya (1983) included the presence of gas assuming it is "perfectly coupled to particles," meaning that the particle velocity ${\boldsymbol{U}}$ exactly equals the gas velocity ${\boldsymbol{u}}$. This becomes true (modulo a gravitationally induced sedimentation velocity) as particles become very small. Sekiya showed that in this limit, gas pressure imposes a very high critical value (e.g., $6\times {{10}^{6}}$ at 5 AU) for the ratio (Φ ≡ ${{\rho }_{{\rm p}}}/{{\rho }_{{\rm g}}}$) of particle to gas density required for gravitational instability. This is easily understood, as described near Equation (32) below, by replacing the sound speed with the effective sound speed in the coupled two-fluid medium. This replacement was suggested by Safronov (1987)¹ and rediscovered by Cuzzi et al. (2008). Sekiya finds that another mode exists whose critical particle density assumes a Roche-like value of $\rho _{{\rm p}}^{*}=0.6\;{{M}_{\odot }}/{{\varpi }^{3}}$, where $\varpi$ is the orbital radius. This density leads to a much lower critical value of Φ = 170 at $\varpi =5$ AU. It is important to note that this mode is incompressible, i.e., it cannot increase the particle density. While such a mode can lead to fragmentation of the midplane layer, it can contribute to planetesimal formation only via (slow) sedimentation toward the center of the fragments. How the constraint imposed by gas pressure relaxes with increasing Stokes number has yet to be quantified for this situation, to our knowledge.

Even if one accepts that the incompressible mode is relevant to planetesimal formation, the midplane layer cannot easily achieve the required density due to self-generated turbulence. Particles at the midplane revolve faster than the gas, which moves slower than the Keplerian speed because it has slight pressure support. This results in an Ekman-like layer with vertical shear, which makes it susceptible to turbulence driven by the Kelvin–Helmholtz instability, which could prevent the critical density from being achieved (Weidenschilling 1980). On the other hand, the layer can be stably stratified due to the presence of particles. The ratio of stabilizing buoyancy to shear is called the Richardson number, ${\rm Ri}$. Ignoring differential rotation together with the Coriolis force, and assuming that particles are perfectly coupled to the gas, one has classical stratified shear flow which is known to be stable provided that the local ${\rm Ri}\geqslant 1/4$ everywhere (Drazin & Reid 1981). Garaud & Lin (2004) find that as the particle layer settles, this (sufficient) criterion for flow stability is violated before Sekiya's incompressible gravitational mode becomes unstable. Cuzzi et al. (1993) solved the Reynolds-averaged two-fluid equations for a turbulent midplane layer using mixing-length-type models for various covariances in the more general case of finite stopping time. They conclude that Sekiya's critical particle density is approached at 10 AU only for particles of about 1 m in size. Likewise, the critical density is not achieved in the simulations of Johansen et al. (2006), which use Lagrangian particles and ignore Keplerian shear. The latter can reduce the strength of the Kelvin–Helmholtz vortices by tilting them away from the vertical shear; however, it can also transiently strengthen them by stretching. Simulations by Barranco (2009) included the Coriolis force and radial shear in the single-fluid limit (i.e., assuming that the particle velocity exactly equals the gas velocity) and concluded that the flow is turbulent. However, these simulations could not obtain the equilibrium particle density distribution because, in the perfectly coupled limit, the particle mass fraction is an advected scalar and retains its maximum value from the initial condition. This is where work on midplane Kelvin–Helmholtz turbulence currently stands.

An interesting mechanism known as the streaming instability was discovered by Youdin & Goodman (2005). It creates particle concentrations in the midplane layer and under certain circumstances can counteract particle diffusion due to the Kelvin–Helmholtz instability. The basic state is the Nakagawa et al. (1986) solution for radial and azimuthal particle drift relative to the gas flow. Youdin & Goodman (2005) perform a local axisymmetric linear stability analysis for this basic state assuming uniform gas and particle density. In particular, vertical shear leading to Kelvin–Helmholtz instability is absent. Non-linear simulations with the same set-up are presented in Johansen & Youdin (2007). The streaming instability does not depend on self-gravity and simultaneously clumps particles both radially (x) and vertically (z), i.e., unstable modes have ${{k}_{x}}\ne 0$ and ${{k}_{z}}\ne 0$; see Figure 2 in Youdin & Goodman (2005). Clumping in the radial direction is explained by Jacquet et al. (2011) as follows: imagine a perturbation mode in which particles are radially compressed. If the stopping time is sufficiently short, the particles will also compress the gas via drag, leading to a local pressure maximum. The pressure maximum then acts as a further attractor for particles (as is well known and described below), which closes the feedback loop. Note that the analysis of Youdin & Goodman (2005) takes the gas to be incompressible; in this limit the increase in gas pressure arises from enforcement of the incompressibility constraint. For the mechanism of radial clumping put forward by Jacquet et al. (2011), one expects optimal instability growth at an intermediate particle size small enough for particles to drag the gas but large enough to undergo radial drift. Finally, we mention two computational studies (Johansen et al. 2009; Bai & Stone 2010) whose set-up allows for both the Kelvin–Helmholtz and streaming instabilities to occur. These studies also include multiple particle sizes in each simulation and therefore have the additional effect of differential radial drift. They conclude that significant clumping takes place when the (vertically averaged) density ratio of dust to gas is supersolar and particles sufficiently large.

Since the Kelvin–Helmholtz and streaming instabilities work against each other in concentrating particles, one may ask under what conditions is the latter sufficiently strong that gravitational instability can take place. In a brief computational study, Johansen et al. (2009) introduced particles with four radii, namely 3, 6, 9, and 12 cm at 5 AU. They found that when nebular dust to gas ratio was Z = 0.01, clumping was insignificant. However at a super-solar value of Z = 0.02, clumping sharply increased.

A local maximum of gas pressure at orbital radius $\varpi ={{\varpi }_{{\rm max} }}$ attracts solids toward it (Haghighipour & Boss 2003). This is because for $\varpi \lt {{\varpi }_{{\rm max} }}$, we have $\partial p/\partial \varpi \gt 0$ leading to faster than Keplerian gas velocity. This boosts the orbital radius of particles, while the opposite is true for particles at $\varpi \lt {{\varpi }_{{\rm max} }}$. A more quantitative rendition of this is as follows. The analysis of Nakagawa et al. (1986) is easily extended to non-small stopping times to give for the radial drift velocity of particles:

$$\begin{eqnarray}&&{{U}_{\varpi }}=-2\eta {{u}_{{\rm K}}}\left[ \frac{{\rm S}{{{\rm t}}_{{\rm K}}}}{(1+{\Phi })+{\rm St}_{{\rm K}}^{2}} \right],\end{eqnarray} \tag{ 1 }$$

where ${{u}_{{\rm K}}}$ is the Keplerian velocity, Φ = ${{\rho }_{{\rm p}}}/{{\rho }_{{\rm g}}}$ is the particle loading, ${\rm S}{{{\rm t}}_{{\rm K}}}={{t}_{{\rm s}}}{{u}_{{\rm K}}}/\varpi$ is the Stokes number based on the characteristic orbital time, and η is the pressure gradient parameter

$$\begin{eqnarray}&&\eta =-\frac{\varpi }{2{{\rho }_{{\rm g}}}u_{{\rm K}}^{2}}\frac{\partial p}{\partial \varpi }.\end{eqnarray} \tag{ 2 }$$

From (1), we see that the radial drift velocity is maximized for ${\rm S}{{{\rm t}}_{{\rm K}}}={{(1+{\Phi })}^{1/2}}\gt 1$. At 3 AU in a minimum-mass nebula this implies particles with radii $a\gtrsim 70$ cm. In other words, this mechanism is most effective for meter-sized bodies. Related to the above phenomenon is the trapping of particles by large-scale anti-cyclonic vortices (Barge & Sommeria 1995), which are pressure highs, or in the vortices created by magneto-rotational turbulence (Johansen et al. 2007). Meter-sized bodies are also the ones that drift toward the central star most rapidly (which is, of course, not coincidental) and Johansen et al. suggest that such bodies can avoid this fate by being trapped in vortices where they can gravitationally collapse on dynamical timescales. This suggestion neglects fragmentation of the boulders in the eddies by collision.

The difficulty with the above mechanisms is that they require either large particles or a nebula with an enhanced solids-to-gas mass fraction. Another scenario for planetesimal formation (Cuzzi et al. 2008, 2010) envisions a direct path from chondrule (mm) sized particles to asteroidal sized bodies. It invokes the phenomenon of turbulent concentration, first articulated by Maxey (1987), whereby particles centrifuge out from regions of high vorticity and accumulate in regions of high strain. Simulations of homogeneous isotropic turbulence confirmed this (Squires & Eaton 1991) and found that the effect was most pronounced for particles having a certain ratio of particle stopping time. Subsequently, Wang & Maxey (1993) refined this result by finding that particles having Stokes number ${\rm S}{{{\rm t}}_{\eta }}\approx {{t}_{{\rm s}}}/{{t}_{\eta }}\approx 1$ are the ones most concentrated. Here ${{t}_{{\rm s}}}$ is the particle stopping or response time and ${{t}_{\eta }}$ is the Kolmogorov time. It so happens that chondrule-sized particles approximately fulfill this condition (Cuzzi et al. 2001). The larger and slower eddies concentrate larger particles (Bec et al. 2007; Pan et al. 2011). The current status of this effort will be provided in the concluding remarks.

Many of the above mechanisms produce dense particle clumps. They include the Kelvin–Helmholtz instability (Johansen et al. 2006), preferential concentration in turbulence (e.g., Cuzzi et al. 2001, 2008, 2010; Pan et al. 2011), and the streaming instability (Johansen & Youdin 2007; Bai & Stone 2010). The present work is concerned with the gravitational collapse of such clumps. Our collapse solutions cover the relevant regimes of the three non-dimensional parameters, namely, the Stokes number, particle loading, and the two-phase Jeans number. In reality, a particle clump will also be subject to dispersive effects from within and without. For the case of preferential concentration by turbulence, the former includes turbulent gas motions within the clump that lead to a particle dispersion velocity. The latter includes straining by eddies of the same size and larger than the clump, and ram gas pressure because dense clumps revolve faster than the gas (Cuzzi et al. 2008). While these dispersal mechanisms should be studied in more detail in the future, in Section 5 we provide an estimate for the rate of erosion by turbulent strain.

The equation of motion of a single particle of mass ${{m}_{{\rm p}}}$ and radius a subject to gas drag is

$$\begin{eqnarray}&&{{m}_{{\rm p}}}\frac{d{\boldsymbol{U}} }{dt}=\frac{4\pi {{a}^{2}}}{3}\;\bar{c}{{\rho }_{{\rm g}}}({\boldsymbol{u}} -{\boldsymbol{U}} ),\end{eqnarray} \tag{ 3 }$$

where ${\boldsymbol{U}}$ is the particle velocity, ${\boldsymbol{u}}$ is the gas velocity, and ${{\rho }_{{\rm g}}}$ is the gas density. The quantity $\bar{c}$ is the mean thermal velocity of gas molecules defined by Epstein (1924) as

$$\begin{eqnarray}&&\bar{c}={{(8/\pi )}^{1/2}}\;{{c}_{{\rm i}}},\end{eqnarray} \tag{ 4 }$$

where ${{c}_{{\rm i}}}$ is the isothermal sound speed. Equation (3) uses the Epstein (1924) drag law, which is valid when $|{\boldsymbol{u}} -{\boldsymbol{U}} |\ll {{c}_{{\rm i}}}$ and $a\lt {{\lambda }_{{\rm g}}}$, where ${{\lambda }_{{\rm g}}}$ is the mean free path of the gas. We will employ the two-fluid treatment for solid particles, which applies mass and momentum conservation to a differential volume containing many particles. Let all particles be identical (the so-called mono-disperse case); this is convenient but not essential. For a differential volume, consider a physical volume large enough to contain many particles but small enough that ${\boldsymbol{u}} ({\boldsymbol{x}} )$ does not change appreciably across it. In a protoplanetary disk the distance between particles (∼1 m) is much smaller than the smallest scale across which changes in velocity occur (the Kolmogorov scale ∼1 km). We also require that ${\boldsymbol{U}}$ not change appreciably across the volume and be a single-valued function of particle position. This assumption disallows "crossing trajectories" such as when two particle streams traveling in different directions and/or at different velocities interpenetrate without colliding. In this case, the continuum treatment averages the particle velocity in the interpenetration zone leading to unphysical behavior. One way to treat this case in the future is to solve a kinetic equation for the particle number density in the position–velocity phase space (Chalons et al. 2012). The final assumption is that particles are sufficiently separated from each other that the drag law (Equation (3)), which is for a single particle in an infinite medium, remains valid.

When solid particles are treated as a continuum phase, the following two-fluid equations of motion result:

$$\begin{eqnarray}&&{{\rho }_{{\rm g}}}_{,t}+{{({{\rho }_{{\rm g}}}{{u}_{i}})}_{,i}}=0,\quad ({\rm gas}\ {\rm mass})\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{({{\rho }_{{\rm g}}}{{u}_{i}})}_{,t}}+{{({{\rho }_{{\rm g}}}{{u}_{i}}{{u}_{j}})}_{,j}} & = & -{{p}_{,i}}+{{\rho }_{{\rm g}}}{{g}_{i}}-{{d}_{i}}, \\ {} & {} & \quad ({\rm gas}\ {\rm momentum}) \\ \end{array}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{e}_{,t}}+{{\left[ {{u}_{j}}(p+e) \right]}_{,j}}={{\rho }_{{\rm g}}}{{g}_{i}}{{u}_{i}}-{{d}_{i}}{{u}_{i}}-{{q}_{j,j}}\quad ({\rm gas}\ {\rm energy})\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{{\rho }_{{\rm p}}}_{,t}+{{({{\rho }_{{\rm p}}}{{U}_{i}})}_{,i}}=0,\quad ({\rm particle}\ {\rm mass})\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{({{\rho }_{{\rm p}}}{{U}_{i}})}_{,t}}+{{({{\rho }_{{\rm p}}}{{U}_{i}}{{U}_{j}})}_{,j}} & = & {{\rho }_{{\rm p}}}{{g}_{i}}+{{d}_{i}}. \\ {} & {} & \quad ({\rm particle}\ {\rm momentum}) \\ \end{array}\end{eqnarray} \tag{ 9 }$$

Here ${{\rho }_{{\rm p}}}$ is the density of the particle phase, q_j is the heat flux vector, and

$$\begin{eqnarray}&&{{d}_{i}}=\beta \bar{c}{{\rho }_{{\rm g}}}{{\rho }_{{\rm p}}}\left( {{u}_{i}}-{{U}_{i}} \right),\ {\rm with}\ \beta \equiv \frac{4\pi {{a}^{2}}}{3{{m}_{{\rm p}}}},\end{eqnarray} \tag{ 10 }$$

is the drag force per unit volume exerted on the particle phase by the fluid phase. The quantity ${{g}_{i}}={{\varphi }_{,i}}$ is the gravitational acceleration, where the potential $\varphi$ satisfies Poisson's equation

$$\begin{eqnarray}&&{{\varphi }_{,ii}}=-4\pi G({{\rho }_{{\rm g}}}+{{\rho }_{{\rm p}}}).\end{eqnarray} \tag{ 11 }$$

The equation of state for an ideal gas of specific heat ratio γ is:

$$\begin{eqnarray}&&e=\frac{p}{\gamma -1}+\frac{1}{2}{{\rho }_{{\rm g}}}|{\boldsymbol{u}} {{|}^{2}}.\end{eqnarray} \tag{ 12 }$$

We have implemented the adiabatic and isothermal cases but have performed calculations for only the latter, in which case the equations of energy (Equation (7)) and state (Equation (12)) are replaced by

$$\begin{eqnarray}&&p={{\rho }_{{\rm g}}}{{c}_{{\rm i}}}^{2}.\end{eqnarray} \tag{ 13 }$$

For the adiabatic case, the heat flux vector is set to zero (q_i = 0).

While the competition between gas pressure and gravity is readily appreciated in the star formation context, it is less obvious in the context of gravitational collapse of solid particles. In this sub-section we introduce a two-phase acoustic speed ${{c}_{{\rm eff}}}$, which was first appreciated by Sekiya (1983); see also the review article by Marble (1970) on dusty gases. This leads to a two-phase Jeans parameter, ${{J}_{{\rm t}}},$ which is the classical Jeans parameter with the acoustic speed replaced by ${{c}_{{\rm eff}}}$.

Let the particle response time ${{t}_{{\rm s}}}$ be much shorter than the flow timescale ${{t}_{{\rm f}}}$ over which the velocity field changes following a gas particle; ${{t}_{{\rm f}}}$ can be obtained from the velocity gradient tensor. Then, particles and gas approximately follow each other and we may adopt the Maxey (1987) expansion

$$\begin{eqnarray}&&{{U}_{i}}={{u}_{i}}+{{t}_{{\rm s}}}{{F}_{i}}({\boldsymbol{x}} ,t)+\mathcal{O}\left( t_{{\rm s}}^{2} \right),\end{eqnarray} \tag{ 14 }$$

where ${{F}_{i}}({\boldsymbol{x}} ,t)$ remains to be determined. In writing Equation (14) we have assumed that each quantity has been non-dimensionalized, in particular, that time has been non-dimensionalized using ${{t}_{{\rm f}}}$. Substituting Equation (14) into the particle momentum equation in acceleration form,

$$\begin{eqnarray}&&\frac{\partial {{U}_{i}}}{\partial t}+{{U}_{j}}\frac{\partial {{U}_{i}}}{\partial {{x}_{j}}}={{g}_{i}}+({{u}_{i}}-{{U}_{i}})/{{t}_{{\rm s}}},\end{eqnarray} \tag{ 15 }$$

gives

$$\begin{eqnarray}&&{{U}_{i}}={{u}_{i}}+{{t}_{{\rm s}}}\left( {{g}_{i}}-\frac{D{{u}_{i}}}{Dt} \right)+\mathcal{O}\left( t_{{\rm s}}^{2} \right).\end{eqnarray} \tag{ 16 }$$

The second term on the right of Equation (16) is simply a sedimentation velocity generalized to include a D'Alembert term due to gas acceleration. Here we need only the leading order result, U_i = u_i, which allows one to add the momentum equations for the gas and particles to give (for the isothermal case)

$$\begin{eqnarray}&&{{\rho }_{{\rm t}}}_{,t}+{{({{\rho }_{{\rm t}}}{{U}_{i}})}_{,i}}=0,\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{{({{\rho }_{{\rm t}}}{{U}_{i}})}_{,t}}+{{({{\rho }_{{\rm t}}}{{U}_{i}}{{U}_{j}})}_{,j}}=-c_{{\rm i}}^{2}{{({{\phi }_{{\rm g}}}{{\rho }_{{\rm t}}})}_{,i}}+{{\rho }_{{\rm t}}}{{g}_{i}}\end{eqnarray} \tag{ 18 }$$

where ${{\rho }_{{\rm t}}}\equiv {{\rho }_{{\rm g}}}+{{\rho }_{{\rm p}}}$ is the total density and ${{\phi }_{{\rm g}}}\equiv {{\rho }_{{\rm g}}}/{{\rho }_{{\rm t}}}$ is the gas mass fraction. Note that the particle loading Φ ≡ ${{\rho }_{{\rm p}}}/{{\rho }_{{\rm g}}}=\phi _{{\rm g}}^{-1}-1$ and conversely ${{\phi }_{{\rm g}}}={{(1+{\Phi })}^{-1}}$. One may then define a two-phase Jeans ${{J}_{{\rm t}}}$ parameter as the ratio of gravitational to pressure force terms in Equation (18):

$$\begin{eqnarray}&&{{J}_{{\rm t}}}=\frac{{{\rho }_{{\rm t}}}G{{\ell }^{2}}}{{{\phi }_{{\rm g}}}c_{{\rm i}}^{2}},\end{eqnarray} \tag{ 19 }$$

where ${{\phi }_{{\rm g}}}$ and ℓ are the characteristic gas mass fraction and size of the region. The criterion for collapse on the dynamical timescale ${{({{\rho }_{{\rm t}}}G)}^{-1/2}}$ is

$$\begin{eqnarray}&&{{J}_{{\rm t}}}\gt {{J}_{{\rm c}}},\end{eqnarray} \tag{ 20 }$$

where ${{J}_{{\rm c}}}$ is a critical value of $\mathcal{O}(1)$. The ratio of characteristic free-fall velocity ${{R}_{0}}{{({{\rho }_{{\rm t}}}G)}^{1/2}}$ to ${{c}_{{\rm i}}}$ is the dynamical Mach number ${\rm M}{{{\rm a}}_{{\rm dyn}}}$ so

$$\begin{eqnarray}&&{{J}_{{\rm t}}}=(1+{\Phi }){\rm Ma}_{{\rm dyn}}^{2}.\end{eqnarray} \tag{ 21 }$$

The criterion (Equation (20)) then becomes

$$\begin{eqnarray}&&{\rm M}{{{\rm a}}_{{\rm dyn}}}\gt {{\left[ {{(1+{\Phi })}^{-1}}{{J}_{{\rm c}}} \right]}^{1/2}}\approx {{(1+{\Phi })}^{-1/2}}.\end{eqnarray} \tag{ 22 }$$

In other words, the gas must have (sufficient) compressibility to allow dynamical collapse of particles that are tightly coupled to it. Otherwise, the best that one may expect is particle sedimentation through an incompressible gas.

To be complete, Equations (17) and (18) need to be supplemented with a transport equation for ${{\phi }_{{\rm g}}}$. We have

$$\begin{eqnarray}&&\frac{1}{{{\rho }_{{\rm p}}}}\frac{D{{\rho }_{{\rm p}}}}{Dt}=\nabla \cdot {\boldsymbol{U}} ,\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&\frac{1}{{{\rho }_{{\rm g}}}}\frac{D{{\rho }_{{\rm g}}}}{Dt}=\nabla \cdot {\boldsymbol{U}} ,\end{eqnarray} \tag{ 24 }$$

where

$$\begin{eqnarray}&&\frac{D}{Dt}\equiv \frac{\partial }{\partial t}+{\boldsymbol{U}} \cdot \nabla .\end{eqnarray} \tag{ 25 }$$

Subtracting Equation (24) from Equation (23) gives

$$\begin{eqnarray}&&\frac{D}{Dt}({\rm ln} \;{\Phi })=0.\end{eqnarray} \tag{ 26 }$$

Thus ${\rm ln} \;{\Phi }$ is a passive scalar and so is any function of it, in particular ${{\phi }_{{\rm g}}}$:

$$\begin{eqnarray}&&\frac{D{{\phi }_{{\rm g}}}}{Dt}=0\end{eqnarray} \tag{ 27 }$$

This means that in the perfectly coupled limit, the maximum value of the particle loading cannot change from its initial value. If, in addition, ${{\phi }_{{\rm g}}}$ is initially uniform² it remains so and may be removed from under the pressure gradient in Equation (18). The equations of motion then become completely equivalent to those of a compressible isothermal gas with effective speed of sound

$$\begin{eqnarray}&&{{c}_{{\rm eff}}}={{c}_{{\rm i}}}\sqrt{{{\phi }_{{\rm g}}}}.\end{eqnarray} \tag{ 28 }$$

The role of the effective speed ${{c}_{{\rm eff}}}$ seems to have been first recognized by Safronov (1969), Marble (1970), and Sekiya (1983) who accounts for the presence of gas in the Goldreich–Ward midplane instability. The equivalence (Equation (28)) allows many known results for isothermal gases to be directly translated to a perfectly coupled particle–gas system. Here are three examples: (i) Jeans' criterion for gravitational collapse translates to

$$\begin{eqnarray}&&\lambda \gt {{\left( \frac{\pi {{\phi }_{{\rm g}}}c_{{\rm i}}^{2}}{G{{\rho }_{{\rm t}}}} \right)}^{1/2}},\end{eqnarray} \tag{ 29 }$$

for instability where λ is the wavelength of a disturbance. (ii) The dispersion relation for an infinitesimally thin Keplerian gas sheet subject to axisymmetric disturbances is (in the local shearing sheet approximation):

$$\begin{eqnarray}&&{{\omega }^{2}}={{k}^{2}}c_{{\rm i}}^{2}+{\Omega }_{{\rm K}}^{2}-2\pi G{{{\Sigma }}_{{\rm g}}}k.\end{eqnarray} \tag{ 30 }$$

To adapt this to a perfectly coupled particle–gas sheet, one makes the replacements $c_{{\rm i}}^{2}\to c_{{\rm i}}^{2}{{\phi }_{{\rm g}}}$ and ${{{\Sigma }}_{{\rm g}}}\to {{{\Sigma }}_{{\rm p}}}$. A range of unstable wave numbers exists when

$$\begin{eqnarray}&&{{\phi }_{{\rm g}}}\lt {{\left( \frac{\pi G{{{\Sigma }}_{{\rm p}}}}{{{{\Omega }}_{{\rm K}}}{{c}_{{\rm i}}}} \right)}^{2}}.\end{eqnarray} \tag{ 31 }$$

For a minimum mass nebula (Hayashi 1981) we have ${{c}_{{\rm i}}}={{10}^{5}}{{\varpi }^{-1/4}}\ {\rm cm}\ {{{\rm s}}^{-1}}$, ${{{\Sigma }}_{{\rm p}}}=7.1{{\varpi }^{-3/2}}\ {\rm g}\ {\rm c}{{{\rm m}}^{-2}}$, and ${{{\Omega }}_{{\rm K}}}=2\times {{10}^{-7}}{{\varpi }^{-3/2}}\;{{{\rm s}}^{-1}}$, where ϖ is the orbital radius in AU. For $\varpi =1$ AU one obtains

$$\begin{eqnarray}&&{\Phi }\gt 1.8\times {{10}^{8}},\end{eqnarray} \tag{ 32 }$$

which is clearly a very tall order and should be compared with the criterion ${\Phi }\gt 2.5\times {{10}^{8}}$ that Sekiya (his Table IIa) obtains using a more complex analysis. (iii) Shu (1977) presented a solution for the rate of infall during self-similar collapse of a singular isothermal gas sphere. If instead one had a gas–particle sphere with an initial total density distribution ${{\rho }_{{\rm t}}}\sim {{r}^{-1}}$ and uniform gas mass fraction, then with the replacement ${{c}_{{\rm i}}}\to {{c}_{{\rm i}}}\sqrt{{{\phi }_{{\rm g}}}}$, Shu's result becomes

$$\begin{eqnarray}&&\dot{M}=.975c_{{\rm i}}^{3}\phi _{{\rm g}}^{3/2}/G.\end{eqnarray} \tag{ 33 }$$

With the above discussion motivating the use of ${{c}_{{\rm eff}}}$ and ${{J}_{{\rm t}}}$, we now describe numerical calculations of the collapse of a spherical clump of particles without a restriction to small stopping times.

With spherical symmetry the governing equations have the form:

$$\begin{eqnarray}&&\frac{\partial {{\rho }_{{\rm g}}}}{\partial t}+\frac{1}{{{r}^{2}}}\frac{\partial }{\partial r}\left( {{r}^{2}}{{\rho }_{{\rm g}}}u \right)=0,\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \frac{\partial }{\partial t}\left( {{\rho }_{{\rm g}}}u \right)+\frac{1}{{{r}^{2}}}\frac{\partial }{\partial r}\left( {{r}^{2}}{{\rho }_{{\rm g}}}{{u}^{2}} \right) & = & -\frac{\partial p}{\partial r}+{{\rho }_{{\rm g}}}g \\ {} & {} & -\beta \bar{c}{{\rho }_{{\rm p}}}{{\rho }_{{\rm g}}}(u-U), \\ \end{array}\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&\frac{\partial {{\rho }_{{\rm p}}}}{\partial t}+\frac{1}{{{r}^{2}}}\frac{\partial }{\partial r}\left( {{r}^{2}}{{\rho }_{{\rm p}}}U \right)=0,\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&\frac{\partial }{\partial t}\left( {{\rho }_{{\rm p}}}U \right)+\frac{1}{{{r}^{2}}}\frac{\partial }{\partial r}\left( {{r}^{2}}{{\rho }_{{\rm p}}}{{U}^{2}} \right)={{\rho }_{{\rm p}}}g+\beta \bar{c}{{\rho }_{{\rm p}}}{{\rho }_{{\rm g}}}(u-U),\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \frac{\partial e}{\partial t}+\frac{1}{{{r}^{2}}}\frac{\partial }{\partial r}\left[ {{r}^{2}}(e+p)u \right] & = & {{\rho }_{{\rm g}}}gu-u\beta \bar{c}{{\rho }_{{\rm p}}}{{\rho }_{{\rm g}}}(u-U) \\ {} & {} & -\frac{1}{{{r}^{2}}}\frac{\partial }{\partial r}\left( {{r}^{2}}{{q}_{r}} \right). \\ \end{array}\end{eqnarray} \tag{ 38 }$$

Here u and U denote the radial velocity of the gas and particles, respectively, and q_r is the radial radiative heat flux. The gravitational acceleration is

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} g(r) & = & -\frac{GM(r)}{{{r}^{2}}}\ {\rm where}\ M(r) \\ {} & = & \int _{0}^{r}4\pi {{r}^{\prime }}^{2}({{\rho }_{{\rm g}}}({{r}^{\prime }})+{{\rho }_{{\rm p}}}({{r}^{\prime }})-{{\rho }_{{\rm g}0}})d{{r}^{\prime }}. \\ \end{array}\end{eqnarray} \tag{ 39 }$$

Note that the density ${{\rho }_{{\rm g}0}}$ of the background gas has been subtracted out in M(r) to ensure that the background gas by itself remains static. This is equivalent to assuming that there is a background pressure gradient that balances the gravity of the background gas. In the present work, we shall assume that compressive and frictional heating are so slow compared with radiative transfer that the gas remains isothermal.

Spatial derivatives on the left-hand sides of Equations (34)–(38) are computed using a second-order accurate finite-volume scheme. Grid points (where the time-evolved quantities reside) are located at the midpoints of computational shells in the radial direction. In a finite-volume scheme, the rate of change of a shell-averaged quantity is the difference of fluxes at the shell faces divided by the shell volume. At second-order, the distinction between a shell-centered and shell-averaged quantity is immaterial; recall that the midpoint rule for quadrature is second-order accurate. We first compute fluxes per unit area (i.e., without the r² factor) at grid-points, linearly interpolate them to shell faces, multiply by r², take the difference, and divide by the shell volume.

The origin requires a slightly special treatment and two schemes are implemented. In the first, a grid-point is placed at the origin and fluxes per unit area are linearly interpolated on to the spherical surface surrounding it. These fluxes determine the rates of change of densities and energy at the origin. The rates of change of radial gas and particle velocities at the origin are explicitly set to zero (as required by spherical symmetry). All results presented this work are obtained using the first scheme. In the second scheme, the origin is taken to be a flux location rather then a grid point. Fluxes can be interpolated there using the symmetry or anti-symmetry of various quantities across the origin. The only non-zero "flux" at the origin is the pressure since it is both an even function across the origin and not multiplied by an r² factor.

At the outer radial boundary, a choice between two boundary conditions for the gas is provided. The first is a non-reflective boundary condition (Thompson 1987) for the gas-dynamic Euler terms, which sets the rate of change of incoming characteristic variables to zero. All spatial derivatives required to implement this condition are calculated using one-sided differences. The second choice is as follows. When the inflow velocity of the gas at the outer boundary is subsonic, which was always the case for the collapse runs, there is one incoming characteristic (for isothermal flow where the characteristic speeds are $u+{{c}_{i}}$ and $u-{{c}_{i}}$) and one may specify one quantity. We chose to set the gas density equal to its ambient value. All the collapse runs were run with this boundary condition. Some tests were conducted to assess differences when the non-reflective boundary condition was used instead. The differences were very small, indicating that the location of the outer boundary was sufficiently far. For the particle equations, a characteristic decomposition is not possible and the fluxes required on the outer surface of the last shell are obtained by linear extrapolation.

Unlike shock-capturing schemes, the scheme just described is free of numerical dissipation. As a result, there is no mechanism to damp oscillations at the Nyquist wavelength, which equals twice the grid spacing. Such waves are numerically generated, even in a linear problem, whenever there are inhomogeneities such as changes in mesh spacing, boundary treatments, or changes in coordinate Jacobians. These waves have zero phase velocity, and once created, remain. Centrally differenced schemes therefore require a small amount of filtering. At the end of every time step we applied the fourth-order Padé filter given by Equation (C.2.1) in Lele (1992) with his coefficients $\beta =d=0$. It has a transfer function that is near unity and drops off sharply near the Nyquist wavenumber $k{\Delta }r=\pi$ where ${\Delta }r$ is the grid spacing. The strength and sharpness of the filter is determined by a parameter ${{\epsilon }_{{\rm filter}}}$ such that Lele's $\alpha =1/2-2{{\epsilon }_{{\rm filter}}}$. A value of ${{\epsilon }_{{\rm filter}}}=0$ gives no filtering. The collapse runs used values of ${{\epsilon }_{{\rm filter}}}$ between 0.02 and 0.05; Figure 1 plots the transfer function of the filter for these two values.

Zoom In Zoom Out Reset image size

The gravitational acceleration at each grid-point is obtained by assuming the density factor ${{\rho }_{{\rm g}}}(r)+{{\rho }_{{\rm p}}}(r)-{{\rho }_{{\rm g}0}}$ to be piecewise linear in r and analytically performing the integration in Equation (39) for M(r). Time advancement is performed using a fourth-order Runge–Kutta scheme.

The purely gas dynamic terms in the code, together with the non-reflective condition at the outer boundary, were validated against the exact solution (starting from arbitrary initial conditions) for spherical acoustic waves beautifully described by Rayleigh (1945, Section 279) accounting for reflection at the origin. The initial condition we used was a Gaussian bump for the density and radial velocity perturbation, made symmetric and anti-symmetric, respectively, about the origin. Figure 2(a) shows the result of using the first treatment of the origin (where the origin is a grid-point). The simulation (lines) is in excellent agreement with the exact solution (symbols). To avoid Nyquist oscillations in the acoustic test, we found that the parameter, ${{\epsilon }_{{\rm filter}}}$, governing the strength of the filter and its cut-off wavenumber needed to be set to 0.03. The second scheme (where the origin is a flux surface) required a lower ${{\epsilon }_{{\rm filter}}}=0.01$ in the acoustic test. In the few comparisons we made for collapse cases, both schemes required comparable smoothing. The largest value used for ${{\epsilon }_{{\rm filter}}}$ in the collapse runs was 0.05.

Zoom In Zoom Out Reset image size

Quadrature for evaluating the gravitational acceleration (Equation (39)) was validated using an analytical (${{{\rm cos} }^{2}}(kr)$) density profile. To test the implementation of the particle transport and mutual drag terms, we performed a calculation for a particle–gas acoustic wave in the tightly coupled limit with uniform gas mass fraction (${{\phi }_{{\rm g}}}$) initially. We would like to suggest this as a test case for other two-fluid codes as well. In this limit (see Section 3.2), we should have the same exact solution as a gaseous acoustic wave except that the speed of sound is replaced by ${{c}_{{\rm i}}}\sqrt{{{\phi }_{{\rm g}}}}$. The stopping time was set to ${{t}_{{\rm s}}}=0.001$ and the gas mass fraction was set to ${{\phi }_{{\rm g}}}=1/4$. To obtain a solution with the same functional form as for the gaseous acoustic wave described in the last paragraph, the speed of sound ${{c}_{{\rm i}}}$ was doubled. Figure 2(b) shows good agreement between the simulation (lines) and exact solution (symbols). However, comparing the cyan line and symbols, we see that there is some spurious reflection of the wave at the outer boundary (r = 6). In the problem of gravitational collapse of particles considered later in this work, the outer region of the domain is free of particles and so this issue does not arise.

For the collapse runs, the mesh consists of a uniformly and finely spaced inner region where the density gradient is strongest. This is followed by a second region where the mesh is stretched. Finally, there is a third region where the mesh is uniform again.

The initial clump consists of a core of radius ${{r}_{{\rm core}}}$ having uniform particle density, and a Gaussian tail of length ${{\sigma }_{{\rm tail}}}$. The particle clump is embedded in a gas of initially uniform density ${{\rho }_{{\rm g}0}}$:

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\rho }_{{\rm p}}}(r,0) & = & {{\rho }_{{\rm p}0}}\left\{ \begin{array}{ccccccccccccccc} 1, & r\leqslant {{r}_{{\rm core}}} \\ {{{\rm e}}^{-{{(r-{{r}_{{\rm core}}})}^{2}}/\sigma _{{\rm tail}}^{2}}} & r\gt {{r}_{{\rm core}}}. \\ \end{array} \right., \\ {} & {} & \qquad \qquad {{\rho }_{{\rm g}}}(r,0)={{\rho }_{{\rm g}0}}. \\ \end{array}\end{eqnarray} \tag{ 40 }$$

The choice of a uniform particle density in the core makes the gas–particle wave speed uniform (in the tightly coupled limit) and therefore the wave dynamics is easier to analyze. Furthermore, this choice is in accord with the assumptions of the ODE model. The necessity for a having a smooth tail is purely numerical: because we do not use shock-capturing methods, discontinuities are not allowed. The size of the tail is set to ${{\sigma }_{{\rm tail}}}=0.2{{r}_{{\rm core}}}$ and the initial clump radius is defined as the sum ${{R}_{0}}\equiv {{r}_{{\rm core}}}+{{\sigma }_{{\rm tail}}}$. A comparison with a purely Gaussian ${{\rho }_{{\rm p}}}(r,0)$ profile will be made in Section 4.8. The code requires computation of the particle velocity U from the momentum ${{\rho }_{{\rm p}}}U$. To avoid division by very small numbers a minimum value of 10⁻⁶ was used for ${{\rho }_{{\rm p}}}(r,0)$. The outermost radius of the computational domain is ${{r}_{{\rm max} }}=4$ and the total number of mesh points is $\approx 1500$. The evolution starts from rest. Corresponding to initial conditions at the origin, we define the dynamical (i.e., free-fall) time and the particle stopping time:

$$\begin{eqnarray}&&{{t}_{{\rm dyn}}}\equiv {{\left( {{\rho }_{{\rm p}0}}G \right)}^{-1/2}},\quad {{t}_{{\rm s}}}\equiv {{\left( \beta \bar{c}{{\rho }_{{\rm g}0}} \right)}^{-1}}.\end{eqnarray} \tag{ 41 }$$

The initial condition is characterized by three non-dimensional parameters defined using conditions at the origin (i) The particle loading

$$\begin{eqnarray}&&{{{\Phi }}_{0}}\equiv {{\rho }_{{\rm p}}}(0,0)/{{\rho }_{{\rm g}0}}.\end{eqnarray} \tag{ 42 }$$

We will often instead use the gas mass fraction ${{\phi }_{{\rm g}0}}={{({{{\Phi }}_{0}}+1)}^{-1}}$. (ii) The two-phase Jeans parameter

$$\begin{eqnarray}&&{{J}_{{\rm t}}}\equiv \frac{{{\rho }_{{\rm t}0}}GR_{0}^{2}(1+{{{\Phi }}_{0}})}{c_{{\rm i}}^{2}},\end{eqnarray} \tag{ 43 }$$

where ${{\rho }_{{\rm t}0}}={{\rho }_{{\rm p}0}}+{{\rho }_{{\rm g}0}}$. Since the wave period

$$\begin{eqnarray}&&{{t}_{{\rm wave}}}={{R}_{0}}/{{c}_{{\rm i}}}\phi _{{\rm g}0}^{1/2},\end{eqnarray} \tag{ 44 }$$

we have

$$\begin{eqnarray}&&{{J}_{{\rm t}}}={{({{t}_{{\rm wave}}}/{{t}_{{\rm dyn}}})}^{2}}.\end{eqnarray} \tag{ 45 }$$

(iii) The Stokes number based on dynamical time:

$$\begin{eqnarray}&&{\rm S}{{{\rm t}}_{{\rm dyn}}}\equiv {{t}_{{\rm s}}}/{{t}_{{\rm dyn}}}.\end{eqnarray} \tag{ 46 }$$

Corresponding to the dimensions of mass, length, and time we may set any three quantities to unity in the code. For these we chose R₀, ${{c}_{{\rm i}}}$, and ${{\rho }_{{\rm g}0}}$. The quantity we use to assess the progress of collapse is the radial centroid of particle density:

$$\begin{eqnarray}&&\bar{r}(t)\equiv \frac{\int _{0}^{\infty }r{{\rho }_{{\rm p}}}(r)dr}{\int _{0}^{\infty }{{\rho }_{{\rm p}}}(r)dr}.\end{eqnarray} \tag{ 47 }$$

The first simulation was for the case ${{{\Phi }}_{0}}=100$ and ${\rm S}{{{\rm t}}_{{\rm dyn}}}=0.02$. Figure 3 plots the time evolution of the radial particle centroid $\bar{r}(t)$ for various Jeans parameters ${{J}_{{\rm t}}}$. As ${{J}_{{\rm t}}}\to 0$ the gas becomes incompressible and with spherical symmetry this implies that it must become static. Indeed, as ${{J}_{{\rm t}}}\to 0$ the curves tend to the behavior (black line) for a static gas, which was computed by setting u = 0 in the code. This is the incompressible self-sedimentation limit and will be discussed in more detail later. Monotonic collapse on the dynamical timescale occurs for ${{J}_{{\rm t}}}=0.40$ (violet line). For intermediate values of ${{J}_{{\rm t}}}$, a combination of bouncing and sedimentation is observed prior to collapse. The first bounce occurs due to the propagation of an expansion wave (in both the gas and particle phases) from the particle-lean exterior to the particle-rich interior. The time at which the first minimum in $\bar{r}(t)$ occurs is consistent with wave propagation at the two-phase acoustic speed ${{c}_{{\rm eff}}}=\phi _{{\rm g}}^{1/2}{{c}_{{\rm i}}}$. The wave is then reflected from the origin as a compression, propagates outward to the edge of the particle core where it reflects as a rarefaction, and so on. This sets up an oscillating standing wave whose behavior (in the absence of gravity) will be discussed in more detail in Section 4.5.

Zoom In Zoom Out Reset image size

We saw in the previous sub-section that, in the incompressible limit ${{J}_{{\rm t}}}\to 0$, the gas remains static with uniform density ${{\rho }_{{\rm g}0}}$. In addition, consider the case where the local ratio of stopping time to dynamical time is everywhere small, i.e.,

$$\begin{eqnarray}&&\frac{{{t}_{{\rm s}}}}{{{t}_{{\rm dyn}}}(r)}=\frac{{{(G{{\rho }_{{\rm p}}}(r))}^{1/2}}}{\beta \bar{c}{{\rho }_{{\rm g}0}}}\ll 1.\end{eqnarray} \tag{ 48 }$$

Note that this condition will always fail in the inner core in the very late stages of collapse. If condition (48) holds then particles will be close to terminal velocity where gravity and drag terms balance each other, acceleration being negligible in comparison. The terminal velocity is then given by

$$\begin{eqnarray}&&{{U}_{{\rm term}}}(r)=g(r){{t}_{{\rm s}}}.\end{eqnarray} \tag{ 49 }$$

The mass M_S enclosed by any collapsing material surface S is constant, which makes the gravity at the surface, ${{g}_{S}}\propto {{M}_{S}}/r_{S}^{2}$, increase as the surface collapses. Hence the sedimentation velocity accelerates with time. It was verified that for the static gas simulation in Section 4.2, the particle velocity is closely given by Equation (49) except at late times in the inner core. The particle dilatation associated with Equation (49) is

$$\begin{eqnarray}&&{{{\Delta }}_{{\rm p}}}\equiv \nabla \cdot {\boldsymbol{U}} =-\alpha {{\rho }_{{\rm p}}},\ {\rm where}\ \alpha \equiv 4\pi G{{t}_{{\rm s}}}.\end{eqnarray} \tag{ 50 }$$

The mass conservation equation D ${{\rho }_{{\rm p}}}$/Dt = −${{\rho }_{{\rm p}}}{{{\Delta }}_{{\rm p}}}$ can then be solved in the (assumed) uniform density particle core to give the evolution of particle density:

$$\begin{eqnarray}&&{{\rho }_{{\rm p}}}(t)={{\rho }_{{\rm p}0}}{{\left( 1-\alpha {{\rho }_{{\rm p}0}}t \right)}^{-1}}.\end{eqnarray} \tag{ 51 }$$

According to Equation (51), the volume of the clump ($\propto {{\rho }_{{\rm p}}}$⁻¹), not its radius, decreases linearly with time. A singularity, which may be termed the sedimentation singularity, occurs at

$$\begin{eqnarray}&&{{t}_{{\rm sed},{\rm sing}}}={{(\alpha {{\rho }_{{\rm p}0}})}^{-1}}={{(4\pi {\rm S}{{{\rm t}}_{{\rm dyn}}})}^{-1}}{{t}_{{\rm dyn}}},\end{eqnarray} \tag{ 52 }$$

which will be referred to as the sedimentation time. Note that previous estimates (Cuzzi et al. 2008, 2010) do not account for the acceleration in sedimentation velocity. As a result they overestimate the sedimentation time by a factor of π.

Figure 4 shows that the behavior predicted by Equation (51) agrees well with the static gas simulation at ${\rm S}{{{\rm t}}_{{\rm dyn}}}=0.02$. Obviously, the result (Equation (52)) cannot hold indefinitely as ${\rm S}{{{\rm t}}_{{\rm dyn}}}$ increases since the collapse time can never be smaller than ${{t}_{{\rm dyn}}}$. Deviations for larger ${\rm S}{{{\rm t}}_{{\rm dyn}}}$ will be studied in Section 4.6.

Zoom In Zoom Out Reset image size

Here we develop a simple model for the behavior of the clump near the origin (for the isothermal case). To leading order, the dependent variables behave as follows near the origin:

$$\begin{eqnarray}&&{{\rho }_{{\rm g}}}(r,t)={{\rho }_{{\rm g}}}(t)+\mathcal{O}\left( {{r}^{2}} \right),\end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray}&&{{\rho }_{{\rm p}}}(r,t)={{\rho }_{{\rm p}}}(t)+\mathcal{O}\left( {{r}^{2}} \right),\end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray}&&u(r,t)={{{\Delta }}_{{\rm g}}}(t)r/3+\mathcal{O}\left( {{r}^{3}} \right),\end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray}&&U(r,t)={{{\Delta }}_{{\rm p}}}(t)r/3+\mathcal{O}\left( {{r}^{3}} \right),\end{eqnarray} \tag{ 56 }$$

where ${{{\Delta }}_{{\rm g}}}(t)$ and ${{{\Delta }}_{{\rm p}}}(t)$ denote gas and particle dilatations, respectively, at the origin. Note that a simplified notation is being used in which an argument of t alone (or no argument at all) indicates the value at the origin. One substitutes the forms (53)–(56) into the governing Equations (34)–(37), equates coefficients of r⁰ in the mass conservation equations, and coefficients of r¹ in the momentum equations. However, note that to get the r¹ coefficient of the pressure gradient requires the density to quadratic order; a similar closure problem would present itself no matter how high the order of the expansion. Therefore, the pressure gradient term, which will require a model, is left in its exact form for now. The result is the ODE system:

$$\begin{eqnarray}&&\frac{d{{\rho }_{{\rm g}}}}{dt}=-{{\rho }_{{\rm g}}}{{{\Delta }}_{{\rm g}}},\end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \frac{d{{{\Delta }}_{{\rm g}}}}{dt} & = & -\frac{1}{3}{\Delta }_{{\rm g}}^{2}-4\pi G({{\rho }_{{\rm p}}}+{{\rho }_{{\rm g}}}-{{\rho }_{{\rm g}0}}) \\ {} & {} & +\beta \bar{c}{{\rho }_{{\rm p}}}\left( {{{\Delta }}_{{\rm p}}}-{{{\Delta }}_{{\rm g}}} \right)+{\Pi }(t), \\ \end{array}\end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray}&&\frac{d{{\rho }_{{\rm p}}}}{dt}=-{{\rho }_{{\rm p}}}{{{\Delta }}_{{\rm p}}},\end{eqnarray} \tag{ 59 }$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \frac{d{{{\Delta }}_{{\rm p}}}}{dt} & = & -\frac{1}{3}{\Delta }_{{\rm p}}^{2}-4\pi G({{\rho }_{{\rm p}}}+{{\rho }_{{\rm g}}}-{{\rho }_{{\rm g}0}}) \\ {} & {} & -\beta \bar{c}{{\rho }_{{\rm g}}}\left( {{{\Delta }}_{{\rm p}}}-{{{\Delta }}_{{\rm g}}} \right). \\ \end{array}\end{eqnarray} \tag{ 60 }$$

In Equation (58) the pressure gradient term (exact) is:

$$\begin{eqnarray}&&{\Pi }(t)\equiv -{{\left[ \frac{3}{{{\rho }_{{\rm g}}}(r,t)r}\frac{\partial p}{\partial r} \right]}_{r=0}},\end{eqnarray} \tag{ 61 }$$

where $p(r,t)={{\rho }_{{\rm g}}}(r,t)c_{{\rm i}}^{2}$ for the isothermal case. Two models for ${\Pi }(t)$ will be introduced. The first one (subscript M1) is

$$\begin{eqnarray}&&{\Pi }(t)\approx {{{\Pi }}_{{\rm M}1}}(t)=-\frac{3({{\rho }_{{\rm g}0}}-{{\rho }_{{\rm g}}}(t))c_{{\rm i}}^{2}}{{{\rho }_{{\rm g}}}(t)f_{{\rm rg}}^{2}r_{{\rm g}}^{2}(t)},\end{eqnarray} \tag{ 62 }$$

where ${{f}_{{\rm rg}}}$ is a modeling constant. The quantity ${{r}_{{\rm g}}}(t)$ is the radius of the region of compressed gas and is obtained by mass conservation:

$$\begin{eqnarray}&&{{r}_{{\rm g}}}(t)={{\left[ {{\rho }_{{\rm g}0}}/{{\rho }_{{\rm g}}}(t) \right]}^{1/3}}{{r}_{{\rm g}}}(0),\end{eqnarray} \tag{ 63 }$$

with its initial value set equal to the initial radius of the particle clump, i.e., ${{r}_{{\rm g}}}(0)={{R}_{0}}$. The model constant ${{f}_{{\rm rg}}}$ in Equation (62) was calibrated to ${{f}_{{\rm rg}}}=0.43$ by obtaining the closest match to one simulation case (${{{\Phi }}_{0}}=100,{{J}_{{\rm t}}}=0.3,{\rm S}{{{\rm t}}_{{\rm dyn}}}=0.02$). Figure 5(a) shows the the time evolution of the inverse particle density at the origin as given by the simulations for various ${{J}_{{\rm t}}}$, while Figure 5(b) shows the corresponding predictions of the model. The ODE model reproduces many features of the simulations, however, the depth of collapse preceding the first bounce is smaller. It was suspected that this difference arises because the pressure gradient model (Equation (62)) fails to account for the time delay in arrival of the pressure wave at the origin. This is confirmed in Figure 6 which shows that the actual (solid line) pressure gradient term ${\Pi }(t)$ at the origin begins to rise after a delay compared to the model (dashed) expressions in Equations (62) and (63). The second model allows for such a delay by tracking a wavefront location ${{r}_{{\rm front}}}(t)$ as follows:

$$\begin{eqnarray}&&\frac{d{{r}_{{\rm front}}}(t)}{dt}=u({{r}_{{\rm front}}},t)-{{c}_{{\rm eff}}},\end{eqnarray} \tag{ 64 }$$

with an initial location set to the particle core boundary (defined in Equation (40)):

$$\begin{eqnarray}&&{{r}_{{\rm front}}}(0)={{r}_{{\rm core}}}.\end{eqnarray} \tag{ 65 }$$

The total wave speed on the right of Equation (64) is that of the inward propagating characteristic with the gas speed given from Equation (55)

$$\begin{eqnarray}&&u({{r}_{{\rm front}}},t)={{{\Delta }}_{{\rm g}}}(t){{r}_{{\rm front}}}/3,\end{eqnarray} \tag{ 66 }$$

and wave speed (relative to the gas) given by

$$\begin{eqnarray}&&{{c}_{{\rm eff}}}(t)={{\left( \frac{{{\rho }_{{\rm g}}}(t)}{{{\rho }_{{\rm g}}}(t)+{{\rho }_{{\rm p}}}(t)} \right)}^{1/2}}{{c}_{{\rm i}}}.\end{eqnarray} \tag{ 67 }$$

Let ${{t}_{{\rm arrival}}}$ be the time that the wave arrives at the origin. Then the pressure gradient term is turned on only for $t\gt {{t}_{{\rm arrival}}}$, i.e.,

$$\begin{eqnarray}&&{\Pi }(t)\approx {{{\Pi }}_{{\rm delay}}}(t)=H(t-{{t}_{{\rm arrival}}}){{{\Pi }}_{{\rm M}1}}(t),\end{eqnarray} \tag{ 68 }$$

where H(t) is the Heaviside function. ${{{\Pi }}_{{\rm delay}}}(t)$ constitutes our second model. Calibration for the same case as before (${{J}_{{\rm t}}}=0.3,{{{\Phi }}_{0}}=100,{\rm S}{{{\rm t}}_{{\rm dyn}}}=0.02$) gave ${{f}_{{\rm rg}}}=0.4579$ for the model constant. The extra digits of precision result from the fact that we tried our best to nail the ${{J}_{{\rm t}}}=0.3$ case, which lies in a sensitive region of the parameter space with a very deep compression followed by a bounce. Figure 5(c) shows that the ODE model with arrival delay gives a better reproduction of the time and depth of the initial compression, however, the bounces are higher than the simulations for ${{J}_{{\rm t}}}\gtrsim 0.3$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It may be noted that the ODE model can be written in a form bearing an analogy to that of two coupled (non-linear) oscillators. For a velocity field of the form in Equation (55) the acceleration of a gas particle is:

$$\begin{eqnarray}&&\frac{{{d}^{2}}{{r}_{{\rm g}}}}{d{{t}^{2}}}=\frac{Du}{Dt}=\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial r}=\frac{{{r}_{{\rm g}}}}{3}\left( \frac{d{{{\Delta }}_{{\rm g}}}}{dt}+\frac{{\Delta }_{{\rm g}}^{2}}{3} \right),\end{eqnarray} \tag{ 69 }$$

and similarly for a solid particle. Hence, Equations (58) and (60) become

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \frac{{{d}^{2}}{{r}_{{\rm g}}}}{d{{t}^{2}}} & = & -\frac{4\pi }{3}G\left( {{\rho }_{{\rm p}}}+{{\rho }_{{\rm g}}}-{{\rho }_{{\rm g}0}} \right){{r}_{{\rm g}}}+\beta \bar{c}{{\rho }_{{\rm p}}}(U-u) \\ {} & {} & +\frac{{\Pi }(t){{r}_{{\rm g}}}}{3}, \\ \end{array}\end{eqnarray} \tag{ 70 }$$

$$\begin{eqnarray}&&\frac{{{d}^{2}}{{r}_{{\rm p}}}}{d{{t}^{2}}}=-\frac{4\pi }{3}G\left( {{\rho }_{{\rm p}}}+{{\rho }_{{\rm g}}}-{{\rho }_{{\rm g}0}} \right){{r}_{{\rm p}}}-\beta \bar{c}{{\rho }_{{\rm g}}}(U-u).\end{eqnarray} \tag{ 71 }$$

The densities on the right of these equations are the following functions of ${{r}_{{\rm g}}}$ and ${{r}_{{\rm p}}}$ due to mass conservation:

$$\begin{eqnarray}&&{{\rho }_{{\rm p}}}(t)={{\left( \frac{{{r}_{{\rm p}}}(0)}{{{r}_{{\rm p}}}(t)} \right)}^{3}}{{\rho }_{{\rm p}}}(0),\quad {{\rho }_{{\rm g}}}(t)={{\left( \frac{{{r}_{{\rm g}}}(0)}{{{r}_{{\rm g}}}(t)} \right)}^{3}}{{\rho }_{{\rm g}}}(0).\end{eqnarray} \tag{ 72 }$$

The terms in Equations (70) and (71) that are proportional to ${{r}_{{\rm g}}}$ and ${{r}_{{\rm p}}}$ represent spring-like (non-linear) restoring forces and the drag terms provide damping.

To develop insight into the nature of bouncing oscillations, gravity was activated for only a short initial period $t/{{t}_{{\rm dyn}}}\leqslant 0.01$ to provide an initial velocity, after which it was switched off. Figure 7(a) shows profiles of gas density for one period of the standing wave.

Zoom In Zoom Out Reset image size

Figure 7(b) shows that the ODE model (dashed line) reproduces the period and damping rate of density fluctuations at the origin quite well. This gives us confidence to use it to develop further insight. Linearizing the ODE model by substituting

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{r}_{{\rm g}}}(t) & = & {{r}_{{\rm g}}}(0)\left[ 1+{{\epsilon }_{{\rm g}}}(t)+\mathcal{O}\left( \epsilon _{{\rm g}}^{2} \right) \right], \\ {{r}_{{\rm p}}}(t) & = & {{r}_{{\rm p}}}(0)\left[ 1+{{\epsilon }_{{\rm p}}}(t)+\mathcal{O}\left( \epsilon _{{\rm p}}^{2} \right) \right] \\ \end{array}\end{eqnarray} \tag{ 73 }$$

into the oscillator Equations (70) and (71) (without the gravity terms and considering times larger than the pressure gradient arrival time) gives to leading order:

$$\begin{eqnarray}&&{{\ddot{\epsilon }}_{{\rm g}}}=t_{{\rm s}}^{-1}{{{\Phi }}_{0}}\left( {{{\dot{\epsilon }}}_{{\rm p}}}-{{{\dot{\epsilon }}}_{{\rm g}}} \right)-{{\omega }^{2}}{{\epsilon }_{{\rm g}}},\end{eqnarray} \tag{ 74 }$$

$$\begin{eqnarray}&&{{\ddot{\epsilon }}_{{\rm p}}}=-t_{{\rm s}}^{-1}\left( {{{\dot{\epsilon }}}_{{\rm p}}}-{{{\dot{\epsilon }}}_{{\rm g}}} \right).\end{eqnarray} \tag{ 75 }$$

The frequency

$$\begin{eqnarray}&&\omega \equiv {{\left( \frac{3c_{{\rm i}}^{2}}{f_{{\rm rg}}^{2}r_{{\rm p}}^{2}(0)} \right)}^{1/2}},\end{eqnarray} \tag{ 76 }$$

is inversely proportional to the sound crossing time and arises from the pressure gradient model. The system of Equations (74) and (75) have eigenvalues ${{\lambda }_{m}}$ ($m=0,\ldots ,3$) each of which corresponds to a mode with time behavior ${\rm exp} ({{\lambda }_{m}}t)$. The first eigenvalue is zero and the other three are roots of the cubic polynomial

$$\begin{eqnarray}&&{{\lambda }^{3}}+t_{{\rm s}}^{-1}{{\phi }_{{\rm g}0}}{{\lambda }^{2}}+{{\omega }^{2}}\lambda +{{\omega }^{2}}t_{{\rm s}}^{-1}=0.\end{eqnarray} \tag{ 77 }$$

In working out the roots of the cubic, the non-dimensional parameter

$$\begin{eqnarray}&&\varepsilon \equiv {{\left( {{\phi }_{{\rm g}0}}\omega {{t}_{{\rm s}}} \right)}^{2}}\end{eqnarray} \tag{ 78 }$$

emerges. It involves the stopping time and will later be seen to set the damping rate of oscillations. It can be expressed as

$$\begin{eqnarray}&&\varepsilon =\frac{3}{f_{{\rm rg}}^{2}}\phi _{{\rm g}0}^{2}{{\left( \frac{{{t}_{{\rm s}}}}{{{t}_{c}}} \right)}^{2}},\end{eqnarray} \tag{ 79 }$$

where ${{t}_{c}}={{r}_{{\rm p}}}(0)/{{c}_{{\rm i}}}$ is the sound crossing time. For the example depicted in Figure 7 we have $\varepsilon \sim {{10}^{-3}}$. To leading order in $\varepsilon$ and the initial gas mass fraction ${{\phi }_{{\rm g}0}}$, the non-zero eigenvalues are

$$\begin{eqnarray}&&{{\lambda }_{1}}=-\omega {{\varepsilon }^{-1/2}},\end{eqnarray} \tag{ 80 }$$

$$\begin{eqnarray}&&{{\lambda }_{2,3}}=\omega \left[ -\frac{1}{2}{{\varepsilon }^{1/2}}\pm i{{\left( {{\phi }_{{\rm g}0}}-\varepsilon /4 \right)}^{1/2}} \right].\end{eqnarray} \tag{ 81 }$$

The first eigenvalue gives a rapidly damped non-oscillating mode. If

$$\begin{eqnarray}&&{{\phi }_{{\rm g}0}}\gt \varepsilon /4,\end{eqnarray} \tag{ 82 }$$

then the second two eigenvalues correspond to a damped oscillation as realized in Figure 7. This is analogous to the underdamped case for a harmonic oscillator. The oscillation frequency is given by

$$\begin{eqnarray}&&{{\omega }_{{\rm osc}}}=\omega {{\left( {{\phi }_{{\rm g}0}}-\varepsilon /4 \right)}^{1/2}}.\end{eqnarray} \tag{ 83 }$$

Without the $\varepsilon /4$ term, the result corresponds to the frequency of waves propagating at speed c_eff. The $\varepsilon /4$ term gives a frequency correction due to damping analogous to that for a damped spring–mass system. If condition (82) is not satisfied then all three modes are damped without oscillation. This is analogous to the overdamped case of the harmonic oscillator. We expect that condition (82) will also roughly determine whether bouncing oscillations will be present in the presence of gravity. For example for ${{J}_{{\rm t}}}=0.20$ and ${{{\Phi }}_{0}}=100$, Equation (82) predicts that bouncing should disappear for ${\rm S}{{{\rm t}}_{{\rm dyn}}}\gt .13$. Indeed, we found no bouncing in the simulations for ${\rm S}{{{\rm t}}_{{\rm dyn}}}=0.15$.

For the case of Figure 7, Equation (81) gives an oscillation period of $0.53{{t}_{{\rm dyn}}}$ and an ${{{\rm e}}^{-1}}$ damping period of $0.70{{t}_{{\rm dyn}}}$. The corresponding values in the simulation are $0.59{{t}_{{\rm dyn}}}$ and $0.96{{t}_{{\rm dyn}}},$ respectively.

Collapse time ${{t}_{{\rm collapse}}}$ was recorded in the code as the instant when the normalized density centroid $\bar{r}/{{R}_{0}}$ falls below 0.03. The plot of $\bar{r}(t)$ becomes very nearly vertical at this time and so the actual instant of the singularity $\bar{r}/{{R}_{0}}\to 0$ is very close to the recorded value.

We begin by investigating dependence on ${\rm S}{{{\rm t}}_{{\rm dyn}}}$, fixing the particle loading at ${{{\Phi }}_{0}}=100$. Figure 8(a) shows that on average (i.e., apart from oscillations in slope, which become quite pronounced for $0.30\lesssim {{J}_{{\rm t}}}\lesssim 0.40$) the curves fan out linearly with respect to ${\rm S}{{{\rm t}}_{{\rm dyn}}}$⁻¹ from the point $({\rm St}_{{\rm dyn}}^{-1}=1.0,{{t}_{{\rm collapse}}}/{{t}_{{\rm dyn}}}=0.57)$. The dashed lines in Figure 8(a) represent a simple linear (with cut-off) formula (Equation (85)) to be developed below. Figures 8(b) and (c) show that both versions of the ODE model provide reasonably good quantitative predictions. Figure 9 shows the simulation result for a smaller particle loading of ${{{\Phi }}_{0}}=10$. Comparing it with Figure 8(a), one concludes that while details of the slope oscillations change, the overall behavior is insensitive to Φ₀ in the range $10\leqslant {{{\Phi }}_{0}}\leqslant 100$ (keeping ${{J}_{{\rm t}}}$ fixed). One also notes that in the incompressible (static gas) limit ${{J}_{{\rm t}}}\to 0$ (solid black line), the terminal velocity approximation ${{t}_{{\rm collapse}}}/{{t}_{{\rm dyn}}}={{(4\pi {\rm S}{{{\rm t}}_{{\rm dyn}}})}^{-1}}$ (red dashed line) is good to within 10% for ${\rm S}{{{\rm t}}_{{\rm dyn}}}\lt 0.05$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next, the dependence on ${{J}_{{\rm t}}}$ is investigated; see Figure 10. Again, there are oscillations in slope as well as a few cases where increasing the strength of gravity (i.e., ${{J}_{{\rm t}}}$) actually increases the collapse time. This is a consequence of bouncing: for these cases, increasing ${{J}_{{\rm t}}}$ causes a deeper initial collapse but the bounce is even stronger, which increases the time to collapse. Apart from such anomalies and oscillations, the behavior is linear. We write an estimate for the collapse time by drawing a straight line from the terminal velocity approximation at the left end (${{J}_{{\rm t}}}=0,{{t}_{{\rm collapse}}}/{{t}_{{\rm dyn}}}={\rm max} [{{(4\pi {\rm S}{{{\rm t}}_{{\rm dyn}}})}^{-1}},0.65]$) to the point $({{J}_{{\rm t}}}=0.43,{{t}_{{\rm collapse}}}/{{t}_{{\rm dyn}}}=0.65)$. The max function has been introduced in writing the coordinates of the left end-point in order to prevent the collapse time from dropping below $0.65{{t}_{{\rm dyn}}}$ when the terminal velocity approximation breaks down for large ${\rm S}{{{\rm t}}_{{\rm dyn}}}$. The straight line is given by

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} F({{J}_{{\rm t}}}) & = & {\rm max} [{{\left( 4\pi {\rm S}{{{\rm t}}_{{\rm dyn}}} \right)}^{-1}},0.65] \\ {} & {} & \times \left[ 1-{{J}_{{\rm t}}}/0.43 \right]+(0.65/0.43){{J}_{{\rm t}}}. \\ \end{array}\end{eqnarray} \tag{ 84 }$$

Finally we write

$$\begin{eqnarray}&&\frac{{{t}_{{\rm collapse}}}}{{{t}_{{\rm dyn}}}}={\rm max} (F({{J}_{{\rm t}}}),0.65),\end{eqnarray} \tag{ 85 }$$

to keep the collapse time from going below (0.65 times) the dynamical time for large ${{J}_{{\rm t}}}$. Figure 10 shows formula (85) using dashed lines. The reader may also refer back to Figures 8(a) and 9 and the dashed lines therein to assess the performance of formula (85) as a function of ${\rm S}{{{\rm t}}_{{\rm dyn}}}$⁻¹. Comparing the results for ${{{\Phi }}_{0}}=10$ and 100 one concludes again that the dependence on particle loading Φ₀ as a separate parameter is weak and Equation (85) continues to work quite well. In other words, the presence of Φ₀ in the definition ${{J}_{{\rm t}}}$ captures most of the influence of Φ₀.

Zoom In Zoom Out Reset image size

Figure 11 shows the predictions of collapse time given by both versions of the ODE model. The results of the simulations, previously shown in Figure 10, are depicted in gray. The overall agreement is satisfactory for both models. The model without delay gives a better prediction for smaller ${{J}_{{\rm t}}}$ while the model with delay is slightly more accurate for ${{J}_{{\rm t}}}\gtrsim 0.4$.

Zoom In Zoom Out Reset image size

Time evolution of the particle density profile was studied for the following set of 12 cases: ${{{\Phi }}_{0}}\in \{10,100\}$ $\;\times {\rm S}{{{\rm t}}_{{\rm dyn}}}\in \{0.005,0.05\}$ $\;\times {{J}_{{\rm t}}}\in \{0,0.2,0.4\}$; see Figure 12. Profiles are compared at the same value of the "progress variable" $\bar{r}(t)$ (the density centroid defined in Equation (47)). Insets are used to show $\bar{r}(t)$. Since $\bar{r}(t)$ sometimes oscillates in time, it is necessary to state that profiles are recorded the first time that $\bar{r}(t)$ goes below the given value. Figures 12(a)–(c) illustrate the ${{J}_{{\rm t}}}$ dependence, other parameters being kept fixed. For ${{J}_{{\rm t}}}=0$ (incompressible static gas), the profile remains flat in the collapsing core. This is in accord with the solution obtained in Section 4.3 where the terminal velocity approximation was invoked. Slight departure from this behavior at the final instant is likely due to breakdown of this approximation. Note that as ${{J}_{{\rm t}}}$ increases (Figures 12(a)–(c)), the profiles become less flat in the core. In particular, when ${{J}_{{\rm t}}}$ is increased from 0.2 to 0.4, the knee in some profiles (e.g., for $\bar{r}/{{R}_{0}}=0.25,$ violet) disappears. In all the other cases this was the only difference observed between the ${{J}_{{\rm t}}}=0.2$ and 0.4 cases. Hence, we omit all other plots for ${{J}_{{\rm t}}}=0.4$.

Zoom In Zoom Out Reset image size

Next compare Figures 12(b) and (d); in the latter case Φ₀ is 10 times larger. Except for a factor of 10 difference in the ordinate, the profiles in the two plots are nearly the same. This was observed to be true for all the other cases in the set.

Between Figures 12(d) and (e), ${\rm S}{{{\rm t}}_{{\rm dyn}}}$ has been increased by 10. Despite the large difference in the collapse time, profiles at the same value of the progress variable are nearly identical. This conclusion holds for all the other cases in the set.

To determine the functional form of the particle density profile about the origin when the singularity forms there, we plot in Figure 12(e) profiles for all 12 cases on log–log axes at the final instant of the simulations, when $\bar{r}/{{R}_{0}}=0.03$. A power law region develops with a slope slightly steeper than −2.

The original (uniform-with-tail) initial particle density profile (Equation (40)) was chosen for analytical advantage: it is consistent with the Taylor approximation (Equations (53)–(56)) of the ODE model over a longer region near the origin and for a longer time. Furthermore, the speed of a particle–gas wave is uniform, making the physics easier to analyze. Here, we explore some differences when the Gaussian initial profile,

$$\begin{eqnarray}&&{{\rho }_{{\rm p}}}(r,0)={{\rho }_{{\rm p}0}}{\rm exp} (-{{r}^{2}}/\sigma _{{\rm G}}^{2}),\end{eqnarray} \tag{ 86 }$$

is used instead. The definitions, Equations (41) and (42), for particle loading, dynamical time, and stopping time remain the same as before. Recall that they are based on conditions at the origin. The definition of the two-phase Jeans parameter, however, requires some comment. Recall the definition (Equation (43)) used for the simulations:

$$\begin{eqnarray}&&{{J}_{{\rm t}}}\equiv \frac{{{\rho }_{{\rm t}0}}GR_{0}^{2}(1+{{{\Phi }}_{0}})}{c_{{\rm i}}^{2}},\end{eqnarray} \tag{ 87 }$$

where ${{R}_{0}}={{R}_{{\rm core}}}+{{\sigma }_{{\rm tail}}}$ for the original profile. We need to rewrite the clump radius R₀ so it has a common definition independent of profile shape. The obvious way to do this (and in retrospect what should have been done in the first place) is to use the total particle mass. For the original profile (where ${{\sigma }_{{\rm tail}}}=0.2\;{{R}_{0}}$) the total particle mass is ${{M}_{{\rm p}}}=4.017{{\rho }_{{\rm p}0}}\;R_{0}^{3}$. Hence we define a mass-effective clump radius

$$\begin{eqnarray}&&{{R}_{0,{\rm eff}}}\equiv {{\left( \frac{{{M}_{{\rm p}}}}{4.017{{\rho }_{{\rm p}0}}} \right)}^{1/3}}.\end{eqnarray} \tag{ 88 }$$

and rewrite the definition of ${{J}_{{\rm t}}}$ as

$$\begin{eqnarray}&&{{J}_{{\rm t}}}\equiv \frac{{{\rho }_{{\rm t}0}}GR_{0,{\rm eff}}^{2}(1+{{{\Phi }}_{0}})}{c_{{\rm i}}^{2}}.\end{eqnarray} \tag{ 89 }$$

The Gaussian profile has ${{M}_{{\rm p}}}={{\rho }_{{\rm p}0}}{{\pi }^{3/2}}\sigma _{{\rm G}}^{3}$ which is substituted into Equation (88) to obtain ${{R}_{0,{\rm eff}}}$. Code units for the Gaussian case are such that ${{\sigma }_{{\rm G}}}={{c}_{{\rm i}}}={{\rho }_{{\rm g}0}}=1$.

Figure 13 shows collapse time as a function of ${{J}_{{\rm t}}}$. Comparing this with Figure 10(a) for the uniform core case, one sees that while in the incompressible (i.e., sedimentation) limit the collapse time is the same (i.e., ${{t}_{{\rm collapse}}}/{{t}_{{\rm dyn}}}$ very nearly equals ${{(4\pi {\rm S}{{{\rm t}}_{{\rm dyn}}})}^{-1}}$), as ${{J}_{{\rm t}}}$ increases, the Gaussian profile collapses more slowly. This is probably because the mass function ${{M}_{{\rm p}}}(r)$ near the origin increases more slowly with r for the Gaussian case. A quadratic function of ${{J}_{{\rm t}}}$ was found necessary to fit the data:

$$\begin{eqnarray}&&\frac{{{t}_{{\rm collapse}}}}{{{t}_{{\rm dyn}}}}={\rm max} ({{F}_{{\rm G}}}({{J}_{{\rm t}}}),0.75),\end{eqnarray} \tag{ 90 }$$

with

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{F}_{{\rm G}}}({{J}_{{\rm t}}}) & \equiv & {\rm max} [{{(4\pi {\rm S}{{{\rm t}}_{{\rm dyn}}})}^{-1}},0.75] \\ {} & {} & \times \;\left[ 1-0.25{{J}_{{\rm t}}}-0.75J_{{\rm t}}^{2} \right]. \\ \end{array}\end{eqnarray} \tag{ 91 }$$

This fit is shown using dashed lines in Figure 13. It should be noted that the Gaussian case also exhibited bouncing in certain cases.

Zoom In Zoom Out Reset image size

Using formulas (85) and (90), collapse time was calculated for clumps located at the midplane of a minimum-mass solar nebula (Hayashi 1981) for which the parameters are:

$$\begin{eqnarray}&&{{c}_{{\rm i}}}={{10}^{5}}{{(\varpi /{\rm AU})}^{-1/4}}\;{\rm cm}\ {{{\rm s}}^{-1}},\end{eqnarray} \tag{ 92 }$$

$$\begin{eqnarray}&&{\rm gas}\ {\rm scale}\ {\rm height}\ {{H}_{{\rm g}}}=\sqrt{2}{{c}_{{\rm i}}}/{{{\Omega }}_{{\rm Kepler}}},\end{eqnarray} \tag{ 93 }$$

$$\begin{eqnarray}&&{\rm surface}\ {\rm density}\ {{{\Sigma }}_{{\rm g}}}=1700{{(\varpi /{\rm AU})}^{-3/2}}\;{\rm gm}\ {\rm c}{{{\rm m}}^{-2}},\end{eqnarray} \tag{ 94 }$$

$$\begin{eqnarray}&&{{\rho }_{{\rm g}0}}={{{\Sigma }}_{{\rm g}}}/({{\pi }^{1/2}}{{H}_{{\rm g}}}).\end{eqnarray} \tag{ 95 }$$

The material density of solids is taken to be ${{\rho }_{{\rm material}}}=3$ gm cm⁻³. Eight cases were selected: particle radius $a\in \{0.5\ {\rm mm},5\ {\rm mm}\}$ × $\varpi \in \{3\ {\rm AU},30\ {\rm AU}\}\times {{{\Phi }}_{0}}{\mkern 1mu} \in \{10,$ $100\}$; validity of the Epstein drag regime was verified for each case. Results are shown in Figure 14 as a function of normalized clump radius ${{R}_{0,{\rm eff}}}/{{H}_{{\rm g}}}$. Solid lines are for the uniform + tail initial profile and dashed lines are used for the Gaussian profile. As ${{R}_{0,{\rm eff}}}/{{H}_{{\rm g}}}\to 0$, the collapse time tends to the sedimentation time which is independent of the clump size. Note that at a fixed radial location, the sedimentation time depends on the product $a{{{\Phi }}_{0}}$. For ${{{\Phi }}_{0}}=100$, the clump radius has to be a few percent of the disk scale height for the Jeans parameter ${{J}_{{\rm t}}}$ to be sufficiently large that the collapse time is significantly reduced below the sedimentation time. On the other hand, for ${{{\Phi }}_{0}}=10$ the clump radius has to be significant fraction of ${{H}_{{\rm g}}}$ for this to happen. For the Gaussian initial profile (dashed lines) the size required to overcome the sedimentation limit is even greater. Currently, no mechanism is known that forms such large dense clumps. We conclude that contraction of clumps containing mm to cm-sized particles will occur on the sedimentation timescale ${{t}_{{\rm sed},{\rm sing}}}={{(4\pi G{{\rho }_{{\rm p}0}}{{t}_{{\rm s}}})}^{-1}}$. Recall that this is a factor of π smaller than earlier estimates (Cuzzi et al. 2008, 2010).

Zoom In Zoom Out Reset image size