Appendix B: Python Quantum Simulation Toolkit

Over the course of 40 chapters, you built a quantum simulation toolkit piece by piece. This appendix consolidates every module into a single reference, documents the public API, and provides a complete workflow example. If you skipped the coding exercises and just want to get everything running, start at Section B.3 (Installation) and then jump to Section B.7 (Complete Workflow).

B.1 Module Overview

The Quantum Simulation Toolkit (qm_toolkit) is organized into the following modules. Each was introduced in a specific chapter and then extended in later chapters as new physics demanded new capabilities.

B.2 Core Classes

B.2.1 Ket

Represents a column vector in a finite-dimensional Hilbert space.

from qm_toolkit.core import Ket

# Create a ket from a list or numpy array
spin_up = Ket([1, 0], label="|↑⟩")
superposition = Ket([1/np.sqrt(2), 1/np.sqrt(2)], label="|+⟩")

# Properties
spin_up.dim          # 2
spin_up.components   # numpy array [1, 0]
spin_up.norm()       # 1.0
spin_up.normalized() # returns normalized copy

# Arithmetic
ket_sum = spin_up + Ket([0, 1])   # vector addition
scaled = 0.5 * spin_up            # scalar multiplication

# Inner product via Bra
bra_down = Bra([0, 1])
overlap = bra_down @ spin_up      # returns complex number 0+0j

# Outer product
projector = spin_up.outer(spin_up)  # returns Operator |↑⟩⟨↑|

# Tensor product
two_spins = spin_up.tensor(Ket([0, 1]))  # |↑⟩ ⊗ |↓⟩, dim=4

B.2.2 Bra

Represents a row vector (dual space). Constructed from a Ket or directly.

from qm_toolkit.core import Bra

bra = Bra([1, 0])
bra_from_ket = spin_up.bra()  # conjugate transpose

result = bra @ some_ket  # inner product (complex number)

B.2.3 Operator

Represents a linear operator on a finite-dimensional Hilbert space as a matrix.

from qm_toolkit.core import Operator
import numpy as np

# Create from a 2D numpy array
sigma_z = Operator(np.array([[1, 0], [0, -1]]), label="σ_z")

# Properties
sigma_z.dim                 # 2
sigma_z.matrix              # numpy array
sigma_z.is_hermitian()      # True
sigma_z.eigenvalues()       # array([1, -1])
sigma_z.eigenstates()       # list of Ket objects
sigma_z.trace()             # 0
sigma_z.det()               # -1

# Arithmetic
product = sigma_x @ sigma_z              # operator multiplication
commutator = sigma_x.commutator(sigma_z) # [σ_x, σ_z]
anticomm = sigma_x.anticommutator(sigma_z)

# Act on a ket
result_ket = sigma_z @ spin_up  # returns Ket with eigenvalue +1

# Expectation value
exp_val = spin_up.expectation(sigma_z)  # ⟨↑|σ_z|↑⟩ = 1.0

# Matrix exponential
U = (1j * theta * sigma_z / 2).expm()  # rotation operator

# Tensor product of operators
two_body_op = sigma_z.tensor(Operator.identity(2))  # σ_z ⊗ I

Built-in operator constructors:

Operator.identity(n)              # n×n identity
Operator.zero(n)                  # n×n zero matrix
Operator.pauli_x()                # σ_x
Operator.pauli_y()                # σ_y
Operator.pauli_z()                # σ_z
Operator.raising(j)               # J+ for spin-j
Operator.lowering(j)              # J- for spin-j
Operator.jx(j), .jy(j), .jz(j)   # angular momentum components
Operator.creation(n_max)          # a† truncated to n_max states
Operator.annihilation(n_max)      # a truncated to n_max states
Operator.number(n_max)            # N = a†a

B.2.4 QuantumState

A higher-level class that bundles a state vector with its basis information and supports both pure and mixed states.

from qm_toolkit.core import QuantumState

