Equations of Stellar Structure

A star sits still for ten billion years, and the reason is a standoff. At every depth the weight of everything above presses down while a pressure gradient pushes back, and the two match to a part in a billion — hydrostatic equilibrium. At the same time the energy that leaks out of the surface as light must be replenished exactly, shell by shell, by nuclear burning far below. The equations of stellar structure are the mathematical statement of these two balances, mechanical and thermal, made local: four coupled first-order ordinary differential equations for how mass, pressure, luminosity, and temperature run from the centre to the surface. Two of them — mass continuity and hydrostatic equilibrium — are pure mechanics; the other two — energy generation and energy transport — are pure thermodynamics. What closes the set is a layer of material physics from this chapter's prerequisite articles: the equation of state $P(\rho,T,\mu)$ , the Rosseland opacity $\kappa(\rho,T)$ , and the nuclear rate $\varepsilon(\rho,T)$ . Rewrite the independent variable from radius to enclosed mass and the star becomes a clean two-point boundary value problem on the fixed domain $0\le m\le M$ . The payoff is enormous: for a given mass and composition the solution is essentially unique — the Vogt–Russell theorem — which is why the main sequence is a one-parameter family and why a whole evolutionary track is just a sequence of static structure solutions computed as the composition slowly burns.

Key takeaways

Four coupled first-order ODEs — mass continuity, hydrostatic equilibrium, energy generation, and energy transport — describe the run of $m$ , $P$ , $L$ , and $T$ through a star.
They are not self-contained: the equation of state $P(\rho,T,\mu)$ , the opacity $\kappa(\rho,T)$ , the nuclear rate $\varepsilon(\rho,T)$ , and a convection criterion supply the missing closure.
Rewritten in the enclosed-mass (Lagrangian) coordinate they become a two-point boundary value problem on the fixed domain $0\le m\le M$ : regularity at the centre, photosphere at the surface.
The Vogt–Russell theorem: mass and composition alone (nearly) fix a star's radius, luminosity, and effective temperature — the physical basis of the main sequence.
A timescale hierarchy — dynamical $\ll$ thermal $\ll$ nuclear — justifies treating evolution as a sequence of static models; the rigorous bound $P_c>GM^2/8\pi R^4$ gives $P_c>4.5\times10^{13}$ Pa for the Sun.
Key idea: four first-order ODEs — continuity, hydrostatic equilibrium, energy generation, and energy transport — closed by the equation of state, opacity, and nuclear rates, and solved as a boundary value problem in the mass coordinate, completely determine a star's structure from its mass and composition.

The core insight

Four first-order ODEs — continuity, hydrostatic equilibrium, energy generation, and energy transport — closed by the equation of state, opacity, and nuclear reaction rates, and solved as a boundary value problem in the mass coordinate, completely determine a star's structure from its mass and composition.

Where this lives today

The four equations, plus modern microphysics, are the engine of every stellar-evolution code — most visibly the open-source instrument MESA (Modules for Experiments in Stellar Astrophysics), which solves the Lagrangian set on an adaptive mesh with current opacities, nuclear networks, and equations of state, from brown dwarfs to supernova progenitors. Their predictions are tested to extraordinary precision by helioseismology: the Sun's measured p-mode frequencies pin the base of the convective envelope to $0.713\pm0.001\,R_\odot$ and the interior sound speed to fractions of a percent, exactly where the standard solar model's radiative/convective transition sits — a direct confirmation of the transport physics in (4) and (5). Asteroseismology with Kepler and TESS now does the same for thousands of other stars, reading off their densities and ages from oscillation spectra. And the neutrino side of the energy equation has closed the loop from the other direction: the Borexino experiment measured the solar pp-chain neutrinos and, in 2020, detected CNO-cycle neutrinos for the first time, directly confirming the $\varepsilon$ term the equations assume. The open frontiers — convective boundaries, overshooting, rotation, mass loss — are all about improving the closure relations, not the four differential balances, which remain as Eddington left them.

The physical picture

