Appendix/Ramblings/ForCohlHoward

From jetwiki
Jump to navigation Jump to search

Discussions With Howard Cohl

These discussions began in late 2021, when Howard 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 Lebovitz & Lifschitz (1996).

Understanding the Dimensionality of EFE Index Symbols

Howard put together a Mathematica script intended to provide — for any specification of the semi-axis length triplet (a,b,c) — very high-precision, numerical evaluations of any of the index symbols, Aijk and Bijk as defined by Eqs. (103 - 104) in §21 of [EFE]. Originally I suggested that, without loss of generality, he should only need to specify the pair of length ratios, (1,b/a,c/a). In response, Howard pointed out that evaluation of all but a few of the lowest-numbered index symbols — as defined by [EFE] — does explicitly depend on specification of (various powers of) the semi-axis length, a.

Joel's response:  Howard is correct! He should leave the explicit dependence of a — to various powers — 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, a; these evaluations should involve combining various index symbols in such a way that the dependence on a disappears. Consider, for example, our accompanying discussion (click to see relevant expressions) of the virial-equilibrium-based determination of the frequency ratio, fζ/Ωf, in equilibrium S-Type Riemann Ellipsoids. Although most of the required index symbols, A1,A2,A3 and B12, are dimensionless parameters, the index symbol A12 has the unit of inverse-length-squared. Notice, however, that when A12 appears along with any of these other dimensionless parameters in the definition of f, it is accompanied by an extra "length-squared" factor, such as a2. Hence, although I strongly agree that Howard should continue to include various powers of a (etc.) in his Mathematica notebook expressions, I suspect that, without loss of generality, in the end we will always be able to set a=1 and only need to specify the pair of length ratios, (1,b/a,c/a).

Evaluation of Index Symbols

Three Lowest-Order Expressions

In our accompanying derivation of expressions for the three lowest-order index symbols Ai, we have used subscripts (,m,s) instead of (1,2,3) in order to identify which associated semi-axis length is (largest, medium-length, smallest). We have derived the following expressions:

Aaamas

=

2a3p2sin3α[F(α,p)E(α,p)];

Amaamas

=

2a3[E(α,p)(1p2)F(α,p)(as/am)p2sinαp2(1p2)sin3α];

Asaamas

=

2a3[(am/as)sinαE(α,p)(1p2)sin3α].

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

A1

=

2a2a3a12msin3(ϕ)[EllipticF[ϕ,m]EllipticE[ϕ,m]];

A2

=

2a2a3a12m(1m)sin3(ϕ)[EllipticE[ϕ,m]cos2θEllipticF[ϕ,m]a3a2sin2θsinϕ];

A3

=

2a2a3a12(1m)sin3(ϕ)[a2a3sin(ϕ)EllipticE[ϕ,m]].

With a little study it should be clear that our derived expressions for Ai precisely match Howard's Mathematica-notebook expressions when =1, m=2, and s=3, that is, in all cases for which a1>a2>a3. But there will be models to consider (for example, in the uppermost region of the so-called "horn-shaped" region for S-Type Riemann Ellipsoids) for which a1>a3>a2, in which case care must be taken in assigning the proper expressions to A2 and A3. Similarly note that most of the Riemann models of Type I, II, and III — see, for example, Figure 16 (p. 161) in Chapter 7 of [EFE] — have either a2>a1 or a3>a1.

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 A1,A2,A3? Specifically …

Index-Symbol Recurrence Relations

Bijk

=

Ajklai2Aijkl

[EFE], §21, p. 54, Eq. (105)

ai2Aiklaj2Ajkl

=

+(ai2aj2)Bijk

[EFE], §21, p. 54, Eq. (106)

AiklAjkl

=

(ai2aj2)Aijk

[EFE], §21, p. 54, Eq. (107)

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

A1A2

=

(a12a22)A12

A12

=

A2A1(a12a22)

a12A12

=

A2A1(1a22/a12);

and, from the first of the relations,

B12

=

A2a12A12

 

=

A2a12[A2A1(a12a22)]

 

=

1(a12a22)[(a12a22)A2a12(A2A1)]

 

=

1(a12a22)[a12A1a22A2]

 

=

1(1a22/a12)[A1a22a12A2].


Also, consider using the set of relations labeled "LEMMA 7" on p. 54 of [EFE].

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:

Figures circa Year 2000

Approximately four years after 📚 Lebovitz & Lifschitz (1996) 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 📚 Lebovitz & Lifschitz (1996).

Image i3.png

Howard's plot of this data — his image i3.png — is shown immediately below; the abscissa is 0b/a1 and the ordinate is 0c/a1.

Howard's "i3.png" image

i3.png

Compare with Figure 2b of 📚 N. R. Lebovitz, & A. Lifschitz (1996, ApJ, Vol. 458, pp. 699 - 713)

Image i5.png

