Section 1.1: From Hydrodynamic Equilibrium to the Lane-Emden Equation

The foundation for modeling a star's structure lies in the equation of hydrostatic equilibrium, which describes the balance of forces within the star. In this context, the inward gravitational force is countered by the outward pressure force. The hydrostatic equilibrium equation is expressed as follows:

P = f(ρ, r, φ)

where P represents pressure, ρ denotes density, r is the radial distance from the star's center, and φ indicates the gravitational potential. The gravitational potential can be described using the Poisson equation, which for a spherically symmetric star is expressed as:

∇²φ = 4πGρ

where G is the gravitational constant.

At this point, we have two differential equations but three unknown functions. To resolve this, we must introduce an additional condition by selecting an equation of state, which serves as a thermodynamic relationship between pressure and density.

As our equation of state, we adopt the polytropic model, which is typically represented using the polytropic index n. Different values of n correspond to various physical scenarios in stars, spanning from entirely degenerate matter to ideal gases.

The Lane-Emden equation is derived by integrating the equations of hydrostatic equilibrium, the Poisson equation, and the polytropic equation of state. By substituting the preceding expressions and introducing dimensionless variables, we obtain:

d²ζ/dz² + (2/z) dζ/dz + ζ^n = 0

where ζ is a dimensionless radial coordinate.

Section 1.2: Numerical Solutions Using Python

Now, let's tackle the solution of the Lane-Emden equation. The initial conditions we have are:

ζ(0) = 1

dζ/dz(0) = 0

Regrettably, analytical solutions exist only for certain values of n that do not accurately represent real stars. For realistic n, we must rely on numerical solutions. We will utilize SciPy's ODE solver, solve_ivp. First, we will reformulate the second-order Lane-Emden equation into a system of two first-order equations:

from scipy.integrate import solve_ivp

def solve_lane_emden(n, z_eval=None):

def f(z, y):

zeta, dzeta_dz = y

return [dzeta_dz, -2/z * dzeta_dz - z**n]

def end_condition(z, y):

return y[0]

end_condition.terminal = True

z_span = (1.E-6, 100)

y0 = [1, 0]

sol = solve_ivp(f, z_span, y0, t_eval=z_eval, events=[end_condition])

return sol

The function solve_lane_emden accepts two parameters: the polytropic index n and an optional parameter z_eval, which indicates the points at which the solution should be evaluated. Inside the function, we define f(z, y), which represents the system of two first-order ordinary differential equations derived from the Lane-Emden equation.

Another nested function, end_condition(z, y), specifies the termination condition for integration. It is set to trigger when ζ reaches zero, which signifies the boundary of the star. This function is marked as terminal, halting the integration when this condition is met.

The solve_ivp function from the SciPy library performs the numerical integration of the equation system over a specified range of z, from a small positive value to 100, with initial conditions of y0 = [1, 0]. If z_eval is provided, the function evaluates the solution at these specified points; otherwise, it integrates until the end condition is reached.

To test the function in a Jupyter notebook for a polytropic index of 3/2:

sol = solve_lane_emden(n=1.5)

The video discusses polytropic models and the Lane-Emden equation, providing further insights into these concepts.

Final Steps: Analyzing Results

After the termination event, we can determine the radius of the star based on the "time"-value of the event in the solver:

z_max = sol.t_events[0][0]

Now, we can request more detailed output:

import numpy as np

sol = solve_lane_emden(n=1.5, z_eval=np.linspace(1.E-6, z_max, 100))

Repeating this for various values of the polytropic index n allows us to visualize the results:

import matplotlib.pyplot as plt

ns = (1, 3/2., 5/3.)

sol = {}

for n in ns:

z_max = solve_lane_emden(n).t_events[0][0]

sol[n] = solve_lane_emden(n, z_eval=np.linspace(1.E-6, z_max, 100))

plt.plot(sol[n].t, sol[n].y[0, :], label=f"n={n:.2f}")

plt.xlabel("z")

plt.ylabel(r"$zeta$")

plt.legend()

plt.grid()

Finally, we can compute the mass of the star with n=3/2 using the integral of the Lane-Emden equation ζ(z):

from scipy.interpolate import interp1d

from scipy.integrate import quad

z = sol[1.5].t

z_max = z[-1]

zeta = sol[1.5].y[0, :]

zeta_interp = interp1d(z, zeta, kind="cubic")

mass_integ, err = quad(lambda z: zeta_interp(z) * z**2, 0, z_max)

The computed mass, mass_integ, yields a dimensionless mass that can be scaled to physical units based on the specific stellar model.