# Pure state
psi = QuantumState.pure(Ket([1, 0]), basis_labels=["|↑⟩", "|↓⟩"])

# Mixed state from density matrix
rho = QuantumState.mixed(density_matrix, basis_labels=["|↑⟩", "|↓⟩"])

# Mixed state from ensemble
rho = QuantumState.ensemble(
    states=[Ket([1,0]), Ket([0,1])],
    probabilities=[0.7, 0.3]
)

# Properties
psi.is_pure()                # True
psi.density_matrix()         # returns Operator (|ψ⟩⟨ψ| for pure states)
psi.von_neumann_entropy()    # 0 for pure states
psi.purity()                 # Tr(ρ²), 1 for pure states

# Measurement
outcome, post_state = psi.measure(sigma_z)  # simulates projective measurement
probs = psi.measurement_probabilities(sigma_z)  # dict of {eigenvalue: probability}

# Time evolution
psi_t = psi.evolve(hamiltonian, t=1.0, hbar=1.0)  # unitary evolution

B.2.5 Wavefunction

Represents a continuous wave function on a discrete grid.

from qm_toolkit.wavefunctions import Wavefunction

# Create from a function
x = np.linspace(-10, 10, 1000)
psi = Wavefunction(x, np.exp(-x**2 / 2), normalize=True)

# Properties
psi.x                       # position grid
psi.values                  # complex array of ψ(x)
psi.probability_density()   # |ψ(x)|²
psi.norm()                  # ∫|ψ|² dx (should be 1)
psi.expectation_x()         # ⟨x⟩
psi.expectation_x2()        # ⟨x²⟩
psi.uncertainty_x()         # Δx

# Momentum space
psi_k = psi.to_momentum_space()  # FFT with correct normalization
psi_k.expectation_p()             # ⟨p⟩
psi_k.uncertainty_p()             # Δp

# Verify uncertainty principle
print(psi.uncertainty_x() * psi_k.uncertainty_p())  # ≥ ℏ/2

# Overlap
overlap = psi.inner_product(phi)  # ⟨φ|ψ⟩

# Time evolution (split-operator method)
psi_evolved = psi.evolve(potential_func, dt=0.01, steps=100, mass=1.0)

B.3 Installation and Dependencies

Requirements

• Python 3.10 or later (3.11 recommended)
• NumPy >= 1.24
• SciPy >= 1.10
• Matplotlib >= 3.7
• QuTiP >= 5.0 (for advanced simulations and verification)
• SymPy >= 1.12 (for symbolic calculations in selected exercises)
• Jupyter >= 1.0 (for interactive notebooks)

Installing the Toolkit

If you have been building the toolkit chapter by chapter, the modules already live in your project directory under qm_toolkit/. To make them importable from anywhere:

# Navigate to the root of your textbook code directory
cd quantum-mechanics/code/

# Install in development mode (editable)
pip install -e .

This requires a minimal setup.py or pyproject.toml at quantum-mechanics/code/:

# pyproject.toml
[build-system]
requires = ["setuptools>=68.0"]
build-backend = "setuptools.backends._legacy:_Backend"

[project]
name = "qm_toolkit"
version = "1.0.0"
requires-python = ">=3.10"
dependencies = [
    "numpy>=1.24",
    "scipy>=1.10",
    "matplotlib>=3.7",
]

[project.optional-dependencies]
full = ["qutip>=5.0", "sympy>=1.12", "jupyter"]

Then install:

pip install -e ".[full]"

Quick Verification

Run the following to verify everything is working:

import numpy as np
from qm_toolkit.core import Ket, Operator
from qm_toolkit.constants import hbar, m_e, e_charge

# Spin-1/2 test
up = Ket([1, 0])
down = Ket([0, 1])
sz = Operator.pauli_z() * hbar / 2

print(f"⟨↑|σ_z|↑⟩ = {up.expectation(sz):.4e} J·s")
print(f"ℏ/2       = {hbar/2:.4e} J·s")