Slice the star into thin concentric shells and do the bookkeeping one shell at a time. A shell cannot fall inward, so the pressure pushing up on its inner face must exceed the pressure on its outer face by exactly the weight of the shell — hydrostatic equilibrium. A shell cannot hoard or leak energy in steady state, so the luminosity entering plus whatever it generates must equal the luminosity leaving — energy generation. These two ledgers applied to every shell are two of the four equations; the other two say how much mass a shell of given thickness holds (continuity) and how steep a temperature gradient pushes the flux through it (transport). The first figure freezes this picture on a single shell.

Figure

Bookkeeping on one spherical shell of thickness dr and mass dm = 4πr²ρ dr: the outward pressure force balances the shell's weight, giving dP/dr = −Gmρ/r², while the luminosity in plus what the shell generates equals the luminosity out, giving dL/dr = 4πr²ρε

The temperature gradient is set by how hard it is to move the heat outward. If radiative diffusion can carry the flux against a gentle gradient, the layer stays radiative. But if the gradient radiation would demand grows steeper than the gradient a rising blob of gas follows adiabatically, the gas becomes buoyantly unstable, overturns, and convection takes over — flattening the profile toward the adiabatic value. The star always transports heat the cheaper way, and a single criterion (Schwarzschild, refined by Ledoux) decides which regime holds at each point. The other essential move is the change of variable from radius to enclosed mass: a parcel keeps its mass label even as the star pulsates or contracts, so the domain never moves and abundances stay pinned to their shells. The second figure shows the finished product — the four variables through a standard model of the Sun.

Figure

A standard solar model: the run of temperature, density, luminosity, and enclosed mass through the Sun versus fractional radius. Central anchor values are T_c = 1.57×10⁷ K, ρ_c = 1.5×10⁵ kg m⁻³, and P_c = 2.5×10¹⁶ Pa; the radiative core gives way to the convective envelope at 0.713 R_⊙, and the surface sits at T_eff = 5772 K

The core equations

The four differential relations are the two mechanical balances plus the two thermal ones. Start with mass continuity, which simply counts the mass in a shell of radius $r$ and thickness $\mathrm{d}r$ : its volume is $4\pi r^2\,\mathrm{d}r$ and its density is $\rho$ .

\frac{\mathrm{d}m}{\mathrm{d}r}=4\pi r^{2}\rho \qquad (1)

where: $m$ — mass enclosed within radius $r$ ; $r$ — radial coordinate measured from the centre; $\rho$ — local mass density. Reading it: the enclosed mass grows outward at a rate fixed entirely by the local density and the surface area of the shell. Hydrostatic equilibrium is Newton's second law for a shell whose acceleration is negligibly small: the net pressure force per unit volume balances the local gravity.

\frac{\mathrm{d}P}{\mathrm{d}r}=-\frac{G m \rho}{r^{2}} \qquad (2)

where: $P$ — total pressure (gas plus radiation); $G$ — Newton's gravitational constant; $m$ — enclosed mass; $\rho$ — local density. Reading it: pressure falls outward, and it falls fastest where gravity $Gm/r^2$ and density are large — deep in the interior. Energy generation is the luminosity ledger: the net power crossing a sphere grows outward at the rate the enclosed shells release energy, from nuclear burning $\varepsilon$ , minus neutrino losses, plus any gravitational release from contraction.

\frac{\mathrm{d}L_r}{\mathrm{d}r}=4\pi r^{2}\rho\left(\varepsilon-\varepsilon_\nu+\varepsilon_{\mathrm{grav}}\right) \qquad (3)

where: $L_r$ — luminosity crossing radius $r$ ; $\varepsilon$ — nuclear energy release per unit mass per unit time; $\varepsilon_\nu$ — specific power lost to neutrinos; $\varepsilon_{\mathrm{grav}}=-T\,\mathrm{d}s/\mathrm{d}t$ — specific power from gravitational contraction or expansion. Reading it: in strict thermal equilibrium ( $\varepsilon_{\mathrm{grav}}=0$ ) the luminosity climbs from zero at the centre to $L$ across the burning core and then stays flat. Energy transport, in its radiative form, is photon diffusion: the flux is proportional to the gradient of the radiation energy density $aT^4$ , throttled by the opacity.

