Perturbation Theory

When exact solutions to the Schrödinger equation are unavailable — as is the case for virtually every system beyond the hydrogen atom — perturbation theory provides a systematic, order-by-order framework for constructing approximate solutions. Starting from a solvable reference Hamiltonian, perturbation theory builds corrections to both energies and wavefunctions, yielding some of the most important quantitative predictions in quantum chemistry and physics.

1. Introduction: When Exact Solutions Don't Exist

The Schrödinger equation can be solved exactly for only a handful of systems: the particle in a box, the harmonic oscillator, the rigid rotor, and the hydrogen atom. For any system with two or more interacting electrons, or for a one-electron system in an external field, we must resort to approximation methods.

Perturbation theory addresses the following scenario: suppose we can write the full Hamiltonian as

$\hat{H} = \hat{H}^{(0)} + \lambda \hat{H}'$

where $\hat{H}^{(0)}$ is a "zeroth-order" Hamiltonian whose eigenstates and eigenvalues are known exactly, $\hat{H}'$ is a "perturbation" representing the additional interaction, and $\lambda$ is a dimensionless parameter that tracks the order of the correction. The key assumption is that $\hat{H}'$ is "small" compared to $\hat{H}^{(0)}$, so that the exact eigenstates and eigenvalues can be expanded as power series in $\lambda$.

The power of perturbation theory lies in its generality: it applies whenever the problem can be decomposed into a solvable part plus a small correction. This includes atoms in external fields (Stark and Zeeman effects), electron-electron repulsion in helium, anharmonic corrections to molecular vibrations, spin-orbit coupling, and the intermolecular interactions that give rise to van der Waals forces.

Central idea: If the perturbation is small, the first few terms of the expansion provide an excellent approximation, and systematic improvement is possible by going to higher order.

2. Non-Degenerate Rayleigh-Schrödinger Perturbation Theory

2.1 Setting Up the Expansion

We assume that the unperturbed Hamiltonian $\hat{H}^{(0)}$ has a complete set of orthonormal eigenstates $|\psi_n^{(0)}\rangle$ with eigenvalues $E_n^{(0)}$:

$\hat{H}^{(0)} |\psi_n^{(0)}\rangle = E_n^{(0)} |\psi_n^{(0)}\rangle$

and that the state of interest $|\psi_k^{(0)}\rangle$ is non-degenerate, meaning$E_k^{(0)} \neq E_n^{(0)}$ for all $n \neq k$. We expand both the exact eigenvalue and eigenstate as power series in $\lambda$:

$E_k = E_k^{(0)} + \lambda E_k^{(1)} + \lambda^2 E_k^{(2)} + \lambda^3 E_k^{(3)} + \cdots$

$|\psi_k\rangle = |\psi_k^{(0)}\rangle + \lambda |\psi_k^{(1)}\rangle + \lambda^2 |\psi_k^{(2)}\rangle + \cdots$

Here $E_k^{(1)}$ is the first-order energy correction, $E_k^{(2)}$ is the second-order correction, and so forth. We adopt the intermediate normalization convention:$\langle \psi_k^{(0)} | \psi_k \rangle = 1$, which implies$\langle \psi_k^{(0)} | \psi_k^{(j)} \rangle = 0$ for all $j \geq 1$.

2.2 Collecting Terms by Powers of λ

