Lane-Emden polytropes: Lane, Ritter, and Emden’s gas spheres

Lane (1870) first treated a star as a self-gravitating polytropic gas sphere, Ritter (1878-1889) generalized the analysis to arbitrary polytropic index \(n\), and Emden (1907) systematized the resulting family of solutions to the dimensionless equation

\[\frac{1}{\xi^2}\frac{d}{d\xi}\left(\xi^2\frac{d\theta}{d\xi}\right) + \theta^n = 0, \qquad \theta(0)=1,\ \theta'(0)=0.\]

lane_emden() integrates exactly this equation; PolytropicStar builds a physical star (radius, mass) from one solution. This example compares the profile shape across several classic polytropic indices, checks two of them against their known closed forms, and builds one physical star.

import matplotlib.pyplot as plt
import numpy as np

from physicskit.astro.stellar_structure import PolytropicStar, lane_emden
from physicskit.astro.visualizers import plot_lane_emden

The Lane-Emden family, across classic polytropic indices

n=0 (uniform density), n=1 (has a closed-form sinc solution), n=1.5 (non-relativistic degenerate matter / fully convective stars), n=3 (Eddington’s “standard model” / relativistic degenerate matter, below), and n=4 (approaching the n=5 index at which the star has infinite radius). Larger n means more centrally concentrated mass.

indices = [0.0, 1.0, 1.5, 3.0, 4.0]
fig1, ax1 = plt.subplots(figsize=(6.5, 5))
for n in indices:
    xi, theta = lane_emden(n)
    ax1.plot(xi, theta, label=f"n={n} (surface at $\\xi_1$={xi[-1]:.4f})")
ax1.axhline(0.0, color="0.7", lw=0.8)
ax1.set_xlabel(r"$\xi$")
ax1.set_ylabel(r"$\theta(\xi)$")
ax1.set_title("The Lane-Emden family: larger n is more centrally concentrated")
ax1.legend(fontsize=8)
fig1.tight_layout()

Checking two indices against their exact closed forms

n=0 has the exact solution \(\theta=1-\xi^2/6\), with its first zero at \(\xi_1=\sqrt6\); n=1 has the exact solution \(\theta=\sin(\xi)/\xi\), with its first zero at \(\xi_1=\pi\).

xi0, theta0 = lane_emden(0.0)
xi1, theta1 = lane_emden(1.0)
print(f"n=0 surface: numeric xi_1={xi0[-1]:.6f}, exact sqrt(6)={np.sqrt(6):.6f}")
print(f"n=1 surface: numeric xi_1={xi1[-1]:.6f}, exact pi={np.pi:.6f}")
print(f"n=0 profile max error vs exact 1-xi^2/6: {np.max(np.abs(theta0 - (1 - xi0**2 / 6))):.2e}")
with np.errstate(divide="ignore", invalid="ignore"):
    exact_n1 = np.where(xi1 > 1e-8, np.sin(xi1) / xi1, 1.0)
print(f"n=1 profile max error vs exact sin(xi)/xi:  {np.max(np.abs(theta1 - exact_n1)):.2e}")

n=0 surface: numeric xi_1=2.449490, exact sqrt(6)=2.449490
n=1 surface: numeric xi_1=3.141593, exact pi=3.141593
n=0 profile max error vs exact 1-xi^2/6: 8.88e-16
n=1 profile max error vs exact sin(xi)/xi:  7.11e-11

One physical star

PolytropicStar turns a single Lane-Emden solution into a physical star, given a polytropic constant K and central density rho_c (here in arbitrary \(G=1\) units).

star = PolytropicStar(n=1.5, K=1.0, rho_c=1.0)
print(f"\nn=1.5 star: xi_1={star.xi1:.6f}, alpha={star.alpha:.6f}, radius={star.radius:.6f}, mass={star.mass:.6f}")

xi_star, theta_star = lane_emden(1.5)
fig2, ax2 = plot_lane_emden(xi_star, theta_star)
ax2.set_title(f"n=1.5 star (radius={star.radius:.3f}, mass={star.mass:.3f} in code units)")
fig2.tight_layout()

plt.show()

n=1.5 star: xi_1=3.653754, alpha=0.446031, radius=1.629688, mass=3.026381