Quantum Mechanics

The quantum domain provides tools for semiclassical quantum mechanics, WKB approximation, and exact solutions for standard potentials.

Overview 

The quantum module implements:

WKB approximation for semiclassical analysis
Bohr-Sommerfeld quantization for bound states
Quantum harmonic oscillator exact solutions
Infinite square well (particle in a box)
Ehrenfest dynamics (quantum-classical correspondence)

Quick Start 

from mechanics_dsl.domains.quantum import (
    QuantumHarmonicOscillator,
    InfiniteSquareWell,
    WKBApproximation
)
import numpy as np

# Quantum harmonic oscillator
qho = QuantumHarmonicOscillator(mass=1.0, omega=1.0, hbar=1.0)

# Energy levels: E_n = ℏω(n + 1/2)
E0 = qho.energy_level(0)  # 0.5 (ground state)
E1 = qho.energy_level(1)  # 1.5

# Wavefunction
x = np.linspace(-5, 5, 1000)
psi_0 = qho.wavefunction(x, n=0)

Classes 

QuantumHarmonicOscillator 

The quantum harmonic oscillator has the Hamiltonian:

\[H = \frac{p^2}{2m} + \frac{1}{2}m\omega^2 x^2\]

with exact energy eigenvalues \(E_n = \hbar\omega(n + \frac{1}{2})\).

class mechanics_dsl.domains.quantum.QuantumHarmonicOscillator(mass: float = 1.0, omega: float = 1.0, hbar: float = 1.0)[source]

Bases: object

Exact quantum harmonic oscillator solution.

H = p²/(2m) + (1/2)mω²x²

Energy levels: E_n = ℏω(n + 1/2)

energy_level(n: int) → float[source]

Exact energy eigenvalue.

E_n = ℏω(n + 1/2)

zero_point_energy() → float[source]: Ground state energy E_0 = ℏω/2.

characteristic_length() → float[source]

Characteristic length scale a = √(ℏ/(mω)).

This is the ground state width.

classical_amplitude(n: int) → float[source]

Classical turning point for energy level n.

x_max = √(2E_n / (mω²))

wavefunction(x: ndarray, n: int) → ndarray[source]

Normalized wavefunction ψ_n(x).

Uses Hermite polynomials.

probability_density(x: ndarray, n: int) → ndarray[source]: Probability density |ψ_n(x)|².

position_expectation(n: int) → float[source]: Expectation value <x> = 0 for all n.

position_variance(n: int) → float[source]: Variance <x²> = a²(n + 1/2).

momentum_variance(n: int) → float[source]: Variance <p²> = (ℏ/a)²(n + 1/2).

uncertainty_product(n: int) → float[source]

Uncertainty product Δx·Δp = ℏ(n + 1/2).

Minimum (ℏ/2) for ground state n=0.

Example: Ground state properties

from mechanics_dsl.domains.quantum import QuantumHarmonicOscillator
import numpy as np
import matplotlib.pyplot as plt

qho = QuantumHarmonicOscillator(mass=1.0, omega=1.0, hbar=1.0)

# Ground state energy (zero-point energy)
E0 = qho.zero_point_energy()  # ℏω/2 = 0.5

# Characteristic length scale
a = qho.characteristic_length()  # √(ℏ/mω) = 1.0

# Classical turning point
x_classical = qho.classical_amplitude(n=0)

# Plot probability density
x = np.linspace(-4, 4, 500)
prob = qho.probability_density(x, n=0)

plt.plot(x, prob)
plt.xlabel('x')
plt.ylabel('|ψ₀(x)|²')
plt.title('Ground State Probability Density')

Example: Uncertainty principle verification

from mechanics_dsl.domains.quantum import QuantumHarmonicOscillator

qho = QuantumHarmonicOscillator(hbar=1.0)

# Ground state has minimum uncertainty
delta_x_delta_p = qho.uncertainty_product(n=0)  # ℏ/2 = 0.5

# Excited states have larger uncertainty
for n in range(5):
    product = qho.uncertainty_product(n)
    print(f"n={n}: ΔxΔp = {product:.1f}ℏ")

InfiniteSquareWell 

Particle in a box with walls at x=0 and x=L:

\[E_n = \frac{n^2\pi^2\hbar^2}{2mL^2}, \quad n = 1, 2, 3, \ldots\]

class mechanics_dsl.domains.quantum.InfiniteSquareWell(length: float = 1.0, mass: float = 1.0, hbar: float = 1.0)[source]

Bases: object

Particle in an infinite square well (1D box).

V(x) = 0 for 0 < x < L V(x) = ∞ otherwise

Energy levels: E_n = n²π²ℏ²/(2mL²)

energy_level(n: int) → float[source]

Energy eigenvalue for quantum number n (n = 1, 2, 3, …).

E_n = n²π²ℏ²/(2mL²)

wavefunction(x: ndarray, n: int) → ndarray[source]: Normalized wavefunction ψ_n(x) = √(2/L) sin(nπx/L).

probability_density(x: ndarray, n: int) → ndarray[source]: Probability density |ψ_n(x)|².

position_expectation(n: int) → float[source]: <x> = L/2 for all n.

