SSC/FreeEnergy/PolytropesEmbedded: Difference between revisions

From jetwiki
Jump to navigation Jump to search
 
Line 9: Line 9:
     <ul>
     <ul>
<li>[[SSC/FreeEnergy/PolytropesEmbedded/Pt3A|Focus on Five-One Free-Energy Expression]]</li>
<li>[[SSC/FreeEnergy/PolytropesEmbedded/Pt3A|Focus on Five-One Free-Energy Expression]]</li>
<li>[[SSC/FreeEnergy/PolytropesEmbedded/Pt3B|Focus on One-One Free-Energy Expression]]</li>
<li>[[SSC/FreeEnergy/PolytropesEmbedded/Pt3B|Focus on Zero-Zero Free-Energy Expression]]</li>
<li>[[SSC/FreeEnergy/PolytropesEmbedded/Pt3C|Overview]]</li>
<li>[[SSC/FreeEnergy/PolytropesEmbedded/Pt3C|Overview]]</li>
     </ul>
     </ul>

Latest revision as of 13:31, 15 October 2023

Chapter Divisions

In October 2023, this very long chapter was subdivided in order to more effectively accommodate edits. Here is a list of the resulting set of shorter chapters:

  1. Free-Energy Synopsis
  2. Free-Energy of Truncated Polytropes
  3. Free-Energy of BiPolytropes

Free-Energy Synopsis

All of the self-gravitating configurations considered below have an associated Gibbs-like free-energy that can be expressed analytically as a power-law function of the dimensionless configuration radius, x. Specifically,

𝔊type*

=

ax1+bx3/n+cx3/j+𝔊0.

Equilibrium Radii and Critical Radii

The first and second (partial) derivatives with respect to x are, respectively,

𝔊type*x

=

ax2(3bn)x3/n1(3cj)x3/j1

 

=

1x2[a(3bn)x(n3)/n(3cj)x(j3)/j],

2𝔊type*x2

=

2ax3+(3bn)(n+3n)x3/n2+(3cj)(j+3j)x3/j2

 

=

1x3{(3bn)(n+3n)x(n3)/n+(3cj)(j+3j)x(j3)/j2a}.

Equilibrium configurations are identified by setting the first derivative to zero. This gives,

0

=

a(3bn)xeq(n3)/n(3cj)xeq(j3)/j

xeq(n3)/n

=

(n3b)[a(3cj)xeq(j3)/j].

bncxeq(n3)/na3c+1jxeq(j3)/j

=

0.

We conclude, as well, that at this equilibrium radius, the second (partial) derivative assumes the value,

[2𝔊type*x2]eq

=

1xeq3{(3bn)(n+3n)x(n3)/n+(3cj)(j+3j)x(j3)/j2a}eq

 

=

1xeq3{(n+3n)[a(3cj)xeq(j3)/j]+(3cj)(j+3j)xeq(j3)/j2a}

 

=

1xeq3{(3cj)[(j+3j)(n+3n)]xeq(j3)/j+(3nn)a}.

Hence, equilibrium configurations for which the second (as well as first) derivative of the free energy is zero are found at "critical" radii given by the expression,

0

=

(3cj)[(j+3j)(n+3n)][xeq(j3)/j]crit+(3nn)a

[xeq(j3)/j]crit

=

[j2a(n3)3c][n(j+3)j(n+3)]1

 

=

a32c[j2(n3)nj].

Examples

Pressure-Truncated Polytropes

For pressure-truncated polytropes of index n, we set, j=1, in which case,

bncxeq(n3)/na3cxeq4

=

0;

[(34π)MtotMSWS](n+1)/nxeq(n3)/n320π(n+1n)(MtotMSWS)2xeq4

=

0;

xeq(n3)/n

=

(n3b)[a+3cxeq4];

 

and

 

[xeq]crit

=

[a(n3)32c(n+1)]1/4.

Case M

More specifically, the expression that describes the "Case M" free-energy surface is,

𝔊K,M*𝔊K,MEnorm

=

3𝒜(RRnorm)1+n(RRnorm)3/n+(4π3)PePnorm(RRnorm)3.


Hence, we have,

a

3𝒜=35𝔣~W𝔣~M2,

b

n=n(4π3)1/n𝔣~A𝔣~M(n+1)/n,

c

4π3(PePnorm),

where the structural form factors for pressure-truncated polytropes are precisely defined here. Therefore, the statement of virial equilibrium is,

0

=

bncxeq(n3)/na3cxeq4

(34π)cxeq4

=

(34π)[bnxeq(n3)/na3]

(PePnorm)xeq4

=

(34π)[(34π)1/n𝔣~A𝔣~M(n+1)/nxeq(n3)/n15𝔣~W𝔣~M2]

 

=

(34π)(n+1)/n𝔣~A𝔣~M(n+1)/nxeq(n3)/n320π𝔣~W𝔣~M2.

And we conclude that,

3c[xeq]crit4

=

(n3)5(n+1)𝔣~W𝔣~M2

(PePnorm)[xeq]crit4

=

120π(n3n+1)𝔣~W𝔣~M2

(34π)(n+1)/n𝔣~A𝔣~M(n+1)/n[xeq]crit(n3)/n320π𝔣~W𝔣~M2

=

120π(n3n+1)𝔣~W𝔣~M2

(34π)(n+1)/n𝔣~A𝔣~M(n+1)/n[xeq]crit(n3)/n

=

120π(n3n+1)𝔣~W𝔣~M2+320π𝔣~W𝔣~M2

 

=

120π(4nn+1)𝔣~W𝔣~M2

[xeq]crit(n3)/n

=

120π(4π3)(n+1)/n(4nn+1)𝔣~W𝔣~A𝔣~M(n1)/n

[xeq]crit

=

[4n15(n+1)(4π3)1/n𝔣~W𝔣~A𝔣~M(n1)/n]n/(n3).

ASIDE:  Let's see what this requires for the case of n=5, where everything is specifiable analytically. We have gathered together:

  • Form factors from here.
  • Hoerdt's equilibrium expressions from here.
  • Conversion from Horedt's units to ours as specified here.

𝔣~M

=

(1+2)3/2

𝔣~W

=

5245[(48321)(1+2)3+tan1()]

𝔣~A

=

3233[tan1()+(21)(1+2)2]

ReqRnorm=ReqRHoredt[4π(n+1)n]1/(n3)

=

{3[(ξe2/3)5(1+ξe2/3)6]}1/2[4π(n+1)n]1/(n3)

 

=

[(1+2)35][π2336]1/2

PePnorm=PePHoredt[(n+1)34π](n+1)/(n3)

=

33[(ξe2/3)3(1+ξe2/3)4]3[(n+1)34π](n+1)/(n3)

 

=

[18(1+2)12][234π]3

So, the radius of the critical equilibrium state should be,

[ReqRnorm]crit4

=

(n3)35(n+1)(322π)(PePnorm)1𝔣~W𝔣~M2

 

=

12235π{(1+2)1218[π234]3}(1+2)3{5245[(48321)(1+2)3+tan1()]}

 

=

π229313{(1+2)1223}{[(48321)+(1+2)3tan1()]};

whereas, each equilibrium configuration has,

[ReqRnorm]4

=

π226312[(1+2)1220].

So the equilibrium state that marks the critical configuration must have a value of that satisfies the relation,

π226312[(1+2)1220]

=

π229313{(1+2)1223}{[(48321)+(1+2)3tan1()]}

2333

=

(48321)+(1+2)3tan1()

[(1+2)3]tan1()

=

1+80324.

The solution is: crit2.223175.


In addition, we know from our dissection of Hoerdt's work on detailed force-balance models that, in the equilibrium state,

(PePnorm)(ReqRnorm)4

=

[θ~n+1(4π)(n+1)(θ~)2]

3cxeq4

=

[θ~n+1(n+1)(θ~)2].

This means that, for any chosen polytropic index, the critical equilibrium state is the equilibrium configuration for which (needs to be checked),

2(92n)θ~n+1

=

3(n3)[(θ~')2θ~(θ~')ξ~].

We note, as well, that by combining the Horedt expression for xeq with our virial equilibrium expression, we find (needs to be checked),

xeqn3

=

4π3[3(n+1)ξ~2+𝔣~W𝔣~M5𝔣~A]n𝔣~M1n.

Case P

First Pass

Alternatively, let's examine the "Case P" free-energy surface. Drawing on Stahler's presentation, we adopt the following radius and mass normalizations:

MSWS=(n+1n)3/2G3/2Kn2n/(n+1)Pe(3n)/[2(n+1)],

RSWS=(n+1n)1/2G1/2Knn/(n+1)Pe(1n)/[2(n+1)].

In terms of these new normalizations, we have,

Rnorm[(GK)nMtot(n1)]1/(n3)

=

(GK)n/(n3)Mtot(n1)/(n3)RSWS(n+1n)1/2G1/2Knn/(n+1)Pe(1n)/[2(n+1)]

 

 

+MSWS(n1)/(n3)[(n+1n)3/2G3/2Kn2n/(n+1)Pe(3n)/[2(n+1)]](n1)/(n3)

 

=

RSWS(MtotMSWS)(n1)/(n3)(n+1n)[3(n1)(n3)]/[2(n3)]G[2n+(n3)3(n1)]/[2(n3)]

 

 

+Knn[2(n1)(n+1)(n3)]/[(n+1)(n3)]Pe(n1)(3n)/[2(n+1)(n3)]Pe(n1)(3n)/[2(n+1)(n3)]

 

=

RSWS(MtotMSWS)(n1)/(n3)(n+1n)n/(n3).

and,

Pnorm[K4nG3(n+1)Mtot2(n+1)]1/(n3)

=

[K4nG3(n+1)]1/(n3)(MtotMSWS)2(n+1)/(n3){(n+1n)3/2G3/2Kn2n/(n+1)Pe(3n)/[2(n+1)]}2(n+1)/(n3)

 

=

(MtotMSWS)2(n+1)/(n3)(n+1n)3(n+1)/(n3)K4n/(n3)G3(n+1)/(n3)

 

 

×G3(n+1)/(n3)Kn4n/(n3){Pe(n3)/[2(n+1)]}2(n+1)/(n3)

 

=

Pe(MtotMSWS)2(n+1)/(n3)(n+1n)3(n+1)/(n3).

Rewriting the expression for the free energy gives,

𝔊K,M*𝔊K,MEnorm

=

3𝒜(RRSWS)1(RnormRSWS)+n(RRSWS)3/n(RnormRSWS)3/n+(4π3)PePnorm(RRSWS)3(RnormRSWS)3

 

=

3𝒜(RRSWS)1[(MtotMSWS)(n1)/(n3)(n+1n)n/(n3)]

 

 

+n(RRSWS)3/n[(MtotMSWS)(n1)/(n3)(n+1n)n/(n3)]3/n

 

 

+(4π3)(MtotMSWS)2(n+1)/(n3)(n+1n)3(n+1)/(n3)(RRSWS)3[(MtotMSWS)(n1)/(n3)(n+1n)n/(n3)]3

 

=

3𝒜(n+1n)n/(n3)(MtotMSWS)(n1)/(n3)(RRSWS)1+n(n+1n)3/(n3)(MtotMSWS)3(n1)/[n(n3)](RRSWS)3/n

 

 

+(4π3)(n+1n)3/(n3)(MtotMSWS)(5n)/(n3)(RRSWS)3.


Therefore, in this case, we have,

a

=

35𝔣~W𝔣~M2(n+1n)n/(n3)(MtotMSWS)(n1)/(n3),

b

=

n(4π3)1/n𝔣~A𝔣~M(n+1)/n(n+1n)3/(n3)(MtotMSWS)3(n1)/[n(n3)],

c

=

4π3(n+1n)3/(n3)(MtotMSWS)(5n)/(n3),

where the structural form factors for pressure-truncated polytropes are precisely defined here. The statement of virial equilibrium is, therefore,

xeq4+α

=

βxeq(n3)/n,

where,

αa3c

=

15𝔣~W𝔣~M2(n+1n)n/(n3)(MtotMSWS)(n1)/(n3){34π(n+1n)3/(n3)(MtotMSWS)(n5)/(n3)}

 

=

320π𝔣~W𝔣~M2(n+1n)(MtotMSWS)2

 

=

(4π35)𝔣~W(n+1n)𝔪2,

βbnc

=

(4π3)1/n𝔣~A𝔣~M(n+1)/n(n+1n)3/(n3)(MtotMSWS)3(n1)/[n(n3)]{34π(n+1n)3/(n3)(MtotMSWS)(n5)/(n3)}

 

=

𝔣~A𝔪(n+1)/n,

𝔪

(34π)1𝔣~M(MtotMSWS).

From a previous derivation, we have,

0

=

bncxeq(n3)/na3cxeq4

 

=

34π(n+1n)3/(n3)(MtotMSWS)(n5)/(n3){(34π)1/n𝔣~A𝔣~M(n+1)/n(n+1n)3/(n3)(MtotMSWS)3(n1)/[n(n3)]}xeq(n3)/n

 

 

34π(n+1n)3/(n3)(MtotMSWS)(n5)/(n3){15𝔣~W𝔣~M2(n+1n)n/(n3)(MtotMSWS)(n1)/(n3)}xeq4

0

=

(34π)(n+1)/n𝔣~A𝔣~M(n+1)/n(MtotMSWS)(n+1)/nxeq(n3)/n320π𝔣~W𝔣~M2(n+1n)(MtotMSWS)2xeq4

 

