SSC/Structure/Polytropes/Numerical

Polytropic Spheres

Lane (1870)

Part I:  Isolated Polytropes

Part II:  Known Analytic Solutions

Part III:  Numerical Solutions

Numerical Solutions

Here we explain how an Excel workbook can be used to numerically solve the Lane-Emden equation, evaluating the Lane-Emden function across a one-dimensional, discrete grid.

Techniques

HSCF Technique

via Self-Consistent Field (SCF) Technique

On the first spreadsheet within the workbook, we establish the following columns of number:

• Column A:   Labeled $r_i/R$ (for i between 1 and N), that represents a discrete radial grid of spacing, $\mathrm{\Delta }=(N-1{)}^{-1}$; each row gives the radial coordinate location of the i^th zone, starting from $r_1/R=0$ and ending at $r_N/R=1$.
• Column B:   Labeled ${\mathrm{r}\mathrm{h}\mathrm{f}}_i$ (for i between 1 and N-1); each row gives the radial coordinate of the mid-point of a grid zone.

${\mathrm{r}\mathrm{h}\mathrm{f}}_i\equiv \frac{1}{2}[\frac{r_i}{R}+\frac{r_{i+1}}{R}].$

• Column C:   Labeled ${\rho }_i$ (for i between 1 and N-1); each row provides an initial guess for the mass-density of the grid zone. Usually it is sufficient to guess, ${\rho }_i=1$ throughout. For an $n=0$ polytrope, this proves also to be the correct final density profile.
• Column D:   Labeled $M_i$ (for i between 1 and N); the i^th row gives the integrated mass enclosed interior to the radial grid coordinate, $r_i/R$. Specifically, $M_1=0$, and thereafter, beginning with zone, $i=2$,

$M_i=M_{i-1}+\frac{4\pi {\rho }_{i-1}}{3}\left[\right(\frac{r_i}{R}{)}^3-\left(\frac{r_{i-1}}{R}{)}^3\right].$

• Note that, $M_{\mathrm{t}\mathrm{o}\mathrm{t}}=M_N.$
• Column E:   Labeled $g_i$ (for i between 2 and N); each row tabulates the inwardly directed gravitational acceleration that is felt at the outer edge of each grid zone. Specifically,

$g_i=\frac{G{M}_i}{(r_i/R{)}^2}.$

• Column F:   Labeled ${\mathrm{\Phi }}_i$ (for i between 1 and N); each row gives the value of the gravitational potential at the mid-point of a grid zone. Here, we start by specifying the value of the potential just (specifically, half a radial grid-zone) outside the surface of the configuration, where it should be, ${\mathrm{\Phi }}_N=-G{M}_{\mathrm{t}\mathrm{o}\mathrm{t}}/(1+\mathrm{\Delta }/2)$. Then, working from the surface, inward and, given that, $g=d\mathrm{\Phi }/dr$, we use the corresponding finite-difference representation of the radial derivative and set,

$\frac{{\mathrm{\Phi }}_i-{\mathrm{\Phi }}_{i-1}}{\mathrm{\Delta }}$	$=$	$g_i$
$\Rightarrow {\mathrm{\Phi }}_{i-1}$	$=$	${\mathrm{\Phi }}_i-g_i\mathrm{\Delta }.$

• Note that the value of the gravitational potential at the surface is not ${\mathrm{\Phi }}_N$ but, rather, must be ${\mathrm{\Phi }}_{\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}}=-G{M}_{\mathrm{t}\mathrm{o}\mathrm{t}}/R$.
• Furthermore, note that a lop-sided Taylor-series expansion about the center of the configuration provides the following good approximation to the gravitational potential at the center:   ${\mathrm{\Phi }}_c\approx (9{\mathrm{\Phi }}_1-{\mathrm{\Phi }}_2)/8$.
• Note as well that all of these numerically determined values of the gravitational potential can be checked against the known analytic expression for the radial profile of the potential in a uniform-density sphere.
• Column G:   Labeled $H_i$ (for i between 1 and N-1); each row provides the value of the fluid enthalpy at the center of a grid cell. Adopting the convention that the enthalpy is zero at the surface of the configuration, and given that the enthalpy and the gravitational potential must sum to zero throughout the configuration, we have,

$H_i={\mathrm{\Phi }}_{\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}}-{\mathrm{\Phi }}_i.$

