Part I: Stellar Physics | Chapter 4

The Hertzsprung-Russell Diagram

The most important diagram in stellar astrophysics: Stefan-Boltzmann law, mass-luminosity relation, spectral classification, the magnitude system, and isochrone fitting

Overview

The Hertzsprung-Russell (HR) diagram, independently developed by Ejnar Hertzsprung (1911) and Henry Norris Russell (1913), plots stellar luminosity against effective temperature (or equivalently, absolute magnitude against spectral type or color). It is the single most important tool for understanding stellar populations, testing theories of stellar structure and evolution, and determining the ages of star clusters.

Stars do not occupy random positions on the HR diagram. Instead, they cluster along well-defined sequences: the main sequence, the red giant branch, the horizontal branch, the asymptotic giant branch, and the white dwarf cooling sequence. Each region corresponds to a specific phase of stellar evolution.

Figure 1. The Hertzsprung-Russell diagram plotting luminosity (in solar units, logarithmic) versus effective temperature (decreasing to the right). The diagonal main sequence band runs from hot, luminous O-stars to cool, faint M-dwarfs. Key regions include red giants (upper right), supergiants (top), and white dwarfs (lower left). Dashed lines show constant stellar radius. Positions of the Sun, Sirius, Betelgeuse, and Proxima Centauri are marked.

1. The Stefan-Boltzmann Law and Lines of Constant Radius

The luminosity of a star is related to its radius and effective temperature through the Stefan-Boltzmann law, which governs the thermal radiation from a blackbody.

1.1 Derivation of the Stefan-Boltzmann Law

A blackbody of temperature $T$ emits radiation with a spectral energy density given by the Planck function. The total energy density of blackbody radiation is obtained by integrating the Planck spectrum over all frequencies:

$$u = \int_0^\infty \frac{8\pi h\nu^3}{c^3} \frac{1}{e^{h\nu/k_BT} - 1} \, d\nu = \frac{8\pi^5 k_B^4}{15 c^3 h^3} T^4 = a T^4$$

where $a = 4\sigma/c$ is the radiation constant and $\sigma = 5.67 \times 10^{-8}$ W m$^{-2}$ K$^{-4}$is the Stefan-Boltzmann constant. The flux emitted per unit area of a blackbody surface is:

$$F = \sigma T^4$$

The factor of $c/4$ relating flux to energy density arises from averaging over the hemisphere of emission directions and accounting for the cosine projection factor.

For a spherical star of radius $R$ radiating as a blackbody at effective temperature $T_{\text{eff}}$, the total luminosity is:

Stefan-Boltzmann Law for Stars

$$\boxed{L = 4\pi R^2 \sigma T_{\text{eff}}^4}$$

This defines the effective temperature: $T_{\text{eff}}$ is the temperature of a blackbody with the same radius and luminosity as the star. For the Sun: $T_{\text{eff},\odot} = 5778$ K.

1.2 Lines of Constant Radius on the HR Diagram

Taking the logarithm of the Stefan-Boltzmann law:

$$\log L = 2\log R + 4\log T_{\text{eff}} + \text{const}$$

In solar units:

$$\log\left(\frac{L}{L_\odot}\right) = 2\log\left(\frac{R}{R_\odot}\right) + 4\log\left(\frac{T_{\text{eff}}}{T_{\text{eff},\odot}}\right)$$

On the HR diagram (where $\log L$ is plotted against $\log T_{\text{eff}}$with temperature increasing to the left), lines of constant radius are straight lines with slope 4:

$$\log L = 4 \log T_{\text{eff}} + \text{const}(R)$$

These diagonal lines run from upper-left (hot, luminous) to lower-right (cool, faint). The main sequence spans roughly $0.1{-}10 R_\odot$, while red giants reach$100{-}1000 R_\odot$ and white dwarfs are at $\sim 0.01 R_\odot$ (Earth-sized).

2. The Mass-Luminosity Relation

One of the most important empirical results in stellar astrophysics is that main-sequence stars obey a tight mass-luminosity relation. We derive this from the equations of stellar structure using homology arguments.

