🎯 Learning Objectives
- Understand the concept of the wavefunction and its probabilistic interpretation
- Learn the time-dependent and time-independent Schrödinger equations
- Master how to compute operators and the expectation values of physical quantities
- Derive the exact solution for a particle in a one-dimensional box
- Compute the eigenfunctions and energy levels of the harmonic oscillator
- Understand the uncertainty principle and commutation relations
- Learn angular momentum operators and spherical harmonics
- Implement simulations of quantum systems in Python
📖 The Wavefunction and the Schrödinger Equation
The Wavefunction and Its Probabilistic Interpretation
In quantum mechanics, the state of a particle is described by the wavefunction \(\psi(\mathbf{r}, t)\):
- \(|\psi(\mathbf{r}, t)|^2\) is the probability density of finding the particle at position \(\mathbf{r}\) at time \(t\)
- Normalization condition: \(\int |\psi(\mathbf{r}, t)|^2 d^3r = 1\)
- The wavefunction is a complex-valued function
Time-dependent Schrödinger equation:
\[
i\hbar \frac{\partial \psi}{\partial t} = \hat{H} \psi
\]
Here, \(\hat{H}\) is the Hamiltonian operator, which represents the total energy of the system.
Time-independent Schrödinger Equation
For a stationary state (an energy eigenstate), we can separate variables as \(\psi(\mathbf{r}, t) = \phi(\mathbf{r}) e^{-iEt/\hbar}\):
\[
\hat{H} \phi = E \phi
\]
This is the time-independent Schrödinger equation (an eigenvalue problem).
In one dimension:
\[
-\frac{\hbar^2}{2m} \frac{d^2\phi}{dx^2} + V(x)\phi = E\phi
\]
💻 Example 1.1: Particle in a One-Dimensional Box
Infinite Potential Well
A box of width \(L\), with potential \(V(x) = 0\) (\(0 < x < L\)) and \(V(x) = \infty\) (elsewhere):
Eigenfunctions:
\[
\phi_n(x) = \sqrt{\frac{2}{L}} \sin\left(\frac{n\pi x}{L}\right), \quad n = 1, 2, 3, \ldots
\]
Energy eigenvalues:
\[
E_n = \frac{n^2 \pi^2 \hbar^2}{2m L^2}
\]
Python Implementation: Particle in a Box
import numpy as np
import matplotlib.pyplot as plt
# Physical constants (atomic units)
hbar = 1.0
m = 1.0
L = 1.0
def energy_level(n, L, m, hbar):
"""Energy level"""
return (n**2 * np.pi**2 * hbar**2) / (2 * m * L**2)
def wavefunction(x, n, L):
"""Normalized eigenfunction"""
return np.sqrt(2/L) * np.sin(n * np.pi * x / L)
def probability_density(x, n, L):
"""Probability density"""
psi = wavefunction(x, n, L)
return np.abs(psi)**2
# Coordinates
x = np.linspace(0, L, 500)
# Different quantum numbers
quantum_numbers = [1, 2, 3, 4]
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# Wavefunctions
ax1 = axes[0, 0]
for n in quantum_numbers:
psi = wavefunction(x, n, L)
E_n = energy_level(n, L, m, hbar)
ax1.plot(x, psi, linewidth=2, label=f'n={n}, E={E_n:.2f}')
ax1.set_xlabel('Position x')
ax1.set_ylabel('Wavefunction ψ_n(x)')
ax1.set_title('Eigenfunctions (Wavefunctions)')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.axhline(0, color='k', linewidth=0.5)
# Probability density
ax2 = axes[0, 1]
for n in quantum_numbers:
rho = probability_density(x, n, L)
ax2.plot(x, rho, linewidth=2, label=f'n={n}')
ax2.set_xlabel('Position x')
ax2.set_ylabel('Probability density |ψ_n(x)|²')
ax2.set_title('Probability Density Distribution')
ax2.legend()
ax2.grid(True, alpha=0.3)
# Energy level diagram
ax3 = axes[1, 0]
n_max = 10
energies = [energy_level(n, L, m, hbar) for n in range(1, n_max+1)]
for i, (n, E) in enumerate(zip(range(1, n_max+1), energies)):
ax3.hlines(E, 0, 1, colors='blue', linewidth=2)
ax3.text(1.1, E, f'n={n}, E={E:.2f}', fontsize=10, va='center')
ax3.set_xlim([-0.1, 2])
ax3.set_ylim([0, max(energies) * 1.1])
ax3.set_ylabel('Energy E_n')
ax3.set_title('Energy Levels')
ax3.grid(True, alpha=0.3, axis='y')
ax3.set_xticks([])
# Expectation value calculation (expectation value of position)
ax4 = axes[1, 1]
expectation_x = []
for n in range(1, n_max+1):
psi = wavefunction(x, n, L)
# = ∫ ψ*(x) x ψ(x) dx
integrand = np.conj(psi) * x * psi
exp_x = np.trapz(integrand, x)
expectation_x.append(exp_x)
ax4.plot(range(1, n_max+1), expectation_x, 'o-', linewidth=2, markersize=8)
ax4.axhline(L/2, color='r', linestyle='--', linewidth=2, label='Classical (L/2)')
ax4.set_xlabel('Quantum number n')
ax4.set_ylabel('⟨x⟩')
ax4.set_title('Expectation Value of Position')
ax4.legend()
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('qchem_particle_in_box.png', dpi=300, bbox_inches='tight')
plt.show()
# Numerical results
print("=== Particle in a One-Dimensional Box ===\n")
print(f"Box width L = {L}")
print(f"Mass m = {m}")
print(f"\nEnergy levels:")
for n in range(1, 6):
E_n = energy_level(n, L, m, hbar)
print(f" n = {n}: E_{n} = {E_n:.4f} (atomic units)")
print(f"\nExpectation value of position (all x = L/2 = {L/2}):")
print(f" By symmetry, ⟨x⟩ = L/2 for every state")
💻 Example 1.2: Harmonic Oscillator
Quantum Harmonic Oscillator
A system with potential \(V(x) = \frac{1}{2}m\omega^2 x^2\):
Energy eigenvalues:
\[
E_n = \hbar\omega \left(n + \frac{1}{2}\right), \quad n = 0, 1, 2, \ldots
\]
Eigenfunctions:
\[
\phi_n(x) = \left(\frac{m\omega}{\pi\hbar}\right)^{1/4} \frac{1}{\sqrt{2^n n!}} H_n(\xi) e^{-\xi^2/2}
\]
Here, \(\xi = \sqrt{m\omega/\hbar} \cdot x\), and \(H_n\) are the Hermite polynomials.
Python Implementation: Quantum Harmonic Oscillator
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import hermite, factorial
# Physical constants
hbar = 1.0
m = 1.0
omega = 1.0
def energy_harmonic(n, omega, hbar):
"""Energy of the harmonic oscillator"""
return hbar * omega * (n + 0.5)
def harmonic_wavefunction(x, n, m, omega, hbar):
"""Eigenfunction of the harmonic oscillator"""
xi = np.sqrt(m * omega / hbar) * x
normalization = (m * omega / (np.pi * hbar))**0.25 / np.sqrt(2**n * factorial(n))
H_n = hermite(n)
psi = normalization * H_n(xi) * np.exp(-xi**2 / 2)
return psi
def classical_turning_point(n, omega, m):
"""Classical turning point (E = V)"""
E_n = energy_harmonic(n, omega, hbar)
return np.sqrt(2 * E_n / (m * omega**2))
# Coordinate range
x_max = 5
x = np.linspace(-x_max, x_max, 500)
# Potential
V = 0.5 * m * omega**2 * x**2
# Different quantum numbers
quantum_numbers = [0, 1, 2, 3, 5, 10]
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# Wavefunctions and potential
ax1 = axes[0, 0]
ax1.plot(x, V, 'k--', linewidth=2, label='Potential V(x)')
for n in [0, 1, 2, 3, 5]:
E_n = energy_harmonic(n, omega, hbar)
psi = harmonic_wavefunction(x, n, m, omega, hbar)
# Shift by the energy level for display
ax1.plot(x, E_n + psi * 2, linewidth=1.5, label=f'n={n}')
ax1.hlines(E_n, -x_max, x_max, colors='gray', linestyles=':', linewidth=0.8)
ax1.set_xlabel('Position x')
ax1.set_ylabel('Energy / Wavefunction')
ax1.set_title('Harmonic Oscillator Eigenfunctions and Energy Levels')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_ylim([0, 12])
# Probability density (low quantum numbers)
ax2 = axes[0, 1]
for n in [0, 1, 2, 3]:
psi = harmonic_wavefunction(x, n, m, omega, hbar)
rho = np.abs(psi)**2
ax2.plot(x, rho, linewidth=2, label=f'n={n}')
ax2.set_xlabel('Position x')
ax2.set_ylabel('Probability density |ψ_n(x)|²')
ax2.set_title('Probability Density Distribution (Low Excited States)')
ax2.legend()
ax2.grid(True, alpha=0.3)
# Classical limit (high quantum numbers)
ax3 = axes[1, 0]
n_high = 10
psi_high = harmonic_wavefunction(x, n_high, m, omega, hbar)
rho_high = np.abs(psi_high)**2
# Classical probability density
E_classical = energy_harmonic(n_high, omega, hbar)
x_turn = classical_turning_point(n_high, omega, m)
rho_classical = np.zeros_like(x)
mask = np.abs(x) < x_turn
# Classically, probability ∝ 1/v ∝ 1/√(E - V)
rho_classical[mask] = 1 / np.sqrt(E_classical - 0.5 * m * omega**2 * x[mask]**2)
rho_classical /= np.trapz(rho_classical, x) # Normalization
ax3.plot(x, rho_high, 'b-', linewidth=2, label=f'Quantum (n={n_high})')
ax3.plot(x, rho_classical, 'r--', linewidth=2, label='Classical')
ax3.axvline(-x_turn, color='k', linestyle=':', linewidth=1, label='Turning points')
ax3.axvline(x_turn, color='k', linestyle=':', linewidth=1)
ax3.set_xlabel('Position x')
ax3.set_ylabel('Probability density')
ax3.set_title(f'Classical Limit (n={n_high})')
ax3.legend()
ax3.grid(True, alpha=0.3)
# Zero-point energy
ax4 = axes[1, 1]
n_range = np.arange(0, 20)
E_quantum = [energy_harmonic(n, omega, hbar) for n in n_range]
E_classical = n_range * hbar * omega # Classical energy (ground state = 0)
ax4.plot(n_range, E_quantum, 'bo-', linewidth=2, markersize=6, label='Quantum')
ax4.plot(n_range, E_classical, 'r--', linewidth=2, label='Classical (E=nℏω)')
ax4.fill_between(n_range, E_classical, E_quantum, alpha=0.3, color='yellow',
label='Zero-point energy')
ax4.set_xlabel('Quantum number n')
ax4.set_ylabel('Energy E_n')
ax4.set_title('Zero-Point Energy (E_0 = ℏω/2)')
ax4.legend()
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('qchem_harmonic_oscillator.png', dpi=300, bbox_inches='tight')
plt.show()
# Numerical results
print("\n=== Quantum Harmonic Oscillator ===\n")
print(f"Angular frequency ω = {omega}")
print(f"Mass m = {m}")
print(f"Zero-point energy E_0 = {energy_harmonic(0, omega, hbar):.4f}\n")
print("Energy levels:")
for n in range(6):
E_n = energy_harmonic(n, omega, hbar)
print(f" n = {n}: E_{n} = {E_n:.4f} = {n + 0.5:.1f}ℏω")
print(f"\nClassical turning points:")
for n in [0, 1, 5, 10]:
x_turn = classical_turning_point(n, omega, m)
print(f" n = {n}: x_turn = ±{x_turn:.4f}")
💻 Example 1.3: Operators and Expectation Values
Operators in Quantum Mechanics
Physical quantities are represented by operators:
- Position operator: \(\hat{x} = x\) (multiplication by position)
- Momentum operator: \(\hat{p} = -i\hbar \frac{\partial}{\partial x}\)
- Energy operator: \(\hat{H} = -\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2} + V(x)\)
Expectation value:
\[
\langle A \rangle = \int \psi^* \hat{A} \psi \, dx
\]
Ehrenfest's theorem:
\[
\frac{d\langle x \rangle}{dt} = \frac{\langle p \rangle}{m}, \quad
\frac{d\langle p \rangle}{dt} = -\left\langle \frac{\partial V}{\partial x} \right\rangle
\]
Python Implementation: Operators and Expectation Values
import numpy as np
import matplotlib.pyplot as plt
hbar = 1.0
m = 1.0
omega = 1.0
def position_expectation(n, omega, m, hbar):
"""Expectation value of position (zero by symmetry)"""
return 0.0 # For the harmonic oscillator, ⟨x⟩ = 0 by symmetry
def momentum_expectation(n, omega, m, hbar):
"""Expectation value of momentum (zero by symmetry)"""
return 0.0
def position_uncertainty(n, omega, m, hbar):
"""Position uncertainty Δx"""
# ⟨x²⟩ = (n + 1/2) ℏ/(mω)
x_squared = (n + 0.5) * hbar / (m * omega)
return np.sqrt(x_squared)
def momentum_uncertainty(n, omega, m, hbar):
"""Momentum uncertainty Δp"""
# ⟨p²⟩ = (n + 1/2) mωℏ
p_squared = (n + 0.5) * m * omega * hbar
return np.sqrt(p_squared)
# Verify the uncertainty relation
n_range = np.arange(0, 20)
Delta_x = np.array([position_uncertainty(n, omega, m, hbar) for n in n_range])
Delta_p = np.array([momentum_uncertainty(n, omega, m, hbar) for n in n_range])
product = Delta_x * Delta_p
heisenberg_limit = hbar / 2
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# Position uncertainty
ax1 = axes[0, 0]
ax1.plot(n_range, Delta_x, 'bo-', linewidth=2, markersize=6)
ax1.set_xlabel('Quantum number n')
ax1.set_ylabel('Δx')
ax1.set_title('Position Uncertainty')
ax1.grid(True, alpha=0.3)
# Momentum uncertainty
ax2 = axes[0, 1]
ax2.plot(n_range, Delta_p, 'ro-', linewidth=2, markersize=6)
ax2.set_xlabel('Quantum number n')
ax2.set_ylabel('Δp')
ax2.set_title('Momentum Uncertainty')
ax2.grid(True, alpha=0.3)
# Uncertainty relation
ax3 = axes[1, 0]
ax3.plot(n_range, product, 'go-', linewidth=2, markersize=6, label='Δx·Δp')
ax3.axhline(heisenberg_limit, color='r', linestyle='--', linewidth=2,
label=f'Heisenberg limit (ℏ/2 = {heisenberg_limit})')
ax3.set_xlabel('Quantum number n')
ax3.set_ylabel('Δx · Δp')
ax3.set_title('Heisenberg Uncertainty Relation')
ax3.legend()
ax3.grid(True, alpha=0.3)
# Action of the momentum operator via numerical differentiation
ax4 = axes[1, 1]
x = np.linspace(-5, 5, 500)
dx = x[1] - x[0]
n = 2
psi = harmonic_wavefunction(x, n, m, omega, hbar)
# Momentum operator p̂ψ = -iℏ dψ/dx
dpsi_dx = np.gradient(psi, dx)
p_psi = -1j * hbar * dpsi_dx
ax4.plot(x, np.real(psi), 'b-', linewidth=2, label='Re(ψ)')
ax4.plot(x, np.real(p_psi), 'r-', linewidth=2, label='Re(p̂ψ)')
ax4.plot(x, np.imag(p_psi), 'g--', linewidth=2, label='Im(p̂ψ)')
ax4.set_xlabel('Position x')
ax4.set_ylabel('Amplitude')
ax4.set_title(f'Action of the Momentum Operator (n={n})')
ax4.legend()
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('qchem_operators_expectation.png', dpi=300, bbox_inches='tight')
plt.show()
# Numerical results
print("\n=== Operators and Expectation Values ===\n")
print("Expectation values of the harmonic oscillator:\n")
for n in [0, 1, 2, 5]:
Dx = position_uncertainty(n, omega, m, hbar)
Dp = momentum_uncertainty(n, omega, m, hbar)
print(f"n = {n}:")
print(f" ⟨x⟩ = {position_expectation(n, omega, m, hbar):.4f} (symmetry)")
print(f" ⟨p⟩ = {momentum_expectation(n, omega, m, hbar):.4f} (symmetry)")
print(f" Δx = {Dx:.4f}")
print(f" Δp = {Dp:.4f}")
print(f" Δx·Δp = {Dx * Dp:.4f} (≥ ℏ/2 = {hbar/2:.4f})\n")
print("Heisenberg uncertainty relation:")
print(f" In the ground state (n=0), Δx·Δp = ℏ/2 (minimum-uncertainty state)")
💻 Example 1.4: Angular Momentum and Spherical Harmonics
Angular Momentum Operators
The angular momentum operator in three-dimensional space:
\[
\hat{\mathbf{L}} = \mathbf{r} \times \hat{\mathbf{p}} = -i\hbar (\mathbf{r} \times \nabla)
\]
Commutation relations:
\[
[\hat{L}_x, \hat{L}_y] = i\hbar \hat{L}_z, \quad [\hat{L}_y, \hat{L}_z] = i\hbar \hat{L}_x, \quad [\hat{L}_z, \hat{L}_x] = i\hbar \hat{L}_y
\]
Since \([\hat{L}^2, \hat{L}_z] = 0\), \(\hat{L}^2\) and \(\hat{L}_z\) share simultaneous eigenstates.
Eigenvalues:
\[
\hat{L}^2 Y_l^m(\theta, \phi) = \hbar^2 l(l+1) Y_l^m(\theta, \phi)
\]
\[
\hat{L}_z Y_l^m(\theta, \phi) = \hbar m Y_l^m(\theta, \phi)
\]
Here, \(Y_l^m(\theta, \phi)\) are the spherical harmonics.
Python Implementation: Spherical Harmonics
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import sph_harm
from mpl_toolkits.mplot3d import Axes3D
def plot_spherical_harmonic(l, m, ax, title):
"""Visualize a spherical harmonic"""
# Spherical coordinates
theta = np.linspace(0, np.pi, 100)
phi = np.linspace(0, 2*np.pi, 100)
Theta, Phi = np.meshgrid(theta, phi)
# Spherical harmonic (scipy convention: Y_l^m(phi, theta))
Y_lm = sph_harm(m, l, Phi, Theta)
# Absolute value (angular dependence of the probability density)
R = np.abs(Y_lm)
# Cartesian coordinates
X = R * np.sin(Theta) * np.cos(Phi)
Y = R * np.sin(Theta) * np.sin(Phi)
Z = R * np.cos(Theta)
# Plot
surface = ax.plot_surface(X, Y, Z, cmap='viridis',
facecolors=plt.cm.viridis(R/R.max()),
alpha=0.9, shade=True)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title(title)
ax.set_box_aspect([1,1,1])
# Spherical harmonics for different (l, m)
fig = plt.figure(figsize=(16, 12))
configs = [
(0, 0, 's orbital (l=0, m=0)'),
(1, 0, 'p_z orbital (l=1, m=0)'),
(1, 1, 'p_x orbital (l=1, m=1)'),
(2, 0, 'd_{z²} orbital (l=2, m=0)'),
(2, 1, 'd_{xz} orbital (l=2, m=1)'),
(2, 2, 'd_{xy} orbital (l=2, m=2)'),
(3, 0, 'f_{z³} orbital (l=3, m=0)'),
(3, 2, 'f orbital (l=3, m=2)'),
]
for idx, (l, m, title) in enumerate(configs):
ax = fig.add_subplot(2, 4, idx+1, projection='3d')
plot_spherical_harmonic(l, m, ax, title)
plt.tight_layout()
plt.savefig('qchem_spherical_harmonics.png', dpi=300, bbox_inches='tight')
plt.show()
# Eigenvalues of angular momentum
print("\n=== Angular Momentum and Spherical Harmonics ===\n")
print("Eigenvalues of angular momentum:\n")
hbar = 1.0
for l in range(4):
L_squared_eigenvalue = hbar**2 * l * (l + 1)
print(f"l = {l}:")
print(f" Eigenvalue of L² = {L_squared_eigenvalue:.4f} = ℏ²·{l}·{l+1}")
print(f" Allowed values of m: from {-l} to {l} ({2*l+1} values)")
for m in range(-l, l+1):
Lz_eigenvalue = hbar * m
print(f" m = {m:2d}: L_z = {Lz_eigenvalue:+.4f} = {m:+d}ℏ")
print()
print("Orthogonality of spherical harmonics:")
print(" ∫ Y_l^m* Y_l'^m' dΩ = δ_{ll'} δ_{mm'}")
📚 Summary
- The wavefunction describes the quantum state of a particle and provides the probability amplitude
- The Schrödinger equation is the fundamental equation of quantum mechanics, governing the time evolution of the wavefunction
- Stationary states (energy eigenstates) are found by solving an eigenvalue problem
- The particle in a box and the harmonic oscillator are important models with exact solutions
- Quantization gives rise to discrete energy levels
- The zero-point energy is a consequence of quantum effects: even in the ground state the particle is in motion
- Operators represent physical quantities, and expectation values give the average of measurements
- The uncertainty relation is an essential constraint of quantum mechanics: position and momentum cannot be determined precisely at the same time
- Angular momentum is quantized and described by spherical harmonics
- Spherical harmonics give the angular dependence of atomic orbitals (s, p, d, f orbitals)
💡 Exercises
- Normalization of a particle: Numerically verify that the wavefunction of the particle in a box (\(n=1\)) is normalized.
- Tunneling effect: Compute the transmission probability through a finite potential barrier and discuss the difference from classical mechanics.
- Operator method for the harmonic oscillator: Construct the eigenfunctions using the creation and annihilation operators \(\hat{a}^\dagger, \hat{a}\).
- Virial theorem: Verify that \(2\langle T \rangle = \langle V \rangle\) holds for the harmonic oscillator.
- Orthogonality of spherical harmonics: Numerically verify by integration that spherical harmonics with different \((l, m)\) are orthogonal.