Appendix/Ramblings/BiPolytropeStability: Difference between revisions

From jetwiki
Jump to navigation Jump to search
Line 503: Line 503:
<tr>
<tr>
   <td align="right">
   <td align="right">
<math>~0</math>
<math>0</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
<math>~=</math>
<math>=</math>
   </td>
   </td>
   <td align="left">
   <td align="left">
<math>~
<math>
\frac{d^2x}{d\xi^2} + \biggl[ 4 - 2 Q_1 \biggr] \frac{1}{\xi} \cdot \frac{dx}{d\xi}   
\frac{d^2x}{d\xi^2} + \biggl[ 4 - 2 Q_1 \biggr] \frac{1}{\xi} \cdot \frac{dx}{d\xi}   
+ 2 \biggl[ \biggl( \frac{\sigma_c^2}{6\gamma_\mathrm{core} } \biggr)  \frac{\xi^2}{\theta}  
+ 2 \biggl[ \biggl( \frac{\sigma_c^2}{6\gamma_\mathrm{core} } \biggr)  \frac{\xi^2}{\theta}  
Line 522: Line 522:
<tr>
<tr>
   <td align="right">
   <td align="right">
<math>~Q_1</math>
<math>Q_1</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
<math>~\equiv</math>
<math>\equiv</math>
   </td>
   </td>
   <td align="left">
   <td align="left">
<math>~- \frac{d\ln\theta}{d\ln\xi} \, .</math>
<math>- \frac{d\ln\theta}{d\ln\xi} \, .</math>
   </td>
   </td>
</tr>
</tr>
Line 548: Line 548:
<tr>
<tr>
   <td align="right">
   <td align="right">
<math>~Q_1</math>
<math>Q_1</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
<math>~=</math>
<math>=</math>
   </td>
   </td>
   <td align="left">
   <td align="left">
<math>~
<math>
- \frac{\xi^2}{\sin\xi} \biggl[ \frac{\cos\xi}{\xi}- \frac{\sin\xi}{\xi^2}\biggr]
- \frac{\xi^2}{\sin\xi} \biggl[ \frac{\cos\xi}{\xi}- \frac{\sin\xi}{\xi^2}\biggr]
</math>
</math>
Line 565: Line 565:
   </td>
   </td>
   <td align="center">
   <td align="center">
<math>~=</math>
<math>=</math>
   </td>
   </td>
   <td align="left">
   <td align="left">
<math>~
<math>
1 -  \xi\cot\xi \, .
1 -  \xi\cot\xi \, .
</math>
</math>
Line 581: Line 581:
<tr>
<tr>
   <td align="right">
   <td align="right">
<math>~0</math>
<math>0</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
<math>~=</math>
<math>=</math>
   </td>
   </td>
   <td align="left">
   <td align="left">
<math>~
<math>
\frac{d^2x}{d\xi^2} + \biggl[ 4 - 2 ( 1 -  \xi\cot\xi ) \biggr] \frac{1}{\xi} \cdot \frac{dx}{d\xi}   
\frac{d^2x}{d\xi^2} + \biggl[ 4 - 2 ( 1 -  \xi\cot\xi ) \biggr] \frac{1}{\xi} \cdot \frac{dx}{d\xi}   
+ 2 \biggl[ \biggl( \frac{\sigma_c^2}{6\gamma_\mathrm{core} } \biggr)  \frac{\xi^3}{\sin\xi}  
+ 2 \biggl[ \biggl( \frac{\sigma_c^2}{6\gamma_\mathrm{core} } \biggr)  \frac{\xi^3}{\sin\xi}  
Line 600: Line 600:
   </td>
   </td>
   <td align="center">
   <td align="center">
<math>~=</math>
<math>=</math>
   </td>
   </td>
   <td align="left">
   <td align="left">
<math>~
<math>
\frac{d^2x}{d\xi^2} + \biggl[ 1 +  \xi\cot\xi  \biggr] \frac{2}{\xi} \cdot \frac{dx}{d\xi}   
\frac{d^2x}{d\xi^2} + \biggl[ 1 +  \xi\cot\xi  \biggr] \frac{2}{\xi} \cdot \frac{dx}{d\xi}   
+ 2 \biggl[ \biggl( \frac{\sigma_c^2}{6\gamma_\mathrm{core} } \biggr)  \frac{\xi^3}{\sin\xi}  
+ 2 \biggl[ \biggl( \frac{\sigma_c^2}{6\gamma_\mathrm{core} } \biggr)  \frac{\xi^3}{\sin\xi}  
Line 617: Line 617:
<tr>
<tr>
   <td align="right">
   <td align="right">
<math>~0</math>
<math>0</math>
   </td>
   </td>
   <td align="center">
   <td align="center">
<math>~=</math>
<math>=</math>
   </td>
   </td>
   <td align="left">
   <td align="left">
<math>~
<math>
\frac{d^2x}{d\xi^2} +  \frac{2}{\xi} \biggl[ 1 +  \xi\cot\xi  \biggr]\frac{dx}{d\xi}   
\frac{d^2x}{d\xi^2} +  \frac{2}{\xi} \biggl[ 1 +  \xi\cot\xi  \biggr]\frac{dx}{d\xi}   
+ \biggl[ \biggl( \frac{\sigma_c^2}{3\gamma_\mathrm{core} } \biggr)  \frac{\xi}{\sin\xi}  
+ \biggl[ \biggl( \frac{\sigma_c^2}{3\gamma_\mathrm{core} } \biggr)  \frac{\xi}{\sin\xi}  
Line 636: Line 636:
   </td>
   </td>
   <td align="center">
   <td align="center">
<math>~=</math>
<math>=</math>
   </td>
   </td>
   <td align="left">
   <td align="left">
<math>~
<math>
\frac{d^2x}{d\xi^2} +  \frac{2}{\xi} \biggl[ 1 +  \xi\cot\xi  \biggr]\frac{dx}{d\xi}   
\frac{d^2x}{d\xi^2} +  \frac{2}{\xi} \biggl[ 1 +  \xi\cot\xi  \biggr]\frac{dx}{d\xi}   
+ \biggl[ \frac{\gamma_g}{\gamma_\mathrm{core}}\biggl( \omega_k^2 \theta_c \biggr)  \frac{\xi}{\sin\xi}  
+ \biggl[ \frac{\gamma_g}{\gamma_\mathrm{core}}\biggl( \omega_k^2 \theta_c \biggr)  \frac{\xi}{\sin\xi}  
Line 647: Line 647:
</tr>
</tr>
</table>
</table>
which matches the expression presented by [http://adsabs.harvard.edu/abs/1985PASAu...6..222M Murphy &amp; Fiedler (1985b)] (see middle of the left column on p. 223 of their article) if we set <math>~\theta_c = 1</math> and <math>~\gamma_g/\gamma_\mathrm{core} = 1</math>.  This LAWE also appears in our [[User:Tohline/SSC/Stability/n1PolytropeLAWE#MurphyFiedler1985b|separate discussion of radial oscillations in n = 1 polytropic spheres]].
which matches the expression presented by [http://adsabs.harvard.edu/abs/1985PASAu...6..222M Murphy &amp; Fiedler (1985b)] (see middle of the left column on p. 223 of their article) if we set <math>~\theta_c = 1</math> and <math>~\gamma_g/\gamma_\mathrm{core} = 1</math>.  This LAWE also appears in our [[SSC/Stability/n1PolytropeLAWE#MurphyFiedler1985b|separate discussion of radial oscillations in n = 1 polytropic spheres]].


==Interface Conditions==
==Interface Conditions==

Revision as of 16:54, 12 April 2022

Ramblings: Marginally Unstable Bipolytropes

Our aim is to determine whether or not there is a relationship between (1) equilibrium models at turning points along bipolytrope sequences and (2) bipolytropic models that are marginally (dynamically) unstable toward collapse (or dynamical expansion).

Overview

file = Dropbox/WorkFolder/Wiki edits/EmbeddedPolytropes/CombinedSequences.xlsx --- worksheet = EqSeqCombined2
file = Dropbox/WorkFolder/Wiki edits/EmbeddedPolytropes/CombinedSequences.xlsx --- worksheet = EqSeqCombined2
Figure 1:  Equilibrium Sequences
of Pressure-Truncated Polytropes

Equilibrium sequences of Pressure-Truncated Polytropes

We expect the content of this chapter — which examines the relative stability of bipolytropes — to parallel in many ways the content of an accompanying chapter in which we have successfully analyzed the relative stability of pressure-truncated polytopes. Figure 1, shown here on the right, has been copied from a closely related discussion. The curves show the mass-radius relationship for pressure-truncated model sequences having a variety of polytropic indexes, as labeled, over the range 1n6. (Another version of this figure includes the isothermal sequence.) On each sequence for which n3, the green filled circle identifies the model with the largest mass. We have shown analytically that the oscillation frequency of the fundamental-mode of radial oscillation is precisely zero for each one of these maximum-mass models. As a consequence, we know that each green circular marker identifies the point along its associated sequence that separates dynamically stable (larger radii) from dynamically unstable (smaller radii) models.


In each case, the fundamental-mode oscillation frequency is precisely zero if, and only if, the adiabatic index governing expansions/contractions is related to the underlying structural polytropic index via the relation, γg=(n+1)/n, and if a constant surface-pressure boundary condition is imposed.



In another accompanying chapter, we have used purely analytic techniques to construct equilibrium sequences of spherically symmetric bipolytropes that have, (nc,ne)=(5,1). For a given choice of μe/μc — the ratio of the mean-molecular weight of envelope material to the mean-molecular weight of material in the core — a physically relevant sequence of models can be constructed by steadily increasing the value of the dimensionless radius at the core/envelope interface, ξi, from zero to infinity. Figure 2, which has been copied from this separate chapter, shows how the fractional core mass, νMcore/Mtot, varies with the fractional core radius, qrcore/R, along sequences having six different values of μe/μc: 1 (blue diamonds), ½ (red squares), 0.345 (dark purple crosses), ⅓ (pink triangles), 0.309 (light green dashes), and ¼ (purple asterisks). Along each of the model sequences, points marked by solid-colored circles correspond to models whose interface parameter, ξi, has one of three values: 0.5 (green circles), 1 (dark blue circles), or 3 (orange circles).

When modeling bipolytropes, the default expectation is that an increase in ξi along a given sequence will correspond to an increase in the relative size — both the radius and the mass — of the core. As Figure 2 illustrates, this expectation is realized along the sequences marked by blue diamonds (μe/μc=1) and by red squares (μe/μc=½). But the behavior is different along the other four illustrated sequences. For sufficiently large ξi, the relative radius of the core begins to decrease. Furthermore, along sequences for which μe/μc<13, eventually the fractional mass of the core reaches a maximum and, thereafter, decreases even as the value of ξi continues to increase. (Additional properties of these equilibrium sequences are discussed in yet another accompanying chapter.)

The principal question is: Along bipolytropic sequences, are maximum-mass models associated with the onset of dynamical instabilities?

Planned Approach

Figure 2: Equilibrium Sequences of Bipolytropes with (nc,ne)=(5,1)

Ideally we would like to answer the just-stated "principal question" using purely analytic techniques. But, to date, we have been unable to fully address the relevant issues analytically, even in what would be expected to be the simplest case:   bipolytropic models that have (nc,ne)=(0,0). Instead, we will streamline the investigation a bit and proceed — at least initially — using a blend of techniques. We will investigate the relative stability of bipolytropic models having (nc,ne)=(5,1) whose equilibrium structures are completely defined analytically; then the eigenvectors describing radial modes of oscillation will be determined, one at a time, by solving the relevant LAWE(s) numerically. We are optimistic that this can be successfully accomplished because we have had experience numerically integrating the LAWE that governs the oscillation of:

A key reference throughout this investigation will be the paper by J. O. Murphy & R. Fiedler (1985b, Proc. Astr. Soc. of Australia, 6, 222). They studied Radial Pulsations and Vibrational Stability of a Sequence of Two Zone Polytropic Stellar Models. Specifically, their underlying equilibrium models were bipolytropes that have (nc,ne)=(1,5). In an accompanying chapter, we describe in detail how Murphy & Fiedler obtained these equilibrium bipolytropic structures and detail some of their equilibrium properties.

Here are the steps we initially plan to take:

  • Governing LAWEs:
  • Determine what surface boundary condition should be imposed on physically relevant LAWE solutions, i.e., on the physically relevant radial-oscillation eigenvectors.
  • Initial Analysis:
    • Choose a maximum-mass model along the bipolytropic sequence that has, for example, μe/μc=1/4. Hopefully, we will be able to identify precisely (analytically) where this maximum-mass model lies along the sequence. Yes! Our earlier analysis does provide an analytic prescription of the model that sits at the maximum-mass location along the chosen sequence.
    • Solve the relevant eigenvalue problem for this specific model, initially for (γc,γe)=(6/5,2) and initially for the fundamental mode of oscillation.

Review of the Analysis by Murphy & Fiedler (1985b)

In the stability analysis presented by Murphy & Fiedler (1985b), the relevant polytropic indexes are, (nc,ne)=(1,5). Structural properties of the underlying equilibrium models have been reviewed in our accompanying discussion.

The Linear Adiabatic Wave Equation (LAWE) that is relevant to polytropic spheres may be written as,

0=d2xdξ2+[4(n+1)Q]1ξdxdξ+(n+1)[(σc26γg)ξ2θαQ]xξ2

where:    Q(ξ)dlnθdlnξ,    σc23ω22πGρc,     and,     α(34γg)

See also …


As we have detailed separately, the boundary condition at the center of a polytropic configuration is,

dxdξ|ξ=0=0;

and the boundary condition at the surface of an isolated polytropic configuration is,

dlnxdlnξ

=

α+ω2γg(14πGρc)ξ(θ')         at         ξ=ξs.

But this surface condition is not applicable to bipolytropes. Instead, let's return to the original, more general expression of the surface boundary condition:

dlnxdlnξ|s

=

α+ω2R3γgGMtot.


Utilizing an accompanying discussion, let's examine the frequency normalization used by Murphy & Fiedler (1985b) (see the top of the left-hand column on p. 223):

Ω2

ω2[R3GMtot]

 

=

ω2[34πGρ¯]=ω2[34πGρc]ρcρ¯=3ω2(nc+1)[(nc+1)4πGρc]ρcρ¯

 

=

3ω2(nc+1)[an2ρcPcθc]ρcρ¯=3γ(nc+1)ρcρ¯[an2ρcPcω2θcγ].

For a given radial quantum number, k, the factor inside the square brackets in this last expression is what Murphy & Fiedler (1985b) refer to as ωk2θc. Keep in mind, as well, that, in the notation we are using,

σc2

3ω22πGρc

σc2

=

(2ρ¯ρc)Ω2=6γ(nc+1)[an2ρcPcω2θcγ]=6γ(nc+1)[ωk2θc].

This also means that the surface boundary condition may be rewritten as,

dlnxdlnξ|s

=

Ω2γgα.


Let's apply these relations to the core and envelope, separately.


Envelope Layers With n = 5

The LAWE for n = 5 structures is, then,

0

=

d2xdη2+[46Q5]1ηdxdη+6[(σc26γenv)η2ϕαenvQ5]xη2

where,

Q5

dlnϕdlnη.

From our accompanying discussion of the underlying equilibrium structure of (nc,ne)=(1,5) bipolytropes, we know that,

ϕ

=

B01sinΔη1/2(32sin2Δ)1/2,

and,

dϕdη

=

B01[3cosΔ3sinΔ+2sin3Δ]2η3/2(32sin2Δ)3/2.

where A0 is a "homology factor," B0 is an overall scaling coefficient, and we have introduced the notation,

Δln(A0η)1/2=12(lnA0+lnη).

Hence,

Q5

=

η[η1/2(32sin2Δ)1/2B01sinΔ]B01[3cosΔ3sinΔ+2sin3Δ]2η3/2(32sin2Δ)3/2

 

=

3sinΔ3cosΔ2sin3Δ2sinΔ(32sin2Δ).

And,

0

=

d2xdη2+[4+3(3cosΔ3sinΔ+2sin3Δ)sinΔ(32sin2Δ)]1ηdxdη+[(σc2γenv)B0η1/2(32sin2Δ)1/2sinΔ+3αenv(3cosΔ3sinΔ+2sin3Δ)η2sinΔ(32sin2Δ)]x

 

=

d2xdη2+[4+3(3cosΔ32sinΔ12sin3Δ)sinΔ(2+cos2Δ)]1ηdxdη+[ωk2θc(γgγenv)B0η1/2(2+cos2Δ)1/2sinΔ+3αenv(3cosΔ32sinΔ12sin3Δ)η2sinΔ(2+cos2Δ)]x,

which matches the expression presented by Murphy & Fiedler (1985b) (see middle of the left column on p. 223 of their article) if we set θc=1 and γg/γenv=1.

Surface Boundary Condition

Next, pulling from our accompanying discussion of the stability of polytropes and an accompanying table that details the properties of (nc,ne)=(1,5) bipolytropes, the surface boundary condition is,

dlnxdlnη|s

=

(γgγenv)α+ω2R3γenvGMtot

dlnxdlnη|s+(γgγenv)α

=

ω2(Rs*)3γenvGMtot*(KcG)3/2(KcG)3/21ρ0

 

=

ω2γenvGρ0[(2π)1/2ξie2(πΔi)]3[(32π)1/2sinξi(3sin2Δi2)1/2e(πΔi)]1(μeμc)

 

=

ω2γenv(2πGρ0)(μeμc)13[ξi2θi](3sin2Δi2)1/2e5(πΔi)

 

=

ω2γenv(2πGρ0)(μeμc)e5π3[ξi2θi]ξi1/2Bθi(ξiA)5/2

 

=

ω2γenv(2πGρ0)(μeμc)Be5π3A5/2

 

=

2ωk2θc(nc+1)(μeμc)Be5π3A5/2.

After acknowledging that, in their specific stability analysis, θc=1, nc=1, and μe/μc=1, this right-hand-side expression matches the equivalent term published by Murphy & Fiedler (1985b) (see the bottom of the left-hand column on p. 223).

Core Layers With n = 1

And for n = 1 structures the LAWE is,

0

=

d2xdξ2+[42Q1]1ξdxdξ+2[(σc26γcore)ξ2θαcoreQ1]xξ2

where,

Q1

dlnθdlnξ.

Given that, for n=1 polytropic structures,

θ(ξ)=sinξξ       and       dθdξ=[cosξξsinξξ2]

we have,

Q1

=

ξ2sinξ[cosξξsinξξ2]

 

=

1ξcotξ.

Hence, the governing LAWE for the core is,

0

=

d2xdξ2+[42(1ξcotξ)]1ξdxdξ+2[(σc26γcore)ξ3sinξαcore(1ξcotξ)]xξ2

 

=

d2xdξ2+[1+ξcotξ]2ξdxdξ+2[(σc26γcore)ξ3sinξαcore(1ξcotξ)]xξ2.

This can be rewritten as,

0

=

d2xdξ2+2ξ[1+ξcotξ]dxdξ+[(σc23γcore)ξsinξ+2αcore(ξcosξsinξ)ξ2sinξ]x

 

=

d2xdξ2+2ξ[1+ξcotξ]dxdξ+[γgγcore(ωk2θc)ξsinξ+2αcore(ξcosξsinξ)ξ2sinξ]x,

which matches the expression presented by Murphy & Fiedler (1985b) (see middle of the left column on p. 223 of their article) if we set θc=1 and γg/γcore=1. This LAWE also appears in our separate discussion of radial oscillations in n = 1 polytropic spheres.

Interface Conditions

Here, we will simply copy the discussion already provided in the context of our attempt to analyze the stability of (nc,ne)=(0,0) bipolytropes; specifically, we will draw from STEP 4: in the Piecing Together subsection. Following the discussion in §§57 & 58 of P. Ledoux & Th. Walraven (1958), the proper treatment is to ensure that fractional perturbation in the gas pressure (see their equation 57.31),

δPP

=

γx(3+dlnxdlnξ),

is continuous across the interface. That is to say, at the interface (ξ=ξi), we need to enforce the relation,

0

=

[γcxcore(3+dlnxcoredlnξ)γexenv(3+dlnxenvdlnξ)]ξ=ξi

 

=

γe[γcγe(3+dlnxcoredlnξ)(3+dlnxenvdlnξ)]ξ=ξi

dlnxenvdlnξ|ξ=ξi

=

3(γcγe1)+γcγe(dlnxcoredlnξ)ξ=ξi.

In the context of this interface-matching constraint (see their equation 62.1), P. Ledoux & Th. Walraven (1958) state the following:   In the static (i.e., unperturbed equilibrium) modeldiscontinuities in ρ or in γ might occur at some [radius]. In the first case — that is, a discontinuity only in density, while γe=γc — the interface conditions imply the continuity of 1xdxdξ at that [radius]. In the second case — that is, a discontinuity in the adiabatic exponent — the dynamical condition may be written as above. This implies a discontinuity of the first derivative at any discontinuity of γ.

The algorithm that Murphy & Fiedler (1985b) used to "… [integrate] through each zone …" was designed "… with continuity in x and dx/dξ being imposed at the interface …" Given that they set γc=γe=5/3, their interface matching condition is consistent with the one prescribed by P. Ledoux & Th. Walraven (1958).

Our Numerical Integration

Let's try to integrate this bipolytrope's LAWE from the center, outward, using as a guideline an accompanying Numerical Integration outline. Generally, for any polytropic index, the relevant LAWE can be written in the form,

θixi

=

[𝒜]xiξi(n+1)6[]xi

where,

𝒜

4θi(n+1)ξi(θ')i=θi[4(n+1)Qi]

σc2γg2α(3θ'ξ)i=𝔉+2α[1(3θ'ξ)i]=𝔉+2α[13θiξi2Qi]

𝔉

[σc2γg2α]=[σc2γg2(34γg)]=[(8+σc2)γg6]

 

σc2=γg(𝔉+6)8.

This leads to a discrete, finite-difference representation of the form,

x+[2θi+δξξi𝒜]

=

x[δξξi𝒜2θi]+xi{4θi(δξ)2(n+1)3}.

This provides an approximate expression for x+xi+1, given the values of xxi1 and xi; this works for all zones, i=3N as long as the center of the configuration is denoted by the grid index, i=1. Note that,

δξ

ξmax(N1)

        and        

ξi

=

(i1)δξ.

In order to kick-start the integration, we will set the displacement function value to x1=1 at the center of the configuration (ξ1=0), then we will draw on the derived power-series expression to determine the value of the displacement function at the first radial grid line, ξ2=δξ, away from the center. Specifically, we will set,

x2

=

x1[1(n+1)𝔉(δξ)260].

Integration Through the n = 1 Core

For an n=1 core, we have,

θi=sinξiξi       and       Qi=1ξicotξi.

Hence,

𝒜core

=

sinξiξi[42(1ξicotξi)]=2sinξiξi[1+ξicotξi]

core

=

𝔉core+2αcore[13θiξi2Qi]=𝔉core+2αcore[13sinξiξi3(1ξicotξi)].

So, first we choose a value of σc2 and γc, which means,

𝔉core

[(8+σc2)γc6]

Then, moving from the center of the configuration, outward to the interface at ξi=ξinterfaceδξ=ξinterface/(N1), we have,

x1

=

1,

x2

=

x1[1𝔉core(δξ)230],

for i=2N,        xi+1[2θi+δξξi𝒜core]

=

xi1[δξξi𝒜core2θi]+xi{4θi(δξ)2(n+1)3core}.

At the interface — that is, when i=N — the logarithmic slope of the displacement function is,

dlnxdlnξ|interface

ξNxN(xN+1xN1)2δξ.

Interface

Keep in mind that, as has been detailed in the accompanying equilibrium structure chapter, for (nc,ne)=(1,5) bipolytropes,

r*

=

(12π)1/2ξ,

          for,

0ξξinterface.

r*

=

[(μeμc)1(32π)1/2]η,

          for,

ηinterfaceηηs;

ηs

=

13(μeμc)ξs=13(μeμc)[ξe2(πΔ)]interface.

         

We now need to determine what the slope is at the interface, viewed from the perspective of the envelope. From above, we deduce that,

dlnydlnη|interface

=

3(γcγe1)+γcγedlnxdlnξ|interface.

Hence, letting the subscript "1" denote the interface location as viewed from the envelope, we have,

η1y1(y2y0)(η2η0)

=

3(γcγe1)+γcγedlnxdlnξ|interface.

y0

=

y22(δη)y1η1{3(γcγe1)+γcγedlnxdlnξ|interface}.

Integration Through the n = 5 Envelope

For an n=5 envelope, we have,

ϕi=B01sinΔη1/2(32sin2Δ)1/2       and       Qi=dlnϕdlnη=3sinΔ3cosΔ2sin3Δ2sinΔ(32sin2Δ),

where A0 is a "homology factor," B0 is an overall scaling coefficient, and we have introduced the notation,

Δln(A0η)1/2=12(lnA0+lnη).

Hence,

𝒜env

ϕi[4(n+1)Qi]

 

=

B01sinΔη1/2(32sin2Δ)1/2{46[3sinΔ3cosΔ2sin3Δ2sinΔ(32sin2Δ)]}

env

σc2γe2αenv(3ϕ'η)i

 

=

1γe[γc(𝔉core+6)8]2[34γe]Qi(3ϕiηi2)

 

=

γcγe[𝔉core+68γc]6[34γe]B01sinΔη5/2(32sin2Δ)1/2[3sinΔ3cosΔ2sin3Δ2sinΔ(32sin2Δ)]

 

=

γcγe[𝔉core+68γc]3B01[34γe][3sinΔ3cosΔ2sin3Δη5/2(32sin2Δ)3/2].

This leads to a discrete, finite-difference representation of the form,

y+[2ϕi+δηηi𝒜env]

=

y[δηηi𝒜env2ϕi]+yi[4ϕi2(δη)2env]

This provides an approximate expression for y+yi+1, given the values of yyi1 and yi; this works for all zones, i=3M as long as the interface between the core and the envelope of the configuration is denoted by the grid index, i=1. Note that,

δη

ηsurfηinterfaceM1

        and        

ηi

=

ηinterface+(i1)δη.

At the interface, we need special treatment in order to ensure that both the amplitude and the first derivative of the displacement function behave properly. Specifically, when i=1, we must set, y1=xN and η1=(μe/μc)ξN/3. Then the value of y2 is obtained from the expression,

y2[2ϕ1+δηη1𝒜env]

=

y0[δηη1𝒜env2ϕ1]+y1{4ϕ12(δη)2env}

 

=

y1[4ϕ12(δη)2env]+[δηη1𝒜env2ϕ1]{y22(δη)y1η1[3(γcγe1)+γcγedlnxdlnξ|interface]}

y2[4ϕ1]

=

y1[4ϕ12(δη)2env]2(δη)y1η1[δηη1𝒜env2ϕ1]{3(γcγe1)+γcγedlnxdlnξ|interface}

y2y1[ϕ1]

=

[ϕ1(δη)22env](δη)2η1[δηη1𝒜env2ϕ1]{3(γcγe1)+γcγedlnxdlnξ|interface}.

Regroup

Foundation

In an accompanying discussion, we derived the so-called,

Adiabatic Wave (or Radial Pulsation) Equation

User:Tohline/Math/EQ RadialPulsation01


whose solution gives eigenfunctions that describe various radial modes of oscillation in spherically symmetric, self-gravitating fluid configurations. Assuming that the underlying equilibrium structure is that of a bipolytrope having (nc,ne)=(1,5), it makes sense to adopt the normalizations used when defining the equilibrium structure, namely,

ρ*

ρ0ρc

;    

r*

r0(Kc/G)1/2

P*

P0Kcρc2

;    

Mr*

M(r0)ρc(Kc/G)3/2

H*

HKcρc

.    

 

We note as well that,

g0

=

GM(r0)r02

 

=

G[Mr*ρc(KcG)3/2][r*(KcG)1/2]2

 

=

Mr*(r*)2[Gρc(KcG)1/2].

Hence, multiplying the LAWE through by (Kc/G) gives,

0

=

d2xdr*2+[4r*(KcG)1/2(g0ρ0P0)]dxdr*+(KcG)(ρ0γgP0)[ω2+(43γg)g0r0]x

 

=

d2xdr*2+{4r*(KcG)1/2(ρcρ*P*Kcρc2)Mr*(r*)2[Gρc(KcG)1/2]}dxdr*+(KcG)(ρ*ρcγgP*Kcρc2){ω2+(43γg)1r*(GKc)1/2Mr*(r*)2[Gρc(KcG)1/2]}x

 

=

d2xdr*2+{4r*(ρ*P*)Mr*(r*)2}dxdr*+(1γgGρc)(ρ*P*){ω2+(43γg)1r*Mr*(r*)2[Gρc]}x

 

=

d2xdr*2+{4r*(ρ*P*)Mr*(r*)2}dxdr*+(ρ*P*){ω2γgGρc+(4γg3)1r*Mr*(r*)2}x

 

=

d2xdr*2+{4(ρ*P*)Mr*(r*)}1r*dxdr*+(ρ*P*){2πσc23γgαgMr*(r*)3}x.

Profile

Now, referencing the derived bipolytropic model profile, we should incorporate the following relations:

Variable

Throughout the Core
0r*ξi2π

Throughout the Envelope
ξi2πr*ξie2(πΔi)2π

Plotted Profiles

ξi=0.5

ξi=1.0

ξi=3.0

 

ξ=2πr*

η=(μeμc)(2π3)1/2r*

 

ρ*

sinξξ

(μeμc)θi[ϕ(η)]5

P*

(sinξξ)2

θi2[ϕ(η)]6

Mr*

(2π)1/2(sinξξcosξ)

(μeμc)2(233π)1/2θi(η2dϕdη)

In order to obtain the various envelope profiles, it is necessary to evaluate ϕ(η) and its first derivative using the information presented in Step 6 of our accompanying discussion.


Therefore, throughout the core we have,

ρ*P*

=

ξsinξ;

Mr*r*

=

2πξ(2π)1/2(sinξξcosξ)=2sinξξ(1ξcotξ).

In which case the governing LAWE throughout the core is,

0

=

d2xdr*2+{42(1ξcotξ)}1r*dxdr*+ξsinξ{2πσc23γgαg4πsinξξ3(1ξcotξ)}x

 

=

d2xdr*2+{42(1ξcotξ)}1r*dxdr*+2πξsinξ{σc23γg+2αgξ3(ξcosξsinξ)}x.

Next, throughout the envelope we have,

ρ*P*

=

(μeμc)θi[ϕ(η)]5{θi2[ϕ(η)]6}1=(μeμc)1θiϕ(η);

Mr*r*

=

(μeμc)2(233π)1/2θi(η2dϕdη){1η(μeμc)(2π3)1/2}=6(μeμc)1θi(ηdϕdη)

So, the governing LAWE throughout the envelope is,

0

=

d2xdr*2+{46(μeμc)1θi(ηdϕdη)(μeμc)1θiϕ(η)}1r*dxdr*+(μeμc)1θiϕ(η){2πσc23γg6αg(μeμc)1θi(ηdϕdη)[1η(μeμc)(2π3)1/2]2}x

 

=

d2xdr*2+[46(dlnϕdlnη)]1r*dxdr*+(μeμc)1θiϕ(η){2πσc23γg6αg(μeμc)θiη(dϕdη)(2π3)}x

 

=

d2xdr*2+[46(dlnϕdlnη)]1r*dxdr*+2π3(μeμc)21η2{(μeμc)1(σc2γg)η2θiϕ(η)6αg(dlnϕdlnη)}x.

Model 10

As we have reviewed in an accompanying discussion, equilibrium Model 10 from Murphy & Fiedler (1985, Proc. Astr. Soc. of Australia, 6, 219) is defined by setting (ξi,m)=(2.5646,1). Drawing directly from our reproduction of their Table 1, we see that a few relevant structural parameters of Model 10 are,

ξs

=

6.5252876

riR=ξiξs

=

0.39302482

ρcρ¯

=

34.346

MenvMtot

=

5.89×104

Here we list a few other model parameter values that will aid in our attempt to correctly integrate the LAWE to find various radial oscillation eigenvectors.

A Sampling of Model 10's Equilibrium Parameter Values

Grid
Line
rR ξ η Δ ϕ dϕdη r* ρ* P* Mr* g0*Mr*(r*)2
25 0.12093071 0.789108         0.31480842 0.89940188 0.80892374 0.122726799 1.23835945
40 0.19651241 1.2823         0.51156369 0.74761972 0.55893525 0.473819194 1.81056130
79 0.393025 2.5646         1.02312737 0.21270605 0.04524386 2.150231108 2.05411964
79 0.393025   1.4806725 2.6746514 1.000000 1.112155 1.02312737 0.21270605 0.04524386 2.15023111 2.0541196
100 0.49883919   1.8793151 2.7938569 0.6505914 0.69070815 1.2985847 0.0247926 0.0034309 2.15127319 1.2757189
150 0.7507782   2.8284641 2.9982701 0.2149684 0.30495637 1.95443562 9.7646E-05 4.4649E-06 2.15149752 0.563246
199 0.9976784   3.7586302 3.1404305 0.00150695 0.17269514 2.59716948 1.653E-15 5.2984E-19 2.15149876 0.31896316

Our chosen (uniform) grid spacing is,

δrR=178(riR)0.00503878;

as a result, the center is at zone 1, the interface is at grid line 79, and the surface is just beyond grid line 199.

Numerical Integration

General Approach

Here, we begin by recognizing that the 2nd-order ODE that must be integrated to obtain the desired eigenvectors has the generic form,

0

=

x+r*x+𝒦x,

where,

x

=

dxdr*

      and      

x

=

d2xd(r*)2.

Adopting the same approach as before when we integrated the LAWE for pressure-truncated polytropes, we will enlist the finite-difference approximations,

x

x+x2δr*

      and      

x

x+2xj+x(δr*)2.

The finite-difference representation of the LAWE is, therefore,

x+2xj+x(δr*)2

=

r*[x+x2δr*]𝒦xj

x+2xj+x

=

δr*2r*[x+x](δr*)2𝒦xj

xj+1[1+(δr*2r*)]

=

[2(δr*)2𝒦]xj[1(δr*2r*)]xj1.

In what follows we will also find it useful to rewrite 𝒦 in the form,

𝒦(σc2γg)𝒦1αg𝒦2.

Case A:   From the above Foundation discussion, the relevant coefficient expressions for all regions of the configuration are,

{4(ρ*P*)Mr*(r*)}

      ,      

𝒦1

2π3(ρ*P*)

      and      

𝒦2

(ρ*P*)Mr*(r*)3.


Case B:   Alternatively, immediately following the above Profile discussion, the relevant coefficient expressions for the core are,

{42(1ξcotξ)}

      ,      

𝒦1

2π3(ξsinξ)

      and      

𝒦2

4πξ2sinξ(sinξξcosξ);

while the coefficient expressions for the envelope are,

=

{46(dlnϕdlnη)}

      ,      

𝒦1

=

2π3(μeμc){1θiϕ(η)}

      and      

𝒦2

=

12π3(μeμc)21η2(dlnϕdlnη).

Grid
Line
rR ξ η Case A Case B
𝒦1 𝒦2 𝒦1 𝒦2
25 0.12093071 0.789108   3.566549 2.328653 4.373676 3.566549 2.328653 4.373676
40 0.19651241 1.2823   2.761112 2.801418 4.734049 2.761112 2.801418 4.734049
79 0.393025 2.5646   -5.880425 9.846430 9.4387879 -5.880424 9.846430 9.438787
79 0.393025   1.4806725 -5.880425 9.846430 9.4387879 -5.880424 9.846430 9.438787
100 0.49883919   1.8793151 -7.971244 15.134659 7.099025 -7.971184 15.134583 7.098989
150 0.7507782   2.8284641 -2.00748E+01 4.58038E+01 6.30260 -2.00749E+01 4.58041E+01 6.30264
199 0.9976784   3.7586302 -2.58045E+03 6.53411E+03 3.83150E+02 -2.58041E+03 6.53401E+03 3.83144E+02

Special Handling at the Center

In order to kick-start the integration, we set the displacement function value to x1=1 at the center of the configuration (ξ1=0), then draw on the derived power-series expression to determine the value of the displacement function at the first radial grid line, ξ2=δξ, away from the center. Specifically, we set,

x2

=

x1[1(n+1)𝔉(δξ)260].

Special Handling at the Interface

Integrating outward from the center, the general approach will work up through the determination of xj+1 when "j+1" refers to the interface location. In order to properly transition from the core to the envelope, we need to determine the value of the slope at this interface location. Let's do this by setting j = i, then projecting forward to what x+ would be — that is, to what the amplitude just beyond the interface would be — if the core were to be extended one more zone. Then, the slope at the interface (as viewed from the perspective of the core) will be,

x'i|core

12δr*{x+xi1}

 

=

xi12δr*+12δr*{[2(δr*)2𝒦]xi[1(δr*2r*)]xi1}[1+(δr*2r*)]1

 

=

12δr*{[2(δr*)2𝒦]xi[1(δr*2r*)]xi1[1+(δr*2r*)]xi1}[1+(δr*2r*)]1

 

=

12δr*{[2(δr*)2𝒦]xi2xi1}[1+(δr*2r*)]1

Conversely, as viewed from the envelope, if we assume that we know xi and x'i, we can determine the amplitude, xi+1, at the first zone beyond the interface as follows:

x

xi+12δr*x'i|env

xi+1[1+(δr*2r*)]

=

[2(δr*)2𝒦]xi[1(δr*2r*)][xi+12δr*x'i|env]

xi+1[1+(δr*2r*)]+[1(δr*2r*)]xi+1

=

[2(δr*)2𝒦]xi+[1(δr*2r*)]2δr*x'i|env

xi+1

=

[112(δr*)2𝒦]xi+[1(δr*2r*)]δr*x'i|env

Eigenvectors

Keep in mind that, for all models, we expect that, at the surface, the logarithmic derivative of each proper eigenfunction will be,

dlnxdlnr*|surf

=

Ω2γα.

Also, keep in mind that, for Model 10 (ξi=2.5646):

riR

=

0.39302482

      ,    

ρcρ¯

=

34.3460405


Our Determinations for Model 10

Mode σc2 Ω2σc22(ρcρ¯) xsurf dlnxdlnr*|surf rR|1 1MrMtot|1 rR|2 1MrMtot|2 rR|3 1MrMtot|3
expected measured
1
(Fundamental)
0.92813095170326 15.93881161 +85.17 8.963286966 8.963085 n/a n/a n/a n/a n/a n/a
2 1.237156768978 21.24571822 - 610 12.14743093 12.147337 0.5724 3.05E-05 n/a n/a n/a n/a
3 1.8656033984 32.0380449 +3225 18.62282676 18.6228 0.4845 1.35E-04 0.787 2.05E-07 n/a n/a
4 2.65901504799 45.66331921 -9410 26.79799153 26.797977 0.4459 2.620E-04 0.7096 1.834E-06 0.8632 1.189E-08
Match Figure 2 from MF85


For Model 17 (ξi=3.0713):

riR

=

0.93276717

      ,    

ρcρ¯

=

3.79693903


Our Determinations for Model 17

Mode σc2 Ω2σc22(ρcρ¯) xsurf dlnxdlnr*|surf rR|1 1MrMtot|1 rR|2 1MrMtot|2 rR|3 1MrMtot|3
expected measured
1
(Fundamental)
1.149837904 2.182932207 +1.275 0.7097593 0.7097550 n/a n/a n/a n/a n/a n/a
2 7.34212930615 13.93880866 - 2.491 7.763285 7.763244 0.7215 0.24006 n/a n/a n/a n/a
3 16.345072567 31.03062198 +4.33 18.01837 18.01826 0.5806 0.5027 0.848 0.0541 n/a n/a
4 27.746934203 52.6767087 -9.1 31.0060 31.0058 0.4859 0.6737 0.7429 0.1974 0.8957 0.0171
Match Figure 3 from MF85


Numerical Values for Some Selected (nc,ne)=(1,5) Bipolytropes
[to be compared with Table 1 of Murphy & Fiedler (1985)]

MODEL Source riR Ω02 Ω12 rR|1 1MrMtot|1
10 MF85 0.393 15.9298 21.2310 0.573 1.00E-03
Here 0.39302 15.93881161 21.24571822 0.5724 3.05E-05
17 MF85 0.933 2.1827 13.9351 0.722 0.232
Here 0.93277 2.182932207 13.93880866 0.7215 0.24006

Try Splitting Analysis Into Separate Core and Envelope Components

Core:

Given that, 2πr*=ξ, lets multiply the LAWE through by (2π)1. This gives,

0

=

d2xdξ2+{4(ρ*P*)Mr*(r*)}1ξdxdξ+12π(ρ*P*){2πσc23γgαgMr*(r*)3}x.

Specifically for the core, therefore, the finite-difference representation of the LAWE is,

x+2xj+x(δξ)2

=

ξ[x+x2δξ][𝒦2π]xj

x+2xj+x

=

δξ2ξ[x+x](δξ)2[𝒦2π]xj

xj+1[1+(δξ2ξ)]

=

[2(δξ)2(𝒦2π)]xj[1(δξ2ξ)]xj1.

This also means that, as viewed from the perspective of the core, the slope at the interface is

[dxdξ]interface

=

12δξ{[2(δξ)2(𝒦2π)]xi2xi1}[1+(δξ2ξ)]1.

Envelope:

Given that,

(μeμc)(2π3)1/2r*=η,

let's multiply the LAWE through by (3/2π)(μe/μc)2. This gives,

0

=

d2xdη2+{4(ρ*P*)Mr*(r*)}1ηdxdη+32π(μeμc)2(ρ*P*){2πσc23γgαgMr*(r*)3}x.

Specifically for the envelope, therefore, the finite-difference representation of the LAWE is,

x+2xj+x(δη)2

=

η[x+x2δη](μeμc)2[3𝒦2π]xj

x+2xj+x

=

δη2η[x+x](δη)2(μeμc)2[3𝒦2π]xj

xj+1[1+(δη2η)]

=

[2(δη)2(μeμc)2(3𝒦2π)]xj[1(δη2η)]xj1.

This also means that, once we know the slope at the interface (see immediately below), the amplitude at the first zone outside of the interface will be given by the expression,

xi+1

=

[112(δη)2(μeμc)2(3𝒦2π)]xi+[1(δη2η)]δη[dxdη]interface.

Interface

If we consider only cases where γe=γc, then at the interface we expect,

dlnxdlnr*

=

dlnxdlnξ=dlnxdlnη

r*dxdr*

=

ξdxdξ=ηdxdη

dxdr*

=

(2π)1/2dxdξ=(μeμc)(2π3)1/2dxdη.

Switching at the interface from ξ to η therefore means that,

[dxdη]interface

=

3(μeμc)1[dxdξ]interface.

Begin Our Analysis

Relevant LAWEs

The LAWE that is relevant to polytropic spheres may be written as,

User:Tohline/Math/EQ RadialPulsation02

Core Layers With n = 5

The LAWE for n = 5 structures is,

0

=

d2xdξ2+[46Q5]1ξdxdξ+6[(σc26γcore)ξ2θ5αcoreQ5]xξ2

where,

Q5

dlnθ5dlnξ.

From our study of the equilibrium structure of (nc,ne)=(5,1) bipolytropes, we have,

θ5

=

[1+13ξ2]1/2;

dθ5dξ

=

ξ3[1+13ξ2]3/2.

Q5=ξθ5dθ5dξ

=

[1+13ξ2]1/2ξ23[1+13ξ2]3/2

 

=

ξ23[1+13ξ2]1.

Hence, for the core the governing LAWE is,

0

=

d2xdξ2+{42ξ2[1+13ξ2]1}1ξdxdξ+6{(σc26γcore)ξ2[1+13ξ2]1/2αcoreξ23[1+13ξ2]1}xξ2

 

=

d2xdξ2+{2ξ2[1+13ξ2]1}2ξdxdξ+{(σc2γcore)[1+13ξ2]1/22αcore[1+13ξ2]1}x

 

=

d2xdξ2+{23ξ2[3+ξ2]1}2ξdxdξ+{(σc23γcore)[3+ξ2]1/26αcore[3+ξ2]1}x

0

=

(3+ξ2)d2xdξ2+{2(3+ξ2)3ξ2}2ξdxdξ+{(σc23γcore)(3+ξ2)3/26αcore}x

 

=

(3+ξ2)d2xdξ2+(6ξ2)2ξdxdξ+[(σc23γcore)(3+ξ2)3/26αcore]x.

This exactly matches our derivation performed in the context of pressure-truncated polytropes. When we insert the eigenfunction obtained via a Eureka Moment on 3/6/2017,

x

=

x0[1ξ215],

dxdξ

=

2x0ξ15

d2xdξ2

=

2x015,

we obtain,

(3+ξ2)d2xdξ2+(6ξ2)2ξdxdξ+[(σc23γcore)(3+ξ2)3/26αcore]x

=

(3+ξ2)2x0152(6ξ2)2x015+x015[(σc23γcore)(3+ξ2)3/26αcore](15ξ2)

 

=

x015{2(3+ξ2)+4(6ξ2)[(σc23γcore)(3+ξ2)3/26αcore](15ξ2)}

 

=

x015{302ξ2+6αcore(15ξ2)[(σc23γcore)(3+ξ2)3/2](15ξ2)}.

The right-hand-side goes to zero if αcore=1/3 and σc=0. Also, notice that,

dlnxdlnξ|core=ξxdxdξ|core

=

2x0ξ215{x0[1ξ215]}1

 

=

{152ξ2[15ξ215]}1

 

=

[15ξ22ξ2]1

 

=

2[115ξ2]1.

So, with γc=6/5 and γe=2 we need,

dlnxenvdlnξ|ξ=ξi

=

3(γcγe1)+γcγe(dlnxcoredlnξ)ξ=ξi

 

=

3(351)+65[115ξi2]1

 

=

65{[ξi215ξi2]11}

 

=

65[15ξi215].

Envelope Layers With n = 1

And for n = 1 structures the LAWE is,

0

=

d2xdη2+[42Q1]1ηdxdη+2[(σc26γenv)η2θ1αenvQ1]xη2

where,

Q1

dlnθ1dlnη.

As has already been pointed out, above, for n = 1 polytropic spheres, this LAWE becomes,

0

=

d2xdη2+2η[1+ηcotη]dxdη+[γgγenv(ωk2θc)ηsinη+2αenv(ηcosηsinη)η2sinη]x.

In a separate chapter, we explain that the analytically defined eigenfunction that satisfies this LAWE — when ωk2=0 and αenv=+1 — is,

xP|n=1

=

3beη2[1ηcotη].

Summary

For a given choice of the equilibrium model parameters, ξi and μe/μc, we can pull the parameters and profiles of the base equilibrium model from our accompanying chapter on (nc,ne)=(5,1) bipolytropes. Note, in particular, that:

r*

=

(32π)1/2ξ,

          for,

0ξξi.

r*

=

[(μeμc)112π(1+ξi23)]η,

          for,

ηiηηs;

ηs

=

ηi+π2+tan1[1ηiξi3].

         


Then, for a choice of the pair of exponents, γcore and γenv, that govern the behavior of adiabatic oscillations, the LAWE to be numerically integrated is:

0

=

(3+ξ2)d2xcoredξ2+(6ξ2)2ξdxcoredξ+[(σc23γcore)(3+ξ2)3/26αcore]xcore,

          for,

0ξξi;

0

=

d2xenvdη2+[1+ηcotη]2ηdxenvdη+2[(σc26γenv)η3sinηαenv(1ηcotη)]xenvη2,

          for,

ηiηηs.

The boundary conditions at the center of the configuration are, xenv=1, while dxenv/dξ=0.

As described above, the two interface boundary conditions are, xenv=xcore, and,

dlnxenvdlnξ|ξ=ξi

=

3(γcoreγenv1)+γcoreγenv(dlnxcoredlnξ)ξ=ξi.

Finally, drawing from a related discussion of the surface boundary condition for isolated n = 3 polytropes, at the surface we need,

dlnxenvdlnη|surface

=

1γenv(43γenv+ω2R3GMtot)

 

=

[3ω2R34πGγenvρ¯αenv]

 

=

12[σc2γenv(ρcρ¯)2αenv],

where, for an (nc,ne)=(5,1) bipolytrope,

ρcρ¯

=

(μeμc)1ηs23Aθi5.

See Also


Tiled Menu

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