\frac{\mathrm{d}T}{\mathrm{d}r}=-\frac{3}{4ac}\frac{\kappa\rho}{T^{3}}\frac{L_r}{4\pi r^{2}}\quad(\text{radiative}) \qquad (4)

where: $T$ — local temperature; $\kappa$ — Rosseland-mean opacity per unit mass; $a=4\sigma/c$ — radiation constant; $c$ — speed of light; $L_r$ — local luminosity. Reading it: carrying a large flux through an opaque, cool region demands a steep temperature gradient — and if the required gradient is too steep, radiation gives way to convection. Which regime holds is decided pointwise by the Schwarzschild criterion, comparing the gradient radiation would need against the adiabatic gradient of a displaced blob.

\nabla_{\mathrm{rad}}>\nabla_{\mathrm{ad}}\;\Longrightarrow\;\text{convection} \qquad (5)

where: $\nabla_{\mathrm{rad}}$ — logarithmic gradient $\mathrm{d}\ln T/\mathrm{d}\ln P$ that pure radiative diffusion would require; $\nabla_{\mathrm{ad}}$ — adiabatic gradient a rising element follows. Reading it: when radiation is overtaxed the layer overturns, and the transport equation switches from (4) to $\nabla\approx\nabla_{\mathrm{ad}}$ . The four ODEs contain more unknown fields than equations, so they are closed by constitutive relations from the prerequisite articles — the equation of state, opacity, and nuclear rate — all functions of the local composition.

P=P(\rho,T,\mu),\quad \kappa=\kappa(\rho,T,\mu),\quad \varepsilon=\varepsilon(\rho,T,\mu) \qquad (6)

where: $\mu$ — mean molecular weight, fixed by the mass fractions $X$ (hydrogen), $Y$ (helium), $Z$ (metals). Reading it: material physics, not the differential structure, makes the count balance and encodes the chemistry. Finally, the system is a two-point boundary value problem: three conditions are set at the centre and two at the surface.

r=0,\ L_r=0 \ \text{at}\ m=0;\qquad \rho\to0,\ T\to T_{\mathrm{eff}} \ \text{at}\ m=M \qquad (7)

where: $M$ — total stellar mass; $T_{\mathrm{eff}}$ — effective (photospheric) temperature, with $L=4\pi R^2\sigma T_{\mathrm{eff}}^4$ . Reading it: regularity forces $r$ and $L_r$ to vanish at the centre, while a photospheric relation closes the surface — solvable only for discrete structures.

Derivation AThe structure equations in the Lagrangian mass coordinate

Our goal is to recast the four Eulerian equations, written with radius $r$ as the independent variable, into the Lagrangian form used by essentially every stellar-evolution code, where the enclosed mass $m$ is independent and the outer boundary sits at the fixed total mass $M$ . We invert mass continuity to get $\mathrm{d}r/\mathrm{d}m$ , chain-rule each remaining equation through that factor, and package the temperature equation through the dimensionless gradient $\nabla$ so one line covers both regimes via the Schwarzschild switch.

Step 1. Invert continuity. Equation (1) says $\mathrm{d}m/\mathrm{d}r=4\pi r^2\rho$ , so its reciprocal gives the geometry equation in the new coordinate:

\frac{\mathrm{d}r}{\mathrm{d}m}=\frac{1}{4\pi r^{2}\rho}.

Why: this factor $\mathrm{d}r/\mathrm{d}m$ is the key that turns every $\mathrm{d}/\mathrm{d}r$ into a $\mathrm{d}/\mathrm{d}m$ .

Step 2. Chain-rule the pressure equation. Multiply (2) by $\mathrm{d}r/\mathrm{d}m$ :

\frac{\mathrm{d}P}{\mathrm{d}m}=\frac{\mathrm{d}P}{\mathrm{d}r}\frac{\mathrm{d}r}{\mathrm{d}m}=-\frac{Gm\rho}{r^{2}}\cdot\frac{1}{4\pi r^{2}\rho}=-\frac{Gm}{4\pi r^{4}}.

