Editing Appendix/Ramblings/ForCohlHoward (section)

=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>.