ψ iℏ ∂ψ/∂t = Ĥψ

Quantum Mechanics

Heisenberg 1925 · Schrödinger 1926 · Dirac 1930 · Born Rule · Bell's Theorem

I. Historical Origins & Quantum Crises

By 1900, classical physics faced three catastrophic failures — phenomena it could not explain. Their resolution forced a complete reinvention of the laws of nature.

The Ultraviolet Catastrophe

Classical statistical mechanics predicted that a blackbody should radiate infinite energy at short wavelengths (the Rayleigh-Jeans law diverges). Planck (1900) resolved this by quantizing energy: radiation is emitted in discrete packets (quanta) of energy \(E = h\nu = \hbar\omega\):

$$B(\nu,T) = \frac{2h\nu^3}{c^2}\cdot\frac{1}{e^{h\nu/k_BT}-1}$$

Planck's constant: \(h = 6.626\times10^{-34}\) J·s, \(\hbar = h/(2\pi) = 1.055\times10^{-34}\) J·s.

Photoelectric Effect (Einstein, 1905)

Light ejects electrons from metals only above a threshold frequency, regardless of intensity. Einstein proposed that light consists of photons, each carrying energy \(E = h\nu\):

$$E_k = h\nu - \phi \qquad (\phi = \text{work function})$$

This earned Einstein the 1921 Nobel Prize and established the particle nature of light.

Bohr Model (1913)

Bohr postulated that electrons orbit the nucleus only on discrete circular orbits where angular momentum is quantized: \(L = n\hbar\). This gave the correct hydrogen spectrum:

$$E_n = -\frac{13.6\text{ eV}}{n^2}, \quad n = 1,2,3,\ldots$$

The model worked for hydrogen but failed for multi-electron atoms and had no theoretical justification.

YearDiscoveryPhysicistKey Idea
1900Blackbody quantizationPlanck\(E=h\nu\)
1905Photoelectric effectEinsteinPhotons
1913Bohr atomBohrQuantized orbits
1923Compton scatteringComptonPhoton momentum \(p=h/\lambda\)
1924Matter wavesde Broglie\(\lambda=h/p\)
1925Matrix mechanicsHeisenbergObservables as matrices
1926Wave mechanicsSchrödingerWave equation for \(\psi\)
1926Statistical interpretationBorn\(|\psi|^2 =\) probability density
1927Uncertainty principleHeisenberg\(\Delta x\Delta p\ge\hbar/2\)
1928Relativistic QMDiracDirac equation; predicts antimatter
1930Transformation theoryDiracBra-ket formalism; unifies SR/WM

II. Wave-Particle Duality

de Broglie Hypothesis (1924)

If light (wave) has particle properties, perhaps particles have wave properties. de Broglie proposed that every particle with momentum \(p\) has an associated wavelength:

$$\boxed{\lambda = \frac{h}{p} = \frac{h}{mv}\cdot\frac{1}{\gamma}} \qquad \text{(de Broglie wavelength)}$$

For an electron at 54 eV (Davisson-Germer experiment): \(\lambda \approx 1.67\) Å — same as crystal lattice spacing. Diffraction was observed. ✓

Double-Slit Experiment

A single electron (or photon, or atom, or even molecule) fired at a double-slit produces an interference pattern on the detector. Each particle goes through "both slits simultaneously." The probability amplitude \(\psi\) interferes; only \(|\psi|^2\) is observable.

Measurement Problem

If you detect which slit the particle passes through, the interference pattern disappears. The act of measurement irreversibly changes the quantum state. This is not a technical limitation — it is fundamental to quantum theory.

III. The Six Postulates of Quantum Mechanics

Postulate I — State Space

The state of a quantum system is completely described by a normalized vector \(\ket{\psi}\) (a "ket") in a complex Hilbert space \(\mathcal{H}\). The wavefunction \(\psi(x) = \braket{x}{\psi}\) is the position-space representation.

Postulate II — Observables

Every physical observable \(A\) corresponds to a Hermitian (self-adjoint) linear operator \(\op{A}\) on \(\mathcal{H}\). The only possible measurement outcomes are the eigenvalues of \(\op{A}\).

Postulate III — Born Rule

If the system is in state \(\ket{\psi}\) and we measure observable \(\op{A}\) with eigenstates \(\ket{a_n}\), the probability of obtaining result \(a_n\) is:

$$P(a_n) = |\braket{a_n}{\psi}|^2$$

For continuous observables (position): \(P(x,x+\d x) = |\psi(x)|^2\d x\).

Postulate IV — Collapse

Immediately after measuring \(a_n\), the state collapses to \(\ket{a_n}\). Subsequent measurements of \(A\) yield \(a_n\) with certainty.

Postulate V — Time Evolution (Schrödinger)

Between measurements, the state evolves unitarily:

$$\ii\hbar\frac{\partial}{\partial t}\ket{\psi} = \op{H}\ket{\psi}$$

where \(\op{H}\) is the Hamiltonian operator (total energy observable).

Postulate VI — Identical Particles

For a system of identical particles, the state must be symmetric (bosons, integer spin) or antisymmetric (fermions, half-integer spin) under exchange of any two particles.

IV. Hilbert Spaces

Definition — Hilbert Space

A Hilbert space \(\mathcal{H}\) is a complete inner product space over \(\mathbb{C}\). For quantum mechanics:

  • \(\mathcal{H} = L^2(\mathbb{R})\): square-integrable functions \(\int_{-\infty}^\infty|\psi(x)|^2\d x < \infty\)
  • Inner product: \(\braket{\phi}{\psi} = \int_{-\infty}^\infty \phi^*(x)\psi(x)\,\d x\)
  • Norm: \(\|\psi\| = \sqrt{\braket{\psi}{\psi}}\)
  • Normalization: \(\braket{\psi}{\psi} = 1\)

Key Properties of the Inner Product

$$\braket{\phi}{\psi} = \braket{\psi}{\phi}^*,\quad \braket{\phi}{\alpha\psi_1+\beta\psi_2} = \alpha\braket{\phi}{\psi_1}+\beta\braket{\phi}{\psi_2}$$

Completeness (Resolution of the Identity)

An orthonormal basis \(\{\ket{n}\}\) satisfies:

$$\braket{n}{m} = \delta_{nm}, \qquad \sum_n\ket{n}\bra{n} = \op{1}$$

Any state can be expanded: \(\ket{\psi} = \sum_n c_n\ket{n}\) where \(c_n = \braket{n}{\psi}\).

For continuous bases (e.g., position): \(\int\ket{x}\bra{x}\d x = \op{1}\), \(\braket{x}{x'} = \delta(x-x')\).

V. Dirac Bra-Ket Notation

Dirac's notation elegantly unifies the wave-mechanics and matrix-mechanics formulations.

SymbolNameMeaning
\(\ket{\psi}\)KetState vector in \(\mathcal{H}\)
\(\bra{\psi}\)BraDual vector (linear functional on \(\mathcal{H}\))
\(\braket{\phi}{\psi}\)Bracket (inner product)\(\int\phi^*\psi\,\d x\)
\(\ket{\psi}\bra{\phi}\)Outer productRank-1 operator: \((\ket{\psi}\bra{\phi})\ket{\chi}=\braket{\phi}{\chi}\ket{\psi}\)
\(\bra{\phi}\op{A}\ket{\psi}\)Matrix element\(\int\phi^*\op{A}\psi\,\d x\)
\(\psi(x) = \braket{x}{\psi}\)WavefunctionProjection onto position eigenstate
\(\tilde\psi(p) = \braket{p}{\psi}\)Momentum-space wavefunctionFourier transform of \(\psi(x)\)

The Fourier Connection

$$\psi(x) = \frac{1}{\sqrt{2\pi\hbar}}\int_{-\infty}^\infty \tilde\psi(p)\,e^{ipx/\hbar}\,\d p, \qquad \tilde\psi(p) = \frac{1}{\sqrt{2\pi\hbar}}\int_{-\infty}^\infty\psi(x)\,e^{-ipx/\hbar}\,\d x$$

VI. Operators, Observables & Hermitian Operators

Hermitian Operator

An operator \(\op{A}\) is Hermitian (self-adjoint) if \(\op{A} = \op{A}^\dagger\), i.e., \(\bra{\phi}\op{A}\ket{\psi} = \bra{\psi}\op{A}\ket{\phi}^* = \bra{\op{A}\phi}\ket{\psi}\) for all \(\ket{\phi},\ket{\psi}\in\mathcal{H}\).

Spectral Theorem for Hermitian Operators

Every Hermitian operator has (i) real eigenvalues and (ii) a complete orthonormal set of eigenstates.