Why: the density cancels cleanly, leaving a right-hand side with no $\rho$ — a real numerical convenience, since $\rho$ is itself a derived quantity from the equation of state.

Step 3. Localise energy generation. Dividing (3) by $\mathrm{d}m/\mathrm{d}r$ strips off the geometric factor $4\pi r^2\rho$ :

\frac{\mathrm{d}L}{\mathrm{d}m}=\varepsilon-\varepsilon_\nu+\varepsilon_{\mathrm{grav}}.

Why: in mass coordinates the luminosity change per unit mass is just the net specific energy release — the ledger becomes purely local, with no reference to where the shell sits.

Step 4. Write transport through the gradient. Define $\nabla\equiv\mathrm{d}\ln T/\mathrm{d}\ln P$ ; then $\mathrm{d}T/\mathrm{d}m=(T/P)\,\nabla\,\mathrm{d}P/\mathrm{d}m$ , and using Step 2:

\frac{\mathrm{d}T}{\mathrm{d}m}=-\frac{Gm\,T}{4\pi r^{4}P}\,\nabla.

Why: one equation now covers both cases — $\nabla=\nabla_{\mathrm{rad}}$ where Schwarzschild is stable, $\nabla\approx\nabla_{\mathrm{ad}}$ where unstable — so transport enters through a single scalar. The four together, $\mathrm{d}r/\mathrm{d}m$ , $\mathrm{d}P/\mathrm{d}m$ , $\mathrm{d}L/\mathrm{d}m$ , $\mathrm{d}T/\mathrm{d}m$ , are the Lagrangian structure equations. Why: multiplying each Lagrangian equation by $\mathrm{d}m/\mathrm{d}r=4\pi r^2\rho$ returns the Eulerian set (1)–(4), confirming the two forms are fully equivalent.

Go deeper · Derivation B · 1 minThe radiative temperature gradient from photon diffusion

Derivation BThe radiative temperature gradient from photon diffusion

The goal is to derive the radiative transport equation (4) from first principles, treating the interior radiation field as a photon gas that diffuses down a temperature gradient with a mean free path set by the opacity. We write the diffusion flux in terms of the radiation energy density $aT^4$ , identify it with the local luminosity per unit area, and solve for the gradient — showing why high opacity or low temperature steepens it and eventually trips convection.

Step 1. Photon mean free path and diffusion coefficient. A photon travels $\ell=1/(\kappa\rho)$ before absorption or scattering, so the radiation diffuses with coefficient $D=\tfrac13 c\ell=c/(3\kappa\rho)$ . Fick's law applied to the radiation energy density $u=aT^4$ gives the radiative flux:

F=-D\,\frac{\mathrm{d}u}{\mathrm{d}r}=-\frac{c}{3\kappa\rho}\frac{\mathrm{d}(aT^{4})}{\mathrm{d}r}.

Why: deep in a star the medium is so opaque that radiation transport is genuinely diffusive — the field's anisotropy is tiny and the diffusion approximation is excellent.

Step 2. Differentiate the source term. Since $\mathrm{d}(aT^4)/\mathrm{d}r=4aT^3\,\mathrm{d}T/\mathrm{d}r$ ,

F=-\frac{4ac\,T^{3}}{3\kappa\rho}\frac{\mathrm{d}T}{\mathrm{d}r}.

Why: the $T^3$ factor makes the flux extremely sensitive to temperature — a small gradient at high $T$ already carries an enormous flux.

Step 3. Identify the flux with the luminosity. The luminosity crossing radius $r$ is spread over a sphere, $F=L_r/4\pi r^2$ . Setting the two expressions equal and solving for the gradient:

\frac{\mathrm{d}T}{\mathrm{d}r}=-\frac{3}{4ac}\frac{\kappa\rho}{T^{3}}\frac{L_r}{4\pi r^{2}}.

