Home Back

2.a The free-fall collapse of a cloud without pressure

The objective of this section is to show and analyze the following facts:

These two facts are inter-related with a third, which is

This fact can be derived very succintly[5] from Newton's theorem, which establishes that the gravitational potential within the cavity inside a constant density homeoid is constant. A homeoid is the shape between two concentric, homotetic ellipsoids. The potential inside an uniform ellipsoid is calculated explicitly in [1].

The problem of collapsing, pressureless gas or "dust" (as sometimes it is called) has one peculiarity not shared by other hydrodynamics problems: without pressure, there are no dynamical constrains which prevent converging or intersecting streamlines, and the path of dust particles can self-intersect. This is not physical because, at least on small scales, we expect pressure to prevent gas parcels from getting too close to each other.

First we derive a property of the contracting motions induced by self-gravity in this context. We stress under the conditions analyzed, the contracting motions are unpposed by any pressure or by any force whatsoever. Consider an homogeneous distribution of fluid mass at rest at \(t=0\), that is, \(\rho(\vec{x},0)=\rho_0\) and \(u(\vec{x},0)=0\). Assume that the dynamics of the fluid does not introduce vorticity, as it is the case if the body forces are conservative, and there is no pressure or viscosity. Then, the density will preserve homogeneity if and only if the velocity field of the fluid is affine linear in space. That is,

We drop the arrow on top of vector quantities when there is no risk of confusion.

\(\boxed{\Longleftarrow}\)

From the conservation of mass (or mass continuity) $$ \frac{D\ln\rho}{Dt}=-\nabla\cdot u = -\mathrm{Tr}\Lambda(t)~~,$$ where \(D/Dt\) is the convective derivative. The initial condition of all streamlines being that the density is \(\rho_0\), we conclude that the density is also the same and independent of \(\vec{x}\) after time \(t\).

\(\boxed{\Longrightarrow}\)

By the equation of continuity, we have that $$ \frac{d\ln\rho(t)}{dt}=-\nabla\cdot u=\xi(t)~~.\tag{1.1}$$ On the other hand, the momentum equation \begin{align*} \partial_tu+(u\cdot \nabla)u&=-\nabla\phi &&/\nabla\cdot\\ -\dot{\xi}(t) +\nabla\cdot\left((u\cdot \nabla)u\right)&=-4\pi G\rho(t) &&/\nabla \tag{1.2}\\ \rightarrow&\nabla\left(\nabla\cdot\left((u\cdot \nabla)u\right)\right)=0~~,&\tag{1.3} \end{align*} where the right side of (1.2) comes from \(\triangle \phi=4\pi G\rho\). Equation (1.3) is equivalent to \begin{align*} \require{cancel} 0=\partial_k\left(\partial_i u_j(\partial_j u_i)\right)&=\partial_k\left((\partial_i u_j)(\partial_j u_i)+u_j(\partial_{ij}u_i)\right)\\ &=\partial_k\left((\partial_i u_j)(\partial_j u_i)+u_j\cancelto{0}{(\partial_j\xi(t))}\right) \\ &=(\partial_{ki}u_j)(\partial_ju_i)+(\partial_iu_j)(\partial_{kj}u_i)\tag{1.4}\\ \end{align*} At this point we use that \(u\) is irrotational, because the initial condition (\(u=0\)) does not have vorticity anywhere and all the forces involved are conservative. Therefore, we can assume that \(u_i=\partial_i w\), where \(w\) is the velocity potential. This implies that \(\partial_iu_j=\partial_ju_i\). Applying this to (1.4), we deduce \begin{equation*} (\partial_{ijk}w)(\partial_{ij}w)=0~~,\tag{1.5} \end{equation*} or, what is the same, that $$ (\partial_{ij}w)(\partial_{ij}w)=\sum_{i,j=1}^3 \left(\frac{\partial^2w}{\partial x_i\partial x_j}\right)^2=\dot{\xi}(t)-4\pi G \rho(t)~~,\tag{1.6}$$ where \( (\dot{\xi}(t)-4\pi G \rho(t)) \) is independent of the spatial coordinates by assumption. Because \(\nabla\cdot u=\nabla \cdot \nabla w=\Delta w =-\xi(t)\), we have that \(\Delta(\partial_i w)=0\). Higher partial derivatives of \(w\) are also harmonic functions[4], including \(\partial_{ij}w\). Therefore, from Equation (1.6): \begin{align*} 0&=\left(\partial_k(\partial_{kij}w)(\partial_{ij}w)\right)\\ &=\cancelto{0}{\partial_{kk}(\partial_{ij}w)}+(\partial_{kij}w)(\partial_{kij}w)=\sum_{i,j,k=1}^n\left(\frac{\partial^3w}{\partial x_k\partial x_i\partial x_j}\right)^2\\ &\rightarrow \frac{\partial^3w}{\partial x_k\partial x_i\partial x_j}=0~~,\tag{1.7} \end{align*} independent of \(i,j,\) and \(k\). That is, the function \(w\) depends at most quadratically on the coordinates, or that $$ w(\vec{x},t)=\frac{1}{2}\vec{x}^\intercal \Lambda(t) \vec{x}+\vec{a}(t)\cdot\vec{x}+\alpha(t)~~.\tag{1.8} $$ Note that \(\Lambda\) can be taken to be symmetric because \(u=\nabla w=\frac{1}{2}(\Lambda^\intercal+\Lambda)x\). Related with this, although expressions in the velocity of the type \(\vec{\Omega}(t)\times \vec{x}\) do preserve the density and are linear with the coordinates, they are not allowed when the velocity comes from a potential because they represent vorticity (or maybe rotation of the system of reference). They do not exist when \(\Lambda(t)\) is symmetric either. Also the constant, solid displacement \(\vec{a}(t)\) preserves the density, but in general we will drop this term as it cannot come from self-gravity. However, it could arise, for example, from an accelerating system of reference. The rest of the dynamics would be unchanged because adding a linear term to the gravitational potential does not affect Poisson's equation.