# Should match to machine precision
assert np.isclose(up.expectation(sz), hbar/2)
print("All checks passed.")

B.4 Key Functions Grouped by Topic

Stationary States and Eigenvalues

from qm_toolkit.solvers import (
    shooting_method,       # 1D eigenvalue solver via shooting (Ch 5)
    matrix_diag_solver,    # discretize Hamiltonian, diagonalize (Ch 5)
    dvr_solver,            # discrete variable representation (Ch 20)
)
from qm_toolkit.hydrogen import (
    hydrogen_radial_wf,    # R_nl(r) for hydrogen (Ch 10)
    hydrogen_energy,       # E_n = -13.6/n² eV (Ch 10)
    hydrogen_wf_3d,        # full ψ_nlm(r,θ,φ) (Ch 10)
)
from qm_toolkit.harmonic import (
    ho_eigenstate,         # ψ_n(x) for harmonic oscillator (Ch 6)
    ho_energy,             # E_n = ℏω(n + 1/2) (Ch 6)
    coherent_state,        # |α⟩ coherent state (Ch 9)
)

Angular Momentum and Spin

from qm_toolkit.angular import (
    spherical_harmonic,        # Y_l^m(θ, φ) (Ch 8)
    clebsch_gordan,            # ⟨j1 m1; j2 m2 | J M⟩ (Ch 14)
    wigner_d_matrix,           # D^j_{m'm}(α,β,γ) rotation matrix (Ch 14)
    couple_angular_momenta,    # returns coupled basis states (Ch 14)
)
from qm_toolkit.spin import (
    spinor,                    # create spin-1/2 state from angles (Ch 12)
    pauli_matrices,            # returns (σ_x, σ_y, σ_z) (Ch 12)
    spin_rotation,             # U(θ, n̂) rotation operator (Ch 13)
    stern_gerlach_simulate,    # simulate SG experiment (Ch 12)
    spin_operators,            # J_x, J_y, J_z for arbitrary j (Ch 14)
)

Perturbation Theory

from qm_toolkit.perturbation import (
    first_order_energy,        # E_n^(1) = ⟨n|H'|n⟩ (Ch 17)
    second_order_energy,       # E_n^(2) = Σ |⟨m|H'|n⟩|²/(E_n-E_m) (Ch 17)
    first_order_state,         # |n^(1)⟩ correction (Ch 17)
    degenerate_perturbation,   # diagonalize H' in degenerate subspace (Ch 18)
    time_dependent_amplitude,  # c_f(t) from Fermi's golden rule (Ch 19)
    transition_rate,           # Γ = (2π/ℏ)|⟨f|H'|i⟩|² ρ(E_f) (Ch 19)
)

Multi-Particle Systems

from qm_toolkit.multiparticle import (
    tensor_product_state,      # |ψ⟩ ⊗ |φ⟩ (Ch 20)
    tensor_product_operator,   # A ⊗ B (Ch 20)
    symmetrize,                # project onto symmetric subspace (Ch 21)
    antisymmetrize,            # project onto antisymmetric subspace (Ch 21)
    slater_determinant,        # build Slater determinant (Ch 22)
    partial_trace,             # Tr_B(ρ_AB) (Ch 30)
)

Scattering

from qm_toolkit.scattering import (
    phase_shift,               # δ_l from potential (Ch 27)
    partial_wave_amplitude,    # f_l(k) (Ch 27)
    differential_cross_section,# dσ/dΩ (Ch 27)
    total_cross_section,       # σ_tot via optical theorem (Ch 27)
    born_amplitude,            # first Born approximation (Ch 28)
    t_matrix_element,          # T-matrix from Lippmann-Schwinger (Ch 28)
)

Quantum Information

