Relativistic Hydrogen Atom

The relativistic_hydrogen.py module implements a complete Dirac equation solver for hydrogen-like atoms, including fine structure, Zitterbewegung, and spin-orbit coupling.

Theory Background

Dirac Equation

The Dirac equation for a relativistic electron in Coulomb potential:

(cα̂·p̂ + βmc² + V(r))ψ = Eψ

where:
  α̂ = (α₁, α₂, α₃) = Dirac alpha matrices
  β = γ⁰ = Dirac gamma matrix
  c ≈ 137.036 a.u. (speed of light)
  V(r) = -Z/r (Coulomb potential)

Fine Structure Energy Levels

From relativistic_hydrogen.py:526-554, the exact Dirac energy:

E(n,κ) = mc² / sqrt(1 + (αZ)^2 / (n - |κ| + sqrt(κ² - (αZ)²))^2)

Binding energy:
E_b = (E - mc²) = E - c²  (in atomic units)

where:
  α = 1/137.036 (fine structure constant)
  κ = relativistic quantum number
       = -(l+1) for j = l + 1/2
       = l     for j = l - 1/2

Installation Check

Verify the module loads correctly:

python -c "from relativistic_hydrogen import Config; print('OK')"

Basic Usage

Energy level calculation

from relativistic_hydrogen import DiracHydrogenAtom, Config

config = Config()
atom = DiracHydrogenAtom(config)

# Calculate energy for n=2, l=1, j=3/2
n, l, j = 2, 1, 1.5

if j == l + 0.5:
    kappa = -(l + 1)  # j = l + 1/2
else:
    kappa = l         # j = l - 1/2

E_dirac = atom.energy_level_dirac(n, kappa)
E_schrodinger = -0.5 / n**2

print(f"n={n}, l={l}, j={j}")
print(f"Dirac energy: {E_dirac:.8f} Ha")
print(f"Schrodinger energy: {E_schrodinger:.8f} Ha")
print(f"Fine structure shift: {E_dirac - E_schrodinger:.10f} Ha")

Generate energy spectrum

spectrum = atom.energy_spectrum(n_max=4)

for level in spectrum:
    print(f"{level['notation']:>8s}: "
          f"E = {level['energy']:.8f} Ha, "
          f"degeneracy = {level['degeneracy']}")

Fine structure splitting

# For n=2, l=1 (2p orbital)
fs = atom.fine_structure_splitting(n=2, l=1)

print(f"2p splitting:")
print(f"  j = {fs['j_upper']}: E = {fs['E_j_upper']:.8f} Ha")
print(f"  j = {fs['j_lower']}: E = {fs['E_j_lower']:.8f} Ha")
print(f"  Splitting: {fs['splitting']:.10f} Ha")

Gamma Matrices

From relativistic_hydrogen.py:112-193, the Dirac gamma matrices in standard representation:

Gamma Matrices
Anticommutation

class GammaMatrices:
    def __init__(self, device='cpu'):
        # γ⁰ (beta)
        self.gamma0 = torch.tensor([
            [1, 0, 0, 0],
            [0, 1, 0, 0],
            [0, 0, -1, 0],
            [0, 0, 0, -1]
        ], dtype=torch.complex128)
        
        # γ¹
        self.gamma1 = torch.tensor([
            [0, 0, 0, 1],
            [0, 0, 1, 0],
            [0, -1, 0, 0],
            [-1, 0, 0, 0]
        ], dtype=torch.complex128)
        
        # γ²
        self.gamma2 = torch.tensor([
            [0, 0, 0, -1j],
            [0, 0, 1j, 0],
            [0, 1j, 0, 0],
            [-1j, 0, 0, 0]
        ], dtype=torch.complex128)
        
        # γ³
        self.gamma3 = torch.tensor([
            [0, 0, 1, 0],
            [0, 0, 0, -1],
            [-1, 0, 0, 0],
            [0, 1, 0, 0]
        ], dtype=torch.complex128)
        
        # Alpha matrices: α_i = γ⁰ γ^i
        self.alpha_x = self.gamma0 @ self.gamma1
        self.alpha_y = self.gamma0 @ self.gamma2
        self.alpha_z = self.gamma0 @ self.gamma3
        
        # Beta = γ⁰
        self.beta = self.gamma0