\(\blacksquare\)

A homotetic velocity field as the one defined in Equation \eqref{eq-A} guarantees the streamlines do not cross, at least as long the matrix \(\Lambda(t)\) is well defined and varies continuously through time. What happens during collapse is that the matrix \(\Lambda(t)\) "explodes" to infinity in finite time. This time correspond to the gravitational collapse time, or free-fall time. At this time, streamlines do cross.

The explicit equations of motion of the fluid are the following \begin{align} & & u=\nabla w& \\ (\partial_t+u\cdot\nabla)\ln(\rho)&=-\nabla\cdot u & \partial_t\rho+\nabla\rho\cdot\nabla w+\rho\Delta w&=0 \label{eq-cont}\\ \partial_tu+(u\cdot\nabla)u&=-\nabla\phi & \partial_tw+\frac{1}{2}\left(\nabla w\right)^2+\phi&=0\label{eq-mom}\\ &\qquad\qquad \Delta \phi=4\pi G \rho~~, &&\label{eq-poisson}\\ \end{align} where in the left side we leave the velocity and in the right side we use the velocity potential.[2]

2a.1 The ellipsoid is the only shape which could collapse homogeneously

Using \eqref{eq-A}, it is rather straightforward to derive this result. An homogeneous density requires the velocity potential to be $$w=\frac{1}{2}x^\intercal \Lambda(t) x~~.$$ Therefore, from Equation \eqref{eq-mom} we deduce that \begin{equation} \phi=\frac{1}{2}x^\intercal (-\dot{\Lambda}(t)-\Lambda^2(t))x~~,\label{eq-ficc} \end{equation} that is, \(\phi\) is squared in the coordinates throughout time and necessarily the potential must arise from an ellipsoidal, uniform distribution of mass. Therefore, other possible shapes with the property that the unopossed free-fall collapse maintains constant density are not possible. However, it is important to verify that these ellipsoids represent a coherent solution to all equations.

2a.2 The contracting ellipsoid is a solution of the hydrodynamic equations

We will look for solutions of ellipsoids with axes aligned with the coordinate axes. Due to the non-rotational initial conditions, we expect that the later ellipsoids would remain also aligned with the coordinate axes. Under these circumstances, we expect matrix \(\Lambda(t)\) to be diagonal (or at least this is the solution we will look for). That is, \begin{equation} \Lambda(t)=-\begin{pmatrix} \lambda_1(t) &0&0\\0& \lambda_2(t) &0\\0&0&\lambda_3(t)\\ \end{pmatrix}\quad\rightarrow\quad-\nabla\cdot u=\sum_{i=1}^3\lambda_i(t)=\frac{d\ln\rho}{dt}~~.\label{eq-contE} \end{equation} Each one of the \(\lambda_i(t)\) characterizes the contraction in the \(\hat{x}_i\) axis of the ellipsoid, which we will assume have initial semimajor axes \(a_i(0)\). The semiaxes in later times are given by \(a_i(t)\), that is, the ellipsoid fills the volume given by \(\lbrace(x_1,x_2,x_3):\frac{x_1^2}{a_1^2}+\frac{x_2^2}{a_2^2}+\frac{x_3^2}{a_3^2}\le1\rbrace\). From \(u=\Lambda x\) we deduce the equation \(\dot{a}_i=a_i\lambda_i\), which, in combination with the initial condition — rest, that is, \(\lambda_i(0)=0\) — it allows us to deduce that \begin{equation} \int_0^t\lambda_i(s)ds=\ln\left(\frac{a_i(0)}{a_i(t)}\right)~~.\label{eq-ai}\end{equation} To start, we will analyze ellipsods with all different semiaxes, which is the most general case. Without losing generality, we will assume hereafter that \( a_1(0) > a_2(0) > a_3(0) \). The equation of continuity \eqref{eq-contE} implies that \begin{equation} \rho(t)=\rho_0\exp\left(\int_0^t\lambda_1(s)+\lambda_2(s)+\lambda_3(s)ds\right)~~. \end{equation}

