Stellar structure: hydrostatic equilibrium, energy transport and nuclear burning

We cannot see inside a star, yet we know the temperature at the centre of the Sun to within a few per hundred. That knowledge comes from physics: a handful of equations that any star must obey, solved with the measured mass, radius, luminosity and composition as boundary conditions, and then tested against neutrinos and vibrations that do reach us from the core. This lesson builds those equations one at a time. They express four ideas: a star is held up by pressure against its own gravity, it is made of gas whose pressure depends on density and temperature, the energy released in its core must be carried out to its surface, and that energy comes from nuclear fusion.

Hydrostatic equilibrium

Lesson 1 showed that the Sun would collapse in about half an hour without pressure, so at every point inside it the pressure must support the weight of the layers above. Consider a thin shell of radius \(r\), thickness \(dr\) and density \(\rho\). Gravity pulls it inward with force per unit area \(G m(r) \rho\,dr / r^2\), where \(m(r)\) is the mass inside radius \(r\). The pressure below it must exceed the pressure above it by exactly this amount, so

These are the equations of hydrostatic equilibrium and mass continuity. They can be turned into estimates at once. Replacing the derivative by \(P_c / R\) and the density by the mean density \(\bar\rho = 3M/4\pi R^3\) gives a central pressure of order \(GM\bar\rho/R\), about \(GM^2/R^4\) up to numerical factors. For the Sun that is about \(10^{14}\) to \(10^{15}\) pascals, a billion atmospheres or more; the full solar model gives \(2.5 \times 10^{16}\) Pa, because the density rises steeply toward the centre.

The virial theorem and the central temperature

What provides this pressure? Stellar interiors are fully ionized gases, and to a good approximation they obey the ideal gas law \(P = \rho k T / \mu m_H\), where \(\mu\) is the mean mass per particle in units of the hydrogen mass, about 0.6 for ionized solar material because each hydrogen atom contributes a proton and an electron. Combining the gas law with hydrostatic equilibrium leads to one of the most useful results in astrophysics, the virial theorem. Multiply hydrostatic equilibrium by \(4\pi r^3\) and integrate over the star; integrating the pressure term by parts and using the gas law, one finds

where \(K\) is the total thermal kinetic energy of the particles, \(U\) the gravitational potential energy and \(\bar T\) a mean temperature. Setting \(3NkT \approx GM^2/R\) with \(N = M/\mu m_H\) gives \(kT \approx GM \mu m_H / 3R\). For the Sun this is about 5 million kelvins as a mass-weighted average, with the centre several times hotter. The virial theorem has a surprising consequence. The total energy is \(E = K + U = -K\): a star that loses energy by radiating gets hotter, not colder, because it contracts and half the gravitational energy released goes into heating the gas. Self-gravitating systems have negative heat capacity, which is why a contracting gas cloud heats up until fusion ignites.

import math
# Eddington's standard model: a polytrope P = K rho^(1 + 1/n) with n = 3. Solve the Lane-Emden equation
# theta'' + (2/xi) theta' + theta^n = 0, theta(0) = 1, theta'(0) = 0, by fourth-order Runge-Kutta.
def lane_emden(n, h=1e-4):
    xi, th, dth = 1e-6, 1.0 - 1e-12 / 6, -1e-6 / 3
    def f(x, y, z): return z, -(max(y, 0.0) ** n) - 2 * z / x
    while th > 0:
        k1 = f(xi, th, dth); k2 = f(xi + h/2, th + h/2*k1[0], dth + h/2*k1[1])
        k3 = f(xi + h/2, th + h/2*k2[0], dth + h/2*k2[1]); k4 = f(xi + h, th + h*k3[0], dth + h*k3[1])
        th_new = th + h/6*(k1[0] + 2*k2[0] + 2*k3[0] + k4[0]); dth_new = dth + h/6*(k1[1] + 2*k2[1] + 2*k3[1] + k4[1])
        if th_new <= 0:
            frac = th / (th - th_new); return xi + frac * h, dth + frac * (dth_new - dth)
        xi, th, dth = xi + h, th_new, dth_new