Why: this is exactly the boxed transport equation (4) — and it makes the physics transparent, since large $\kappa$ , large $\rho$ , large $L_r$ , or small $T$ all demand a steeper gradient, precisely the conditions that make $\nabla_{\mathrm{rad}}$ exceed $\nabla_{\mathrm{ad}}$ and switch on convection. *Why: this reproduces the boxed transport equation (4) exactly, with $a=4\sigma/c$ , confirming the diffusion derivation is consistent with the Rosseland-mean treatment.*

Go deeper · Derivation C · 1 minA rigorous lower bound on the central pressure

Derivation CA rigorous lower bound on the central pressure

The goal is to prove — using only the two mechanical equations, with no thermodynamics whatsoever — that the pressure at the centre of any star in hydrostatic equilibrium must exceed $GM^2/8\pi R^4$ . The trick is to add a cleverly chosen auxiliary term to the pressure so the combination becomes a total derivative in the mass coordinate with a definite sign, then evaluate that monotone combination at the two ends.

Step 1. Build a perfect derivative. Add $Gm^2/8\pi r^4$ to $P$ and differentiate with respect to $m$ , using $\mathrm{d}P/\mathrm{d}m=-Gm/4\pi r^4$ from the Lagrangian pressure equation:

\frac{\mathrm{d}}{\mathrm{d}m}\!\left(P+\frac{Gm^{2}}{8\pi r^{4}}\right)=-\frac{Gm}{4\pi r^{4}}+\frac{Gm}{4\pi r^{4}}-\frac{Gm^{2}}{2\pi r^{5}}\frac{\mathrm{d}r}{\mathrm{d}m}.

Why: the auxiliary term is engineered so its explicit piece cancels $\mathrm{d}P/\mathrm{d}m$ exactly, leaving only the term from the shrinking radius.

Step 2. Fix the sign. The surviving term uses $\mathrm{d}r/\mathrm{d}m=1/(4\pi r^2\rho)>0$ , so

\frac{\mathrm{d}}{\mathrm{d}m}\!\left(P+\frac{Gm^{2}}{8\pi r^{4}}\right)=-\frac{Gm^{2}}{2\pi r^{5}}\cdot\frac{1}{4\pi r^{2}\rho}<0.

Why: the combination decreases monotonically outward, so its value at the centre must exceed its value at the surface — an inequality with no free parameters.

Step 3. Evaluate at the two ends. At the centre $m=0$ , so the auxiliary term vanishes and the combination equals $P_c$ . At the surface $m=M$ , $r=R$ , and $P_{\mathrm{surf}}\approx0$ :

P_c > P_{\mathrm{surf}}+\frac{GM^{2}}{8\pi R^{4}}\;\ge\;\frac{GM^{2}}{8\pi R^{4}}.

Why: dropping the tiny positive surface pressure only weakens the bound, so the result is rigorous for every star in hydrostatic equilibrium, regardless of composition or equation of state. *Why: for the Sun this bound gives $P_c>4.5\times10^{13}$ Pa, comfortably below the modelled $P_c\approx2.5\times10^{16}$ Pa — the two are compared in detail in Worked example 1.*

Limiting cases

Case	Result	Meaning
Mechanical limit (polytrope)	$\dfrac{\mathrm{d}P}{\mathrm{d}r}=-\dfrac{Gm\rho}{r^{2}},\ P=K\rho^{1+1/n}$	With a polytropic closure the two mechanical equations decouple from the thermal pair into the Lane–Emden equation, solvable without opacity or nuclear physics.
Fully radiative envelope	$\dfrac{\mathrm{d}T}{\mathrm{d}r}=-\dfrac{3}{4ac}\dfrac{\kappa\rho}{T^{3}}\dfrac{L_r}{4\pi r^{2}}$	Everywhere Schwarzschild-stable, $\nabla_{\mathrm{rad}}<\nabla_{\mathrm{ad}}$ ; transport is pure photon diffusion.
Efficient convection	$\nabla\approx\nabla_{\mathrm{ad}}=\dfrac{P\delta}{T\rho c_P}$	In deep, dense zones convection is so efficient the actual gradient sits a hair above adiabatic; the transport equation is replaced by $\nabla=\nabla_{\mathrm{ad}}$ .
Homology scaling	$L\propto\mu^{4}M^{3}\kappa_0^{-1}$	For an ideal gas with constant opacity, dimensional analysis of the four equations reproduces the observed mass–luminosity slope $L\propto M^{3\text{–}4}$ .