Proof — Real Eigenvalues

Let \(\op{A}\ket{a} = a\ket{a}\). Then \(\bra{a}\op{A}\ket{a} = a\braket{a}{a} = a\). Also, by Hermiticity: \(\bra{a}\op{A}\ket{a} = (\bra{a}\op{A}\ket{a})^* = a^*\). So \(a = a^*\), hence \(a \in \mathbb{R}\).

Proof — Orthogonality of Eigenstates

Let \(\op{A}\ket{a}\!=\!a\ket{a}\) and \(\op{A}\ket{b}\!=\!b\ket{b}\). Then \(\bra{b}\op{A}\ket{a} = a\braket{b}{a}\). But also \(\bra{b}\op{A}\ket{a} = (\bra{a}\op{A}\ket{b})^* = (b\braket{a}{b})^* = b^*\braket{b}{a} = b\braket{b}{a}\) (since \(b\) is real). So \((a-b)\braket{b}{a} = 0\). For \(a\ne b\): \(\braket{b}{a} = 0\).

Fundamental Operators in QM

ObservableOperator (position space)Eigenvalues
Position\(\op{x} = x\cdot\)Any real \(x\)
Momentum\(\op{p} = -\ii\hbar\partial/\partial x\)Any real \(p\)
Kinetic energy\(\op{T} = -\hbar^2\nabla^2/(2m)\)\(p^2/2m \ge 0\)
Hamiltonian\(\op{H} = \op{T}+V(\op{x})\)Energy eigenvalues \(E_n\)
Angular momentum\(\op{L}_z = -\ii\hbar\partial/\partial\phi\)\(m\hbar\), \(m\in\mathbb{Z}\)
Parity\(\op{\Pi}\psi(x) = \psi(-x)\)\(\pm 1\)

Commutators

$$\comm{\op{A}}{\op{B}} \equiv \op{A}\op{B} - \op{B}\op{A}$$
Canonical Commutation Relation $$\boxed{\comm{\op{x}}{\op{p}} = \ii\hbar}$$
Proof

Acting on an arbitrary \(\psi(x)\):

$$\comm{\op{x}}{\op{p}}\psi = x\left(-\ii\hbar\pd{\psi}{x}\right) - \left(-\ii\hbar\right)\pd{(x\psi)}{x} = -\ii\hbar x\psi' + \ii\hbar(\psi + x\psi') = \ii\hbar\psi$$

Key commutators for angular momentum: \(\comm{L_i}{L_j} = \ii\hbar\varepsilon_{ijk}L_k\), \(\comm{L^2}{L_i} = 0\).

Unitary Time Evolution

The formal solution to the Schrödinger equation is:

$$\ket{\psi(t)} = e^{-\ii\op{H}t/\hbar}\ket{\psi(0)} \equiv \op{U}(t)\ket{\psi(0)}$$

where \(\op{U}(t)\) is unitary: \(\op{U}^\dagger\op{U} = \op{1}\) — preserves norms (probability is conserved).

VII. The Heisenberg Uncertainty Principle — Complete Proof

Robertson Uncertainty Relation

For any two Hermitian operators \(\op{A}\) and \(\op{B}\) in state \(\ket{\psi}\):

$$\boxed{\Delta A \cdot \Delta B \ge \frac{1}{2}\left|\expect{\comm{\op{A}}{\op{B}}}\right|}$$

where \(\Delta A = \sqrt{\expect{\op{A}^2} - \expect{\op{A}}^2}\) is the standard deviation.

Complete Proof (Robertson 1929)

Step 1. Define shifted operators: \(\delta\op{A} = \op{A} - \expect{A}\), \(\delta\op{B} = \op{B} - \expect{B}\). Then \((\Delta A)^2 = \expect{(\delta\op{A})^2}\) and similarly for \(B\).

Step 2. Define \(\ket{f} = \delta\op{A}\ket{\psi}\) and \(\ket{g} = \delta\op{B}\ket{\psi}\). By Cauchy-Schwarz:

$$(\Delta A)^2(\Delta B)^2 = \braket{f}{f}\braket{g}{g} \ge |\braket{f}{g}|^2$$

Step 3. Decompose \(\braket{f}{g}\) into symmetric and antisymmetric parts:

$$\braket{f}{g} = \underbrace{\frac{\braket{f}{g}+\braket{g}{f}}{2}}_{\equiv Z/2,\;\text{real}} + \underbrace{\frac{\braket{f}{g}-\braket{g}{f}}{2}}_{\equiv W/2,\;\text{imaginary}}$$ $$Z = \expect{\acom{\delta\op{A}}{\delta\op{B}}},\quad W = \expect{\comm{\delta\op{A}}{\delta\op{B}}} = \expect{\comm{\op{A}}{\op{B}}}$$

Step 4. Since \(\braket{f}{g} = Z/2 + W/2\), and \(W\) is purely imaginary (anti-Hermitian commutator):

$$|\braket{f}{g}|^2 = \frac{Z^2}{4} + \frac{|W|^2}{4} \ge \frac{|\expect{\comm{\op{A}}{\op{B}}}|^2}{4}$$

Step 5. Combining Steps 2 and 4:

$$(\Delta A)^2(\Delta B)^2 \ge \frac{1}{4}\left|\expect{\comm{\op{A}}{\op{B}}}\right|^2$$ $$\implies \boxed{\Delta A\cdot\Delta B \ge \frac{1}{2}\left|\expect{\comm{\op{A}}{\op{B}}}\right|}$$

Heisenberg Uncertainty for Position and Momentum

Since \(\comm{\op{x}}{\op{p}} = \ii\hbar\):

$$\boxed{\Delta x \cdot \Delta p \ge \frac{\hbar}{2}}$$

The minimum uncertainty state (Gaussian wave packet) saturates this inequality. Similarly, for energy and time: \(\Delta E \cdot \Delta t \ge \hbar/2\).

Physical Interpretation

The uncertainty principle is NOT about measurement disturbance — it is an intrinsic property of quantum states. A particle that has a definite position (Δx=0) must have completely undefined momentum, and vice versa. It is physically impossible for both Δx and Δp to be simultaneously zero.

VIII. Position and Momentum Representations

Position RepresentationMomentum Representation
State\(\psi(x) = \braket{x}{\psi}\)\(\tilde\psi(p) = \braket{p}{\psi}\)
\(\op{x}\)Multiply by \(x\)\(\ii\hbar\partial/\partial p\)
\(\op{p}\)\(-\ii\hbar\partial/\partial x\)Multiply by \(p\)
\(\op{H}\)\(-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}+V(x)\)\(\frac{p^2}{2m}+V(\ii\hbar\partial/\partial p)\)
Inner product\(\int\psi^*\phi\,\d x\)\(\int\tilde\psi^*\tilde\phi\,\d p\)
Normalization\(\int|\psi|^2\d x = 1\)\(\int|\tilde\psi|^2\d p = 1\)

The two representations are related by the Fourier transform, which is itself a unitary operation on \(L^2\).

IX. The Schrödinger Equation — Derivation & Analysis

Motivation from de Broglie and Planck

A free particle with momentum \(p\) and energy \(E\) has a plane-wave wavefunction:

$$\psi(x,t) = A\,e^{\ii(kx-\omega t)} = A\,e^{\ii(px-Et)/\hbar}$$

Acting with \(\op{p} = -\ii\hbar\partial/\partial x\) gives \(\op{p}\psi = p\psi\) ✓. Acting with \(\ii\hbar\partial/\partial t\) gives \(\ii\hbar\partial\psi/\partial t = E\psi\). The classical energy \(E = p^2/2m + V\) becomes the operator equation:

Time-Dependent Schrödinger Equation (TDSE) $$\boxed{\ii\hbar\frac{\partial\Psi}{\partial t} = \op{H}\Psi = \left[-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2} + V(x,t)\right]\Psi(x,t)}$$

Time-Independent Schrödinger Equation (TISE)

If \(V\) does not depend on time, separate variables: \(\Psi(x,t) = \psi(x)\,e^{-\ii Et/\hbar}\). The spatial part satisfies:

$$\boxed{\op{H}\psi(x) = E\psi(x) \quad\Leftrightarrow\quad -\frac{\hbar^2}{2m}\psi'' + V(x)\psi = E\psi}$$

General Solution by Superposition

$$\Psi(x,t) = \sum_n c_n\psi_n(x)\,e^{-\ii E_n t/\hbar}, \qquad c_n = \int_{-\infty}^\infty\psi_n^*(x)\Psi(x,0)\,\d x$$