position_variance(n: int) → float[source]: <x²> - <x>² for level n.

Example: Energy level spacing

from mechanics_dsl.domains.quantum import InfiniteSquareWell

well = InfiniteSquareWell(length=1.0, mass=1.0, hbar=1.0)

# First few energy levels
for n in range(1, 6):
    E = well.energy_level(n)
    print(f"E_{n} = {E:.3f}")

# Note: E_n ∝ n²

WKBApproximation 

The WKB (Wentzel-Kramers-Brillouin) approximation is valid when the de Broglie wavelength varies slowly compared to the potential.

Bohr-Sommerfeld quantization condition:

\[\oint p \, dx = \left(n + \frac{1}{2}\right) h\]

class mechanics_dsl.domains.quantum.WKBApproximation(potential: Callable[[float], float], mass: float = 1.0, hbar: float = 1.0)[source]

Bases: object

Implements the WKB (Wentzel-Kramers-Brillouin) approximation.

Valid in the semiclassical limit where the de Broglie wavelength varies slowly compared to the potential.

The WKB wavefunction is:: ψ(x) ≈ C/√p(x) * exp(±i/ℏ ∫p(x’)dx’)

where p(x) = √(2m(E-V(x))) is the classical momentum.

Example

>>> wkb = WKBApproximation(potential=lambda x: 0.5*x**2, mass=1.0)
>>> levels = wkb.bohr_sommerfeld_levels(n_max=10)

classical_momentum(x: float, E: float) → float[source]

Calculate classical momentum p(x) = √(2m(E-V(x))).

Returns 0 if E < V(x) (classically forbidden region).

find_turning_points(E: float, x_range: Tuple[float, float], n_points: int = 1000) → List[float][source]

Find classical turning points where E = V(x).

Parameters:

E – Total energy
x_range – (x_min, x_max) search range
n_points – Number of grid points

Returns:

List of turning point positions

action_integral(E: float, x1: float, x2: float) → float[source]

Compute action integral ∫p(x)dx between turning points.

Parameters:

E – Total energy
x1 – Left turning point
x2 – Right turning point

Returns:

Action integral value

bohr_sommerfeld_condition(E: float, x_range: Tuple[float, float]) → float[source]

Evaluate Bohr-Sommerfeld quantization condition.

For bound states: ∮p dx = (n + 1/2)h

Returns the value that should equal (n + 1/2) for valid energies.

Parameters:

E – Trial energy
x_range – Search range for turning points

Returns:

Quantization condition value

find_energy_level(n: int, E_range: Tuple[float, float], x_range: Tuple[float, float]) → float[source]

Find the n-th energy level using Bohr-Sommerfeld quantization.

Parameters:

n – Quantum number (0, 1, 2, …)
E_range – (E_min, E_max) search range
x_range – Spatial range for turning points

Returns:

Energy eigenvalue

bohr_sommerfeld_levels(n_max: int, E_range: Tuple[float, float], x_range: Tuple[float, float]) → List[EnergyLevel][source]

Compute multiple energy levels.

Parameters:

n_max – Maximum quantum number
E_range – Energy search range
x_range – Spatial range

Returns:

List of EnergyLevel objects

Example: Finding energy levels with WKB

from mechanics_dsl.domains.quantum import WKBApproximation

# Harmonic oscillator potential
def V(x):
    return 0.5 * x**2

wkb = WKBApproximation(potential=V, mass=1.0, hbar=1.0)

# Find ground state energy
E0 = wkb.find_energy_level(n=0, E_range=(0.1, 2.0), x_range=(-5, 5))
print(f"WKB ground state: E_0 = {E0:.3f}")  # ≈ 0.5

# Compare with exact: E_n = (n + 0.5)ℏω = 0.5

EhrenfestDynamics 

Ehrenfest theorem states that expectation values follow classical equations for quadratic potentials:

\[\frac{d\langle x \rangle}{dt} = \frac{\langle p \rangle}{m}\]

\[\frac{d\langle p \rangle}{dt} = -\left\langle \frac{dV}{dx} \right\rangle\]

class mechanics_dsl.domains.quantum.EhrenfestDynamics(potential: Callable[[float], float], potential_derivative: Callable[[float], float], mass: float = 1.0)[source]

Bases: object

Ehrenfest theorem: quantum-classical correspondence.

d<x>/dt = /m d/dt = -<dV/dx>

Expectation values follow classical equations for quadratic potentials.

classical_force(x: float) → float[source]: Classical force F = -dV/dx.

equations_of_motion(t: float, y: ndarray) → ndarray[source]

Classical equations for expectation values.

Parameters:

t – Time
y – State vector [<x>, ]

Returns:

Derivatives [d<x>/dt, d/dt]

propagate(x0: float, p0: float, t_span: Tuple[float, float], n_points: int = 100) → Dict[str, ndarray][source]

Propagate expectation values in time.

Parameters:

x0 – Initial <x>
p0 – Initial
t_span – (t_start, t_end)
n_points – Number of output points

Returns:

Dictionary with t, x_exp, p_exp arrays

Example: Wave packet evolution