Common trapThe four equations are self-contained

They are not. The four ODEs involve $m$ , $P$ , $L$ , $T$ and also $\rho$ , $\kappa$ , $\varepsilon$ , and the gradient $\nabla$ — strictly more unknowns than equations. The system closes only when you attach the constitutive relations $P(\rho,T,\mu)$ , $\kappa(\rho,T,\mu)$ , and $\varepsilon(\rho,T,\mu)$ plus the composition profile $\mu(m)$ . A subtler trap: the transport equation is not one fixed law. You must test the Schwarzschild (or, with composition gradients, Ledoux) criterion at every point to decide whether $\nabla=\nabla_{\mathrm{rad}}$ or $\nabla\approx\nabla_{\mathrm{ad}}$ applies — writing (4) everywhere gives a wrong, convectively unstable structure in a massive-star core or a cool-star envelope.

Worked examples

1 · The central pressure of the Sun

Use the rigorous bound of Derivation C to estimate the Sun's central pressure, and compare with a detailed solar model. Take $G=6.674\times10^{-11}\ \mathrm{N\,m^2\,kg^{-2}}$ , $M_\odot=1.989\times10^{30}$ kg, $R_\odot=6.96\times10^{8}$ m.

Step 1. State the bound: $P_c>GM_\odot^2/(8\pi R_\odot^4)$ , which needs no thermodynamics — only that the Sun is in hydrostatic equilibrium.

Step 2. Numerator: $GM_\odot^2=6.674\times10^{-11}\times(1.989\times10^{30})^2=\mathbf{2.64\times10^{50}}$ (SI).

Step 3. Denominator: $8\pi R_\odot^4=8\pi\times(6.96\times10^{8})^4=\mathbf{5.90\times10^{36}}$ (SI).

Step 4. Divide: $P_c>2.64\times10^{50}/5.90\times10^{36}=\mathbf{4.5\times10^{13}\ Pa}$ . A standard solar model gives $P_c\approx2.5\times10^{16}$ Pa, so the bound holds and undershoots the true value by roughly 560 — expected, since real stars are far more centrally concentrated than the uniform-density limit the crude bound compares against. The bound $P_c>4.5\times10^{13}$ Pa is rigorous; the modelled value is $2.5\times10^{16}$ Pa.

2 · The central temperature of the Sun

Combine hydrostatic equilibrium with the ideal-gas law to estimate the Sun's central temperature to order of magnitude, using mean molecular weight $\mu\approx0.6$ for ionized solar-composition gas.

Step 1. Dimensional balance. Hydrostatic equilibrium sets the pressure scale $P_c\sim GM_\odot\rho_c/R_\odot$ , and the ideal gas law gives $P_c=\rho_c kT_c/(\mu m_H)$ . Eliminating $\rho_c$ leaves the compact estimate

kT_c\sim\frac{\mu m_H\,GM_\odot}{R_\odot}\;\Longrightarrow\;T_c\sim\frac{\mu m_H\,GM_\odot}{k R_\odot}.

Step 2. Evaluate $GM_\odot/R_\odot=6.674\times10^{-11}\times1.989\times10^{30}/6.96\times10^{8}=\mathbf{1.91\times10^{11}\ m^2\,s^{-2}}$ .

Step 3. Insert $\mu m_H=0.6\times1.673\times10^{-27}=1.00\times10^{-27}$ kg and $k=1.381\times10^{-23}$ J K⁻¹:

T_c\sim\frac{1.00\times10^{-27}\times1.91\times10^{11}}{1.381\times10^{-23}}=\mathbf{1.4\times10^{7}\ K}.