xi1, dth1 = lane_emden(3)
print(f"n = 3: surface at xi1 = {xi1:.4f}, -xi1^2 theta'(xi1) = {-xi1**2 * dth1:.4f}")
conc = xi1 / (3 * -dth1)
W = 1 / (4 * math.pi * 4 * dth1**2)          # P_c = W G M^2 / R^4 for n = 3
G, Msun, Rsun, k, mH = 6.67430e-11, 1.989e30, 6.957e8, 1.380649e-23, 1.6735575e-27
rho_bar = Msun / (4 / 3 * math.pi * Rsun**3)
rho_c = conc * rho_bar
P_c = W * G * Msun**2 / Rsun**4
mu = 0.62
T_c = mu * mH * P_c / (k * rho_c)
print(f"central density = {conc:.2f} x mean = {rho_c:.3g} kg/m^3")
print(f"central pressure = {P_c:.3g} Pa, central temperature = {T_c / 1e6:.1f} million K")
print("standard solar model: 1.5e5 kg/m^3, 2.5e16 Pa, 15.7 million K")

n = 3: surface at xi1 = 6.8968, -xi1^2 theta'(xi1) = 2.0182
central density = 54.18 x mean = 7.64e+04 kg/m^3
central pressure = 1.25e+16 Pa, central temperature = 12.3 million K
standard solar model: 1.5e5 kg/m^3, 2.5e16 Pa, 15.7 million K

Arthur Eddington's 1926 book The Internal Constitution of the Stars used exactly this kind of polytropic model, in which pressure depends on density as a power law. With no knowledge of nuclear physics it gives a central temperature for the Sun within about 20 per hundred of the modern value of 15.7 million kelvins. The real Sun is more centrally concentrated than the polytrope, partly because 4.6 billion years of fusion have turned much of the core's hydrogen into helium, so the model's central density and temperature come out too low, but it shows how much follows from equilibrium alone.

Radiation pressure and the Eddington limit

Light carries momentum, and in hot stars the pressure of radiation, \(P_\text{rad} = aT^4/3\) with \(a = 4\sigma/c\), becomes comparable to gas pressure. Radiation pushing outward on electrons through Thomson scattering, with cross-section \(\sigma_T = 6.65 \times 10^{-29}\ \mathrm{m^2}\), can balance gravity pulling on protons. Setting the outward force on an electron, \(\sigma_T L / 4\pi r^2 c\), equal to the gravitational force on a proton, \(GMm_p/r^2\), gives a maximum luminosity for a star of mass \(M\), the Eddington luminosity:

Since luminosity grows faster than mass, the most massive stars approach this limit, and above about 100 to 150 solar masses radiation pressure drives off their outer layers in powerful winds. The same limit governs how fast black holes can swallow gas and how bright quasars can be.

Getting the energy out

The energy released in the core must reach the surface, and it travels in one of three ways: radiation, convection or conduction. In most of the Sun it travels as radiation, but not in a straight line. A photon in the core travels only about a millimetre before being absorbed or scattered, and each time it sets off in a random direction. A random walk of \(N\) steps of length \(\ell\) goes a distance of only \(\ell\sqrt N\), so to cross the Sun's radius \(R\) a photon needs \(N = (R/\ell)^2\) steps and a total path of \(R^2/\ell\). With a mean free path of about a millimetre the energy takes over 100,000 years to diffuse out, while the light from the surface reaches the Earth in eight minutes. The rate of radiative flow is set by the temperature gradient and the opacity \(\kappa\), how strongly the gas absorbs:

When the required temperature gradient becomes too steep, the gas begins to churn. A blob of gas displaced upward expands adiabatically; if it ends up hotter and less dense than its new surroundings, it keeps rising, and the layer is unstable to convection. Karl Schwarzschild stated this criterion in 1906: convection sets in when the actual temperature gradient is steeper than the adiabatic gradient. Convection happens where the opacity is high, as in the cool outer layers of the Sun where hydrogen is partly ionized, or where the energy flux is very concentrated, as in the cores of massive stars. The Sun has a radiative interior out to about 0.71 of its radius and a convective zone above that, whose tops of rising cells are visible as granulation on the solar surface. Stars below about 0.35 solar masses are convective throughout, and stars above about 1.3 solar masses have convective cores and radiative envelopes.