=

𝔣~A[(34π)1𝔣~M(MtotMSWS)](n+1)/nxeq(n3)/n15(4π3)𝔣~W(n+1n)[(34π)1𝔣~M(MtotMSWS)]2xeq4

which, thankfully, matches! We conclude as well that the transition from stable to dynamically unstable configurations occurs at,

[xeq]crit4

=

[(n3)3(n+1)]α.

When combined with the statement of virial equilibrium at this critical point, we find,

{[(n3)3(n+1)]+1}αβ

=

[xeq]crit(n3)/n

 

=

{[(n3)3(n+1)]α}(n3)/(4n)

[4n3(n+1)]4n(αβ)4n

=

[(n3)3(n+1)](n3)α(n3)

[3n(n3)(n+1n)](3n)[34(n+1n)]4n

=

α3(n+1)β4n

 

=

{(4π35)𝔣~W(n+1n)𝔪2}3(n+1){𝔣~A𝔪(n+1)/n}4n

 

=

𝔣~A4n[(4π35)𝔣~W(n+1n)]3(n+1)𝔪2(n+1)

𝔪2(n+1)

=

[(4π35)𝔣~W(n+1n)]3(n+1)[3n(n3)(n+1n)](3n)[3𝔣~A4(n+1n)]4n

 

=

[(354π)1𝔣~W]3(n+1)[3n(n3)](3n)[3𝔣~A4]4n

 

=

[325n4π(n3)1𝔣~W](3n)[(32524π)𝔣~A𝔣~W]4n.

This also means that the critical radius is,

[xeq]crit4

=

[(n3)3(n+1)](4π35)𝔣~W(n+1n)𝔪2

 

=

[325n4π(n3)1𝔣~W]1𝔪2

[xeq]crit4(n+1)

=

[325n4π(n3)1𝔣~W](n+1)[325n4π(n3)1𝔣~W](3n)[(32524π)𝔣~A𝔣~W]4n

 

=

[4π(n3)325n𝔣~W]2(n1)[(32524π)𝔣~A𝔣~W]4n

[xeq]crit2(n+1)

=

[n(n3)(3254π𝔣~W)](1n)[(3254π𝔣~W)𝔣~A4]2n

 

=

(nn3)(1n)(3254π𝔣~W)(n+1)(𝔣~A4)2n.


The following parallel derivation was done independently. [Note that a factor of 2n/(n-1) appears to correct a mistake made during the original derivation.] Beginning with the virial expression,

βxeq(n3)/n

=

α+xeq4

(34π)(n+1)/n𝔣~A𝔣~M(n+1)/n(MtotMSWS)(n+1)/n[xeq]crit(n3)/n

=

320π𝔣~W𝔣~M2(n+1n)(MtotMSWS)2+(n3)20πn(MtotMSWS)2𝔣~W𝔣~M2

 

=

(n1)10πn(MtotMSWS)2𝔣~W𝔣~M2[2n(n1)]

[xeq]crit(n3)/n

=

2(n1)15n(4π3)1/n𝔣~W𝔣~A𝔣~M(n1)/n(MtotMSWS)(n1)/n[2n(n1)]

[xeq]crit(n3)

=

[2(n1)15n]n(4π3)𝔣~Wn𝔣~An𝔣~M(n1)(MtotMSWS)(n1)[2n(n1)]n

 

=

[2(n1)15n]n(4π3)𝔣~Wn𝔣~An𝔣~M(n1){[20πn(n3)](n1)/2(𝔣~M2𝔣~W)(n1)/2[xeq]crit2(n1)}[2n(n1)]n

 

=

[2(n1)15n]n(4π3)𝔣~W(n+1)/2𝔣~An{[20πn(n3)](n1)/2[xeq]crit2(n1)}[2n(n1)]n

[xeq]crit(n+1)

=

(34π)[15n2(n1)]n[(n3)20πn](n1)/2𝔣~An𝔣~W(n+1)/2[(n1)2n]n

[xeq]crit(n+1)

=

(34π)[1522]n[(n3)20πn](n1)/2𝔣~An𝔣~W(n+1)/2


Also from Stahler's work we know that the equilibrium mass and radius are,

MtotMSWS

=

