Editing SSC/Stability/BiPolytropes/RedGiantToPN/Pt4 (section)

==Numerical Integration Through Envelope==

===Finite-Difference Expressions===

The discussion  in this subsection is guided by our [[SSC/Stability/Polytropes/Pt3#Numerical_Integration_from_the_Center,_Outward|previous attempt at numerical integration]].

Here, we focus on the LAWE that is relevant to the envelope, namely,
<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[4 - (n+1) Q_1 \biggr] \frac{1}{r_0} \cdot\frac{dx}{dr_0} 
- \biggl[(n+1)\alpha Q_1 \biggr]  \frac{x}{r_0^2} 
\, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{d^2x}{dr_0^2} + \biggl[4 - 2 Q_1 \biggr] \frac{1}{r_0} \cdot\frac{dx}{dr_0} 
- \biggl[2 Q_1 \biggr]  \frac{x}{r_0^2} 
\, ,
</math>
  </td>
</tr>
</table>
<span id="FD">where we have plugged</span> in the values, <math>(n,\alpha) = (1, 1)</math>.  Using the [[Appendix/Ramblings/NumericallyDeterminedEigenvectors#General_Approach|general finite-difference approach described separately]], we make the substitutions,

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

<tr>
  <td align="right">
<math>\biggl[\frac{dx}{dr_0}\biggr]_i</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>\biggl[\frac{d^2x}{dr_0^2}\biggr]_i</math>
  </td>
  <td align="center">
<math>\approx</math>
  </td>
  <td align="left">
<math>\frac{x_+ - 2x_i + x_-}{\Delta_r^2} \, ;</math>
  </td>
</tr>
</table>
which will provide 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>.

<font color="orange"><b>A:</b></font> &nbsp; &nbsp; Pick <math>\xi_\mathrm{int}</math>; this will give analytic expressions for <math>\eta_\mathrm{int}</math>, <math>B</math>, and for <math>\eta_\mathrm{surf}</math>, as well as analytic expressions for <math>(r_0)_\mathrm{int}</math> and <math>(r_0)_\mathrm{surf}</math>.

<font color="orange"><b>B:</b></font> &nbsp; &nbsp; Divide the radial coordinate grid into 99 spherical shells <math>\Rightarrow~ \Delta_r = [(r_0)_\mathrm{surf} - (r_0)_\mathrm{int}]/99.</math>  Then tabulate 100 values of <math>(r_0)_i, \eta_i, (Q_1)_i = [1 -  \eta\cot(\eta-B) ]_i</math>.

Generally speaking, after multiplying through by <math>r_0^2</math>, the finite-difference  representation of the envelope's LAWE takes the 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>
r_0^2\biggl[\frac{x_+ - 2x_i + x_-}{\Delta_r^2}\biggr] 
+ \biggl[4 - 2 Q_1 \biggr] r_0 \biggl[\frac{x_+ - x_-}{2 \Delta_r}\biggr] 
- \biggl[2 Q_1 \biggr]  x_i 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
x_+ \biggl\{ \frac{r_0^2}{\Delta_r^2} + (4-2Q_1)\frac{r_0}{2 \Delta_r} \biggr\}
+
x_i \biggl\{- \frac{2r_0^2}{\Delta_r^2} - 2Q_1  \biggr\}
+
x_- \biggl\{ \frac{r_0^2}{\Delta_r^2} - (4-2Q_1) \frac{r_0}{2 \Delta_r} \biggr\}
</math>
  </td>
</tr>
</table>
Multiplying through by <math>(\Delta_r^2/r_0^2)</math> and solving for <math>x_+</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>
x_+ \biggl\{ 1 + (4-2Q_1)\frac{\Delta_r}{2r_0} \biggr\}
-
2x_i \biggl\{1 + Q_1\biggl( \frac{\Delta_r^2}{r_0^2} \biggr)  \biggr\}
+
x_- \biggl\{ 1 - (4-2Q_1) \frac{\Delta_r}{2r_0} \biggr\}
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow ~~~ x_+ \biggl\{ 1 + (4-2Q_1)\frac{\Delta_r}{2r_0} \biggr\}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
2x_i \biggl\{1 + Q_1\biggl( \frac{\Delta_r^2}{r_0^2} \biggr)  \biggr\}
-
x_- \biggl\{ 1 - (4-2Q_1) \frac{\Delta_r}{2r_0} \biggr\}
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow ~~~ x_+ </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl\{~
2x_i \biggl[1 + Q_1\biggl( \frac{\Delta_r^2}{r_0^2} \biggr)  \biggr]
-
x_- \biggl[ 1 - (4-2Q_1) \frac{\Delta_r}{2r_0} \biggr]
~\biggr\}~\biggl[ 1 + (4-2Q_1)\frac{\Delta_r}{2r_0} \biggr]^{-1} \, .
</math>
  </td>
</tr>
</table>
Now, at the interface &#8212; as viewed from the perspective of both the core and the envelope &#8212; we know the value of <math>x_i =x_\mathrm{int}</math>, but we don't know the value of <math>x_-</math> as viewed from the envelope.  However &#8212; [[#STEPS|see <font color="maroon">STEP #4</font> below]] &#8212; we know analytically the value of the first derivative at the interface as viewed from the perspective of the envelope, namely, 
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>\biggl[ \frac{dx}{dr_0}\biggr]_\mathrm{int}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\frac{x_\mathrm{int}}{r_0} \cdot \biggl\{ \frac{d\ln x}{d\ln r_0}\biggr|_\mathrm{int} \biggr\}_\mathrm{env}</math>
  </td>
</tr>
</table>
Therefore, from the [[#FD|above-specified finite-difference representation]] of the first derivative, we deduce that,
<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>
x_+ - 2\Delta_r \cdot \biggl[ \frac{dx}{dr_0}\biggr]_\mathrm{int}
</math>
  </td>
</tr>
</table>
Hence, at the interface &#8212; and only ''at'' the interface &#8212; the finite-difference representation of the envelope's LAWE can be written as,

<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_+ \biggl[ 1 + (4-2Q_1)\frac{\Delta_r}{2r_0} \biggr]
-
2x_i \biggl[1 + Q_1\biggl( \frac{\Delta_r^2}{r_0^2} \biggr)  \biggr]
+
\biggl\{x_+ - 2\Delta_r \cdot \biggl[ \frac{dx}{dr_0}\biggr]_\mathrm{int} \biggr\}\cdot \biggl[ 1 - (4-2Q_1) \frac{\Delta_r}{2r_0} \biggr]
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
x_+ \biggl[ 1 + (4-2Q_1)\frac{\Delta_r}{2r_0} \biggr]
-
2x_i \biggl[1 + Q_1\biggl( \frac{\Delta_r^2}{r_0^2} \biggr)  \biggr]
+
x_+  \biggl[ 1 - (4-2Q_1) \frac{\Delta_r}{2r_0} \biggr]
- 2\Delta_r \cdot \biggl[ \frac{dx}{dr_0}\biggr]_\mathrm{int} \cdot \biggl[ 1 - (4-2Q_1) \frac{\Delta_r}{2r_0} \biggr]
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow~~~
x_+ \biggl[ 1 + (4-2Q_1)\frac{\Delta_r}{2r_0} \biggr]
+
x_+  \biggl[ 1 - (4-2Q_1) \frac{\Delta_r}{2r_0} \biggr]
</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
2x_i \biggl[1 + Q_1\biggl( \frac{\Delta_r^2}{r_0^2} \biggr)  \biggr]
+ 2\Delta_r \cdot \biggl[ \frac{dx}{dr_0}\biggr]_\mathrm{int} \cdot \biggl[ 1 - (4-2Q_1) \frac{\Delta_r}{2r_0} \biggr]
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow~~~
x_+ 
</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
x_i \biggl[1 + Q_1\biggl( \frac{\Delta_r^2}{r_0^2} \biggr)  \biggr]
+ \Delta_r \cdot \biggl[ \frac{dx}{dr_0}\biggr]_\mathrm{int} \cdot \biggl[ 1 - (4-2Q_1) \frac{\Delta_r}{2r_0} \biggr]
\, .
</math>
  </td>
</tr>
</table>

===Steps===
<table border="1" align="center" width="80%" cellpadding="8"><tr><td align="left">
<font color="maroon">STEP 1:</font>&nbsp;&nbsp;Specify the interface location from the perspective of the core; that is, specify <math>\xi_\mathrm{int}</math>, in which case,
<table border="0" align="center" cellpadding="8">
<tr>
  <td align="right">
<math>(r_0)_\mathrm{int} = a_5\cdot \xi_\mathrm{int}</math>
  </td>
  <td align="center">=</td>
  <td align="left">
<math>
\biggl[ K_5 G^{-1}\rho_c^{-4/5} \biggr]^{1 / 2}\biggl(\frac{3}{2\pi}\biggr)^{1 / 2} \xi_\mathrm{int} \, .
</math>
  </td>
</tr>
</table>

<font color="maroon">STEP 2:</font>&nbsp; &nbsp;Adopting the normalization <math>\phi_\mathrm{int} = 1</math>, determine numerous additional [[SSC/Structure/BiPolytropes/Analytic51#Step_6:_Envelope_Solution|equilibrium properties]] at the interface, such as &hellip;
<table border="0" align="center" cellpadding="8" width="80%">
<tr><td align="center" colspan="4">
<table border="1" align="center" cellpadding="8"><tr><td align="center"><font color="darkgreen">Example numerical values inside parentheses assume <math>(\mu_e/\mu_c) = 1</math> and <math>\xi_\mathrm{int} = 1.668646016</math><br /><math>\Rightarrow~~~(r_0)_\mathrm{int}[ K_5^{-1} G\rho_c^{4/5} ]^{1 / 2} = 1.153014872 \, .</math></td></tr></table>
  </td>
</tr>

<tr>
  <td align="right">
<math>\theta_\mathrm{int}</math>
  </td>
  <td align="center">=</td>
  <td align="left">
<math>
\biggl[1 + \frac{\xi^2_\mathrm{int}}{3}\biggr]^{-1 / 2} \, ;
</math>
  </td>
  <td align="right">(0.720165375)</td>
</tr>

<tr>
  <td align="right">
<math>\biggl( \frac{d\theta}{d\xi} \biggr)_\mathrm{int}</math>
  </td>
  <td align="center">=</td>
  <td align="left">
<math>
-~\frac{\xi_\mathrm{int}}{3}\biggl[1 + \frac{\xi^2_\mathrm{int}}{3}\biggr]^{-3 / 2} \, ;
</math>
  </td>
  <td align="right">(- 0.207749350)</td>
</tr>

<tr>
  <td align="right">
<math>\eta_\mathrm{int}</math>
  </td>
  <td align="center">=</td>
  <td align="left">
<math>
3^{1 / 2}\biggl(\frac{\mu_e}{\mu_c}\biggr)\theta^2_\mathrm{int}\cdot \xi_\mathrm{int} \, ;
</math>
  </td>
  <td align="right">(1.498957494)</td>
</tr>

<tr>
  <td align="right">
<math>\biggl( \frac{d\phi}{d\eta} \biggr)_\mathrm{int}</math>
  </td>
  <td align="center">=</td>
  <td align="left">
<math>
3^{1 / 2} \theta^{-3}_\mathrm{int}\cdot \biggl( \frac{d\theta}{d\xi} \biggr)_\mathrm{int} \, ;
</math>
  </td>
  <td align="right">(- 0.963393227)</td>
</tr>

<tr>
  <td align="right">
<math>\Lambda_\mathrm{int}</math>
  </td>
  <td align="center">=</td>
  <td align="left">
<math>
\frac{1}{\eta_\mathrm{int}} + \biggl( \frac{d\phi}{d\eta} \biggr)_\mathrm{int} \, ;
</math>
  </td>
  <td align="right">(- 0.296262902)</td>
</tr>

<tr>
  <td align="right">
<math>A</math>
  </td>
  <td align="center">=</td>
  <td align="left">
<math>
\eta_\mathrm{int}(1 + \Lambda^2_\mathrm{int})^{1 / 2} \, ;
</math>
  </td>
  <td align="right">(1.563357124)</td>
</tr>

<tr>
  <td align="right">
<math>B</math>
  </td>
  <td align="center">=</td>
  <td align="left">
<math>
\eta_\mathrm{int} - \frac{\pi}{2} + \tan^{-1}(\Lambda_\mathrm{int}) \, .
</math>
  </td>
  <td align="right">(- 0.359863580)</td>
</tr>

<tr>
  <td align="right">
<math>\eta_\mathrm{surf}</math>
  </td>
  <td align="center">=</td>
  <td align="left">
<math>
B + \pi \, .
</math>
  </td>
  <td align="right">(2.781729074)</td>
</tr>
</table>

<font color="maroon">STEP 3:</font>&nbsp;&nbsp; Throughout the core &#8212; that is, at all radial positions, <math>0 \le r_0 \le (r_0)_\mathrm{int}</math> &#8212; the displacement amplitude and its first and second derivatives, respectively, are given by the expressions,
<table border="0" cellpadding="5" align="center" width="80%">
<tr>
  <td align="right">
<math>x</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl[1 - \frac{r_0^2}{15a_5^2} \biggr] = \biggl[1 - \frac{\xi^2}{15} \biggr]
\, ;
</math>
  </td>
  <td align="right">(0.814374698)</td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow ~~~r_0\cdot \frac{dx}{dr_0}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
-~\frac{2r_0^2}{15a_5^2} = -~\frac{2\xi^2}{15}
\, ;
</math>
  </td>
  <td align="right">(- 0.371250604)</td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow ~~~ r_0^2 \cdot \frac{d^2x}{dr_0^2}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
-~\frac{2r_0^2}{15a_5^2} = -~\frac{2\xi^2}{15}
\, ;
</math>
  </td>
  <td align="right">(- 0.371250604)</td>
</tr>

<tr>
  <td align="right">
also &hellip; &nbsp; &nbsp; 
<math>
\biggl\{ \frac{d\ln x}{d\ln \xi} \biggr\}_\mathrm{core}
= 
\biggl\{ \frac{d\ln x}{d\ln r_0} \biggr\}_\mathrm{core} 
= \frac{r_0}{x} \cdot \frac{dx}{dr_0}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl[\frac{15}{15 - \xi^2} \biggr] \cdot \biggl[-~\frac{2\xi^2}{15 }\biggr]
= 
\biggl[\frac{2\xi^2}{\xi^2 - 15} \biggr]
\, .
</math>
  </td>
  <td align="right">(-0.455871977)<sup>&dagger;</sup></td>
</tr>
</table>

<font color="maroon">STEP #4:</font>&nbsp; &nbsp; From the determination of the logarithmic slope of the displacement function at the edge of the core &#8212; <i>i.e.,</i> at the core-envelope interface &#8212; determine the slope as viewed from the perspective of the envelope.
<table border="0" cellpadding="5" align="center" width="80%">

<tr>
  <td align="right">
<math>~
\biggl\{ \frac{d\ln x}{d\ln r_0}\biggr|_\mathrm{int} \biggr\}_\mathrm{env}
=
\biggl\{ \frac{d\ln x}{d\ln \eta}\biggr|_\mathrm{int} \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|_\mathrm{int} \biggr\}_\mathrm{core} \, .</math>
  </td>
  <td align="right">(-1.473523186)<sup>&dagger;</sup></td>
</tr>
</table>

----

<sup>&dagger;</sup>This analytically determined value matches the [[SSC/Stability/BiPolytropes/Pt3#Eigenfunction_Details|previous determination]] that was obtained via numerical integration of the LAWE.
</td></tr></table>

[[SSC/Structure/BiPolytropes/Analytic51#Step_8:_Throughout_the_envelope_(ηi_≤_η_≤_ηs)|Throughout the envelope]] &#8212; that is, over the range, <math>(\eta_\mathrm{int} \le \eta \le \eta_\mathrm{surf})</math> &#8212; the radial coordinate, <math>r_0</math>, is a linear function of <math>\eta</math> and takes on values given by the expression,
<table border="0" cellpadding="5" align="center" width="80%">

<tr>
  <td align="right">
<math>
r_0 [K_5^{-1} G \rho_c^{4/5}]^{1 / 2}
</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl[ \biggl(\frac{\mu_e}{\mu_c}\biggr)^{-1} \theta^{-2}_\mathrm{int} (2\pi)^{-1 / 2} \biggr]\cdot \eta
</math>
  </td>
  <td align="right">(0.769211186 &times; &eta;)</td>
</tr>

<tr>
  <td align="right">
<math>
\Rightarrow ~~~1.153014872
</math>
  </td>
  <td align="center">
<math>\leq r_0 \leq</math>
  </td>
  <td align="left">
<math>2.139737121
\, .
</math>
  </td>
  <td align="right">&nbsp;</td>
</tr>
</table>



[[SSC/Stability/BiPolytropes/Pt3#Eigenfunction_Details|From our earlier discussions]],


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>