The third mechanism, conduction by electrons, is unimportant in ordinary stars, where electrons collide too often, but it dominates in the degenerate matter of white dwarfs and the cores of red giants. With the four equations, of hydrostatic equilibrium, mass continuity, energy transport and energy generation, plus the gas law, the opacity and the nuclear reaction rates, a computer can integrate a model star from centre to surface. The Vogt–Russell theorem, a rule of thumb rather than a proof, states that a star's structure is determined by its mass and composition. That is why the main sequence is a line rather than a region of the HR diagram.

Why the mass–luminosity relation is steep

The structure equations explain the mass–luminosity relation of lesson 6 with a short scaling argument. From the virial theorem, \(T \propto M/R\). The radiative transport equation, written as scalings with \(dT/dr \sim T/R\), gives \(L \propto R^4 T^4/(\kappa \rho R) \cdot (1/R)\) up to constants, which with \(\rho \propto M/R^3\) becomes \(L \propto R^4 T^4/(\kappa M)\). Substituting \(T \propto M/R\), the radius cancels: \(L \propto M^3/\kappa\). With roughly constant opacity, as for electron scattering in hot massive stars, \(L \propto M^3\); with the Kramers opacity of cooler stars, which grows as density rises and temperature falls, the relation steepens toward \(M^{3.5}\) to \(M^5\). The radius drops out: luminosity is fixed by how fast radiation can leak out, and the nuclear furnace adjusts its rate to match.

Timescales

Three timescales describe a star. The dynamical timescale, the free-fall time of about half an hour for the Sun, is how fast it responds to a loss of balance. The Kelvin–Helmholtz timescale, \(t_\text{KH} \approx GM^2/RL\), is how long gravitational contraction could supply its luminosity: about 30 million years for the Sun. The nuclear timescale is how long fusion can supply it: fusing hydrogen to helium releases 0.7 per hundred of the rest mass, and a star burns about the inner tenth of its hydrogen on the main sequence, so \(t_\text{nuc} \approx 0.007 \times 0.1 \times Mc^2/L\), about 10 billion years for the Sun. The ordering, dynamical much shorter than thermal much shorter than nuclear, means that a main-sequence star is always in mechanical and thermal equilibrium, slowly changing its composition.

Nuclear fusion in stars

At 15 million kelvins the average proton has a kinetic energy of about 2 keV, far less than the roughly 1 MeV Coulomb barrier between two protons. Classically fusion would be impossible. As you saw in modern physics, protons tunnel through the barrier, with a probability that rises steeply with energy, while the number of protons with high energy falls steeply in the Maxwell–Boltzmann tail. Their product peaks in a narrow window, the Gamow peak, near 6 keV in the Sun's core. Fusion therefore proceeds slowly, which is why the Sun lasts billions of years.

In stars like the Sun hydrogen fuses mainly through the proton–proton chain, worked out by Hans Bethe and Charles Critchfield in 1938. Its main branch is

with net result \(4p \to {}^4\mathrm{He} + 2e^+ + 2\nu_e\) and 26.73 MeV released, about 2 per hundred of it carried away by neutrinos. The first step is extraordinarily slow, because it needs a weak interaction to turn a proton into a neutron while the two protons are briefly close: a given proton in the Sun's core waits about 10 billion years on average. In stars more massive than about 1.3 solar masses, with hotter cores, hydrogen fuses mainly through the CNO cycle, proposed by Bethe and independently by Carl Friedrich von Weizsäcker in 1938 and 1939, in which carbon, nitrogen and oxygen nuclei act as catalysts, capturing four protons in turn and releasing a helium nucleus. Because the CNO cycle must penetrate the larger Coulomb barriers of carbon and nitrogen, its rate rises much more steeply with temperature, roughly as \(T^{17}\) near 20 million kelvins compared with about \(T^4\) for the proton–proton chain. Bethe received the 1967 Nobel Prize in Physics for this work.