X. Probability Density & Continuity Equation

Probability Density $$\rho(x,t) = |\Psi(x,t)|^2 = \Psi^*\Psi \ge 0, \qquad \int_{-\infty}^\infty\rho\,\d x = 1$$
Theorem — Conservation of Probability

If \(\Psi\) satisfies the TDSE, then \(\partial\rho/\partial t + \partial J/\partial x = 0\) (continuity equation), where:

$$J(x,t) = \frac{\hbar}{2m\ii}\left(\Psi^*\frac{\partial\Psi}{\partial x} - \Psi\frac{\partial\Psi^*}{\partial x}\right) = \frac{\hbar}{m}\,\mathrm{Im}\left(\Psi^*\frac{\partial\Psi}{\partial x}\right)$$
Proof

Compute \(\partial\rho/\partial t = \Psi^*\partial_t\Psi + \Psi\partial_t\Psi^*\). Use TDSE to replace \(\partial_t\Psi = (\ii\hbar/2m)\partial_{xx}\Psi - (\ii/\hbar)V\Psi\) and its complex conjugate:

$$\frac{\partial\rho}{\partial t} = \frac{\ii\hbar}{2m}\left(\Psi^*\Psi'' - \Psi\Psi^{*''}\right) = -\frac{\partial J}{\partial x}$$

where we recognized \(J = \frac{\hbar}{2m\ii}(\Psi^*\Psi' - \Psi\Psi^{*'})\).

XI. Expectation Values & Standard Deviations

$$\expect{\op{A}}_\psi = \bra{\psi}\op{A}\ket{\psi} = \int_{-\infty}^\infty\Psi^*(x,t)\,\op{A}\,\Psi(x,t)\,\d x$$

Key examples:

$$\expect{x} = \int x|\Psi|^2\d x, \qquad \expect{p} = \int\Psi^*\left(-\ii\hbar\pd{}{x}\right)\Psi\,\d x$$ $$\expect{H} = \int\Psi^*\left(-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}+V\right)\Psi\,\d x = E \text{ (for energy eigenstate)}$$

Ehrenfest's Theorem

Ehrenfest Theorem $$\frac{\d\expect{x}}{\d t} = \frac{\expect{p}}{m}, \qquad \frac{\d\expect{p}}{\d t} = -\left\langle\frac{\partial V}{\partial x}\right\rangle$$
Proof (first equation)

Compute \(\frac{\d}{\d t}\expect{x} = \int x\,\partial_t|\Psi|^2\d x\). Use the continuity equation \(\partial_t|\Psi|^2 = -\partial_x J\) and integrate by parts:

$$= -\int x\,\partial_x J\,\d x = \int J\,\d x = \frac{\hbar}{m}\mathrm{Im}\int\Psi^*\Psi'\,\d x = \frac{\expect{p}}{m}$$

Ehrenfest's theorem says the expectation values obey Newton's laws — QM reduces to CM for macroscopic objects where the wave packet is narrow.

XII. Particle in an Infinite Square Well

A particle confined to \(0 \le x \le L\) with \(V=0\) inside and \(V=\infty\) outside. Boundary conditions: \(\psi(0) = \psi(L) = 0\).

Solution

Inside: \(-(\hbar^2/2m)\psi'' = E\psi \implies \psi'' = -k^2\psi\), \(k = \sqrt{2mE}/\hbar\). General solution: \(\psi = A\sin kx + B\cos kx\). Boundary condition at \(x=0\): \(B=0\). At \(x=L\): \(A\sin kL = 0 \implies kL = n\pi\), \(n=1,2,3,\ldots\).

Particle in a Box — Exact Solution $$\psi_n(x) = \sqrt{\frac{2}{L}}\sin\!\left(\frac{n\pi x}{L}\right), \qquad E_n = \frac{n^2\pi^2\hbar^2}{2mL^2} = n^2 E_1$$ $$E_1 = \frac{\pi^2\hbar^2}{2mL^2} = \frac{h^2}{8mL^2} \qquad \text{(ground state energy)}$$

Key features: (1) Discrete energy levels (\(E \propto n^2\)); (2) Zero-point energy \(E_1 > 0\) — a consequence of the uncertainty principle; (3) \(n-1\) nodes in \(\psi_n\); (4) Orthonormality: \(\int_0^L\psi_m^*\psi_n\,\d x = \delta_{mn}\).

3D Infinite Square Well

$$E_{n_x n_y n_z} = \frac{\pi^2\hbar^2}{2m}\left(\frac{n_x^2}{L_x^2}+\frac{n_y^2}{L_y^2}+\frac{n_z^2}{L_z^2}\right),\quad \psi = \psi_{n_x}(x)\psi_{n_y}(y)\psi_{n_z}(z)$$

For a cube \(L_x=L_y=L_z\), the state \(E_{211}=E_{121}=E_{112}\) is 3-fold degenerate.

XIII. Quantum Harmonic Oscillator — Analytic Method

Potential: \(V(x) = \frac{1}{2}m\omega^2 x^2\). TISE:

$$-\frac{\hbar^2}{2m}\psi'' + \frac{1}{2}m\omega^2 x^2\psi = E\psi$$

Solution via Hermite Polynomials

Introduce dimensionless variable \(\xi = x\sqrt{m\omega/\hbar}\). The normalizable solutions are:

QHO Eigenstates $$\psi_n(x) = \left(\frac{m\omega}{\pi\hbar}\right)^{1/4}\frac{1}{\sqrt{2^n n!}}H_n(\xi)\,e^{-\xi^2/2}, \quad E_n = \hbar\omega\!\left(n+\tfrac{1}{2}\right), \quad n=0,1,2,\ldots$$

Hermite polynomials \(H_n(\xi) = (-1)^n e^{\xi^2}\frac{\d^n}{\d\xi^n}e^{-\xi^2}\):

$$H_0=1,\;H_1=2\xi,\;H_2=4\xi^2-2,\;H_3=8\xi^3-12\xi,\;H_4=16\xi^4-48\xi^2+12,\ldots$$

Ladder Operators (Algebraic Method)

Define:

$$\op{a} = \sqrt{\frac{m\omega}{2\hbar}}\left(\op{x}+\frac{\ii}{\:m\omega}\op{p}\right), \qquad \op{a}^\dagger = \sqrt{\frac{m\omega}{2\hbar}}\left(\op{x}-\frac{\ii}{\:m\omega}\op{p}\right)$$
Ladder Operator Properties $$\comm{\op{a}}{\op{a}^\dagger} = 1, \qquad \op{H} = \hbar\omega\!\left(\op{a}^\dagger\op{a}+\tfrac{1}{2}\right) = \hbar\omega\!\left(\op{N}+\tfrac{1}{2}\right)$$ $$\op{a}\ket{n} = \sqrt{n}\,\ket{n-1}, \qquad \op{a}^\dagger\ket{n} = \sqrt{n+1}\,\ket{n+1}$$ $$\op{N}\ket{n} = n\ket{n} \qquad\text{(number operator)}$$
Proof — Hamiltonian in terms of ladder operators

Compute \(\op{a}^\dagger\op{a} = \frac{m\omega}{2\hbar}\left(\op{x}^2 + \frac{\op{p}^2}{m^2\omega^2} + \frac{\ii}{m\omega}[\op{x},\op{p}]\right) = \frac{m\omega}{2\hbar}\left(\op{x}^2 + \frac{\op{p}^2}{m^2\omega^2}\right) - \frac{1}{2}\)

Therefore \(\op{H} = \frac{\op{p}^2}{2m}+\frac{1}{2}m\omega^2\op{x}^2 = \hbar\omega\left(\op{a}^\dagger\op{a}+\frac{1}{2}\right)\).

Proof — Energy Ladder (all eigenvalues are (n+½)ℏω)

If \(\op{H}\ket{E} = E\ket{E}\), consider \(\op{a}\ket{E}\). Using \(\op{H}\op{a} = \op{a}(\op{H}-\hbar\omega)\):

$$\op{H}(\op{a}\ket{E}) = (E-\hbar\omega)(\op{a}\ket{E})$$

So \(\op{a}\ket{E}\) is an eigenstate with energy \(E-\hbar\omega\). Repeated lowering must terminate since \(\braket{n}{\op{N}}{n} = n\braket{n}{n} \ge 0\). Thus there exists a ground state \(\ket{0}\) with \(\op{a}\ket{0} = 0\). Solving: \(\op{H}\ket{0} = \frac{1}{2}\hbar\omega\ket{0}\). Then \(E_n = \hbar\omega(n+\frac{1}{2})\).

Matrix Elements