from qm_toolkit.entanglement import (
    bell_state,                # |Φ+⟩, |Φ-⟩, |Ψ+⟩, |Ψ-⟩ (Ch 31)
    concurrence,               # entanglement measure for 2 qubits (Ch 31)
    chsh_correlator,           # compute CHSH value (Ch 31)
    quantum_teleportation,     # simulate teleportation protocol (Ch 33)
)
from qm_toolkit.density import (
    von_neumann_entropy,       # S = -Tr(ρ ln ρ) (Ch 30)
    purity,                    # Tr(ρ²) (Ch 30)
    fidelity,                  # F(ρ, σ) (Ch 30)
    decohere,                  # apply decoherence channel (Ch 33)
)

Visualization

from qm_toolkit.visualization import (
    plot_wavefunction,         # plot ψ(x) and |ψ(x)|² (Ch 3)
    plot_potential_and_levels, # V(x) with energy eigenvalues overlaid (Ch 5)
    plot_probability_current,  # j(x,t) for wave packets (Ch 9)
    plot_spherical_harmonic,   # 3D surface plot of |Y_l^m|² (Ch 8)
    plot_hydrogen_orbital,     # 3D density plot of |ψ_nlm|² (Ch 10)
    plot_bloch_sphere,         # single qubit state on Bloch sphere (Ch 12)
    plot_wigner_function,      # Wigner quasi-probability W(x,p) (Ch 34)
    plot_energy_diagram,       # energy level diagram with transitions (Ch 17)
    animate_wavepacket,        # time evolution animation (Ch 9)
    plot_cross_section,        # dσ/dΩ polar plot (Ch 27)
)

B.5 Chapter-by-Chapter Module Index

This table tells you exactly which chapter introduced or extended each piece of the toolkit, so you can find the pedagogical context for any function.

B.6 QuTiP Interoperability

The toolkit provides seamless conversion to and from QuTiP objects for cross-validation and access to QuTiP's advanced solvers.

from qm_toolkit.core import Ket, Operator

# Toolkit -> QuTiP
ket = Ket([1, 0])
qutip_ket = ket.to_qutip()       # returns qutip.Qobj (ket type)

op = Operator.pauli_z()
qutip_op = op.to_qutip()         # returns qutip.Qobj (oper type)

# QuTiP -> Toolkit
import qutip
q_state = qutip.basis(2, 0)
toolkit_ket = Ket.from_qutip(q_state)

q_op = qutip.sigmax()
toolkit_op = Operator.from_qutip(q_op)

When to use QuTiP directly: - Master equation solvers (qutip.mesolve) for open quantum systems (Chapter 33) - Lindblad dynamics with multiple collapse operators - Large Hilbert spaces where QuTiP's sparse matrix optimizations shine - Bloch sphere animations via qutip.Bloch

When to use the toolkit: - Learning exercises where you want to see the math explicitly - Custom visualizations tuned to the textbook's plotting style - Small-system calculations where readability of code matters more than speed

B.7 Complete Workflow Example

Here is a complete example that starts from scratch, creates a quantum system, solves for eigenstates, applies perturbation theory, and visualizes the results. This example studies the Stark effect in hydrogen (linear electric field perturbation on the $n=2$ states).

import numpy as np
import matplotlib.pyplot as plt
from qm_toolkit.hydrogen import hydrogen_energy, hydrogen_wf_3d
from qm_toolkit.perturbation import degenerate_perturbation
from qm_toolkit.angular import spherical_harmonic
from qm_toolkit.core import Ket, Operator
from qm_toolkit.visualization import plot_energy_diagram
from qm_toolkit.constants import e_charge, a_0, hbar

# Step 1: Unperturbed n=2 energy
E2 = hydrogen_energy(n=2)
print(f"Unperturbed E_2 = {E2:.4f} eV")

# Step 2: Build the perturbation matrix H' = eEz in the n=2 subspace
# Basis: |2,0,0⟩, |2,1,0⟩, |2,1,1⟩, |2,1,-1⟩
# Only ⟨2,0,0|H'|2,1,0⟩ is nonzero (selection rule Δl=±1, Δm=0)

