Editing Aps/MaclaurinSpheroidFreeFall

__FORCETOC__
<!-- __NOTOC__ will force TOC off -->
=Free-Fall Collapse of an Homogeneous Spheroid=
{| class="FreeFallSpheroid" style="float:left; margin-right: 20px; border-style: solid; border-width: 3px border-color: black"
|- 
! style="height: 125px; width: 125px; background-color:#ffff99;" |[[H_BookTiledMenu#Nonlinear_Dynamical_Evolution_2|<b>Free-Fall<br />Collapse<br />of an<br />Homogeneous<br />Spheroid</b>]]
|}
<table border="0" cellpadding="3" align="center" width="80%">
<tr><td align="left">
<font color="darkgreen">
"What is the form of the collapse under gravitational forces of a uniformly rotating spheroidal gas cloud?  In the special case where initially the gas is absolutely cold and of uniform density within the spheroid, we show that the collapse proceeds through a series of uniform, uniformly rotating spheroids until a disk is formed."
</font>
</td></tr>
<tr><td align="right">
&#8212; Drawn from the first paragraph of {{ LB62 }}
<!--[https://www-cambridge-org.libezp.lib.lsu.edu/core/journals/mathematical-proceedings-of-the-cambridge-philosophical-society/article/on-the-gravitational-collapse-of-a-cold-rotating-gas-cloud/EA358305189B0F305818C79AAEB0709F D. Lynden-Bell (1962)], Math. Proc. Cambridge Phil. Soc., Vol. 58, Issue 4, pp. 709 - 711 
-->
</td></tr></table>

==Simplified Governing Relations==

When studying the dynamical evolution of strictly axisymmetric configurations, it proves useful to write the spatial operators in our overarching set of [[PGE|principal governing equations]] in terms of cylindrical coordinates, <math>(\varpi, \phi, z)</math>, and to simplify the individual equations as described in our [[AxisymmetricConfigurations/PGE#Governing_Equations_.28CYL..29|accompanying discussion]].  The resulting set of simplified governing relations is &hellip;

<div align="center">
<span id="Continuity"><font color="#770000">'''Equation of Continuity'''</font></span><br />

<math>\frac{d\rho}{dt} + \frac{\rho}{\varpi} \frac{\partial}{\partial\varpi} \biggl[ \varpi \dot\varpi \biggr] 
+ \rho \frac{\partial}{\partial z} \biggl[ \rho \dot{z} \biggr] = 0 </math><br />


<span id="PGE:Euler">
<font color="#770000">'''Euler Equation'''</font>
</span><br />

<math>
{\hat{e}}_\varpi \biggl[ \frac{d \dot\varpi}{dt} -  \frac{j^2}{\varpi^3} \biggr]  + {\hat{e}}_z \biggl[ \frac{d \dot{z}}{dt} \biggr] = -
{\hat{e}}_\varpi \biggl[ \frac{1}{\rho}\frac{\partial P}{\partial\varpi} + \frac{\partial \Phi}{\partial\varpi}\biggr] -  {\hat{e}}_z \biggl[ \frac{1}{\rho}\frac{\partial P}{\partial z} + \frac{\partial \Phi}{\partial z} \biggr] 
</math><br />
where, the specific angular momentum, <math>j(\varpi,z) \equiv \varpi^2 \dot\phi =  \mathrm{constant} ~(\mathrm{i.e.,}~\mathrm{independent~of~time})
</math><br />



<span id="PGE:AdiabaticFirstLaw">Adiabatic Form of the<br />
<font color="#770000">'''First Law of Thermodynamics'''</font></span><br />

{{Math/EQ_FirstLaw02}}


<span id="PGE:Poisson"><font color="#770000">'''Poisson Equation'''</font></span><br />

<math>
\frac{1}{\varpi} \frac{\partial }{\partial\varpi} \biggl[ \varpi \frac{\partial \Phi}{\partial\varpi} \biggr] + \frac{\partial^2 \Phi}{\partial z^2} = 4\pi G \rho . 
</math><br />
</div>

This study is closely tied to our [[SSC/Dynamics/FreeFall#Free-Fall|separate discussion of the free-fall collapse of uniform-density ''spheres'']].  For example, by definition, an element of fluid is in "free fall" if its motion in a gravitational field is unimpeded by pressure gradients.  The most straightforward way to illustrate how such a system evolves is to set <math>P = 0</math> in all of the governing equations.  In doing this, the continuity equation and the Poisson equation remain unchanged; the equation formulated by the first law of thermodynamics becomes irrelevant; and the two components of the Euler equation become,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right"><math>\hat{\mathbf{e}}_\varpi</math>:</td>
  <td align="right">
<math>\frac{d\dot{\varpi}}{dt} - \frac{j^2}{\varpi^3}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>- \frac{\partial\Phi}{\partial\varpi}    \, ,</math>
  </td>
</tr>

<tr>
  <td align="right"><math>\hat{\mathbf{e}}_z</math>:</td>
  <td align="right">
<math>\frac{d\dot{z}}{dt} </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>- \frac{\partial\Phi}{\partial z}    \, .</math>
  </td>
</tr>
</table>

==Gravitational Potential==
Here, our specific interest is in modeling the free-fall collapse of a uniform-density spheroid. Ignoring, for the moment, the time-dependent nature of this problem, we appreciate from [[ThreeDimensionalConfigurations/HomogeneousEllipsoids#Gravitational_Potential|a separate, detailed derivation]] that the gravitational potential inside (or on the surface) of an homogeneous, triaxial ellipsoid with semi-axes <math>(x, y, z) = (a_1, a_2, a_3)</math> is given, to within an arbitrary additive constant, by the expression,

<div align="center">
<math>
~\Phi(\vec{x}) = \pi G \rho \biggl[ A_1 x^2 + A_2 y^2 +A_3 z^2 \biggr],
</math>
</div>
where the three, spatially independent coefficients, <math>A_1, A_2,</math> and <math>A_3</math> are functions of the chosen lengths of the three semi-axes.  When deriving mathematical expressions for each of the three <math>A_i</math> coefficients, in our [[ThreeDimensionalConfigurations/HomogeneousEllipsoids#Derivation_of_Expressions_for_Ai|accompanying discussion]] we have found it useful to initially attach a subscript, <math>(\ell, m, ~\mathrm{or}~ s)</math>  &#8212; indicating whether the coefficient is associated with the (largest, medium-length, or smallest) semi-axis &#8212; before specifying how, for a given problem, <math>(\ell, m, s)</math> are appropriately associated with the three <math>(x, y, z)</math> coordinate axes.  

===Oblate Spheroids===
For example, for an oblate-spheroidal mass distribution, by definition the "largest" and "medium-length" semi-axes are equal to one another.  Hence, <math>a_\ell = a_m</math> and, according to [[ThreeDimensionalConfigurations/HomogeneousEllipsoids#When_am_.3D_a.E2.84.93|our associated derivation]],

<table border="0" cellpadding="5" align="center">
<tr>
  <td align="right">
<math>A_\ell = A_m</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>-
\frac{(1-e^2)}{e^2} + 
~\frac{(1-e^2)^{1 / 2}}{ e^3 } \biggl[
\sin^{-1}e
\biggr] \, ,
</math>
  </td>
</tr>
</table>

where, <math>e \equiv (1 - a_s^2/a_\ell^2)^{1 / 2}</math>; and,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>A_s</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{2}{e^2} - 
~\frac{2(1-e^2)^{1 / 2}}{ e^3 } \biggl[
\sin^{-1}e
\biggr] \, .
</math>
  </td>
</tr>
</table>
Conventionally, the <math>z</math>-axis is aligned with the symmetry (in this case, shortest) axis of the mass distribution, so we set <math>A_1 = A_2 = A_\ell</math>, and <math>A_3 = A_s</math>.  Therefore &#8212; see also our [[ThreeDimensionalConfigurations/HomogeneousEllipsoids#Oblate_Spheroids_.7F.27.22.60UNIQ--postMath-00000039-QINU.60.22.27.7F|parallel discussion]] &#8212; we appreciate that, for oblate-spheroidal mass distributions,

<div align="center">
<math>
\Phi(\vec{x}) = \pi G \rho \biggl[A_\ell \varpi^2 + A_s z^2 \biggr].
</math>
</div>
These same coefficient expressions may also be found in, for example: &nbsp;Chapter 3, Eq. (36) of [<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>]; &sect;4.5, Eqs. (48) &amp; (49) of [<b>[[Appendix/References#T78|<font color="red">T78</font>]]</b>]; and the first column of Table 2-1 (p. 57) of [<b>[[Appendix/References#BT87|<font color="red">BT87</font>]]</b>].  Note that, as we have pointed out in a [[Apps/MaclaurinSpheroids#Gravitational_Potential|separate discussion of Maclaurin Spheroids]], the expressions for <math>A_\ell</math> and <math>A_s</math> have the following values in the limit of a sphere <math>(e=0)</math> or in the limit of an infinitesimally thin disk <math>(e=1)</math>:


<table align="center" border=1 cellpadding="8">
<tr>
  <td colspan="3" align="center">
'''Table 1:''' &nbsp;[[Appendix/Ramblings/PowerSeriesExpressions#Maclaurin_Spheroid_Index_Symbols|Limiting Values]] (for oblate spheroids)
  </td>
</tr>
<tr>
  <td align="center">
&nbsp;
  </td>
  <td align="center">
<b><math>e \rightarrow 0</math></b>
  </td>
  <td align="center">
<b><math>\frac{a_3}{a_1} \rightarrow 0</math></b>
  </td>
</tr>
<tr>
  <td align="center">
<b><math>~\frac{\sin^{-1}e}{e}</math></b>
  </td>
  <td align="center">
<math>~1 + \frac{e^2}{6} + \mathcal{O}\biggl(e^4\biggr)</math>
  </td>
  <td align="center">
<math>~\frac{\pi}{2} - \biggl(\frac{a_3}{a_1}\biggr) +\frac{\pi}{4}\biggl(\frac{a_3}{a_1}\biggr)^2
- \mathcal{O}\biggl(\frac{a_3^3}{a_1^3}\biggr)</math>
  </td>
</tr>
<tr>
  <td align="center">
<b><math>A_\ell</math></b>
  </td>
  <td align="center">
<math>\frac{2}{3}\biggl[1 - \frac{e^2}{5} - \mathcal{O}\biggl(e^4\biggr)\biggr]</math>
  </td>
  <td align="center">
<math>\frac{\pi}{2} \biggl( \frac{a_3}{a_1}\biggr) - 2\biggl(\frac{a_3}{a_1}\biggr)^2+ \mathcal{O}\biggl(\frac{a_3^3}{a_1^3}\biggr)</math>
  </td>
</tr>
<tr>
  <td align="center">
<b><math>A_s</math></b>
  </td>
  <td align="center">
<math>\frac{2}{3}\biggl[1 + \frac{2e^2}{5} + \mathcal{O}\biggl(e^4\biggr)\biggr]</math>
  </td>
  <td align="center">
<math> 2 - \pi \biggl( \frac{a_3}{a_1}\biggr) + 4\biggl(\frac{a_3}{a_1}\biggr)^2
- \mathcal{O}\biggl(\frac{a_3^3}{a_1^3}\biggr)</math>
  </td>
</tr>
</table>

===Prolate Spheroids===
For a prolate-spheroidal mass distribution, by definition the "smallest" and "medium-length" semi-axes are equal to one another.  Hence, <math>a_m = a_s</math> and, according to [[ThreeDimensionalConfigurations/HomogeneousEllipsoids#When_am_.3D_as|our associated derivation]],

<table border="0" cellpadding="5" align="center">
<tr>
  <td align="right">
<math>A_\ell</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{ (1-e^2)}{e^3} \cdot 
\ln \biggl[ \frac{1+e}{ 1-e }   \biggr]
- \frac{2(1-e^2)}{e^2 } 
\, ,
</math>
  </td>
</tr>
</table>

where, as above, <math>e \equiv (1 - a_s^2/a_\ell^2)^{1 / 2}</math>; and,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>A_s = A_m</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{1}{e^2 } 
- \frac{ (1-e^2)}{2e^3} 
\cdot
\ln \biggl( \frac{1+e}{ 1-e }   \biggr) \, .
</math>
  </td>
</tr>
</table>
If the symmetry (in this case, longest) axis of this prolate-spheroidal mass distribution is aligned with the <math>z</math>-axis of the coordinate system, then we should set <math>A_1 = A_2 = A_s</math>, and <math>A_3 = A_\ell</math>.  This means that the expression for the gravitational potential is,

<div align="center">
<math>
\Phi(\vec{x}) = \pi G \rho \biggl[A_s \varpi^2 + A_\ell z^2 \biggr].
</math>
</div>
These same coefficient expressions may also be found in, for example, the second column of Table 2-1 (p. 57) of [<b>[[Appendix/References#BT87|<font color="red">BT87</font>]]</b>].


<font color="red">NOTE:</font>&nbsp; If, following [<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>], we instead align the longest (and, in this case, symmetry) axis of the prolate mass distribution with the <math>x</math>-axis, then <math>A_2 = A_3 = A_s</math> and <math>A_1 = A_\ell</math>.  This matches the coefficient expressions presented in our [[ThreeDimensionalConfigurations/HomogeneousEllipsoids#Prolate_Spheroids_.7F.27.22.60UNIQ--postMath-0000004C-QINU.60.22.27.7F|parallel discussion]] of the potential inside and on the surface of a prolate-spheroidal mass distribution.

===Consider a Time-Varying Eccentricity===

If the eccentricity of an homogeneous spheroid varies with time &#8212; that is, if <math>e \rightarrow e(t)</math> &#8212; while it remains homogeneous, the result will be a potential of the form,

<div align="center">
<math>
\Phi(\vec{x}, t) = A(t) \varpi^2 + C(t) z^2 \, ,
</math>
</div>

whether the spheroid is oblate or the spheroid is prolate.

==Lynden-Bell's (1962) Insight==
Let's examine the analysis by {{ LB62full }} &#8212; hereafter, {{ LB62hereafter }} &#8212; of the "Gravitational Collapse of a Cold Rotating Gas Cloud."

===Motion of a Single Particle===
Consider a particle that, at time <math>t=0</math>, is at position <math>(\varpi_0, \phi_0, z_0)</math> and is moving about the <math>z</math>-axis with velocity, <math>\varpi_0\Omega ~\Rightarrow~ j_0 = \varpi_0^2 \Omega</math>.  Consider furthermore that its acceleration is subject to the force arising from an axisymmetric gravitational potential of the form,

<div align="center">
<math>
\Phi(\vec{x}, t) = A(t) \varpi^2 + C(t) z^2 .
</math>
</div>

[This is the gravitational potential adopted by {{ LB62hereafter }} &#8212; see his equation (1) &#8212; except he adopted a different sign convention to ours.  He would therefore have also attached a sign to the gradient of the potential that is the opposite of the sign that appears on the right-hand side of our Euler equation expression.]  In this case, the two components of the Euler equation that govern the particle's motion are,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right"><math>\hat{\mathbf{e}}_\varpi</math>:</td>
  <td align="right">
<math>\ddot\varpi - \frac{j_0^2}{\varpi^3}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>- 2A\varpi    \, ,</math>
  </td>
</tr>

<tr>
  <td align="right"><math>\hat{\mathbf{e}}_z</math>:</td>
  <td align="right">
<math>\ddot{z} </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>- 2Cz    \, ,</math>
  </td>
</tr>
</table>
where we have adopted the familiar shorthand notation, <math>d\dot\varpi/dt \rightarrow \ddot\varpi</math> and <math>d\dot{z}/dt \rightarrow \ddot{z}</math>.  If we divide the first of these relations by <math>\varpi_0</math> and the second by <math>z_0</math>, then adopt the dimensionless variables, <math>R \equiv \varpi/\varpi_0</math> and <math>Z \equiv z/z_0</math>, we can write,
<table border="0" cellpadding="3" align="center">

<tr>
  <td align="right">
<math>\frac{\ddot\varpi}{\varpi_0}  - \frac{j_0^2}{\varpi_0 \varpi^3}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align=left">
<math>- 2A \biggl( \frac{\varpi}{\varpi_0} \biggr)</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow ~~~ \ddot{R}  - \frac{\Omega^2}{R^3}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align=left">
<math>- 2A R \, ,</math>
  </td>
</tr>
<tr>
  <td align="center" colspan="3">
{{ LB62hereafter }}, p. 710, Eq. (10)
  </td>
</tr>
</table>
and,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>\frac{\ddot{z}}{z_0} </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>- 2C \biggl(\frac{z}{z_0}\biggr)</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow ~~~ \ddot{Z} </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>- 2C Z \, .</math>
  </td>
</tr>
<tr>
  <td align="center" colspan="3">
{{ LB62hereafter }}, p. 710, Eq. (11)
  </td>
</tr>
</table>
Finally, we use <math>\varphi(t)</math> to represent the particle's time-varying angular-coordinate position ''relative to'' its initial position &#8212; that is, we adopt the definition, <math>\varphi(t) \equiv \phi(t) - \phi_0</math>. Then, conservation of angular momentum implies that, at any moment, the particle's rotation frequency about the symmetry axis will be,
<table border="0" cellpadding="0" align="center">

<tr>
  <td align="right">
<math>\dot\varphi = \dot\phi</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>\frac{\Omega}{R^2} \, .</math>
  </td>
</tr>
<tr>
  <td align="center" colspan="3">
{{ LB62hereafter }}, p. 710, Eq. (9)
  </td>
</tr>
</table>
This governing set of evolutionary equations has been set up such that at time, <math>t=0</math>: <math>R = 1</math>, <math>Z = 1</math>, <math>\varphi = 0</math>, <math>\dot{R} = 0</math>, <math>\dot{Z} = 0</math>, and <math>\dot\varphi = \Omega</math>. With this set of initial conditions in hand, along with an appropriate specification of the two time-dependent coefficients, <math>A(t)</math> and <math>C(t)</math>, the set of governing relations can be integrated (numerically) to give <math>R(t), \varphi(t)</math>, and <math>Z(t)</math>.  This is the result that {{ LB62 }} established for the motion of one particle.

===Evolution of the Spheroid===
Following his examination of the motion of an individual particle, {{ LB62hereafter }} recognized that, <font color="darkgreen">"<math>R(t), \varphi(t)</math>, and <math>Z(t)</math> are all independent of <math>\varpi_0, \phi_0</math>, and <math>z_0</math> because</font> [none of the three evolutionary equations] <font color="darkgreen">nor the above initial conditions mention them."</font>  You only need to integrate the coupled set of governing relations once then &#8212; assuming that the functions, <math>A(t)</math> and <math>C(t)</math>, are the same in all cases &#8212; the time-dependent coordinates of any particle are given by <math>(\varpi_0 R, \phi_0 + \varphi, z_0 Z)</math>, where <math>(\varpi_0 , \phi_0, z_0)</math> are the initial coordinates of that particle.  <font color="darkgreen">"Thus the result of the motion is merely a change of scales."</font>

Consider then, as did {{ LB62hereafter }}, the evolution of a spheroid that is initially uniformly filled with free particles and whose only motion, initially, is uniform rotation, <math>\Omega</math>, about the z-axis.  As {{ LB62hereafter }} puts it, since the motion of each particle can be described merely via a change of scales: <font color="darkgreen">"&hellip; the distribution of the particles remains uniform, and the boundary remains spheroidal"</font>; and, while the angular frequency of each particle, <math>\dot\phi</math>, varies with time, <font color="darkgreen">"&hellip; since <math>\dot\phi = \dot\varphi</math> the rotation remains uniform in space."</font>

It should be clear, as well, that the eccentricity of the evolving spheroid will vary with time.  Specifically in the case of an oblate spheroid, the time-dependent semi-axes are <math>(a_\ell R(t),a_\ell R(t), a_s Z(t))</math>; hence,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>e(t)</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl[ 1 - \biggl( \frac{a_s Z}{a_\ell R} \biggr)^2 \biggr]^{1 / 2}
=
\biggl[ 1 - ( 1 - e_0^2 )\biggl( \frac{Z}{R} \biggr)^2 \biggr]^{1 / 2} \, ,
</math>
  </td>
</tr>
</table>
where, <math>e_0 = (1 - a_s^2/a_\ell^2)^{1 / 2}</math> is the eccentricity of the spheroid initially, and the time-variation enters via the pair of functions, <math>Z(t)</math> and <math>R(t)</math>.  In the case of an prolate spheroid, the time-dependent semi-axes are <math>(a_s R(t),a_s R(t), a_\ell Z(t))</math>; hence,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>e(t)</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl[ 1 - \biggl( \frac{a_s R}{a_\ell Z} \biggr)^2 \biggr]^{1 / 2}
=
\biggl[ 1 - ( 1 - e_0^2 )\biggl( \frac{R}{Z} \biggr)^2 \biggr]^{1 / 2} \, .
</math>
  </td>
</tr>
</table>
<span id="Table1">In turn,</span> the time-dependent behavior of the coefficients in the expression for the gravitational potential, <math>A(t)</math> and <math>C(t)</math>, is drawn from <math>e(t)</math> as detailed in Table 1, immediately below.

<table border="1" align="center" cellpadding="3">
<tr>
  <td align="center" colspan="3">
Table 1: &nbsp; Time-Dependent Coefficients of the Gravitational Potential<br />
<div align="center"><math>\Phi(\vec{x}, t) = A(t) \varpi^2 + C(t) z^2\, ,</math></div>
where, it is understood that the eccentricity of the spheroid, <math>e(t)</math>, varies with time.
  </td>
</tr>
<tr>
  <td align="center">Oblate Spheroid</td>
  <td align="center">Prolate Spheroid</td>
</tr>
<tr>
  <td align="center">
<table border="0" cellpadding="5" align="center">
<tr>
  <td align="right">
<math>\frac{A(t)}{\pi G \rho(t)}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>-
\frac{(1-e^2)}{e^2} + 
~\frac{(1-e^2)^{1 / 2}}{ e^3 } \biggl[
\sin^{-1}e
\biggr] \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\frac{C(t)}{\pi G \rho(t)}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{2}{e^2} - 
~\frac{2(1-e^2)^{1 / 2}}{ e^3 } \biggl[
\sin^{-1}e
\biggr] \, ,
</math>
  </td>
</tr>
<tr>
  <td align="left" colspan="3">where,</td>
</tr>

<tr>
  <td align="right">
<math>e(t)</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl\{ 1 - ( 1 - e_0^2 )\biggl[ \frac{Z(t)}{R(t)} \biggr]^2 \biggr\}^{1 / 2} \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\frac{\rho(t)}{\rho_0}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{1}{R^2(t)Z(t)} \, .
</math>
  </td>
</tr>
</table>
  </td>
  <td align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>\frac{A(t)}{\pi G \rho(t)}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{1}{e^2 } 
- \frac{ (1-e^2)}{2e^3} 
\cdot
\ln \biggl( \frac{1+e}{ 1-e }   \biggr) \, ,
</math>
  </td>
</tr>
<tr>
  <td align="right">
<math>\frac{C(t)}{\pi G \rho(t)}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{ (1-e^2)}{e^3} \cdot 
\ln \biggl[ \frac{1+e}{ 1-e }   \biggr]
- \frac{2(1-e^2)}{e^2 } 
\, ,
</math>
  </td>
</tr>
<tr>
  <td align="left" colspan="3">where,</td>
</tr>

<tr>
  <td align="right">
<math>e(t)</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl\{ 1 - ( 1 - e_0^2 )\biggl[ \frac{R(t)}{Z(t)} \biggr]^2 \biggr\}^{1 / 2} \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\frac{\rho(t)}{\rho_0}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{1}{R^2(t)Z(t)} \, .
</math>
  </td>
</tr>
</table>

  </td>
</tr>
</table>

==Our Numerical Integration==

===Case of the Oblate Spheroid===
====Governing, Dimensionless Differential Equations====
We are primarily interested in determining, for various initial values of the parameter pair <math>(e_0, \Omega)</math>, how rapidly an oblate spheroid collapses to an infinitesimally thin disk <math>(e \rightarrow 1)</math>, and what the radius of this disk is at the instant it forms.  As has been detailed above, according to {{ LB62hereafter }} the time-dependent behavior of <math>R</math> and <math>Z</math> &#8212; and, hence, of <math>e</math> &#8212; can be obtained by integrating the following pair of governing relations:
<table border="0" cellpadding="3" align="center">

<tr>
  <td align="right">
<math>\ddot{R}  - \frac{\Omega^2}{R^3}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align=left">
<math>- 2A R \, ,</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\ddot{Z} </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>- 2C Z \, .</math>
  </td>
</tr>
</table>

Turning both of these 2<sup>nd</sup>-order ODEs into a pair of 1<sup>st</sup>-order ODEs, then multiplying through by <math>\tau_\mathrm{ff} \equiv [3\pi/(32G\rho_0)]^{1 / 2}</math> or, as appropriate, by <math>\tau_\mathrm{ff}^2</math>, we have,
<table border="0" cellpadding="3" align="center">

<tr>
  <td align="right">
<math>\tau_\mathrm{ff} \cdot \frac{d(\tau_\mathrm{ff}\dot{R})}{dt}  </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align=left">
<math>\frac{3\pi}{32G\rho_0} \biggl[ \frac{\Omega^2}{R^3}- 2A R \biggr] </math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align=left">
<math>\frac{3\pi^2}{16} \biggl\{ 
\frac{\Omega^2/(2\pi G \rho_0)}{R^3} - \biggl[\frac{\rho(t)}{\rho_0}  \biggr]\frac{A(t)}{\pi G \rho(t)} \cdot R 
\biggr\} \, ,</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\tau_\mathrm{ff} \cdot \frac{dR}{dt}  </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align=left">
<math>\tau_\mathrm{ff}\dot{R} \, ,</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\tau_\mathrm{ff} \cdot \frac{d(\tau_\mathrm{ff}\dot{Z})}{dt}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
- \frac{3\pi}{32G\rho_0} \biggl[ 2C Z \biggr]
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
- \frac{3\pi^2}{16} \biggl\{ \biggl[\frac{\rho(t)}{\rho_0}  \biggr] \frac{C(t) }{\pi G \rho(t)}\cdot Z \biggr\} 
\, ,</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\tau_\mathrm{ff} \cdot \frac{dZ}{dt}  </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align=left">
<math>\tau_\mathrm{ff}\dot{Z} \, .</math>
  </td>
</tr>
</table>
Adopting the shorthand notation,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>\sigma_0^2</math>
  </td>
  <td align="center">
<math>\equiv</math>
  </td>
  <td align="left">
<math>\frac{\Omega^2}{2\pi G\rho_0}   \, ,</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\tau</math>
  </td>
  <td align="center">
<math>\equiv</math>
  </td>
  <td align="left">
<math>\frac{t}{\tau_\mathrm{ff}}   \, ,</math>
  </td>
</tr>
</table>
and pulling [[#Table1|expressions for the oblate spheroid from Table 1]] above, we have,
<table border="0" cellpadding="3" align="center">

<tr>
  <td align="right">
<math>\frac{d(\tau_\mathrm{ff}\dot{R})}{d\tau}  </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align=left">
<math>\frac{3\pi^2}{16} \biggl\{ 
\frac{\sigma_0^2}{R^3} - \frac{1}{R Z} \biggl[ 
-\frac{(1-e^2)}{e^2} + 
~\frac{(1-e^2)^{1 / 2}}{ e^3 } \cdot
\sin^{-1}e
\biggr] \biggr\} \, ,</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\frac{dR}{d\tau}  </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align=left">
<math>\tau_\mathrm{ff}\dot{R} \, ,</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\frac{d(\tau_\mathrm{ff}\dot{Z})}{d\tau}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
- \frac{3\pi^2}{16} \biggl\{ \frac{1}{R^2} \biggl[ 
\frac{2}{e^2} - 
~\frac{2(1-e^2)^{1 / 2}}{ e^3 } \cdot \sin^{-1}e \biggr] \biggr\} 
\, ,</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\frac{dZ}{d\tau}  </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align=left">
<math>\tau_\mathrm{ff}\dot{Z} \, .</math>
  </td>
</tr>
</table>

====Discrete Representations====
Along a uniformly segmented, discrete time grid &#8212; time step, <math>\Delta = (\tau_{n+1}-\tau_{n})</math> &#8212; let's define the following discretized variables:
<ul>
  <li>At integer values of <math>\Delta</math>: &nbsp; <math>R</math> and <math>Z</math></li>
  <li>At half-integer values of <math>\Delta</math>: &nbsp; <math>F \equiv (\tau_\mathrm{ff} \dot{R}) \, ; K \equiv (\tau_\mathrm{ff} \dot{Z})</math></li>
</ul>
The four finite-difference relations are, then:
<table border="0" cellpadding="3" align="center">

<tr>
  <td align="right">
<math>\frac{F_{n+1/2} - F_{n-1/2}}{\Delta}  </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align=left">
<math>\frac{3\pi^2}{16} \biggl\{ 
\frac{\sigma_0^2}{R^3} - \frac{1}{R Z} \biggl[ 
-\frac{(1-e^2)}{e^2} + 
~\frac{(1-e^2)^{1 / 2}}{ e^3 } \cdot
\sin^{-1}e
\biggr] \biggr\}_n \, ,</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\frac{R_{n+1} - R_n}{\Delta}  </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align=left">
<math>F_{n+1/2} \, ,</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\frac{K_{n+1/2} - K_{n-1/2}}{\Delta}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
- \frac{3\pi^2}{16} \biggl\{ \frac{1}{R^2} \biggl[ 
\frac{2}{e^2} - 
~\frac{2(1-e^2)^{1 / 2}}{ e^3 } \cdot \sin^{-1}e \biggr] \biggr\}_n 
\, ,</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\frac{Z_{n+1} - Z_{n}}{\Delta}  </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align=left">
<math>K_{n+1/2} \, .</math>
  </td>
</tr>
</table>

<font color="red">STEP #0:</font> 
<ul>
  <li>Choose time-independent values of the parameters, <math>(c_0/a_0)</math> and <math>\beta_0</math>, which means that, <math>\sigma_0^2 = 3\pi^2\beta_0/8</math> and <math>e_0 = [1 - (c_0/a_0)^2]^{1 / 2}</math>.  Also set the (uniform) integration time step, <math>\Delta</math>; for 201 time steps, for example, set <math>\Delta \approx 0.005</math>.<br />
The example values shown below (inside parentheses) assume that <math>(c_0/a_0, \beta_0) = (0.90, 0.0) ~~\Rightarrow (e_0, \sigma_0^2) = (0.43588989, 0)</math>, which corresponds to one of the model-evolutions presented by {{ LMS65 }}.</li>
</ul>

<font color="red">STEP #1:</font> Initially, that is, at integration time step, <math>n = 0</math>
<ul>
  <li>Set <math>R_0 = 1</math>, <math>Z_0=1</math>, and <math>e = e_0</math>; this means that the RHS of the first and third discrete evolution equations is known.</li>
  <li>Given that the configuration is collapsing ''from rest'', we want to set <math>\dot{R}_0</math> and <math>\dot{Z}_0</math> to zero. This is accomplished by establishing reflection symmetry through the (time) origin, that is, by setting <math>F_{+1/2}= - F_{-1/2}</math> and <math>K_{+1/2}= - K_{-1/2}</math>.  Initially, then, the LHS of the first and third discrete equations are, respectively, <math>2F_{+1/2}/\Delta</math> and <math>2K_{+1/2}/\Delta</math>.</li>
  <li>Use the first and third discrete equations to solve for the "R" and "Z" velocities at time step <math>n=0</math>, namely,</li>
</ul>
<table border="0" width="80%" cellpadding="3" align="center">

<tr>
  <td align="right">
<math>F_{+1/2}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align=left">
<math>\frac{\Delta}{2} \cdot \frac{3\pi^2}{16} \biggl\{ 
\sigma_0^2 - \biggl[ 
-\frac{(1-e_0^2)}{e_0^2} + 
~\frac{(1-e_0^2)^{1 / 2}}{ e_0^3 } \cdot
\sin^{-1}e_0
\biggr] \biggr\} \, ,</math>
  </td>
  <td align="right" width="25%">&nbsp;&nbsp;&nbsp; <math>(-2.952452 \times 10^{-3})</math></td>
</tr>

<tr>
  <td align="right">
<math>K_{+1/2}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
- \frac{\Delta}{2} \cdot \frac{3\pi^2}{16} \biggl\{ \biggl[ 
\frac{2}{e_0^2} - 
~\frac{2(1-e_0^2)^{1 / 2}}{ e_0^3 } \cdot \sin^{-1}e_0 \biggr] \biggr\} 
\, .</math>
  </td>
  <td align="right" width="25%">&nbsp;&nbsp;&nbsp; <math>(-3.347849\times 10^{-3})</math></td>
</tr>
</table>
<ul>
  <li>From the second and fourth discrete relations, determine the advanced coordinate positions.</li>
</ul>

<table border="0" width="80%" cellpadding="3" align="center">

<tr>
  <td align="right">
<math>R_{+1} </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align=left">
<math>1 + \Delta \cdot F_{+1/2} \, ,</math>
  </td>
  <td align="right" width="25%">&nbsp;&nbsp;&nbsp; <math>(1.0-1.476226 \times 10^{-5})</math></td>
</tr>

<tr>
  <td align="right">
<math>Z_{+1}  </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align=left">
<math>1 + \Delta \cdot K_{+1/2} \, .</math>
  </td>
  <td align="right" width="25%">&nbsp;&nbsp;&nbsp; <math>(1.0-1.673924\times 10^{-5})</math></td>
</tr>
</table>
<ul>
  <li>Determine the eccentricity at this advanced time.</li>
</ul>
<table border="0" width="80%" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>e_{+1}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl[ 1 - ( 1 - e_0^2 )\biggl( \frac{Z_{+1}}{R_{+1}} \biggr)^2 \biggr]^{1 / 2} \, .
</math>
  </td>
  <td align="right" width="25%">&nbsp;&nbsp;&nbsp; <math>(1 - 0.809996797)^{1 / 2} = (0.435893568)</math></td>
</tr>
</table>

<font color="red">STEP #2:</font> &nbsp; Repeat, in sequence for all values of <math>n > 1</math> until <math>Z</math> passes through zero.
<ul>
  <li>Set <math>n \rightarrow n+1</math>.</li>
  <li>Given knowledge of the various variable values at time-step, <math>(n-1/2)</math> and <math>n</math>, use the first and third discrete evolution relations to determine <math>F_{n+1/2}</math> and  <math>K_{n+1/2}</math>; specifically,</li>
</ul>
<table border="0" width="80%" cellpadding="3" align="center">

<tr>
  <td align="right">
<math>F_{n+1/2} </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align=left">
<math>F_{n-1/2} + \Delta \cdot \frac{3\pi^2}{16} \biggl\{ 
\frac{\sigma_0^2}{R^3} - \frac{1}{R Z} \biggl[ 
-\frac{(1-e^2)}{e^2} + 
~\frac{(1-e^2)^{1 / 2}}{ e^3 } \cdot
\sin^{-1}e
\biggr] \biggr\}_n \, ,</math>
  </td>
  <td align="right" width="25%">&nbsp;&nbsp;&nbsp; <math>(F_{+1/2} -5.905085\times 10^{-3})</math></td>
</tr>

<tr>
  <td align="right">
<math>K_{n+1/2} </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>K_{n-1/2} 
- \Delta \cdot \frac{3\pi^2}{16} \biggl\{ \frac{1}{R^2} \biggl[ 
\frac{2}{e^2} - 
~\frac{2(1-e^2)^{1 / 2}}{ e^3 } \cdot \sin^{-1}e \biggr] \biggr\}_n 
\, .</math>
  </td>
  <td align="right" width="25%">&nbsp;&nbsp;&nbsp; <math>(K_{+1/2} - 6.695907\times 10^{-3})</math></td>
</tr>
</table>
<ul>
  <li>Similarly, from the second and fourth discrete relations, determine the advanced coordinate positions.</li>
</ul>
<table border="0" width="80%" cellpadding="3" align="center">

<tr>
  <td align="right">
<math>R_{n+1}  </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align=left">
<math>R_n + \Delta\cdot F_{n+1/2} \, ,</math>
  </td>
  <td align="right" width="25%">&nbsp;&nbsp;&nbsp; <math>(0.99994095)</math></td>
</tr>

<tr>
  <td align="right">
<math>Z_{n+1}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align=left">
<math>Z_n + \Delta \cdot K_{n+1/2} \, .</math>
  </td>
  <td align="right" width="25%">&nbsp;&nbsp;&nbsp; <math>(0.999933042)</math></td>
</tr>
</table>
<ul>
  <li>Determine the eccentricity at this advanced time.</li>
</ul>
<table border="0" width="80%" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>e_{n+1}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl[ 1 - ( 1 - e_0^2 )\biggl( \frac{Z_{n+1}}{R_{n+1}} \biggr)^2 \biggr]^{1 / 2} \, .
</math>
  </td>
  <td align="right" width="25%">&nbsp;&nbsp;&nbsp; <math>(1 - 0.809987188)^{1 / 2} = (0.43590459)</math></td>
</tr>
</table>

====Results====

<table border="1" align="center" cellpadding="3">
<tr>
  <td align="center" colspan="6">
Table 2:  Collapse of a Nonrotating, Pressure-Free Oblate Spheroid
  </td>
</tr>
<tr>
  <td align="center"><math>\frac{c_0}{a_0}</math></td>
  <td align="center"><math>e_0</math></td>
  <td align="center"><math>\tau_c</math></td>
  <td align="center"><math>R_c</math></td>
  <td align="center"><math>\dot{R}_c</math></td>
  <td align="center"><math>\dot{Z}_c</math></td>
</tr>
<tr>
  <td align="center"><math>0.99</math></td>
  <td align="center"><math>0.141</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
</tr>
<tr>
  <td align="center"><math>0.95</math></td>
  <td align="center"><math>0.312</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
</tr>
<tr>
  <td align="center"><math>0.90</math></td>
  <td align="center"><math>0.436</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
</tr>
<tr>
  <td align="center"><math>0.80</math></td>
  <td align="center"><math>0.600</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
</tr>
<tr>
  <td align="center"><math>0.70</math></td>
  <td align="center"><math>0.714</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
</tr>
<tr>
  <td align="center"><math>0.60</math></td>
  <td align="center"><math>0.800</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
</tr>
<tr>
  <td align="center"><math>0.10</math></td>
  <td align="center"><math>0.995</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
  <td align="center"><math>~</math></td>
</tr>
</table>

=Key References=

* {{ LB62full }}:  ''On the gravitational collapse of a cold rotating gas cloud''  <br />NOTE &hellip; according to the [https://ui.adsabs.harvard.edu/abs/1962PCPS...58..709L/abstract new ADS listing], the authors associated with this paper should be, D. Lynden-Bell &amp; C. T. C. Wall ([https://en.wikipedia.org/wiki/C._T._C._Wall Charles Terence Clegg "Terry" Wall]); however, the archived article, itself, lists Lynden-Bell as the sole author while indicating that the paper was simply being ''communicated'' by Wall.

[[File:CommentButton02.png|right|100px|Comment by J. E. Tohline:  In &sect;II of this "1964" article, Lynden-Bell references his 1962 article with an incorrect year (Lynden-Bell 1963); within his list of REFERENCES, the year (1962) is correct, but the journal volume is incorrectly identified as 50 (it should be vol. 58).]]&nbsp;<br />

* [https://ui.adsabs.harvard.edu/abs/1964ApJ...139.1195L/abstract D. Lynden-Bell (1964)], ApJ, 139, 1195 - 1216:  ''On Large-Scale Instabilities during Gravitational Collapse and the Evolution of Shrinking Maclaurin Spheroids''
* Classic paper by C. C. Lin, Leon Mestel, and Frank Shu [https://ui.adsabs.harvard.edu/abs/1965ApJ...142.1431L/abstract (1965, ApJ, 142, 1431 - 1446)] titled, "The Gravitational Collapse of a Uniform Spheroid."

=See Also=


{{ SGFfooter }}