Editing Appendix/Ramblings/ForCohlHoward

__FORCETOC__ 
<!-- __NOTOC__ will force TOC off -->
=Discussion Regarding Stability Studies Performed by Lebovitz=
==Motivation==
These discussions began in late 2021, when [https://www.nist.gov/people/howard-cohl Howard Cohl] asked if I would be interested in working with him on establishing a better understanding of the stability of Riemann S-Type Ellipsoids.  This discussion relates directly to our study of the work by {{ LL96full }}.

We are motivated to pursue this discussion, in part, because our own research group at LSU has previously carried out some nonlinear dynamical simulations that relate to this topic. See &hellip;
<ul>
<li>
[https://www.youtube.com/watch?v=Fi8Fp26uc-g YouTube Simulation] that draws from the following ApJ publication.
</li>
<li>
<b>Relevant Publication:</b> &nbsp; {{ OTM2007full }}, titled, ''Further Evidence for an Elliptical Instability in Rotating Fluid Bars and Ellipsoidal Stars.''
<table border="1" align="center" cellpadding="8" width="90%">
<tr><td align="left" rowspan="2">
<div align="center"><b>Publication Abstract</b></div><br />
Using a three-dimensional nonlinear hydrodynamic code, we examine the dynamical stability of more than 20 self-gravitating, compressible, ellipsoidal fluid configurations that initially have the same velocity structure as Riemann S-type ellipsoids. Our focus is on ''adjoint'' configurations, in which internal fluid motions dominate over the collective spin of the ellipsoidal figure; Dedekind-like configurations are among this group. We find that, although some models are stable and some are moderately unstable, the majority are violently unstable toward the development of m=1, m=3, and higher-order azimuthal distortions that destroy the coherent, m=2 barlike structure of the initial ellipsoidal configuration on a dynamical timescale. The parameter regime over which our models are found to be unstable generally corresponds with the regime over which incompressible Riemann S-type ellipsoids have been found to be susceptible to an elliptical strain instability. We therefore suspect that an elliptical instability is responsible for the destruction of our compressible analogs of Riemann ellipsoids. The existence of the elliptical instability raises concerns regarding the final fate of neutron stars that encounter the secular bar-mode instability and regarding the spectrum of gravitational waves that will be radiated from such systems.
</td>
<td align="center">[[File:OuSimulation2007.png|300px|OTM Simulation (2007)]]<br />
Initial frame of a [https://www.youtube.com/watch?v=Fi8Fp26uc-g YouTube animation] that shows one model's evolution
</td>
</tr>
<tr>
  <td align="center">[[File:OTM2007Fig14.png|250px|{{ OTM2007 }}, Fig. 14]]</td>
</tr>
</table>
</li>
</ul>

==Understanding the Dimensionality of EFE Index Symbols==
Howard put together a Mathematica script intended to provide &#8212; for any specification of the semi-axis length triplet <math>(a, b, c)</math> &#8212; very high-precision, numerical evaluations of any of the index symbols, <math>A_{ijk\ldots}</math> and <math>B_{ijk\ldots}</math> as defined by Eqs. (103 - 104) in &sect;21 of [<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>].  Originally I suggested that, without loss of generality, he should only need to specify the ''pair'' of length ratios, <math>(1, b/a, c/a)</math>.  In response, Howard pointed out that evaluation of all but a few of the lowest-numbered index symbols &#8212; as defined by [<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>] &#8212; ''does'' explicitly depend on specification of (various powers of) the semi-axis length, <math>a</math>.

<font color="red">Joel's response:</font>&nbsp; Howard is correct!  He should leave the explicit dependence of <math>a</math> &#8212; to various powers &#8212; in his Mathematica notebook's determination of all the EFE index symbols.

Instead, what we should expect is that the evaluation of various ''physically relevant'' parameters will produce results that are independent of the semi-axis length, <math>a</math>; these evaluations should involve combining various index symbols in such a way that the dependence on <math>a</math> disappears.  Consider, for example, our [[ThreeDimensionalConfigurations/RiemannStype#Based_on_Virial_Equilibrium|accompanying discussion]] (click to see relevant expressions) of the virial-equilibrium-based determination of the frequency ratio, <math>f \equiv \zeta/\Omega_f</math>, in equilibrium S-Type Riemann Ellipsoids.  Although most of the required index symbols, <math>A_1, A_2, A_3</math> and <math>B_{12}</math>, are dimensionless parameters, the index symbol <math>A_{12}</math> has the unit of inverse-length-squared.  Notice, however, that when <math>A_{12}</math> appears along with any of these other ''dimensionless'' parameters in the definition of <math>f</math>, it is accompanied by an extra "length-squared" factor, such as <math>a^2</math>.  Hence, although I strongly agree that Howard should continue to include various powers of <math>a</math> (etc.) in his Mathematica notebook expressions, I suspect that, without loss of generality, in the end we will always be able to set <math>a=1</math> and only need to specify the ''pair'' of length ratios, <math>(1, b/a, c/a)</math>.

==Evaluation of Index Symbols==
===Three Lowest-Order Expressions===
In our [[ThreeDimensionalConfigurations/HomogeneousEllipsoids#Derivation_of_Expressions_for_Ai|accompanying derivation of expressions]] for the three lowest-order index symbols <math>A_i</math>, we have used subscripts <math>(\ell, m, s)</math> instead of <math>(1, 2, 3)</math> in order to identify which associated semi-axis length is (largest, medium-length, smallest).  We have derived the following expressions:

<table border="1" align="center" cellpadding="8"><tr><td align="left">
<table border="0" cellpadding="5" align="center">
<tr>
  <td align="right">
<math>\frac{A_\ell}{a_\ell a_m a_s}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{2}{a_\ell^3 ~p^2 \sin^3\alpha} 
\biggl[ F(\alpha, p) - E(\alpha, p) \biggr] \, ;
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\frac{A_m}{a_\ell a_m a_s} </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{ 2}{a_\ell^3 }
\biggl[ \frac{
E(\alpha, p) 
-~(1-p^2)
F(\alpha, p)
-~(a_s/a_m)p^2\sin\alpha}{p^2 (1-p^2)\sin^3\alpha}
\biggr] \, ;
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\frac{A_s}{a_\ell a_m a_s}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{ 2}{a_\ell^3 }
\biggl[\frac{
(a_m/a_s) \sin\alpha
- E(\alpha, p)}{ (1-p^2) \sin^3\alpha } 
\biggr] \, .
</math>
  </td>
</tr>

</table>

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

The corresponding expressions that appear in Howard's Mathematica notebook are:

<table border="1" align="center" cellpadding="8"><tr><td align="left">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>A_1</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{2a_2 a_3}{a_1^2 m \sin^3(\phi) } \biggl[
\mathrm{EllipticF}[\phi, m] - \mathrm{EllipticE}[\phi, m] 
\biggr]
\, ;
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>A_2</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{2a_2 a_3}{a_1^2 m(1-m) \sin^3(\phi) } \biggl[
\mathrm{EllipticE}[\phi, m] - \cos^2\theta \cdot \mathrm{EllipticF}[\phi, m] - \frac{a_3}{a_2}\cdot\sin^2\theta \sin\phi
\biggr]
\, ;
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>A_3</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{2a_2 a_3}{a_1^2 (1-m) \sin^3(\phi) } \biggl[
\frac{a_2}{a_3}\cdot \sin(\phi) - \mathrm{EllipticE}[\phi, m] 
\biggr]
\, .
</math>
  </td>
</tr>
</table>
</td></tr></table>
With a little study it should be clear that our derived expressions for <math>A_i</math> precisely match Howard's Mathematica-notebook expressions when <math>\ell = 1</math>, <math>m = 2</math>, and <math>s = 3</math>, that is, in all cases for which <math>a_1 > a_2 > a_3</math>.  But there will be models to consider (for example, in the uppermost region of the so-called [[ThreeDimensionalConfigurations/Stability/RiemannEllipsoids#Background|"horn-shaped" region for S-Type Riemann Ellipsoids]]) for which <math>a_1 > a_3 > a_2</math>, in which case care must be taken in assigning the proper expressions to <math>A_2</math> and <math>A_3</math>.  Similarly note that most of the Riemann models of [[ThreeDimensionalConfigurations/RiemannTypeI#Riemann_Type_1_Ellipsoids|Type I]], II, and III &#8212; see, for example, Figure 16 (p. 161) in Chapter 7 of [<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>] &#8212; have either <math>a_2 > a_1</math> or <math>a_3 > a_1</math>.

===Determination of Higher-Order Expressions===
Howard's Mathematica notebook performs brute-force integrations to evaluate various higher-order index-symbol expressions.  Why doesn't he instead use recurrence relations, which point back to the elliptic-integral-based expressions for <math>A_1, A_2, A_3</math>? Specifically &hellip;

<table border="0" cellpadding="5" align="center">
<tr>
  <td align="center" colspan="3">
<font color="#770000">'''Index-Symbol Recurrence Relations'''</font>
  </td>
</tr>
<tr>
  <td align="right">
<math>B_{ijk\ell\ldots}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> 
A_{jkl\ldots} - a_i^2 A_{ijkl\ldots}
</math>
  </td>
</tr>
<tr>
  <td align="center" colspan="3">
[<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>], §21, p. 54, Eq. (105)
  </td>
</tr>

<tr>
  <td align="right">
<math>a_i^2A_{ikl\ldots} - a_j^2 A_{jkl\ldots}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> 
+ (a_i^2 - a_j^2) B_{ijk\ell\ldots}
</math>
  </td>
</tr>
<tr>
  <td align="center" colspan="3">
[<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>], §21, p. 54, Eq. (106)
  </td>
</tr>

<tr>
  <td align="right">
<math>A_{ikl\ldots} - A_{jkl\ldots}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> 
- (a_i^2 - a_j^2) A_{ijk\ell\ldots}
</math>
  </td>
</tr>
<tr>
  <td align="center" colspan="3">
[<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>], §21, p. 54, Eq. (107)
  </td>
</tr>
</table>

For example, setting <math>i = 1</math> and <math>j = 2</math> in the third of these expressions gives,

<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>A_{1} - A_{2}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> 
- (a_1^2 - a_2^2) A_{12}
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow ~~~ A_{12}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> 
\frac{A_{2} - A_{1}}{(a_1^2 - a_2^2)} 
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow ~~~ a_1^2 A_{12}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> 
\frac{A_{2} - A_{1}}{(1 - a_2^2/a_1^2)} \, ;
</math>
  </td>
</tr>
</table>
and, from the first of the relations,

<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>B_{12}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> 
A_{2} - a_1^2 A_{12}
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> 
A_{2} - a_1^2 \biggl[ \frac{A_{2} - A_{1}}{(a_1^2 - a_2^2)} \biggr]
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> 
\frac{1}{(a_1^2 - a_2^2)} \biggl[(a_1^2 - a_2^2)A_2
- a_1^2 ( A_{2} - A_{1} )\biggr]
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> 
\frac{1}{(a_1^2 - a_2^2)} \biggl[a_1^2 A_1 - a_2^2A_2\biggr]
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> 
\frac{1}{(1 - a_2^2/a_1^2)} \biggl[A_1 - \frac{a_2^2}{a_1^2} \cdot A_2\biggr]\, .
</math>
  </td>
</tr>
</table>


Also, consider using the set of relations labeled "LEMMA 7" on p. 54 of [<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>].

===Example Test Evaluations===

Some of Howard's 20-digit-precision evaluations of various index symbols have been recorded, for comparison with our separate lower-precision evaluations, as follows:
<ul>
  <li>Values of <math>(A_1, A_2, A_3)</math> are recorded for a model with <math>(a_1, a_2, a_3) = (1, 0.9, 0.641)</math> in the [[ThreeDimensionalConfigurations/RiemannStype#TestPart1|table titled, ''TEST (part 1)'', near the top of our chapter on Riemann S-Type ellipsoids]].</li>
  <li>Values of <math>(A_{12}, B_{12})</math> are recorded for a model with <math>(a_1, a_2, a_3) = (1, 0.9, 0.641)</math> in the [[ThreeDimensionalConfigurations/RiemannStype#TestPart2|table titled, ''TEST (part 2)'' in our chapter on Riemann S-Type ellipsoids]].</li>
</ul>

==Figures ''circa'' Year 2000==

Approximately four years after {{ LL96 }} was published, Norman Lebovitz gave a copy of his stability-analysis (FORTRAN) code to Howard Cohl.  Using this code, Howard was able to generate a large set of growth-rate data that essentially allowed him to reproduce Figure 2b from {{ LL96 }}. 

===Image i3.png===
Howard's plot of this data &#8212; his image i3.png &#8212; is shown immediately below; the abscissa is <math>0 \le b/a \le 1</math> and the ordinate is <math>0 \le c/a \le 1</math>.
<table border="1" align="center" cellpadding="8">
<tr>
  <td align="center">Howard's "i3.png" image</td>
  <td align="center" rowspan="3">[[File:HighResSelfAdjointi3.png|400px|High Resolution]]</td>
</tr>
<tr><td align="center">
[[File:I3 FromCohl.png|450px|i3.png]]
</td></tr>
<tr>
<td align="left">
Compare with Figure 2b of {{ LL96full }}
</td>
</tr>
</table>

===Image i5.png===

In an effort to better examine growth-rate trends in the lower-left quadrant of this {{ LL96 }} figure, Howard plotted the same set of stability-analysis data on an axis pair where the abscissa is still <math>0 \le b/a \le 1</math>, but where, for each value of <math>b/a</math>, the ordinate extends from the lower self-adjoint sequence to the upper self-adjoint sequence &#8212; labeled, respectively, <math>x = +1</math> and <math>x = -1</math> in the classic EFE diagram ([[Appendix/References#EFE|<font color="red">EFE</font>]], &sect;49, p. 147, Fig. 15 or, see [[ThreeDimensionalConfigurations/RiemannStype#Fig2|our accompanying discussion]]).  This is displayed immediately below as Howard's "i5.png" image.

<table border="1" align="center" cellpadding="8">
<tr>
  <td align="center">Howard's "i5.png" image</td>
  <td align="center" rowspan="3">[[File:HighResSelfAdjointi5.png|400px|High Resolution]]</td>
</tr>
<tr><td align="center">
[[File:I5 FromCohl.png|450px|i5.png]]
</td></tr>
<tr>
<td align="left">
Notice &hellip;
</td>
</tr>
</table>
In generating his "i5.png" image, precisely how did Howard "stretch" the ordinate from <math>c/a</math> (as used in his "i3.png" image) to an ordinate ranging from the lower to the upper self-adjoint sequences?  Drawing from {{ LL96 }} I presume that, for a given point in the EFE diagram <math>(b/a, c/a)</math>, Howard used the expression,
<table border="0" cellpadding="5" align="center">
<tr>
  <td align="right">
<math>x_\pm</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> - C \pm \sqrt{C^2 - 1} \, ,</math>
  </td>
</tr>
<tr>
  <td align="center" colspan="3">
{{ LL96 }}, §2, p. 701, Eq. (8)
  </td>
</tr>
</table>
where,
<table border="0" cellpadding="5" align="center">
<tr>
  <td align="right">
<math>C</math>
  </td>
  <td align="center">
<math>\equiv</math>
  </td>
  <td align="left">
<math> \frac{abB_{12}}{c^2 A_3 - b^2 a^2 A_{12}} \, ,</math>
  </td>
</tr>
<tr>
  <td align="center" colspan="3">
{{ LL96 }}, §2, p. 701, Eq. (6)
  </td>
</tr>
</table>
Then I presume that the ordinate, <math>y</math> &#8212; which runs from zero to unity in the "i5.png" image &#8212; is determined from the expression,
<table border="0" cellpadding="5" align="center">
<tr>
  <td align="right">
<math>y</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>0.5(1 - x_\pm)  \, .</math>
  </td>
</tr>
</table>
<font color="red">Is this the way Howard generated "i5.png"?</font>


<table border="1" align="center" width="80%" cellpadding="10"><tr><td align="left">
In an email dated 26 January 2022, Howard provided the following answer to this question &#8212;

<ul>
  <li>
I had numerical data that I think Norman provided me for the lower and upper self-adjoint sequences.
  </li>
  <li>
They were simply two sets of curvilinear data in the (b/a,c/a) diagram. Call ay=b/a and then az=c/a.
  </li>
  <li>
Then lsa and usa are both functions of ay. Note that when I encountered points which didn't lie on Norman's data, I interpolated using a 9th degree polynomial  to the lower and upper self-adjoint sequence data.
  </li>
  <li>
Now consider that you have some data "g" which gives you points in the (ay,az) plane, then g=g(ay,az).
  </li>
  <li>
Take for instance data "g" which are points in the (ay,az) plane where the growth rate are above some critical value such as 10e-5. 
  </li>
  <li>
For every data point g, there is a fixed ay coordinate value.  Normally you would plot that point at (ay,az).
  </li>
  <li>
The remapping that I did now plots it instead at a point on the ordinate given by some az'= (az-lsa(ay))/(usa(ay)-lsa(ay))
  </li>
  <li>
So if az=lsa(ay) then it appears at the ordinate value of zero.
  </li>
  <li>
and if az=usa(ay) then it appears at the ordinate value of unity.
  </li>
  <li>
So the whole horn shaped region is mapped into the unit square.
  </li>
</ul>
</td></tr></table>

===Image i4.png===
Howard's "i4.png" image, immediately below, presents a magnification of the upper-right-hand portion (identified, by hand, as the "E-group") of his "i5.png" image. The abscissa spans the parameter range, <math>0.3 \le b/a \le 1.0</math> while the ordinate spans the parameter range, <math>0.96 \le y \le 1</math>. 
<table border="1" align="center" cellpadding="8">
<tr>
  <td align="center">Howard's "i4.png" image</td>
  <td align="center" rowspan="3">[[File:HighResSelfAdjointi4.png|400px|High Resolution]]</td>
</tr>
<tr><td align="center">
[[File:I4 FromCohl.png|450px|i4.png]]
</td></tr>
<tr>
<td align="left">
Notice &hellip;
</td>
</tr>
</table>

===Summary===
Howard is interested in understanding &#8212; in greater detail than appears in {{ LL96 }} &#8212; what gives rise to, and what is the extent of these various bands of instability in the classic EFE diagram. Explicit comments/questions:
<ol>
<li>Notice in "i4.png" that the bands labeled E4, E6, and E8 appear to extend all the way to, and intersect, the Maclaurin spheroid sequence.</li>
</ol>

==Self-Adjoint Sequences==

<table border="1" align="center" width="80%" cellpadding="10"><tr><td align="left">
In an email dated 26 January 2022, Howard asked, "Do you have analytic curves for the lower and upper self-adjoint sequences? Otherwise, do you have very accurate data for the lower and upper self-adjoint sequences?"
</td></tr></table>

On the same day, I sent the following response to Howard:


I have added a subsection to my online chapter discussion of {{ LL96 }} in which I derive an expression whose solution/root should map out the **upper** boundary (x = -1) of the horned-shape region.  [[ThreeDimensionalConfigurations/Stability/RiemannEllipsoids#Upper_Boundary|Click here]] to see the entire derivation; this derivation ends with the following recommended strategy:

<table border="1" align="center" width="80%" cellpadding="8"><tr><td align="left">
<font color="red">STRATEGY</font> for finding the locus of points that define the upper boundary of the horned-shape region &hellip; &nbsp; &nbsp;Set <math>a = 1</math>, and pick a value for <math>0 < b < 1</math>; then, using an iterative technique, vary <math>c</math> until the following expression is satisfied:

<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>
\biggl[c^2(a+b)+ ab^2 \biggr]
\biggl[  \frac{(b/a) \sin\theta - (c/a)E(\theta,k)}{(1-k^2) \sin^3\theta} \biggr]
</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> a^2b + bc(a  - b)
\biggl[  \frac{F(\theta,k) - E(\theta,k)}{k^2 \sin^3\theta} \biggr]
\, .</math>
  </td>
</tr>
</table>
Choose another value of <math>0 < b < 1</math>, then iterate again to find the value of <math>c</math> that corresponds to this new, chosen value of <math>b</math>.  Repeat!
</td></tr></table>


Related remarks:
<ol>
<li>I have not actually plugged in numbers -- that is, (b,c) pairs -- to see if it works, but I am pretty confident in the result because the derivation was pretty straightforward.  Would you mind trying it out for me, since you have working elliptic integral routines?</li>
<li>It would be wise to start by trying to duplicate -- then improve upon -- the set of (b, c) coordinate-pairs that were derived by Chandrasekhar and presented in EFE Table VI (section 48, p. 142).</li>
<li>Shortly, I will derive the complementary expression that maps out the "lower" boundary (x = +1).</li>
</ol>

<!--
<table border="1" align="center" cellpadding="10"><tr><td align="left">
<div align="center"><b>Upper Self-Adjoint</b><br />(Howard's email on 1/27/2022)</div>
<div align="center">
[[File:EmailFromHowardJan27Yr2022.png|650px|FromHowardJan27Yr2022]]
</div>
</td></tr></table>
-->

On 27 January 2022, Joel added a subsection to the online chapter discussion of {{ LL96 }} in which he derives an expression whose solution/root should map out the **lower** boundary (x = +1) of the horned-shape region.  [[ThreeDimensionalConfigurations/Stability/RiemannEllipsoids#Lower_Boundary|Click here]] to see the entire derivation; this derivation ends with the following recommended strategy:

<table border="1" align="center" width="80%" cellpadding="8"><tr><td align="left">
<font color="red">STRATEGY</font> for finding the locus of points that define the lower boundary of the horned-shape region &hellip; &nbsp; &nbsp;Set <math>a = 1</math>, and pick a value for <math>0 < b < 1</math>; then, using an iterative technique, vary <math>c</math> until the following expression is satisfied:

<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>
\biggl[ c^2(a-b)+ ab^2 \biggr]
\biggl[  \frac{(b/a) \sin\theta - (c/a)E(\theta,k)}{(1-k^2) \sin^3\theta} \biggr]
</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>a^2 b - bc(b + a) 
\biggl[  \frac{F(\theta,k) - E(\theta,k)}{k^2 \sin^3\theta} \biggr] \, .
</math>
  </td>
</tr>
</table>
Choose another value of <math>0 < b < 1</math>, then iterate again to find the value of <math>c</math> that corresponds to this new, chosen value of <math>b</math>.  Repeat!
</td></tr>
<!--
<tr><td align="center">
<div align="center"><b>Lower Self-Adjoint</b><br />(Howard's derivation 1/28/2022)</div>
<div align="center">
[[File:EmailFromHowardJan28Yr2022.png|650px|Lower Self-Adjoint Definition]]
</div>
</td></tr>
-->
</table>
&nbsp;<br />

<div id="HowardHighResolution">
<table border="1" align="center" cellpadding="10" width="80%">
<tr><td align="center">
<div align="center"><b>Upper (USA) and Lower (LSA) Self-Adjoint</b><br />(Howard's Compact Expressions Plus Plot 1/28/2022)</div>
<div align="center">
[[File:SelfAdjointDefinitions.png|650px|Self Adjoint Definitions]]<br />[[ThreeDimensionalConfigurations/RiemannStype#SAdata|(tabulated data here)]]
</div>
</td></tr>
<tr><td align="center" bgcolor="black">
[[File:SelfAdjointPlot.png|600px|USA and LSA Plot]]
</td></tr>
</table>
</div>

==COLLADA 3D Animations==
<!-- 
When you have a chance (this is certainly not urgent), please go to the following web page on the lsu physics server:
[http://phys.lsu.edu/tohline/COLLADA/ http://phys.lsu.edu/tohline/COLLADA/]

The yellow boxes in Table 1 on this web page contain links to eight separate 3D animated models of Riemann S-type ellipsoids. (Actually there are four different models, but for each there is a model as viewed from the inertial frame and there is a second model as viewed from a rotating frame.)  In each case, the (binary) file type is ".glb".  If you have a "3D Viewer" application on your PC (hopefully it is a standard installation on your PC), you should be able to open any one of these files with "3D Viewer" and immediately see the three-dimensional structure and flow in action.  Some additional helpful information is provided, both preceding and following Table 1, on this "COLLADA" web page.
-->


<table border="3" align="center" cellpadding="10"><tr><td align="center" bgcolor="red">
<table border="0" cellpadding="8" align="center">
<tr>
  <td align="center">
[[Appendix/Ramblings/COLLADA/RiemannSType|<font color="white" size="+1">DOWNLOADABLE 3D MODELS</font>]]
  </td>
</tr>
</table>
</td></tr></table>

==Riemann Type I Ellipsoids==

[[ThreeDimensionalConfigurations/RiemannTypeI#Figure3|Click here.]]

==Excerpts From Various Lebovitz Publications Regarding Stability Studies==

===Why is ''Lagrangian'' Perturbation Method Preferred?===
<table border="0" cellpadding="3" align="center" width="80%">
<tr><td align="left">
<font color="darkgreen">
"… a preliminary version</font> &#8212; see Lebovitz (1987) &#8212; <font color="darkgreen">of the methods of the present paper [have shown that] &hellip; they enjoy certain advantages over other methods that have been employed in related problems (like Cartan's method and the virial method)."
</font>
</td></tr>
<tr><td align="right">
— Drawn from the last paragraph on p. 222 of {{ Lebovitz89a }}. 
</td></tr></table>

<table border="0" cellpadding="3" align="center" width="80%">
<tr><td align="left">
<font color="darkgreen">
[<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>] "&hellip; utilizes the virial equations, which require a separate derivation for modes associated with ellipsoidal harmonics of various orders; it is accordingly restricted to those of low order (at most four).."
</font>
</td></tr>
<tr><td align="right">
— Drawn from the first paragraph of &sect;1 (p. 225) in {{ Lebovitz89b }}. 
</td></tr></table>

==Contrast with Radial Pulsation Eigenvectors==

<table border="1" align="center" cellpadding="5">
<tr>
  <td align="center" colspan="2">
[[SSC/Stability/n3PolytropeLAWE#Schwarzschild_(1941)|Schwarzschild's 1941 Analysis of Radial Oscillations in n = 3 Polytropes]]
  </td>
</tr>
<tr>
  <td align="center">
{{ Math/EQ_RadialPulsation02 }}
  </td>
  <td align="center">
[[File:Schwarzschild1941movie.gif|300px|link=SSC/Stability/n3PolytropeLAWE#SchwarzschildMovie|Schwarzschild's Modal Analysis]]
  </td>
</tr>
</table>


By contrast, in {{ LL96hereafter }}, <font color="green">"the basic equation"</font> appears in the form,

<table border="0" align="center" cellpadding="5">
<tr>
  <td align="right"><math>0</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\boldsymbol{\xi}_{tt} +
A {\boldsymbol\xi}_t + B \boldsymbol\xi + \rho^{-1} \nabla ( \Delta p ) \, .
</math>
  </td>
</tr>

<tr>
  <td align="center" colspan="3">
{{ LL96hereafter }}, &sect;3.1, p. 701, Eq. (10)
  </td>
</tr>
</table>
This means that the matrix operators, <math>M</math> &amp; <math>\Lambda</math>, found in {{ Lebovitz89b }} and [[ThreeDimensionalConfigurations/Stability/RiemannEllipsoids#Stability_Equations|re-derived herein]], have simply been renamed in {{ LL96hereafter }}.  That is to say, in {{ LL96hereafter }},

<table border="0" align="center" cellpadding="5">
<tr>
  <td align="right"><math>A\boldsymbol\xi</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
2\{ (\mathbf{u}\cdot \nabla)\boldsymbol\xi 
+ \boldsymbol\omega \boldsymbol\times \boldsymbol\xi \} \, ,
</math>
  </td>
  <td align="center">&nbsp; &nbsp; &nbsp; and, &nbsp; &nbsp; &nbsp;</td>
  <td align="right"><math>B\boldsymbol\xi</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\{ (\mathbf{u}\cdot \nabla)^2\boldsymbol\xi 
+ 2\boldsymbol\omega \boldsymbol\times [(\mathbf{u}\cdot \nabla)\boldsymbol\xi ] \} 
- L \boldsymbol\xi \, ,
</math>
  </td>
</tr>
</table>
where,

<table border="0" align="center" cellpadding="5">
<tr>
  <td align="right"><math>L \boldsymbol\xi</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\rho^{-1} \nabla ( \boldsymbol\xi \cdot \nabla p )
+
\nabla \delta \Phi_\mathrm{L89}
+ ( \boldsymbol\xi \cdot \mathbf\nabla ) [ 
(\mathbf{u} \cdot \nabla)\mathbf{u} + 2\boldsymbol\omega \boldsymbol\times \mathbf{u} 
] \, .
</math>
  </td>
</tr>

<tr>
  <td align="center" colspan="3">
{{ LL96hereafter }}, &sect;3.1, p. 703, Eq. (17)
  </td>
</tr>
</table>

==Determining Oscillation Frequencies and Growthrates==

6 February 2022: &nbsp; Howard asked me to clarify how the solution to the "basic linear perturbation equation" provides values for the growth rate of an unstable mode.

===Spherically Symmetric Situations===

In a [[SSC/Perturbations#The_Eigenvalue_Problem|separate chapter]] I have explained how the eigenvalue problem is set up for one-dimensional (spherically symmetric) configurations.  Traditionally, the assumption is that you are starting from an equilibrium configuration for which you know how pressure <math>(P_0)</math>, density <math>(\rho_0)</math>, and radius <math>(r_0)</math> vary with mass shell <math>(m)</math>,
everywhere inside (and on the surface) of the configuration. As the following three equations illustrate, you "perturb" the configuration by assuming that there are low-amplitude fluctuations to each variable, namely, <math>P_1</math>, <math>\rho_1</math>, and <math>r_1</math>, and that these fluctuations each vary with spatial location &#8212; that is, with mass shell &#8212; inside the configuration as well a s with time. 

<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~P(m,t)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~P_0(m) + P_1(m,t) = P_0(m) \biggl[1 + p(m) e^{i\omega t}  \biggr] \, ,</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\rho(m,t)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\rho_0(m) + \rho_1(m,t) = \rho_0(m) \biggl[1 + d(m) e^{i\omega t}  \biggr] \, ,</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~r(m,t)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~r_0(m) + r_1(m,t) = r_0(m) \biggl[1 + x(m) e^{i\omega t}  \biggr] \, ,</math>
  </td>
</tr>
</table>
</div>
The convention is to assume, as shown, that the spatial portion of the variation can be separated from the time-dependent portion, and that the time-dependent portion is of the form, <math>e^{i\omega t}</math>.  Then, when these quantities are introduced into the governing linear perturbation equation, any term involving <b>one</b> partial time-derivative &#8212; for example, <math>\partial r/\partial t</math> &#8212; will generate a term of the form, <math>i\omega \cdot e^{i\omega t}</math>; while any term involving <b>a second</b> partial time-derivative will generate a term of the form, <math>-\omega^2 \cdot e^{i\omega t}</math>.  Usually, every term in the governing linear perturbation equation will contain a factor of <math>e^{i\omega t}</math>, so it can be divided out.  This leaves you with a perturbation equation that has no explicit time dependence; it only contains spatial derivatives of variables along with a few <math>\omega</math> or <math>\omega^2</math> terms.

The ''very nature'' of an eigenvalue problem is to see if the perturbation equation can be solved to give you the square of the characteristic oscillation frequency, <math>\omega^2</math> for one (or, hopefully more) mode ''along with'' a specification of how the spatial part of the mode varies with location inside the configuration.  If <math>\omega^2</math> is positive, then the exponent of <math>e^{i\omega t}</math> is imaginary, which identifies a periodic oscillation; if <math>\omega^2</math> is negative, then the exponent of <math>e^{i\omega t}</math> is real, which identifies an exponentially growing (or, formally as well, decaying) mode.

===Riemann S-Type Ellipsoids===

Notice that Eq. (25) of {{ Lebovitz89ahereafter }} assumes that the three-dimensional, time-dependent perturbation, <math>\boldsymbol\xi</math>, of the Riemann S-Type ellipsoid is of the form,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>\boldsymbol\xi(\mathbf{x},t)</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>\sum_{k=1}^n a_k(t) \xi_k(\mathbf{x})\, .</math>
  </td>
</tr>
</table>
Notice that time-dependence and spatial-dependence have been separated.  And I presume that the time-dependent coefficients, <math>a_k(t)</math>, will include factors of <math>e^{i\omega t}</math>.  Notice as well that, in {{ LL96hereafter }}, <font color="green">"the basic equation"</font> appears in the form,

<table border="0" align="center" cellpadding="5">
<tr>
  <td align="right"><math>0</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\boldsymbol{\xi}_{tt} +
A {\boldsymbol\xi}_t + B \boldsymbol\xi + \rho^{-1} \nabla ( \Delta p ) \, .
</math>
  </td>
</tr>

<tr>
  <td align="center" colspan="3">
{{ LL96hereafter }}, &sect;3.1, p. 701, Eq. (10)
  </td>
</tr>
</table>
The subscript <math>tt</math> means that the perturbation will be twice differentiated in time, while the subscript <math>t</math> means a single time-differentiation.  The oscillation frequency (or growth rate) probably will appear in the basic perturbation equation via these two terms.

==Index Symbols for Type I Riemann Ellipsoids==

2 March 2022: &nbsp; Joel requests Howard's assistance in verifying the correct expressions to use for the index symbols, <math>(A_1, A_2, A_3)</math>, when constructing "Type I" Riemann ellipsoids.

As has been [[#Three_Lowest-Order_Expressions|summarized above]], when I have derived expressions for the three lowest-order index symbols, I have used <math>(\ell, m, s)</math> instead of <math>(1, 2, 3)</math> as the subscript notation.  In the context of our discussion of Riemann S-type ellipsoids, it is usually the case that <math>a_1</math> is the longest semi-axis and <math>a_3</math> is the shortest semi-axis, so we can adopt the association, <math>(\ell, m, s) \rightarrow (1, 2, 3)</math>. This is consistent with the set of expressions for <math>A_1</math>,  <math>A_2</math>, and  <math>A_3</math>, that Howard has defined inside of his Mathematica notebook, namely &hellip;

<table border="1" align="center" cellpadding="8"><tr><td align="left">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>A_1</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{2a_2 a_3}{a_1^2 m \sin^3(\phi) } \biggl[
\mathrm{EllipticF}[\phi, m] - \mathrm{EllipticE}[\phi, m] 
\biggr]
\, ;
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>A_2</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{2a_2 a_3}{a_1^2 m(1-m) \sin^3(\phi) } \biggl[
\mathrm{EllipticE}[\phi, m] - \cos^2\theta \cdot \mathrm{EllipticF}[\phi, m] - \frac{a_3}{a_2}\cdot\sin^2\theta \sin\phi
\biggr]
\, ;
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>A_3</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{2a_2 a_3}{a_1^2 (1-m) \sin^3(\phi) } \biggl[
\frac{a_2}{a_3}\cdot \sin(\phi) - \mathrm{EllipticE}[\phi, m] 
\biggr]
\, .
</math>
  </td>
</tr>
</table>
</td></tr></table>

But I am also interested in developing a better understanding of the structural properties of Type I Riemann Ellipsoids.  For these systems, it is the case that <math>a_2</math> is the longest semi-axis and <math>a_3</math> is the shortest semi-axis, so we should adopt the association, <math>(\ell, m, s) \rightarrow (2, 1, 3)</math>.  In <font color="red">STEP #2</font> of the chapter that I am writing [[3Dconfigurations/DescriptionOfRiemannTypeI#Description_of_Riemann_Type_I_Ellipsoids|about Type I Riemann Ellipsoids]], I claim that the correct expressions for the three lowest-order index symbols are &hellip;
<div align="center">
<table align="center" border=0 cellpadding="3">
<tr>
  <td align="right">
<math>
A_2
</math>
  </td>
  <td align="center">
<math>
=
</math>
  </td>
  <td align="left">
<math>2 \biggl( \frac{a_1}{a_2} \biggr)\biggl( \frac{a_3}{a_2} \biggr)
\biggl[  \frac{F(\theta,k) - E(\theta,k)}{k^2 \sin^3\theta} \biggr] \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>
A_3
</math>
  </td>
  <td align="center">
<math>
=
</math>
  </td>
  <td align="left">
<math>
2\biggl( \frac{a_1}{a_2}\biggr) \biggl[  \frac{(a_1/a_2) \sin\theta - (a_3/a_2)E(\theta,k)}{(1-k^2) \sin^3\theta} \biggr]  \, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>
A_1 = 2 - (A_2 + A_3)
</math>
  </td>
  <td align="center">
<math>
=
</math>
  </td>
  <td align="left">
<math>
\frac{ 2a_1 a_3}{a_2^2 }
\biggl[ \frac{
E(\theta, k) 
-~(1-k^2)
F(\theta, k)
-~(a_3/a_1)k^2\sin\theta}{k^2 (1-k^2)\sin^3\theta}
\biggr] \, ,
</math>
  </td>
</tr>

</table>
</div>
where, the arguments of the incomplete elliptic integrals are,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>\theta = \cos^{-1} \biggl(\frac{a_3}{a_2} \biggr)</math>
  </td>
  <td align="center">
&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp; 
  </td>
  <td align="left">
<math>k = \biggl[\frac{1 - (a_1/a_2)^2}{1 - (a_3/a_2)^2} \biggr]^{1/2} \, .</math>
  </td>
</tr>
</table>
</div>
I would like to know if Howard agrees with this claim.  In particular, I would like to know what values of <math>(A_1, A_2, A_3)</math> he obtains using Mathematica's tools, for the following pairs of model axis ratios.

<table border="1" align="center" cellpadding="5">
<tr>
  <td align="center" colspan="5">
Example Models Drawn ''Primarily''<sup>&dagger;</sup> from Table 4 (p. 858) of  {{ Chandrasekhar66_XXVIII }}
  </td>
</tr>
<tr>
  <td align="center" rowspan="2"><math>\frac{a_2}{a_1}</math></td>
  <td align="center" rowspan="2"><math>\frac{a_3}{a_1}</math></td>
  <td align="center" colspan="3">Howard Cohl's 20-Digit Precision Evaluation of Index Symbols</td>
</tr>
<tr>
  <td align="center" bgcolor="red" width="25%"><font color="white">A<sub>1</sub></font></td>
  <td align="center" bgcolor="red" width="25%"><font color="white">A<sub>2</sub></font></td>
  <td align="center" bgcolor="red" width="25%"><font color="white">A<sub>3</sub></font></td>
</tr>
<tr>
  <td align="right">1.05263</td>
  <td align="right">0.41667</td>
  <td align="left">0.43008853859109863146</td>
  <td align="left">0.40190463270335853286</td>
  <td align="left">1.1680068287055428357</td>
</tr>
<tr>
  <td align="right">1.25000</td>
  <td align="right">0.50000</td>
  <td align="left">0.50824409128926609544</td>
  <td align="left">0.37944175746381924039</td>
  <td align="left">1.1123141512469146642</td>
</tr>
<tr>
  <td align="right">1.44065</td>
  <td align="right">0.49273</td>
  <td align="left">0.52404662956445493304</td>
  <td align="left">0.32351600902859843592</td>
  <td align="left">1.1524373614069466310</td>
</tr>
<tr>
  <td align="right">1.66666</td>
  <td align="right">0.33333</td>
  <td align="left">0.41804279087942201679</td>
  <td align="left">0.20718018294728778618</td>
  <td align="left">1.3747770261732901970</td>
</tr>
<tr>
  <td align="right">1.36444</td>
  <td align="right">0.09518</td>
  <td align="left">0.14378468450707924448</td>
  <td align="left">0.091526900363135453419</td>
  <td align="left">1.7646884151297853021</td>
</tr>
<tr>
  <td align="right">1.69351</td>
  <td align="right">0.11813</td>
  <td align="left">0.18177227461540501422</td>
  <td align="left">0.084646944102441440238</td>
  <td align="left">1.7335807812821535455</td>
</tr>
<tr>
  <td align="right">1.52303</td>
  <td align="right">0.05315</td>
  <td align="left">0.085943612594485367787</td>
  <td align="left">0.046184257303756474717</td>
  <td align="left">1.8678721301017581575</td>
</tr>
<tr>
  <td align="right">1.78590</td>
  <td align="right">0.06233</td>
  <td align="left">0.10259594310295396216</td>
  <td align="left">0.043583894016884923661</td>
  <td align="left">1.8538201628801611142</td>
</tr>
<tr>
  <td align="right" bgcolor="yellow">1.2500</td>
  <td align="right" bgcolor="yellow">0.4703</td>
  <td align="left">0.48955940032702523984</td>
  <td align="left">0.36486593343389634429</td>
  <td align="left">1.1455746662390784159</td>
</tr>
<tr>
  <td align="left" colspan="5">
<sup>&dagger;</sup>NOTE: &nbsp; The model whose axis-ratios are highlighted with a yellow background has been drawn from Table 6a (p. 871) of  {{ Chandrasekhar66_XXVIII }}
  </td>
</tr>
</table>

If you want to know how I derived the relevant index-symbol expressions, go to my MediaWiki chapter that discusses ''The Gravitational Potential (A<sub>i</sub> coefficients)'' in the context of <b>Ellipsoidal &amp; Ellipsoidal-Like</b> equilibrium structures.

<table border="1" align="center" cellpadding="5">
<tr>
  <td align="center">[[File:TiledMenuIndexSymbolEvaluations.png|600px|Tiled Menu Chapter titled "Index Symbol Evaluations"]]
</tr>
</table>

More specifically, go to the subsection of this chapter titled, "[[ThreeDimensionalConfigurations/HomogeneousEllipsoids#Derivation_of_Expressions_for_Ai|Derivation of Expressions for A<sub>i</sub>]]," where I evaluate <math>A_\ell, A_m, A_s</math>.  These expressions will always remain the same.  Then the question is, for your specific problem of interest, which ellipsoidal axis <math>a_1, a_2,</math> or <math>a_3</math> is the "largest, medium-sized, or smallest"?  When considering Riemann Type-I Ellipsoids, <math>a_2</math> is the largest axis, so <math>A_2 \leftrightarrow A_\ell</math>; <math>a_3</math> is the smallest axis, so <math>A_3 \leftrightarrow A_s</math>; and <math>a_1</math> is the medium-sized axis, so <math>A_1 \leftrightarrow A_m</math>.

=Considering Utility of Finite-Element Techniques=

==Joel's Opening Comments==
Among the introductory chapters of our [[H_BookTiledMenu|MediaWiki-based discussion of ''Self-Gravitating Fluids'']], you will find a formal [[PGE|list of the relevant "Principal Governing Equations]]" and an accompanying discussion of, for example, the [[PGE/Euler|Euler equation]].  Here, we pull from that discussion the so-called,

<div align="center">
<span id="ConservingMomentum:Conservative"><font color="#770000">'''Conservative Form'''</font></span><br />
of the Euler Equation,

<math>\frac{\partial(\rho\vec{v})}{\partial t} + \underbrace{\nabla\cdot [(\rho\vec{v})\vec{v}]}_{\mathrm{advection}} 
= 
\overbrace{- \nabla P - \rho \nabla \Phi}^{\mathrm{source}}</math>

[<b>[[Appendix/References#BLRY07|<font color="red">BLRY07</font>]]</b>], p. 8, Eq. (1.31)
</div>

which contains the divergence of the "[https://en.wikipedia.org/wiki/Dyadics dyadic product]" or "[https://en.wikipedia.org/wiki/Outer_product outer product]" of the vector momentum density and the velocity vector, and is sometimes written as, <math>\nabla\cdot [(\rho \vec{v}) \otimes \vec{v}]</math>.  How do we rewrite this ''continuum'' equation in a ''discrete'' form that will allow us (using digital computers) to accurately integrate this equation &#8212; in tandem with the other "Principal Governing Equations" &#8212; to determine how the fluid momentum at every point in space evolves in time?  Answering this question is a particularly challenging task because the so-called ''advection'' term on the left-hand-side is inherently nonlinear in the velocity.

Over the years, applied mathematicians &#8212; including, for example, engineers, physicists, astronomers, aerodynamicists &#8212; have devised a large array of schemes, the most popular falling under one of the following three broad (Wikipedia-acknowledged!) categories:
<ul>
  <li>
[https://en.wikipedia.org/wiki/Finite_difference_method Finite-Difference Method] (FDM)
  </li>
  <li>
[https://en.wikipedia.org/wiki/Finite_volume_method Finite-Volume Method] (FVM) &nbsp; &nbsp; <math>\Leftarrow</math> &nbsp; Our category of choice, historically (see immediately below)
  </li>
  <li>
[https://en.wikipedia.org/wiki/Finite_element_method Finite-Element Method] (FEM)
  </li>
</ul>

When I googled "finite element vs finite difference " on 13 October 2022, I was pointed to, for example,
<ul>
  <li>
[https://www.sciencedirect.com/science/article/pii/B9780122532504500191 Comparison of Finite-Element and Finite-Difference Methods]
  </li>
  <li>
[https://www.machinedesign.com/3d-printing-cad/fea-and-simulation/article/21832072/whats-the-difference-between-fem-fdm-and-fvm What's the Difference Between FEM, FDM, and FVM?]
  </li>
</ul>

===Our Choice, Historically===

Over the past four decades, the nonrelativistic fluid simulations that have been performed by LSU's astrophysics group have been carried out using various renditions of an ''explicit'' <b>Finite-Volume Method</b>.  The explicit FVM method was chosen primarily because &hellip;
<ol>
  <li>
It was the method to which I was introduced by my [[Appendix/Ramblings/MyDoctoralStudents#Years_1976_-_1978|doctoral dissertation advisors]], Peter Bodenheimer and David Black, that was considered appropriate to the type of fluid-flows problems that were the focus of my dissertation research.  Note that Bodenheimer is the first author of the book that is [[#Joel%27s_Opening_Comments|referenced above]] as [[Appendix/References#BLRY07|<font color="red">BLRY07</font>]].
  </li>
  <li>
Time-integration is fairly straightforward because all terms (except the time derivative) are specified entirely in terms of the current (as opposed to advanced) time.
  </li>
  <li>
The nonlinear ''advection'' term on the LHS is written in a way that guarantees global momentum conservation, if the ''source'' terms on the RHS of the Euler equation sum to zero.
  </li>
  <li>
The amount of computer time that is required to take each step forward in time is typically much smaller when using an ''explicit'' as opposed to ''implicit'' integration scheme.
  </li>
</ol>

A major disadvantage of our decades-old, explicit-FVM approach is that, in order for the scheme to be numerically stable, roughly speaking each integration time-step <math>\Delta t</math> must be smaller than the smallest value of the ratio, <math>\Delta x/c_s</math>, as determined across the entire grid, where <math>\Delta x</math> is the grid spacing and <math>c_s</math> is the fluid sound-speed.  (Note that <math>\Delta t</math> is essentially the sound-crossing time across the smallest grid cell.)  This constraint on the integration time-step is computationally burdensome when attempting to employ a high-resolution (small <math>\Delta x</math>) grid and/or when attempting to model environments where the fluid is relatively hot and therefore characterized by a high sound speed.  Even worse, <b>it is ''impossible'' to model the evolution of ''incompressible'' fluids using an explicit-FVM scheme</b> because, in effect in such fluids, sound travels at an infinite speed; that is, <math>c_s = \infty</math> so <math>\Delta t</math> goes to zero.

===Quandary===

Because we have chosen over the past four decades to employ an explicit-FVM approach in our examination of the stability of rotating, self-gravitating fluids, we have been constrained to examine the behavior of ''compressible'' rather than ''incompressible'' fluids.  This is perfectly acceptable &#8212; and even desirable &#8212; in the context of studies of star-formation, for example, because we are confident that star-forming gas clouds and young protostars are composed of fluids that obey compressible equations of state.  But it is difficult to critique the validity &#8212; and therefore relevance &#8212; of the results of such compressible-fluid simulations because <b>the published literature contains virtually no linear-stability analyses of compressible fluid systems</b> against which our nonlinear simulations can be compared.

On the other hand, as Chandrasekhar details in [[Appendix/References#EFE|<font color="red">EFE</font>]], the literature is packed full of studies that have quantitatively assessed the relative stability of rotating, self-gravitating, ''incompressible'' fluid systems.  (And, associated with the arrival of Stephen Sorokanich at NIST, we expect this literature to expand significantly.)  So, the validity &#8212; and relevance &#8212; of our numerical simulations would be on a much firmer foundation if we had the ability to model the evolution of incompressible fluid systems.  As an additional reward, a numerical simulation tool of this type would allow us to follow the ''nonlinear'' development of instabilities (in all Riemann ellipsoids, for example) whose onset has been predicted only from ''linear'' stability analyses.

===Perhaps We Should Consider the FEM===

Although I have had virtually no experience developing or using numerical algorithms that fall into the category of the Finite-Element Method (FEM), it is my understanding that such algorithms have the following features:
<ul>
  <li>
They can be used to follow the evolution of ''incompressible'' &#8212; and, for example, uniform-density &#8212; fluid systems.
  </li>
  <li>
They generally (always?) incorporate an ''implicit'' time-integration scheme in which case time steps are not constrained by &#8212; and can be much larger than &#8212; the often severely limiting sound-crossing times encountered in ''explicit'' FVM schemes.  
  </li>
  <li>
They are particularly good at identifying ''surfaces'' and modeling the time-evolution of surface distortions.
  </li>
</ul>

If you (Howard) and Sorokanich decide to build a hydro-code based on the Finite-Element Method, <font color="red">a potentially helpful reference</font> is:  
<div align="center">
{{ Meier99figure }}
</div>
I interacted with David Meier at a couple of different meetings back at the time he was developing this FEM algorithm. I'm not sure what astrophysics problems he tackled with this code or, ultimately, what he did with it; for example, I have never seen "the second paper in this series" to which he refers in the abstract.

In the context of our examination of the stability of Riemann ellipsoids, an FEM might permit us to reduce the dimensionality of the problem.  For example, ignore details associated with the interior of the 3D configuration and focus on modeling the distortion of its 2D surface.  

We might even be able to model (for the first time!) the topological transformation that is experienced by the surface when a rapidly spinning, incompressible ellipsoid ''fissions'' into a pair of disconnected tear-drop shaped objects in orbit about one another.  See, for example, the dissertation research of [https://www.semanticscholar.org/paper/A-hybrid-variational-level-set-approach-to-handle-Walker/7ee624ab9ffe45241cb0e5a0ce0898a6da201e7b Shawn W. Walker (2007)] (LSU Mathematics &amp; CCT), or more recent papers that Walker has coauthored with his dissertation advisor, Ricardo H. Nochetto (U. Maryland, Mathematics).

==Question Regarding Pressure-Poisson Equation==

<table border="1" align="center" width="80%" cellpadding="10"><tr><td align="left">
In an email dated 07 January 2023, Howard asked, "Do you think the pressure Poisson formulation of the Euler equations coupled with a 3D Poisson solver for gravity and evolution of the momentum equation might be our best strategy to evolve unstable Riemann S-type ellipsoids? We’ve already identified two different pieces of software (Fenics and COMSOL) which look like they might have this ability using finite-element. We also now have access to software which should be able to build good grids for FEM. Thanks! "
</td></tr></table>

Joel's response (11 January 2023) &hellip;

Your suggested "coupled Poisson formulation" sounds like a smart way to proceed.  Note, however, that I say this with limited authority; I have never written a code to solve the Euler equations that is based on a pressure-Poisson formulation.  My confidence in offering support for this approach comes from the following historical note.


In the early 1990s, in collaboration with Priya Vashishta, I helped the LSU Physics & Astronomy department secure approximately $0.8 million from the Louisiana Board of Regents to purchase [[Appendix/Ramblings/MyDoctoralStudents#Years_1988_-_1994|an 8K-node SIMD architecture MasPar MP1]].  As you know, our astrophysics group developed a 3D Poisson solver that was extraordinally well suited to the 8K-node hardware layout of the MasPar.  Around the same time that we received this ''state'' funding, I heard that the U.S. Department of Energy (DOE) was expecting to fund a limited number of large-scale fluid-flow (and heat-transfer) simulations in support of the aircraft industry's efforts to design better turbine engines.  LSU's research office organized a meeting of (primarily engineering) faculty to see if a group could be put together to submit a viable research proposal to the DOE.  I attended this meeting, just for grins; but it turned out to be a very smart decision.  At the meeting, I met [https://www.lsu.edu/eng/mie/people/formerfaculty/acharya.php Sumanta Acharya] (and [https://www.lsu.edu/eng/mie/people/faculty/nikitopoulos.php Dimitris Nikitopoulos]) and learned about his multidimensional simulations of incompressible fluids; he was looking for funding opportunities that would support his efforts to scale from 2D to 3D simulations.  I also learned that he usually solved the Euler equations using a pressure-Poisson approach.  We both realized that my group's successful implementation of a 3D Poisson solver (for gravity) on the MasPar could serve as an illustration of how his group could likely efficiently solve the pressure-Poisson equation in 3D.  We wrote the proposal to DOE, and it was funded at $0.5 million.  Thus began a very successful, multi-decade collaboration between my group and Sumanta's mechanical engineering group.


So &hellip; I am fairly confident that the approach you take to solve the gravitational-Poisson equation will help you understand how to solve the pressure-Poisson equation (or vise versa) on the finite-element grid.  Note, however, that to my knowledge, Sumanta's simulations always involved a finite-difference, rather than a finite-element, methodology.

=Recommended Riemann S-Type Ellipsoids for Modeling=
==Introductory Remarks==

<table border="1" align="center" width="65%" cellpadding="5">
<tr>
<td align="left">
On 4 January 2023, Shawn Walker reintroduced the following strategy:

"<font color="darkgreen">But if I remember correctly, I think the idea was to start with some baby steps.  Perhaps have a fixed domain (the star) and solve the Euler equations in this with perhaps a fixed “fake” gravitational potential.  Just to get something working.  Then add the Poisson equation for the gravity, which would couple with the fluid.</font>"
</td>
</tr>
</table>

Joel's response:  I am completely with you regarding the idea of starting with baby steps.  But if we totally focus on Riemann S-Type ellipsoids, even the first few baby steps can be very realistic.  Immediately below is a table that itemizes my recommendation regarding which ellipsoids should be modeled fist, second, third, etc. and why.  What follows is a short preview.

Every S-Type configuration is a uniform-density ellipsoid (semi-axes: a, b, c) whose gravitational potential -- inside, outside, and on the surface of the ellipsoid -- is known analytically in terms of incomplete elliptic integrals of the 1st and 2nd kind (or simpler, trigonometric functions).  When the chosen ellipsoid is viewed from a frame that is rotating (about its c-axis) with a characteristic frequency, <math>\Omega</math>, the ellipsoidal configuration appears stationary.  As a result, the (analytically specified) gravitational potential does not vary with time.  For our purposes, this can naturally serve as the "fixed fake" potential to which you refer.  

As viewed from the proper (<math>\Omega</math>) rotating frame, we also know the exact solution to the Euler equations throughout every S-Type ellipsoid.  It is a steady-state flow in which all Lagrangian fluid elements move along closed elliptical orbits that &hellip;
<ul>
  <li>lie in planes parallel to the system's (''a-b'') equatorial plane;</li>
  <li>have identical ellipticities (that is, identical axis ratios, <math>b/a</math>);</li>
  <li>have identical orbital frequencies.</li>
</ul>  

So, as we try to build a code that can solve the Euler equations, we can (A) build a coordinate grid that is rotating with the desired (<math>\Omega</math>) frequency; (B) "guess" a velocity flow-field that corresponds to what we know to be the correct, steady-state solution; (C) hold fixed the analytically known gravitational potential; (D) then use the code to integrate the Euler equations forward in time and see how well the configuration maintains a steady-state flow.  We will have successfully debugged the "Euler equations" code if steady-state has been achieved for a variety of different Riemann S-Type ellipsoids.

Once the "Euler equations" code has been properly debugged, we can then add a Poisson solver; this will allow us to examine non-steady-state flows that result from natural instabilities among the set of Riemann S-Type ellipsoids.

----


Other useful attributes of Riemann S-Type ellipsoids &hellip;
<ul>
<li>Each is a (incompressible) uniform-density ellipsoid with semi-axes, <math>(a, b, c)</math>; ''usually'' without loss of generality we are able to set <math>a = 1</math>.</li>
<li>The steady-state velocity flow-field that is associated with each equilibrium configuration can be characterized by the values of two physical parameters: a frequency, <math>\Omega</math>, that is associated with the rate of the ellipsoidal figure's spin about its <math>c</math>-axis; and, a vorticity, <math>\zeta</math>, associated with the elliptical orbital motion of Lagrangian fluid elements ''inside and on the surface of'' the ellipsoid, when viewed from a frame that is rotating with the frequency, <math>\Omega</math>.</li>
</ul>

<!--
Details &hellip;
<ul>
<li>A pair of equilibrium configurations exists if the specified set of axis ratios <math>(b/a, c/a)</math> falls between the "upper-self-adjoint" (USA; <math>x = -1</math>) and "lower-self-adjoint" (LSA; <math>x=+1</math>) sequences; this "horn-shaped" sub-domain is identified in both figure panels (1a) and (1b), immediately below. </li>
<li>At each point in this <math>(b/a, c/a)</math> sub-domain, the pair of equilibrium models are distinguished from one other by the behavior of the fluid flow;</li>
<li>Each is steady-state when viewed from ... </li>
</ul>
-->

<!--
<table border="1" align="center" cellpadding="8">
<tr>
  <td align="center" colspan="3">Interesting Domains and/or Equilibrium Sequences in the <math>(b/a, c/a)</math> Diagram</td>
</tr>
<tr>
  <td align="center">
Figure 2 extracted from &hellip;<br />
{{ Chandrasekhar65_XXVfigure }}
  </td>
  <td align="center">
Boundary defined by USA and LSA sequences
  </td>
  <td align="center">
Jacobi Sequence:<br />(blue) Points defined by data in Table IV of [<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>], <font color="#00CC00">Chapter 6, &sect;39 (p. 103)</font>; <br />(red) points generated here from [[#Roots_of_the_Governing_Relation|above-defined roots of the governing relation]].
  </td>
</tr>
<tr>
  <td align="center">
<b>(1a)</b><br />
[[File:ChandrasekharFig2annotated.png|400px|Chandra Diagram]]
  </td>
  <td align="center">
<b>(1b)</b><br />
[[File:SelfAdjointPlot.png|400px|USA and LSA Plot]]
  </td>
  <td align="center">
<b>(1c)</b><br />
[[File:EFEdiagram02.png|400px|Jacobi/Dedekind Sequence]]
  </td>
</tr>
</table>
-->

==Initial Equilibrium Configurations Recommended for Dynamical Analysis==

===(Axisymmetric) Maclaurin Spheroids===
The following table lists four "critical points" drawn from Table 4 (Appendix A, p. 611) of {{ HTE87 }}; go [[Appendix/Ramblings/MacSphCriticalPoints#HTE87|here]] for additional details.

<span id="RRSTEMtable1">&nbsp;</span>
<table border="1" align="center" cellpadding="5" width="85%">
<tr>
  <td align="center" colspan="9">
<b>RRSTEM Table 1</b><br />
Established Critical-Point Models Along the (Axisymmetric) Maclaurin Spheroid Sequence
  </td>
</tr>
<tr>
  <td align="center" rowspan="2">
Model
  </td>
  <td align="center" rowspan="2">
<math>e</math>
  </td>
  <td align="center" rowspan="2">
<math>\Omega^2</math>
  </td>
  <td align="center" rowspan="2">
<math>\tau</math>
  </td>
  <td align="center" rowspan="2">
<math>j^2 = \frac{1}{3}\biggl(\frac{4\pi}{3}\biggr)^{-4 / 3} L_*^2</math>
  </td>
  <td align="center" colspan="3">
Bifurcation Characteristics &hellip;
  </td>
  <td align="center" rowspan="6">
[[File:CrossSectionsAnnotated.png|300px|Meridional Cross Sections]]
  </td>
</tr>

<tr>
  <td align="center"><sup>&dagger;</sup>Geometric Distortion</td>
  <td align="center">Angular Mom. Profile</td>
  <td align="center">Instability Type</td>
</tr>

<tr>
  <td align="center" rowspan="1">
<font color="red"><b>A</b></font>
  </td>
  <td align="center"><math>0.812670</math></td>
  <td align="center"><math>0.093557</math></td>
  <td align="center"><math>0.1375</math></td>
  <td align="center"><math>4.555\times 10^{-3}</math></td>
  <td align="left"><math>P_2^2(\eta)\cos(2\phi)</math></td>
  <td align="center">Uniform <math>\omega_0</math></td>
  <td align="center">Secular</td>
</tr>

<tr>
  <td align="center" rowspan="1">
<font color="red"><b>B</b></font>
  </td>
  <td align="center"><math>0.95289</math></td>
  <td align="center"><math>0.11006</math></td>
  <td align="center"><math>0.2738</math></td>
  <td align="center"><math>1.280\times 10^{-2}</math></td>
  <td align="left"><math>P_2^2(\eta)\cos(2\phi)</math></td>
  <td align="center"><math>n' = 0</math></td>
  <td align="center">Dynamical</td>
</tr>

<tr>
  <td align="center" rowspan="1">
<font color="red"><b>C</b></font>
  </td>
  <td align="center"><math>0.98522</math></td>
  <td align="center"><math>0.087271</math></td>
  <td align="center"><math>0.3589</math></td>
  <td align="center"><math>2.174\times 10^{-2}</math></td>
  <td align="left"><math>P_4(\eta)</math></td>
  <td align="center">Uniform <math>\omega_0</math></td>
  <td align="center">Secular</td>
</tr>

<tr>
  <td align="center" rowspan="1">
<font color="red"><b>D</b></font>
  </td>
  <td align="center"><math>0.998556</math></td>
  <td align="center"><math>0.036820</math></td>
  <td align="center"><math>0.4511</math></td>
  <td align="center"><math>4.304\times 10^{-2}</math></td>
  <td align="left"><math>P_6(\eta)</math></td>
  <td align="center"><math>n' = 0</math></td>
  <td align="center">Dynamical</td>
</tr>

<tr>
  <td align="left" colspan="9">
<sup>&dagger;</sup>Following {{ EH85full }} &#8212; see especially their Appendix and their Table 2 &#8212; <math>P_{2n}</math> is the Legendre polynomial.  See also, p. 429 of {{ Bardeen71full }}.
  </td>
</tr>

<tr>
  <td align="left" colspan="9">
Given the value of the eccentricity, <math>e</math>, we immediately know from the [[Apps/MaclaurinSpheroidSequence#Maclaurin_Spheroid_Sequence|accompanying detailed discussion]] that &hellip;

<table align="center" border="0" cellpadding="5">

<tr>
  <td align="right">
<math>
\frac{c}{a }
</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
(1 - e^2)^{1 / 2}
\, ;
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>
\Omega^2 \equiv \frac{\omega_0^2}{4\pi G \rho }
</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>\frac{1}{2e^2} \biggl[
(3-2e^2)(1-e^2)^{1 / 2} \cdot \frac{\sin^{-1}e}{e} 
- 3(1-e^2)
\biggr]
\, ;
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\tau \equiv \frac{T_\mathrm{rot}}{|W_\mathrm{grav}|}</math>
  </td>
  <td align="center">
<math>= </math>
  </td>
  <td align="left">
<math>~ \frac{3}{2e^2}\biggl[ 1 - \frac{e(1-e^2)^{1 / 2}}{\sin^{-1} e}\biggr] - 1 \, ;</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>L_*^2 \equiv \frac{L^2}{(GM^3\bar{a})}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>\frac{6}{5^2} \biggl[ (3-2e^2)(1-e^2)^{1 / 2} \cdot \frac{\sin^{-1}e}{e^3} - \frac{3(1-e^2)}{e^2}\biggr](1 - e^2)^{-2 / 3}  \, .</math>
  </td>
</tr>
</table>
 </td>
</tr>
</table>

====Second-Harmonic Distortions====

<table border="0" cellpadding="3" align="center" width="80%">
<tr><td align="left">
<font color="darkgreen">
"The disturbances of the Maclaurin spheroids that belong to the second harmonics lead to instability for smaller values of the angular momentum than do those that belong to any higher harmonic &hellip; [Of] the five oscillation frequencies associated with these disturbances &hellip; one is able to distinguish two critical points &hellip;  The first, when <math>e</math> is 0.8127</font> &#8212; our <font color="red">Model A</font> &#8212; <font color="darkgreen">is the point of bifurcation where the Maclaurin and Jacobi sequences intersect.  Although &sigma;<sup>2</sup> vanishes at this critical point, it is positive on either side, so that instability does not set in here, at least in the absence of further effects.  The second, when <math>e</math> is 0.9529 </font>&#8212; our <font color="red">Model B</font> &#8212; <font color="darkgreen"> represents the onset of </font> [dynamical]<font color="darkgreen"> instability, and also the most flattened Maclaurin spheroid that can lie on a Riemann sequence of type S &hellip;</font>"
</font>
</td></tr>
<tr><td align="right">
— Drawn from pp. 475 - 476 of {{ Lebovitz67_XXXIV }} 
</td></tr></table>

<font color="red"><b>Model A</b></font>:
<table border="0" align="center" width="95%"><tr><td align="left">
This is the point along the Maclaurin spheroid sequence where the Jacobi sequence bifurcates.  Some of the quantitative characteristics of this critical axisymmetric configuration are identified in Table IV (Chapter 6, &sect;39, p. 103) of [<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>].  See also &hellip; the first line in Table B1 (p. 446) of {{ Bardeen71 }}; and Appendices D.3 & D.4 (pp. 485-486) of [<b>[[Appendix/References#T78|<font color="red">T78</font>]]</b>].
</td></tr></table>

<span id="RRSTEMfigure1">&nbsp;</span>
<table border="1" align="center" cellpadding="5" width="70%">
<tr><th align="center" colspan="2">RRSTEM Figure 1</th></tr>
<tr>
  <td align="center">
<b>EFE Diagram</b><br />

[[File:OurEFEannotated2.png|300px|OurEFE]]
  </td>
  <td align="center">
<b>&Omega;<sup>2</sup> vs. j<sup>2</sup> Diagram</b><br />

[[File:OurHE84Fig1annotated2.png|425px|OurHE84Fig1]]
  </td>
</tr>

<tr>
  <td align="left" colspan="2">
''Left panel:''  This version of the [[ThreeDimensionalConfigurations/RiemannStype#Summary|familiar EFE Diagram]] displays the horn-shaped region &#8212; bounded by the lower (LSA) and upper (USA) self-adjoint sequences &#8212; where all equilibrium, S-type Riemann ellipsoids reside.  Axisymmetric <math>(b/a = 1)</math>, uniformly rotating equilibrium configurations belong to the Maclaurin spheroid sequence which, as indicated, functions as the right-hand vertical boundary of the EFE diagram. Maclaurin spheroids that have an eccentricity less than that of <font color="red">Model A</font> <math>(e = 0.812670)</math> lie along the ''blue'' segment of the sequence; the ''orange'' segment is populated by Maclaurin spheroids with eccentricities larger than that of Model A but less than that of <font color="red">Model B</font> <math>(e = 0.95289)</math>; all other Maclaurin spheroids <math>(0.95289 < e \le 1)</math> lie along the ''black'' segment.

''Right panel:''  The solid, multi-colored curve shows how the (square of the) dimensionless rotation frequency, <math>\Omega^2</math>, varies with the (square of the) dimensionless total angular momentum, <math>j^2</math>, along the Maclaurin spheroid sequence.  As in the accompanying  (''left panel'') EFE diagram, <font color="red">Model A</font> and <font color="red">Model B</font> mark the ends of the differently colored curve segments.  

''Both panels:''  A pair of small yellow circular markers identify the points where, respectively, the USA and LSA sequences intersect the Maclaurin spheroid sequence; and the small green square marker identifies where <math>\Omega^2</math> has its maximum value along the Maclaurin spheroid sequence. The sequence of uniformly rotating Jacobi ellipsoids is identified by the series of small solid purple markers; <font color="red">Model A</font> is the axisymmetric equilibrium configuration at the point where the Jacobi ellipsoid sequence intersect (bifurcates from) the Maclaurin spheroid sequence.  Similarly, <font color="red">Model B</font> lies at the intersection of the LSA with the Maclaurin spheroid sequence.
  </td>
</tr>
</table>

<font color="red"><b>Model B</b></font>:
<table border="0" align="center" width="95%"><tr><td align="left">
From p. 141 of [<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>], we find that this model, <font color="darkgreen">"&hellip; the Maclaurin spheroid on the verge of dynamical instability, is the first member of the</font> [Riemann S-type ellipsoid] <font color="darkgreen"> self-adjoint sequence <math>x = + 1</math>."</font>  Some of the quantitative characteristics of this critical (axisymmetric) Maclaurin spheroid are identified in Table VI (Chapter 7, &sect;48, p. 142) of [<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>]; for example, <math>c/a = a_3/a_2 = 0.30333 \Rightarrow e = 0.95289</math>.  In Table B1 (p. 446) of {{ Bardeen71 }}, this is referred to as the "First nonaxisymmetric dynamical instability."
</td></tr></table>

<table border="0" cellpadding="3" align="center" width="80%">
<tr><td align="left">
<font color="darkgreen">
"&hellip; to what kind of configuration does [this] instability lead?  {{ Rossner67_XXXVIII }} has answered this question with the aid of Riemann's formulation of the nonlinear, ordinary differential equations describing the motion of a liquid ellipsoid.  Integrating the equations numerically, Rossner found that the configuration neither gets disrupted nor finds its way to a new steady state, but performs a complicated unsteady motion."</font>  &#8212; See [[#Finite-Amplitude_Oscillations_of_the_Maclaurin_Spheroid|further reference to Rossner's work, below]].
</td></tr>
<tr><td align="right">
— Drawn from pp. 475 - 476 of {{ Lebovitz67_XXXIV }} 
</td></tr></table>

====Ring-Like Distortions====


<font color="red"><b>Model C</b></font>: &nbsp; ([[Appendix/Ramblings/MacSphCriticalPoints#One-Ring_(Dyson-Wong)_Sequence|additional supporting discussion]])
<table border="0" align="center" width="95%"><tr><td align="left">
{{ ES81 }} claim that the one-ring (Dyson-Wong toroid) sequence bifurcates from the Maclaurin sequence precisely at the point where the spheroid has an eccentricity, <math>e = e_\mathrm{cr} = 0.98523</math> &#8212; in which case, also, <math>\Omega^2 = 0.08726</math> and <math>j^2 = 0.02174</math>.  In support of this conjecture, they point out that, {{ Chandrasekhar67_XXXfull }} &#8212; and {{ Bardeen71 }} have shown that this is <font color="darkgreen">&hellip; a neutral point on the Maclaurin sequence against the perturbation of <math>P_4(\eta)</math> displacement at the surface where <math>\eta</math> is one of the spheroidal coordinates."</font>  This is also the "neutral point" on the Maclaurin sequence labeled "F" in Table I of {{ HE82 }}; and the "bifurcation point" along the Maclaurin sequence that is labeled by the quantum numbers, <math>(n, m) = (4, 0)</math> in Table 1 of {{ HE84 }} as well as in the inset box of the ''left panel'' of our [[#RRSTEMfigure2|RRSTEM Figure 2]] (immediately below).

See also &hellip; {{ AKM2003full }}
</td></tr></table>
<span id="RRSTEMfigure2">&nbsp;</span>
<table border="1" align="center" cellpadding="5">
<tr><th align="center" colspan="2">RRSTEM Figure 2</th></tr>
<tr>
  <td align="center" colspan="1" rowspan="3">
[[File:ES81OneRingWithInsetBox3.png|600px|One-Ring Sequence]]
  </td>
  <td align="center" colspan="1" rowspan="1">
Figure 1 extracted from &sect;2.2, p. 488 of &hellip;<br />{{ CKST95bfigure }}
  </td>
</tr>
<tr>
  <td align="center" colspan="1">[[File:CKST95bFig1annotated3.png|300px|CKST95b Figure 1]]</td>
</tr>
<tr>
  <td align="center">An analogous illustration appears as Figure 1 (p. 585) of &hellip;<br />{{ HE83figure }}</td>
</tr>

<tr>
  <td align="left" colspan="2">
''Left panel (primary plot):''  As in the ''right panel'' of [[#RRSTEMfigure1|RRSTEM Figure 1, above]], the solid, multi-colored curve shows how <math>\Omega^2</math> varies with <math>0 \le j^2 \le 0.04</math> along the Maclaurin spheroid sequence.  Models A and B are not explicitly labeled, but the plot still shows the small solid circular markers (purple and yellow, respectively) that identify the locations of these models on the spheroid sequence.  The point where the Jacobi sequence bifurcates from the Maclaurin spheroid sequence (Model A) is labeled by the quantum numbers, <math>(n, m) = (2, 2)</math>, to indicate what geometric distortion, <math>P_n^m</math>, that is associated with this particular bifurcation.  Drawing from Table 1 of {{ HE84 }}, five other neutral points are marked by red crosses; three carry labels in this ''primary plot'' &#8212; <math>(3, 3), (4, 4), (3, 1)</math> &#8212; and the remaining two are labeled in the ''inset box'' &#8212; <math>(4, 2), (4, 0)</math>.

''Left panel (inset box):''  An ''inset box'' is used to magnify the segment of the Maclaurin spheroid sequence where bifurcation to the so-called one-ring (Dyson-Wong toroid) sequence occurs.  The neutral point that is believed to be associated with the bifurcation point, itself, carries the geometric distortion label, <math>(n, m) = (4, 0)</math>, and is identified as our <font color="red">Model C</font>.  As has been detailed in [[Appendix/Ramblings/MacSphCriticalPoints#Model-Sequence_Details|our separate chapter discussion]], the smooth (pink) curve that connects the spheroid sequence to the one-ring sequence has been defined by the set of 18 equilibrium models presented by {{ ES81 }}; the set of small green square markers identify nine equilibrium models obtained from Table I of the separate study by {{ HES82 }}; and the small purple triangular markers identify eight equilibrium models obtained from Table Ia of {{ Hachisu86a }}.

''Right panel:''  Figure 1 (plus caption) from {{ CKST95b }} has been reprinted here to emphasize its similarity to, and overlap with our ''inset box''.  According to the caption of this reprinted figure, the filled circular marker labeled "A" identifies the bifurcation point on the Maclaurin spheroid sequence, where <math>e = 0.985226</math>. Accordingly, we have annotated the reprinted figure to indicate that the axisymmetric equilibrium model associated with point "A" is exactly our <font color="red">Model C</font>.  As is stated in the caption of this reprinted figure, the dotted line XBC denotes the (hypothesized) onset of a secular instability that &#8212; in the nonlinear regime and conserving total angular momentum (vertical dotted line) &#8212; should deform the spheroid into a ring-like configuration.  
  </td>
</tr>
</table>


<font color="red"><b>Model D</b></font>: 
<table border="0" align="center" width="95%"><tr><td align="left">
First ''dynamical'' ring mode instability and bifurcation point to the Maclaurin toroid sequence. Also identified at <math>(e, j^2) = (0.998556, 0.04305)</math> as a <math>P_6</math> bifurcation point in Table 2 (p. 292) of {{ EH85 }}.
</td></tr></table>
<span id="RRSTEMfigure3">&nbsp;</span>
<table border="1" align="center" cellpadding="5">
<tr><th align="center" colspan="2">RRSTEM Figure 3</th></tr>
<tr>
  <td align="center" colspan="1" rowspan="2">
[[File:EH85MacToroidWithInsetBox2.png|600px|Maclaurin Toroid Sequence]]
  </td>
  <td align="center" colspan="1" rowspan="1">
Figure 6 extracted from &sect;4.2, p. 493 of &hellip;<br />{{ CKST95bfigure }}
  </td>
</tr>
<tr>
  <td align="center" colspan="1">[[File:CKST95bFig6annotated2.png|400px|CKST95b Figure 1]]</td>
</tr>

<tr>
  <td align="left" colspan="2">
''Left panel (primary plot):''  The solid, multi-colored curve shows how <math>\tau \equiv T_\mathrm{rot}/|W_\mathrm{grav}|</math> (rather than <math>\Omega^2</math>) varies with <math>0 \le j^2 \le 0.07</math> along the Maclaurin spheroid sequence.  The positions of Models A, B, C and, now also, <font color="red">Model D</font> along this sequence are labeled. The colors of the curve segments and the positions of the (yellow and purple) small circular markers have the same meanings as in [[#RRSTEMfigure1|RRSTEM Figure 1]].  The (pink curve) one-ring sequence from [[#RRSTEMfigure2|RRSTEM Figure 2]] is also displayed.

''Left panel (inset box):''  An ''inset box'' is used to magnify the segment of the Maclaurin spheroid sequence where bifurcation to the so-called Maclaurin toroid sequence occurs.  Three neutral points along the Maclaurin spheroid sequence are identified by light-green circular markers, and are labeled according to their respective geometric distortions <math>(P_6, P_8, P_4)</math> &#8212; see Table 2 of {{ EH85 }}.  They have argued that the neutral point that is associated with bifurcation, itself, carries the geometric distortion label, <math>P_6</math>; it has been identified here as our <font color="red">Model D</font>.  As has been detailed in [[Apps/MaclaurinToroid#Maclaurin_Toroid_(EH85)|our separate chapter discussion]], the smooth (violet) curve that connects the spheroid sequence to the Maclaurin toroid sequence has been defined by the set of 15 equilibrium models presented by {{ EH85 }} in their Table 2.  [This data supersedes the modeling of the Maclaurin toroid sequence presented by {{ MPT77 }} in their discovery paper.] 

''Right panel:''  Figure 6 (plus caption) from {{ CKST95b }} has been reprinted here to emphasize its similarity to, and overlap with our ''inset box''.  According to the caption of this reprinted figure, the filled circular marker labeled "A" identifies the point along the Maclaurin spheroid sequence where <math>e = 0.985226</math>; according to {{ Bardeen71 }}, this neutral point is associated with a <font color="red">P<sub>4</sub></font> geometric distortion and should be associated with the onset of dynamical axisymmetric instability. {{ EH85 }} agree that this is the eccentricity at which the P<sub>4</sub> neutral point resides, but they argue that bifurcation &#8212; and onset of dynamical axisymmetric instability &#8212; should be associated instead with the <font color="red">P<sub>6</sub></font> neutral point, as labeled in our ''inset box'', because "P<sub>6</sub>" is encountered ''earlier'' than "P<sub>4</sub>" along the Maclaurin spheroid sequence.  [As they argue, {{ Bardeen71 }} misidentified the bifurcation point because he only investigated models undergoing a P<sub>4</sub> distortion.]  In accordance with the arguments of {{ EH85 }}, our choice for <font color="red">Model D</font> is the Maclaurin spheroid whose eccentricity is the same as the P<sub>6</sub> neutral point, that is, <math>e = 0.998556</math>.
  </td>
</tr>
</table>

====Finite-Amplitude Oscillations of the Maclaurin Spheroid====

In &sect;53 (pp. 172 - 184) of [<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>] we find a discussion titled, "''A class of finite-amplitude oscillations of the Maclaurin spheroid''".  The primary focus is on simulations presented by {{ Rossner67_XXXVIIIfull }} &#8212; which may provide excellent points of comparison for our own investigations into the oscillatory motions of initially stable &#8212; but perturbed &#8212; Maclaurin spheroids,  and the nonlinear evolution of initially unstable Maclaurin spheroids.


===Jacobi Ellipsoids===
<span id="RRSTEMtable2">&nbsp;</span>
<table border="1" align="center" cellpadding="5" width="85%">
<tr>
  <td align="center" colspan="9">
<b>RRSTEM Table 2</b><br />
Established Critical-Point Models Along the Jacobi Ellipsoid Sequence
  </td>
</tr>
<tr>
  <td align="center" rowspan="2">
Model
  </td>
  <td align="center" rowspan="2">
<math>\frac{b}{a}</math>
  </td>
  <td align="center" rowspan="2">
<math>\frac{c}{a}</math>
  </td>
  <td align="center" rowspan="2">
<math>\Omega^2</math>
  </td>
  <td align="center" rowspan="2">
<math>\tau</math>
  </td>
  <td align="center" rowspan="2">
<math>j^2 = \frac{1}{3}\biggl(\frac{4\pi}{3}\biggr)^{-4 / 3} L_*^2</math>
  </td>
  <td align="center" colspan="3">
Bifurcation Characteristics &hellip;
  </td>
</tr>

<tr>
  <td align="center"><sup>&dagger;</sup>Geometric Distortion</td>
  <td align="center">Angular Mom. Profile</td>
  <td align="center">Instability Type</td>
</tr>

<tr>
  <td align="center" rowspan="1">
<font color="red"><b>E</b></font>
  </td>
  <td align="center"><math>0.2972</math></td>
  <td align="center"><math>0.2575</math></td>
  <td align="center"><math>0.053286</math></td>
  <td align="center"><math>0.1863</math></td>
  <td align="center"><math>1.1507\times 10^{-2}</math></td>
  <td align="center">Dumbbell-shaped</td>
  <td align="center">Uniform <math>\omega_0</math></td>
  <td align="center">Secular</td>
</tr>

<tr>
  <td align="center" rowspan="1">
<font color="red"><b>F</b></font>
  </td>
  <td align="center"><math>0.432232</math></td>
  <td align="center"><math>0.345069</math></td>
  <td align="center"><math>0.07101</math></td>
  <td align="center"><math>0.1628</math></td>
  <td align="center"><math>7.491\times 10^{-3}</math></td>
  <td align="center">Pear-shaped</td>
  <td align="center">Uniform <math>\omega_0</math></td>
  <td align="center">Secular</td>
</tr>

<tr>
  <td align="left" colspan="9">
Given the value of the square of the angular velocity, <math>\Omega^2 \equiv \omega_0^2/(4\pi G \rho)</math>, we immediately know from the [[ThreeDimensionalConfigurations/JacobiEllipsoids#Angular_Momentum_Constraint|accompanying detailed discussion]] that &hellip;

<table align="center" border="0" cellpadding="5">

<tr>
  <td align="right"><math>L_*^2 \equiv \frac{L^2}{GM^3 \bar{a}}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>\frac{3 \Omega^2}{5^2} \cdot \frac{(a^2 + b^2)^2}{\bar{a}^4} \, .</math>
  </td>
</tr>
<tr><td align="center" colspan="3">[<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>], <font color="#00CC00">&sect;39, p. 103, Eq. (16)</font> <br /><font size="-1">(Note different definition of <math>\Omega^2</math> in EFE)</font></td></tr>
</table>
 </td>
</tr>
</table>

====Pear-Shaped Distortion====

<font color="red"><b>Model F</b></font>:
<table border="0" align="center" width="95%"><tr><td align="left">

According to <font color="#00CC00">Chapter 6, &sect;40</font> of [<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>], a pear-shaped configuration bifurcates from the Jacobi Sequence at <math>(b/a, c/a) = (0.432232, 0.345069)</math>, where <math>\omega_0^2/(\pi G \rho) = 0.284030</math>;  this result has been obtained via the virial relations &#8212; see Eq. (28) on p. 106 &#8212; as well as via a direct perturbation analysis &#8212; see Eq. (57) on p. 110.  From the information provided in the first row of Table I in {{ EHS82 }}, we appreciate as well that <math>\omega_0^2/(4\pi G \rho) = 0.07101; ~j^2 = 0.07821; ~\tau \equiv T/|W| = 0.1628</math>. (This information also has been recorded near the end of our [[ThreeDimensionalConfigurations/JacobiEllipsoids#Bifurcation_from_Maclaurin_to_Jacobi_Sequence|accompanying discussion of the Jacobi ellipsoid sequence]].)
</td></tr></table>

<span id="RRSTEMfigure4">&nbsp;</span>
<table border="1" align="center" cellpadding="5" width="80%">
<tr><th align="center" colspan="2">RRSTEM Figure 4</th></tr>
<tr>
  <td align="center" colspan="1" rowspan="1">
[[File:PearAndDumbbellModelF2.png|350px|Pear-Shaped Sequence]]
  </td>
  <td align="left" colspan="1" rowspan="1">
''Primary plot:''  As in the ''right panel'' of [[#RRSTEMfigure1|RRSTEM Figure 1, above]], the solid, multi-colored curve shows how <math>\Omega^2</math> varies with <math>0 \le j^2 \le 0.04</math> along the Maclaurin spheroid sequence; for reference, <font color="red">Model A</font> and <font color="red">Model B</font> are labeled. Also as in [[#RRSTEMfigure1|RRSTEM Figure 1]], starting at the Model A bifurcation point, the purple dashed curve identifies the equilibrium sequence of Jacobi ellipsoids. A red cross labeled <font color="red">F</font> identifies the position along the Jacobi ellipsoid sequence that is a "neutral point against the 3<sup>rd</sup>-order harmonic perturbation."  A so-called pear-shaped sequence branches off at this point; here, we have displayed the behavior of this sequence by drawing the properties of equilibrium models from Table I (p. 1071) of {{ EHS82 }}. 

''Inset box:'' Especially Because the pear-shaped sequence is  very short, an ''inset box'' is used to magnify the plotted sequence; this was also done by {{ EHS82 }} &#8212; see their Figure 1 (p. 1072).  A red arrow, that accompanies our <font color="red">Model F</font> label, points to the location along the Jacobi sequence where the relevant neutral point lies.  As is suggested by {{ EHS82 }}, the presumption is that this pear-shaped sequence bifurcates from the Jacobi sequence at the neutral point.  But this is not the case in our plot.  Upon closer inspection, we see that the pear-shaped sequence ''does'' appear to bifurcate from the Jacobi sequence in Figure 1 of {{ EHS82 }} and that this happens because the Jacobi sequence is shifted a bit from ''our'' depiction of the Jacobi sequence location.  Curious!
  </td>
</tr>

</table>

====Dumbbell-Shaped Distortion====

<font color="red"><b>Model E</b></font>:

<ul><li>
According to the last pair of equations on p. 128 of [<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>], <font color="#00CC00">Chapter 6, &sect;45</font>, a neutral point belonging to the fourth harmonic (a dumbbell-shaped) distortion arises on the Jacobi Sequence at <math>(b/a) = 0.2972</math> and <math>\cos^{-1}(c/a) = 75\overset{\circ}{.}081~~~~\Rightarrow ~~~~ (c/a) = 0.2575</math>.  Chronologically, this result for <math>(b/a, c/a)</math> appears first in Eq. (93) on p. 635 of {{ Chandrasekhar67_XXXIIfull }}.  Then, in Eq. (66) on p. 302 of {{ Chandrasekhar68_XXXVfull }} &#8212; we find <math>\cos^{-1}(c/a) = 75\overset{\circ}{.}068</math>, along with a footnote [5] which states, <font color="darkgreen">"The value <math>\cos^{-1}(c/a) = 75\overset{\circ}{.}081</math> found earlier differs slightly; but the difference is not outside the limits of accuracy of the numerical evaluation."</font>
</li>
<li>
According to the first row of properties in Table II of {{ EHS82 }}, we find that <font color="red">Model E</font> is characterized by the properties &hellip;  <math>\Omega^2/(4\pi G \rho) = 0.0532</math>; <math>j^2 = ( 3\cdot 2^{-8} \pi^{-4} )^{1/3}  L^2/(GM^3\bar{a}) = 0.01157</math>; and <math>\tau \equiv T/|W| = 0.1863 </math>.  &nbsp;I have not (yet) found the corresponding value of <math>\Omega^2</math> in any of Chandrasekhar's publications, but if we combine the value of <math>\Omega^2</math> obtained from {{ EHS82hereafter }} with the values of <math>(b/a, c/a)</math> obtained from [<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>], we find &hellip; <math>\Omega^2/(\pi G \rho) = 0.2128</math>; from the [[#Angular_Momentum_Constraint|above expression]], <math>L_* = 0.48242</math>; and &nbsp; <math>j^2 = ( 3\cdot 2^{-8} \pi^{-4})^{1/3}  L_*^2 = 0.01149</math>.  &nbsp; This value of <math>j^2</math> is very close to the value obtained by {{ EHS82hereafter }}.</li>
<li>
In the paragraph at the top of the right-hand column of p. 467 of {{Hachisu86bfull }}, we find &hellip; <math>\Omega^2/(4\pi G\rho) = 0.0535</math>; <math>j^2 = 0.01157</math>.</li>
<li>
{{ CKST95bfull }} grab parameter values from a variety of sources.  In subsection "B" (''Jacobi Ellipsoid to Binary'') of their Table 1 (p. 494) and in the first paragraph of their &sect;3.2 (p. 492), they state &hellip; <math>(b/a, c/a) = (0.29720, 0.25746)</math>; &nbsp;<math>\Omega^2/(4\pi G \rho) = 0.0532790</math>; &nbsp; and, <math>j^2 = 0.0115082</math>.</li>
<li><font color="red">NOTE (23 May 2023):</font> &nbsp; Plugging the axis values from p. 128 of [<b>[[Appendix/References#EFE|<font color="red">EFE</font>]]</b>]  &#8212; that is, <math>(b/a, c/a) = (0.2972, 0.2575)</math> &#8212; into our "Riemann01.for" application, we find <math>(A_1, A_2, A_3) = 0.171772322973, 0.844742744895, 0.983484932132)</math>, and <math>\Omega^2/(4\pi G \rho) = 0.053286</math>, and <math>j^2 = 0.011507</math>.</li>
</ul>


<span id="RRSTEMfigure5">&nbsp;</span>
<table border="1" align="center" cellpadding="5">
<tr><th align="center" colspan="2">RRSTEM Figure 5</th></tr>
<tr>
  <td align="center" colspan="1" rowspan="2">
[[File:PearAndDumbbellModelE4.png|600px|Pear and Dumbbell Sequences]]
  </td>
  <td align="center" colspan="1" rowspan="1">
Figure 4 extracted from &sect;3.2, p. 493 of &hellip;<br />{{ CKST95bfigure }}
  </td>
</tr>
<tr>
  <td align="center" colspan="1">[[File:CKST95bFig4annotatedBetter.png|300px|CKST95b Figure 4]]</td>
</tr>

<tr>
  <td align="left" colspan="2">
''Left panel (primary plot):''  Same as the ''primary plot'' displayed in [[#RRSTEMfigure4|RRSTEM Figure 4, immediately above]].  A red cross labeled <font color="red">E</font> identifies the position along the Jacobi ellipsoid sequence that is a neutral point against a 4<sup>th</sup>-order harmonic perturbation.  A so-called dumbbell-shaped sequence branches off at this point; it, in turn, transitions to a sequence of equal-mass binaries.

''Left panel (inset box):''  An ''inset box'' shows more clearly the sequence of equilibrium models that make up the dumbbell (green markers and curve) and binary (blue markers and curve) sequences.  Here, the dumbbell sequence is defined by a set of equilibrium models drawn from Table II (p. 1073) of {{ EHS82 }} while the binary sequence is defined by a set of equilibrium models drawn from the subsection (p. 243) of Table 1 labeled "N = 0" in {{ HE84b }}.

''Right panel:''  Figure 4 (plus caption) from {{ CKST95b }} has been reprinted here to emphasize its similarity to, and overlap with our ''inset box''.  According to the caption of this reprinted figure, the filled circular marker labeled "A" identifies the bifurcation point on the Jacobi ellipsoid sequence, where <math>(b/a, c/a) = (0.2972, 0.2575)</math>.  Accordingly, we have annotated the reprinted figure to indicate that the Jacobi ellipsoid model associated with point "A" is exactly our <font color="red">Model E</font>.  As is stated in the caption of this reprinted figure, the dotted line XBC denotes the (hypothesized) onset of a secular instability that &#8212; in the nonlinear regime and conserving total angular momentum (vertical dotted line) &#8212; should lead to fission of the ellipsoid into an equal-mass binary system.  
  </td>
</tr>
</table>

==Details==

<ul>
  <li>Maclaurin Spheroid Sequence:
  <ol>
    <li>[[Apps/MaclaurinSpheroidSequence#Maclaurin_Spheroid_Sequence|Standard Presentation]]</li>
    <li>[[Apps/MaclaurinSpheroidSequence#Alternate_Sequence_Diagrams|Alternate Sequence Diagrams]]</li>
    <li>[[Apps/MaclaurinSpheroidSequence#Oblate_Spheroidal_Coordinates|Oblate-Spheroidal Coordinates]]</li>
    <li>Various ''[[AxisymmetricConfigurations/SolutionStrategies#SRPtable|Simple Rotation Profiles]]''
    <li>[[Apps/MaclaurinSpheroidSequence#Bifurcation_Points_Along_Maclaurin-Spheroid_Sequence|Two Particularly Relevant Angular Momentum Distributions]]
      <ul type="-">
        <li>Uniform Rotation &nbsp; &rarr; &nbsp; <math>(n, m) = (4, 0)</math> bifurcation leads to Dyson-Wong "one-ring" sequence</li>
        <li><math>n' = 0</math> &nbsp; &rarr; &nbsp; <math>P_6</math> bifurcation leads to {{ MPT77 }} "Maclaurin Toroid" sequence</li>
      </ul>
    </li>
  </ol>
  </li>
  <li>[[Apps/EriguchiHachisu/Models|Eriguchi, Hachisu, and their various Colleagues]]</li>
  <li>[[Apps/MaclaurinToroid|Maclaurin Toroid Sequence (MPT77 &amp; EH85)]]</li>
  <li>[[Appendix/Ramblings/MacSphCriticalPoints#PhaseTransition|Phase Transition]] to Dyson-Wong (one-ring) sequence</li>
</ul>

<table border="1" cellpadding="8" align="center" width="100%"><tr><td align="center" bgcolor="lightblue">
See also our associated, more detailed discussion of [[Appendix/Ramblings/MacSphCriticalPoints|Critical Points along the Maclaurin Spheroid Sequence]].
</td></tr></table>

=Lebovitz and Lifschitz=
==Invitation==

<table border="1" align="center" width="65%" cellpadding="8">
<tr>
<td align="left">
On 8 May 2023, Stephen Sorokanich sent the following email to Joel:

"<font color="darkgreen">Dr. Cohl and I are wondering if you would like to join our discussion of &hellip;"</font>
<div align="center">{{ LL96bfigure }}</div>

"<font color="darkgreen">This paper is very interesting to us because the authors consider Eulerian perturbations to the steady-state S-type ellipsoid velocity and pressure fields. The resulting linear equations in the Eulerian frame seem ideal for us to compare with our current finite element simulations of the nonlinear system. This is the only paper I've found by Lebovitz and Lifschitz where the stability analysis is carried out in an Eulerian frame. 

There's some nice analysis in the paper too. The linearized equations are solved using the "geometrical-optics approximation," which reduces the 3-dimensional linearized Euler equations to a system of ordinary differential equations. These solutions are used to draw conclusions for the stability of the S-type ellipsoids. Maybe there are other approximations to this system?</font>"

NOTE from Joel:  The referenced paper is obviously related to the companion publication &hellip;
<div align="center">{{ LL96figure }}</div>We have already prepared a [[ThreeDimensionalConfigurations/Stability/RiemannEllipsoids#Euler_Equation|MediaWiki chapter]] whose focus is on this "{{ LL96hereafter }}" publication.

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

==WKB (geometrical-optics) Method==

As Stephen has pointed out, in {{ LL96b }}, the linearized equations are solved using the "geometrical-optics approximation."  I am not familiar with this so-called "geometrical-optics" method, but in the last couple of paragraphs of &sect;1 of their paper, {{ LL96bhereafter }} refer to it also as the "WKB" approximation/method.  I ''have'' been exposed to this alternate terminology; the WKB approximation has been used by astronomers in connection with efforts to understand the onset and development of "spiral-arm" structures in disk galaxies.  Here is an excerpt from [<b>[[Appendix/References#BT87|<font color="red">BT87</font>]]</b>]. 

<table border="1" align="center" width="80%" cellpadding="8">
<tr><td align="left">
<div align="center">'''The Tight-Winding Approximation'''<br />&sect;6.2.2 (p. 352) of [<b>[[Appendix/References#BT87|<font color="red">BT87</font>]]</b>] </div>
<font color="darkgreen">"&hellip; In the early 1960s a number of workers, notably A. J. Kalnajs, C. C. Lin, and A. Toomre, realized that for tightly wound density waves (''i.e.,'' waves whose radial wavelength is much less than the radius</font> [of the galaxy's disk]) <font color="darkgreen">the long-range coupling</font> [due to gravity] <font color="darkgreen">is negligible, the response is determined locally, and the relevant solutions are analytic.  As we shall see, this '''tight-winding''' or '''WKB approximation'''<sup>&dagger;</sup> is an indispensable tool for understanding the origin and evolution of density waves in galaxies."

----

<sup>&dagger;</sup>Named  after the closely related [https://en.wikipedia.org/wiki/WKB_approximation Wentzel-Kramers-Brillouin approximation] of quantum mechanics.</font>
</td></tr></table>

==Governing Linearized Equations==
In our [[ThreeDimensionalConfigurations/Stability/RiemannEllipsoids#EulerRotating|accompanying discussion]] we present a derivation of the &hellip;

<table border="1" align="center" cellpadding="8" width="80%">
<tr><td align="center" bgcolor="lightgreen">{{ Lebovitz89afigure }}</td></tr>
<tr><td align="left">
<div align="center">
<font color="#770000">'''Mixed Lagrangian/Eulerian Representation'''</font><br />
of the Euler Equation <br />
<font color="#770000">'''as viewed from a Rotating Reference Frame'''</font>
</div>
<table border="0" cellpadding="3" align="center">
<tr>
  <td align="right">
<math>\frac{D\mathbf{u}}{Dt} + 2\boldsymbol\omega \boldsymbol\times \mathbf{u}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> 
-~ \rho^{-1} \nabla p + \mathbf\nabla \{ \Phi_\mathrm{L89} + \tfrac{1}{2} |\boldsymbol\omega \boldsymbol\times \mathbf{x}|^2 \}
\, .
</math>
  </td>
</tr>
<tr>
  <td align="center" colspan="3">
{{ Lebovitz89a }}, §2, p. 223, Eq. (2)
  </td>
</tr>
</table>
</td></tr>
</table>

When rewritten with an alternate set of variable notations &#8212; specifically, <math>(\mathbf{u}, p, \Phi_\mathrm{L96}) \rightarrow (\mathbf{V}, P, \Psi)</math> &#8212; and recognizing that <math>D/Dt = (\partial/\partial t + \mathbf{V}\cdot \nabla)</math>, this key relation becomes,
<table border="1" align="center" cellpadding="8" width="80%">
<tr><td align="center" bgcolor="yellow">{{ LL96bfigure }}</td></tr>
<tr><td align="left">
<table border="0" cellpadding="3" align="center">
<tr>
  <td align="right">
<math>\biggl[\frac{\partial}{\partial t} + \mathbf{V}\cdot \nabla\biggr]\mathbf{V} + 2\boldsymbol\omega \boldsymbol\times \mathbf{V}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> 
-~ \rho^{-1} \nabla P + \mathbf\nabla \{ \Psi + \tfrac{1}{2} |\boldsymbol\omega \boldsymbol\times \mathbf{x}|^2 \}
\, .
</math>
  </td>
</tr>
<tr>
  <td align="center" colspan="3">
{{ LL96b }}, &sect;2, p. 929, Eq. (2.1)
  </td>
</tr>
</table>
</td></tr>
</table>

Next, following {{ LL96bhereafter }}, we let each variable be represented by the sum of an unperturbed, steady-state component (subscript zero) and a small perturbation,
<table border="0" align="center" cellpadding="8">
<tr>
  <td align="right"><math>\mathbf{V}</math></td>
  <td align="center"><math>\rightarrow</math></td>
  <td align="left"><math>\mathbf{V}_0 + a_0(\pi G \rho_0)^{1 / 2}\mathbf{u} \, ,</math></td>
</tr>

<tr>
  <td align="right"><math>\rho</math></td>
  <td align="center"><math>\rightarrow</math></td>
  <td align="left"><math>\rho_0 </math> &nbsp; &nbsp; &nbsp; (''i.e.,'' remains unchanged)</td>
</tr>

<tr>
  <td align="right"><math>P</math></td>
  <td align="center"><math>\rightarrow</math></td>
  <td align="left"><math>P_0 + (\pi G \rho_0^2 a_0^2) p \, ,</math></td>
</tr>

<tr>
  <td align="right"><math>\Psi</math></td>
  <td align="center"><math>\rightarrow</math></td>
  <td align="left"><math>\Psi_0 + (\pi G a_0^2 \rho_0)\psi \, ,</math></td>
</tr>
</table>

where <math>a_0</math> is a characteristic length scale.  As a result, we have,
<table border="0" cellpadding="3" align="center">
<tr>
  <td align="right">
<math>\biggl[ \frac{\partial}{\partial t} + [\mathbf{V}_0 + a_0(\pi G \rho_0)^{1 / 2}\mathbf{u}]\cdot 
\nabla\biggr][\mathbf{V}_0 + a_0(\pi G \rho_0)^{1 / 2}\mathbf{u}] 
+ 2\boldsymbol\omega \boldsymbol\times [\mathbf{V}_0 + a_0(\pi G \rho_0)^{1 / 2}\mathbf{u}]</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> 
-~ \rho_0^{-1} \nabla [P_0 + (\pi G \rho_0^2 a_0^2)p] + \mathbf\nabla \{ [\Psi_0 + (\pi G \rho_0^2 a_0)\psi ] + \tfrac{1}{2} |\boldsymbol\omega \boldsymbol\times \mathbf{x}|^2 \}
\, .
</math>
  </td>
</tr>
</table>

Then, subtracting off the unperturbed steady-state component, namely,
<table border="0" cellpadding="3" align="center">
<tr>
  <td align="right">
<math>(\mathbf{V}_0 \cdot \nabla) \mathbf{V}_0 + 2\boldsymbol\omega \boldsymbol\times \mathbf{V}_0</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> 
-~ \rho_0^{-1} \nabla P_0 + \mathbf\nabla \{ \Psi_0 + \tfrac{1}{2} |\boldsymbol\omega \boldsymbol\times \mathbf{x}|^2 \}
\, ,
</math>
  </td>
</tr>
</table>
we obtain,
<table border="0" cellpadding="3" align="center">

<tr>
  <td align="right">
<math>
a_0(\pi G\rho_0)^{1 / 2}
\frac{\partial \mathbf{u}}{\partial t}  
+
\biggl[ \mathbf{V}_0 \cdot \nabla\biggr]a_0(\pi G\rho_0)^{1 / 2}\mathbf{u} 
+
\biggl[ a_0(\pi G\rho_0)^{1 / 2}\mathbf{u}\cdot \nabla\biggr]\mathbf{V}_0 
+
(\pi G \rho_0 a_0^2)
\cancelto{\mathrm{small}}{\biggl[ \mathbf{u}\cdot \nabla\biggr]\mathbf{u} }
+ 2\boldsymbol\omega \boldsymbol\times \mathbf{u} (a_0^2 \pi G\rho_0)^{1 / 2}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> 
-~ (\pi G \rho_0^2 a_0^2)\rho_0^{-1} \nabla p + (\pi G \rho_0 a_0^2)\mathbf\nabla \psi
\, .
</math>
  </td>
</tr>
</table>
Finally, dividing thru by <math>(\pi G \rho_0 a_0)</math>, 
<table border="0" cellpadding="3" align="center">

<tr>
  <td align="right">
<math>
(\pi G\rho_0)^{-1 / 2}
\frac{\partial \mathbf{u}}{\partial t}  
+
\biggl[ \mathbf{V}_0 \cdot \nabla\biggr](\pi G\rho_0)^{-1 / 2}\mathbf{u} 
+
\biggl[ (\pi G\rho_0)^{-1 / 2}\mathbf{u}\cdot \nabla\biggr]\mathbf{V}_0 
+ 2\boldsymbol\omega \boldsymbol\times \mathbf{u} (\pi G\rho_0)^{-1 / 2}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> 
-~ a_0 \nabla p + a_0\mathbf\nabla \psi
\, ,
</math>
  </td>
</tr>
</table>
and adopting the characteristic timescale, <math>t_0 \equiv (\pi G \rho_0)^{-1 / 2}</math>, we have the dimensionless Euler equation,
<table border="0" cellpadding="3" align="center">

<tr>
  <td align="right">
<math>
t_0 \underbrace{\biggl[
\frac{\partial }{\partial t}  
+
\mathbf{V}_0 \cdot \nabla \biggr]}_{D/Dt}\mathbf{u} 
+
t_0 \overbrace{[ \mathbf{u}\cdot \nabla ]\mathbf{V}_0}^{L\mathbf{u}} 
+ 2t_0 ~\overbrace{\boldsymbol\omega \boldsymbol\times \mathbf{u}}^{\Omega\mathbf{u} }</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math> 
-~ a_0 \nabla p + a_0\mathbf\nabla \psi
\, .
</math>
  </td>
</tr>
<tr>
  <td align="center" colspan="3">
{{ LL96b }}, &sect;2, p. 929, Eqs. (2.4) &amp; (2.7)
  </td>
</tr>
</table>

=See Also=
<ul>
<li>
MediaWiki chapter on [[ThreeDimensionalConfigurations/RiemannStype#Riemann_S-Type_Ellipsoids|Riemann S-Type Ellipsoids]]  &nbsp; &nbsp; <== &nbsp; &nbsp; includes [[ThreeDimensionalConfigurations/RiemannStype#SAdata|Self-Adjoint Data Table]] from EFE and from Howard's high-precision (20-digit accuracy) evaluation.
</li>
<li>
MediaWiki chapter on [[ThreeDimensionalConfigurations/Stability/RiemannEllipsoids#Lebovitz_&_Lifschitz_(1996)|Lebovitz &amp; Lifschitz (1996)]]&nbsp; &nbsp; <== &nbsp; &nbsp; Analytic expressions that give the location of the Self-Adjoint Sequences in the EFE Diagram is presented here. 
</li>
<li>
{{ OTM2007full }}, titled, ''Further Evidence for an Elliptical Instability in Rotating Fluid Bars and Ellipsoidal Stars.''
</li>
</ul>

{{ SGFfooter }}