$$\op{x} = \sqrt{\frac{\hbar}{2m\omega}}(\op{a}+\op{a}^\dagger), \qquad \op{p} = \ii\sqrt{\frac{m\omega\hbar}{2}}(\op{a}^\dagger-\op{a})$$ $$\bra{m}\op{x}\ket{n} = \sqrt{\frac{\hbar}{2m\omega}}\left(\sqrt{n}\,\delta_{m,n-1}+\sqrt{n+1}\,\delta_{m,n+1}\right)$$

XIV. Free Particle & Gaussian Wave Packets

For \(V=0\), energy eigenstates are plane waves \(\psi_k = e^{ikx}\), \(E = \hbar^2k^2/2m\). These are not normalizable (not in \(L^2\)) — the physical state is a wave packet:

$$\Psi(x,t) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty\phi(k)\,e^{\ii(kx-\omega(k)t)}\,\d k, \qquad \omega(k) = \frac{\hbar k^2}{2m}$$

Gaussian Wave Packet

Initial state \(\Psi(x,0) = (2\pi\sigma^2)^{-1/4}\exp\!\left(-\frac{x^2}{4\sigma^2}+\ii k_0 x\right)\):

$$\phi(k) = \left(\frac{2\sigma^2}{\pi}\right)^{1/4}\exp\!\left(-\sigma^2(k-k_0)^2\right)$$

Time evolution (exact for free particle):

$$\Psi(x,t) = \left(\frac{1}{2\pi\sigma^2(t)^2}\right)^{1/4}\exp\!\left(-\frac{(x-v_g t)^2}{4\sigma^2(t)}\right)\cdot\text{(phase)}$$ $$\sigma(t) = \sigma\sqrt{1+\left(\frac{\hbar t}{2m\sigma^2}\right)^2}\quad \text{(spreading width)}$$

Group velocity: \(v_g = \hbar k_0/m\). Phase velocity: \(v_\phi = \hbar k_0/(2m) = v_g/2\). The packet spreads as \(\sigma(t)\sim \hbar t/(2m\sigma)\) for large \(t\).

XV. Quantum Tunneling & the WKB Method

Rectangular Barrier

A particle of energy \(E < V_0\) hitting a barrier of height \(V_0\) and width \(a\):

$$T = \left[1 + \frac{V_0^2\sinh^2(\kappa a)}{4E(V_0-E)}\right]^{-1}, \qquad \kappa = \frac{\sqrt{2m(V_0-E)}}{\hbar}$$

For thick barriers \(\kappa a \gg 1\): \(T \approx 16\frac{E(V_0-E)}{V_0^2}e^{-2\kappa a}\).

WKB Approximation

For a slowly varying potential, the wavefunction takes the semiclassical form. The tunneling transmission coefficient through a general barrier from \(x_1\) to \(x_2\) is:

$$\boxed{T \approx \exp\!\left(-\frac{2}{\hbar}\int_{x_1}^{x_2}\sqrt{2m(V(x)-E)}\,\d x\right)}$$

Applications: alpha decay (Gamow factor), Josephson junctions, scanning tunneling microscopy, Fowler-Nordheim field emission.

Gamow Factor for Alpha Decay

$$\Gamma \propto e^{-G}, \qquad G = \frac{2}{\hbar}\int_{R}^{r_0}\sqrt{2m_\alpha\!\left(\frac{2Ze^2}{4\pi\varepsilon_0 r}-E\right)}\,\d r \approx \frac{4Ze^2}{\hbar v}\left(\arccos\sqrt{\frac{E}{V_0}}-\sqrt{\frac{E}{V_0}\left(1-\frac{E}{V_0}\right)}\right)$$

XVI. Orbital Angular Momentum

In 3D: \(\bv{L} = \bv{r}\times\bv{p}\), components \(L_x = yp_z-zp_y\), etc. In operator form:

$$L_z = -\ii\hbar\frac{\partial}{\partial\phi}, \qquad L^2 = -\hbar^2\!\left[\frac{1}{\sin\theta}\pd{}{\theta}\!\left(\sin\theta\pd{}{\theta}\right)+\frac{1}{\sin^2\theta}\frac{\partial^2}{\partial\phi^2}\right]$$
Angular Momentum Eigenvalue Equations $$L^2\ket{l,m} = l(l+1)\hbar^2\ket{l,m}, \quad l = 0,1,2,\ldots$$ $$L_z\ket{l,m} = m\hbar\ket{l,m}, \quad m = -l,-l+1,\ldots,+l$$

Spherical Harmonics \(Y_l^m(\theta,\phi)\)

The simultaneous eigenfunctions of \(L^2\) and \(L_z\):

$$Y_l^m(\theta,\phi) = \varepsilon\sqrt{\frac{(2l+1)}{4\pi}\frac{(l-|m|)!}{(l+|m|)!}}\,P_l^{|m|}(\cos\theta)\,e^{\ii m\phi}$$
\(l\)\(m\)Name\(Y_l^m(\theta,\phi)\)
00s\(\frac{1}{2\sqrt{\pi}}\)
10p\(_z\)\(\frac{1}{2}\sqrt{\frac{3}{\pi}}\cos\theta\)
1±1p\(_x\),p\(_y\)\(\mp\frac{1}{2}\sqrt{\frac{3}{2\pi}}\sin\theta\,e^{\pm\ii\phi}\)
20d\(_{z^2}\)\(\frac{1}{4}\sqrt{\frac{5}{\pi}}(3\cos^2\theta-1)\)
2±1d\(_{xz}\),d\(_{yz}\)\(\mp\frac{1}{2}\sqrt{\frac{15}{2\pi}}\sin\theta\cos\theta\,e^{\pm\ii\phi}\)
2±2d\(_{xy}\),d\(_{x^2-y^2}\)\(\frac{1}{4}\sqrt{\frac{15}{2\pi}}\sin^2\theta\,e^{\pm 2\ii\phi}\)

XVII. Spin-½ & Pauli Matrices

Spin is an intrinsic angular momentum with no classical analog. For spin-½ (electrons, protons, neutrons): \(s = 1/2\), \(m_s = \pm 1/2\).

Pauli Matrices $$\sigma_x = \begin{pmatrix}0&1\\1&0\end{pmatrix},\quad \sigma_y = \begin{pmatrix}0&-\ii\\\ii&0\end{pmatrix},\quad \sigma_z = \begin{pmatrix}1&0\\0&-1\end{pmatrix}$$ $$\op{S} = \frac{\hbar}{2}\bv{\sigma}, \qquad S^2 = \frac{3\hbar^2}{4}\op{I} = \hbar^2 s(s+1)\op{I}$$

Pauli Matrix Algebra

$$\sigma_i\sigma_j = \delta_{ij}\op{I} + \ii\varepsilon_{ijk}\sigma_k$$ $$\comm{\sigma_i}{\sigma_j} = 2\ii\varepsilon_{ijk}\sigma_k, \qquad \acom{\sigma_i}{\sigma_j} = 2\delta_{ij}\op{I}$$

Spin-up and spin-down eigenstates of \(S_z\):

$$\ket{\uparrow} = \ket{+} = \begin{pmatrix}1\\0\end{pmatrix},\quad \ket{\downarrow} = \ket{-} = \begin{pmatrix}0\\1\end{pmatrix}$$

Spin in an Arbitrary Direction

Spin operator along unit vector \(\hat{n} = (\sin\theta\cos\phi, \sin\theta\sin\phi, \cos\theta)\):

$$\bv{S}\cdot\hat{n} = \frac{\hbar}{2}\begin{pmatrix}\cos\theta & e^{-\ii\phi}\sin\theta \\ e^{\ii\phi}\sin\theta & -\cos\theta\end{pmatrix}$$

Eigenvalues: \(\pm\hbar/2\). Eigenstates: \(\ket{+,\hat{n}} = \cos(\theta/2)\ket{+}+e^{\ii\phi}\sin(\theta/2)\ket{-}\).

Addition of Angular Momenta — Clebsch-Gordan

Two systems with angular momenta \(j_1\) and \(j_2\). Total angular momentum \(J = j_1+j_2\) takes values \(|j_1-j_2|\le J\le j_1+j_2\) (integer steps). The coupled basis is related to the uncoupled basis by:

$$\ket{J,M} = \sum_{m_1,m_2}\braket{j_1,m_1;j_2,m_2}{J,M}\ket{j_1,m_1}\ket{j_2,m_2}$$

The coefficients \(\braket{j_1,m_1;j_2,m_2}{J,M}\) are the Clebsch-Gordan coefficients. For two spin-½ particles:

$$\ket{1,1} = \ket{\uparrow\uparrow},\quad \ket{1,0} = \frac{1}{\sqrt{2}}(\ket{\uparrow\downarrow}+\ket{\downarrow\uparrow}),\quad \ket{1,-1} = \ket{\downarrow\downarrow}$$ $$\ket{0,0} = \frac{1}{\sqrt{2}}(\ket{\uparrow\downarrow}-\ket{\downarrow\uparrow}) \qquad\text{(singlet state)}$$

XVIII. The Hydrogen Atom

The Coulomb potential \(V(r) = -e^2/(4\pi\varepsilon_0 r)\) in 3D gives the TISE:

$$-\frac{\hbar^2}{2\mu}\nabla^2\psi - \frac{e^2}{4\pi\varepsilon_0 r}\psi = E\psi$$

In spherical coordinates, separate: \(\psi_{nlm}(r,\theta,\phi) = R_{nl}(r)\,Y_l^m(\theta,\phi)\).

Radial Equation

Substituting \(R_{nl} = u_{nl}(r)/r\):

$$-\frac{\hbar^2}{2\mu}\frac{\d^2 u}{\d r^2} + \underbrace{\left[\frac{\hbar^2 l(l+1)}{2\mu r^2} - \frac{e^2}{4\pi\varepsilon_0 r}\right]}_{V_\text{eff}(r)} u = E\,u$$

Introduce the Bohr radius \(a_0 = 4\pi\varepsilon_0\hbar^2/(\mu e^2) = 0.529\) Å and \(\rho = 2r/(na_0)\). The normalizable solutions involve associated Laguerre polynomials:

$$R_{nl}(r) = -\sqrt{\left(\frac{2}{na_0}\right)^3\frac{(n-l-1)!}{2n[(n+l)!]^3}}\,e^{-\rho/2}\rho^l\,L_{n-l-1}^{2l+1}(\rho)$$

Energy Spectrum & Quantum Numbers

Hydrogen Atom Energy Levels $$E_n = -\frac{\mu e^4}{2(4\pi\varepsilon_0)^2\hbar^2}\cdot\frac{1}{n^2} = -\frac{13.606\text{ eV}}{n^2}, \quad n = 1,2,3,\ldots$$

Quantum number constraints: \(n \ge 1\), \(0 \le l \le n-1\), \(-l \le m \le l\). Degeneracy of level \(n\): \(g_n = n^2\) (ignoring spin), \(2n^2\) (with spin).

SeriesLower levelTransitionsSpectral region
Lymann=1n=2,3,4,…→1UV (λ = 91–122 nm)
Balmern=2n=3,4,5,…→2Visible (λ = 365–656 nm)
Paschenn=3n=4,5,6,…→3IR (λ = 820–1875 nm)
Brackettn=4n=5,6,7,…→4IR

Selection Rules (Electric Dipole)

$$\Delta l = \pm 1, \qquad \Delta m = 0, \pm 1, \qquad \Delta n = \text{any}$$

These arise from the parity of spherical harmonics and conservation of angular momentum (photon carries \(j=1\)).

XIX. Time-Independent Perturbation Theory