The momentum equations, on the other hand, are given by $$ -\dot{\lambda}_i(t)+(\lambda_i(t))^2=-\partial_i\phi/x_i ~~,\text{(no sum in }i)~~. $$ The expression for the potential due to an uniform ellipsoid is given by[1] \begin{equation} \phi=\pi G \rho x^\intercal \begin{pmatrix} A_1(t)&& \\ & A_2(t) & \\ & & A_3(t)\\ \end{pmatrix} x\quad\rightarrow\quad \partial_i\phi/x_i=2\pi G\rho A_i(t)~~.\label{eq-potE} \end{equation} Expressions for the \(A_i\) are given in terms of the following function $$ \mathcal{A}(\alpha,\beta)=\int_0^\infty\frac{d\eta}{(1+\eta)^{3/2}\sqrt{(1+\alpha^2 \eta)(1+\beta^2\eta)}}~~,$$ for \(\alpha,\beta > 0\). Clearly this function is symmetric in its arguments. Then \begin{equation} A_1(t)=\mathcal{A}\left(\frac{a_1}{a_2},\frac{a_1}{a_3}\right)\quad A_2(t)=\mathcal{A}\left(\frac{a_2}{a_1},\frac{a_2}{a_3}\right)\quad A_3(t)=\mathcal{A}\left(\frac{a_3}{a_1},\frac{a_3}{a_2}\right)~~. \end{equation} There are elementary expressions for the \(A_i\) in the case of spheroids, that is, for a pair of equal semiaxes. In the spherical case, \(A_i=\mathcal{A}(1,1)=2/3\) for each \(i\). An important general property is[1] \begin{equation}A_1+A_2+A_3=2~~.\label{eq-propA}\end{equation}

We conclude then that the momentum equations are given by \begin{align} \dot{\lambda}_i(t)&=(\lambda_i(t))^2+2\pi G A_i(t)\rho_0\exp\left(\int_0^t \lambda_1(s)+\lambda_2(s)+\lambda_3(s)ds \right)~~, \label{eq-lmom}\\ \lambda_i(0)&=0~~,\label{eq-l0}\\ \end{align} for \(i=1,2,3\). Despite the complications, the solutions to these equations solve all the hydrodynamic equations and describe exactly what we expected: a homogeneous, contracting ellipsoid whose axes remain parallel throughout the entire free-fall collapse. The velocity field remains homotetic respect to the origin during the collapse. The collapse, however, ends in finite time, as it is shown next.

2a.3 The free-fall time of ellipsoids, spheroids, the sphere.

The closed system of equations \eqref{eq-lmom} with initial conditions \eqref{eq-l0} define the evolution of the collapsing ellipsoid. We define \(\beta_i=\int_0^t\lambda_i(s)ds=\ln(a_i(0)/a_i(t))\) in order to get rid of the integral, and we scale the time by defining \(d\xi=\sqrt{2\pi G\rho_0}dt\). The equivalent system is \begin{align} \ddot{\beta_i}(\xi)&=\left(\dot{\beta_i}(\xi)\right)^2+ A_i(\xi) \exp\left(\beta_1(\xi)+\beta_2(\xi)+\beta_3(\xi)\right)~~,\label{eq-bmom}\\ \beta_i(0)&=\dot{\beta_i}(0)=0\label{eq-b0}~~,\\ \end{align} where \(\dot{(\,\,)}\) represent derivative respect to \(\xi\). In explicit terms, \begin{equation} A_1(\xi)=\mathcal{A}\left(e^{\beta_2-\beta_1}\frac{a_1(0)}{a_2(0)},e^{\beta_3-\beta_1}\frac{a_1(0)}{a_3(0)}\right)~~, \end{equation} with analogous expressions for \(A_2\) and \(A_3\). The function \(\mathcal{A}(\alpha,\beta)\) is a symmetric, decreasing function in both parameters. If \(\alpha \le \beta\), then \(\mathcal{A}(\beta,\beta)\le\mathcal{A}(\alpha,\beta)\le\mathcal{A}(\alpha,\alpha)\). An explicit formula when \(\alpha=\beta\) is \begin{equation}\mathcal{A}(\alpha,\alpha)=\frac{2}{\alpha^2-1}\left(\frac{\alpha\,\text{acosh}(\alpha)}{\sqrt{\alpha^2-1}}-1\right)=\frac{2}{\alpha^2+\alpha+1}+\epsilon(\alpha)~~,\end{equation} where \(|\epsilon(\alpha)| < 1/5\) and \(\epsilon(0)=\epsilon(1)=\epsilon(\infty)=0\). An even better approximation can be obtained if the first term in the rightmost side of the previos equation is multiplied by \((1+\ln(\alpha/2+1/2))\) for \(\alpha \ge 1\), with the error being less than 25% and 0.15 in relative and absolute terms, respectively. Despite the explicit appearance of the equation, therere is no discontinuity in \(\alpha=1\). Another region where the function is simple is when \(\beta=0\), in which case \(\mathcal{A}(\alpha,0)=2/(\alpha+1)\).