(n34π)1/2[θ~n(n3)/2ξ~2(θ~')],

ReqRSWS

=

(n4π)1/2[ξ~θ~n(n1)/2].

Additional details in support of an associated PowerPoint presentation can be found here.

Reconcile

[ReqRSWS]crit4

=

[(n3)20π(n+1)](n+1n)(MtotMSWS)2𝔣~W𝔣~M2

(ReqRnorm)crit4

=

120π(n3n+1)(PePnorm)1𝔣~W𝔣~M2

Taking the ratio, the RHS is,

(MtotMSWS)2(PePnorm)

=

PeMtot2[G3(n+1)Mtot2(n+1)K4n]1/(n3)[(n+1n)3/2G3/2Kn2n/(n+1)Pe(3n)/[2(n+1)]]2(n+1n)

 

=

(n+1n)2PeMtot2[G3Mtot2](n+1)/(n3)Kn4n/(n3)[G3Kn4n/(n+1)Pe(n3)/(n+1)]

 

=

(n+1n)2[G3Mtot2][(n3)+(n+1)]/(n3)[Kn[(n+1)+(n3)]/[(n+1)(n3)]]4nPe[(n+1)+(n3)]/(n+1)

 

=

(n+1n)2Mtot4(n1)/(n3)G[6(n1)]/(n3)Kn8(n1)/[(n+1)(n3)]Pe2(n1)/(n+1);

while the LHS is,

(RnormRSWS)4

=

[(GK)nMtot(n1)]4/(n3){(n+1n)1/2G1/2Knn/(n+1)Pe(1n)/[2(n+1)]}4

 

=

(n+1n)2Mtot4(n1)/(n3)G[6(n1)]/(n3)K8n(n1)/[(n3)(n+1)]Pe2(n1)/(n+1).

Q.E.D.

Now, given that,

MSWS4(n1)/(n3)

=

[(n+1n)3/2G3/2Kn2n/(n+1)Pe(3n)/[2(n+1)]]4(n1)/(n3)

 

=

(n+1n)6(n1)/(n3)G6(n1)/(n3)Kn8n(n1)/[(n+1)(n3)]Pe2(n1)/(n+1)

we have,

(RnormRSWS)4

=

(n+1n)2(MtotMSWS)4(n1)/(n3)(n+1n)6(n1)/(n3)

 

=

(MtotMSWS)4(n1)/(n3)(n+1n)4n/(n3)

(RnormRSWS)n3

=

(MtotMSWS)n1(n+1n)n


By inspection, in the specific case of n=5 (see above), this critical configuration appears to coincide with one of the "turning points" identified by Kimura. Specifically, it appears to coincide with the "extremal in r1" along an M1 sequence, which satisfies the condition,

[n3n1]n=5

=

[ξ~θ~n(θ~')]n=5

12

=

31/2[(1+2)1/2]5[31/2(1+2)3/2]1

 

=

3(1+2)1

=

51/2.

Hence, according to Kimura, the turning point associated with the configuration with the largest equilibrium radius, corresponds to the equilibrium configuration having,

|Rmax=52.2360680.

This is, indeed, very close to — but decidedly different from — the value of crit determined, above!


Streamlined

Let's copy the expression for the "Case P" free energy derived above, then factor out a common term:


𝔊K,MEnorm

=

3𝒜(n+1n)n/(n3)(MtotMSWS)(n1)/(n3)(RRSWS)1+n(n+1n)3/(n3)(MtotMSWS)3(n1)/[n(n3)](RRSWS)3/n

 

 

+(4π3)(n+1n)3/(n3)(MtotMSWS)(5n)/(n3)(RRSWS)3.

 

=

(MtotMSWS)(5n)/(n3)(n+1n)3/(n3){3𝒜(n+1n)(MtotMSWS)2(RRSWS)1+n(MtotMSWS)(n+1)/n(RRSWS)3/n+4π3(RRSWS)3}

Defining a new normalization energy,

ESWS

Enorm(MtotMSWS)(5n)/(n3)(n+1n)3/(n3)

 

=

(n+1n)3/2K3n/(n+1)G3/2Pe(5n)/[2(n+1)],

we can write,

𝔊K,M*𝔊K,MESWS

=

3𝒜(n+1n)(MtotMSWS)2(RRSWS)1+n(MtotMSWS)(n+1)/n(RRSWS)3/n+4π3(RRSWS)3,

in which case the coefficients of the generic free-energy expression are,

a

=

35𝔣~W𝔣~M2(n+1n)(MtotMSWS)2=35(4π3)2(n+1n)𝔣~W𝔪2

b

=

n(34π)1/n𝔣~A𝔣~M(n+1)/n(MtotMSWS)(n+1)/n=(4πn3)𝔣~A𝔪(n+1)/n

c

=

4π3,

where, as above,

𝔪

(34π)1𝔣~M(MtotMSWS).

Now, if we define the pair of parameters,

α

a3c

β

bnc,

then the statement of virial equilibrium is,

xeq4+α

=

βxeq(n3)/n,

and the boundary between dynamical stability and instability occurs at,

[xeq]crit4

=

[n33(n+1)]α.

Combining these last two expressions means that the boundary between dynamical stability and instability is associated with the parameter condition,

[xeq]crit(n3)/n

=

[n33(n+1)+1]αβ

{[n33(n+1)]α}(n3)/(4n)

=

[4n3(n+1)]αβ

βα3(n+1)/(4n)

=

[4n3(n+1)][n3nn3(n+1)](3n)/(4n)

 

=

4[n3(n+1)]3(n+1)/(4n)[n3n](3n)/(4n)

(β4)4nα3(n+1)

=

[n3(n+1)]3(n+1)[n3n](3n)

(β4)4n

=

[nα3(n+1)]3(n+1)[nn3]n3.

Case M

a

3𝒜=35𝔣~W𝔣~M2,

b

n=n(4π3)1/n𝔣~A𝔣~M(n+1)/n,

c

4π3(PePnorm).

Hence,

α

=

(4π15)𝔣~W(34π𝔣~M)2(PePnorm)1

β

=

𝔣~A(34π𝔣~M)(n+1)/n(PePnorm)1.

So the dynamical stability conditions are:

(PePnorm)(nn3)[xeq]crit4

=

[(nn+1)4π𝔣~W325](34π𝔣~M)2;

and,

(𝔣~A4)4n(34π𝔣~M)4(n+1)(PePnorm)4n

=

[n3(n+1)]3(n+1)[nn3]n3(4π𝔣~W15)3(n+1)(34π𝔣~M)6(n+1)(PePnorm)3(n+1)

(𝔣~A4)4n

=

[(nn+1)4π𝔣~W325]3(n+1)(34π𝔣~M)2(n+1)[nn3(PePnorm)]n3

[nn3(PePnorm)]n3

=

(𝔣~A4)4n[(nn+1)4π𝔣~W325]3(n+1)(34π𝔣~M)2(n+1).

Together, then,

[xeq]crit4(n3)

=

[(nn+1)4π𝔣~W325]n3(34π𝔣~M)2(n3)[(PePnorm)nn3](n3)

 

=

[(nn+1)4π𝔣~W325]n3(34π𝔣~M)2(n3)(𝔣~A4)4n[(nn+1)4π𝔣~W325]3(n+1)(34π𝔣~M)2(n+1)

 

=

[(nn+1)4π𝔣~W325]4n(34π𝔣~M)4(n1)(𝔣~A4)4n

[xeq]crit(n3)

=

[4𝔣~A(nn+1)4π𝔣~W325]n(34π𝔣~M)(n1).

Case P

a

=

35(4π3)2(n+1n)𝔣~W𝔪2

b

=

(4πn3)𝔣~A𝔪(n+1)/n

c

=

4π3,

where, as above,

𝔪

(34π)1𝔣~M(MtotMSWS).

Hence,

α

=

15(4π3)(n+1n)𝔣~W𝔪2

β

=

𝔣~A𝔪(n+1)/n.

So the dynamical stability conditions are:

[xeq]crit4

=

[n3(n+1)][n3n]15(4π3)(n+1n)𝔣~W𝔪2

 

=

[n3n](4π𝔣~W325)𝔪2

and,

(𝔣~A4)4n𝔪4(n+1)

=

[(4π325)𝔣~W𝔪2]3(n+1)[nn3]n3

𝔪2(n+1)

=

[(4π325)𝔣~W]3(n+1)[nn3](n3)(𝔣~A4)4n.

Together, then,

[xeq]crit4(n+1)

=

[n3n](n+1)(4π𝔣~W325)(n+1)(4π𝔣~W325)3(n+1)[n3n](n3)(𝔣~A4)4n

 

=

[n3n]2(n1)(4π𝔣~W325)2(n+1)(𝔣~A4)4n.

Compare

Let's see if the two cases, in fact, provide the same answer.

(RnormRSWS)n3=[xPxM]critn3

=

{[n3n](4π𝔣~W325)𝔪2}(n3)/4{[4𝔣~A(nn+1)4π𝔣~W325]n(34π𝔣~M)(n1)}1

 

=

{[n3n](4π𝔣~W325)𝔪2}(n3)/4{[4𝔣~A(nn+1)4π𝔣~W325]n(34π𝔣~M)(n1)}1

Five-One Bipolytropes

For analytically prescribed, "five-one" bipolytropes, n=5 and j=1, in which case,

xeq2/5

=

(53b)[a3cxeq2];

 

and

 

[xeq]crit

=

[18ca]1/2.


More specifically, the expression that describes the free-energy surface is,

𝔊51*24(qν2)χeq[𝔊51Enorm]

=

1i2[X3/5(5𝔏i)+X3(4𝔎i)X1(3𝔏i+12𝔎i)].

Hence, we have,

a

3χeq(𝔏i+4𝔎i),

b

5𝔏iχeq3/5,

c

4𝔎iχeq3,

and conclude that,

[χeq]crit

=

[18(4𝔎iχeq3)3χeq(𝔏i+4𝔎i)]crit1/2

 

=

[χeq]crit[24𝔎i(𝔏i+4𝔎i)]1/2

[24𝔎i(𝔏i+4𝔎i)]crit

=

1

[𝔏i𝔎i]crit

=

20.

Also, from our detailed force balance derivations, we know that,

χeqReqRnorm

=

(π23)1/2ν2q(1+i2)333i5.

Zero-Zero Bipolytropes

General Form

In this case, we retain full generality making the substitutions, nnc and jne, to obtain,

xeq(nc3)/nc

=

nc3b[a(3cne)xeq(ne3)/ne];

 

and

 

[xeq(ne3)/ne]crit

=

{ne2(nc3)3[nc(ne+3)ne(nc+3)]}ac

 

=

[ne2(nc3)32(ncne)]ac.


And here, the expression that describes the free-energy surface is,

𝔊00*5(qν2)χeq[𝔊00Enorm]

=

52q3[ncA2X3/nc+neB2X3/ne3C2X1].

Hence, we have,

a3χeq(52q3)C2

=

3fχeq,

bncχeq3/nc(52q3)A2

ncχeq3/nc[1+q3(f1𝔉)],

cneχeq3/ne(52q3)B2

neχeq3/ne(52q3)[25q3fA2]

 

=

neχeq3/ne{f[1+q3(f1𝔉)]},

where the definitions of f and 𝔉 are given below. We immediately deduce that the critical equilibrium state is identified by,

[xeq(ne3)/ne]crit

=

{fne(nc3)3(ncne)}[χeq(ne3)/ne]crit{f[1+q3(f1𝔉)]}1

1f[1+q3(f1𝔉)]

=

1[ne(nc3)3(ncne)]

 

=

nc(3ne)3(ncne).


From our associated detailed-force-balance derivation, we know that the associated equilibrium radius is,

χeq

=

{(π3)22ncνnc1q3nc[1+25q3(f1𝔉)]nc}1/(nc3).


Compare with Five-One

It is worthwhile to set nc=5 and ne=1 in this expression and compare the result to the comparable expression shown above for the "Five-One" Bipolytrope. Here we have,

[χeq]51

=

{(π3)23ν4q2[1+25q3(f1𝔉)]5}1/2

 

=

(π23)1/2ν2q13[1+25q3(f1𝔉)]5/2;

whereas, rewriting the above relation gives,

χeq|51

=

(π23)1/2ν2q13[(1+i2)6/53i2]5/2.

And, here, we should conclude that the critical equilibrium configuration is associated with,

1f[1+q3(f1𝔉)]

=

56.

Free-Energy of Truncated Polytropes

In this case, the Gibbs-like free energy is given by the sum of three separate energies,

𝔊

=

Wgrav+𝔖therm+PeV

 

=

3𝒜[GM2R]+n[KM(n+1)/nR3/n]+4π3PeR3,

where the constants,

𝒜15𝔣~W𝔣~M2

      and     

(4π3)1/n𝔣~A𝔣~M(n+1)/n,

and, as derived elsewhere,

Structural Form Factors for Pressure-Truncated Polytropes (n5)

𝔣~M

=

(3θ~'ξ~)

𝔣~W

=

35(5n)ξ~2[θ~n+1+3(θ~')2𝔣~Mθ~]

𝔣~A

=

1(5n){6θ~n+1+(n+1)[3(θ~')2𝔣~Mθ~]}

As we have shown separately, for the singular case of n=5,

𝔣M

=

(1+2)3/2

𝔣W

=

5245[(48321)(1+2)3+tan1()]

𝔣A

=

3233[tan1()+(21)(1+2)2]

where, ξ~/3


In general, then, the warped free-energy surface drapes across a four-dimensional parameter "plane" such that,

𝔊

=

𝔊(R,K,M,Pe).

In order to effectively visualize the structure of this free-energy surface, we will reduce the parameter space from four to two, in two separate ways: First, we will hold constant the parameter pair, (K,M); giving a nod to Kimura's (1981b) nomenclature, we will refer to the resulting function, 𝔊K,M(R,Pe), as a "Case M" free-energy surface because the mass is being held constant. Second, we will hold constant the parameter pair, (K,Pe), and examine the resulting "Case P" free-energy surface, 𝔊K,Pe(R,M).

Virial Equilibrium and Dynamical Stability

The first (partial) derivative of 𝔊 with respect to R is,

𝔊R

=

1R[3𝒜GM2R13KM(n+1)/nR3/n+4πPeR3];

and the second (partial) derivative is,

2𝔊R2

=

1R2[6𝒜GM2R1+(n+3n)3KM(n+1)/nR3/n+8πPeR3].

The virial equilibrium radius is identified by setting the first derivative to zero. This means that,

3KM(n+1)/nReq3/n

=

3𝒜GM2Req1+4πPeReq3.

This expression can be usefully rewritten in the following forms:

Virial Equilibrium Condition
Case 1:

3(n+3)KM(n+1)/nReq3/n

=

3(n+3)𝒜GM2Req1+4π(n+3)PeReq3

Case 2:

6n𝒜GM2Req1

=

8πnPeReq36nKM(n+1)/nReq3/n

Case 3:

8πnPeReq3

=

6nKM(n+1)/nReq3/n6n𝒜GM2Req1

Dynamical stability is determined by the sign of the second derivative expression evaluated at the equilibrium radius; setting the second derivative to zero identifies the transition from stable to unstable configurations. The criterion is,

0

=

[6n𝒜GM2R1+3(n+3)KM(n+1)/nR3/n+8πnPeR3]Req

Case 1 Stability Criterion

Using the "Case 1" virial expression to define the equilibrium radius means that the stability criterion is,

0

=

6n𝒜GM2Req1+3(n+3)𝒜GM2Req1+4π(n+3)PeReq3+8πnPeReq3

 

=

𝒜GM2Req1[3(n+3)6n]+4πPeReq3[(n+3)+2n]

4πPeReq3[3(n+1)]

=

𝒜GM2Req1[3(n3)]

4πPeReq4(n+1)

=

𝒜GM2(n3)

Case 2 Stability Criterion

Using the "Case 2" virial expression to define the equilibrium radius means that the stability criterion is,

0

=

8πnPeReq36nKM(n+1)/nReq3/n+3(n+3)KM(n+1)/nReq3/n+8πnPeReq3

 

=

16πnPeReq3[3(n3)]KM(n+1)/nReq3/n

16πnPeReq3

=

[3(n3)]KM(n+1)/nReq3/n

16πnPeReq3(n+1)/n

=

[3(n3)]KM(n+1)/n

Case 3 Stability Criterion

Using the "Case 3" virial expression to define the equilibrium radius means that the stability criterion is,

0

=

6n𝒜GM2Req1+3(n+3)KM(n+1)/nReq3/n+6nKM(n+1)/nReq3/n6n𝒜GM2Req1

 

=

12n𝒜GM2Req1+[6n+3(n+3)]KM(n+1)/nReq3/n

9(n+1)KM(n+1)/nReq3/n

=

12n𝒜GM2Req1

Reqn3

=

[4n𝒜3(n+1)]n(GK)nMn1

Case M

Now, in our discussion of "Case M" sequence analyses, the configuration's radius is normalized to,

Rnorm

[GnKnMn1]1/(n3).

Our "Case 3" stability criterion directly relates. We conclude that the transition from stability to dynamical instability occurs when,

[ReqRnorm]critn3

=

[4n𝒜3(n+1)]n

[ReqRnorm]crit(n3)/n

=

4n15(n+1)(4π3)1/n𝔣~W𝔣~A𝔣~M(n1)/n

Also in the "Case M" discussions, the external pressure is normalized to,

Pnorm

[G3(n+1)K4nM2(n+1)]1/(n3).

If we raise the "Case 1" stability criterion expression to the (n3) power, then divide it by the "Case 3" stability criterion expression raised to the fourth power, we find,

[Pe]critn3

=

[𝒜GM2(n3)4π(n+1)]n3{[4n𝒜3(n+1)]n(GK)nMn1}4

 

=

[𝒜(n3)4π(n+1)]n3Gn3M2(n3)[3(n+1)4n𝒜]4n(KG)4nM4(1n)

 

=

[(n3)4π(n+1)]n3[3(n+1)4n]4n𝒜3(n+1)4nK4nM2(n+1)G3(n+1)

[PePnorm]critn3

=

[(n3)4π(n+1)]n3[3(n+1)4n]4n𝒜3(n+1)4n

 

=

[(n3)4π(n+1)]n3[3(n+1)4n]4n[5𝔣~M2𝔣~W]3(n+1)(34π)4[𝔣~A𝔣~M(n+1)/n]4n

 

=

(34π)4[(n3)4π(n+1)]n3[3(n+1)4n]4n[5𝔣~W]3(n+1)𝔣~M2(n+1)𝔣~A4n

Case P

Flipping around this expression for [Pe]crit, we also can write,

[M]crit2(n+1)

=

[(n3)4π(n+1)]n3[3(n+1)4n]4n𝒜3(n+1)4nK4nG3(n+1)Pe3n.

Now, in our "Case P" discussions we normalized the mass to

MSWS

(n+1n)3/2G3/2K2n/(n+1)Pe(3n)/[2(n+1)].

Hence, we have,

[MMSWS]crit2(n+1)

=

[(n3)4π(n+1)]n3[3(n+1)4n]4n𝒜3(n+1)4n(n+1n)3(n+1)

 

=

[(n3)4πn]n3(34)4n𝒜3(n+1)4n,

where the constants,

𝒜15𝔣~W𝔣~M2

      and     

(4π3)1/n𝔣~A𝔣~M(n+1)/n.

So we can furthermore conclude that,

[MMSWS]crit2(n+1)

=

[(n3)4πn]n3(34)4n{15𝔣~W𝔣~M2}3(n+1){(4π3)1/n𝔣~A𝔣~M(n+1)/n}4n

 

=

(34π)4[(n3)4πn]n3(3𝔣~A4)4n[53𝔣~M2𝔣~W3](n+1).


Our expression for [M]crit2(n+1) can also be combined with the "Case 2 stability criterion" to eliminate the mass entirely, giving,

{16πnPeReq3(n+1)/n}2n

=

{[3(n3)]K}2n[(n3)4π(n+1)]n3[3(n+1)4n]4n𝒜3(n+1)4nK4nG3(n+1)Pe3n

Req6(n+1)

=

[3(n3)16πn]2n[(n3)4π(n+1)]n3[3(n+1)4n]4n𝒜3(n+1)6nK6nG3(n+1)Pe3(1n)

Req2(n+1)

=

{[(n3)4πn]2n[(n3)4πn]n3[(n+1)n]4n+(3n)(34)6n}1/3𝒜(n+1)2nK2nG(n+1)Pe(1n)

 

=

[(n3)4πn](n1)[(n+1)n](n+1)(34)2n𝒜(n+1)2nK2nG(n+1)Pe(1n).

Finally, recognizing that in our "Case P" discussions we normalized the radius to

RSWS

(n+1n)1/2G1/2Kn/(n+1)Pe(1n)/[2(n+1)],

we have,

[Req]crit2(n+1)

=

[(n3)4πn](n1)(n+1n)(n+1)(34)2n𝒜(n+1)2n{RSWS(n+1n)1/2}2(n+1)

[ReqRSWS]crit2(n+1)

=

[(n3)4πn](n1)(34)2n𝒜(n+1)2n

 

=

[nn3](1n)(4π)1n(34)2n[15𝔣~W𝔣~M2](n+1)[(4π3)1/n𝔣~A𝔣~M(n+1)/n]2n

 

=

[nn3](1n)(4π)1n232n+242n[5𝔣~M2𝔣~W](n+1)[𝔣~A2n𝔣~M2(n+1)]

 

=

[nn3](1n)[3254π𝔣~W](n+1)[𝔣~A4]2n.

Case M Free-Energy Surface

It is useful to rewrite the free-energy function in terms of dimensionless parameters. Here we need to pick normalizations for energy, radius, and pressure that are expressed in terms of the gravitational constant, G, and the two fixed parameters, K and M. We have chosen to use,

Rnorm

[(GK)nMtotn1]1/(n3),

Pnorm

[K4nG3(n+1)Mtot2(n+1)]1/(n3),

which, as is detailed in an accompanying discussion, are similar but not identical to the normalizations used by Horedt (1970) and by Whitworth (1981). The self-consistent energy normalization is,

Enorm

PnormRnorm3.

As we have demonstrated elsewhere, after implementing these normalizations, the expression that describes the "Case M" free-energy surface is,

𝔊K,M*𝔊K,MEnorm=3𝒜(RRnorm)1+n(RRnorm)3/n+(4π3)PePnorm(RRnorm)3,

Given the polytropic index, n, we expect to obtain a different "Case M" free-energy surface for each choice of the dimensionless truncation radius, ξ~; this choice will imply corresponding values for θ~ and θ~' and, hence also, corresponding (constant) values of the coefficients, 𝒜 and .


Case P Free-Energy Surface

Again, it is useful to rewrite the free-energy function in terms of dimensionless parameters. But here we need to pick normalizations for energy, radius, and mass that are expressed in terms of the gravitational constant, G, and the two fixed parameters, K and Pe. As is detailed in an accompanying discussion, we have chosen to use the normalizations defined by Stahler (1983), namely,

RSWS

(n+1nG)1/2Kn/(n+1)Pe(1n)/[2(n+1)],

MSWS

(n+1nG)3/2K2n/(n+1)Pe(3n)/[2(n+1)].

The self-consistent energy normalization is,

ESWS(nn+1)GMSWS2RSWS

=

(n+1n)3/2G3/2K3n/(n+1)Pe(5n)/[2(n+1)].

After implementing these normalizations — see our accompanying analysis for details — the expression that describes the "Case P" free-energy surface is,

𝔊K,Pe*𝔊K,PeESWS

=

3𝒜(n+1n)(MMSWS)2(RRSWS)1+n(MMSWS)(n+1)/n(RRSWS)3/n+4π3(RRSWS)3.

Given the polytropic index, n, we expect to obtain a different "Case P" free-energy surface for each choice of the dimensionless truncation radius, ξ~; this choice will imply corresponding values for θ~ and θ~' and, hence also, corresponding (constant) values of the coefficients, 𝒜 and .

Summary

  DFB Equilibrium Onset of Dynamical Instability
Case M:

[ReqRnorm]n3

=

[4π(n+1)n]ξ~(n3)(ξ~2θ'~)(1n)

[ReqRnorm]critn3

=

[4n15(n+1)]n(4π3)𝔣~Wn𝔣~An𝔣~M(n1)

[PePnorm]n3

=

[(n+1)34π](n+1)θ~(n+1)(n3)(ξ~2θ'~)2(n+1)

[PePnorm]critn3

=

(34π)4[(n3)4πn]n3(n+1n)3(n+1)[53𝔣~M2𝔣~W3]n+1(3𝔣~A4)4n

Case P:

[ReqRSWS]2

=

(n4π)ξ~2θ~n1

[ReqRSWS]crit2

=

[3254π𝔣~W][n3n](n1)/(n+1)[𝔣~A4]2n/(n+1)

[MtotMSWS]2

=

(n34π)θ~n3(ξ~2θ'~)2

[MtotMSWS]crit2

=

[53𝔣~M2𝔣~W3](34π)4/(n+1)[(n3)4πn](n3)/(n+1)(3𝔣~A4)4n/(n+1)

In all four cases, the expression on right intersects (is equal to) the expression on the left when the following condition applies:

For n5:      

2(92n)θ~n+1

=

3(n3)[(θ~')2(θ~θ~'ξ~)];

For n=5:      

[2453]3

=

(1+2)3tan1+(41).


If (for n5) we adopt the shorthand notation,

Υ

[3(θ~')2𝔣~Mθ~]=3[(θ~')2(θ~θ~'ξ~)],

and

τ

θ~n+1,

then the critical condition becomes,

(n3)Υ

=

2(92n)τ,

and at the critical state, the expressions for the structural form-factors become,

𝔣~A

=

1(5n)[6τ+(n+1)Υ]

 

=

1(5n){6+(n+1)[2(92n)n3]}τ

 

=

1(5n)[6(n3)+2(92n)(n+1)n3]τ

 

=

1(5n)[4n(5n)n3]τ

 

=

4nτ(n3);

𝔣~W

=

35(5n)ξ~2[τ+Υ]

 

=

35(5n)ξ~2{1+[2(92n)n3]}τ

 

=

35(5n)ξ~2[3(5n)n3]τ

 

=

325τ(n3)ξ~2

53𝔣~M2𝔣~W3

=

[(n3)ξ~232τ]3(3θ~'ξ~)2=32[(n3)ξ~232τ]3(ξ~2θ~'ξ~3)2

 

=

[(n3)334τ3](ξ~2θ~')2.

Hence (1),

[ReqRnorm]critn3

=

[4π(n+1)n][4n15]n(13)[3254nξ~2]n𝔣~M1n

 

=

[4π(n+1)n][1ξ~2n](ξ~2θ~'ξ~3)1n

 

=

[4π(n+1)n]ξ~n3(ξ~2θ~')1n

Q.E.D.


And (2),

[PePnorm]critn3

=

(34π)4[(n3)4πn]n3(n+1n)3(n+1)[53𝔣~M2𝔣~W3]n+1(3𝔣~A4)4n

 

=

34(4π)(n+1)(n3n)n3(n+1n)3(n+1)[(n3)334τ3]n+1(ξ~2θ~')2(n+1)[3nτn3]4n

 

=

34(4π)(n+1)(n3n)n3(n+1n)3(n+1)n3(n+1)[n3n]3(n+1)(ξ~2θ~')2(n+1)[nn3]4nτ4n3(n+1)34n4(n+1)

 

=

[(n+1)34π]n+1(ξ~2θ~')2(n+1)τn3

[PePnorm]critn3

=

[(n+1)34π]n+1(ξ~2θ~')2(n+1)θ~(n+1)(n3).

Q.E.D.


And (3),

[ReqRSWS]crit2(n+1)

=

[3254π𝔣~W]n+1[n3n](n1)[𝔣~A4]2n

 

=

[3254π]n+1[n3n](n1)[nτn3]2n[(n3)ξ~2325τ]n+1

 

=

[14π]n+1[nn3]n+1[(n3)ξ~2]n+1τn1

[ReqRSWS]crit2

=

(n4π)ξ~2θ~n1

Q.E.D.

And (4),

[MtotMSWS]crit2(n+1)

=

[53𝔣~M2𝔣~W3]n+1(34π)4[(n3)4πn](n3)(3𝔣~A4)4n

 

=

{[(n3)334τ3](ξ~2θ~')2}n+134(4π)(n+1)(n3n)(n3)[3nτn3]4n

 

=

[n34π]n+1(ξ~2θ~')2(n+1)τn3

[MtotMSWS]crit2

=

[n34π](ξ~2θ~')2θ~n3

Q.E.D.

Free-Energy of Bipolytropes

In this case, the Gibbs-like free energy is given by the sum of four separate energies,

𝔊

=

[Wgrav+𝔖therm]core+[Wgrav+𝔖therm]env.

In addition to specifying (generally) separate polytropic indexes for the core, nc, and envelope, ne, and an envelope-to-core mean molecular weight ratio, μe/μc, we will assume that the system is fully defined via specification of the following five physical parameters:

  • Total mass, Mtot;
  • Total radius, R;
  • Interface radius, Ri, and associated dimensionless interface marker, qRi/R;
  • Core mass, Mc, and associated dimensionless mass fraction, νMc/Mtot;
  • Polytropic constant in the core, Kc.

In general, the warped free-energy surface drapes across a five-dimensional parameter "plane" such that,

𝔊

=

𝔊(R,Kc,Mtot,q,ν).

Order of Magnitude Derivation

Let's begin by providing very rough, approximate expressions for each of these four terms, assuming that nc=5 and ne=1.

Wgrav|core

𝔞c[GMtotMc(Ri/2)]=2𝔞c[GMtot2R(νq)];

Wgrav|env

𝔞e[GMtotMe(Ri+R)/2]=2𝔞e[GMtot2R(1ν1+q)];

𝔖therm|core=Uint|core

𝔟cncKcMc(ρ¯c)1/nc=5𝔟cKcMtotν[3Mc4πRi3]1/5

 

=

𝔟c(35522π)1/5Kc(Mtotν)6/5(Rq)3/5;

𝔖therm|env=Uint|env

𝔟eneKeMenv(ρ¯e)1/ne=𝔟eKeMtot(1ν)[3Menv4π(R3Ri3)]

 

=

𝔟e(322π)Ke[Mtot(1ν)]2[R3(1q3)]1.

In writing this last expression, it has been necessary to (temporarily) introduce a sixth physical parameter, namely, the polytropic constant that characterizes the envelope material, Ke. But this constant can be expressed in terms of Kc via a relation that ensures continuity of pressure across the interface while taking into account the drop in mean molecular weight across the interface, that is,

Ke(ρ¯e)(ne+1)/ne

Kc(ρ¯c)(nc+1)/nc

Ke[(μeμc)ρ¯c]2

Kc(ρ¯c)6/5

KeKc(μeμc)2

[3Mtotν4π(Rq)3]4/5.

Hence, the fourth energy term may be rewritten in the form,

𝔖therm|env=Uint|env

𝔟e(322π)(μeμc)2Kc[3Mtotν4π(Rq)3]4/5[Mtot(1ν)]2[R3(1q3)]1

 

=

𝔟e(322π)1/5(μeμc)2KcMtot6/5R3/5[q3ν]4/5(1ν)2(1q3).

Putting all the terms together gives,

𝔊

2𝔞c[GMtot2R(νq)]2𝔞e[GMtot2R(1ν1+q)]+𝔟c(35522π)1/5Kc(Mtotν)6/5(Rq)3/5

 

 

+𝔟e(322π)1/5(μeμc)2KcMtot6/5R3/5[q3ν]4/5(1ν)2(1q3)

 

=

2𝒜biP[GMtot2R]+biPKc[(νMtot)2qR]3/5

𝔊Enorm

=

2𝒜biP[GMtot2R](G3Kc5)1/2+biP(ν2q)3/5Kc[Mtot2R]3/5(G3Kc5)1/2

 

=

2𝒜biP[RnormR]+biP(ν2q)3/5[RnormR]3/5,

where,

𝒜biP

[𝔞c(νq)+𝔞e(1ν1+q)],

biP

(322π)1/5[5𝔟c+𝔟e(μeμc)2q3(1ν)2ν2(1q3)].

Equilibrium Radius

Order of Magnitude Estimate

This means that,

𝔊R

=

+2𝒜biP[GMtot2R2]35biPKc[ν2q]3/5Mtot6/5R8/5.

Hence, because equilibrium radii are identified by setting 𝔊/R=0, we have,

ReqRnorm

=

(253)5/2[𝒜biPbiP]5/2(qν2)3/2.

Reconcile With Known Analytic Expression

From our earlier derivations, it appears as though,

χeqReqRnorm

=

(3825π)1/2(324)(qi)5(νq3)2(1+i2)3

 

=

(253)5/2(qν2)3/2[(π28355)1/2(ν2q)5/2(1+i2)3i5].

This implies that,

𝒜biPbiP

[(π28355)1/2(ν2q)5/2(1+i2)3i5]2/5

 

=

(ν2q)(π28355)1/5(1+i2)6/5i2

[𝔞c(νq)+𝔞e(1ν1+q)]

1225(ν2q)(1+i2)6/5i2[5𝔟c+𝔟e(μeμc)2q3(1ν)2ν2(1q3)]

[𝔞c+𝔞eq(1ν)ν(1+q)]

ν225(1+i2)6/5i2[5𝔟c+𝔟e(μeμc)2q3(1ν)2ν2(1q3)]

Focus on Five-One Free-Energy Expression

Approximate Expressions

Let's plug this equilibrium radius back into each term of the free-energy expression.

WgravEnorm|core

2𝔞c(G3Kc5)1/2[GMtot2Req(νq)]

 

=

2𝔞c(νq)[RnormReq];

WgravEnorm|env

2𝔞e(G3Kc5)1/2[GMtot2Req(1ν1+q)]

 

=

2𝔞e(1ν1+q)[RnormReq];

ScoreEnorm=[3(γc1)2]UintEnorm|core

[325]𝔟c(35522π)1/5(G3Kc5)1/2Kc(Mtotν)6/5(Reqq)3/5

 

=

[325]𝔟c(35522π)1/5(ν2q)3/5(RnormReq)3/5;

SenvEnorm=[3(γe1)2]UintEnorm|env

[32]𝔟e(322π)1/5(μeμc)2(G3Kc5)1/2KcMtot6/5Req3/5[q3ν]4/5(1ν)2(1q3)

 

=

[32]𝔟e(322π)1/5(μeμc)2[q3ν]4/5(1ν)2(1q3)(RnormReq)3/5.

From Detailed Force-Balance Models

In the following derivations, we will use the expression,

χeqReqRnorm

=

(μeμc)3(π23)1/21A2ηs=(π23)1/2ν2q(1+i2)333i5.

Keep in mind, as well — as derived in an accompanying discussion — that,

νMcoreMtot

=

(m32i3)(1+i2)1/2[1+(1m3)2i2]1/2[m3i+(1+i2)(π2+tan1Λi)]1,

where,

m33(μeμc).

From the accompanying Table 1 parameter values, we also can write,

q

=

ηiηs=ηi{π2+ηi+tan1[1ηii]}1

 

=

ηi{ηi+cot1[i1ηi]}1,

where,

ηi

=

m3[i(1+i2)].

Let's also define the following shorthand notation:

𝔏i

(i41)i2+(1+i2)3i3tan1i;

𝔎i

(1+Λi2)ηi[π2+tan1Λi]+Λiηi.


Gravitational Potential Energy of the Core

Pulling from our detailed derivations,

[WcoreEnorm]eq

=

(3825π)1/2[i(i483i21)(1+i2)3+tan1(i)].

χeq[WcoreEnorm]eq

=

(3825π)1/2[i(i483i21)(1+i2)3+tan1(i)](π23)1/2ν2q(1+i2)333i5

 

=

(324)ν2q1i5[i(i483i21)+(1+i2)3tan1(i)]

Out of equilibrium, then, we should expect,

WcoreEnorm

=

χ1(324)ν2q1i5[i(i483i21)+(1+i2)3tan1(i)]

 

=

χ1(324)ν2q1i2[𝔏i83],

which, in comparison with our above approximate expression, implies,

𝔞c

=

(325)νi5[i(i483i21)+(1+i2)3tan1(i)].

Thermal Energy of the Core

Again, pulling from our detailed derivations,

[ScoreEnorm]eq

=

12(3825π)1/2[i(i41)(1+i2)3+tan1(i)]

χeq3[ScoreEnorm]eq5

=

125(3825π)5/2[i(i41)(1+i2)3+tan1(i)]5[(π23)1/2ν2q(1+i2)333i5]3

 

=

1π(322)11(ν2q)3[i(i41)(1+i2)3+tan1(i)]5[(1+i2)9i15].

Out of equilibrium, we should then expect,

ScoreEnorm

=

(322π)1/5[χ1(ν2q)1(1+i2)2]3/5(322)2𝔏i.

In comparison with our above approximate expression, we therefore have,

[(325)𝔟c(35522π)1/5(ν2q)3/5]5

=

1π(322)11(ν2q)3[i(i41)(1+i2)3+tan1(i)]5[(1+i2)9i15]

𝔟c

=

323i3(1+i2)6/5[i(i41)+(1+i2)3tan1(i)].


Gravitational Potential Energy of the Envelope

Again, pulling from our detailed derivations and appreciating, in particular, that (see, for example, our notes on equilibrium conditions),

A

=

ηisin(ηiB),

(ηsB)

=

π,

ηiB

=

π2tan1(Λi),

sin(ηiB)=(1+Λi2)1/2

     and    

sin[2(ηiB)]=2Λi(1+Λi2)1 ,

we have,

[WenvEnorm]eq

=

(123π)1/2(μeμc)3A2{[6(ηsB)3sin[2(ηsB)]4ηssin2(ηsB)+4B]

 

 

[6(ηiB)3sin[2(ηiB)]4ηisin2(ηiB)+4B]}

 

=

(123π)1/2(μeμc)3[ηisin(ηiB)]2{6π[6(ηiB)3sin[2(ηiB)]4ηisin2(ηiB)]}

 

=

(123π)1/2(μeμc)3ηi2(1+Λi2){6π6[π2tan1(Λi)]+6[Λi(1+Λi2)]+4ηi[1(1+Λi2)]}

 

=

(322π)1/2(μeμc)3ηi2{(1+Λi2)[π2+tan1(Λi)]+Λi+23ηi}.

So, in equilibrium we can write,

χeq[WenvEnorm]eq

=

(322π)1/2(μeμc)3ηi2{(1+Λi2)[π2+tan1(Λi)]+Λi+23ηi}(π23)1/2ν2q(1+i2)333i5

 

=

322(ηim3)3{(1+Λi2)ηi[π2+tan1(Λi)]+Λiηi+23}ν2q(1+i2)3i5

 

=

322(ν2q)1i2{(1+Λi2)ηi[π2+tan1(Λi)]+Λiηi+23}.

And out of equilibrium,

WenvEnorm

=

χ1322(ν2q)1i2[𝔎i+23].

This, in turn, implies that both in and out of equilibrium,

𝔞e

=

323[ν2(1+q)q(1ν)]1i2{(1+Λi2)ηi[π2+tan1(Λi)]+Λiηi+23}.

Thermal Energy of the Envelope

Again, pulling from our detailed derivations,

[SenvEnorm]eq

=

(125π)1/2(μeμc)3A2{[6(ηsB)3sin[2(ηsB)]][6(ηiB)3sin[2(ηiB)]]}

 

=

(125π)1/2(μeμc)3[ηisin(ηiB)]2{6π6(ηiB)+3sin[2(ηiB)]}

 

=

(125π)1/2(μeμc)3ηi2(1+Λi2){6[π2+tan1(Λi)]+6[Λi(1+Λi2)1]}

 

=

12(322π)1/2(μeμc)3ηi2{(1+Λi2)[π2+tan1(Λi)]+Λi}.

So, in equilibrium we can write,

χeq3[SenvEnorm]eq

=

12(322π)1/2(μeμc)3ηi2{(1+Λi2)[π2+tan1(Λi)]+Λi}[(π23)1/2ν2q(1+i2)333i5]3

 

=

(ν2q)3(32π2212)1/2(μeμc)3ηi3{(1+Λi2)ηi[π2+tan1(Λi)]+Λiηi}[(1+i2)939i15]

 

=

(ν2q)3(π2635)[(1+i2)6i12]{(1+Λi2)ηi[π2+tan1(Λi)]+Λiηi}.

And, out of equilibrium,

[SenvEnorm]eq

=

χ3(ν2q)3(π2635)[(1+i2)6i12]𝔎.

Combined in Equilibrium

Notice that, in combination,

[2Senv+WenvEnorm]eq

=

23(322π)1/2(μeμc)3ηi3

 

=

23(322π)1/2(μeμc)3[3(μeμc)i(1+i2)1]3

 

=

(236π)1/2[i3(1+i2)3].

Also, from above,

[2Score+WcoreEnorm]eq

=

(3825π)1/2[i(83i2)(1+i2)3]

 

=

+(236π)1/2[i3(1+i2)3].

So, in equilibrium, these terms from the core and envelope sum to zero, as they should.

Out of Equilibrium

And now, in combination out of equilibrium,

𝔊Enorm

=

(χχeq)1{[WcoreEnorm]eq+[WenvEnorm]eq}+(χχeq)3/5(2nc3)[ScoreEnorm]eq+(χχeq)3(2ne3)[SenvEnorm]eq.

Hence, quite generally out of equilibrium,

χ[𝔊Enorm]

=

χ1(χχeq)1{[WcoreEnorm]eq+[WenvEnorm]eq}35χ1(χχeq)3/5(103)[ScoreEnorm]eq3χ1(χχeq)3(23)[SenvEnorm]eq.

Let's see what the value of this derivative is if the dimensionless radius, χ, is set to the value that has been determined, via a detailed force-balanced analysis, to be the equilibrium radius, namely, χ=χeq. In this case, we have,

{χ[𝔊Enorm]}χχeq

=

χeq1{[WcoreEnorm]eq+[WenvEnorm]eq+2[ScoreEnorm]eq+2[SenvEnorm]eq}.

But, according to the virial theorem — and, as we have just demonstrated — the four terms inside the curly braces sum to zero. So this demonstrates that the derivative of our out-of-equilibrium free-energy expression does go to zero at the equilibrium radius, as it should!

Summary51

In summary, the desired out of equilibrium free-energy expression is,

𝔊Enorm

=

WcoreEnorm+WenvEnorm+(2nc3)ScoreEnorm+(2ne3)SenvEnorm

 

=

χ1(324)ν2q1i2[𝔏i83]χ1322(ν2q)1i2[𝔎i+23]

 

 

+(253)(322π)1/5[χ1(ν2q)1(1+i2)2]3/5(322)2𝔏i+(23)χ3(ν2q)3(π2635)[(1+i2)6i12]𝔎

 

=

(324)[χ1ν2q1i2][𝔏i+4𝔎i]+(322π)1/5(3523)[χ1(ν2q)1(1+i2)2]3/5𝔏i

 

 

+(π2536)[χ1(ν2q)(1+i2)2i4]3𝔎.

Or, in terms of the ratio,

Xχχeq,

and pulling from the above expressions,

[WcoreEnorm]eq

=

(3825π)1/2[i(i483i21)(1+i2)3+tan1(i)]

 

=

(3825π)1/2[i(1+i2)]3[𝔏i83]

[WenvEnorm]eq

=

(322π)1/2(μeμc)3ηi2{(1+Λi2)[π2+tan1(Λi)]+Λi+23ηi}

 

=

(3825π)1/2[i(1+i2)]3[4𝔎i+83]

[ScoreEnorm]eq

=

12(3825π)1/2[i(i41)(1+i2)3+tan1(i)]

 

=

12(3825π)1/2[i(1+i2)]3𝔏i

[SenvEnorm]eq

=

12(322π)1/2(μeμc)3ηi2{(1+Λi2)[π2+tan1(Λi)]+Λi}

 

=

12(3825π)1/2[i(1+i2)]3(4𝔎i),

we have the streamlined,

(25π36)1/2[(1+i2)i]3[𝔊Enorm]

=

+X3/5(5𝔏i)+X3(4𝔎i)X1(3𝔏i+12𝔎i)

or, better yet,

Out-of-Equilibrium, Free-Energy Expression for BiPolytropes with (nc,ne)=(5,1)

24(qi2ν2)χeq[𝔊Enorm]

=

X3/5(5𝔏i)+X3(4𝔎i)X1(3𝔏i+12𝔎i)


where,

𝔏i

(i41)i2+(1+i2)3i3tan1i,

𝔎i

Λiηi+(1+Λi2)ηi[π2+tan1Λi],

Λi

1ηii,

ηi

=

3(μeμc)[i(1+i2)].

From the accompanying Table 1 parameter values, we also can write,

1q

=

ηsηi=1+1ηi[π2+tan1Λi],

ν

=

iq(1+Λi2)1/2.

Radial Derivatives

𝔊*X

=

X8/5(3𝔏i)X4(12𝔎i)+X2(3𝔏i+12𝔎i)

2𝔊*X2

=

35[X13/5(8𝔏i)+X5(80𝔎i)X1(10𝔏i+40𝔎i)]

Consistent with our generic discussion of the stability of bipolytropes and the specific discussion of the stability of bipolytropes having (nc,ne)=(5,1), it can straightforwardly be shown that 𝔊/χ=0 is satisfied by setting X=1; that is, the equilibrium condition is,

χ=χeq

=

(π23)1/2ν2q(1+i2)333i5.

Furthermore, the equilibrium configuration is unstable whenever 2𝔊/χ2<0, that is, it is unstable whenever,

𝔏i𝔎i

>

20.

Table 1 of an accompanying chapter — and the red-dashed curve in the figure adjacent to that table — identifies some key properties of the model that marks the transition from stable to unstable configurations along equilibrium sequences that have various values of the mean-molecular weight ratio, μe/μc.

Focus on Zero-Zero Free-Energy Expression

Here, we will draw heavily from the following accompanying chapters:


From Detailed Force-Balance Models

Equilibrium Radius

First View

In an accompanying chapter we find,

P0Req4GMtot2

=

(323π)(νq3)2[q2+25q5(f1𝔉)]

where,

f

1+52(ρeρc)(1q21)+(ρeρc)2[1q51+52(11q2)],

𝔉

52(ρeρc)1q5[(2q2+3q3q5)+35(ρeρc)(1+5q25q3+q5)],

ρeρc

=

q3(1ν)ν(1q3).

Here, we prefer to normalize the equilibrium radius to Rnorm. So, let's replace the central pressure with its expression in terms of Kc. Specifically,

P0

=

Kcρcγc=Kc[3Mcore4πRi3]γc=Kc[3νMtot4πq3Req3](nc+1)/ncP0Pnorm=[34π(νq3)1χeq3](nc+1)/nc

Kc[3νMtot4πq3Req3](nc+1)/ncReq4GMtot2

=

(323π)(νq3)2[q2+25q5(f1𝔉)]

Req(nc3)/nc

=

(GKc)Mtot(nc1)/nc[3ν4πq3](nc+1)/nc(323π)(νq3)2[q2+25q5(f1𝔉)]

χeq(nc3)/nc[ReqRnorm](nc3)/nc

=

12(4π3)1/nc(νq3)(nc1)/nc[q2+25q5(f1𝔉)].

Or, in terms of γc,

χeq43γc

=

12(34π)1γc(νq3)2γc[q2+25q5(f1𝔉)].

Second View

Alternatively, from our derivation and discussion of analytic detailed force-balance models,

[R4GMtot2]P0

  = 

(323π)ν2g2q4,

where,

[g(ν,q)]2

1+(ρeρ0)[2(1ρeρ0)(1q)+ρeρ0(1q21)].

In order to show that this expression is the same as the other one, above, we need to show that,

(323π)(νq3)2[q2+25q5(f1𝔉)]

=

(323π)ν2g2q4

f1𝔉

=

52q3[g21]

 

=

52q3(ρeρ0)[2(1ρeρ0)(1q)+ρeρ0(1q21)]

 

=

52q5(ρeρ0){2(q2q3)+ρeρ0[13q2+2q3]}.

Let's see …

f1𝔉

=

52(ρeρc)(1q21)+(ρeρc)2[1q51+52(11q2)]52(ρeρc)1q5[(2q2+3q3q5)+35(ρeρc)(1+5q25q3+q5)]

 

=

52(ρeρc)(1q21)52(ρeρc)1q5[(2q2+3q3q5)]

 

 

52(ρeρc)1q5[35(ρeρc)(1+5q25q3+q5)]+(ρeρc)2[1q51+52(11q2)]

 

=

52(ρeρc)1q5{(q3q5)+(2q23q3+q5)}

 

 

+12(ρeρc)21q5[3(15q2+5q3q5)]+12(ρeρc)21q5[22q5+5(q5q3)]

 

=

52(ρeρc)1q5[(q3q5)+(2q23q3+q5)]+12(ρeρc)21q5[3(15q2+5q3q5)+22q5+5(q5q3)]

 

=

52(ρeρc)1q5[2q22q3]+52q5(ρeρc)2[13q2+2q3].

Q.E.D.

Hence, the equilibrium radius can also be written as,

χeq43γc

=

12(34π)1γc(νq3)2γcq2g2;

or, in terms of the polytropic index,

χeqnc3

=

12nc(4π3)(νq3)nc1(qg)2nc.

Gravitational Potential Energy

Also from our accompanying discussion, we have,

WgravEnorm

=

X1(35)(νq3)2q5[12nc(4π3)(νq3)nc1(qg)2nc]1/(nc3)f(ν,q)

 

=

X1(65)q5f[2nc(nc3)(34π)(νq3)(1nc)+2(nc3)bξnc]1/(nc3)

 

=

X1(65)q5f[(6π)(νq3)nc5bξnc]1/(nc3).

Internal Energy Components

First View

Before writing out the expressions for the internal energy of the core and of the envelope, we note from our separate detailed derivation that, in either case,

[Piχ3γPnorm]eqχ33γ

=

[(PiP0)(P0Pnorm)χ3]eq[χχeq]33γ

 

=

{(PiP0)[34π(νq3)]γcχ33γc}eqX33γ,

where, in equilibrium,

(PiP0)eq

=

1bξq2

bξq2

=

{25q3f+[125q3(1+𝔉)]}1

 

=

[1+25q3(f1𝔉)]1

So, copying from our accompanying detailed derivation, we have,

(𝔖AEnorm)core

=

4π/3(γc1){(PiP0)[34π(νq3)]γcχ33γc}eqX33γc{(P0Pic)[q3(3bξ5)q5]}

 

=

1(γc1)[(4π3)1γc(νq3)γcχeq33γc]X33γcq3[1(35)bξq2],

(𝔖AEnorm)env

=

4π/3(γe1){(PiP0)[34π(νq3)]γcχ33γc}eqX33γe{(1q3)+bξ(P0Pie)[25q5𝔉]}

 

=

1(γe1)[(4π3)1γc(νq3)γcχeq33γc]X33γe(PiP0){(1q3)+bξ(P0Pie)[25q5𝔉]}

 

=

1(γe1)[(4π3)1γc(νq3)γcχeq33γc]X33γe{(1bξq2)(1q3)+bξ[25q5𝔉]}

 

=

1(γe1)[(4π3)1γc(νq3)γcχeq33γc]X33γe(1q3){1[125(q31q3)𝔉]bξq2}.

Furthermore,

[(4π3)1γc(νq3)γcχeq33γc]

=

(34π)γc1(νq3)γc{χeq43γc}(33γc)/(43γc)

 

=

(34π)γc1(νq3)γc{12(34π)1γc(νq3)2γc[q2+25q5(f1𝔉)]}(33γc)/(43γc)

 

=

(34π)(γc1)/(43γc)(νq3)(65γc)(43γc){q22[1+25q3(f1𝔉)]}(33γc)/(43γc)

 

=

(34π)1/(nc3)(νq3)(nc5)(nc3){q22[1+25q3(f1𝔉)]}3/(nc3)

 

=

[(23π)(νq3)(nc5)bξ3]1/(nc3).

Hence, we have,

(𝔖AEnorm)core

=

nc[(4π3)1γc(νq3)γcχeq33γc]X3/ncq3[1(35)bξq2]

 

=

nc[(23π)(νq3)(nc5)bξ3]1/(nc3)X3/ncq3[1(35)bξq2],

(𝔖AEnorm)env

=

ne[(23π)(νq3)(nc5)bξ3]1/(nc3)X3/ne(1q3){1[125(q31q3)𝔉]bξq2}.


Second View

In our accompanying discussion of energies associated with detailed force balance models, we used the notation,

Π

(323π)GMtot2R4(νq3)2=Pnormχ4(323π)(νq3)2,

which allows us to rewrite the above quoted relationship between the central pressure and the radius of the bipolytrope as,

P0=Π(qg)2.

We also showed that, in equilibrium, the relationship between the central pressure and the interface pressure is,

P0=Pi+Πeqq2.

This means that, in equilibrium, the ratio of the interface pressure to the central pressure is,

(PiP0)eq

=

1Πeqq2P0=11g2,

or given that (see above),

52q3[g21]

=

f1𝔉

g2

=

1+25q3(f1𝔉),

we have,

(PiP0)eq

=

1Πeqq2P0=1[1+25q3(f1𝔉)]1.

This is exactly the pressure-ratio expression presented in our "first view" and unveils the notation association,

bξq2

1g2.

From our separate derivation, we have, in equilibrium,

𝔊core=(2nc3)Score

=

(2nc3)(4π5)Req3q5(5Pi2q2+Π)eq

 

=

(q5nc5)Req3(23π3)Πeq[52q2(PiΠ)eq+1]

 

=

(nc5)[Rnorm3Pnorm]χeq1(ν2q)[52q2(PiP0)eq(P0Π)eq+1]

[𝔊coreEnorm]eq

=

(nc5)(ν2q)[52q2(11g2)(q2g2)+1]χeq1

 

=

(nc2)(ν2q)[g235]{12nc(4π3)(νq3)nc1(qg)2nc}1/(nc3)

 

=

nc[1(35)1g2](12)(ν2q)g2{2nc(34π)(νq3)1nc(qg)2nc}1/(nc3)

 

=

nc[1(35)1g2]{2nc2(3nc)(34π)(νq3)1nc(νq3)2(nc3)q5(nc3)q2ncg2ncg2(nc3)}1/(nc3)

 

=

nc[1(35)1g2]{(23π)(νq3)nc5q3nc15g6}1/(nc3).

Finally, switching from the g notation to the bξ notation gives,

[𝔊coreEnorm]eq

=

nc[1(35)bξq2]{(23π)(νq3)nc5q3nc15bξ3q6}1/(nc3)

 

=

ncq3[1(35)bξq2]{(23π)(νq3)nc5bξ3}1/(nc3),

which, after setting X=1, precisely matches the above, "first view" expression. Also from our previous derivation, we can write,

𝔊env=(2ne3)Senv

=

2π(2ne3)Req3Πeq{(1q3)(PiΠ)eq+(ρeρ0)[(2q2+3q3q5)+35(ρeρ0)(1+5q25q3+q5)]}

 

=

2π(2ne3)Req3[Pnormχ4(323π)(νq3)2]eq{(1q3)q2(g21)+(25)q5𝔉}

 

=

[PnormRnorm3]ne2(ν2q4)(1q3){(g21)+25(q31q3)𝔉}χeq1

[𝔊envEnorm]eq

=

ne(1q3){(g21)+25(q31q3)𝔉}q22(νq3)2[12nc(4π3)(νq3)nc1(qg)2nc]1/(nc3)

 

=

ne(1q3){(g21)+25(q31q3)𝔉}[2[nc(nc3)](34π)(νq3)(1nc)+2(nc3)q2(nc3)2ncg2nc]1/(nc3)

 

=

ne(1q3){(g21)+25(q31q3)𝔉}[(23π)(νq3)nc5q6g2nc]1/(nc3).

And, finally, switching from the g notation to the bξ notation gives,

[𝔊envEnorm]eq

=

ne(1q3)(bξq2)1{1[125(q31q3)𝔉]bξq2}[(23π)(νq3)nc5q6(bξq2)nc]1/(nc3)

 

=

ne(1q3){1[125(q31q3)𝔉]bξq2}[(23π)(νq3)nc5q62(nc3)+2ncbξ3nc+nc]1/(nc3)

 

=

ne[(23π)(νq3)nc5bξ3]1/(nc3)(1q3){1[125(q31q3)𝔉]bξq2},

which, after setting X=1, precisely matches the above, "first view" expression.


Summary00

In summary, the desired out of equilibrium free-energy expression is,

𝔊Enorm

=

A0X3/nc+B0X3/neC0X1

where,

A0(𝔖coreEnorm)eq

=

ncbξ[(23π)(νq3)(nc5)bξnc]1/(nc3)q3[1(35)bξq2],

B0(𝔖envEnorm)eq

=

nebξ[(23π)(νq3)(nc5)bξnc]1/(nc3)(1q3){1[125(q31q3)𝔉]bξq2},

C0(WgravEnorm)eq

=

(65)q5f[(23π)(νq3)nc5bξnc]1/(nc3).

Or, in a more compact form,

𝔊*[(23π)(νq3)(nc5)bξnc]1/(nc3)[𝔊Enorm]

=

ncA1X3/nc+neB1X3/ne3C1X1

where,

A1

1bξ(q3)[1(35)bξq2],

B1

1bξ(1q3){1[125(q31q3)𝔉]bξq2},

C1

(25)q5f.

Let's examine the behavior of the first radial derivative.

𝔊*X

=

3X[A1X3/ncB1X3/ne+C1X1].

Let's see whether the sum of terms inside the square brackets is zero at the derived equilibrium radius, that is, when X=1 and, hence, when

χ=χeq

=

[12nc(4π3)(νq3)nc1(qg)2nc]1/(nc3)

 

=

[12nc(4π3)(νq3)nc1bξnc]1/(nc3).

C1A1B1

=

(25)q5f1bξ(q3)[1(35)bξq2]1bξ(1q3){1[125(q31q3)𝔉]bξq2}

 

=

(25)q5f1bξ{1[125(q31q3)𝔉]bξq2}+q3bξ{1[125(q31q3)𝔉]bξq2}q3bξ[1(35)bξq2]

 

=

(25)q5f1bξ+[125(q31q3)𝔉]q2+q3bξ[125(q31q3)𝔉]q5q3bξ+(35)q5

 

=

q2{(25)q3f1bξq2+[125(q31q3)𝔉](1q3)+(35)q3}

 

=

q2{(25)q3f[1+25q3(f1𝔉)]+[(1q3)25q3𝔉]+(35)q3}

 

=

q2{0}.

Q.E.D.

Even slightly better:

1q2[(π23)(νq3)(5nc)bξnc]1/(nc3)[𝔊Enorm]

=

ncA2X3/nc+neB2X3/ne3C2X1,

or, better yet,

Out-of-Equilibrium, Free-Energy Expression for BiPolytropes with Structural (nc,ne)=(0,0)

2(q2ν)2χeq[𝔊Enorm]

=

ncA2X3/nc+neB2X3/ne3C2X1

where, keeping in mind that,

1(bξq2)

=

[1+25q3(f1𝔉)],

we have,

A2A1q2

q3(bξq2)[1(35)bξq2]

 

=

q3{[1+25q3(f1𝔉)](35)}

 

=

25q3[1+q3(f1𝔉)],

B2B1q2

1(bξq2)(1q3){1[125(q31q3)𝔉]bξq2}

 

=

(1q3){1(bξq2)1+25(q31q3)𝔉}

 

=

(1q3){[1+25q3(f1𝔉)]1+25(q31q3)𝔉}

 

=

25q3{(1q3)(f1𝔉)+𝔉}

 

=

25q3{f[1+q3(f1𝔉)]}

 

=

25q3fA2,

C2C1q2

25q3f.

As before, the equilibrium system is dynamically unstable if 2𝔊/X2<0. We have deduced that the system is unstable if,

ne3[3nencne]

<

A2C2=1f[1+q3(f1𝔉)].


Overview

BiPolytrope51

Key Analytic Expressions

Out-of-Equilibrium, Free-Energy Expression for BiPolytropes with (nc,ne)=(5,1)

𝔊51*24(qν2)χeq[𝔊51Enorm]

=

1i2[X3/5(5𝔏i)+X3(4𝔎i)X1(3𝔏i+12𝔎i)]

where,

𝔏i

(i41)i2+(1+i2)3i3tan1i,

𝔎i

Λiηi+(1+Λi2)ηi[π2+tan1Λi],

Λi

1ηii,

ηi

=

3(μeμc)[i(1+i2)].

From the accompanying Table 1 parameter values, we also can write,

1q

=

ηsηi=1+1ηi[π2+tan1Λi],

ν

=

iq(1+Λi2)1/2.

Consistent with our generic discussion of the stability of bipolytropes and the specific discussion of the stability of bipolytropes having (nc,ne)=(5,1), it can straightforwardly be shown that 𝔊51*/χ=0 is satisfied by setting X=1; that is, the equilibrium condition is,

χ=χeq

=

(π23)1/2ν2q(1+i2)333i5

 

=

{(π3)22ncνnc1q3nc[(1+i2)6/5(3i2)]nc}1/(nc3),

where the last expression has been cast into a form that more clearly highlights overlap with the expression, below, for the equilibrium radius for zero-zero bipolytropes. Furthermore, the equilibrium configuration is unstable whenever,

[2𝔊51*χ2]X=1<0,

that is, it is unstable whenever,

𝔏i𝔎i

>

20.

Table 1 of an accompanying chapter — and the red-dashed curve in the figure adjacent to that table — identifies some key properties of the model that marks the transition from stable to unstable configurations along equilibrium sequences that have various values of the mean-molecular weight ratio, μe/μc.

Behavior of Equilibrium Sequence

Here we reprint Figure 1 from an accompanying chapter wherein the structure of five-one bipolytropes has been derived. It displays detailed force-balance sequences in the qν plane for a variety of choices of the ratio of mean-molecular-weights, μe/μc, as labeled.

Five-One Bipolytropic Equilibrium Sequences for Various ratios of the mean molecular weight
Five-One Bipolytropic Equilibrium Sequences for Various ratios of the mean molecular weight
Limiting Values

Each sequence begins (i=0) at the origin, that is, at (q,ν)=(0,0). As i, however, the sequences terminate at different coordinate locations, depending on the value of m33(μe/μc). In deriving the various limits, it will be useful to note that,

1ηi

=

(1+i2)m3i,

Λi

=

(1+i2)m3ii

 

=

1m3i+[(1m3)m3]i

 

=

1m3i[1(m31)i2]

 

=

(m31)im3[11(m31)i2],

1+Λi2

=

1+1m32i2[1+(1m3)i2]2

 

=

1m32i2{m32i2+[1+(1m3)i2]2}

 

=

1m32i2{1+(22m3+m32)i2+(1m3)2i4}

Examining the three relevant parameter regimes, we see that:

  • For μe/μc<13, that is, m3<1

tan1Λi|i

tan1[(1m3)m3]i

 

π2[m3(1m3)i]

1q|i

1+(1+i2)m3i[πm3(1m3)i]

 

m3+πim3

q|i

11+(πi/m3)0.

 

and

 

(νq)2

=

i21+Λi2

 

=

m32i4{1+(22m3+m32)i2+(1m3)2i4}1

 

=

m32{i4+(22m3+m32)i2+(1m3)2}1

νq|i

m31m3

ν|i

[m31m3]11+(πi/m3)0.

  • For μe/μc=13, that is, m3=1

tan1Λi

=

tan1(1i)

tan1Λi|i

1i

1q|i

1+(1+i2)i[π2+1i]

 

(π2)i

q|i

2πi0

 

and

 

(νq)

=

i(1+1/i2)1/2

ν|i

i(2πi)=2π0.63662


  • For μe/μc>13, that is, m3>1

tan1Λi|i

tan1[(m31m3)i]

 

π2+[m3(m31)i]

1q|i

1+(1+i2)m3i[m3(m31)i]

 

=

1+(1+1/i2)(m31)

 

1+1(m31)=m3(m31)

q|i

=

(m31)m3

 

and

 

(νq)2

=

i21+Λi2

 

=

m32i4{1+(22m3+m32)i2+(m31)2i4}1

 

=

m32{i4+(22m3+m32)i2+(m31)2}1

νq|i

m3m31

ν|i

=

m3m31[m31m3]1.

Summarizing:

  • For μe/μc<13, that is, m3<1       …       (q,ν)i=(0,0).


  • For μe/μc=13, that is, m3=1       …       (q,ν)i=(0,2π).


  • For μe/μc>13, that is, m3>1       …       (q,ν)i=[(m31)/m3,1].
Turning Points

Let's identify the location of two turning points along the ν(q) sequence — one defines qmax and the other identifies νmax. They occur, respectively, where,

dlnqdlni=0

      and      

dlnνdlni=0.

In deriving these expressions, we will use the relations,

dηidi

=

m3(1i2)(1+i2)2,

dΛidi

=

1m3i2[1i2(1m3)],

where,

m33(μeμc).


Given that,

q

=

{1+1ηi[π2+tan1Λi]}1,

we find,

dlnqdlni

=

iq(q2)ddi{1ηi[π2+tan1Λi]}

 

=

qi{1ηi2[π2+tan1Λi]dηidi+1ηi(1+Λi2)dΛidi}

 

=

qi{(1i2)m3i2[π2+tan1Λi]+(1+i2)m32i3(1+Λi2)[1i2(1m3)]}

 

=

qm32i2{m3i(1i2)[π2+tan1Λi]+(1+i2)(1+Λi2)[1i2(1m3)]}.

And, given that,

ν

=

iq(1+Λi2)1/2.

we find,

dlnνdlni

=

iν{q(1+Λi2)1/2+q(1+Λi2)1/2dlnqdlniiqΛi(1+Λi2)3/2dΛidi}

 

=

qiν(1+Λi2)1/2{1+dlnqdlni+Λim3i(1+Λi2)[1i2(1m3)]}

In summary, then, the qmax turning point occurs where,

0

=

(1+Λi2)[π2+tan1Λi]+(1+i2)m3i(1i2)[1i2(1m3)];

and the νmax turning point occurs where,

0

=

1+Λim3i(1+Λi2)[1i2(1m3)]+qi3(1i2)m3[π2+tan1Λi]+qi2m32(1+i2)(1+Λi2)[1i2(1m3)]

 

=

1+qi3(1i2)m3[π2+tan1Λi]+[Λim3i(1+Λi2)+qi2m32(1+i2)(1+Λi2)][1i2(1m3)]

 

=

1+qi3(1i2)m3[π2+tan1Λi]+1m3i[Λi(1+Λi2)+qi3m3(1+i2)(1+Λi2)][1i2(1m3)].

NOTE:  As we show above, for the special case of m3=1 — that is, when μe/μc=13, precisely — the equilibrium sequence (as i) intersects the q=0 axis at precisely the value, ν=2/π. As is illustrated graphically in Figure 1 of an accompanying chapter, no νmax turning point exists for values of m3>1.

For the record, we repeat, as well, that the transition from stable to dynamically unstable configurations occurs along the sequence when,

(i41)i2+(1+i2)3i3tan1i

=

20{Λiηi+(1+Λi2)ηi[π2+tan1Λi]}

 

=

20(1+Λi2)(1+i2)m3i{Λi(1+Λi2)+[π2+tan1Λi]}

m3i(i41)+m3(1+i2)3tan1i

=

20i2(1+Λi2)(1+i2){Λi(1+Λi2)+[π2+tan1Λi]}

m3i(i41)+m3(1+i2)3tan1i20i2(1+i2)

=

Λi+(1+Λi2)[π2+tan1Λi].


In order to clarify what equilibrium sequences do not have any turning points, let's examine how the qmax turning-point expression behaves as i.

(1+i2)(1+Λi2)[1i2(1m3)]

=

m3i(i21)[π2+tan1Λi]

(1+i2)m3i(i21)[1+i2(m31)]

=

{1+1m32i2[(m31)i21]2}{π2+[π21Λi+13Λi3+𝒪(Λi5)]}

(1+i2)i2(m31)m3i(i21)[1+1i2(m31)]

=

{1+(m31)2i2m32[11(m31)i2]2}1(Λi)[113Λi2+𝒪(Λi4)0]

(1+i2)i(i21)m3(m31)[1+1i2(m31)]

=

[12(m31)i2+m32(m31)2i2+1(m31)2i4]m3(m31)i[11(m31)i2]1{1m323(m31)2i2[11(m31)i2]2}

(1+1i2)[1+1i2(m31)][11(m31)i2]

=

(11i2)[12(m31)i2+m32(m31)2i2+1(m31)2i4]{1m323(m31)2i2[11(m31)i2]2}

The leading-order term is unity on both sides of this expression, so they cancel; let's see what results from keeping terms i2.

1i2[1+1(m31)1(m31)]

=

1i2[12(m31)+m32(m31)2m323(m31)2]

2

=

2(m31)+2m323(m31)2

6(m31)2

=

6(m31)+2m32

6m3212m3+6

=

6m3+6+2m32

m3

=

32.

We therefore conclude that the qmax turning point does not appear along any sequence for which,

m3

>

32

μeμc

>

12.


Five-One Bipolytrope Equilibrium Sequences in qν Plane

Full Sequences for Various μeμc

Magnified View with Turning Points and Stability Transition-Points Identified

Five-One Sequences

Graphical Depiction of Free-Energy Surface

Figure 1:   Free-Energy Surface for (nc,ne)=(5,1) and μeμc=1
Free-Energy surface for 5_1 bipolytrope
Free-Energy surface for 5_1 bipolytrope

Left Panel: The free energy (vertical, blue axis) is plotted as a function of the radial interface location, ξi (red axis), and the normalized configuration radius, Xχ/χeq (green axis). Right Panel: Same as the left panel, but animated in order to highlight undulations of the surface. The value of the free energy is indicated by color as well as by the height of the warped surface — red identifies the lowest depicted energies while blue identifies the highest depicted energies; these same colors have been projected down onto the z=0 plane to present a two-dimensional, color-contour plot. A multi-colored line segment drawn parallel to the ξi axis at the value, X=1, identifies the configuration's equilibrium radius for each value of the interface location. Equilibrium configurations marked in white lie at the bottom of the principal free-energy "valley" and are stable, while configurations marked in blue lie at the top of a free-energy "hill," indicating that they are unstable; the red dot identifies the location along this equilibrium sequence where the transition from stable to unstable configurations occurs.

For purposes of reproducibility, it is incumbent upon us to clarify how the values of the free energy were normalized in order to produce the free-energy surface displayed in Figure 1. The normalization steps are explicitly detailed within the fortran program, below, that generated the data for plotting purposes; here we provide a brief summary. We evaluated the normalized free energy, 𝔊51*, across a 200×200 zone grid of uniform spacing, covering the following (x,y)=(i,X) domain:

13

i

33

0.469230769

X

2.0

(With this specific definition of the y-coordinate grid, X=1 is associated with zone 70.) After this evaluation, a constant, Efudge=10, was added to 𝔊* in order to ensure that the free energy was negative across the entire domain. Then (inorm = 1), for each specified interface location, x=i, employing the equilibrium value of the free energy,

E0=𝔊51*(i,X=1)+Efudge,

the free energy was normalized across all values of y=X via the expression,

fe=(𝔊51*+Efudge)(E0)i|E0|i.

Finally, for plotting purposes, at each grid cell vertex ("vertex") — as well as at each grid cell center ("cell") — the value of the free energy, fe, was renormalized as follows,

vertex=femin(fe)max(fe)min(fe).

Via this last step, the minimum "vertex" energy across the entire domain was 0.0 while the maximum "vertex" energy was 1.0.


FORTRAN Program Documentation Example Evaluations

(See also associated Table 1)
Coord. Axis Coord. Name Associated Physical Quantity μeμc=1 μeμc=0.305
x-axis bsize iξi3 2.4163=1.395 8.19383=4.7307 14.3893=8.3076
y-axis csize Xχχeq 1 1 1
Relevant Lines of Code  
      eta = 3.0d0*muratio*bsize/(1.0d0+bsize**2)
      Gami = 1.0d0/eta-bsize
      frakL = (bsize**4-1.0d0)/bsize**2 + &
     &        DATAN(bsize)*((1.0d0+bsize**2)/bsize)**3
      frakK = Gami/eta + ((1.0d0+Gami**2)/eta)*(pii/2.0d0+DATAN(Gami))
      E0    = ((5.0d0*frakL) + (4.0d0*frakK)&
     &        - (3.0d0*frakL+12.0d0*frakK))/bsize**2+Efudge
          fescalar(j,k) = (csize**(-0.6d0)*(5.0d0*frakL)&
     &        + csize**(-3.0d0)*(4.0d0*frakK)&
     &   - (3.0d0*frakL+12.0d0*frakK)/csize)/bsize**2 + Efudge
          if(inorm.eq.1)fescalar(j,k)=fescalar(j,k)/DABS(E0) &
     &                  - E0/DABS(E0)
Variable Represents Value calculated via the expression …
eta ηi

3(μeμc)[i(1+i2)]

1.421 0.1851 0.1086
Gami Λi 1ηii 0.691 0.6705 0.9033
frakL 𝔏i (i41)i2+[1+i2i]3tan1i 10.37 186.80 937.64
frakK 𝔎i Λiηi+(1+Λi2)ηi[π2+tan1Λi] 0.518 20.544 46.882
  𝔏i𝔎i   20 9.093 20
E0 - Efudge 𝔊51*(i,X=1)

1i2[5𝔏i+4𝔎i(3𝔏i+12𝔎i)]=2(𝔏i4𝔎i)i2

8.525 9.3496 21.737
Figure 2:   Free-Energy Surface for (nc,ne)=(5,1) and μeμc=0.305
Free-Energy surface for 5_1 bipolytrope
Free-Energy surface for 5_1 bipolytrope

Left Panel: The free energy (vertical, blue axis) is plotted as a function of the radial interface location, ξi (red axis), and the normalized configuration radius, Xχ/χeq (green axis). Right Panel: Same as the left panel, but animated in order to highlight undulations of the surface. The value of the free energy is indicated by color as well as by the height of the warped surface — red identifies the lowest depicted energies while blue identifies the highest depicted energies; these same colors have been projected down onto the z=0 plane to present a two-dimensional, color-contour plot. A multi-colored line segment drawn parallel to the ξi axis at the value, X=1, identifies the configuration's equilibrium radius for each value of the interface location. Equilibrium configurations marked in white lie at the bottom of the principal free-energy "valley" and are stable, while configurations marked in blue lie at the top of a free-energy "hill," indicating that they are unstable; the red dot identifies the location along this equilibrium sequence where the transition from stable to unstable configurations occurs.

BiPolytrope00

Out-of-Equilibrium, Free-Energy Expression for BiPolytropes with Structural (nc,ne)=(0,0)

𝔊00*5(qν2)χeq[𝔊00Enorm]

=

52q3[ncA2X3/nc+neB2X3/ne3C2X1]

where,

A2

25q3[1+q3(f1𝔉)],

B2

25q3fA2,

C2

25q3f,

f

1+52(ρeρc)(1q21)+(ρeρc)2[1q51+52(11q2)],

𝔉

52(ρeρc)1q5[(2q2+3q3q5)+35(ρeρc)(1+5q25q3+q5)],

ρeρc

=

q3(1ν)ν(1q3).

The associated equilibrium radius is,

χeq

=

{(π3)22ncνnc1q3nc[1+25q3(f1𝔉)]nc}1/(nc3).

We have deduced that the system is unstable if,

ne3[3nencne]

<

A2C2=1f[1+q3(f1𝔉)].

Fortran Code

This is the program that generated the free-energy data for the "five-one" bipolytrope that is displayed in the above, Figure 1 image/animation.

      PROGRAM BiPolytrope
      real*8 pii
      real*8 bmin,bmax,cmin,cmax,db,dc
      real*8 c(200),b(200),chalf(199),bhalf(199)
      real*8 bsize,csize,emin,emax
      real*8 fepoint(200,200),fescalar(199,199)
      real*8 ell(200),ellhalf(199)
      real*8 muratio,eta,Gami,frakK,frakL
      real*8 eshift,ediff
      real   xx(200),yy(200),cell(199,199),vertex(200,200)
      real   valuemin,minlog,valufudge
      real*8 q,nu,chiEq,Enorm,E0,Efudge
      integer j,k,n,nmax,inorm
  101 format(4x,'bsize',7x,'csize',8x,'xi',10x,'A',12x,'B',12x,&
     &'fM',13x,'fW',11x,'fA',/)
! 102 format(1p8d12.4)
  103 format(2i5,1p3d14.6)
  201 format(5x,'valuemin = ',1pe15.5,//)
  205 format(5x,'For Cell-center ... emin, emax = ',1p2d14.6,/)
  206 format(5x,'For Cell-vertex ... emin, emax = ',1p2d14.6,/)
  701 format(5x,1p10d10.2)  
  700 format(8x,'xi',9x,'ell',8x,'eta',8x,'Lambda',5x,'frakK',&
     &       5x,'frakL',8x,'q',5x,'nu',5x,'chiEq',8x,'E0',/)
!!!!!!!!
!!!!!!!!
      inorm=1
      pii = 4.0d0*datan(1.0d0) 
      muratio = 1.0d0
      bsize = 0.0d0
      csize = 0.0d0
      Efudge = -10.0d0
      write(*,101)
!     write(*,102)bsize,csize,xival,coefA,coefB,fM,fW,fA
!!!!!!!!!!!
!
!  In this free-energy routine, c = X = chi/chi_eq and b = xi_i
!
!!!!!!!!!!!
      nmax = 200
      bmin = 1.0d0
      bmax = 3.0d0
      db   = (bmax-bmin)/dfloat(nmax-1)
      b(1) = bmin
      ell(1) = b(1)/dsqrt(3.0d0)
!  These values of cmin and cmax ensure that X=1 occurs at zone 70
      cmin = 0.469230769d0
      cmax = 2.00d0
      dc   = (cmax-cmin)/dfloat(nmax-1)
      c(1) = cmin
      do n=2,nmax
        b(n) = b(n-1)+db
        c(n) = c(n-1)+dc
        ell(n) = b(n)/dsqrt(3.0d0)
      enddo
      do n=1,nmax-1
        bhalf(n) = 0.5d0*(b(n)+b(n+1))
        chalf(n) = 0.5d0*(c(n)+c(n+1))
        ellhalf(n) = bhalf(n)/dsqrt(3.0d0)
      enddo
!
!  BEGIN LOOP to evaluate free energy (cell center)
!
      emin = 0.0d0
      emax = 0.0d0
      write(*,700)
      do j=1,nmax-1
        bsize = ellhalf(j)
      eta = 3.0d0*muratio*bsize/(1.0d0+bsize**2)
      Gami = 1.0d0/eta-bsize
      frakL = (bsize**4-1.0d0)/bsize**2 + &
     &        DATAN(bsize)*((1.0d0+bsize**2)/bsize)**3
      frakK = Gami/eta + ((1.0d0+Gami**2)/eta)*(pii/2.0d0+DATAN(Gami))
      q = 1.0d0/(1.0d0 + (0.5d0*pii+DATAN(Gami))/eta)
      nu = bsize*q/dsqrt(1.0d0+Gami**2)
      chiEq = dsqrt(pii/8.0d0)*(nu**2/(q*bsize**2))&
     &        *((1.0d0+bsize**2)/(3.0d0*bsize))**3
      Enorm = 16.0d0*(q/nu**2)*chiEq
      E0    = ((5.0d0*frakL) + (4.0d0*frakK)&
     &        - (3.0d0*frakL+12.0d0*frakK))/bsize**2+Efudge
      write(*,701)b(j),bsize,eta,Gami,frakK,frakL,q,nu,chiEq,E0
        do k=1,nmax-1
          csize=chalf(k)
          fescalar(j,k) = (csize**(-0.6d0)*(5.0d0*frakL)&
     &        + csize**(-3.0d0)*(4.0d0*frakK)&
     &   - (3.0d0*frakL+12.0d0*frakK)/csize)/bsize**2 + Efudge
          if(inorm.eq.1)fescalar(j,k)=fescalar(j,k)/DABS(E0) &
     &                  - E0/DABS(E0)
        if(fescalar(j,k).gt.0.5d0)fescalar(j,k)=0.5d0
          if(fescalar(j,k).lt.emin)emin=fescalar(j,k)
          if(fescalar(j,k).gt.emax)emax=fescalar(j,k)
!         write(*,103)j,k,bsize,csize,fescalar(j,k)
        enddo
      enddo
      write(*,205)emin,emax
!
!  BEGIN LOOP to evaluate free energy (cell vertex)
!
      emin = 0.0d0
      emax = 0.0d0
      do j=1,nmax
        bsize = ell(j)
      eta = 3.0d0*muratio*bsize/(1.0d0+bsize**2)
      Gami = 1.0d0/eta-bsize
      frakL = (bsize**4-1.0d0)/bsize**2 + &
     &        DATAN(bsize)*((1.0d0+bsize**2)/bsize)**3
      frakK = Gami/eta + ((1.0d0+Gami**2)/eta)*(pii/2.0d0+DATAN(Gami))
      q = 1.0d0/(1.0d0 + (0.5d0*pii+DATAN(Gami))/eta)
      nu = bsize*q/dsqrt(1.0d0+Gami**2)
      chiEq = dsqrt(pii/8.0d0)*(nu**2/(q*bsize**2))&
     &        *((1.0d0+bsize**2)/(3.0d0*bsize))**3
      Enorm = 16.0d0*(q/nu**2)*chiEq
      E0    = ((5.0d0*frakL) + (4.0d0*frakK)&
     &        - (3.0d0*frakL+12.0d0*frakK))/bsize**2 + Efudge
        do k=1,nmax
          csize=c(k)
          fepoint(j,k) = (csize**(-0.6d0)*(5.0d0*frakL)&
     &        + csize**(-3.0d0)*(4.0d0*frakK)&
     &   - (3.0d0*frakL+12.0d0*frakK)/csize)/bsize**2 + Efudge
          if(inorm.eq.1)fepoint(j,k)=fepoint(j,k)/DABS(E0) &
     &                  - E0/DABS(E0)
        if(fepoint(j,k).gt.0.5d0)fepoint(j,k)=0.5d0
          if(fepoint(j,k).lt.emin)emin=fepoint(j,k)
          if(fepoint(j,k).gt.emax)emax=fepoint(j,k)
!         write(*,103)j,k,bsize,csize,fepoint(j,k)
        enddo
      enddo
      write(*,206)emin,emax
!
!  Now fill single-precision arrays for plotting.
!
      do n=1,nmax
!       xx(n)=b(n)/b(nmax)
!       yy(n)=c(n)/c(nmax)
        xx(n)=b(n)-bmin
        yy(n)=c(n)-cmin
!       xx(n)=b(n)
!       yy(n)=c(n)
      enddo
      valuemin= -emin
      valufudge = 1.0d0/(emax-emin)
      do k=1,nmax
        do j=1,nmax
          vertex(j,k)=fepoint(j,k)+valuemin
          vertex(j,k)=vertex(j,k)*valufudge
        enddo
      enddo
      do k=1,nmax-1
        do j=1,nmax-1
          cell(j,k)=fescalar(j,k)+valuemin
          cell(j,k)=cell(j,k)*valufudge
        enddo
      enddo
      call XMLwriter01(nmax,xx,yy,cell,vertex)
      stop
      END PROGRAM BiPolytrope
      Subroutine XMLwriter01(imax,x,y,cell_scalar,point_scalar)
      real x(200),y(200),z(1)
      real cell_scalar(199,199),point_scalar(200,200)
      integer imax
      integer extentX,extentY,extentZ
      integer ix0,iy0,iz0
      integer norm(200,3)
!     imax=200
      ix0=0
      iy0=0
      iz0=0
      extentX=imax-1
      extentY=imax-1
      extentZ=0
      z(1) = 0.0
! Set normal vector 1D array
      do i=1,imax
        norm(i,1)=0
        norm(i,2)=0
        norm(i,3)=1
      enddo
  201 format('<?xml version="1.0"?>')
  202 format('<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian">')
  302 format('</VTKFile>')
  203 format(2x,'<RectilinearGrid WholeExtent="',6I4,'">')
  303 format(2x,'</RectilinearGrid>')
  204 format(4x,'<Piece Extent="',6I4,'">')
  304 format(4x,'</Piece>')
  205 format(6x,'<CellData Scalars="cell_scalars" Normals="magnify">')
  305 format(6x,'</CellData>')
  206 format(8x,'<DataArray type="Float32" Name="magnify" NumberOfComponents="3" format="ascii">')
  207 format(8x,'<DataArray type="Float32" Name="cell_scalars" format="ascii">')
  399 format(8x,'</DataArray>')
  208 format(6x,'<PointData Scalars="colorful" Normals="direction">')
  308 format(6x,'</PointData>')
  209 format(8x,'<DataArray type="Float32" Name="colorful" format="ascii">')
  210 format(6x,'<Coordinates>')
  310 format(6x,'</Coordinates>')
  211 format(8x,'<DataArray type="Float32" format="ascii" RangeMin="0" RangeMax="5">')
  212 format(8x,'<DataArray type="Float32" format="ascii">')
  213 format(8x,'<DataArray type="Float32" Name="direction" NumberOfComponents="3" format="ascii">')
  501 format(10f9.5)
  502 format(10f9.5)
  503 format(5x,9(1x,3I2))
  504 format(10f9.5)
  505 format(5x,10(1x,3I2))
!!!!!
!
!  Begin writing out XML tags.
!
!!!!!
      write(*,201)                              !<?xml
      write(*,202)                              !VTKFile
        write(*,203)ix0,extentX,iy0,extentY,iz0,extentZ        !  RectilinearGrid
          write(*,204)ix0,extentX,iy0,extentY,iz0,extentZ      !    Piece
            write(*,205)                        !      CellData
              write(*,207)                      !        DataArray(cell_scalars)
      do j=1,imax-1
      write(*,501)(cell_scalar(i,j),i=1,imax-1)
      enddo
              write(*,399)                      !        /DataArray
              write(*,206)                      !        DataArray(cell_scalars)
      do j=1,imax-1
      write(*,503)(norm(i,1),norm(i,2),norm(i,3),i=1,imax-1)
      enddo
              write(*,399)                      !        /DataArray
            write(*,305)                        !      /CellData
            write(*,208)                        !      PointData
              write(*,209)                      !        DataArray(points)
      write(*,502)((point_scalar(i,j),i=1,imax),j=1,imax)
              write(*,399)                      !        /DataArray
              write(*,213)                      !        DataArray(cell_scalars)
      do j=1,imax
      write(*,505)(norm(i,1),norm(i,2),norm(i,3),i=1,imax)
      enddo
              write(*,399)                      !        /DataArray
            write(*,308)                        !      /PointData
            write(*,210)                        !      Coordinates
              write(*,212)                      !        DataArray(x-direction)
        write(*,504)(x(i),i=1,imax)
              write(*,399)                      !        /DataArray
              write(*,212)                      !        DataArray(y-direction)
        write(*,504)(y(i),i=1,imax)
              write(*,399)                      !        /DataArray
              write(*,212)                      !        DataArray(z-direction)
        write(*,504)z(1)
              write(*,399)                      !        /DataArray
            write(*,310)                        !      /Coordinates
          write(*,304)                          !    /Piece
        write(*,303)                            !  /RectilinearGrid
      write(*,302)                              !/VTKFile
      return
      end

Nonstandard Examination

In our introductory remarks, above, we said the warped free-energy surface drapes across a five-dimensional parameter "plane" such that,

𝔊51

=

𝔊(R,Kc,Mtot,q,ν).

From a more pragmatic point of view, we should have said that the "five-one" free-energy surface drapes across a five-dimensional parameter "plane" such that,

𝔊51

=

𝔊(R,Kc,Mtot,i,μeμc).

In our initial, standard examination of the structure of this warped free-energy surface, we held three parameters fixed — namely, Kc,Mtot,μeμc — and plotted 𝔊51(i,XR/Req). Now, let's fix X=1 and plot 𝔊51(i,μeμc). The following plot shows how a portion of the (i,μe/μc) grid maps onto the traditional (q,ν) plane. The numerical labels identify lines of constant ξi=3i — 7 (light green), 9 (pink), and 12 (orange) — and lines of constant μe/μc — 0.330 (purple), 0.315 (dark green), and 0.305 (white/blue).

xi-ell grid drawn on q-nu grid
xi-ell grid drawn on q-nu grid

See Also


Tiled Menu

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