Verify gamma matrix anticommutation relations:

gamma = GammaMatrices()

# {γμ, γν} = 2gμν
for mu in range(4):
    for nu in range(4):
        anticomm = (
            gamma.gammas[mu] @ gamma.gammas[nu] +
            gamma.gammas[nu] @ gamma.gammas[mu]
        )
        
        # Metric: g = diag(1, -1, -1, -1)
        expected = 2 * (1 if mu == nu and mu == 0 else
                      -1 if mu == nu else 0)
        
        result = torch.allclose(
            anticomm, 
            expected * torch.eye(4, dtype=torch.complex128)
        )
        
        print(f"{{γ^{mu}, γ^{nu}}} = {expected}: {result}")

Dirac Hamiltonian Operator

From relativistic_hydrogen.py:199-304:

class DiracHamiltonianOperator:
    """H_Dirac = cα̂·p̂ + βmc² + V(r)"""
    
    def apply_dirac_hamiltonian(self, spinor, potential=None):
        """
        Apply Dirac Hamiltonian to 4-component spinor.
        
        Args:
            spinor: [4, H, W] or [batch, 4, H, W]
            potential: Optional scalar potential V(r)
        
        Returns:
            Hψ with same shape
        """
        # Kinetic term: cα̂·p̂
        for c in range(4):
            psi_c = spinor[:, c, :, :]
            
            # Momentum in Fourier space: p̂ = -i∇
            psi_c_fft = torch.fft.fft2(psi_c)
            px_psi = torch.fft.ifft2(self.kx_grid * psi_c_fft)
            py_psi = torch.fft.ifft2(self.ky_grid * psi_c_fft)
            
            # Apply alpha matrices
            for d in range(4):
                result[:, c, :, :] += self.c * (
                    self.gamma.alpha_x[c, d] * px_psi +
                    self.gamma.alpha_y[c, d] * py_psi
                )
        
        # Mass term: βmc²
        for c in range(4):
            for d in range(4):
                result[:, c, :, :] += (
                    self.gamma.beta[c, d] * self.mass * self.c**2 *
                    spinor[:, d, :, :]
                )
        
        # Potential term: V(r)ψ
        if potential is not None:
            for c in range(4):
                result[:, c, :, :] += potential * spinor[:, c, :, :]
        
        return result
    
    def time_evolution(self, spinor, dt, potential=None):
        """First-order split-step: ψ(t+dt) ≈ (1 - iH dt)ψ(t)"""
        H_psi = self.apply_dirac_hamiltonian(spinor, potential)
        result = spinor - 1j * dt * H_psi
        
        # Normalize to preserve probability
        norm_original = torch.norm(spinor.view(spinor.shape[0], -1), dim=1)
        norm_evolved = torch.norm(result.view(result.shape[0], -1), dim=1)
        result = result * (norm_original / norm_evolved).unsqueeze(-1).unsqueeze(-1)
        
        return result

Zitterbewegung Simulation

Zitterbewegung is the “trembling motion” of a relativistic electron caused by interference between positive and negative energy states.

Theory

Frequency: ω = 2mc²/ℏ ≈ 2c² ≈ 2 × 137² ≈ 37,500 a.u.⁻¹ Amplitude: Δx ∼ ℏ/(2mc) ∼ 1/(2c) ≈ 0.004 a.u. ≈ 2 × 10⁻³ Å

Implementation

From relativistic_hydrogen.py:639-827:

class ZitterbewegungSimulator:
    def create_gaussian_wave_packet(self, sigma=0.1, momentum=0.0):
        """Create Gaussian wave packet for free particle."""
        grid_size = self.config.GRID_SIZE
        x = torch.linspace(-np.pi, np.pi, grid_size)
        y = torch.linspace(-np.pi, np.pi, grid_size)
        X, Y = torch.meshgrid(x, y, indexing='ij')
        
        # Gaussian envelope
        gaussian = torch.exp(-(X**2 + Y**2) / (2 * sigma**2))
        
        # 4-component spinor
        spinor = torch.zeros((4, grid_size, grid_size), 
                           dtype=torch.complex128)
        
        # Upper components (positive energy)
        spinor[0] = gaussian          # spin up
        spinor[1] = gaussian * 0.5j   # small spin down
        
        # Lower components (negative energy admixture for ZBW)
        spinor[2] = gaussian * 0.01   # small lower component
        spinor[3] = gaussian * 0.01j
        
        # Add momentum
        if momentum != 0:
            phase = torch.exp(1j * momentum * X)
            for i in range(4):
                spinor[i] = spinor[i] * phase
        
        # Normalize
        norm = torch.sqrt(torch.sum(torch.abs(spinor)**2))
        spinor = spinor / norm
        
        return spinor

Visualization

Energy Spectrum
Orbital Visualization
Zitterbewegung

from relativistic_hydrogen import (
    DiracHydrogenAtom, DiracVisualizer, Config
)

config = Config()
atom = DiracHydrogenAtom(config)
visualizer = DiracVisualizer(config)

# Generate spectrum
spectrum = atom.energy_spectrum(n_max=4)

# Visualize
visualizer.visualize_energy_spectrum(
    spectrum,
    save_path="dirac_energy_spectrum.png"
)

Output shows:

Fine structure splitting for each n, l
Comparison with non-relativistic energies
Degeneracies

from relativistic_hydrogen import (
    DiracMonteCarloSampler, DiracModelWrapper, Config
)

config = Config()
model = DiracModelWrapper(config)
sampler = DiracMonteCarloSampler(config, model)
visualizer = DiracVisualizer(config)

# Sample 2p_{3/2} orbital
data = sampler.sample_orbital(
    n=2, l=1, j=1.5,
    num_samples=100000
)

# Visualize
visualizer.visualize_orbital(
    data,
    save_path="orbital_2p_3half.png"
)

from relativistic_hydrogen import ZitterbewegungSimulator

zbw_sim = ZitterbewegungSimulator(config, model)

zbw_data = zbw_sim.simulate(
    duration=1.0,   # a.u.
    dt=0.001,       # a.u.
    sigma=0.1       # wave packet width
)

visualizer.visualize_zitterbewegung(
    zbw_data,
    save_path="zitterbewegung.png"
)

print(f"ZBW frequency: {zbw_data['zitterbewegung_frequency']:.1f} a.u.")
print(f"Expected: {zbw_data['zitterbewegung_frequency_expected']:.1f} a.u.")
print(f"Amplitude: {zbw_data['zitterbewegung_amplitude']:.6f} a.u.")

Model Loading

The module can use trained Dirac spectral networks (relativistic_hydrogen.py:391-510):

class DiracModelWrapper:
    def _load_model(self):
        # Look for checkpoint in checkpoints_dirac_phase4/
        checkpoint_path = self._find_best_checkpoint()
        
        if checkpoint_path:
            model = DiracSpectralNetwork(
                grid_size=16,
                hidden_dim=32,
                expansion_dim=64,
                num_spectral_layers=2,
                spinor_components=4
            )
            
            checkpoint = torch.load(checkpoint_path)
            model.load_state_dict(checkpoint['model_state_dict'])
            model.eval()
            
            self.model = model
            self.is_loaded = True
        else:
            # Use analytical Dirac operator
            self.is_loaded = False

If no trained model is found, the analytical Dirac operator is used.

Configuration

From relativistic_hydrogen.py:40-89:

@dataclass
class Config:
    # Model
    GRID_SIZE: int = 16
    HIDDEN_DIM: int = 32
    EXPANSION_DIM: int = 64
    NUM_SPECTRAL_LAYERS: int = 2
    SPINOR_COMPONENTS: int = 4
    
    # Checkpoints
    CHECKPOINT_DIR: str = 'checkpoints_dirac_phase4'
    CHECKPOINT_BEST_ALPHA: str = 'latest.pth'
    
    # Physical constants (atomic units)
    HBAR: float = 1.0
    ELECTRON_MASS: float = 1.0
    C_LIGHT: float = 137.035999084
    ALPHA_FS: float = 1.0 / 137.035999084
    
    # Visualization
    FIGURE_DPI: int = 150
    FIGURE_SIZE_X: int = 24
    FIGURE_SIZE_Y: int = 20
    
    # Monte Carlo
    MONTE_CARLO_BATCH_SIZE: int = 100000
    MONTE_CARLO_MAX_PARTICLES: int = 2000000
    
    # Zitterbewegung
    ZBW_TIME_STEPS: int = 1000
    ZBW_DT: float = 0.001
    ZBW_PACKET_WIDTH: float = 0.1

Key Results

Fine Structure Constant

α = 1/137.036 is THE fundamental constant determining:

Fine structure splitting magnitude
Lamb shift
Relativistic corrections

Energy Level Formula

For hydrogen (Z=1):

# n=2, l=1 (2p orbital)
j_upper = 3/2  # 2p_{3/2}
j_lower = 1/2  # 2p_{1/2}

E_nonrel = -0.5 / 2**2 = -0.125 Ha

# Dirac energies
kappa_upper = -(l+1) = -2
kappa_lower = l = 1

E_2p_3half = -0.12499590 Ha  # j=3/2
E_2p_1half = -0.12499820 Ha  # j=1/2

Splitting = E_2p_3half - E_2p_1half = 2.3 × 10⁻⁶ Ha ≈ 0.05 μeV

This matches experimental fine structure splitting.

Command Line Usage

# Energy spectrum
python -c "
from relativistic_hydrogen import DiracHydrogenAtom, Config
atom = DiracHydrogenAtom(Config())
for level in atom.energy_spectrum(n_max=3):
    print(f\"{level['notation']:>8s}: E = {level['energy']:.8f} Ha\")
"

# Orbital visualization (requires full script)
python relativistic_hydrogen.py

Get Started

Core Concepts

Guides

Advanced

Relativistic Hydrogen

Relativistic Hydrogen Atom

Theory Background

Dirac Equation

Fine Structure Energy Levels

Installation Check

Basic Usage

Gamma Matrices

Dirac Hamiltonian Operator

Zitterbewegung Simulation

Theory

Implementation

Visualization

Model Loading

Configuration

Key Results

Fine Structure Constant

Energy Level Formula

Command Line Usage

Next Steps

Visualization

Quantum Algorithms

Build docs developers (and LLMs) love

Get Started

Core Concepts

Guides

Advanced

​Relativistic Hydrogen Atom

​Theory Background

​Dirac Equation

​Fine Structure Energy Levels

​Installation Check

​Basic Usage

​Gamma Matrices

​Dirac Hamiltonian Operator

​Zitterbewegung Simulation

​Theory

​Implementation

​Visualization

​Model Loading

​Configuration

​Key Results

​Fine Structure Constant

​Energy Level Formula

​Command Line Usage

​Next Steps

Visualization

Quantum Algorithms

Build docs developers (and LLMs) love

Relativistic Hydrogen Atom

Theory Background

Dirac Equation

Fine Structure Energy Levels

Installation Check

Basic Usage

Gamma Matrices

Dirac Hamiltonian Operator

Zitterbewegung Simulation

Theory

Implementation

Visualization

Model Loading

Configuration

Key Results

Fine Structure Constant

Energy Level Formula

Command Line Usage

Next Steps