Stellar Interior Structure Model

A numerical simulation of stellar interior dynamics, solving the coupled system of stellar structure ODEs from a star's surface inward to its core.

Overview

This project models the physical state of a star (pressure, mass, luminosity, temperature, density, opacity, and energy generation) as a function of radial depth, under the following simplifying assumptions:

• Fixed composition throughout the star (constant mean molecular weight $\mu$)
• Spherical symmetry — all quantities depend only on radius $r$
• Two energy transport modes — radiation and convection (whichever dominates at a given depth)

Physics

Stellar Structure Equations

The four coupled first-order ODEs solved are:

Equation	Description
$\frac{dP}{dr} = -\frac{GM_r\rho}{r^2}$	Hydrostatic equilibrium
$\frac{dM_r}{dr} = 4\pi r^2\rho$	Mass continuity
$\frac{dL_r}{dr} = 4\pi r^2\rho\varepsilon$	Luminosity gradient
$\frac{dT}{dr} = \min\left(\nabla_\text{rad},, \nabla_\text{adia}\right)$	Temperature gradient

Derived Quantities

At each integration step, the following are computed from the current state $(P, T)$:

• Density (ideal gas law): $$\rho = \frac{P,\mu, m_H}{k_B T}$$

• Opacity (Kramers' law, bound-free absorption): $$\kappa = 0.32,\rho\left(\frac{T}{10^6}\right)^{-3.5}$$

• Energy generation (simplified CNO cycle): $$\varepsilon = \varepsilon_0,\rho\times 0.70\left(\frac{T}{T_0}\right)^{20}$$

Convection vs. Radiation

The temperature gradient switches between:

• Radiative: $\nabla_\text{rad} = \frac{3\kappa\rho}{16\pi ac}\frac{L_r}{r^2 T^3}$
• Adiabatic (convective): $\nabla_\text{adia} = \frac{1}{\gamma_\text{rat}}\frac{GM_r}{r^2}\frac{\mu m_H}{k_B}$

Convection dominates when $\nabla_\text{rad} > \nabla_\text{adia}$ (typically at the core and outer envelope).

Numerical Methods

Fixed-Step RK4

A standard 4th-order Runge-Kutta integrator. Shown to be insufficient for this stiff ODE system — large step sizes cause numerical blow-up.

Adaptive RK45 (Runge-Kutta-Fehlberg)

The primary integration method. Uses the Butcher tableau to compute both a 4th- and 5th-order estimate at each step, then adapts the step size $h$ based on the local truncation error:

$$h_\text{new} = 0.9, h_\text{old} \left(\frac{\varepsilon_\text{tol}}{\varepsilon_\text{current}}\right)^{1/4}$$

This dramatically reduces computation time and memory use versus fixed-step RK4, and avoids the instability caused by the system's stiffness.

The state vector is represented as a NumPy array $\mathbf{y} = [P, M_r, L_r, T]$, enabling vectorised matrix operations via the Butcher tableau for efficiency.

Units

All quantities are in CGS units:

Constant	Value	Unit
$G$	$6.673 \times 10^{-8}$	cm³ g⁻¹ s⁻²
$k_B$	$1.381 \times 10^{-16}$	erg K⁻¹
$c$	$2.998 \times 10^{10}$	cm s⁻¹
$R_\odot$	$6.960 \times 10^{10}$	cm
$M_\odot$	$1.989 \times 10^{33}$	g
$L_\odot$	$3.826 \times 10^{33}$	erg s⁻¹

Boundary Conditions (Surface)

Integration starts at the stellar surface $r = R_\odot$ and proceeds inward:

Variable	Value
Pressure $P$	$0.1$ dyn cm⁻²
Enclosed mass $M_r$	$8,M_\odot$
Luminosity $L_r$	$10{,}000,L_\odot$
Temperature $T$	$7500$ K
Mean molecular weight $\mu$	$1.28$

Output

The integrator returns radial profiles of:

• Pressure
• Enclosed mass
• Luminosity
• Temperature
• Density
• Opacity
• Energy generation rate

These are visualised using matplotlib subplots.

Requirements

• Python 3.12+
• numpy
• matplotlib