Substituting the expansions into the full Schrödinger equation$(\hat{H}^{(0)} + \lambda \hat{H}')|\psi_k\rangle = E_k |\psi_k\rangle$ and collecting terms order by order in $\lambda$:

Order $\lambda^0$: $\hat{H}^{(0)} |\psi_k^{(0)}\rangle = E_k^{(0)} |\psi_k^{(0)}\rangle$ (the unperturbed equation)

Order $\lambda^1$: $\hat{H}^{(0)} |\psi_k^{(1)}\rangle + \hat{H}' |\psi_k^{(0)}\rangle = E_k^{(0)} |\psi_k^{(1)}\rangle + E_k^{(1)} |\psi_k^{(0)}\rangle$

Order $\lambda^2$: $\hat{H}^{(0)} |\psi_k^{(2)}\rangle + \hat{H}' |\psi_k^{(1)}\rangle = E_k^{(0)} |\psi_k^{(2)}\rangle + E_k^{(1)} |\psi_k^{(1)}\rangle + E_k^{(2)} |\psi_k^{(0)}\rangle$

2.3 First-Order Energy Correction

From the first-order equation, project onto $\langle \psi_k^{(0)} |$ from the left:

$\langle \psi_k^{(0)} | \hat{H}^{(0)} | \psi_k^{(1)} \rangle + \langle \psi_k^{(0)} | \hat{H}' | \psi_k^{(0)} \rangle = E_k^{(0)} \langle \psi_k^{(0)} | \psi_k^{(1)} \rangle + E_k^{(1)}$

Since $\hat{H}^{(0)}$ is Hermitian,$\langle \psi_k^{(0)} | \hat{H}^{(0)} | \psi_k^{(1)} \rangle = E_k^{(0)} \langle \psi_k^{(0)} | \psi_k^{(1)} \rangle$. The terms involving $\langle \psi_k^{(0)} | \psi_k^{(1)} \rangle$ cancel on both sides (by intermediate normalization, this inner product is zero). We obtain the celebrated result:

$\boxed{E_k^{(1)} = \langle \psi_k^{(0)} | \hat{H}' | \psi_k^{(0)} \rangle}$

The first-order energy correction is the expectation value of the perturbation in the unperturbed state.

2.4 First-Order Wavefunction Correction

To find $|\psi_k^{(1)}\rangle$, we project the first-order equation onto$\langle \psi_n^{(0)} |$ with $n \neq k$. Since$\langle \psi_n^{(0)} | \hat{H}^{(0)} = E_n^{(0)} \langle \psi_n^{(0)} |$:

$E_n^{(0)} \langle \psi_n^{(0)} | \psi_k^{(1)} \rangle + \langle \psi_n^{(0)} | \hat{H}' | \psi_k^{(0)} \rangle = E_k^{(0)} \langle \psi_n^{(0)} | \psi_k^{(1)} \rangle$

Solving for the expansion coefficient $\langle \psi_n^{(0)} | \psi_k^{(1)} \rangle$:

$\langle \psi_n^{(0)} | \psi_k^{(1)} \rangle = \frac{\langle \psi_n^{(0)} | \hat{H}' | \psi_k^{(0)} \rangle}{E_k^{(0)} - E_n^{(0)}}$

Since $|\psi_k^{(1)}\rangle$ can be expanded in the complete basis of unperturbed states, and using intermediate normalization ($\langle \psi_k^{(0)} | \psi_k^{(1)} \rangle = 0$):

$\boxed{|\psi_k^{(1)}\rangle = \sum_{n \neq k} \frac{\langle \psi_n^{(0)} | \hat{H}' | \psi_k^{(0)} \rangle}{E_k^{(0)} - E_n^{(0)}} |\psi_n^{(0)}\rangle}$

Note: this requires $E_k^{(0)} \neq E_n^{(0)}$ — the non-degeneracy condition.

2.5 Second-Order Energy Correction

From the second-order equation, projecting onto $\langle \psi_k^{(0)} |$:

$\langle \psi_k^{(0)} | \hat{H}' | \psi_k^{(1)} \rangle = E_k^{(2)}$

Substituting the expression for $|\psi_k^{(1)}\rangle$ derived above:

$\boxed{E_k^{(2)} = \sum_{n \neq k} \frac{|\langle \psi_n^{(0)} | \hat{H}' | \psi_k^{(0)} \rangle|^2}{E_k^{(0)} - E_n^{(0)}}}$

The second-order energy requires knowledge of all unperturbed states, not just the state of interest.

Important properties of $E_k^{(2)}$:

• For the ground state ($k = 0$), every denominator$E_0^{(0)} - E_n^{(0)} < 0$ (since $E_0^{(0)}$ is the lowest energy), so $E_0^{(2)} \leq 0$. The second-order correction always lowers the ground-state energy.
• States that are close in energy to state $k$ contribute most strongly (small denominator).
• States not connected by the perturbation ($\langle \psi_n^{(0)} | \hat{H}' | \psi_k^{(0)} \rangle = 0$) do not contribute, regardless of their energy.

3. Degenerate Perturbation Theory

3.1 The Problem with Degeneracy

When two or more unperturbed states share the same energy, $E_k^{(0)} = E_l^{(0)}$, the denominator $E_k^{(0)} - E_l^{(0)}$ in the first-order wavefunction vanishes, and non-degenerate perturbation theory breaks down. The physical reason is clear: when states are degenerate, any linear combination of them is also an eigenstate of $\hat{H}^{(0)}$, and the perturbation must "choose" the correct linear combinations.

3.2 The Secular Equation

Suppose states $|\psi_1^{(0)}\rangle, |\psi_2^{(0)}\rangle, \ldots, |\psi_d^{(0)}\rangle$ are$d$-fold degenerate with energy $E^{(0)}$. We seek the correct zeroth-order states as linear combinations:

$|\phi\rangle = \sum_{i=1}^{d} c_i |\psi_i^{(0)}\rangle$

Substituting into the first-order equation and projecting onto each $\langle \psi_j^{(0)} |$within the degenerate subspace yields the secular equation:

$\sum_{i=1}^{d} \left( H'_{ji} - E^{(1)} \delta_{ji} \right) c_i = 0$

where $H'_{ji} = \langle \psi_j^{(0)} | \hat{H}' | \psi_i^{(0)} \rangle$. This is an eigenvalue problem: diagonalize the perturbation matrix within the degenerate subspace. The eigenvalues give the first-order corrections $E^{(1)}$, and the eigenvectors give the correct zeroth-order states.

$\boxed{\det \left( \mathbf{H}' - E^{(1)} \mathbf{I} \right) = 0}$

The secular determinant: a $d \times d$ eigenvalue problem within the degenerate subspace.

3.3 Procedure

Identify the degenerate subspace: all states with $E_i^{(0)} = E^{(0)}$.
Construct the $d \times d$ matrix of the perturbation $H'_{ij}$ within this subspace.
Diagonalize $\mathbf{H}'$ to find eigenvalues $E_\alpha^{(1)}$ and eigenvectors $|\phi_\alpha\rangle$.
The correct zeroth-order states are the eigenvectors $|\phi_\alpha\rangle$.
The first-order energies are $E_\alpha = E^{(0)} + \lambda E_\alpha^{(1)}$.
If the degeneracy is fully lifted, higher-order corrections follow from non-degenerate PT using the new basis.

4. The Stark Effect in Hydrogen (n = 2)

4.1 Setup

The Stark effect — the splitting of spectral lines in an external electric field — is the classic application of degenerate perturbation theory. Consider the $n = 2$ level of hydrogen, which is 4-fold degenerate (ignoring spin):

$|2,0,0\rangle, \quad |2,1,0\rangle, \quad |2,1,1\rangle, \quad |2,1,-1\rangle$

All four states have energy $E_2^{(0)} = -13.6/4 = -3.4$ eV.

The perturbation due to a uniform electric field $\mathcal{E}$ along the $z$-axis is:

$\hat{H}' = e\mathcal{E}z = e\mathcal{E}r\cos\theta$

4.2 Selection Rules and Matrix Elements

Since $z = r\cos\theta$ is odd under parity, the perturbation connects only states of opposite parity. Furthermore, $\cos\theta \propto Y_1^0$, so the selection rules are:

• $\Delta l = \pm 1$ (parity selection rule)
• $\Delta m = 0$ (since $\hat{H}'$ has no $\phi$-dependence)

The only nonzero matrix element within the $n = 2$ subspace is:

$\langle 2,0,0 | \hat{H}' | 2,1,0 \rangle = e\mathcal{E} \langle 2,0,0 | r\cos\theta | 2,1,0 \rangle = -3e\mathcal{E}a_0$

This matrix element is computed by evaluating the radial integral$\int_0^\infty R_{20}(r) \, r \, R_{21}(r) \, r^2 \, dr = -3\sqrt{6} \, a_0$ (using the explicit hydrogen radial wavefunctions) and the angular integral$\int Y_0^0 \cos\theta \, Y_1^0 \, d\Omega = 1/\sqrt{3}$.

4.3 The Perturbation Matrix and Diagonalization

In the basis $\{|2,0,0\rangle, |2,1,0\rangle, |2,1,1\rangle, |2,1,-1\rangle\}$, the perturbation matrix is:

$\mathbf{H}' = \begin{pmatrix} 0 & -3e\mathcal{E}a_0 & 0 & 0 \\ -3e\mathcal{E}a_0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}$

The secular equation $\det(\mathbf{H}' - E^{(1)}\mathbf{I}) = 0$ yields:

$(E^{(1)})^2 \left[ (E^{(1)})^2 - 9e^2\mathcal{E}^2 a_0^2 \right] = 0$

The four eigenvalues are:

$\boxed{E^{(1)} = +3e\mathcal{E}a_0, \quad 0, \quad 0, \quad -3e\mathcal{E}a_0}$

The linear Stark effect: splitting proportional to the electric field strength.

The correct zeroth-order states corresponding to the nonzero eigenvalues are:

$|\phi_{\pm}\rangle = \frac{1}{\sqrt{2}} \left( |2,0,0\rangle \mp |2,1,0\rangle \right)$

The states $|2,1,\pm1\rangle$ remain degenerate and unshifted to first order. The 4-fold degenerate level splits into 3 levels, with the middle level remaining doubly degenerate. This is the linear Stark effect, characteristic of hydrogen where states of different$l$ but the same $n$ are degenerate.

5. Ground State Polarizability

5.1 Second-Order Energy in an Electric Field

For the non-degenerate ground state of an atom in a uniform electric field, the first-order energy correction vanishes by parity ($\langle 0 | z | 0 \rangle = 0$ for states of definite parity). The leading effect comes from the second-order correction:

$E^{(2)} = e^2 \mathcal{E}^2 \sum_{n \neq 0} \frac{|\langle n | z | 0 \rangle|^2}{E_0 - E_n}$

This energy shift is proportional to $\mathcal{E}^2$, which defines the static electric polarizability $\alpha$ through the relation$E^{(2)} = -\frac{1}{2}\alpha \mathcal{E}^2$:

$\boxed{\alpha = 2e^2 \sum_{n \neq 0} \frac{|\langle n | z | 0 \rangle|^2}{E_n - E_0} = 2 \sum_{n \neq 0} \frac{|\langle n | \mu_z | 0 \rangle|^2}{E_n - E_0}}$

where $\mu_z = ez$ is the $z$-component of the dipole moment operator.

5.2 Closure Approximation for Hydrogen

For the hydrogen atom ground state, the exact sum over all excited states (including the continuum) is difficult to evaluate directly. The closure approximation replaces all energy denominators by an average excitation energy $\Delta E$:

$\alpha \approx \frac{2e^2}{\Delta E} \sum_{n \neq 0} |\langle n | z | 0 \rangle|^2 = \frac{2e^2}{\Delta E} \left( \langle 0 | z^2 | 0 \rangle - |\langle 0 | z | 0 \rangle|^2 \right)$

where we used completeness. For the hydrogen 1s state, $\langle 0 | z | 0 \rangle = 0$ by parity and$\langle 0 | z^2 | 0 \rangle = \langle 0 | r^2 \cos^2\theta | 0 \rangle = \frac{1}{3}\langle r^2 \rangle_{1s} = a_0^2$. Using $\Delta E \approx E_2 - E_1 = 10.2$ eV as a lower bound for the average excitation energy:

$\alpha \approx \frac{2e^2 a_0^2}{\Delta E} \approx 4 a_0^3 \cdot \frac{4\pi\epsilon_0}{e^2/(2a_0)} \cdot \frac{2 \cdot 13.6}{10.2} \approx 5.3 \, a_0^3 \cdot 4\pi\epsilon_0$

The exact result, obtained by Epstein (1926) using parabolic coordinates, is:

$\boxed{\alpha = \frac{9}{2} a_0^3 \cdot 4\pi\epsilon_0 = 4.5 \, a_0^3 \cdot 4\pi\epsilon_0}$

Exact polarizability of the hydrogen atom ground state.

The closure approximation overestimates the polarizability because $\Delta E$ was taken as the minimum excitation energy rather than a proper average. Nevertheless, it provides a useful estimate and correctly captures the $a_0^3$ scaling.

6. Applications

Molecular Polarizabilities

Second-order perturbation theory predicts how molecules respond to external electric fields. The polarizability tensor $\alpha_{ij}$ governs the induced dipole moment$\mu_i^{\text{ind}} = \sum_j \alpha_{ij} \mathcal{E}_j$ and determines refractive indices, dielectric constants, and Raman scattering intensities. Larger molecules with low-lying excited states tend to have larger polarizabilities.

Van der Waals Interactions

London dispersion forces arise from second-order perturbation theory applied to the instantaneous dipole-dipole interaction between two atoms or molecules. The resulting attractive potential scales as $V(R) \sim -C_6/R^6$, where $C_6$ depends on the polarizabilities and ionization energies of both species. These ubiquitous forces are responsible for the condensation of noble gases and much of the cohesion in molecular crystals.

Zeeman Effect

The interaction of an atom with an external magnetic field, $\hat{H}' = -\hat{\boldsymbol{\mu}} \cdot \mathbf{B}$, splits degenerate $m_l$ and $m_s$ sublevels. In the weak-field (anomalous) Zeeman effect, degenerate perturbation theory within the $J$ manifold yields energies $E_m = g_J \mu_B B m_J$, where $g_J$ is the Landé g-factor.

Spin-Orbit Coupling

The relativistic interaction between an electron's spin and its orbital angular momentum,$\hat{H}_{\text{SO}} = \xi(r)\hat{\mathbf{L}} \cdot \hat{\mathbf{S}}$, is treated as a perturbation to the non-relativistic Hamiltonian. This lifts the degeneracy of states with the same $L$ and $S$ but different $J$, producing the fine structure of atomic spectra. The splitting follows the Landé interval rule.

Møller-Plesset Perturbation Theory in Quantum Chemistry

In modern computational chemistry, Møller-Plesset (MP) perturbation theory uses the Hartree-Fock Hamiltonian as $\hat{H}^{(0)}$ and the difference between the exact electron-electron repulsion and the mean-field HF potential as $\hat{H}'$. The second-order correction (MP2) captures a large fraction of the electron correlation energy at a cost scaling as $\mathcal{O}(N^5)$ with system size, making it one of the most widely used post-Hartree-Fock methods for computing molecular energies, geometries, and properties.

7. Historical Context

1894

Lord Rayleigh develops perturbation methods for vibrating systems in classical mechanics, establishing the mathematical framework that would later be adapted for quantum theory. His approach of expanding solutions in power series of a small parameter laid the foundation for what would become Rayleigh-Schrödinger perturbation theory.

1913

Johannes Stark observes the splitting of hydrogen spectral lines in an electric field, providing the experimental phenomenon that would become one of perturbation theory's greatest triumphs when explained quantum mechanically.

1926

Erwin Schrödinger formulates the quantum mechanical version of perturbation theory in his foundational papers on wave mechanics. He derives the formulas for energy and wavefunction corrections and applies them to the Stark effect, reproducing Epstein's earlier old-quantum-theory results within the new framework.

1932

Brillouin and Wigner independently develop an alternative formulation of perturbation theory in which the exact energy appears in the denominators rather than the unperturbed energy. Brillouin-Wigner perturbation theory is not size-extensive but is conceptually simpler and converges differently from the Rayleigh-Schrödinger form.

1934

Christian Møller and Milton Plesset apply Rayleigh-Schrödinger perturbation theory to many-electron atoms and molecules using the Hartree-Fock solution as the reference state. Their formulation (Møller-Plesset perturbation theory, or MP theory) remains one of the most important methods in computational quantum chemistry, with the MP2 approximation becoming a standard tool for including electron correlation.

8. Python Simulation: Anharmonic Oscillator

The following simulation computes energy corrections for an anharmonic oscillator with potential$V(x) = \frac{1}{2}x^2 + \lambda x^4$ using perturbation theory (first and second order) and compares them against exact diagonalization in a truncated harmonic oscillator basis. The matrix elements of $x^4$ are computed exactly using the creation and annihilation operator representation. Observe how perturbation theory becomes progressively less accurate as $\lambda$ increases.

Anharmonic Oscillator: Perturbation Theory vs Exact Diagonalization

Python

Compares first-order and second-order perturbation theory energy corrections with exact eigenvalues for the anharmonic oscillator V = x^2/2 + lambda*x^4.

anharmonic_oscillator_pt.py100 lines

import numpy as np

# ──────────────────────────────────────────────────────────────
# Perturbation Theory vs Exact Diagonalization
# Anharmonic oscillator: V = 0.5*x^2 + lambda*x^4
# Compare PT energy corrections with exact eigenvalues
# ──────────────────────────────────────────────────────────────

N = 60  # basis size (harmonic oscillator eigenstates)

# Build matrix elements of x^4 in the HO basis
# x = sqrt(1/2) * (a + a^dag), so x^4 involves products of (a + a^dag)^4
# Matrix elements of a and a^dag:
#   <m|a|n> = sqrt(n) * delta(m, n-1)
#   <m|a^dag|n> = sqrt(n+1) * delta(m, n+1)

def build_x_matrix(N):
    """Build the position operator x in the HO basis."""
    x = np.zeros((N, N))
    for n in range(N - 1):
        x[n, n + 1] = np.sqrt((n + 1) / 2.0)
        x[n + 1, n] = np.sqrt((n + 1) / 2.0)
    return x

x_mat = build_x_matrix(N)
x2 = x_mat @ x_mat
x4 = x2 @ x2

# Unperturbed Hamiltonian (diagonal): E_n = n + 0.5
H0 = np.diag(np.arange(N) + 0.5)

# Sweep over lambda values
lambda_vals = np.linspace(0.0, 0.3, 50)

# Store results for ground state
E_exact_ground = []
E_pt1_ground = []
E_pt2_ground = []

# Store results for first excited state
E_exact_first = []
E_pt1_first = []
E_pt2_first = []

for lam in lambda_vals:
    # --- Exact diagonalization ---
    H = H0 + lam * x4
    eigvals = np.sort(np.linalg.eigvalsh(H))
    E_exact_ground.append(eigvals[0])
    E_exact_first.append(eigvals[1])

# --- Perturbation theory for ground state (n=0) ---
    # E0 = 0.5
    # First-order: E1 = <0|x^4|0>
    E0_0 = 0.5
    E1_0 = x4[0, 0]  # <0|x^4|0> = 3/4
    # Second-order: E2 = sum_{n!=0} |<n|x^4|0>|^2 / (E_0 - E_n)
    E2_0 = 0.0
    for n in range(1, N):
        matrix_el = x4[n, 0]
        if abs(matrix_el) > 1e-15:
            E2_0 += matrix_el**2 / (E0_0 - (n + 0.5))
    E_pt1_ground.append(E0_0 + lam * E1_0)
    E_pt2_ground.append(E0_0 + lam * E1_0 + lam**2 * E2_0)

# --- Perturbation theory for first excited state (n=1) ---
    E0_1 = 1.5
    E1_1 = x4[1, 1]  # <1|x^4|1> = 15/4
    E2_1 = 0.0
    for n in range(N):
        if n == 1:
            continue
        matrix_el = x4[n, 1]
        if abs(matrix_el) > 1e-15:
            E2_1 += matrix_el**2 / (E0_1 - (n + 0.5))
    E_pt1_first.append(E0_1 + lam * E1_1)
    E_pt2_first.append(E0_1 + lam * E1_1 + lam**2 * E2_1)

E_exact_ground = np.array(E_exact_ground)
E_exact_first = np.array(E_exact_first)
E_pt1_ground = np.array(E_pt1_ground)
E_pt2_ground = np.array(E_pt2_ground)
E_pt1_first = np.array(E_pt1_first)
E_pt2_first = np.array(E_pt2_first)

print("lambda        E0_exact    E0_PT1      E0_PT2      E1_exact    E1_PT1      E1_PT2")
for i in range(0, len(lambda_vals), 5):
    lam = lambda_vals[i]
    print(f"{lam:8.4f}      {E_exact_ground[i]:10.6f}  {E_pt1_ground[i]:10.6f}  {E_pt2_ground[i]:10.6f}  {E_exact_first[i]:10.6f}  {E_pt1_first[i]:10.6f}  {E_pt2_first[i]:10.6f}")

print()
print("--- Ground State Error Analysis ---")
print("lambda    |E_PT1 - E_exact|   |E_PT2 - E_exact|   PT2/PT1 improvement")
for i in range(0, len(lambda_vals), 5):
    lam = lambda_vals[i]
    err1 = abs(E_pt1_ground[i] - E_exact_ground[i])
    err2 = abs(E_pt2_ground[i] - E_exact_ground[i])
    ratio = err1 / err2 if err2 > 1e-15 else float('inf')
    print(f"{lam:8.4f}    {err1:14.8f}      {err2:14.8f}      {ratio:8.2f}x")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

9. Fortran Simulation: Stark Effect

This Fortran program computes the energy splittings of the hydrogen $n = 2$ states as a function of the external electric field strength. The 4×4 perturbation matrix is constructed using the analytically known matrix element $\langle 2,0,0 | eEz | 2,1,0 \rangle = -3ea_0\mathcal{E}$and diagonalized using Jacobi rotations. The output confirms the linear Stark splitting into three distinct energy levels.

Stark Effect: Hydrogen n=2 Energy Splittings

Fortran

Computes energy levels of the hydrogen n=2 states as a function of electric field strength, demonstrating the linear Stark effect via degenerate perturbation theory.

stark_effect_hydrogen.f90172 lines

program stark_effect_hydrogen
    ! ──────────────────────────────────────────────────────────────
    ! Stark Effect: Energy splittings of hydrogen n=2 states
    ! as a function of external electric field strength.
    !
    ! The n=2 level of hydrogen has 4 degenerate states:
    !   |2,0,0>, |2,1,0>, |2,1,1>, |2,1,-1>
    !
    ! The perturbation H' = eEz = eEr*cos(theta) connects only
    ! states differing by Delta_l = +/-1 and Delta_m = 0.
    ! The only nonzero matrix element is:
    !   <2,0,0|H'|2,1,0> = -3 e a0 E_field
    !
    ! This produces a linear Stark splitting into 3 levels.
    ! ──────────────────────────────────────────────────────────────
    implicit none

integer, parameter :: dp = selected_real_kind(15)
    integer, parameter :: n_fields = 60
    integer, parameter :: n_states = 4

real(dp), parameter :: a0 = 5.29177210903e-11_dp   ! Bohr radius (m)
    real(dp), parameter :: e_charge = 1.602176634e-19_dp ! electron charge (C)
    real(dp), parameter :: eV = 1.602176634e-19_dp       ! J per eV
    real(dp), parameter :: E_hartree = 27.211386245988_dp ! Hartree in eV

real(dp) :: E_field, E_max, dE
    real(dp) :: H(n_states, n_states)
    real(dp) :: eigenvalues(n_states)
    real(dp) :: work(20)
    real(dp) :: V_12  ! off-diagonal matrix element
    real(dp) :: E0    ! unperturbed n=2 energy
    integer :: i, info

! Unperturbed energy of n=2 in eV: E = -13.6/4 = -3.4 eV
    E0 = -13.6_dp / 4.0_dp

! Maximum electric field (V/m) - typical lab fields
    E_max = 5.0e8_dp
    dE = E_max / real(n_fields - 1, dp)

write(*,'(A)') "# Stark Effect: Hydrogen n=2 Energy Levels"
    write(*,'(A)') "# E_field(V/m)      E1(eV)          E2(eV)          E3(eV)          E4(eV)"

do i = 1, n_fields
        E_field = real(i - 1, dp) * dE

! Build the 4x4 perturbation matrix in the basis:
        ! |1> = |2,0,0>,  |2> = |2,1,0>,  |3> = |2,1,1>,  |4> = |2,1,-1>
        H = 0.0_dp

! Diagonal elements: all equal to E0
        H(1,1) = E0
        H(2,2) = E0
        H(3,3) = E0
        H(4,4) = E0

! Off-diagonal: <2,0,0|eEz|2,1,0> = -3*a0*e*E_field
        ! Convert to eV: multiply by (1/eV)
        V_12 = -3.0_dp * a0 * e_charge * E_field / eV
        H(1,2) = V_12
        H(2,1) = V_12
        ! All other off-diagonal elements vanish by selection rules

! Diagonalize using a simple 4x4 eigenvalue routine
        ! For this symmetric matrix, we can solve analytically:
        ! Eigenvalues: E0 + 3*a0*e*E, E0, E0, E0 - 3*a0*e*E
        ! But let us use direct diagonalization for generality
        call diag_symmetric_4x4(H, eigenvalues)

write(*,'(ES16.6, 4(2X, F14.8))') E_field, &
            eigenvalues(1), eigenvalues(2), eigenvalues(3), eigenvalues(4)
    end do

! Print summary
    write(*,'(A)') ""
    write(*,'(A)') "# --- Summary at maximum field ---"
    write(*,'(A, ES12.4, A)') "# Electric field: ", E_max, " V/m"
    write(*,'(A, F12.6, A)') "# Splitting: ", &
        eigenvalues(4) - eigenvalues(1), " eV"
    write(*,'(A, F12.6, A)') "# = ", &
        (eigenvalues(4) - eigenvalues(1)) * 8065.54, " cm^-1"
    write(*,'(A)') "# Linear Stark effect confirmed: splitting proportional to E_field"

contains

subroutine diag_symmetric_4x4(A, evals)
        ! Simple Jacobi-like diagonalization for 4x4 symmetric matrix
        implicit none
        real(dp), intent(inout) :: A(4,4)
        real(dp), intent(out) :: evals(4)
        real(dp) :: B(4,4), V(4,4), off_diag, theta, c, s, t, tau
        real(dp) :: Bip, Bjp, Vpi, Vpj
        integer :: p, q, r, sweep, max_sweeps
        real(dp) :: threshold

B = A
        V = 0.0_dp
        do p = 1, 4
            V(p,p) = 1.0_dp
        end do

max_sweeps = 100
        threshold = 1.0e-14_dp

do sweep = 1, max_sweeps
            ! Check convergence
            off_diag = 0.0_dp
            do p = 1, 3
                do q = p+1, 4
                    off_diag = off_diag + B(p,q)**2
                end do
            end do
            if (off_diag < threshold) exit

! Sweep over all off-diagonal elements
            do p = 1, 3
                do q = p+1, 4
                    if (abs(B(p,q)) < 1.0e-15_dp) cycle

! Compute rotation angle
                    if (abs(B(p,p) - B(q,q)) < 1.0e-15_dp) then
                        theta = atan(1.0_dp)  ! pi/4
                    else
                        t = 2.0_dp * B(p,q) / (B(p,p) - B(q,q))
                        theta = 0.5_dp * atan(t)
                    end if

c = cos(theta)
                    s = sin(theta)

! Apply Jacobi rotation
                    do r = 1, 4
                        Bip = B(r,p)
                        Bjp = B(r,q)
                        B(r,p) = c * Bip + s * Bjp
                        B(r,q) = -s * Bip + c * Bjp
                    end do
                    do r = 1, 4
                        Bip = B(p,r)
                        Bjp = B(q,r)
                        B(p,r) = c * Bip + s * Bjp
                        B(q,r) = -s * Bip + c * Bjp
                    end do
                    do r = 1, 4
                        Vpi = V(r,p)
                        Vpj = V(r,q)
                        V(r,p) = c * Vpi + s * Vpj
                        V(r,q) = -s * Vpi + c * Vpj
                    end do
                end do
            end do
        end do

! Extract and sort eigenvalues
        do p = 1, 4
            evals(p) = B(p,p)
        end do
        ! Simple bubble sort
        do p = 1, 3
            do q = p+1, 4
                if (evals(q) < evals(p)) then
                    off_diag = evals(p)
                    evals(p) = evals(q)
                    evals(q) = off_diag
                end if
            end do
        end do
    end subroutine

end program stark_effect_hydrogen

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Key Equations Summary

Quantity	Expression
First-order energy	$E_k^{(1)} = \langle \psi_k^{(0)} \| \hat{H}' \| \psi_k^{(0)} \rangle$
First-order wavefunction	$\|\psi_k^{(1)}\rangle = \sum_{n \neq k} \frac{\langle \psi_n^{(0)} \| \hat{H}' \| \psi_k^{(0)} \rangle}{E_k^{(0)} - E_n^{(0)}} \|\psi_n^{(0)}\rangle$
Second-order energy	$E_k^{(2)} = \sum_{n \neq k} \frac{\|\langle \psi_n^{(0)} \| \hat{H}' \| \psi_k^{(0)} \rangle\|^2}{E_k^{(0)} - E_n^{(0)}}$
Degenerate PT secular eq.	$\det(\mathbf{H}' - E^{(1)} \mathbf{I}) = 0$
Stark splitting (n=2)	$E^{(1)} = \pm 3e\mathcal{E}a_0, \; 0$
Polarizability	$\alpha = 2 \sum_{n \neq 0} \frac{\|\langle n \| \mu_z \| 0 \rangle\|^2}{E_n - E_0}$

← Part I Overview Course Home →