The detailed solar model gives $T_c=1.57\times10^{7}$ K — a one-line estimate lands within about 10%. A few-line hydrostatic argument predicts $T_c\sim1.4\times10^{7}$ K, essentially the fifteen-million-kelvin core temperature required to fuse hydrogen.

Physics of the portrait

The equations were assembled over sixty years, but one figure turned them from a set of relations into a working theory of what stars are made of: Arthur Eddington, who insisted the interior of a star is governed by the balance of gravity, gas pressure, and radiation pressure, and computed the consequences before anyone knew where the energy came from.

Origins · 1 minArthur Eddington and how this was discovered

Arthur Eddington

1882–1944 · the theory of the internal constitution of the stars

Arthur Stanley Eddington was the astrophysicist who made the inside of a star calculable. Culminating in his 1926 masterwork The Internal Constitution of the Stars, he combined hydrostatic equilibrium with the then-new understanding of radiation pressure to argue that a star is held up by the combined pressure of gas and photons, its luminosity throttled by the opacity of its interior. From this he derived the mass–luminosity relation — that a star's brightness is fixed mainly by its mass — years before the nuclear energy source was identified, a prediction that matched the data so well it convinced astronomers the framework was right while the energy source was still a mystery (his retort to the doubters: they should "go and find a hotter place"). He introduced the Eddington standard model, an $n=3$ polytrope in which radiation pressure plays a fixed fractional role, and the Eddington luminosity, the limit at which radiation pressure would blow a star apart. His deeper legacy is the very idea that the structure of a star follows deductively from a handful of local physical laws, closed by material properties and solved as a boundary value problem — the programme every stellar-evolution code still runs today.

History

1870

— Lane models the Sun as a self-gravitating gas sphere and derives the first polytropic relation

1907

— Emden's Gaskugeln systematizes polytropes into the Lane–Emden equation

1926

— Eddington publishes The Internal Constitution of the Stars and the mass–luminosity relation; the Vogt–Russell theorem is stated

1939

— Bethe works out the proton–proton chain and CNO cycle, supplying the missing energy source

1959

— Henyey introduces the relaxation method, making full numerical solution of the structure equations routine

Going further

Texts. Carroll & Ostlie, An Introduction to Modern Astrophysics, 2nd ed., Ch. 10 (The Interiors of Stars) · Kippenhahn, Weigert & Weiss, Stellar Structure and Evolution, 2nd ed., Ch. 9–11 · Hansen, Kawaler & Trimble, Stellar Interiors, 2nd ed., Ch. 1–5 · Prialnik, An Introduction to the Theory of Stellar Structure and Evolution, 2nd ed., Ch. 5–6. Key numerical anchors: the central-pressure bound $P_c>GM^2/8\pi R^4$ evaluates to $4.5\times10^{13}$ Pa for the Sun using $G=6.674\times10^{-11}$ N m² kg⁻², $M_\odot=1.989\times10^{30}$ kg, $R_\odot=6.96\times10^{8}$ m; the central-temperature estimate $T_c\sim\mu m_H GM_\odot/(kR_\odot)=1.4\times10^{7}$ K uses $\mu=0.6$ ; and the Lagrangian set reduces to the Eulerian set on multiplication by $4\pi r^2\rho$ . The convection-zone base $0.713\,R_\odot$ is the helioseismic value; the solar $T_c$ , $\rho_c$ , and $P_c$ are standard-solar-model values.

Prerequisites. Hydrostatic Equilibrium in Stars · Mass Continuity Equation · Equation of State for Stellar Matter · Virial Theorem for Stars. Leads to. Energy Generation Equation · Radiative Temperature Gradient · Convective Energy Transport · Schwarzschild and Ledoux Criteria for Convection · Mixing-Length Theory.

Ask the tutor: Ask me to walk through how these four equations, with a real opacity and nuclear network, actually build the main sequence and predict a star's radius and luminosity from its mass alone.

Why are there four differential equations but the system still needs more input?

