Editing Appendix/Ramblings/BiPolytropeStability (section)

==Our Numerical Integration==

Let's try to integrate this bipolytrope's LAWE from the center, outward, using as a guideline an [[SSC/Stability/Polytropes#Numerical_Integration_from_the_Center.2C_Outward|accompanying ''Numerical Integration'' outline]].  Generally, for any polytropic index, the relevant LAWE can be written in the form,

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

<tr>
  <td align="right">
<math>~\theta_i {x_i''}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~- \biggl[\mathcal{A} \biggr] \frac{x_i'}{\xi_i}  
- \frac{(n+1)}{6} \biggl[ \mathcal{B} \biggr]  x_i </math>
  </td>
</tr>
</table>
where,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\mathcal{A}</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~
4\theta_i - (n+1)\xi_i (- \theta^')_i 
=
\theta_i [ 4 - (n+1)Q_i] 
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\mathcal{B}</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~
\frac{\sigma_c^2}{\gamma_g}  - 2\alpha \biggl(- \frac{3\theta^'}{\xi} \biggr)_i 
=
\mathfrak{F} + 2\alpha \biggl[ 1 - \biggl(- \frac{3\theta^'}{\xi} \biggr)_i  \biggr]
=
\mathfrak{F} + 2\alpha \biggl[ 1 - \frac{3\theta_i}{\xi_i^2} \cdot Q_i  \biggr]
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~
\mathfrak{F}
</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~
\biggl[ \frac{\sigma_c^2}{\gamma_g}  - 2\alpha\biggr] 
=
\biggl[ \frac{\sigma_c^2}{\gamma_g}  - 2\biggl(3 - \frac{4}{\gamma_g} \biggr) \biggr] 
=
\biggl[ \frac{(8 + \sigma_c^2)}{\gamma_g}  - 6\biggr] 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~\Rightarrow~~~</math>
  </td>
  <td align="left">
<math>~
\sigma_c^2   =
\gamma_g (\mathfrak{F} + 6) -8
\, .
</math>
  </td>
</tr>
</table>

This leads to a discrete, finite-difference representation of the form,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~x_+ \biggl[2\theta_i + \frac{\delta\xi}{\xi_i} \cdot \mathcal{A}\biggr] </math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
x_- \biggl[\frac{\delta\xi }{\xi_i} \cdot \mathcal{A} - 2\theta_i\biggr] 
+  x_i\biggl\{4\theta_i - \frac{(\delta\xi)^2(n+1)}{3}\cdot \mathcal{B} \biggr\} \, .</math>
  </td>
</tr>
</table>
</div>

This provides an approximate expression for <math>~x_+ \equiv x_{i+1}</math>, given the values of <math>~x_- \equiv x_{i-1}</math> and <math>~x_i</math>; this works for all zones, <math>~i = 3 \rightarrow N</math>  as long as the center of the configuration is denoted by the grid index, <math>~i=1</math>.  Note that,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\delta\xi</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~\frac{\xi_\mathrm{max}}{(N - 1)} </math>
  </td>
  <td align="center">
&nbsp; &nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp; &nbsp; 
  </td>
  <td align="right">
<math>~\xi_i</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~(i-1)\delta\xi \, . </math>
  </td>
</tr>
</table>

In order to kick-start the integration, we will set the displacement function value to <math>~x_1 = 1</math> at the center of the configuration <math>~(\xi_1 = 0)</math>, then we will 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 will 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] \, .</math>
  </td>
</tr>
</table>
</div>

===Integration Through the n = 1 Core===

For an <math>~n = 1</math> core, we have,
<div align="center">
<math>
\theta_i = \frac{\sin\xi_i}{\xi_i}
</math>
&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp; 
<math>
Q_i = 1 -  \xi_i \cot\xi_i \, .
</math>
</div>
Hence,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\mathcal{A}_\mathrm{core}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{\sin\xi_i}{\xi_i} \biggl[ 4 - 2(1 -  \xi_i \cot\xi_i) \biggr] 
=
\frac{2\sin\xi_i}{\xi_i} \biggl[ 1 + \xi_i \cot\xi_i \biggr] 
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\mathcal{B}_\mathrm{core}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\mathfrak{F}_\mathrm{core} + 2\alpha_\mathrm{core} \biggl[ 1 - \frac{3\theta_i}{\xi_i^2} \cdot Q_i  \biggr]
=
\mathfrak{F}_\mathrm{core} + 2\alpha_\mathrm{core} \biggl[ 1 - \frac{3\sin\xi_i}{\xi_i^3} \biggl( 1 - \xi_i \cot\xi_i \biggr)\biggr] \, .
</math>
  </td>
</tr>
</table>

So, first we choose a value of <math>~\sigma_c^2</math> and <math>~\gamma_c</math>, which means,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~
\mathfrak{F}_\mathrm{core}
</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~
\biggl[ \frac{(8 + \sigma_c^2)}{\gamma_c}  - 6\biggr] 
</math>
  </td>
</tr>
</table>

Then, moving from the center of the configuration, outward to the interface at <math>~\xi_i = \xi_\mathrm{interface} ~~ \Rightarrow ~~\delta\xi = \xi_\mathrm{interface}/(N-1)</math>, we have,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~x_1</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~1 \, ,</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~
x_2 
</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
x_1 \biggl[ 1 - \frac{\mathfrak{F}_\mathrm{core}  (\delta\xi)^2}{30} \biggr] \, ,</math>
  </td>
</tr>

<tr>
  <td align="right">
for <math>~i = 2 \rightarrow N \, ,</math>&nbsp; &nbsp; &nbsp; &nbsp;
<math>~x_{i+1} \biggl[2\theta_i + \frac{\delta\xi}{\xi_i} \cdot \mathcal{A}_\mathrm{core} \biggr] </math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
x_{i-1} \biggl[\frac{\delta\xi }{\xi_i} \cdot \mathcal{A}_\mathrm{core} - 2\theta_i\biggr] 
+  x_i\biggl\{4\theta_i - \frac{(\delta\xi)^2(n+1)}{3}\cdot \mathcal{B}_\mathrm{core} \biggr\} \, .</math>
  </td>
</tr>
</table>
At the interface &#8212; that is, when <math>~i=N</math> &#8212; the logarithmic slope of the displacement function is,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{d\ln x}{d\ln\xi}\biggr|_\mathrm{interface}</math>
  </td>
  <td align="center">
<math>~\approx</math>
  </td>
  <td align="left">
<math>~
\frac{\xi_N}{x_N} \cdot \frac{(x_{N+1} - x_{N-1})}{2\delta\xi} \, .
</math>
  </td>
</tr>
</table>

===Interface===

Keep in mind that, as has been [[SSC/Structure/BiPolytropes/Analytic15#Parameter_Values|detailed in the accompanying ''equilibrium structure'' chapter]], for <math>~(n_c, n_e) = (1, 5)</math> bipolytropes,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~r^*</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\biggl( \frac{1}{2\pi}\biggr)^{1 / 2} \xi \, ,</math>
  </td>
  <td align="left">
&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; for, 
  </td>
  <td align="left">
<math>~0 \le \xi \le \xi_\mathrm{interface} \, .</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~r^*</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\biggl[ \biggl( \frac{\mu_e}{\mu_c}\biggr)^{-1} \biggl( \frac{3}{2\pi}\biggr)^{1 / 2} \biggr] \eta \, ,</math>
  </td>
  <td align="left">
&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; for, 
  </td>
  <td align="left">
<math>~\eta_\mathrm{interface} \le \eta \le \eta_s \, ;</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\eta_s</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{1}{\sqrt{3}}\biggl( \frac{\mu_e}{\mu_c}\biggr)  \xi_s = 
\frac{1}{\sqrt{3}}\biggl( \frac{\mu_e}{\mu_c}\biggr) \biggl[ \xi e^{2(\pi - \Delta)} \biggr]_\mathrm{interface}
\, .
</math>
  </td>
  <td align="left" colspan="2">
&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;  
  </td>
</tr>
</table>

We now need to determine what the slope is at the interface, viewed from the perspective of the envelope.  From [[#Interface_Conditions|above]], we deduce that,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{d\ln y}{d\ln \eta} \biggr|_\mathrm{interface}</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} \cdot \frac{d\ln x}{d\ln\xi}\biggr|_\mathrm{interface} \, .</math>
  </td>
</tr>
</table>

Hence, letting the subscript "1" denote the interface location as viewed from the envelope, we have,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{\eta_1}{y_1} \cdot \frac{(y_2 - y_0)}{ (\eta_2 - \eta_0)}</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} \cdot \frac{d\ln x}{d\ln\xi}\biggr|_\mathrm{interface} \, .</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow ~~~ y_0</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~y_2 - \frac{2 (\delta\eta) y_1}{\eta_1} \biggl\{ 3\biggl(\frac{\gamma_c}{\gamma_e}  -1\biggr) + \frac{\gamma_c}{\gamma_e} \cdot \frac{d\ln x}{d\ln\xi}\biggr|_\mathrm{interface} \biggr\} \, .</math>
  </td>
</tr>
</table>

===Integration Through the n = 5 Envelope===

For an <math>~n = 5</math> envelope, we have,
<div align="center">
<math>
\phi_i = \frac{B_0^{-1}\sin\Delta}{\eta^{1/2}(3-2\sin^2\Delta)^{1/2}}
</math>
&nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp; 
<math>
Q_i = - \frac{d\ln\phi}{d\ln\eta} = \frac{ 3\sin\Delta - 3\cos\Delta - 2\sin^3\Delta }{2 \sin\Delta (3-2\sin^2\Delta)}\, ,
</math>
</div>

where <math>~A_0</math> is a "homology factor," <math>~B_0</math> is an overall scaling coefficient, and we have introduced the notation,
<div align="center">
<math>~\Delta \equiv \ln(A_0\eta)^{1/2} = \frac{1}{2} (\ln A_0 + \ln\eta) \, .</math>
</div>

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

<tr>
  <td align="right">
<math>~\mathcal{A}_\mathrm{env}</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~
\phi_i [ 4 - (n+1)Q_i] 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{B_0^{-1}\sin\Delta}{\eta^{1/2}(3-2\sin^2\Delta)^{1/2}}\biggl\{ 4 
~-~ 6 \biggl[ \frac{ 3\sin\Delta - 3\cos\Delta - 2\sin^3\Delta }{2 \sin\Delta (3-2\sin^2\Delta)}  \biggr] \biggr\} 
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\mathcal{B}_\mathrm{env}</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~
\frac{\sigma_c^2}{\gamma_e}  - 2\alpha_\mathrm{env} \biggl(- \frac{3\phi^'}{\eta} \biggr)_i 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{1}{\gamma_e}\biggl[ \gamma_c (\mathfrak{F}_\mathrm{core} + 6) -8 \biggr]  - 2\biggl[ 3 - \frac{4}{\gamma_e}\biggr] Q_i \biggl( \frac{3\phi_i }{\eta_i^2} \biggr)
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{\gamma_c}{\gamma_e}\biggl[ \mathfrak{F}_\mathrm{core} + 6  -\frac{8}{\gamma_c} \biggr]  - 6\biggl[ 3 - \frac{4}{\gamma_e}\biggr] 
\frac{B_0^{-1}\sin\Delta}{\eta^{5/2}(3-2\sin^2\Delta)^{1/2}}
\biggl[ \frac{ 3\sin\Delta - 3\cos\Delta - 2\sin^3\Delta }{2 \sin\Delta (3-2\sin^2\Delta)} \biggr]
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\frac{\gamma_c}{\gamma_e}\biggl[ \mathfrak{F}_\mathrm{core} + 6  -\frac{8}{\gamma_c} \biggr]  - 3B_0^{-1}\biggl[ 3 - \frac{4}{\gamma_e}\biggr] 
\biggl[ \frac{ 3\sin\Delta - 3\cos\Delta - 2\sin^3\Delta }{ \eta^{5/2}(3-2\sin^2\Delta)^{3 / 2}} \biggr]
\, .
</math>
  </td>
</tr>
</table>

This leads to a discrete, finite-difference representation of the form,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~y_+ \biggl[2\phi_i + \frac{\delta\eta}{\eta_i} \cdot \mathcal{A}_\mathrm{env} \biggr] </math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
y_- \biggl[\frac{\delta\eta }{\eta_i} \cdot \mathcal{A}_\mathrm{env} - 2\phi_i\biggr] 
+  y_i\biggl[4\phi_i - 2(\delta\eta)^2 \cdot \mathcal{B}_\mathrm{env} \biggr] </math>
  </td>
</tr>
</table>
</div>

This provides an approximate expression for <math>~y_+ \equiv y_{i+1}</math>, given the values of <math>~y_- \equiv y_{i-1}</math> and <math>~y_i</math>; this works for all zones, <math>~i = 3 \rightarrow M</math>  as long as the interface between the core and the envelope of the configuration is denoted by the grid index, <math>~i=1</math>.  Note that,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\delta\eta</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~\frac{\eta_\mathrm{surf}- \eta_\mathrm{interface} }{M - 1} </math>
  </td>
  <td align="center">
&nbsp; &nbsp; &nbsp; &nbsp; and &nbsp; &nbsp; &nbsp; &nbsp; 
  </td>
  <td align="right">
<math>~\eta_i</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\eta_\mathrm{interface} + (i-1)\delta\eta \, . </math>
  </td>
</tr>
</table>

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 <math>~i = 1</math>, we must set, <math>~y_1 = x_N</math> and <math>~\eta_1 = (\mu_e/\mu_c)\xi_N/\sqrt{3}</math>.  Then the value of <math>~y_2</math> is obtained from the expression,

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

<tr>
  <td align="right">
<math>~y_2 \biggl[2\phi_1 + \frac{\delta\eta}{\eta_1} \cdot \mathcal{A}_\mathrm{env} \biggr] </math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
y_0 \biggl[\frac{\delta\eta }{\eta_1} \cdot \mathcal{A}_\mathrm{env} - 2\phi_1\biggr] 
+  y_1\biggl\{4\phi_1 - 2(\delta\eta)^2 \cdot \mathcal{B}_\mathrm{env} \biggr\} </math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
y_1\biggl[ 4\phi_1 - 2(\delta\eta)^2 \cdot \mathcal{B}_\mathrm{env} \biggr] 
+  
\biggl[\frac{\delta\eta }{\eta_1} \cdot \mathcal{A}_\mathrm{env} - 2\phi_1\biggr] 
\biggl\{
y_2 ~-~ \frac{2 (\delta\eta) y_1}{\eta_1} \biggl[ 3\biggl(\frac{\gamma_c}{\gamma_e}  -1\biggr) + \frac{\gamma_c}{\gamma_e} \cdot \frac{d\ln x}{d\ln\xi}\biggr|_\mathrm{interface} \biggr]
\biggr\}
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow~~~y_2 \biggl[4\phi_1  \biggr] </math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
y_1\biggl[ 4\phi_1 - 2(\delta\eta)^2 \cdot \mathcal{B}_\mathrm{env} \biggr] 
~-~  
\frac{2 (\delta\eta) y_1}{\eta_1} \biggl[\frac{\delta\eta }{\eta_1} \cdot \mathcal{A}_\mathrm{env} - 2\phi_1\biggr] 
\biggl\{3\biggl(\frac{\gamma_c}{\gamma_e}  -1\biggr) + \frac{\gamma_c}{\gamma_e} \cdot \frac{d\ln x}{d\ln\xi}\biggr|_\mathrm{interface} 
\biggr\} 
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>~\Rightarrow~~~\frac{y_2}{y_1} \biggl[\phi_1  \biggr] </math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\biggl[ \phi_1 - \frac{(\delta\eta)^2}{2} \cdot \mathcal{B}_\mathrm{env} \biggr] 
~-~  
\frac{ (\delta\eta) }{2\eta_1} \biggl[\frac{\delta\eta }{\eta_1} \cdot \mathcal{A}_\mathrm{env} - 2\phi_1\biggr] 
\biggl\{3\biggl(\frac{\gamma_c}{\gamma_e}  -1\biggr) + \frac{\gamma_c}{\gamma_e} \cdot \frac{d\ln x}{d\ln\xi}\biggr|_\mathrm{interface} 
\biggr\} \, .
</math>
  </td>
</tr>
</table>