• At the center of the configuration, we have, $H_c={\mathrm{\Phi }}_{\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}}-{\mathrm{\Phi }}_c$.
• Column H:   Labeled $H_{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}$ (for i between 1 and N-1); each row provides the value of the fluid enthalpy, renormalized to the central value, specifically,

$[H_{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}{]}_i=\frac{H_i}{H_c}.$

The second spreadsheet within the workbook should be initially created by generating a copy of the first spreadsheet. Then:

• Column C:   Labeled ${\rho }_i$ (for i between 1 and N-1); generate a new, improved guess for the normalized mass-density at each grid zone based on the corresponding value of the normalized enthalpy from the previous spreadsheet/iteration. Specifically, given that the relationship between the density and enthalpy in a polytrope of index, $n$, is, $\rho \propto H^n$, we should set,

$\{\frac{{\rho }_i}{{\rho }_c}{\}}_{\mathrm{s}\mathrm{h}\mathrm{e}\mathrm{e}\mathrm{t}\mathrm{2}}=\{[H_{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}{]}_i^n{\}}_{\mathrm{s}\mathrm{h}\mathrm{e}\mathrm{e}\mathrm{t}\mathrm{1}}.$

Straight Numerical Integration

via Direct Numerical Integration

The above governing relation may be rewritten as,