Count the unknowns. The four ODEs govern $m$ , $P$ , $L$ , and $T$ , but their right-hand sides also contain $\rho$ , $\kappa$ , $\varepsilon$ , and the gradient $\nabla$ — more fields than equations. The missing relations are supplied by material physics: the equation of state $P(\rho,T,\mu)$ ties pressure to density and temperature, the opacity $\kappa(\rho,T,\mu)$ sets radiative transport, the nuclear rate $\varepsilon(\rho,T,\mu)$ sets energy generation, and the convection criterion selects $\nabla$ . With those closures and a composition profile $\mu(m)$ the count balances and a unique structure exists. This is why an evolution code is mostly microphysics tables, not the four lines of differential structure.

Why rewrite the equations in the mass coordinate instead of radius?

Because matter carries its mass label with it. As a star pulsates, contracts, or swells into a giant, a parcel keeps its enclosed-mass value $m$ even though its radius $r$ changes wildly. In the Lagrangian formulation the domain is always the fixed interval $0\le m\le M$ , so the outer boundary never moves and the mesh need not chase an expanding surface. Just as importantly, abundances are properties of the matter, so tracking $X(m)$ , $Y(m)$ , $Z(m)$ as burning proceeds is natural in mass and awkward in radius. A bonus: the pressure equation loses its explicit $\rho$ , which is numerically cleaner.

Is the Vogt–Russell theorem actually a theorem?

Not in the strict mathematical sense — it is better called the Vogt–Russell conjecture. It says mass and composition profile essentially determine a star's equilibrium structure. In practice it holds remarkably well and underlies the main sequence as a one-parameter family, but counterexamples exist: some stars admit more than one equilibrium structure for the same mass and composition (parts of evolutionary tracks, certain shell-burning configurations). So treat it as a powerful organizing principle that is true generically but not universally — enough to justify computing evolution as a sequence of unique static models, while remembering that bifurcations can occur.

How does the code decide whether a layer is radiative or convective?

At every mesh point it computes two gradients. The first, $\nabla_{\mathrm{rad}}$ , is the logarithmic temperature gradient that pure radiative diffusion would need to carry the local flux — read straight off equation (4) and the opacity. The second, $\nabla_{\mathrm{ad}}$ , is the adiabatic gradient a displaced blob of gas would follow, read off the equation of state. If $\nabla_{\mathrm{rad}}>\nabla_{\mathrm{ad}}$ the layer is buoyantly unstable and convects (Schwarzschild); the transport equation then uses $\nabla\approx\nabla_{\mathrm{ad}}$ in deep zones or a mixing-length value near the surface. If a composition gradient is present, the sterner Ledoux criterion replaces Schwarzschild. This pointwise test is what carves a star into its radiative and convective zones.

Why can stellar evolution be computed as a sequence of static models?

Because of a huge separation of timescales. The Sun's dynamical (sound-crossing) time is under an hour, its thermal (Kelvin–Helmholtz) time about thirty million years, and its nuclear time about ten billion years. Mechanical equilibrium is re-established almost instantly and thermal equilibrium far faster than composition changes, so at any instant the star is in essentially perfect hydrostatic and thermal balance relative to the glacial drift of its chemistry. Evolution then proceeds by taking a static solution, advancing the composition by one nuclear timestep, and re-solving — a quasi-static sequence. The $\varepsilon_{\mathrm{grav}}$ term handles the rare phases (pre-main-sequence contraction, shell flashes) where thermal balance briefly breaks.

What makes the surface boundary condition tricky?

The centre is easy — regularity forces $r\to0$ and $L_r\to0$ as $m\to0$ . The surface is subtle because a real star has no sharp edge: density and temperature fall smoothly through a photosphere. You cannot simply demand $P=0$ and $T=0$ , because the transport equation is singular there and the outer layers carry information about opacity and radiative transfer. In practice one attaches a model atmosphere or a photospheric relation linking $P$ , $T$ , and $L$ (often via $\tau=2/3$ and $L=4\pi R^2\sigma T_{\mathrm{eff}}^4$ ) and integrates inward to match the interior. Because conditions are set at both ends, it is a genuine two-point boundary value problem, solved by shooting or, more commonly, by the Henyey relaxation method.