2.1 Homologous Stellar Models

Two stars are homologous if one can be obtained from the other by a simple rescaling of mass and radius: $r = (R/R_0) r_0$ for corresponding points. Under this assumption, all dimensionless variables (fractional mass, fractional radius, etc.) have the same profiles.

From the equation of hydrostatic equilibrium, $dP/dr = -Gm\rho/r^2$, we can extract scaling relations. Dimensionally:

$$\frac{P}{R} \sim \frac{G M}{R^2} \frac{M}{R^3} \implies P \sim \frac{G M^2}{R^4}$$

2.2 Derivation for Radiative Stars

From the ideal gas law: $P = \rho k_B T / (\mu m_H)$, and using $\rho \sim M/R^3$:

$$\frac{GM^2}{R^4} \sim \frac{M}{R^3} \frac{k_B T}{\mu m_H} \implies T \sim \frac{G \mu m_H M}{k_B R}$$

From the radiative transport equation:

$$\frac{T}{R} \sim \frac{\kappa \rho L}{T^3 R^2} \implies L \sim \frac{T^4 R^3}{\kappa \rho R^2} = \frac{T^4 R}{\kappa (M/R^3)}$$

Substituting the temperature scaling $T \propto M/R$:

$$L \sim \frac{(M/R)^4 R}{\kappa M / R^3} = \frac{M^4 R^3}{\kappa M R^4} \cdot R = \frac{M^3}{\kappa}$$

For electron-scattering opacity ($\kappa = \text{const}$), this gives:

Mass-Luminosity Relation (Homology)

$$\boxed{L \propto M^3}$$

For Thomson scattering opacity (constant $\kappa$), the exponent is 3. Including Kramers opacity ($\kappa \propto \rho T^{-3.5}$) steepens this to approximately $L \propto M^{3.5}$ for intermediate-mass main-sequence stars. Very massive stars approach $L \propto M$ as they near the Eddington limit.

2.3 Empirical Mass-Luminosity Relation

Observationally, the mass-luminosity relation for main-sequence stars is well described by:

$$\frac{L}{L_\odot} \approx \begin{cases} 0.23\left(\frac{M}{M_\odot}\right)^{2.3} & M < 0.43 M_\odot \\ \left(\frac{M}{M_\odot}\right)^{4.0} & 0.43 M_\odot < M < 2 M_\odot \\ 1.5\left(\frac{M}{M_\odot}\right)^{3.5} & 2 M_\odot < M < 20 M_\odot \\ 3200\left(\frac{M}{M_\odot}\right) & M > 55 M_\odot \end{cases}$$

The breaking of the relation into different regimes reflects the changing opacity, equation of state, and energy transport mechanism with stellar mass.

3. Spectral Classification and the Saha Equation

Stars are classified into spectral types OBAFGKM (and the cooler L, T, Y types for brown dwarfs) based on the absorption lines in their spectra. The sequence is one of decreasing temperature, and the appearance of spectral lines is governed by ionization equilibrium.

3.1 The Saha Equation

The degree of ionization of an atom in a stellar atmosphere is determined by the balance between ionization (by photon absorption or electron collision) and recombination. In thermal equilibrium, this balance is described by the Saha equation, derived from statistical mechanics.

Consider the ionization reaction $X_r \rightleftharpoons X_{r+1} + e^-$, where$X_r$ is an atom in ionization state $r$. The number densities satisfy:

Saha Equation

$$\boxed{\frac{n_{r+1} n_e}{n_r} = \frac{2 U_{r+1}}{U_r} \left(\frac{2\pi m_e k_BT}{h^2}\right)^{3/2} \exp\left(-\frac{\chi_r}{k_BT}\right)}$$

where $U_r$ and $U_{r+1}$ are the partition functions,$\chi_r$ is the ionization potential from stage $r$ to $r+1$,$n_e$ is the electron density, and the factor of 2 accounts for the two spin states of the free electron.

The derivation proceeds from the chemical potential equality in thermal equilibrium. The chemical potential of a species with number density $n$ and partition function $Z$ in a non-degenerate gas is:

$$\mu = -k_BT \ln\left(\frac{Z}{n}\right) + \varepsilon_0$$

where $\varepsilon_0$ is the ground-state energy. The equilibrium condition $\mu_r = \mu_{r+1} + \mu_e$ (ionization equilibrium) directly yields the Saha equation after substituting the translational partition function for the electron and noting that $\varepsilon_0^{(r+1)} + \varepsilon_0^{(e)} - \varepsilon_0^{(r)} = \chi_r$.

3.2 Application to Hydrogen

The visibility of hydrogen Balmer lines (the strongest optical absorption features) requires hydrogen atoms in the $n = 2$ excited state. The fraction of atoms in this state depends on:

Excitation (Boltzmann factor): The population of the $n=2$ level relative to $n=1$ scales as $\exp(-10.2 \text{ eV}/k_BT)$. This increases with temperature.
Ionization (Saha factor): As temperature increases further, hydrogen becomes ionized ($\chi = 13.6$ eV) and the neutral fraction drops.

The product of these two effects creates a peak in Balmer line strength at$T \approx 9500$ K (spectral type A0). Below this temperature, too few atoms are in the $n=2$ state; above it, too few neutral atoms remain.

3.3 The Spectral Sequence

Type	T_eff (K)	Color	Key Features
O	30,000–50,000	Blue	He II lines, weak H
B	10,000–30,000	Blue-white	He I, H lines strengthen
A	7,500–10,000	White	Strongest H Balmer lines
F	6,000–7,500	Yellow-white	Ca II, Fe I appear
G	5,200–6,000	Yellow	Strong Ca II H&K, Fe I
K	3,700–5,200	Orange	Strong metals, CH, CN bands
M	2,400–3,700	Red	TiO molecular bands

4. The Magnitude System

The astronomical magnitude system, rooted in Hipparchus's 2nd-century classification, provides a logarithmic measure of stellar brightness.

4.1 Apparent and Absolute Magnitude

The apparent magnitude $m$ measures the observed brightness of a star. A difference of 5 magnitudes corresponds to a factor of 100 in flux:

$$m_1 - m_2 = -2.5 \log_{10}\left(\frac{F_1}{F_2}\right)$$

The absolute magnitude $M$ is defined as the apparent magnitude a star would have at a distance of 10 parsecs. Since flux scales as $F \propto d^{-2}$:

$$m - M = -2.5 \log_{10}\left(\frac{F(d)}{F(10\text{pc})}\right) = -2.5 \log_{10}\left(\frac{10^2}{d^2}\right)$$

Distance Modulus

$$\boxed{m - M = 5\log_{10}\left(\frac{d}{10\text{ pc}}\right) = 5\log_{10}(d) - 5}$$

This is the distance modulus, one of the most fundamental relations in observational astronomy. For the Sun at 10 pc: $M_V = +4.83$. For Sirius at 2.64 pc:$m_V = -1.46$, $M_V = +1.42$.

4.2 Bolometric Corrections

The bolometric magnitude $M_{\text{bol}}$ measures the total luminosity integrated over all wavelengths. The bolometric correction$\text{BC}$ relates it to the visual magnitude:

$$M_{\text{bol}} = M_V + \text{BC}$$

The bolometric correction is always negative (convention: $\text{BC} \leq 0$) and is smallest ($\sim -0.07$) for F-type stars where most of the radiation falls in the visual band. For very hot (O-type) or very cool (M-type) stars, the bolometric correction can be several magnitudes because most of the radiation is in the UV or IR.

The relation between bolometric magnitude and luminosity is:

$$M_{\text{bol}} - M_{\text{bol},\odot} = -2.5 \log_{10}\left(\frac{L}{L_\odot}\right)$$

with $M_{\text{bol},\odot} = +4.74$.

4.3 Color Index

The color index is the difference in magnitudes measured through two different filters. The most common is $B - V$ (blue minus visual):

$$B - V = m_B - m_V = -2.5 \log_{10}\left(\frac{F_B}{F_V}\right) + \text{const}$$