$\xi {\theta }^{{\text{'}}^{\prime }}+2{\theta }^{\text{'}}$

$=$

$-\xi {\theta }^n.$

We'll adopt the following finite-difference approximations for the first and second derivatives on a grid of radial spacing, ${\mathrm{\Delta }}_{\xi }$:

${{\theta }_i}^{\prime }$

$\approx$

$\frac{{\theta }_{+}-{\theta }_{-}}{2{\mathrm{\Delta }}_{\xi }}$

and,

${{\theta }_i}^{″}$

$\approx$

$\frac{{\theta }_{+}-2{\theta }_i+{\theta }_{-}}{{\mathrm{\Delta }}_{\xi }^2}.$

Our finite-difference approximation of the governing equation is, then,

${\xi }_i\left[\frac{{\theta }_{+}-2{\theta }_i+{\theta }_{-}}{{\mathrm{\Delta }}_{\xi }^2}\right]+2\left[\frac{{\theta }_{+}-{\theta }_{-}}{2{\mathrm{\Delta }}_{\xi }}\right]$	$=$	$-{\xi }_i{\theta }_i^n$
$\Rightarrow {\xi }_i[{\theta }_{+}-2{\theta }_i+{\theta }_{-}]+{\mathrm{\Delta }}_{\xi }[{\theta }_{+}-{\theta }_{-}]$	$=$	$-{\mathrm{\Delta }}_{\xi }^2{\xi }_i{\theta }_i^n$
$\Rightarrow {\theta }_{+}$	$=$	$\frac{2{\xi }_i{\theta }_i+{\theta }_{-}({\mathrm{\Delta }}_{\xi }-{\xi }_i)-{\mathrm{\Delta }}_{\xi }^2{\xi }_i{\theta }_i^n}{{\mathrm{\Delta }}_{\xi }+{\xi }_i}.$

Now, for the first two steps away from the center — where, ${\theta }_i={\theta }_0=1$ and ${\xi }_i={\xi }_0=0$ — we will use the following power-series expansion (see, for example, eq. 62 from §5 in Chapter IV of [C67]) to determine the value of ${\theta }_i$:

${\theta }_1$

$=$

$1-\frac{{\mathrm{\Delta }}_{\xi }^2}{6}+\frac{n{\mathrm{\Delta }}_{\xi }^4}{120}-\frac{n}{378}(\frac{n}{5}-\frac{1}{8}){\mathrm{\Delta }}_{\xi }^6,$

and,

${\theta }_2$

$=$

$1-\frac{(2{\mathrm{\Delta }}_{\xi }{)}^2}{6}+\frac{n(2{\mathrm{\Delta }}_{\xi }{)}^4}{120}-\frac{n}{378}(\frac{n}{5}-\frac{1}{8})(2{\mathrm{\Delta }}_{\xi }{)}^6.$

Results

Tabulated Global Properties

Here, drawing from tables that have been previously published by other authors, we record numerically determined properties of polytropic models having a fairly wide range of polytropic indexes. First, we draw from Table 4 (p.96) of [C67] . To convert from his tabulated variables to our desired 3 structural form-factors, our normalized equilibrium radius (see earlier ASIDE), and the "virial" (drawn from a more general overview), note that for isolated polytropes,

${\mathfrak{𝔣}}_{\mathrm{M}}$	$=$	$(\frac{{\rho }_c}{\overline{\rho }}{)}^{-1}=[-\frac{3{\theta }^{\text{'}}}{\xi }{]}_{{\xi }_1},$
$(5-n){\mathfrak{𝔣}}_{\mathrm{W}}$	$=$	$5{\mathfrak{𝔣}}_{\mathrm{M}}^2,$
$(5-n){\mathfrak{𝔣}}_{\mathrm{A}}$	$=$	$\left[\right(\frac{4\pi }{3})W_n{]}^{-1}=3(n+1)(-{\theta }^{\text{'}}{)}_{{\xi }_1}^2,$
$x_{\mathrm{e}\mathrm{q}}\equiv \frac{R_{\mathrm{e}\mathrm{q}}}{R_{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}}$	$=$	$[\frac{4\pi }{(n+1{)}^n}{]}^{1/(n-3)}{\xi }_1(-{\xi }^2{\theta }^{\text{'}}{)}_{{\xi }_1}^{(1-n)/(n-3)},$
Virial	$=$	$(5-n)[\frac{b}{n}\cdot x_{\mathrm{e}\mathrm{q}}^{(n-3)/n}-\frac{a}{3}]$
	$=$	$(\frac{3}{4\pi }{)}^{1/n}\frac{(5-n){\mathfrak{𝔣}}_A}{{\mathfrak{𝔣}}_M^{(n+1)/n}}\cdot x_{\mathrm{e}\mathrm{q}}^{(n-3)/n}-1.$

From Table 4 of [C67]
Copied Directly from Table (1st 5 columns)     …     Implied Values of 3 Structural Form Factors, $x_{\mathrm{e}\mathrm{q}}$, and Virial (last 5 columns)
n xi_1 "mass" rho_c/rho_avg W_n f_M (5-n)f_W (5-n)f_A x_eq Virial 0 2.4494 4.8988 1 0.119366 1 5 2.00000 0.620335 --- 0.5 2.7528 3.7871 1.8361 0.26227 0.544632645 1.483123592 0.910254 0.831089 -0.19009 1 3.14159 3.14159 3.28987 0.392699 0.303963378 0.461968677 0.607927 1.253313 3.0E-06 1.5 3.65375 2.71406 5.99071 0.77014 0.166925122 0.139319982 0.309857 2.357285 2.4E-06 2 4.35287 2.41105 11.40254 1.63818 0.087699758 0.038456238 0.145730 7.517481 1.4E-06 2.5 5.35528 2.18720 23.40646 3.90906 0.042723248 0.00912638 0.061072 186.3666 1.6E-08 3 6.89685 2.01824 54.1825 11.05066 0.018456144 0.001703146 0.0216035 --- --- 3.25 8.01894 1.94980 88.153 20.365 0.011343913 0.000643422 0.0117227 3.3265E-06 2.9E-06 3.5 9.53581 1.89056 152.884 40.9098 0.006540907 0.000213917 0.00583558 0.00166854 2.2E-06 4 14.97155 1.79723 622.408 247.558 0.001606663 1.29068E-05 0.00096435 0.051854 6.1E-06 4.5 31.83646 1.73780 6189.47 4922.125 0.000161565 1.30516E-07 4.8502E-05 0.284868 -5.9E-05 4.9 169.47 1.73205 9.348E+05 3.693E+06 1.06975E-06 5.7218E-12 6.4645D-08 2.129056 1.5E-04
The column labeled "mass" contains the tabulated quantity, $(-{\xi }^2{\theta }^{\text{'}}{)}_{{\xi }_1}$.

From Table 2.5.2 (p. 77) of Horedt (2004)   —   "Polytropic Spheres (N = 3)"
Copied Directly from Table (1st 3 columns)     …     Implied Values (last 7 columns)
n xi_1 theta' "mass" rho_c/rho_avg W_n f_M (5-n)f_A x_eq Virial 0 2.44948974 -8.164966E-01 4.898980 1.000000 1.193662E-01 1.000000E+00 2.000000E+00 6.2035049E-01 0.5 2.75269805 -4.999971E-01 3.788651 1.835143 2.122091E-01 5.449168E-01 1.124987E+00 8.3099030E-01 0.0E+00 1 3.14159265 -3.183099E-01 3.141593 3.289868 3.926990E-01 3.039636E-01 6.079272E-01 1.2533141E+00 0.0E+00 1.5 3.65375374 -2.033013E-01 2.714055 5.990704 7.701402E-01 1.669253E-01 3.099856E-01 2.3572860E+00 0.0E+00 2 4.35287460 -1.272487E-01 2.411047 1.140254E+01 1.638182E+00 8.769977E-02 1.457301E-01 7.5164793E+00 0.0E+00 2.5 5.35527546 -7.626491E-02 2.187199 2.340646E+01 3.909062E+00 4.272324E-02 6.107153E-02 1.8636634E+02 0.0E+00 3 6.89684862 -4.242976E-02 2.018236 5.418248E+01 1.105068E+01 1.845615E-02 2.160341E-02 3.5 9.53580534 -2.079098E-02 1.890557 1.528837E+02 4.090983E+01 6.540920E-03 5.835575E-03 1.6685566E-03 0.0E+00 4 1.49715463E+01 -8.018079E-03 1.797230 6.224079E+02 2.475594E+02 1.606664E-03 9.643439E-04 5.1854394E-02 0.0E+00 4.5 3.18364632E+01 -1.714549E-03 1.737799 6.189473E+03 4.921842E+03 1.615646E-04 4.850469E-05 2.8486849E-01 0.0E+00 4.99 1.75818915E+03 -5.598955E-07 1.730765 1.046736E+09 4.237887E+10 9.553503E-10 5.633289E-12 2.3460204E+01 2.2E-15
The column labeled "mass" contains the tabulated quantity, $(-{\xi }^2{\theta }^{\text{'}}{)}_{{\xi }_1}$.

Plotted Structural Profiles

Using the just-described numerical techniques, we have solved the polytropic Lane-Emden equation on a 200-zone, uniform grid for a variety of values of the polytropic index. In each case we have recorded how the dimensionless enthalpy, ${\theta }_n(\xi )$, and its first radial derivative, ${\theta }_n^{\text{'}}(\xi )\equiv d{\theta }_n/d\xi$, vary with $\xi$, from the center of the polytropic configuration to its surface. For the record, these tabulated results reside in the following DropBox files:

• n = 2.5: (10 SCF iterations)   WorkFolder/Wiki edits/HSCF/n25.xlsx
• n = 3: (19 SCF iterations)   WorkFolder/Wiki edits/HSCF/n300.xlsx
• n = 3.005: (15 SCF iterations)   WorkFolder/Wiki edits/HSCF/n3005.xlsx
• n = 3.05: (15 SCF iterations)   WorkFolder/Wiki edits/HSCF/n305.xlsx
• n = 3.5: (18 SCF iterations)   WorkFolder/Wiki edits/HSCF/n25.xlsx
• n = 6: (direct integration)   WorkFolder/Wiki edits/EmbeddedPolytropes/N6.xlsx

For each of these models, as indicated (n = 2.5, 3, 3.005, 3.05, 3.5, 6), Figure 4 illustrates how the normalized mass, $M/M_{\mathrm{S}\mathrm{W}\mathrm{S}}$, varies with the normalized radius, $R/R_{\mathrm{S}\mathrm{W}\mathrm{S}}$, where the definition of these two functions,

$\frac{M}{M_{\mathrm{S}\mathrm{W}\mathrm{S}}}$	$\equiv$	$\left(\frac{n^3}{4\pi }{)}^{1/2}{\theta }_n^{(n-3)/2}{\xi }^2\right|\frac{d{\theta }_n}{d\xi }|,$
$\frac{R}{R_{\mathrm{S}\mathrm{W}\mathrm{S}}}$	$\equiv$	$(\frac{n}{4\pi }{)}^{1/2}\xi {\theta }_n^{(n-1)/2},$

has been drawn from an accompanying discussion of pressure-truncated polytropic configurations. In four of the Figure 4 panels, we have compared the profile of our numerically determined polytropic function (curve defined by $\sim 200$ small, black circular markers) to results (7 - 9 larger, blue circular markers) taken from Table 2.5.1 of Horedt (2004) — see, specifically the segment of his table on pp. 74 - 75 that applies to polytropic spheres — in an effort to demonstrate that our numerically determined solutions are accurate.

Figure 4:   Numerically Determined Solutions to the Polytropic Lane-Emden Equation
↲	↲	↲
↲	↲	↲
Data examples from Table 2.5.1 (pp. 74 - 75) of Horedt (2004): $n$ $\xi$ ${\theta }_n$ $\frac{d{\theta }_n}{d\xi }$ $\frac{R}{R_{\mathrm{S}\mathrm{W}\mathrm{S}}}$ $\frac{M}{M_{\mathrm{S}\mathrm{W}\mathrm{S}}}$ 2.5 4.00000 0.1376807 - 0.1340534 0.4032551 3.926310 3.0 5.00000 0.1108198 - 0.08012604 0.2707342 2.936234 3.5 5.00000 0.1786843 - 0.07362030 0.3065541 2.210326 6.0 5.00000 0.3973243 - 0.05113662 0.3437981 1.327430

Emden's (1907) Tabulated Data

From Table 13 (p. 84) of Emden (1907)   —   "Global Properties"
n xi_1 - theta' 2nd deriv. "mass" rho_c/rho_avg 0 2.4494 0.81647 -0.33333 4.8988 1 0.5 2.7528 0.49975 0.36309 3.7871 1.8361 1 3.14159 0.31831 0.20264 3.14159 3.2899 1.5 3.6571 0.20316 0.11355 2.7176 6.0003 2 4.3518 0.12729 0.06262 2.4107 11.396 2.5 5.4172 0.075 0.02795 2.201 24.076 3 6.9011 0.04231 0.01282 2.015 54.36 4 14.999 0.00803 0.00107 1.8064 623.4 4.5 32.14 0.00168 0.000104 1.7354 6377.7 4.9 169.47 6.04E-05 4.208E-07 1.73554 9.485E+05 5 infinity 0 0 sqrt(3) infinity
The column labeled "mass" contains the tabulated quantity, $(-{\xi }^2{\theta }^{\text{'}}{)}_{{\xi }_1}$.

Horedt's (1986) Tabulated Data

G. P. Horedt (1986), Astrophysics and Space Science, Vol. 126, Issue 2, pp. 357 - 408: Seven-digit tables of Lane-Emden functions

In Table I we present seven digit numerical solutions of the Lane-Emden equation for the plane-parallel (N = 1), cylindrical (N = 2), and spherical (N = 3) case for polytropic indices of $n=-10,-5,-4,-3,-2,-1.5,-1.01,-0.9,-0.5,0,0.5,1,1.5,2,3,4,5,6,10,20,\pm \mathrm{\infty }$, supplemented by $n=2.5,3.5,4.5,$ and $4.99$ for the spherical case. In Table II some finite boundary values of polytropic slabs, cylinders, and spheres are summarized. For polytropic spheres (N = 3) we have also quoted boundary values near the minimum of the dimensionless mass $-{\xi }_1^2{\theta }_1$, occurring at n ≈ 4.823 (Z. F. Seidov and R. Kh. Kuzakhmedov, 1978).

Focusing specifically on the spherically symmetric (N = 3) configurations, we list here the page number(s) on which the table associated with each individual polytropic index can be found in Horedt (1986).

Spherical (N = 3) Configurations
n	page(s)	${\xi }_1$
- 10	386 → 387	∞
- 5	387 → 388	∞
- 4	388 → 389	∞
- 3	389	∞
- 2	390	∞
- 1.5	390 → 391	∞
- 1.01	391 → 392	∞
- 0.9	392 → 393	2.05040073E+00
- 0.5	393	2.20858842E+00
0	393 → 394	√6 = 2.44948974E+00
0.5	394	2.75269805E+00
1	394 → 395	π = 3.14159265E+00
1.5	395	3.65375374E+00
2	395 → 396	4.35287460E+00
2.5	396 → 397	5.35527546E+00
3	397 → 398	6.89684862E+00
3.5	398 → 399	9.53580534E+00
4	399	1.49715463E+01
4.5	399 → 400	3.18364632E+01
4.99	400 → 401	1.75818915+03
5	401 → 402	∞
6	402 → 403	∞
10	403 → 404	∞
20	404 → 405	∞
$\pm \mathrm{\infty }$	405 → 406	∞