Write \(\op{H} = \op{H}^{(0)} + \lambda\op{H}'\) where \(\op{H}^{(0)}\) is exactly solvable. Expand: \(E_n = E_n^{(0)} + \lambda E_n^{(1)} + \lambda^2 E_n^{(2)} + \ldots\)

First & Second Order Corrections $$E_n^{(1)} = \bra{n^{(0)}}\op{H}'\ket{n^{(0)}}$$ $$E_n^{(2)} = \sum_{m\ne n}\frac{|\bra{m^{(0)}}\op{H}'\ket{n^{(0)}}|^2}{E_n^{(0)}-E_m^{(0)}}$$ $$\ket{n^{(1)}} = \sum_{m\ne n}\frac{\bra{m^{(0)}}\op{H}'\ket{n^{(0)}}}{E_n^{(0)}-E_m^{(0)}}\ket{m^{(0)}}$$
Derivation (First Order)

Write the eigenvalue equation \((\op{H}^{(0)}+\lambda\op{H}')(|n^{(0)}\rangle+\lambda|n^{(1)}\rangle+\ldots) = (E_n^{(0)}+\lambda E_n^{(1)}+\ldots)(|n^{(0)}\rangle+\ldots)\). Collect terms at order \(\lambda^1\):

$$\op{H}^{(0)}\ket{n^{(1)}}+\op{H}'\ket{n^{(0)}} = E_n^{(0)}\ket{n^{(1)}}+E_n^{(1)}\ket{n^{(0)}}$$

Project onto \(\bra{n^{(0)}}\): \(\bra{n^{(0)}}\op{H}'\ket{n^{(0)}} = E_n^{(1)}\) (using \(\bra{n^{(0)}}\op{H}^{(0)}\ket{n^{(1)}} = E_n^{(0)}\braket{n^{(0)}}{n^{(1)}}\) and \(\braket{n^{(0)}}{n^{(1)}} = 0\) by normalization).

Example: Anharmonic Oscillator

For \(\op{H}' = \lambda x^4\), using \(x = \sqrt{\hbar/(2m\omega)}(\op{a}+\op{a}^\dagger)\):

$$E_n^{(1)} = 3\lambda\left(\frac{\hbar}{2m\omega}\right)^2(2n^2+2n+1)$$

XX. Variational Principle

Variational Theorem

For any normalized trial state \(\ket{\psi_\text{trial}}\):

$$E_\text{gs} \le \bra{\psi_\text{trial}}\op{H}\ket{\psi_\text{trial}} \equiv E[\psi_\text{trial}]$$
Proof

Expand in energy eigenstates: \(\ket{\psi} = \sum_n c_n\ket{n}\). Then \(\bra{\psi}\op{H}\ket{\psi} = \sum_n|c_n|^2 E_n \ge E_\text{gs}\sum_n|c_n|^2 = E_\text{gs}\cdot 1 = E_\text{gs}\).

Helium atom: With trial wavefunction \(\psi(r_1,r_2) = (Z_\text{eff}^3/\pi a_0^3)e^{-Z_\text{eff}(r_1+r_2)/a_0}\), optimize over \(Z_\text{eff}\): \(Z_\text{eff} = 27/16 = 1.6875\) (screening). Gives \(E_\text{gs} \approx -77.5\) eV vs exact \(-79.0\) eV (0.9% error).

XXI. Time-Dependent Perturbation Theory & Fermi's Golden Rule

For a sinusoidal perturbation \(\op{H}'(t) = \op{V}\,e^{-\ii\omega t} + \op{V}^\dagger\,e^{+\ii\omega t}\) applied to state \(\ket{i}\), the first-order transition rate to state \(\ket{f}\) (Fermi's Golden Rule):

Fermi's Golden Rule $$\boxed{\Gamma_{i\to f} = \frac{2\pi}{\hbar}\left|\bra{f}\op{V}\ket{i}\right|^2\rho(E_f)}$$

where \(\rho(E_f)\) is the density of final states. Applications: spontaneous emission, Auger effect, beta decay, scattering cross-sections.

XXII. Identical Particles & the Pauli Exclusion Principle

For two identical particles, swapping them cannot change any observable. Define the exchange operator \(\op{P}_{12}\). Then \(\op{P}_{12}^2 = \op{1}\), so eigenvalues are \(\pm 1\).

Spin-Statistics Theorem (Pauli, 1940)

Bosons (integer spin): \(\psi(\bv{r}_2,\bv{r}_1) = +\psi(\bv{r}_1,\bv{r}_2)\) — symmetric state.

Fermions (half-integer spin): \(\psi(\bv{r}_2,\bv{r}_1) = -\psi(\bv{r}_1,\bv{r}_2)\) — antisymmetric state.

Pauli Exclusion Principle

Two identical fermions cannot occupy the same quantum state. If we attempt to put both in state \(\phi_a\), the antisymmetric combination \(\psi = \frac{1}{\sqrt{2}}[\phi_a(\bv{r}_1)\phi_a(\bv{r}_2) - \phi_a(\bv{r}_2)\phi_a(\bv{r}_1)] = 0\).

Slater determinant for \(N\) fermions in states \(\phi_1,\ldots,\phi_N\):

$$\Psi = \frac{1}{\sqrt{N!}}\begin{vmatrix}\phi_1(\bv{r}_1) & \phi_1(\bv{r}_2) & \cdots & \phi_1(\bv{r}_N) \\ \phi_2(\bv{r}_1) & \phi_2(\bv{r}_2) & \cdots & \phi_2(\bv{r}_N) \\ \vdots & & \ddots & \vdots \\ \phi_N(\bv{r}_1) & \cdots & & \phi_N(\bv{r}_N)\end{vmatrix}$$

XXIII. Quantum Entanglement & Bell's Theorem

EPR Paradox (1935)

Einstein, Podolsky, and Rosen considered the singlet state \(\ket{\Psi^-} = \frac{1}{\sqrt{2}}(\ket{\uparrow\downarrow}-\ket{\downarrow\uparrow})\). Measuring spin of particle 1 along \(z\) instantly determines the spin of particle 2, no matter how far apart they are. EPR argued QM must be incomplete — there should be hidden variables predetermining outcomes.

Bell's Theorem (1964)

Bell proved that any local hidden variable theory must satisfy the inequality:

$$|E(a,b) - E(a,c)| \le 1 + E(b,c)$$

where \(E(a,b)\) is the correlation of spin measurements along directions \(a\) and \(b\). Quantum mechanics predicts \(E(a,b) = -\cos(a-b)\), which violates Bell's inequality for certain angle choices. Experiments (Aspect 1982, Hensen et al. 2015) conclusively confirm QM and rule out local hidden variables.

Entanglement is Not Faster-Than-Light Communication

Although measuring particle 1 instantly correlates with particle 2, no information can be transmitted this way — the outcome of measurement 1 is random, and particle 2's owner cannot determine what measurement was made on particle 1 until they receive a classical signal.

XXIV. Density Matrix & Mixed States

For a pure state \(\ket{\psi}\): \(\op{\rho} = \ket{\psi}\bra{\psi}\). For a mixed state (statistical ensemble with probability \(p_i\) for state \(\ket{\psi_i}\)):

$$\op{\rho} = \sum_i p_i\ket{\psi_i}\bra{\psi_i}, \qquad \sum_i p_i = 1$$

Properties: \(\mathrm{Tr}(\op{\rho}) = 1\), \(\op{\rho}^\dagger = \op{\rho}\), \(\op{\rho}\ge 0\). For pure states: \(\op{\rho}^2 = \op{\rho}\) (\(\mathrm{Tr}(\op\rho^2)=1\)). For mixed: \(\mathrm{Tr}(\op\rho^2)<1\).

Von Neumann entropy: \(S = -\mathrm{Tr}(\op\rho\ln\op\rho) = -\sum_i\lambda_i\ln\lambda_i\).

Time evolution: \(\ii\hbar\,\partial_t\op\rho = [\op{H},\op\rho]\) (Liouville–von Neumann equation).

XXV. GNU Octave Examples

1 — Particle in a Box: Wavefunctions, Energies & Expectation Values

Infinite Square WellGNU Octave
%% Particle in an infinite square well [0,L]
%% Computes eigenstates, energies, expectation values, and uncertainty
clear; clc;

%% Constants (SI)
hbar = 1.0546e-34;  m = 9.109e-31;  L = 1e-9;  % electron, 1 nm box
E1 = pi^2*hbar^2/(2*m*L^2);   % ground state energy in Joules
printf('Ground state energy E1 = %.4f eV\n', E1/1.602e-19);

%% Wavefunctions and probability densities
x = linspace(0, L, 1000);
figure(1); clf;

for n = 1:4
  psi_n = sqrt(2/L) .* sin(n*pi*x/L);
  E_n   = n^2 * E1 / 1.602e-19;       % energy in eV

  %% Expectation values via numerical integration
  dx = x(2)-x(1);
  x_exp  = trapz(x, x .* psi_n.^2);
  x2_exp = trapz(x, x.^2 .* psi_n.^2);
  % ⟨p⟩ = 0 by symmetry; ⟨p²⟩ = n²π²ℏ²/L²
  p2_exp = n^2*pi^2*hbar^2/L^2;
  Delta_x = sqrt(x2_exp - x_exp^2);
  Delta_p = sqrt(p2_exp);           % since ⟨p⟩=0
  sigma_xp= Delta_x * Delta_p / hbar; % should be ≥ 1/2

  printf('n=%d: E=%.3f eV, ⟨x⟩=%.3f nm, Δx·Δp/ℏ=%.4f ≥ 0.5\n', ...
         n, E_n, x_exp/1e-9, sigma_xp);

  subplot(2,4,n);
  plot(x/1e-9, psi_n*sqrt(1e-9), 'b-', 'LineWidth', 1.8); grid on;
  title(sprintf('ψ_%d,  E=%.2f eV',n,E_n));
  xlabel('x (nm)'); ylabel('ψ_n (nm^{-1/2})');

  subplot(2,4,n+4);
  plot(x/1e-9, psi_n.^2*1e-9, 'r-', 'LineWidth', 1.8); grid on;
  title(sprintf('|ψ_%d|²',n));
  xlabel('x (nm)'); ylabel('|ψ|² (nm^{-1})');
end

%% Time evolution of superposition |ψ⟩ = (|1⟩ + |2⟩)/√2
omega1 = E1/hbar; omega2 = 4*E1/hbar;
T_beat = 2*pi/(omega2-omega1);   % beat period
t_arr = linspace(0, T_beat, 100);
figure(2);
psi1 = sqrt(2/L)*sin(pi*x/L);
psi2 = sqrt(2/L)*sin(2*pi*x/L);
for k = 1:5:50
  t = t_arr(k);
  PSI = (1/sqrt(2))*(psi1*exp(-1i*omega1*t) + psi2*exp(-1i*omega2*t));
  plot(x/1e-9, abs(PSI).^2*1e-9); hold on;
end
title('Time evolution of |ψ⟩=(|1⟩+|2⟩)/√2 — probability density');
xlabel('x (nm)'); ylabel('|Ψ|² (nm^{-1})'); grid on;

2 — Finite Difference Method: General 1D Schrödinger Solver

Finite Difference Schrödinger SolverGNU Octave
%% Numerically solve the 1D TISE using the finite difference method
%% H = T + V,  T_ij = -ℏ²/(2m dx²) * (δ_{i,j-1} - 2δ_{ij} + δ_{i,j+1})
%% Works for any potential V(x)

function [E, psi] = solve_schrodinger(x, V, m, hbar)
  N = length(x);
  dx = x(2) - x(1);
  coeff = -hbar^2/(2*m*dx^2);

  %% Build Hamiltonian as sparse tridiagonal matrix
  diag_main = -2*coeff*ones(N,1) + V(:);
  diag_off  =    coeff*ones(N-1,1);
  H = diag(diag_main) + diag(diag_off,1) + diag(diag_off,-1);

  %% Diagonalise (eigenvalues in ascending order)
  [psi, Emat] = eig(H);
  E = diag(Emat);
  [E, idx] = sort(real(E));
  psi = real(psi(:, idx));

  %% Normalise each eigenstate
  for k = 1:size(psi,2)
    norm_k = sqrt(trapz(x, psi(:,k).^2));
    psi(:,k) /= norm_k;
  end
end

clear; clc;
hbar = 1.0546e-34;  m_e = 9.109e-31;

%% ─── Example 1: Infinite square well (verify analytic) ───────
L = 1e-9;  N = 500;
x1 = linspace(0, L, N);
V1 = zeros(1,N);  V1([1 end]) = 1e10;  % hard walls
[E1, psi1] = solve_schrodinger(x1, V1, m_e, hbar);
E_analytic = (((1:5).^2)*pi^2*hbar^2/(2*m_e*L^2)) / 1.602e-19;
printf('\nInfinite well — comparison:\n');
printf('n  Numeric(eV)  Analytic(eV)  Error\n');
for n=1:5
  printf('%d  %10.4f  %12.4f  %8.2e\n', n, E1(n)/1.602e-19, E_analytic(n), ...
         abs(E1(n)/1.602e-19-E_analytic(n)));
end

%% ─── Example 2: Harmonic oscillator ──────────────────────────
omega = 1e14;  x_sc = sqrt(hbar/(m_e*omega));
x2 = linspace(-6*x_sc, 6*x_sc, 600);
V2 = 0.5*m_e*omega^2*x2.^2;
[E2, psi2] = solve_schrodinger(x2, V2, m_e, hbar);
printf('\nHO: ℏω = %.4f eV\n', hbar*omega/1.602e-19);
for n=0:4
  E_exact = hbar*omega*(n+0.5)/1.602e-19;
  printf('n=%d: num=%.6f eV, exact=%.6f eV, err=%.1e\n', ...
         n, E2(n+1)/1.602e-19, E_exact, abs(E2(n+1)/1.602e-19-E_exact));
end

%% ─── Example 3: Double-well potential ────────────────────────
x3 = linspace(-3e-9, 3e-9, 800);
V0 = 5*1.602e-19;   % 5 eV barrier
d  = 0.3e-9;
V3 = zeros(1,length(x3));
V3(abs(x3) < d) = V0;  % central barrier
[E3, psi3] = solve_schrodinger(x3, V3, m_e, hbar);
printf('\nDouble well: splitting E2-E1 = %.4f meV\n', (E3(2)-E3(1))*1e3/1.602e-19);

3 — Quantum Harmonic Oscillator: Ladder Operators & Coherent States

QHO — Hermite Polynomials & Coherent StatesGNU Octave
%% Quantum harmonic oscillator: analytic eigenstates and coherent states
%% Uses dimensionless units: ξ = x√(mω/ℏ),  ε_n = E_n/(ℏω) = n+1/2
clear; clc;

%% Hermite polynomials by recurrence H_{n+1} = 2ξH_n - 2nH_{n-1}
function H = hermite(n, xi)
  if n == 0,   H = ones(1,length(xi));
  elseif n==1, H = 2*xi;
  else
    H0 = ones(1,length(xi)); H1 = 2*xi;
    for k=2:n
      Hn = 2*xi.*H1 - 2*(k-1)*H0;
      H0=H1; H1=Hn;
    end; H=H1;
  end
end

function psi = qho_state(n, xi)
  Hn = hermite(n, xi);
  psi = (1/pi)^(1/4) / sqrt(2^n * factorial(n)) .* Hn .* exp(-xi.^2/2);
end

xi = linspace(-5, 5, 1000);
V_dim = 0.5*xi.^2;   % dimensionless potential ξ²/2

figure(1); clf; hold on;
cols = {'b','g','r','m','c','y'};
for n = 0:5
  psi_n = qho_state(n, xi);
  E_n = n + 0.5;
  plot(xi, 0.4*psi_n + E_n, cols{n+1}, 'LineWidth', 1.5);
  plot(get(gca,'XLim'), [E_n E_n], ':', 'Color', cols{n+1}, 'LineWidth', 0.8);
end
plot(xi, V_dim, 'k-', 'LineWidth', 2);  % potential
axis([-5 5 -0.5 7]); grid on;
title('QHO: ψ_n(ξ) + E_n (offset plot)');
xlabel('ξ = x√(mω/ℏ)'); ylabel('Energy (units of ℏω)');
legend({'n=0','','n=1','','n=2','','n=3','','n=4','','n=5','','V(ξ)'});

%% Coherent state: α = 2, |α⟩ = e^{-|α|²/2} Σ αⁿ/√n! |n⟩
alpha = 2;  N_max = 20;
t_arr = linspace(0, 4*pi, 200);  % 2 oscillation periods
figure(2);
for kt = 1:5:50
  t = t_arr(kt);
  PSI = zeros(1,length(xi));
  for n = 0:N_max
    cn = exp(-abs(alpha)^2/2) * alpha^n / sqrt(factorial(n));
    PSI += cn * exp(-1i*(n+0.5)*t) * qho_state(n, xi);
  end
  plot(xi, abs(PSI).^2); hold on;
end
title(sprintf('Coherent state α=%d: |⟨x|α,t⟩|² at different times',alpha));
xlabel('ξ'); ylabel('|Ψ|²'); grid on;

4 — Hydrogen Atom: Radial Wavefunctions & Probability

Hydrogen Atom Radial WavefunctionsGNU Octave
%% Hydrogen atom: radial wavefunctions R_nl(r) and associated quantities
%% Uses exact normalized forms with associated Laguerre polynomials
clear; clc;

a0 = 5.292e-11;  % Bohr radius (m)
E_ryd = 13.606;   % Rydberg energy (eV)

%% Associated Laguerre polynomial L_q^p(x)
function L = assoc_laguerre(q, p, x)
  % L_q^p via explicit sum: L_q^p(x) = Σ_{j=0}^{q} C(q+p,q-j)(-x)^j/j!
  L = zeros(1,length(x));
  for j = 0:q
    L += nchoosek(q+p,q-j) * (-x).^j / factorial(j);
  end
end

%% Normalised radial wavefunction R_nl(r)
function R = radial_wf(n, l, r, a0)
  rho = 2*r/(n*a0);
  Norm = -sqrt((2/(n*a0))^3 * factorial(n-l-1) / (2*n*factorial(n+l)^3));
  Lag  = assoc_laguerre(n-l-1, 2*l+1, rho);
  R = Norm .* exp(-rho/2) .* rho.^l .* Lag;
end

r = linspace(0, 25*a0, 2000);

states = {1,0; 2,0; 2,1; 3,0; 3,1; 3,2};
names   = {'1s', '2s', '2p', '3s', '3p', '3d'};
figure(3); clf;

for k = 1:6
  n = states{k,1}; l = states{k,2};
  R = radial_wf(n, l, r, a0);
  P = r.^2 .* R.^2;   % radial probability density P(r) = r²|R|²

  %% Most probable radius: where dP/dr = 0
  [~,idx] = max(P);
  r_mp = r(idx)/a0;
  r_mean = trapz(r, r .* P) / a0;   % ⟨r⟩
  r_rms  = sqrt(trapz(r, r.^2 .* P)) / a0; % √⟨r²⟩

  printf('%s: r_mp=%.2fa₀, ⟨r⟩=%.2fa₀, √⟨r²⟩=%.2fa₀, E=%.3feV\n', ...
         names{k}, r_mp, r_mean, r_rms, -E_ryd/n^2);

  subplot(2,3,k);
  plot(r/a0, P*a0, 'LineWidth', 1.8); hold on;
  plot([r_mp r_mp],[0 max(P*a0)],'r--');
  title(sprintf('%s: P(r) = r²|R|²',names{k}));
  xlabel('r/a₀'); ylabel('P(r)·a₀'); grid on;
end

5 — Spin-½ Matrix Mechanics & Stern-Gerlach

Spin-½ Algebra & Measurement SimulationGNU Octave
%% Spin-1/2 system: Pauli matrices, measurements, Bloch sphere, spin precession
clear; clc;
hbar = 1.0546e-34;

%% Pauli matrices
sx = [0 1; 1 0]/2;  sy = [0 -1i; 1i 0]/2;  sz = [1 0; 0 -1]/2;
Sx = hbar*sx; Sy = hbar*sy; Sz = hbar*sz;

%% Commutator verification: [Sx,Sy] = iℏSz
comm_xy = Sx*Sy - Sy*Sx;
printf('[Sx,Sy]/(iℏSz) = %.6f (should be 1)\n', ...
       real(comm_xy(1,2))/(1i*hbar*Sz(1,2)));

%% Eigenvalues and eigenvectors of S_x, S_y, S_z
for [S, name] in {Sx,'Sx'; Sy,'Sy'; Sz,'Sz'}'
  [V,D] = eig(S);
  printf('%s eigenvalues: ±%.4e J·s\n', name, abs(D(1,1)));
end

%% Spin precession in magnetic field B = Bz ẑ
% H = -γ B·S, γ_e = e/m_e (gyromagnetic ratio)
gamma_e = 1.761e11;  % rad/(T·s)
B0 = 1.0;             % Tesla
omega_L = gamma_e * B0; % Larmor frequency
printf('\nLarmor frequency: %.4e rad/s = %.2f GHz\n', omega_L, omega_L/2/pi/1e9);

% Initial state: spin-up along x: |+x⟩ = (|↑⟩ + |↓⟩)/√2
chi0 = [1; 1]/sqrt(2);
t_arr = linspace(0, 4*pi/omega_L, 1000);

Sx_exp = zeros(1,length(t_arr));
Sy_exp = Sx_exp; Sz_exp = Sx_exp;

for k = 1:length(t_arr)
  t = t_arr(k);
  % U(t) = exp(-iHt/ℏ) = exp(iω_L t σ_z/2)
  U = [exp(1i*omega_L*t/2), 0; 0, exp(-1i*omega_L*t/2)];
  chi_t = U * chi0;
  Sx_exp(k) = real(chi_t' * Sx * chi_t) / hbar;
  Sy_exp(k) = real(chi_t' * Sy * chi_t) / hbar;
  Sz_exp(k) = real(chi_t' * Sz * chi_t) / hbar;
end

figure(4);
plot(t_arr*1e9, Sx_exp, t_arr*1e9, Sy_exp, t_arr*1e9, Sz_exp, 'LineWidth',2);
legend({'⟨Sx⟩/ℏ','⟨Sy⟩/ℏ','⟨Sz⟩/ℏ'}); grid on;
title('Larmor precession of ⟨S⟩ in B = 1 T (z-axis)');
xlabel('t (ns)'); ylabel('⟨S_i⟩/ℏ');

6 — Perturbation Theory: Stark Effect

Linear Stark Effect — Hydrogen in Electric FieldGNU Octave
%% Linear Stark effect: hydrogen n=2 states in electric field E_field
%% H' = eE_field·z = eE_field·r cosθ
%% Demonstrates degenerate perturbation theory
clear; clc;
a0 = 5.292e-11; e = 1.602e-19; E_ryd = 13.606*e;

%% n=2 states: |2,0,0⟩ (2s), |2,1,-1⟩, |2,1,0⟩ (2p₀), |2,1,1⟩
%% Only non-zero matrix element: ⟨2,0,0|z|2,1,0⟩ = -3√3 a₀
%% (all others vanish by selection rules Δm=0, Δl=±1)

z_matrix_elem = -3*sqrt(3)*a0;  % ⟨2,0,0|r cosθ|2,1,0⟩
E_fields = linspace(0, 5e8, 200);  % V/m

%% In the {|2,0,0⟩, |2,1,0⟩} subspace (Δm=0), the perturbation matrix is:
%% H'_sub = e·E·[[0, z_01]; [z_10, 0]] — off-diagonal (linear Stark)
E_shifts = zeros(2, length(E_fields));

for k = 1:length(E_fields)
  Ef = E_fields(k);
  W = e * Ef * z_matrix_elem;   % coupling element
  H_sub = [0, W; W, 0];        % (E0 subtracted)
  ev = eig(H_sub);
  E_shifts(:,k) = sort(ev);
end

figure(5);
plot(E_fields/1e8, E_shifts(1,:)/e*1000, 'b-', ...
     E_fields/1e8, E_shifts(2,:)/e*1000, 'r-', 'LineWidth', 2);
legend({'E_{-} (|2s⟩-|2p⟩)/√2','E_{+} (|2s⟩+|2p⟩)/√2'});
title('Linear Stark Effect: n=2 Hydrogen');
xlabel('Electric field (10⁸ V/m)'); ylabel('Energy shift (meV)'); grid on;
printf('At E=10⁸ V/m: ΔE = ±%.3f meV\n', abs(E_shifts(1,100)/e*1000));

7 — WKB Tunneling: Transmission Coefficient

WKB Tunneling Through Arbitrary BarrierGNU Octave
%% WKB transmission coefficient T ≈ exp(-2/ℏ ∫√(2m(V-E))dx)
%% Compare with exact transfer-matrix result for rectangular barrier
clear; clc;
hbar = 1.0546e-34;  m_e = 9.109e-31;  eV = 1.602e-19;

%% ─── Rectangular barrier: exact vs WKB ────────────────────
V0 = 5*eV;   a = 1e-9;    % 5 eV, 1 nm
E_arr = linspace(0.01, 0.99, 200) * V0;

T_exact = zeros(1,length(E_arr));
T_WKB   = T_exact;

for k = 1:length(E_arr)
  E = E_arr(k);
  kappa = sqrt(2*m_e*(V0-E))/hbar;
  %% Exact: T = [1 + V0²sinh²(κa)/(4E(V0-E))]⁻¹
  T_exact(k) = 1/(1 + V0^2*sinh(kappa*a)^2/(4*E*(V0-E)));
  %% WKB: T ≈ exp(-2κa)
  T_WKB(k) = exp(-2*kappa*a);
end

figure(6); clf;
subplot(1,2,1);
semilogy(E_arr/eV, T_exact, 'b-', E_arr/eV, T_WKB, 'r--', 'LineWidth',2);
legend({'Exact','WKB'}); grid on;
title('Rectangular Barrier T(E): exact vs WKB');
xlabel('E (eV)'); ylabel('Transmission T');

%% ─── WKB for Gaussian barrier V(x) = V₀ exp(-x²/2σ²) ──────
sigma_b = 0.5e-9;  V0g = 3*eV;
x_b = linspace(-4*sigma_b, 4*sigma_b, 2000);
E_arr2 = (0.1:0.05:0.95)*V0g;
T_gauss = zeros(1,length(E_arr2));

for k = 1:length(E_arr2)
  E = E_arr2(k);
  Vx = V0g*exp(-x_b.^2/(2*sigma_b^2));
  classically_forbidden = Vx > E;
  integrand = sqrt(2*m_e*max(Vx-E,0));
  G = (2/hbar)*trapz(x_b(classically_forbidden), integrand(classically_forbidden));
  T_gauss(k) = exp(-G);
end

subplot(1,2,2);
semilogy(E_arr2/eV, T_gauss, 'g-', 'LineWidth',2); grid on;
title('WKB: Gaussian barrier V₀e^{-x²/2σ²}');
xlabel('E (eV)'); ylabel('T_{WKB}');

8 — Split-Operator FFT Method: TDSE Time Evolution

Split-Operator FFT — Real-Time TDSEGNU Octave
%% Split-operator Fourier method for real-time Schrödinger evolution
%% Ψ(t+dt) ≈ exp(-iV dt/2ℏ) · FFT⁻¹[exp(-ik²ℏdt/2m) · FFT[Ψ]] · exp(-iV dt/2ℏ)
%% Spectral accuracy O(dt²) with exact treatment of kinetic operator

clear; clc;
hbar = 1.0546e-34;  m_e = 9.109e-31;

%% Grid setup
L = 20e-9;  N = 1024;  % 20 nm box, 1024 points
x = linspace(-L/2, L/2, N);
dx = x(2)-x(1);

%% Momentum grid (FFT-ordered)
dk = 2*pi/L;
k = [0:N/2-1, -N/2:-1] * dk;

%% Potential: harmonic oscillator well
omega = 2e13;
V = 0.5*m_e*omega^2*x.^2;

%% Initial Gaussian wave packet (displaced from equilibrium)
x0 = 3e-9;  sigma = 0.5e-9;  k0 = 0;
psi = (2*pi*sigma^2)^(-0.25) * exp(-(x-x0).^2/(4*sigma^2) + 1i*k0*(x-x0));
psi /= sqrt(trapz(x, abs(psi).^2));  % normalise

%% Time evolution parameters
dt = 1e-16;  N_steps = 5000;  N_save = 50;
save_every = round(N_steps/N_save);
x_expect = zeros(1,N_save);  t_saved = x_expect;

%% Precompute phase factors
phase_V  = exp(-1i*V*dt/(2*hbar));   % V half-step
phase_T  = exp(-1i*hbar*k.^2*dt/(2*m_e));  % T full step

save_k = 0;
for step = 1:N_steps
  psi = phase_V .* psi;                % half-step in V
  psi = ifft(phase_T .* fft(psi));  % full-step in T (spectral)
  psi = phase_V .* psi;                % half-step in V

  if mod(step, save_every) == 0
    save_k++;
    rho = abs(psi).^2;
    x_expect(save_k) = trapz(x, x .* rho);
    t_saved(save_k)  = step*dt;
  end
end

%% ⟨x⟩(t) should oscillate at frequency ω (Ehrenfest)
T_osc = 2*pi/omega;
printf('Classical period: %.4f fs\n', T_osc/1e-15);

figure(7);
subplot(2,1,1);
plot(t_saved/1e-15, x_expect/1e-9, 'g-', 'LineWidth',2); grid on;
title('Ehrenfest: ⟨x⟩(t) for displaced Gaussian in QHO');
xlabel('t (fs)'); ylabel('⟨x⟩ (nm)');

%% Check norm conservation
norm_check = abs(trapz(x, abs(psi).^2) - 1);
printf('Norm conservation error: %.2e\n', norm_check);

XXVI. Interactive Simulations

Simulation I — Particle in an Infinite Square Well

Infinite Square Well — ψ_n(x) & |ψ|² & Energy Levels
1
n=1, E₁=0.376 eV (electron in 1 nm box)

Simulation II — Quantum Harmonic Oscillator

QHO — Eigenstates on Potential & Classical Turning Points
0
n=0, E=½ℏω (zero-point energy)

Simulation III — Gaussian Wave Packet Evolution

Free Particle Gaussian Wave Packet — Spreading & Phase
10
8
|ψ|² (green), Re[ψ] (cyan), Im[ψ] (pink)

Simulation IV — Double Slit: Quantum Interference Build-Up

Double Slit — Single-Particle Interference Pattern
12
3
N = 0 particles

Simulation V — Hydrogen Orbital Probability Density

Hydrogen ψ_nlm — |ψ(x,z)|² Cross-Section
Click Render to compute |ψ_nlm(x,z)|² heatmap