Hotter stars have bluer spectra ($B - V < 0$), while cooler stars are redder ($B - V > 1$). The Sun has $B - V = +0.65$. The color index is a proxy for effective temperature, with the advantage that it does not require knowing the distance to the star.

5. Color-Magnitude Diagrams and Isochrone Fitting

A color-magnitude diagram (CMD) plots apparent magnitude versus color index for stars in a cluster. Since all stars in a cluster are at essentially the same distance and have the same age and initial composition, the CMD provides a powerful tool for determining cluster ages.

5.1 Isochrones

An isochrone (from Greek “equal time”) is a line on the HR diagram connecting the positions of stars that all have the same age but different masses. Isochrones are computed from stellar evolution models by evolving stars of different masses to the same age and plotting their $(\log T_{\text{eff}}, \log L)$ positions.

The key feature used for age determination is the main-sequence turnoff point (MSTO): the brightest and hottest point on the main sequence. Stars above the turnoff have exhausted their core hydrogen and evolved off the main sequence. The turnoff mass $M_{\text{TO}}$ satisfies:

$$t_{\text{age}} = t_{\text{MS}}(M_{\text{TO}}) \approx 10^{10} \left(\frac{M_{\text{TO}}}{M_\odot}\right)^{-2.5} \text{ years}$$

Inverting this gives the turnoff mass for a cluster of known age:

Main Sequence Turnoff Mass

$$\boxed{M_{\text{TO}} = M_\odot \left(\frac{t_{\text{age}}}{10^{10} \text{ yr}}\right)^{-0.4}}$$

For the Pleiades ($\sim 100$ Myr): $M_{\text{TO}} \approx 6 M_\odot$ (B-type stars). For 47 Tucanae ($\sim 12$ Gyr): $M_{\text{TO}} \approx 0.9 M_\odot$ (G-type stars).

5.2 Age Dating from CMDs

The procedure for isochrone fitting involves:

Observe the CMD of a star cluster (apparent magnitude vs color)
Apply a distance modulus to convert to absolute magnitudes
Correct for interstellar reddening ($E(B-V)$) and extinction ($A_V$)
Compare the observed CMD with a grid of theoretical isochrones at different ages
The best-fitting isochrone gives the cluster age

Additional age-sensitive features include the luminosity of the red giant branch tip, the morphology of the horizontal branch, and the luminosity function of the white dwarf cooling sequence.

5.3 Metallicity Effects

The chemical composition (metallicity) affects the isochrones through opacity: metal-rich stars have higher opacity, are more opaque, and thus appear cooler and fainter at a given mass. The main-sequence shifts redward with increasing metallicity. This degeneracy between age and metallicity is a major challenge in stellar population studies, partly broken by using multiple CMD features simultaneously.

6. Numerical Simulation: HR Diagram and Isochrones

We construct a comprehensive HR diagram showing the main sequence, giants, and white dwarfs, and compute isochrones for clusters of different ages.

HR Diagram with Main Sequence, Giants, White Dwarfs, and Isochrones

Python

script.py335 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ============================================================
# Stellar property functions
# ============================================================

def ms_luminosity(M):
    """Main-sequence luminosity for mass M (solar units)."""
    if hasattr(M, '__len__'):
        L = np.zeros_like(M, dtype=float)
        for i, m in enumerate(M):
            L[i] = ms_luminosity(m)
        return L
    if M < 0.43:
        return 0.23 * M**2.3
    elif M < 2.0:
        return M**4.0
    elif M < 20.0:
        return 1.5 * M**3.5
    else:
        return 3200.0 * M

def ms_temperature(M):
    """Main-sequence effective temperature for mass M (solar units)."""
    if hasattr(M, '__len__'):
        T = np.zeros_like(M, dtype=float)
        for i, m in enumerate(M):
            T[i] = ms_temperature(m)
        return T
    if M < 0.7:
        return 5780.0 * M**0.5
    elif M < 2.0:
        return 5780.0 * M**0.55
    else:
        return 5780.0 * M**0.5