In order to estimate very roughly the free-falling times of these ellipsoids we note that at the initial time, the largest \(A_i\) is associated with the smallest semi-axis (\(a_3\) in our case). The evolution of the system will take \(\beta_3\) to dominate over \(\beta_1\) and \(\beta_2\), and finally blowing up before them. This type of organized collapse is somewhat described as a progression into a "pancake" and then a "spindle" (and a 0-dimensional singularity, presumably). It is possible to formalize these expectation using differential inequalities, but the derivation is somewhat involved and I will skip it. In [3], Lin et al. give estimations of the free-fall times. Our rough free-fall time estimations start with the following implicit solution of a system related with \eqref{eq-bmom}

The implicit solution given by Equation \eqref{eq-diff} can be obtained first transforming the equation using the substitution[6] $$ \varpi(b)=\left(\frac{db}{dt}(t(b))\right)^2~~,\tag{17.1}$$ which implies that $$ \frac{d}{dt}\left(\frac{db}{dt}\right)=\frac{d}{dt}\left(\sqrt{\varpi}\right)=\frac{1}{2\sqrt{\varpi}}\frac{d\varpi}{db}\frac{db}{dt}=\frac{1}{2\cancel{\sqrt{\varpi}}}\frac{d\varpi}{db}\cancel{\sqrt{\varpi}}=\frac{\varpi'}{2}~~.\tag{17.2}$$ On the other hand, the initial condition is evaluated in \(\varpi(b(0))=\varpi(0)=(\dot{b}(0))^2=0\). Therefore, the system \eqref{eq-diff} becomes a linear, first order system for \(\varpi\) in the new independent variable \(b\), that is, \begin{equation*} \begin{aligned}\varpi'&=2\varpi+2g(b)e^b\\\varpi(0)&=0\end{aligned}~~\rightarrow~~\varpi(b)=2e^{2b}\int_0^be^{-s}g(s)ds~~.\tag{17.3}\end{equation*} From (17.3) we derive the autonomous equation $$ \frac{db}{dt}=\sqrt{2}e^b\left(\int_0^be^{-s}g(s)ds\right)^{1/2}~~,\tag{17.4}$$ from where an implicit solution can be obtained (using also \(b(0)=0\)) $$ \int_0^b\frac{dx}{\sqrt{2}e^x\left(\int_0^xe^{-s}g(s)ds\right)^{1/2}}=t~~,\tag{17.5}$$ which is \eqref{eq-diff}.

\(\blacksquare\)

The implicit solution inmediately gives us a criterion for blow-up in finite time (or free-fall time): \(b(t)\) diverges in finite time if and only if the integral above, integrated between \(0\) and \(\infty\), converges to a finite value. The finite value of the integral is the blow-up time.

By replacing \(A_3(\xi)\) by \(A_3(0)\) in \eqref{eq-bmom}, we obtain the following estimations of the free-fall \(\xi\) \begin{equation} \xi_\text{ff,3} \begin{cases} \le\sqrt{\frac{2}{A_3(0)}}& a_3 < a_2,a_1, \text{(includes oblate case); } \beta_{1,2}\approx0& (2>A_3(0)>\frac{2}{3})\\ \le\sqrt{\frac{\pi}{2A_3(0)}}& a_3 = a_2 < a_1, \text{(prolate case); } \beta_3=\beta_2\text{ and }\beta_1\approx0 & (1>A_3(0)>\frac{2}{3})\\ =\sqrt{\frac{3}{4}}\frac{\pi}{2} & \text{sphere; } \beta_{i}=\beta, ~i=1,2,3& (A_i(t)=\frac{2}{3})\\ \end{cases}~~,\label{eq-xiff} \end{equation} from where we obtain time using \(\xi=\sqrt{2\pi G \rho_0}t\). This justifies the well-known result\(t_\text{ff}=\sqrt{\frac{3\pi}{32 G \rho_0}}\) in the case of spherical collapse. Because in general \(\beta_{1,2}\le\beta_3\), we have that the collapse time is bounded by the spherical case from below. In a way, the spherical case is the fastest way of collapse. Note that we are following the ellipsoid evolution only up to the "first" collapse, that is, the one that occurs along the smallest semi-axis.