Editing SSC/Stability/BiPolytropes/Pt3

=Radial Oscillations of (n<sub>c</sub>, n<sub>e</sub>) = (5, 1) Models=
==Foundation==
In an [[SSC/Perturbations#2ndOrderODE|accompanying discussion]], we derived the so-called,

<div align="center" id="2ndOrderODE">
<font color="#770000">'''Adiabatic Wave''' (or ''Radial Pulsation'') '''Equation'''</font><br />

{{Math/EQ_RadialPulsation01}}
</div>

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 <math>~(n_c, n_e) = (5, 1)</math>, it makes sense to adopt the normalizations used when defining the equilibrium structure, namely,
<div align="center">
<table border="0" cellpadding="3">
<tr>
  <td align="right">
<math>~\rho^*</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~\frac{\rho_0}{\rho_c}</math>
  </td>

  <td align="center">; &nbsp;&nbsp;&nbsp;</td>

  <td align="right">
<math>~r^*</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~\frac{r_0}{[K_c^{1/2}/(G^{1/2}\rho_c^{2/5})]}</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~P^*</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~\frac{P_0}{K_c\rho_c^{6/5}}</math>
  </td>

  <td align="center">; &nbsp;&nbsp;&nbsp;</td>

  <td align="right">
<math>~M_r^*</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~\frac{M_r}{[K_c^{3/2}/(G^{3/2}\rho_c^{1/5})]}</math>
  </td>
</tr>

</table>
</div>

We note as well that,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~g_0</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\frac{GM(r_0)}{r_0^2}</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
G \biggl[ M_r^* \rho_c^{-1 / 5} \biggl( \frac{K_c}{G}\biggr)^{3 / 2} \biggr] \biggl[ r^* \rho_c^{-2 / 5}\biggl( \frac{K_c}{G}\biggr)^{1 / 2} \biggr]^{-2}
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{G M_r^*}{(r^*)^2} \biggl[ \rho_c^{3 / 5} \biggl( \frac{K_c}{G}\biggr)^{1 / 2} \biggr] \, .
</math>
  </td>
</tr>
</table>

Hence, multiplying the LAWE through by <math>~(K_c/G)\rho_c^{-4 / 5}</math> gives,

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

<tr>
  <td align="right">
<math>~0</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{d^2x}{dr_0^2} + \biggl[\frac{4}{r_0} - \biggl(\frac{g_0 \rho_0}{P_0}\biggr) \biggr] \frac{dx}{dr_0} 
+ \biggl(\frac{\rho_0}{\gamma_\mathrm{g} P_0} \biggr)\biggl[\omega^2 + (4 - 3\gamma_\mathrm{g})\frac{g_0}{r_0} \biggr]  x 
</math>
  </td>
</tr>
</table>



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

<tr>
  <td align="right">
<math>~0</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{d^2x}{dr*^2} + \biggl[\frac{4}{r^*} ~-~ \rho_c^{-2 / 5} \biggl( \frac{K_c}{G} \biggr)^{1 / 2}\biggl(\frac{g_0 \rho_0}{P_0}\biggr) \biggr] \frac{dx}{dr*} 
~+~ \rho_c^{-4 / 5}\biggl( \frac{K_c}{G} \biggr)\biggl(\frac{\rho_0}{\gamma_\mathrm{g} P_0} \biggr)\biggl[\omega^2 ~+~ (4 - 3\gamma_\mathrm{g})\frac{g_0}{r_0} \biggr]  x
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{d^2x}{dr*^2} ~+~ \biggl\{\frac{4}{r^*} ~-~ \rho_c^{-2 / 5} \biggl( \frac{K_c}{G} \biggr)^{1 / 2}  \frac{G M_r^*}{(r^*)^2} \biggl[ \rho_c^{3 / 5} \biggl( \frac{K_c}{G}\biggr)^{1 / 2} \biggr] \biggl[ \frac{\rho_c \rho^*}{P^* K_c \rho_c^{6/5}}\biggr]  \biggr\} \frac{dx}{dr*} 
~+~ \rho_c^{-4 / 5}\biggl( \frac{K_c}{G} \biggr)\biggl[ \frac{\rho_c \rho^*}{\gamma_\mathrm{g}P^* K_c \rho_c^{6/5}}\biggr] 
 \biggl\{\omega^2 ~+~ (4 - 3\gamma_\mathrm{g}) \frac{G M_r^*}{(r^*)^2} \biggl[ \rho_c^{3 / 5} \biggl( \frac{K_c}{G}\biggr)^{1 / 2} \biggr] \frac{\rho_c^{2 / 5}}{r^*}\biggl( \frac{G}{K_c}\biggr)^{1 / 2} \biggr\}  x
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{d^2x}{dr*^2} ~+~ \biggl\{\frac{4}{r^*} ~-~ \frac{M_r^*}{(r^*)^2} \biggl[ \frac{\rho^*}{P^* }\biggr]  \biggr\} \frac{dx}{dr*} 
~+~ \biggl( \frac{1}{G\rho_c} \biggr)\biggl[ \frac{ \rho^*}{\gamma_\mathrm{g}P^* }\biggr] 
 \biggl\{\omega^2 ~+~ (4 - 3\gamma_\mathrm{g}) \frac{G\rho_c M_r^*}{(r^*)^3} \biggr\}  x
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{d^2x}{dr*^2} + \biggl\{ \frac{4}{r^*} -\biggl(\frac{\rho^*}{P^*}\biggr)\frac{ M_r^*}{(r^*)^2}\biggr\} \frac{dx}{dr*} 
+ \biggl(\frac{\rho^*}{ P^* } \biggr)\biggl\{ \frac{\omega^2}{\gamma_\mathrm{g} G\rho_c} + \biggl(\frac{4}{\gamma_\mathrm{g}} - 3\biggr)\frac{ M_r^*}{(r^*)^3}\biggr\}  x 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{d^2x}{dr*^2} + \biggl\{ 4 -\biggl(\frac{\rho^*}{P^*}\biggr)\frac{ M_r^*}{(r^*)}\biggr\}\frac{1}{r^*} \frac{dx}{dr*} 
+ \biggl(\frac{\rho^*}{ P^* } \biggr)\biggl\{ \frac{2\pi \sigma_c^2}{3\gamma_\mathrm{g}}  ~-~\frac{\alpha_\mathrm{g} M_r^*}{(r^*)^3}\biggr\}  x \, .
</math>
  </td>
</tr>
</table>

==Profile==

Now, referencing the [[SSC/Structure/BiPolytropes/Analytic51#Profile|derived bipolytropic model profile]], we should incorporate the following relations:


<div align="center">
<b>Table 1:  Radial Profile of Various Physical Variables</b>
<table border="1" cellpadding="6">
<tr>
  <td align="center" rowspan="2">
Variable
  </td>
  <td align="center" rowspan="2">
Throughout the Core<br>
<math>~0 \le \xi \le \xi_i</math>
  </td>
  <td align="center" rowspan="2">
Throughout the Envelope<sup>&dagger;</sup><br>
<math>~\eta_i \le \eta \le \eta_s</math>
  </td>

  <td align="center" colspan="3">
Plotted Profiles
  </td>
</tr>

<tr>
  <td align="center">
<math>~\xi_i = 0.5</math>
  </td>
  <td align="center">
<math>~\xi_i = 1.0</math>
  </td>
  <td align="center">
<math>~\xi_i = 3.0</math>
  </td>
</tr>

<tr>
  <td align="center">
<math>~r^*</math>
  </td>
  <td align="center">
<math>\biggl( \frac{3}{2\pi} \biggr)^{1/2} \xi</math>
  </td>
  <td align="center">
<math>\biggl( \frac{\mu_e}{\mu_c} \biggr)^{-1} \theta^{-2}_i (2\pi)^{-1/2}\eta</math>
  </td>

  <td align="center" colspan="3">
&nbsp;
  </td>
</tr>

<tr>
  <td align="center">
<math>~\rho^*</math>
  </td>
  <td align="center">
<math>\biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{-5/2}</math>
  </td>
  <td align="center">
<math>\biggl( \frac{\mu_e}{\mu_c} \biggr) \theta^{5}_i \phi(\eta)</math>
  </td>

  <td align="center">
<!--  [[File:PlotDensity_xi_0.5.jpg|thumb|75px]] -->
[[Image:DenXi05.jpg|thumb|75px]]
  </td>
  <td align="center">
[[Image:DenXi10.jpg|thumb|75px]]
  </td>
  <td align="center">
[[Image:DenXi30.jpg|thumb|75px]]
  </td>
</tr>

<tr>
  <td align="center">
<math>~P^*</math>
  </td>
  <td align="center">
<math>\biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{-3}</math>
  </td>
  <td align="center">
<math>\theta^{6}_i [\phi(\eta)]^{2}</math>
  </td>

  <td align="center">
<!-- [[File:PlotPressure_xi_0.5.jpg|thumb|75px]] -->
[[Image:PresXi05.jpg|thumb|75px]]
  </td>
  <td align="center">
[[Image:PresXi10.jpg|thumb|75px]]
  </td>
  <td align="center">
[[Image:PresXi30.jpg|thumb|75px]]
  </td>
</tr>

<tr>
  <td align="center">
<math>~M_r^*</math>
  </td>
  <td align="center">
<math>\biggl( \frac{2\cdot 3}{\pi } \biggr)^{1/2} \biggl[ \xi^3 \biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{-3/2} \biggr]</math>
  </td>
  <td align="center">
<math>\biggl( \frac{\mu_e}{\mu_c} \biggr)^{-2} \theta^{-1}_i \biggl( \frac{2}{\pi} \biggr)^{1/2} \biggl(-\eta^2 \frac{d\phi}{d\eta} \biggr)</math>
  </td>

  <td align="center">
<!-- [[File:PlotPressure_xi_0.5.jpg|thumb|75px]] -->
[[Image:MassXi05.jpg|thumb|75px]]
  </td>
  <td align="center">
[[Image:MassXi10.jpg|thumb|75px]]
  </td>
  <td align="center">
[[Image:MassXi30.jpg|thumb|75px]]
  </td>
</tr>

<tr>
  <td align="left" colspan="6">
<sup>&dagger;</sup>In order to obtain the various envelope profiles, it is necessary to evaluate <math>\phi(\eta)</math> and its first derivative using the information presented in Step 6, above.
  </td>
</tr>
</table>
</div>


Therefore, throughout the core we have,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{\rho^*}{P^*}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{1 / 2}  \, ;</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\frac{M_r^*}{r^*}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\biggl( \frac{2\cdot 3}{\pi } \biggr)^{1/2} \biggl[ \xi^3 \biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{-3/2} \biggr]\biggl( \frac{2\pi}{3}\biggr)^{1 / 2} \frac{1}{\xi} 
= 2 \xi^2 \biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{-3/2} 
\, ;
</math>
  </td>
</tr>
</table>

and, throughout the envelope we have,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{\rho^*}{P^*}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\biggl( \frac{\mu_e}{\mu_c} \biggr) \theta^{-1}_i \phi(\eta)^{-1}
\, ;
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\frac{M_r^*}{r^*}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\biggl( \frac{\mu_e}{\mu_c} \biggr)^{-2} \theta^{-1}_i \biggl( \frac{2}{\pi} \biggr)^{1/2} \biggl(-\eta^2 \frac{d\phi}{d\eta} \biggr)
\biggl[ \biggl( \frac{\mu_e}{\mu_c} \biggr)^{-1} \theta^{-2}_i (2\pi)^{-1/2}\eta \biggr]^{-1}
=
2 \biggl( \frac{\mu_e}{\mu_c} \biggr)^{-1} \frac{\theta_i}{\eta} \biggl(-\eta^2 \frac{d\phi}{d\eta} \biggr)
\, .
</math>
  </td>
</tr>
</table>

<table border="1" align="center" width="80%" cellpadding="8"><tr><td align="left">
<font color="red">'''NOTE on 15 May 2019:'''</font>&nbsp; Prior to this date the last RHS expression had an incorrect exponent on <math>~\eta</math>.  It previously (incorrectly) read,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{M_r^*}{r^*}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
2 \biggl( \frac{\mu_e}{\mu_c} \biggr)^{-1} \theta_i  \eta \biggl(-\eta^2 \frac{d\phi}{d\eta} \biggr)
\, .
</math>
  </td>
</tr>
</table>

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

==Numerical Integration==

===General Approach===
Here, we begin by recognizing that the 2<sup>nd</sup>-order ODE that must be integrated to obtain the desired eigenvectors has the generic form,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~0</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
x'' + \frac{\mathcal{H}}{r^*} x' + \mathcal{K}x \, ,
</math>
  </td>
</tr>
</table>
where,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~x'</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\frac{dx}{dr^*}</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp; </td>
  <td align="right">
<math>~x''</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\frac{d^2x}{d(r^*)^2} \, .</math>
  </td>
</tr>
</table>
Adopting the same approach [[SSC/Stability/Polytropes#Numerical_Integration_from_the_Center.2C_Outward|as before when we integrated the LAWE for pressure-truncated polytropes]], we will enlist the finite-difference approximations,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~x'</math>
  </td>
  <td align="center">
<math>~\approx</math>
  </td>
  <td align="left">
<math>~
\frac{x_+ - x_-}{2\delta r^*}
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp; </td>
  <td align="right">
<math>~x''</math>
  </td>
  <td align="center">
<math>~\approx</math>
  </td>
  <td align="left">
<math>~
\frac{x_+ -2x_j + x_-}{(\delta r^*)^2} \, .
</math>
  </td>
</tr>
</table>

The finite-difference representation of the LAWE is, therefore,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{x_+ -2x_j + x_-}{(\delta r^*)^2}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
-~ \frac{\mathcal{H}}{r^*} \biggl[ \frac{x_+ - x_-}{2\delta r^*} \biggr] ~-~ \mathcal{K}x_j 
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~ x_+ -2x_j + x_-</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
-~ \frac{\delta r^*}{2r^*} \biggl[ x_+ - x_- \biggr]\mathcal{H} ~-~ (\delta r^*)^2\mathcal{K}x_j 
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~ x_{j+1} \biggl[1 + \biggl( \frac{\delta r^*}{2r^*}\biggr) \mathcal{H} \biggr]</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\biggl[ 2 - (\delta r^*)^2\mathcal{K}\biggr] x_j ~-~\biggl[ 1 - \biggl( \frac{\delta r^*}{2r^*} \biggr) \mathcal{H} \biggr]x_{j-1} \, .
</math>
  </td>
</tr>
</table>

In what follows we will also find it useful to rewrite <math>~\mathcal{K}</math> in the form,
<div align="center">
<math>~\mathcal{K} ~\rightarrow ~\biggl(\frac{\sigma_c^2}{\gamma_\mathrm{g}}\biggr) \mathcal{K}_1 - \alpha_\mathrm{g} \mathcal{K}_2 \, .</math>
</div>

The relevant coefficient expressions for ''all'' regions of the configuration are,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\mathcal{H}</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~
\biggl\{ 4 -\biggl(\frac{\rho^*}{P^*}\biggr)\frac{ M_r^*}{(r^*)}\biggr\}
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; , &nbsp; &nbsp; &nbsp; </td>
  <td align="right">
<math>~\mathcal{K}_1</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~
\frac{2\pi }{3}\biggl(\frac{\rho^*}{ P^* } \biggr) 
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp; </td>
  <td align="right">
<math>~\mathcal{K}_2</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~
\biggl(\frac{\rho^*}{ P^* } \biggr)\frac{M_r^*}{(r^*)^3} \, .
</math>
  </td>
</tr>
</table>

===Special Handling at the Center===
In order to kick-start the integration, we set the displacement function value to <math>~x_1 = 1</math> at the center of the configuration <math>~(\xi_1 = 0)</math>, then draw on the [[Appendix/Ramblings/PowerSeriesExpressions#PolytropicDisplacement|derived power-series expression]] to determine the value of the displacement function at the first radial grid line, <math>~\xi_2 = \delta\xi</math>, away from the center.  Specifically, we set,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~
x_2 
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
x_1 \biggl[ 1 - \frac{(n+1) \mathfrak{F}  (\delta\xi)^2}{60} \biggr]_{n=5} 
=
x_1 \biggl[ 1 - \frac{\mathfrak{F}  (\delta\xi)^2}{10} \biggr]
\, .</math>
  </td>
</tr>
</table>
</div>

===Special Handling at the Interface===

Integrating outward from the center, the ''general approach'' will work up through the determination of <math>~x_{j+1}</math> 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 <math>~x_+</math> ''would be'' &#8212; that is, to what the amplitude just beyond the interface ''would be'' &#8212; 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,

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

<tr>
  <td align="right">
<math>~x'_i\biggr|_\mathrm{core}</math>
  </td>
  <td align="center">
<math>~\approx</math>
  </td>
  <td align="left">
<math>~
\frac{1}{2\delta r^*} \biggl\{
x_+ - x_{i-1}
\biggr\}
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
-\frac{x_{i-1}}{2\delta r^*} 
+
\frac{1}{2\delta r^*} \biggl\{
\biggl[ 2 - (\delta r^*)^2\mathcal{K}\biggr] x_i ~-~\biggl[ 1 - \biggl( \frac{\delta r^*}{2r^*} \biggr) \mathcal{H} \biggr]x_{i-1} 
\biggr\}\biggl[1 + \biggl( \frac{\delta r^*}{2r^*}\biggr) \mathcal{H} \biggr]^{-1}
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{1}{2\delta r^*} \biggl\{
\biggl[ 2 - (\delta r^*)^2\mathcal{K}\biggr] x_i ~-~\biggl[ 1 - \biggl( \frac{\delta r^*}{2r^*} \biggr) \mathcal{H} \biggr]x_{i-1} ~-~\biggl[1 + \biggl( \frac{\delta r^*}{2r^*}\biggr) \mathcal{H} \biggr]x_{i-1} 
\biggr\}\biggl[1 + \biggl( \frac{\delta r^*}{2r^*}\biggr) \mathcal{H} \biggr]^{-1}
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{1}{2\delta r^*} \biggl\{
\biggl[ 2 - (\delta r^*)^2\mathcal{K}\biggr] x_i ~-~2x_{i-1} 
\biggr\}\biggl[1 + \biggl( \frac{\delta r^*}{2r^*}\biggr) \mathcal{H} \biggr]^{-1}
</math>
  </td>
</tr>
</table>

Conversely, as viewed from the ''envelope'', if we assume that we know <math>~x_i</math> and <math>~x'_i</math>, we can determine the amplitude, <math>~x_{i+1}</math>, at the first zone beyond the interface as follows:
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~x_-</math>
  </td>
  <td align="center">
<math>~\approx</math>
  </td>
  <td align="left">
<math>~
x_{i+1} - 2\delta r^*\cdot x'_i\biggr|_\mathrm{env}
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~ x_{i+1} \biggl[1 + \biggl( \frac{\delta r^*}{2r^*}\biggr) \mathcal{H} \biggr]</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\biggl[ 2 - (\delta r^*)^2\mathcal{K}\biggr] x_i ~-~\biggl[ 1 - \biggl( \frac{\delta r^*}{2r^*} \biggr) \mathcal{H} \biggr]
\biggl[ x_{i+1} - 2\delta r^*\cdot x'_i\biggr|_\mathrm{env}
\biggr]
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~ x_{i+1} \biggl[1 + \biggl( \frac{\delta r^*}{2r^*}\biggr) \mathcal{H} \biggr]
~+~
\biggl[ 1 - \biggl( \frac{\delta r^*}{2r^*} \biggr) \mathcal{H} \biggr]
x_{i+1} 
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\biggl[ 2 - (\delta r^*)^2\mathcal{K}\biggr] x_i 
~+~
\biggl[ 1 - \biggl( \frac{\delta r^*}{2r^*} \biggr) \mathcal{H} \biggr]
2\delta r^*\cdot x'_i\biggr|_\mathrm{env} 
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~ x_{i+1} 
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\biggl[ 1 - \tfrac{1}{2}(\delta r^*)^2\mathcal{K}\biggr] x_i 
~+~
\biggl[ 1 - \biggl( \frac{\delta r^*}{2r^*} \biggr) \mathcal{H} \biggr]
\delta r^*\cdot x'_i\biggr|_\mathrm{env} 
</math>
  </td>
</tr>
</table>

==Splitting Analysis Into Separate Core and Envelope Components==
===Core:===
Given that, <math>~\sqrt{2\pi/3}~r^* = \xi</math>, lets multiply the LAWE through by <math>~3/(2\pi)</math>.  This gives,

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

<tr>
  <td align="right">
<math>~0</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{d^2x}{d\xi^2} + \biggl\{ 4 -\biggl(\frac{\rho^*}{P^*}\biggr)\frac{ M_r^*}{(r^*)}\biggr\}\frac{1}{\xi} \cdot \frac{dx}{d\xi} 
+ \frac{3}{2\pi}\biggl(\frac{\rho^*}{ P^* } \biggr)\biggl\{ \frac{2\pi \sigma_c^2}{3\gamma_\mathrm{g}}  ~-~\frac{\alpha_\mathrm{g} M_r^*}{(r^*)^3}\biggr\}  x \, .
</math>
  </td>
</tr>
</table>

Specifically for the core, therefore, the finite-difference representation of the LAWE is, 
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{x_+ -2x_j + x_-}{(\delta \xi)^2}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
-~ \frac{\mathcal{H}}{\xi} \biggl[ \frac{x_+ - x_-}{2\delta \xi} \biggr] ~-~ \biggl[ \frac{3\mathcal{K}}{2\pi} \biggr]x_j 
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~ x_+ -2x_j + x_-</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
-~ \frac{\delta \xi}{2\xi} \biggl[ x_+ - x_- \biggr]\mathcal{H} ~-~ (\delta \xi)^2 \biggl[ \frac{3\mathcal{K}}{2\pi} \biggr] x_j 
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~ x_{j+1} \biggl[1 + \biggl( \frac{\delta \xi}{2\xi}\biggr) \mathcal{H} \biggr]</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\biggl[ 2 - (\delta \xi)^2\biggl( \frac{3\mathcal{K}}{2\pi} \biggr) \biggr] x_j ~-~\biggl[ 1 - \biggl( \frac{\delta \xi}{2\xi} \biggr) \mathcal{H} \biggr]x_{j-1} \, .
</math>
  </td>
</tr>
</table>
This also means that, as viewed from the perspective of the core, the slope at the interface is

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

<tr>
  <td align="right">
<math>~\biggl[ \frac{dx}{d\xi}\biggr]_\mathrm{interface}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{1}{2\delta \xi} \biggl\{
\biggl[ 2 - (\delta \xi)^2 \biggl( \frac{3\mathcal{K}}{2\pi} \biggr)\biggr] x_i ~-~2x_{i-1} 
\biggr\}\biggl[1 + \biggl( \frac{\delta \xi}{2\xi}\biggr) \mathcal{H} \biggr]^{-1} \, .
</math>
  </td>
</tr>
</table>

===Envelope:===
Given that,
<div align="center">
<math>~\biggl( \frac{\mu_e}{\mu_c} \biggr) \theta_i^2 (2\pi)^{1 / 2}~r^* = \eta \, ,</math>
</div>
let's multiply the LAWE through by <math>~(2\pi)^{-1} \theta_i^{-4}( \mu_e/\mu_c)^{-2} </math>.  This gives,

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

<tr>
  <td align="right">
<math>~0</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{d^2x}{d\eta^2} + \biggl\{ 4 -\biggl(\frac{\rho^*}{P^*}\biggr)\frac{ M_r^*}{(r^*)}\biggr\}\frac{1}{\eta} \cdot \frac{dx}{d\eta} 
+ \frac{1}{2\pi \theta_i^4} \biggl( \frac{\mu_e}{\mu_c} \biggr)^{-2} \biggl(\frac{\rho^*}{ P^* } \biggr)\biggl\{ \frac{2\pi \sigma_c^2}{3\gamma_\mathrm{g}}  ~-~\frac{\alpha_\mathrm{g} M_r^*}{(r^*)^3}\biggr\}  x \, .
</math>
  </td>
</tr>
</table>

Specifically for the envelope, therefore, the finite-difference representation of the LAWE is, 
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{x_+ -2x_j + x_-}{(\delta \eta)^2}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
-~ \frac{\mathcal{H}}{\eta} \biggl[ \frac{x_+ - x_-}{2\delta \eta} \biggr] ~-~  \biggl( \frac{\mu_e}{\mu_c} \biggr)^{-2}\biggl[ \frac{\mathcal{K}}{2\pi \theta_i^4} \biggr]x_j 
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~ x_+ -2x_j + x_-</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
-~ \frac{\delta \eta}{2\eta} \biggl[ x_+ - x_- \biggr]\mathcal{H} ~-~ (\delta \eta)^2 \biggl( \frac{\mu_e}{\mu_c} \biggr)^{-2}\biggl[ \frac{\mathcal{K}}{2\pi \theta_i^4} \biggr] x_j 
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~ x_{j+1} \biggl[1 + \biggl( \frac{\delta \eta}{2\eta}\biggr) \mathcal{H} \biggr]</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\biggl[ 2 - (\delta \eta)^2 \biggl( \frac{\mu_e}{\mu_c} \biggr)^{-2} \biggl( \frac{\mathcal{K}}{2\pi \theta_i^4} \biggr) \biggr] x_j ~-~\biggl[ 1 - \biggl( \frac{\delta \eta}{2\eta} \biggr) \mathcal{H} \biggr]x_{j-1} \, .
</math>
  </td>
</tr>
</table>

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,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~x_{i+1} 
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\biggl[ 1 - \tfrac{1}{2}(\delta \eta)^2 \biggl( \frac{\mu_e}{\mu_c} \biggr)^{-2} \biggl( \frac{\mathcal{K}}{2\pi \theta_i^4} \biggr)\biggr] x_i 
~+~
\biggl[ 1 - \biggl( \frac{\delta \eta}{2\eta} \biggr) \mathcal{H} \biggr]
\delta \eta \cdot \biggl[ \frac{dx}{d\eta} \biggr]_\mathrm{interface} \, .
</math>
  </td>
</tr>
</table>

===Interface===
If we consider only cases where <math>~\gamma_e = \gamma_c</math>, then at the interface we expect,

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

<tr>
  <td align="right">
<math>~\frac{d\ln x}{d\ln r^*}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\frac{d\ln x}{d\ln \xi} = \frac{d\ln x}{d\ln \eta}</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~ r^*\frac{dx}{d r^*}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\xi \frac{dx}{d \xi} = \eta \frac{d x}{d \eta}</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~ \frac{dx}{dr^*}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\biggl(\frac{2\pi}{3} \biggr)^{1 / 2}\frac{dx}{d\xi} = \biggl(\frac{\mu_e}{\mu_c}\biggr) \theta_i^2 (2\pi)^{1 / 2} \frac{dx}{d\eta} \, .</math>
  </td>
</tr>
</table>
Switching at the interface from <math>~\xi</math> to <math>~\eta</math> therefore means that,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~ \biggl[ \frac{dx}{d\eta}\biggr]_\mathrm{interface}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\frac{1}{\sqrt{3}} \biggl(\frac{\mu_e}{\mu_c}\biggr)^{-1} \theta_i^{-2} \biggl[ \frac{dx}{d\xi}\biggr]_\mathrm{interface} \, .</math>
  </td>
</tr>
</table>

If, however, we want to consider values for the adiabatic index that are different in the two regions, we have to follow the [[#Interface_Conditions|above-outlined guidelines]], that is,

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

<tr>
  <td align="right">
<math>~ \frac{d\ln x}{d\ln \eta} \biggr|_i</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~3\biggl(\frac{\gamma_c}{\gamma_e}  -1\biggr) + \frac{\gamma_c}{\gamma_e} \biggl( \frac{d\ln x}{d\ln \xi} \biggr)_i </math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~ \Rightarrow ~~~ \biggl[ \frac{dx}{d\eta} \biggr]_i</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\frac{x_i}{\eta_i} \biggl\{
3\biggl(\frac{\gamma_c}{\gamma_e}  -1\biggr) + \frac{\gamma_c}{\gamma_e} \cdot \frac{\xi_i}{x_i}\biggl[ \frac{dx}{d\xi} \biggr]_i 
\biggr\}</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{3x_i}{\eta_i} \biggl(\frac{\gamma_c}{\gamma_e}  -1\biggr) + \frac{\gamma_c}{\gamma_e} \cdot \frac{1}{\sqrt{3}} \biggl(\frac{\mu_e}{\mu_c}\biggr)^{-1} \theta_i^{-2} \biggl[ \frac{dx}{d\xi} \biggr]_i \, .
</math>
  </td>
</tr>
</table>
</div>

==Eigenvectors for Marginally Unstable Models with (&gamma;<sub>c</sub>, &gamma;<sub>e</sub>) = (6/5, 2)==

We now have the tools in hand to identify the eigenvectors &#8212; that is, various radial eigenfunctions and the corresponding eigenfrequency for each &#8212; associated with various modes of oscillation in <math>~(n_c, n_e) = (5,1)</math> bipolytropes.  Which models should we examine?  


In our [[SSC/Stability/MurphyFiedler85#Review_of_the_BiPolytrope_Stability_Analysis_by_Murphy_.26_Fiedler_.281985b.29|accompanying review of the bipolytrope stability analysis presented by Murphy &amp; Fiedler (1983b)]], our primary objective was to show that we were able to match their results quantitatively.  We therefore set <math>~\mu_e/\mu_c</math> = 1 &#8212; the only <math>~\mu</math>-ratio that they considered &#8212; and picked values of the core-envelope interface radius, <math>~\xi_i</math>, that were listed among their set of chosen models.  For a fixed value of <math>~\xi_i</math>, we integrated the relevant LAWE from the center toward the surface for many different eigenfrequency <math>~(\sigma_c^2)</math> ''guesses'' until an eigenfunction was found whose behavior at the surface matched with high precision the physically justified surface boundary condition.


With the [[#Virial_Stability_Evaluation|above virial stability analysis]] in mind (see especially [[#Virial_Stability_Evaluation|Figure 3]]), here we have chosen to focus on models that reside along six of the equilibrium sequences that have already been analytically identified, above &#8212; specifically, the sequences for which <math>~\mu_e/\mu_c</math> = 1, ½, 0.345, ⅓, 0.309, and ¼ &#8212; and to examine oscillation modes under the assumption that,
<div align="center">
<math>~\gamma_c = \frac{(n_c+1)}{n_c} = \frac{6}{5} \, ,</math> &nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp; <math>~\gamma_e = \frac{(n_e+1)}{n_e} =2 \, .</math>
</div>

===Fundamental Modes===

We decided to examine, first, whether any model along each sequence marks a transition from dynamically stable to dynamically unstable configurations.  We accomplished this by setting <math>~\sigma_c^2</math> = 0, then integrating the relevant LAWE from the center toward the surface for many different ''guesses'' of the core-envelope interface radius until an eigenfunction with no radial nodes &#8212; ''i.e.,'' an eigenfunction associated with the fundamental mode of radial oscillation &#8212; was found whose behavior at the surface matched with high precision the physically desired surface boundary condition.  We were successful in this endeavor.  A marginally unstable model was identified on each of the six separate equilibrium sequences.  

====Equilibrium Properties of Marginally Unstable Models====
Table 2 summarizes some of the equilibrium properties of these six models.  For example, the second column of the table gives the value of the core-envelope interface radius, <math>~\xi_i</math>, associated with each marginally unstable model.  The table also lists:  the value of the model's dimensionless radius, <math>~R^*_\mathrm{surf}</math>, the key structural parameters, <math>~q</math> &amp; <math>~\nu</math>, and the central-to-mean density associated with each model; and in each case the dimensionless thermal energy <math>~(\mathfrak{s})</math> and dimensionless gravitational potential energy <math>~(\mathfrak{w})</math> associated, separately, with the core and the envelope.  Note that, once the pair of parameters, <math>~(\mu_e/\mu_c, \xi_i)</math>, has been specified, we can legitimately assign high-precision values to all of the other model parameters because they are [[SSC/Structure/BiPolytropes/Analytic51#Parameter_Values|analytically prescribed]].

<div align="center">
<table border="1" cellpadding="8" align="center">
<tr>
  <th align="center" colspan="10">
'''Table 2:'''  Properties of Marginally Unstable Bipolytropes Having<br /><br /><math>~(n_c, n_e) = (5, 1)</math> and <math>~(\gamma_c, \gamma_e) = (\tfrac{6}{5}, 2)</math><br /><br />Determined from Integration of the LAWE
  </th>
</tr>
<tr>
  <td align="center"><math>~\frac{\mu_e}{\mu_c}</math></td>
  <td align="center"><math>~\xi_i</math></td>
  <td align="center"><math>~R^*_\mathrm{surf}</math></td>
  <td align="center"><math>~q \equiv \frac{r_\mathrm{core}}{R_\mathrm{surf}}</math></td>
  <td align="center"><math>~\nu \equiv \frac{M_\mathrm{core}}{M_\mathrm{tot}}</math></td>
  <td align="center"><math>~\frac{\rho_c}{\bar\rho}</math></td>
  <td align="center"><math>~\mathfrak{s}_\mathrm{core}</math></td>
  <td align="center"><math>~\mathfrak{w}_\mathrm{core}</math></td>
  <td align="center"><math>~\mathfrak{s}_\mathrm{env}</math></td>
  <td align="center"><math>~\mathfrak{w}_\mathrm{env}</math></td>
</tr>
<tr>
  <td align="center">1</td>
  <td align="left">1.6686460157</td>
  <td align="right">2.139737</td>
  <td align="center">0.53885819</td>
  <td align="center">0.497747626</td>
  <td align="center">8.51704656</td>
  <td align="center">3.021916335</td>
  <td align="center">-3.356583022</td>
  <td align="center">1.47780476</td>
  <td align="center">-5.642859167</td>
</tr>
<tr>
  <td align="center"><math>~\tfrac{1}{2}</math></td>
  <td align="left">2.27925811317</td>
  <td align="right">5.146499</td>
  <td align="center">0.306021732</td>
  <td align="center">0.401776274</td>
  <td align="center">63.29514949</td>
  <td align="center">4.241287117</td>
  <td align="center">-6.074241035</td>
  <td align="center">4.284931508</td>
  <td align="center">-10.97819621</td>
</tr>
<tr>
  <td align="center">0.345</td>
  <td align="left">2.560146865247</td>
  <td align="right">9.554041</td>
  <td align="center">0.185160563</td>
  <td align="center">0.234302525</td>
  <td align="center">209.7739052</td>
  <td align="center">4.639705843</td>
  <td align="center">-7.125754184</td>
  <td align="center">11.72861751</td>
  <td align="center">-25.61089252</td>
</tr>
<tr>
  <td align="center"><math>~\tfrac{1}{3}</math></td>
  <td align="left">2.582007485476</td>
  <td align="right">10.120558</td>
  <td align="center">0.176288391</td>
  <td align="center">0.218241608</td>
  <td align="center">230.4125398</td>
  <td align="center">4.667042505</td>
  <td align="center">-7.200966267</td>
  <td align="center">13.15887139</td>
  <td align="center">-28.45086152</td>
</tr>
<tr>
  <td align="center">0.309</td>
  <td align="left">2.6274239687695</td>
  <td align="center">11.464303</td>
  <td align="center">0.158362807</td>
  <td align="center">0.184796947</td>
  <td align="center">279.0788798</td>
  <td align="center">4.722277318</td>
  <td align="center">-7.354156963</td>
  <td align="center">17.1374434</td>
  <td align="center">-36.36528446</td>
</tr>
<tr>
  <td align="center"><math>~\tfrac{1}{4}</math></td>
  <td align="left">2.7357711469398</td>
  <td align="center">15.895632</td>
  <td align="center">0.118924863</td>
  <td align="center">0.11071211</td>
  <td align="center">430.0444648</td>
  <td align="center">4.84592201</td>
  <td align="center">-7.70305421</td>
  <td align="center">37.84289623</td>
  <td align="center">-77.67458196</td>
</tr>
</table>
</div>

<!--
OLD VERSION; IGNORE!

<div align="center">
<font color="red">'''OLD Table 2'''</font>

<table border="1" cellpadding="8" align="center">
<tr>
  <th align="center" colspan="9">
Properties of Marginally Unstable Bipolytropes Having<br /><br /><math>~(n_c, n_e) = (5, 1)</math> and <math>~(\gamma_c, \gamma_e) = (\tfrac{6}{5}, 2)</math><br /><br />Determined from Integration of the LAWE
  </th>
</tr>
<tr>
  <td align="center"><math>~\frac{\mu_e}{\mu_c}</math></td>
  <td align="center"><math>~\xi_i</math></td>
  <td align="center"><math>~q \equiv \frac{r_\mathrm{core}}{R_\mathrm{surf}}</math></td>
  <td align="center"><math>~\nu \equiv \frac{M_\mathrm{core}}{M_\mathrm{tot}}</math></td>
  <td align="center"><math>~\frac{\rho_c}{\bar\rho}</math></td>
  <td align="center"><math>~\mathfrak{s}_\mathrm{env}</math></td>
  <td align="center"><math>~\mathfrak{w}_\mathrm{env}</math></td>
  <td align="center"><math>~\mathfrak{s}_\mathrm{core}</math></td>
  <td align="center"><math>~\mathfrak{w}_\mathrm{core}</math></td>
</tr>
<tr>
  <td align="center">1</td>
  <td align="left">2.467359668</td>
  <td align="center">0.597684036</td>
  <td align="center">0.692367564</td>
  <td align="center">17.09749847</td>
  <td align="center">4.518031091</td>
  <td align="center">-6.79580606</td>
  <td align="center">0.857904827</td>
  <td align="center">-3.956065776</td>
</tr>
<tr>
  <td align="center"><math>~\tfrac{1}{2}</math></td>
  <td align="left">2.65925887</td>
  <td align="center">0.314759926</td>
  <td align="center">0.481811565</td>
  <td align="center">95.04044773</td>
  <td align="center">4.759771212</td>
  <td align="center">-7.459080087</td>
  <td align="center">3.20977926</td>
  <td align="center">-8.480020858</td>
</tr>
<tr>
  <td align="center">0.345</td>
  <td align="left">2.77457272408</td>
  <td align="center">0.183201946</td>
  <td align="center">0.25725514</td>
  <td align="center">281.7482802</td>
  <td align="center">4.887554727</td>
  <td align="center">-7.822410223</td>
  <td align="center">10.22665208</td>
  <td align="center">-22.40600339</td>
</tr>
<tr>
  <td align="center"><math>~\tfrac{1}{3}</math></td>
  <td align="left">2.7843424754</td>
  <td align="center">0.17412463</td>
  <td align="center">0.237677946</td>
  <td align="center">305.4881577</td>
  <td align="center">4.897826446</td>
  <td align="center">-7.852004624</td>
  <td align="center">11.6114224</td>
  <td align="center">-25.16649306</td>
</tr>
<tr>
  <td align="center">0.309</td>
  <td align="left">2.8050375512</td>
  <td align="center">0.156021514</td>
  <td align="center">0.197918988</td>
  <td align="center">359.3276918</td>
  <td align="center">4.91930991</td>
  <td align="center">-7.914090174</td>
  <td align="center">15.49005222</td>
  <td align="center">-32.90463409</td>
</tr>
</table>
</div>

OLD VERSION; IGNORE!
-->

As was expected from our [[#What_to_Expect_for_Equilibrium_Configurations|above discussion of virial equilibrium conditions]], we found that to high precision for each of these equilibrium models,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>
(\mathfrak{w}_\mathrm{core}  
~+~\mathfrak{w}_\mathrm{env})
~+~2(\mathfrak{s}_\mathrm{core}  
~+~\mathfrak{s}_\mathrm{env})  
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~0 \, .</math>
  </td>
</tr>
</table>
<span id="Figure4">However, contrary to expectations,</span> in no case did we find that <math>~\mathfrak{s}_\mathrm{core}/\mathfrak{s}_\mathrm{env} = 5</math>.  That is to say, we found that ''none'' of the models lies on the (red-dashed) curve in the <math>~q-\nu</math> parameter space that separates stable from unstable models as defined by our [[#What_to_Expect_for_Equilibrium_Configurations|above free-energy-based stability analysis]]. The left-hand panel of Figure 4 shows this (red-dashed) demarcation curve; for all intents and purposes, it is a reproduction of the right-hand panel of [[#Virial_Stability_Evaluation|Figure 3, above]] &#8212; turning-point markers have been removed to minimize clutter, the equilibrium sequences have been labeled, and the horizontal axis has been extended to unity in order to include a longer portion of the <math>~\mu_e/\mu_c = 1</math> sequence.  The orange triangular markers that appear in the right-hand panel of Figure 4 pinpoint where each of the Table 2 "marginally unstable" models resides in this <math>~q-\nu</math> plane.  Clearly, all six of the orange triangles lie well off of &#8212; and to the ''stable'' side of &#8212; the red-dashed demarcation curve.  This discrepancy, which has resulted from our use of two separate approaches to stability analysis, will be discussed further and gratifyingly resolved, below.

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

<tr>
  <th align="center">Figure 4</th>
</tr>
<tr>
  <td align="center" colspan="10">[[File:NEWCompositeDlabeled.png|800px|Marginally unstable models]]</td>
</tr>
</table>

====Eigenfunction Details====

Here we examine some of the properties of the fundamental-mode eigenfunctions that we have found are associated with marginally unstable, <math>~(n_c, n_e) = (5,1)</math> bipolytropes.
<table border="0" align="right" width="40%">
<tr>
  <th align="center">Figure 5</th>
</tr>
<tr><td align="center">
[[File:Mod0MuRatio100.png|450px|Example eigenvector]]
</td></tr>
</table>
Consider the model on the <math>~\mu_e/\mu_c = 1</math> sequence for which <math>~\sigma_c^2=0~</math>; key properties of this specific equilibrium model are enumerated in the first row of numbers provided in [[#Equilibrium_Properties_of_Marginally_Unstable_Models|Table 2, above]].  Figure 5 shows how our numerically derived, fundamental-mode eigenfunction, <math>~x = \delta r/r_0</math>, varies with the fractional radius over the entire range, <math>~0 \le r/R \le 1</math>.  By prescription, the eigenfunction has a value of unity and a slope of zero at the center  <math>~(r/R = 0)</math>.  Integrating the LAWE outward from the center, through the model's core (blue curve segment), <math>~x</math> drops smoothly to the value <math>~x_i = 0.81437</math> at the interface <math>~(\xi_i = 1.6686460157 ~\Rightarrow~ q = r_\mathrm{core}/R_\mathrm{surf} = 0.53885819)</math>.  Our numerical integration of the LAWE showed that, at the interface, the logarithmic slope of the ''core'' (blue) segment of the eigenfunction is,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~
\biggl\{ \frac{d\ln x}{d\ln r}\biggr|_i \biggr\}_\mathrm{core}
=
\biggl\{ \frac{d\ln x}{d\ln \xi}\biggr|_i \biggr\}_\mathrm{core}
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~- 0.455872 \, .</math>
  </td>
</tr>
</table>
Next, following the [[#Interface|above discussion of matching conditions at the interface]], we determined that, from the perspective of the envelope, the slope of the eigenfunction at the interface must therefore be,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~
\biggl\{ \frac{d\ln x}{d\ln r}\biggr|_i \biggr\}_\mathrm{env}
=
\biggl\{ \frac{d\ln x}{d\ln \eta}\biggr|_i \biggr\}_\mathrm{env}
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~3\biggl(\frac{\gamma_c}{\gamma_e}  -1\biggr) + \frac{\gamma_c}{\gamma_e} \biggl\{ \frac{d\ln x}{d\ln \xi}\biggr|_i \biggr\}_\mathrm{core} = -1.47352 \, .</math>
  </td>
</tr>
</table>

Adopting this "env" slope along with the amplitude, <math>~x_i = 0.81437</math>, as the appropriate ''interface'' boundary conditions, we integrated the LAWE from the interface to the surface, obtaining the green-colored segment of the eigenfunction that is shown in Figure 5.  The amplitude continued to steadily decrease, reaching a value of <math>~x_s = 0.38203</math>, at the model's surface <math>~(r/R = 1)</math>.  At the surface, this ''envelope'' (green) segment of the eigenfunction exhibits a logarithmic slope that matches to eight significant digits the value that is [[#SurfaceCondition|expected from astrophysical arguments]] for this marginally unstable <math>~(\sigma_c^2=0)</math> model, namely,
<div align="center">
<math>~
\frac{d\ln x}{d\ln \eta}\biggr|_s = \biggl[ \biggl( \frac{\rho_c}{\bar\rho} \biggr)\frac{\cancelto{0}{\sigma_c^2}}{2\gamma_e} - \biggl(3 - \frac{4}{\gamma_e}\biggr)\biggr] = -1 \, .
</math>
</div>


<table border="1" cellpadding="10" width="80%" align="center"><tr><td align="left">
<font color="red">'''Key Reminder:'''</font>  We were able to find an eigenfunction whose surface boundary condition matched the desired value &#8212; in this particular case, a logarithmic slope of negative one &#8212; to this high level of precision only by iterating many times and, at each step, fine-tuning our choice of the equilibrium model's radial interface location, <math>~\xi_i</math> before performing a numerical integration of the LAWE.
</td></tr></table>


The discontinuous jump that occurs in the ''slope'' of the eigenfunction at the interface results from our assumption that the effective adiabatic index of material in the core <math>~(\gamma_c = 6/5)</math> is different from the effective adiabatic index of the envelope material <math>~(\gamma_e = 2)</math>.   In an effort to emphasize and more clearly illustrate the behavior of this fundamental-mode eigenfunction as it crosses the core/envelope interface, we have added a pair of dashed line segments to the Figure 5 plot.  The red-dashed line segment touches, and is tangent to, the blue segment of the eigenfunction at the location of the core/envelope interface; it has a slope,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~
\frac{dx}{d(r/R)}\biggr|_i
=
\frac{x_i}{(r_i/R)}\biggl\{ \frac{d\ln x}{d\ln r}\biggr|_i\biggr\}_\mathrm{core}
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~- 0.455872 \biggl(\frac{ 0.81437 }{ 0.53885819 }\biggr) = - 0.68895\, .</math>
  </td>
</tr>
</table>
On the other hand, the purple-dashed line segment  touches, and is tangent to, the green segment of the eigenfunction at the location of the core/envelope interface; it has a slope,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~
\frac{dx}{d(r/R)}\biggr|_i
=
\frac{x_i}{(r_i/R)}\biggl\{ \frac{d\ln x}{d\ln r}\biggr|_i\biggr\}_\mathrm{env}
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~- 1.47352 \biggl(\frac{ 0.81437 }{ 0.53885819 }\biggr) = - 2.22691\, .</math>
  </td>
</tr>
</table>
For comparison purposes, the eigenfunction shown in Figure 5 has been presented again in Figure 6, along with several other of our numerically derived eigenfunctions, but in Figure 6 the plotted amplitude has been renormalized to give a surface value &#8212; rather than a central value &#8212; of unity.


In Figure 6 we show the behavior of the fundamental-mode eigenfunction for each of the marginally unstable models identified in Table 2.  In the top figure panel, each curve shows &#8212; on a linear-linear plot &#8212; how the amplitude varies with radius; in the bottom figure panel, the amplitude is plotted on a logarithmic scale.  On each curve, the black plus sign marks the radial location of the core-envelope interface; in the bottom panel, these markers are accompanied by the values of <math>~\xi_i</math> that are associated with each corresponding model (see also the second column of Table 2).  Each eigenfunction has been normalized such that the surface amplitude is unity.  In the top panel, the value of the central amplitude of the eigenfunction that results from this normalization is recorded near the point where each eigenfunction touches the vertical axis.   (In each case, the value provided on the plot is simply the inverse of the value of <math>~x_s</math> given in Table 3, below.)

<table border="0" align="center" cellpadding="8">
<tr>
  <th align="center">
[[File:DataFileButton02.png|right|60px|file = Dropbox/WorkFolder/Wiki edits/BiPolytrope/LinearPerturbation/FaulknerBipolytrope1.xlsx --- worksheet = Mode0Ensemble]]'''Figure 6:  Eigenfunctions Associated with the Fundamental-Mode of Radial Oscillation'''<br /> 
'''in Marginally Unstable Models having Various''' <math>~\mu_e/\mu_c</math>
  </th>
</tr>
<tr>
  <td align="center">
[[File:Mode0EigenfunctionsCombinedSmall.png|800px|Eigenfunctions for Marginally Unstable Models]]
  </td>
</tr>
</table>

Notice that, especially as they approach the surface, the "envelope" segments of these six marginally unstable eigenfunction appear to merge into the same curve, irrespective of their value of the ratio of mean molecular weights.  Note as well that the discontinuous jump that occurs in the ''slope'' of each eigenfunction at the radial location of the core/envelope interface &#8212; resulting from our choice to adopt a different adiabatic index, <math>~\gamma_g</math>, in the core from the one in the envelope &#8212; becomes less and less noticeable for smaller and smaller values of the ratio of mean molecular weights.

====Is There an Analytic Expression for the Eigenfunction?====
After noticing that, in Figure 6,  the ''envelope'' segments of all of the marginally unstable eigenfunctions merge into the same curve, we began to wonder whether a single expression &#8212; and, even better, an ''analytically defined'' expression &#8212; would perfectly describe the eigenfunction.  We had reason to believe that this might actually be possible because, in [[SSC/Stability/InstabilityOnsetOverview#Analyses_of_Radial_Oscillations|pressure-truncated polytropic configurations, we have derived analytic expressions for the marginally unstable, fundamental-mode eigenfunctions]] of both <math>~n = 5</math> and <math>~n=1</math> systems.

Very quickly, we convinced ourselves that a parabolic function does indeed perfectly match the "core" segment of each displayed eigenfunction.  Specifically, throughout the core <math>~(0 \le \xi \le \xi_i)</math>,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~x_P\biggr|_\mathrm{core}</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~1 - \frac{\xi^2}{15}</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow~~~\frac{dx_P}{d\xi}\biggr|_\mathrm{core}</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~- \frac{2\xi}{15} </math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow~~~\frac{d\ln x_P}{d\ln \xi}\biggr|_\mathrm{core}</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~- \frac{2\xi^2}{15} \biggl[ \frac{(15 - \xi^2)}{15} \biggr]^{-1} = - \frac{2\xi^2}{(15 - \xi^2)}  \, .</math>
  </td>
</tr>
</table>


The envelope segment posed a much greater challenge.  In the context of our [[SSC/Stability/n1PolytropeLAWE#Radial_Oscillations_of_n_.3D_1_Polytropic_Spheres|discussion of ''Radial Oscillations of n = 1 Polytropic Spheres'']], and in an [[Appendix/Ramblings/BiPolytrope51AnalyticStability#Is_There_an_Analytic_Expression_for_the_Eigenfunction.3F|accompanying ''Ramblings Appendix'' chapter]] we have detailed some trial derivations that are mostly blind alleyways.  Twice &#8212; once in [[SSC/Stability/n1PolytropeLAWE#tagJanuary2019|January, 2019]] and again (independently) in [[Appendix/Ramblings/BiPolytrope51AnalyticStability#Attempt_4B|April 2019]] &#8212; we have analytically demonstrated that the following appears to work for the envelope:

Given that the ''Structural Properties'' of the envelope are described by the Lane-Emden function,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\phi</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~a_0 \biggl[ \frac{\sin(\eta - b_0)}{\eta} \biggr] </math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow~~~ Q \equiv - \frac{d\ln \phi}{d\ln\eta}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\biggl[1 - \eta \cot(\eta - b_0) \biggr] \, ,</math>
  </td>
</tr>
</table>
the relevant LAWE is satisfied by the fractional displacement function,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~x_P</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\frac{3c_0 Q}{\eta^2} \, ,</math>
  </td>
</tr>
</table>
where, <math>~c_0</math> is an arbitrary scale factor.  

<table border="1" align="center" cellpadding="8" width="60%"><tr><td align="left">
Note that,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{dx_P}{d\eta}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{3c_0}{\eta^2}\biggl[\eta -\cot(\eta - b_0)
+\eta\cot^2(\eta - b_0) 
\biggr] 
- \frac{6c_0 Q}{\eta^3} 
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~\frac{d\ln x_P}{d\ln\eta} = \frac{\eta}{x_P}\cdot \frac{dx_P}{d\eta}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\frac{\eta^3}{3c_0 Q} \biggl\{
\frac{3c_0}{\eta^2}\biggl[\eta -\cot(\eta - b_0)
+\eta\cot^2(\eta - b_0) 
\biggr] 
- \frac{6c_0 Q}{\eta^3} 
\biggr\}
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\frac{1}{Q} \biggl\{
\biggl[\eta^2 -\eta \cot(\eta - b_0)
+\eta^2\cot^2(\eta - b_0) 
\biggr] 
- 2 Q 
\biggr\}
</math>
  </td>
</tr>

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

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\frac{1}{Q} 
\biggl[\eta^2  - Q + Q^2\biggr] 
- 2 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{\eta^2}{Q} + Q - 3 \, .
</math>
  </td>
</tr>
</table>

</td></tr></table>
But, as far as we have been able to determine (as of 16 April 2019), this analytic displacement function does not match the displacement function that has been generated through numerical integration of the LAWE (see the light-green segment of the eigenfunction displayed [[#Eigenfunction_Details|above in Figure 5]]).  It remains unclear whether (a) the numerical integration is at fault, (b) we are imposing an incorrect slope at the core-envelope interface, or ( c) we are misinterpreting how to compare the two separately derived (one, numerical, and the other, analytic) envelope eigenfunctions.

===Other Modes===
<!--
Keep in mind that, for all models, we ''expect'' that, at the surface, the logarithmic derivative of each proper eigenfunction will be,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{d\ln x}{d\ln r^*}\biggr|_\mathrm{surf}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\frac{\Omega^2}{\gamma} - \alpha \, .</math>
  </td>
</tr>
</table>

Also, keep in mind that, for Model 10 <math>~(\xi_i = 2.5646)</math>:
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{r_i}{R}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~0.39302482</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; , &nbsp;&nbsp;&nbsp;</td>
  <td align="right">
<math>~\frac{\rho_c}{\bar\rho}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~34.3460405</math>
  </td>
</tr>
</table>
-->

<table border="1" align="center" cellpadding="8">
<tr>
  <td align="center" colspan="15">
[[File:DataFileButton02.png|right|60px|file = Dropbox/WorkFolder/Wiki edits/BiPolytrope/LinearPerturbation/FaulknerBipolytrope1.xlsx --- worksheet = mu100plot]]'''Our Determinations for Marginally Unstable Model Having <math>~\mu_e/\mu_c = 1</math>'''<br />&nbsp;<br />
<math>~[ ~\xi_i = 1.6686460157 ~]</math>
  </td>
</tr>
<tr>
  <td align="center" rowspan="2">Mode</td>
  <td align="center" rowspan="2"><math>~\sigma_c^2</math></td>
  <td align="center" rowspan="2"><math>~\Omega^2 \equiv \frac{\sigma_c^2}{2} \biggl( \frac{\rho_c}{\bar\rho}\biggr)</math></td>
  <td align="center" rowspan="2"><math>~x_i</math></td>
  <td align="center" colspan="2"><math>~\frac{d\ln x}{d\ln r^*}\biggr|_i</math></td>
  <td align="center" rowspan="2"><math>~x_\mathrm{surf}</math></td>
  <td align="center" colspan="2"><math>~\frac{d\ln x}{d\ln r^*}\biggr|_\mathrm{surf}</math></td>
  <td align="center" rowspan="2"><math>~\frac{r}{R}\biggr|_1</math></td>
  <td align="center" rowspan="2"><math>~1 - \frac{M_r}{M_\mathrm{tot}}\biggr|_1</math></td>
  <td align="center" rowspan="2"><math>~\frac{r}{R}\biggr|_2</math></td>
  <td align="center" rowspan="2"><math>~1 - \frac{M_r}{M_\mathrm{tot}}\biggr|_2</math></td>
  <td align="center" rowspan="2"><math>~\frac{r}{R}\biggr|_3</math></td>
  <td align="center" rowspan="2"><math>~1 - \frac{M_r}{M_\mathrm{tot}}\biggr|_3</math></td>
</tr>
<tr>
  <td align="center">core</td>
  <td align="center">env</td>
  <td align="center">''expected''</td>
  <td align="center">measured</td>
</tr>
<tr>
  <td align="center">1<br /><font size="-1">(Fundamental)</font></td>
  <td align="right">0.00</td>
  <td align="right">0.00</td>
  <td align="right">+0.81437470</td>
  <td align="right">-0.455872</td>
  <td align="right">-1.473523</td>
  <td align="right">+0.3820</td>
  <td align="right">-1</td>
  <td align="right">-0.999999992</td>
  <td align="center">n/a</td>
  <td align="center">n/a</td>
  <td align="center">n/a</td>
  <td align="center">n/a</td>
  <td align="center">n/a</td>
  <td align="center">n/a</td>
</tr>
<tr>
  <td align="center">2</td>
  <td align="right">2.51513333</td>
  <td align="right">10.7107538</td>
  <td align="right">0.20482050</td>
  <td align="right">-7.09124</td>
  <td align="right">-5.4547441</td>
  <td align="right">- 0.9962</td>
  <td align="right">4.355376917</td>
  <td align="right">4.35537692</td>
  <td align="right">0.64133</td>
  <td align="right">0.3502</td>
  <td align="center">n/a</td>
  <td align="center">n/a</td>
  <td align="center">n/a</td>
  <td align="center">n/a</td>
</tr>
<tr>
  <td align="center">3</td>
  <td align="right">5.72371888</td>
  <td align="right">24.3745901</td>
  <td align="right">-0.14269277</td>
  <td align="right">+8.046019</td>
  <td align="right">+3.627611</td>
  <td align="right">+0.9308</td>
  <td align="right">11.18729505</td>
  <td align="right">11.18729506</td>
  <td align="right">0.4837</td>
  <td align="right">0.5864</td>
  <td align="right">0.842</td>
  <td align="right">0.0854</td>
  <td align="center">n/a</td>
  <td align="center">n/a</td>
</tr>
<tr>
  <td align="center">4</td>
  <td align="right">10.3458476</td>
  <td align="right">44.0622916</td>
  <td align="right">-0.20845197</td>
  <td align="right">-0.6949966</td>
  <td align="right">-1.61699793</td>
  <td align="right">-1.1443</td>
  <td align="right">21.03114578</td>
  <td align="right">21.03114577</td>
  <td align="right">0.3939</td>
  <td align="right">0.7154</td>
  <td align="right">0.6902</td>
  <td align="right">0.2777</td>
  <td align="center">0.9115</td>
  <td align="center">0.0284</td>
</tr>
<tr>
  <td align="center" colspan="15">
[[File:Mod0MuRatio100.png|550px|Our determination of eigenvector for mu_ratio = 1]]&nbsp;&nbsp;[[File:FourModesMuRatio100.png|550px|Our determination of multiple eigenvectors for mu_ratio = 1]]
  </td>
</tr>
</table>


{{ SGFfooter }}