def ms_lifetime(M):
    """Main-sequence lifetime in years."""
    L = ms_luminosity(M)
    return 1e10 * M / L  # proportional to fuel/luminosity

def bv_from_teff(T):
    """Approximate B-V color from effective temperature."""
    # Ballesteros formula (approximate)
    if hasattr(T, '__len__'):
        return np.array([bv_from_teff(t) for t in T])
    if T < 1000:
        return 2.0
    bv = -3.684 * np.log10(T) + 14.551
    return max(-0.4, min(2.0, bv))

def compute_isochrone(age_yr, n_points=200):
    """Compute an isochrone for a given age.
    Returns arrays of (log_Teff, log_L) for stars from low mass up."""
    # Turnoff mass
    M_to_arr = np.logspace(np.log10(0.1), np.log10(100), 500)
    t_ms_arr = np.array([ms_lifetime(m) for m in M_to_arr])

# Find turnoff mass
    M_to = 0.1
    for i in range(len(M_to_arr)):
        if t_ms_arr[i] < age_yr:
            M_to = M_to_arr[max(0, i-1)]
            break

log_T = []
    log_L = []

# Main sequence stars (below turnoff)
    ms_masses = np.logspace(np.log10(0.1), np.log10(max(M_to, 0.11)), n_points // 2)
    for m in ms_masses:
        L = ms_luminosity(m)
        T = ms_temperature(m)
        # Slight evolution along MS based on fraction of lifetime spent
        frac = min(age_yr / max(ms_lifetime(m), 1), 0.99)
        L *= (1.0 + 0.5 * frac)
        T *= (1.0 - 0.03 * frac)
        log_T.append(np.log10(T))
        log_L.append(np.log10(L))

# Post-MS evolution: subgiant, RGB, HB
    if M_to > 0.5:
        L_to = ms_luminosity(M_to) * 1.5
        T_to = ms_temperature(M_to) * 0.97

# Subgiant branch
        n_sg = n_points // 6
        for i in range(n_sg):
            frac = i / max(n_sg - 1, 1)
            T_sg = T_to * (1.0 - 0.35 * frac)
            L_sg = L_to * (1.0 + 3.0 * frac)
            log_T.append(np.log10(max(T_sg, 3000)))
            log_L.append(np.log10(max(L_sg, 0.001)))

# RGB
        n_rgb = n_points // 4
        T_rgb_base = 4200 if M_to > 2 else 3800
        L_rgb_base = L_to * 4
        L_rgb_tip = L_to * (30 if M_to > 2 else 80)
        for i in range(n_rgb):
            frac = i / max(n_rgb - 1, 1)
            T_rgb = T_rgb_base + 300 * (1.0 - frac)
            L_rgb = L_rgb_base * (L_rgb_tip / L_rgb_base)**frac
            log_T.append(np.log10(max(T_rgb, 3000)))
            log_L.append(np.log10(max(L_rgb, 0.001)))

return np.array(log_T), np.array(log_L), M_to

# ============================================================
# Create figure
# ============================================================
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
fig.patch.set_facecolor('#0a0a0a')

for ax in axes.flat:
    ax.set_facecolor('#111111')
    ax.tick_params(colors='white')
    ax.xaxis.label.set_color('white')
    ax.yaxis.label.set_color('white')
    ax.title.set_color('white')
    for spine in ax.spines.values():
        spine.set_color('#333')

# --- Plot 1: Full HR Diagram ---
ax1 = axes[0, 0]

# Main sequence
ms_m = np.logspace(np.log10(0.08), np.log10(100), 300)
ms_L = np.array([ms_luminosity(m) for m in ms_m])
ms_T = np.array([ms_temperature(m) for m in ms_m])
ms_logL = np.log10(ms_L)
ms_logT = np.log10(ms_T)
ax1.plot(ms_logT, ms_logL, color='#ffa500', linewidth=3, label='Main Sequence', zorder=3)

# Red giants (schematic)
rng = np.random.RandomState(42)
n_rg = 150
rg_logT = 3.5 + 0.15 * rng.randn(n_rg)
rg_logL = 1.0 + 1.5 * rng.rand(n_rg)
ax1.scatter(rg_logT, rg_logL, c='#ff4444', s=15, alpha=0.6, label='Red Giants', zorder=2)

# Supergiants
n_sg = 40
sg_logT = 3.55 + 0.3 * rng.randn(n_sg)
sg_logL = 3.5 + 1.5 * rng.rand(n_sg)
ax1.scatter(sg_logT, sg_logL, c='#ff6600', s=30, alpha=0.7, label='Supergiants', zorder=2, marker='*')

# White dwarfs
n_wd = 100
wd_logT = 3.8 + 0.5 * rng.rand(n_wd)
wd_logL = -2.0 - 2.5 * rng.rand(n_wd)
ax1.scatter(wd_logT, wd_logL, c='#aaaaff', s=10, alpha=0.6, label='White Dwarfs', zorder=2)

# Lines of constant radius
for log_R, label_R in [(-2, '0.01 R_sun'), (-1, '0.1 R_sun'), (0, '1 R_sun'),
                        (1, '10 R_sun'), (2, '100 R_sun'), (3, '1000 R_sun')]:
    T_line = np.logspace(3.3, 5.3, 100)
    # L = 4pi R^2 sigma T^4 => log L = 2*log_R + 4*log_T + const
    log_L_line = 2 * log_R + 4 * (np.log10(T_line) - np.log10(5778))
    ax1.plot(np.log10(T_line), log_L_line, 'w--', linewidth=0.5, alpha=0.3)
    # Label
    idx = len(T_line) // 3
    if -5 < log_L_line[idx] < 7:
        ax1.text(np.log10(T_line[idx]), log_L_line[idx] + 0.2, label_R,
                 color='#666', fontsize=7, rotation=-40)

ax1.set_xlabel('log T_eff (K)', fontsize=12)
ax1.set_ylabel('log L/L_sun', fontsize=12)
ax1.set_title('Hertzsprung-Russell Diagram', fontsize=14, fontweight='bold')
ax1.legend(fontsize=8, facecolor='#1a1a1a', edgecolor='#444', labelcolor='white', loc='lower left')
ax1.set_xlim(5.3, 3.3)
ax1.set_ylim(-5, 7)
ax1.grid(True, alpha=0.15)

# --- Plot 2: Isochrones for different ages ---
ax2 = axes[0, 1]
ages = [1e7, 5e7, 1e8, 5e8, 1e9, 5e9, 1.2e10]
age_labels = ['10 Myr', '50 Myr', '100 Myr', '500 Myr', '1 Gyr', '5 Gyr', '12 Gyr']
iso_colors = ['#ff2222', '#ff6633', '#ff9944', '#ffcc00', '#88dd44', '#44bbff', '#8844ff']

for i, (age, label) in enumerate(zip(ages, age_labels)):
    logT, logL, M_to = compute_isochrone(age)
    ax2.plot(logT, logL, color=iso_colors[i], linewidth=2, label=f'{label} (M_TO={M_to:.1f})', alpha=0.9)
    # Mark turnoff
    if len(logT) > 0:
        # Find the turnoff point (highest L on MS part)
        ms_end = len(logT) // 2
        if ms_end > 0:
            ax2.plot(logT[ms_end-1], logL[ms_end-1], 'o', color=iso_colors[i], markersize=5)

# ZAMS for reference
ax2.plot(ms_logT, ms_logL, 'w--', linewidth=1, alpha=0.4, label='ZAMS')

ax2.set_xlabel('log T_eff (K)', fontsize=12)
ax2.set_ylabel('log L/L_sun', fontsize=12)
ax2.set_title('Isochrones at Different Ages', fontsize=14, fontweight='bold')
ax2.legend(fontsize=7, facecolor='#1a1a1a', edgecolor='#444', labelcolor='white', loc='lower left')
ax2.set_xlim(4.8, 3.4)
ax2.set_ylim(-2, 6)
ax2.grid(True, alpha=0.15)

# --- Plot 3: Mass-Luminosity Relation ---
ax3 = axes[1, 0]
m_plot = np.logspace(np.log10(0.08), np.log10(100), 300)
L_plot = np.array([ms_luminosity(m) for m in m_plot])

ax3.loglog(m_plot, L_plot, color='#ffa500', linewidth=3)

# Power law references
m_ref = np.logspace(-0.5, 1.5, 100)
ax3.loglog(m_ref, m_ref**3.5, 'w--', linewidth=1, alpha=0.4, label=r'$L propto M^{3.5}$')
ax3.loglog(m_ref, m_ref**4.0, color='#ff4444', linestyle='--', linewidth=1, alpha=0.4, label=r'$L propto M^{4}$')

# Mark specific stars
star_data = [(1.0, 'Sun'), (0.12, 'Proxima Cen'), (2.0, 'Sirius A'),
             (17.0, 'Spica'), (50.0, 'Zeta Pup')]
for m, name in star_data:
    L = ms_luminosity(m)
    ax3.plot(m, L, 'o', color='#ffdd00', markersize=8, zorder=5)
    ax3.text(m * 1.15, L * 1.3, name, color='white', fontsize=8)

ax3.set_xlabel('Mass (M_sun)', fontsize=12)
ax3.set_ylabel('Luminosity (L_sun)', fontsize=12)
ax3.set_title('Mass-Luminosity Relation', fontsize=14, fontweight='bold')
ax3.legend(fontsize=10, facecolor='#1a1a1a', edgecolor='#444', labelcolor='white')
ax3.set_xlim(0.05, 200)
ax3.set_ylim(1e-5, 1e7)
ax3.grid(True, alpha=0.15)

# --- Plot 4: Saha equation - ionization fraction of hydrogen ---
ax4 = axes[1, 1]

T_range = np.linspace(3000, 30000, 500)
# Saha equation for hydrogen
# n_HII * n_e / n_HI = (2 * U_II / U_I) * (2pi me kT / h^2)^(3/2) * exp(-chi/kT)
# For H: U_I = 2 (ground state), U_II = 1, chi = 13.6 eV

k_B = 1.381e-23
m_e = 9.109e-31
h = 6.626e-34
eV = 1.602e-19
chi_H = 13.6 * eV

# Assume electron pressure P_e ~ 10 Pa (typical photosphere)
P_e_values = [1.0, 10.0, 100.0]
pe_colors = ['#44bbff', '#ffa500', '#ff4444']
pe_labels = ['P_e = 1 Pa', 'P_e = 10 Pa', 'P_e = 100 Pa']

for j, P_e in enumerate(P_e_values):
    x_ion = np.zeros_like(T_range)
    for i, T in enumerate(T_range):
        n_e = P_e / (k_B * T)
        saha_rhs = (2.0 / 2.0) * (2 * np.pi * m_e * k_B * T / h**2)**1.5 * np.exp(-chi_H / (k_B * T))
        # x^2 * n_total / (1-x) = saha_rhs  where n_e = x * n_total
        # For simplicity, n_total ~ n_e / x, solve quadratic
        # saha_rhs = x * n_e / (1 - x)
        # x = saha_rhs / (saha_rhs + n_e)
        x_ion[i] = saha_rhs / (saha_rhs + n_e)

ax4.plot(T_range, x_ion, color=pe_colors[j], linewidth=2.5, label=pe_labels[j])

# Mark spectral type boundaries
spec_temps = [(30000, 'O'), (10000, 'B'), (7500, 'A'), (6000, 'F'),
              (5200, 'G'), (3700, 'K'), (3000, 'M')]
for T_s, sp in spec_temps:
    if T_s <= 30000:
        ax4.axvline(x=T_s, color='#555', linestyle=':', linewidth=0.5)
        ax4.text(T_s, 1.05, sp, color='#aaa', fontsize=10, ha='center')

# Balmer line strength (schematic)
balmer = np.zeros_like(T_range)
for i, T in enumerate(T_range):
    # Fraction in n=2: Boltzmann * neutral fraction
    n_e = 10.0 / (k_B * T)
    saha_rhs = (2.0 / 2.0) * (2 * np.pi * m_e * k_B * T / h**2)**1.5 * np.exp(-chi_H / (k_B * T))
    x_neutral = n_e / (saha_rhs + n_e)
    # Boltzmann: n_2/n_1 = (g2/g1) * exp(-10.2 eV / kT)
    f_excited = 4.0 * np.exp(-10.2 * eV / (k_B * T))
    balmer[i] = x_neutral * f_excited

if np.max(balmer) > 0:
    balmer_norm = balmer / np.max(balmer)
else:
    balmer_norm = balmer
ax4.plot(T_range, balmer_norm, color='#ff69b4', linewidth=2, linestyle='--', label='Balmer strength (norm)')

ax4.set_xlabel('Temperature (K)', fontsize=12)
ax4.set_ylabel('Ionization fraction / Line strength', fontsize=12)
ax4.set_title('Hydrogen Ionization (Saha Equation)', fontsize=14, fontweight='bold')
ax4.legend(fontsize=9, facecolor='#1a1a1a', edgecolor='#444', labelcolor='white')
ax4.set_xlim(3000, 30000)
ax4.set_ylim(0, 1.15)
ax4.grid(True, alpha=0.15)

# Print results
print("=" * 65)
print("HERTZSPRUNG-RUSSELL DIAGRAM: KEY RESULTS")
print("=" * 65)

print("\nStefan-Boltzmann Law: L = 4*pi*R^2*sigma*T_eff^4")
print(f"  Sun: L = {3.828e26:.3e} W, R = {6.96e8:.2e} m, T_eff = 5778 K")

print("\nMass-Luminosity Relation (Main Sequence):")
print(f"{'Mass (Msun)':>12} {'L (Lsun)':>12} {'T_eff (K)':>12} {'Spec Type':>10}")
print("-" * 50)
spec_types = {0.1: 'M7', 0.3: 'M3', 0.5: 'M0', 0.8: 'K2', 1.0: 'G2',
              1.5: 'F5', 2.0: 'A5', 5.0: 'B5', 10.0: 'B1', 30.0: 'O7', 60.0: 'O4'}
for m, sp in sorted(spec_types.items()):
    L = ms_luminosity(m)
    T = ms_temperature(m)
    print(f"{m:>12.1f} {L:>12.1f} {T:>12.0f} {sp:>10}")

print("\nIsochrone turnoff masses:")
for age, label in zip(ages, age_labels):
    _, _, M_to = compute_isochrone(age)
    L_to = ms_luminosity(M_to)
    print(f"  Age = {label:>8s}: M_TO = {M_to:>5.1f} Msun, L_TO = {L_to:>8.0f} Lsun")

print("\nDistance modulus examples:")
objects = [('LMC', 50000), ('Andromeda (M31)', 778000), ('Virgo Cluster', 16.5e6)]
for name, d_pc in objects:
    dm = 5 * np.log10(d_pc / 10)
    print(f"  {name:>20s}: d = {d_pc:>12.0f} pc, m-M = {dm:>6.2f}")

print("\nSpectral type temperature scale:")
for T_s, sp in spec_temps:
    print(f"  {sp}: T_eff ~ {T_s:>6d} K")

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
print("\nPlot saved to output.png")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Chapter Summary

The HR diagram organizes stars by luminosity and temperature, revealing fundamental relationships. The Stefan-Boltzmann law $L = 4\pi R^2\sigma T_{\text{eff}}^4$ defines lines of constant radius, while the mass-luminosity relation $L \propto M^{3.5}$(derived from homology) explains why the main sequence is a mass sequence.

Spectral classification (OBAFGKM) is fundamentally a temperature sequence, with line strengths governed by the Saha ionization equation. The magnitude system and distance modulus $m - M = 5\log(d/10\text{pc})$ provide the observational framework.

Isochrone fitting to color-magnitude diagrams of star clusters is the primary method for determining stellar ages, with the main-sequence turnoff point being the key age indicator. The turnoff mass relates directly to age through $M_{\text{TO}} \propto t_{\text{age}}^{-0.4}$.

← Stellar Evolution Part II: Compact Objects →