# The matrix element is -3eEa_0 (derived in Chapter 19)
E_field = 1e6  # V/m (modest lab field)
matrix_element = -3 * e_charge * E_field * a_0

H_prime = np.array([
    [0, matrix_element, 0, 0],
    [matrix_element, 0, 0, 0],
    [0, 0, 0, 0],
    [0, 0, 0, 0],
])

H_prime_op = Operator(H_prime, label="H'_Stark")

# Step 3: Diagonalize the perturbation in the degenerate subspace
energies, states = degenerate_perturbation(
    H_prime_op,
    unperturbed_energy=E2,
)

print("\nStark-shifted n=2 levels:")
for i, (e, s) in enumerate(zip(energies, states)):
    print(f"  State {i}: E = {e:.6f} eV")
    print(f"    Components: {s.components}")

# Step 4: Visualize
fig, ax = plot_energy_diagram(
    unperturbed_levels={"n=2": E2},
    perturbed_levels={"n=2": energies},
    title="Stark Effect in Hydrogen (n=2)",
)
plt.tight_layout()
plt.savefig("stark_effect_n2.png", dpi=150)
plt.show()

# Step 5: Verify — the splitting should be 2|⟨200|H'|210⟩|
splitting = abs(energies[0] - energies[1])
predicted = 2 * abs(matrix_element)
print(f"\nMeasured splitting:  {splitting:.6e} eV")
print(f"Predicted splitting: {predicted:.6e} eV")
assert np.isclose(splitting, predicted, rtol=1e-10)
print("Verification passed.")

Expected output:

Unperturbed E_2 = -3.4000 eV

Stark-shifted n=2 levels:
  State 0: E = -3.400024 eV
    Components: [ 0.707+0.j -0.707+0.j  0.   +0.j  0.   +0.j]
  State 1: E = -3.399976 eV
    Components: [ 0.707+0.j  0.707+0.j  0.   +0.j  0.   +0.j]
  State 2: E = -3.400000 eV
    Components: [ 0.   +0.j  0.   +0.j  1.   +0.j  0.   +0.j]
  State 3: E = -3.400000 eV
    Components: [ 0.   +0.j  0.   +0.j  0.   +0.j  1.   +0.j]

Measured splitting:  4.800000e-05 eV
Predicted splitting: 4.800000e-05 eV
Verification passed.

States 0 and 1 are the linear combinations $\frac{1}{\sqrt{2}}(|2,0,0\rangle \mp |2,1,0\rangle)$, which are the "parabolic" eigenstates. States 2 and 3 ($m = \pm 1$) are unaffected by the linear Stark effect because the selection rule $\Delta m = 0$ forbids mixing.

B.8 Coding Style Conventions

Throughout the textbook, all code follows these conventions:

1. SI units by default. All functions accept and return quantities in SI unless the docstring says otherwise. Use the units module for conversions.
2. NumPy arrays everywhere. All state vectors and matrices are stored as numpy.ndarray under the hood. The classes provide a friendlier interface but you can always access .components or .matrix directly.
3. Immutability. Ket, Bra, and Operator objects are immutable. Methods like normalized() return new objects rather than modifying in place.
4. Validation. Dimension mismatches raise ValueError with informative messages. Trying to multiply a 2-dimensional operator by a 3-dimensional ket tells you exactly what went wrong.
5. Labels are optional but encouraged. Pass label="..." to any constructor. Labels propagate through operations when possible and appear in __repr__ output and plot legends.
6. Plotting style. All visualization functions accept an optional ax parameter (a matplotlib Axes object). If not provided, they create a new figure. This makes it easy to compose multi-panel figures.
7. Random seeds. Any function that involves randomness (e.g., measure()) accepts an optional seed parameter for reproducibility in tests and grading.

B.9 Troubleshooting

The complete source code for every module, including docstrings and unit tests, lives in the code/qm_toolkit/ directory. Each module header includes the chapter of origin and a brief description. If you find a bug, the issue tracker for this textbook is the best place to report it.