In an effort to better examine growth-rate trends in the lower-left quadrant of this 📚 Lebovitz & Lifschitz (1996) figure, Howard plotted the same set of stability-analysis data on an axis pair where the abscissa is still 0b/a1, but where, for each value of b/a, the ordinate extends from the lower self-adjoint sequence to the upper self-adjoint sequence — labeled, respectively, x=+1 and x=1 in the classic EFE diagram (EFE, §49, p. 147, Fig. 15 or, see our accompanying discussion). This is displayed immediately below as Howard's "i5.png" image.

Howard's "i5.png" image

i5.png

Notice …

In generating his "i5.png" image, precisely how did Howard "stretch" the ordinate from c/a (as used in his "i3.png" image) to an ordinate ranging from the lower to the upper self-adjoint sequences? Drawing from 📚 Lebovitz & Lifschitz (1996) I presume that, for a given point in the EFE diagram (b/a,c/a), Howard used the expression,

x±

=

C±C21,

📚 Lebovitz & Lifschitz (1996), §2, p. 701, Eq. (8)

where,

C

abB12c2A3b2a2A12,

📚 Lebovitz & Lifschitz (1996), §2, p. 701, Eq. (6)

Then I presume that the ordinate, y — which runs from zero to unity in the "i5.png" image — is determined from the expression,

y

=

0.5(1x±).

Is this the way Howard generated "i5.png"?


In an email dated 26 January 2022, Howard provided the following answer to this question —

  • I had numerical data that I think Norman provided me for the lower and upper self-adjoint sequences.
  • They were simply two sets of curvilinear data in the (b/a,c/a) diagram. Call ay=b/a and then az=c/a.
  • 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.
  • Now consider that you have some data "g" which gives you points in the (ay,az) plane, then g=g(ay,az).
  • 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.
  • For every data point g, there is a fixed ay coordinate value. Normally you would plot that point at (ay,az).
  • 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))
  • So if az=lsa(ay) then it appears at the ordinate value of zero.
  • and if az=usa(ay) then it appears at the ordinate value of unity.
  • So the whole horn shaped region is mapped into the unit square.

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, 0.3b/a1.0 while the ordinate spans the parameter range, 0.96y1.

Howard's "i4.png" image

i4.png

Notice …

Summary

Howard is interested in understanding — in greater detail than appears in 📚 Lebovitz & Lifschitz (1996) — what gives rise to, and what is the extent of these various bands of instability in the classic EFE diagram. Explicit comments/questions:

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

Self-Adjoint Sequences

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?"

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


I have added a subsection to my online chapter discussion of 📚 Lebovitz & Lifschitz (1996) in which I derive an expression whose solution/root should map out the **upper** boundary (x = -1) of the horned-shape region. Click here to see the entire derivation; this derivation ends with the following recommended strategy:

STRATEGY for finding the locus of points that define the upper boundary of the horned-shape region …    Set a=1, and pick a value for 0<b<1; then, using an iterative technique, vary c until the following expression is satisfied:

[c2(a+b)+ab2][(b/a)sinθ(c/a)E(θ,k)(1k2)sin3θ]

=

a2b+bc(ab){[F(θ,k)E(θ,k)k2sin3θ]}.

Choose another value of 0<b<1, then iterate again to find the value of c that corresponds to this new, chosen value of b. Repeat!


Related remarks:

  1. 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?
  2. 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).
  3. Shortly, I will derive the complementary expression that maps out the "lower" boundary (x = +1).
Howard's email response on 1/27/2022

Dear Joel:

I worked with your expression on your website. I proved that this condition for the upper self- adjoint (USA) sequence to be true is equivalent to there being a root of the following equation

F(b, c) := −bc2 (b + 2) c 2 − b 2 + c(b 2 + c 2 (1 + b(b + 3)) (c 2 − b 2 ) √ 1 − c 2 (b + 1) E

cos−1 (c), 1 − b 2 1 − c 2


− bc √ 1 − c 2 (b + 1) F

cos−1 (c), 1 − b 2 1 − c 2


Note that this equation seems to have two roots as a function of c at a fixed <math>b. One corresponds to the USA sequence and the other occurs at c = 0. This can be seen by examining the Maclaurin series approximation as c → 0, namely F(c, b) = − c b + 1

bK( √ 1 − b 2 ) + E( √ 1 − b 2

+ O(c 2 ), (2)

where K and E are the complete elliptic integrals of the first kind. Another interesting limit is the evaluation of the F function when b = c. However, I am unable to do this at the moment. For this k → 1 and the leading order term which comes from E cancels with the constant term. However, it seems as though one then needs the first order term to get an accurate expression for F(b, b) since the remaining term from the third term does not give the correct answer (I checked this numerically). I know it’s not important because this won’t correspond to the USA sequence, but I just thought it was interesting. In any case to use (1) one has to avoid the root at 0 to find the second root which corresponds to the USA sequence. I used this equation to compute the USA to high precision greater than 20 decimal places.

COLLADA 3D Animations

Tiled Menu

Appendices: | VisTrailsEquations | VisTrailsVariables | References | Ramblings | VisTrailsImages | myphys.lsu | ADS |