from mechanics_dsl.domains.quantum import EhrenfestDynamics
import numpy as np

# Harmonic oscillator
def V(x): return 0.5 * x**2
def dV(x): return x

dynamics = EhrenfestDynamics(potential=V, potential_derivative=dV, mass=1.0)

# Initial conditions: x₀ = 1, p₀ = 0
result = dynamics.propagate(x0=1.0, p0=0.0, t_span=(0, 10))

# <x>(t) oscillates like classical oscillator
print(f"Initial <x> = {result['x_exp'][0]:.2f}")
print(f"Final <x> = {result['x_exp'][-1]:.2f}")

Physical Constants 

from mechanics_dsl.domains.quantum import HBAR, PLANCK_H

print(f"ℏ = {HBAR:.6e} J·s")        # 1.054571817e-34
print(f"h = {PLANCK_H:.6e} J·s")     # 6.62607015e-34

Convenience Functions 

from mechanics_dsl.domains.quantum import (
    de_broglie_wavelength,
    compton_wavelength,
    heisenberg_minimum
)

# de Broglie wavelength: λ = h/p = 2πℏ/p
lam = de_broglie_wavelength(momentum=1e-24)  # For electron-like particle

# Compton wavelength: λ_C = h/(mc)
lam_C = compton_wavelength(mass=9.11e-31)  # Electron Compton wavelength

# Heisenberg uncertainty minimum
# Heisenberg uncertainty minimum
min_uncertainty = heisenberg_minimum(hbar=1.0)  # ℏ/2

Quantum Tunneling 

The QuantumTunneling class provides exact and WKB-based tunneling calculations.

Rectangular Barrier (Exact)

For a barrier of height V₀ and width a:

\[T = \frac{1}{1 + \frac{V_0^2 \sinh^2(\kappa a)}{4E(V_0 - E)}}\]

where \(\kappa = \sqrt{2m(V_0 - E)}/\hbar\).

from mechanics_dsl.domains.quantum import QuantumTunneling

tunneling = QuantumTunneling(mass=1.0, hbar=1.0)

# Transmission through rectangular barrier
T = tunneling.rectangular_barrier(E=1.0, V0=2.0, width=1.0)
print(f"Transmission probability: {T:.4f}")

# WKB for arbitrary potentials
def gaussian_barrier(x):
    return 2.0 * np.exp(-x**2)

T_wkb = tunneling.wkb_transmission(E=0.5, potential=gaussian_barrier,
                                    x1=-2, x2=2)

Alpha Decay (Gamow Factor)

from mechanics_dsl.domains.quantum import alpha_decay_rate

# Uranium-238 alpha decay (E_alpha ≈ 4.2 MeV)
E_alpha = 4.2e6 * 1.602e-19  # Convert MeV to Joules
Z_thorium = 90  # Daughter nucleus

decay_rate = alpha_decay_rate(E_alpha, Z_daughter=Z_thorium)

Finite Square Well 

The FiniteSquareWell solves the transcendental eigenvalue equations:

\[\text{Even parity:} \quad k \tan(ka/2) = \kappa\]

\[\text{Odd parity:} \quad -k \cot(ka/2) = \kappa\]

where \(k = \sqrt{2m(E+V_0)}/\hbar\) and \(\kappa = \sqrt{-2mE}/\hbar\).

from mechanics_dsl.domains.quantum import FiniteSquareWell

well = FiniteSquareWell(depth=10.0, width=2.0, mass=1.0, hbar=1.0)

# Find all bound states
states = well.find_bound_states()
for state in states:
    print(f"n={state.n}: E = {state.energy:.4f}")

# Scattering transmission (E > 0)
T = well.transmission_coefficient(E=5.0)

Hydrogen Atom 

Exact energy levels for hydrogen-like atoms:

\[E_n = -\frac{Z^2 \times 13.6 \text{ eV}}{n^2}\]

from mechanics_dsl.domains.quantum import HydrogenAtom

H = HydrogenAtom(Z=1)

# Ground state energy
E1 = H.energy_level(n=1)  # -13.6 eV

# Bohr radius
r1 = H.bohr_radius_n(n=1)  # 5.29e-11 m

# Lyman-alpha wavelength (2→1)
lam_alpha = H.transition_wavelength(n_initial=2, n_final=1)  # 121.6 nm

# Balmer series (visible spectrum)
balmer = H.spectral_series(n_final=2, n_max=7)

# Helium ion (Z=2)
He_plus = HydrogenAtom(Z=2)
E1_He = He_plus.energy_level(n=1)  # -54.4 eV

Step Potential & Delta Barrier 

from mechanics_dsl.domains.quantum import StepPotential, DeltaFunctionBarrier

# Step potential
step = StepPotential(height=5.0)
R, T = step.reflection_transmission(E=10.0)
print(f"R = {R:.3f}, T = {T:.3f}")

# Delta function barrier: V(x) = λδ(x)
barrier = DeltaFunctionBarrier(strength=1.0)
T = barrier.transmission(E=0.5)

# Attractive delta well has one bound state
well = DeltaFunctionBarrier(strength=-1.0)
E_bound = well.bound_state_energy()  # -0.5