import math
# Energy budget of hydrogen fusion and the solar neutrino flux
MeV, c, Msun, Lsun, AU, yr = 1.602176634e-13, 2.99792458e8, 1.989e30, 3.828e26, 1.496e11, 3.156e7
Q = 26.73 * MeV                       # energy per helium nucleus made
m_p, m_He = 1.67262192e-27, 6.6446573e-27
eff = Q / (4 * m_p * c**2)
print(f"fraction of rest mass released: {eff:.4f}")
rate_He = Lsun / (Q * 0.98)           # about 2 per hundred of the energy leaves as neutrinos
print(f"helium nuclei made per second: {rate_He:.3g}; hydrogen consumed {4 * m_p * rate_He / 1e9:.0f} billion kg per second")
nu_flux = 2 * rate_He / (4 * math.pi * AU**2)
print(f"solar neutrinos at Earth: {nu_flux:.3g} per m^2 per second ({nu_flux / 1e4:.3g} per cm^2 per second)")
t_nuc = 0.1 * eff * Msun * c**2 / Lsun / yr
print(f"nuclear timescale (inner tenth of the mass burned): {t_nuc / 1e9:.1f} billion years")
# temperature sensitivity: rate ~ exp(-3 (E_G / 4kT)^(1/3)) for the Gamow peak
k = 8.617333e-11                      # MeV per K
def gamow(T, Z1, Z2, A_red):
    EG = (math.pi * Z1 * Z2 / 137.036)**2 * 2 * A_red * 938.272
    tau = 3 * (EG / (4 * k * T))**(1 / 3)
    return tau**2 * math.exp(-tau), tau / 3
for name, Z1, Z2, Ar in (("p + p", 1, 1, 0.5), ("p + C-12", 1, 6, 12 / 13), ("p + N-14", 1, 7, 14 / 15)):
    r1, nu = gamow(15e6, Z1, Z2, Ar)
    print(f"{name:9s}: rate rises as about T^{nu - 2 / 3:.1f} near 15 million K")

fraction of rest mass released: 0.0071
helium nuclei made per second: 9.12e+37; hydrogen consumed 610 billion kg per second
solar neutrinos at Earth: 6.49e+14 per m^2 per second (6.49e+10 per cm^2 per second)
nuclear timescale (inner tenth of the mass burned): 10.5 billion years
p + p    : rate rises as about T^3.9 near 15 million K
p + C-12 : rate rises as about T^17.8 near 15 million K
p + N-14 : rate rises as about T^19.9 near 15 million K

The Sun converts about 4 million tonnes of mass into energy each second, burning roughly 600 million tonnes of hydrogen. The temperature exponents come from the Gamow factor: the steeper Coulomb barrier for carbon and nitrogen makes their reaction rates far more sensitive to temperature. The proton–proton exponent here is for the tunnelling step alone; the full chain behaves as about \(T^4\). The CNO cycle is limited by its slowest step, proton capture on nitrogen-14, whose exponent is close to the commonly quoted 17 to 20.

Testing the model: neutrinos and vibrations

Two kinds of messenger let us check the solar model directly. Neutrinos escape from the core almost unhindered, about 65 billion of them crossing each square centimetre of the Earth every second. Raymond Davis began detecting them in 1968 with 615 tonnes of cleaning fluid in the Homestake gold mine in South Dakota, and found only about a third of the predicted number. This solar neutrino problem persisted for three decades until the Sudbury Neutrino Observatory, 2 km underground in a nickel mine near Sudbury, Ontario, using 1000 tonnes of heavy water, showed in 2001 and 2002 that the total number of neutrinos matches the solar model but that about two thirds have changed type on the way to Earth. Neutrino oscillation requires neutrinos to have mass, and Arthur McDonald of SNO shared the 2015 Nobel Prize in Physics with Takaaki Kajita for this discovery. The Borexino experiment in Italy detected neutrinos from the Sun's CNO cycle in 2020, confirming that it supplies about 1 per hundred of the Sun's energy.

The second check is helioseismology. The Sun rings with millions of sound waves, excited by convection, and they make its surface rise and fall by a few hundred metres in periods near five minutes. The frequencies of these oscillations depend on the speed of sound throughout the interior, and measuring them from the ground and from space has mapped the Sun's internal sound speed to better than 1 part in 1000, confirming the depth of the convection zone and the rotation of the interior. The standard solar model passes both tests, and it is the foundation on which the theory of stellar evolution in the next lesson is built.