Linearity of the Schrödinger equation means $A\Psi_{\text{raw}}$ is a state for any constant $A$; Born's rule fixes $|A|$. Demand

So the normalized state is $\Psi=\Psi_{\text{raw}}/\sqrt{\int|\Psi_{\text{raw}}|^2dx}$ (the overall phase of $A$ is physically irrelevant and left free). By construction $\int|\Psi|^2dx=\int|\Psi_{\text{raw}}|^2dx\big/\int|\Psi_{\text{raw}}|^2dx=1$. This is exactly what normalize does — divide by $\sqrt{\int|\Psi_{\text{raw}}|^2dx}$ — so total_probability(normalize(raw, x), x) returns $1.0$.

The Born density is the normal density $|\Psi|^2=\frac{1}{\sqrt{2\pi}\,\sigma}e^{-(x-x_0)^2/2\sigma^2}$ with $x_0=2,\ \sigma=1$. The interval $[1,3]$ is exactly $x_0\pm\sigma$, so substituting $u=(x-x_0)/\sigma$,

using $\frac{1}{\sqrt{2\pi}}\int_{-a}^{a}e^{-u^2/2}du=\operatorname{erf}(a/\sqrt2)$. Over the whole line the integral is the normalization, $=1$. These match prob_between(psi, x, 1, 3) $\approx0.68269$ and prob_between(psi, x, x[0], x[-1]) $\approx1.0$.

For any mean $x_0$ and width $\sigma$, the substitution $u=(x-x_0)/\sigma$ removes both parameters:

a pure number depending only on $L$. Evaluating: $\operatorname{erf}(1/\sqrt2)=0.6827$, $\operatorname{erf}(2/\sqrt2)=0.9545$, $\operatorname{erf}(3/\sqrt2)=0.9973$ — the familiar 68–95–99.7% rule. The loop prints these alongside math.erf(L/math.sqrt(2)), agreeing to grid accuracy.

The phase $e^{ik_0x}$ cancels in $|\Psi|^2=\frac{1}{\sqrt{2\pi}\,\sigma}e^{-(x-x_0)^2/2\sigma^2}$, a normal density of mean $x_0$ and variance $\sigma^2$. Its first moment is the mean,

By the variance theorem $\sigma_x^2=\langle x^2\rangle-\langle x\rangle^2=\sigma^2$, so $\langle x^2\rangle=x_0^2+\sigma^2=4+1=5$ and $\sigma_x=\sqrt{5-4}=\sigma=1$. These reproduce expectation_x(psi, x) $\approx2$, expectation_x2(psi, x) $\approx5$, sigma_x(psi, x) $\approx1$.

Write $\Psi=g(x)\,e^{ik_0x}$ with $g$ the real Gaussian envelope. Then $\hat p\Psi=-i\hbar\,\partial_x\Psi=(\hbar k_0\,g-i\hbar\,g')\,e^{ik_0x}$, so

Integrating: $\int g^2dx=1$ gives $\hbar k_0$, while $\int g g'dx=\tfrac12\int(g^2)'dx=0$ (boundary terms vanish), so $\langle p\rangle=\hbar k_0=5$. For the spread, $\langle p^2\rangle=\hbar^2\!\int|\partial_x\Psi|^2dx=\hbar^2\!\big(k_0^2+\tfrac{1}{4\sigma^2}\big)$, hence $\sigma_p^2=\hbar^2/4\sigma^2$ and $\sigma_p=\hbar/2\sigma=0.5$. This matches expectation_p_fft(psi, x) $=5.0$ (the finite-difference expectation_p $\approx4.998$) and sigma_p(psi, x) $\approx0.499$.

Combine the two closed forms from P4 and P5, $\sigma_x=\sigma$ and $\sigma_p=\hbar/2\sigma$:

The $\sigma$ cancels, so every Gaussian — no matter how wide or narrow — has the same product $\hbar/2$, the exact floor of the Heisenberg inequality $\sigma_x\sigma_p\ge\hbar/2$ (~QM-07). Squeezing the packet in $x$ (small $\sigma$) inflates $\sigma_p=\hbar/2\sigma$ in lockstep. Numerically sigma_x(psi, x)*sigma_p(psi, x) $\approx0.4994\approx\hbar/2=0.5$ (the tiny shortfall is finite-difference/grid error).

Insert $\Psi=(2\pi\sigma^2)^{-1/4}e^{-(x-x_0)^2/4\sigma^2}e^{ik_0x}$ into $\Phi(p)=(2\pi\hbar)^{-1/2}\!\int\Psi\,e^{-ipx/\hbar}dx$. Collecting exponentials, the integral is the Fourier transform of a Gaussian at wavenumber $q=p/\hbar-k_0$, which (completing the square) gives $\Phi(p)\propto\exp\!\big[-\sigma^2(p/\hbar-k_0)^2\big]=\exp\!\big[-\sigma^2(p-\hbar k_0)^2/\hbar^2\big]$. Therefore

a Gaussian centred at $p_0=\hbar k_0=5$ with width $\hbar/2\sigma$. Plancherel's theorem makes $\int|\Phi|^2dp=\int|\Psi|^2dx=1$. The code confirms the peak at p[argmax|phi|^2] $\approx5.0$ and np.trapz(|phi|^2, p) $\approx1.0$.

Free evolution acts in momentum space as a pure phase: $\widetilde\Psi(k,t)=e^{-i\hbar k^2 t/2m}\,\widetilde\Psi(k,0)$. Since $|e^{-i\hbar k^2t/2m}|=1$, every mode keeps its modulus, so by Plancherel

for all $t$ — normalization is conserved. The phase does, however, dephase the modes and broaden the packet: a Gaussian spreads as $\sigma_x(t)=\sigma_0\sqrt{1+(\hbar t/2m\sigma_0^2)^2}$. At $t=2,\ m=\hbar=\sigma_0=1$ this is $\sqrt{1+1}=\sqrt2\approx1.4142$, matching total_probability(pt, xx) $\approx1.0$ and sigma_x(pt, xx) $\approx1.41421$.

Ehrenfest's theorem gives $\frac{d}{dt}\langle x\rangle=\frac{\langle p\rangle}{m}$ and $\frac{d}{dt}\langle p\rangle=\langle-\partial_xV\rangle$. With $V=0$ there is no force, so $\langle p\rangle=\hbar k_0$ stays constant; integrating the first equation,

i.e. the centroid drifts at the group velocity $\hbar k_0/m$. For $x_0=0,\ k_0=4,\ m=\hbar=1$, at $t=1$ this is $\langle x\rangle=0+4\cdot1=4$, while $\langle p\rangle=\hbar k_0=4$ is unchanged. These match expectation_x(free_propagate(p0, xx, 1.0), xx) $\approx4.0$ and expectation_p_fft(...) $\approx4.0$.

Insert $\Psi=\psi(x)\varphi(t)$ into $i\hbar\,\partial_t\Psi=\hat H\Psi$. The left side is $i\hbar\,\psi\dot\varphi$ and the right is $\varphi\big(-\tfrac{\hbar^2}{2m}\psi''+V\psi\big)$; divide through by $\psi\varphi$:

The left depends only on $t$, the right only on $x$, so for the equality to hold for all $x,t$ both must equal one constant $E$. The time half integrates to $\dot\varphi=-\tfrac{iE}{\hbar}\varphi\Rightarrow\varphi(t)=e^{-iEt/\hbar}$, and the space half is the TISE $\hat H\psi=E\psi$. Every separable state shares this phase: with $\hbar=1$, $t=\pi/E$ gives $e^{-i\pi}=-1$ (returns $-\psi$) and $t=2\pi/E$ gives $e^{-2\pi i}=1$ (returns $+\psi$) — exactly what stationary_state(psi[:,0], E[0], ...) produces.

Inside the well $V=0$, so the TISE is $-\tfrac{\hbar^2}{2m}\psi''=E\psi$, i.e. $\psi''=-k^2\psi$ with $k=\sqrt{2mE}/\hbar$. The general solution $A\sin kx+B\cos kx$ must vanish at $x=0$ (so $B=0$) and at $x=L$, forcing $\sin kL=0$, i.e. $k_nL=n\pi$. Hence

For $n=1\ldots5$ this is $4.9348,\,19.7392,\,44.4132,\,78.9568,\,123.3701$. The 400-point finite-difference grid returns $4.9348,\,19.7388,\,44.4112,\,78.9504,\,123.354$ — agreeing to $<0.02\%$, and E[0] $\approx$ infinite_well_energy(1,1) $\approx 4.9348$.

Normalizing $\psi_n=A\sin(n\pi x/L)$: $\int_0^L A^2\sin^2(n\pi x/L)\,dx=A^2\tfrac L2=1$ gives $A=\sqrt{2/L}$, so $\psi_n=\sqrt{2/L}\sin(n\pi x/L)$. Orthogonality uses the product-to-sum identity:

since each cosine integrates to zero unless $m=n$. Thus $\langle\psi_m|\psi_n\rangle=\delta_{mn}$. The ground density $|\psi_1|^2=(2/L)\sin^2(\pi x/L)$ peaks at $x=L/2$ where $\sin^2=1$, giving $2/L=2$. The grid confirms max prob_density(psi[:,0]) = 2.0 and inner_product(psi[:,0], psi[:,1], dx) $\approx-1.4\times10^{-16}$.

A single separable state is $\Psi_n(x,t)=\psi_n(x)\,e^{-iE_nt/\hbar}$. Its Born density factorizes the modulus:

because $|e^{-iE_nt/\hbar}|^2=e^{-iE_nt/\hbar}\,e^{+iE_nt/\hbar}=1$ ($E_n$ is real). So the density — and hence every probability and every expectation of a time-independent operator — is frozen, even though $\Psi$ itself keeps rotating in phase at rate $E_n/\hbar$. That is the sense of "stationary." The code confirms prob_density(stationary_state(psi[:,0],E[0],t)) is $t$-independent to $\sim10^{-16}$.

Each term evolves with its own phase, so $\Psi(x,t)=\tfrac1{\sqrt2}\big(\psi_1e^{-iE_1t/\hbar}+\psi_2e^{-iE_2t/\hbar}\big)$. With real $\psi_1,\psi_2$ the density is

the cross-term reducing to a cosine because only the relative phase $e^{-i(E_2-E_1)t/\hbar}$ survives. The beat angular frequency is $\omega=(E_2-E_1)/\hbar=(4-1)\pi^2/2=3\pi^2/2\approx14.80$, so $\langle x\rangle$ oscillates about $L/2$ with period $T=2\pi/\omega\approx0.4244$. This matches superposition_period(E[1],E[0]) $\approx0.4244$ and measure_period on the sampled $\langle x\rangle(t)$.

Fourier's trick projects onto the orthonormal eigenbasis: multiply $\Psi_0=\sum_m c_m\psi_m$ by $\psi_n^\ast$ and integrate, using $\langle\psi_n|\psi_m\rangle=\delta_{nm}$:

For $\Psi_0=0.6\,\psi_1+0.8\,\psi_5$ only the matching terms survive, so $c_1=0.6$, $c_5=0.8$ (all others $0$); note $0.6^2+0.8^2=1$, so $\Psi_0$ is already normalized. The Born rule for energy gives $|c_1|^2=0.36$ and $|c_5|^2=0.64$, summing to $1$, and $\langle H\rangle=0.36\,E_1+0.64\,E_5$ — constant, since the $c_n$ never change. The code returns exactly c1=0.6, c5=0.8, sum|c_n|^2 = 1.0, and psi @ c rebuilds $\Psi_0$.

The Hamiltonian assembly $H=-\tfrac{\hbar^2}{2m}\,\mathrm{tridiag}(1,-2,1)/(\Delta x)^2+\mathrm{diag}\,V_i$ is agnostic to $V$, so swapping the flat $V=0$ for the parabola $V=\tfrac12 m\omega^2x^2$ and re-diagonalizing solves the oscillator TISE $-\tfrac{\hbar^2}{2m}\psi''+\tfrac12 m\omega^2x^2\psi=E\psi$. Its exact spectrum (derived analytically with ladder operators in ~QM-09) is

With $\hbar=m=\omega=1$ this is $0.5,1.5,2.5,3.5,4.5$ — equally spaced by $\hbar\omega=1$. The grid run on $[-8,8]$ returns $0.5,1.4999,2.4998,3.4997,4.4995$ ($\sim0.1\%$), so Eh[0] $\approx0.5$ and Eh[4] $\approx4.5$.

Solve each equation for its time derivative,

and assemble $\partial_t\rho=\Psi^\partial_t\Psi+\Psi\,\partial_t\Psi^$:

The potential piece $-\tfrac{i}{\hbar}(V-V^)|\Psi|^2$ vanishes iff $V=V^$ — that is the one and only place reality of $V$ is used (a complex $V$ would leave a source). With the identity $\Psi^\nabla^2\Psi-\Psi\nabla^2\Psi^=\nabla\!\cdot(\Psi^\nabla\Psi-\Psi\nabla\Psi^)$,

Numerically the residual $\partial_t\rho+\partial_x j$ on a free-evolved packet is $1.6\times10^{-4}$, far below the $\sim0.12$ magnitude of $\partial_t\rho$ itself — the continuity equation holds to discretization error.

For any complex number $z$, $\mathrm{Im}\,z=\dfrac{z-z^}{2i}$. Take $z=\Psi^\nabla\Psi$; since $(\nabla\Psi)^=\nabla\Psi^$ its conjugate is $z^=\Psi\nabla\Psi^$, so

using $1/i=-i$ — exactly the second form. Reality is automatic: $\tfrac{\hbar}{m}\mathrm{Im}(\cdots)$ is $\hbar/m$ times a real number; equivalently the bracket $\Psi\nabla\Psi^-\Psi^\nabla\Psi=z^*-z=-2i\,\mathrm{Im}\,z$ is purely imaginary, and $i\times(\text{imaginary})$ is real. This is why the two arrays agree to $10^{-12}$ while the imaginary part of the second form is $\sim10^{-12}$ — one real current wearing two algebraic faces.

For $\Psi=A\,e^{i(kx-\omega t)}$ the density is $\rho=|\Psi|^2=|A|^2$ — uniform and time-independent. The derivative is $\partial_x\Psi=ik\Psi$, so

and therefore

So $v$ is the de Broglie velocity: the probability fluid drifts rigidly at $v$, the same $\rho\mathbf v$ flux carried by a mass current (~CM-22). Sending $k\to-k$ reverses $v$, hence $j$. With $k=1.5$, prob_current(...).mean()/prob_density(...).mean() $\approx1.4991\approx k=v$ (the $\sim1\%$ deficit is central-difference error), and it changes sign for $k\to-k$ — exactly $j=\rho v$.

Write $\Psi(x)=R(x)\,e^{i\alpha}$ with $R$ real and $\alpha$ a constant global phase. Then

which is real, so $\mathrm{Im}(\Psi^\nabla\Psi)=0$ and $\mathbf j=0$ everywhere. For the well mode $\psi_n\propto\sin(n\pi x/a)$ we have $R=\sin(n\pi x/a)$ real ($\alpha=0$), hence $\mathbf j\equiv0$. This does not contradict "motion": a real eigenstate is a standing wave, the equal superposition $\sin(kx)\propto e^{ikx}-e^{-ikx}$ of right- and left-movers whose currents are equal and opposite and cancel. There is real momentum content — $\langle p^2\rangle\neq0$, i.e. nonzero kinetic energy — but $\langle p\rangle=0$ and zero net flux, just as a standing wave on a string stores energy yet transports none. Thus prob_current(np.cos(1.5x), dx) is $0$ to $10^{-12}$ (in fact identically zero, since $\mathrm{Im}$ of a real array vanishes).

Integrate continuity over all space and pull the derivative out of the integral:

By the divergence theorem the right-hand side is a flux through the surface at infinity:

because any normalizable $\Psi$ — and with it $\mathbf j$ — decays to zero at infinity, killing the surface integral. Hence $\tfrac{d}{dt}\int|\Psi|^2\,dV=0$: normalized once, normalized forever, which is exactly what makes Schrödinger evolution consistent with the Born rule (~QM-02). Numerically total_probability holds at $1.000000$ across $50$ free-evolution steps (drift $<10^{-9}$) — global conservation confirmed.

Stationary states evolve by a pure phase, so (taking $c_1,c_2$ real WLOG — a phase in $c_n$ only shifts the sine)

Form $\Psi^*\partial_x\Psi$. The diagonal pieces $c_1^2\psi_1\psi_1'+c_2^2\psi_2\psi_2'$ are real (the $\psi_n$ are real), so they drop out of $\mathrm{Im}$ — each state alone gives $j=0$, recovering P4. Only the cross terms remain:

Taking the imaginary part with $\mathrm{Im}\,e^{\pm i\omega_{21}t}=\pm\sin\omega_{21}t$,

The current sloshes sinusoidally at the Bohr frequency $\omega_{21}=(E_2-E_1)/\hbar$, though each $\psi_n$ alone carries none; its spatial moment $\frac{d}{dt}\langle x\rangle\propto\int j\,dx$ is the oscillating dipole that radiates at the transition frequency (~QM-16). Building two real well-states with relative phase $e^{-i(E_2-E_1)t/\hbar}$ and evaluating prob_current shows exactly this back-and-forth slosh, of period $2\pi\hbar/(E_2-E_1)$.

Every one of these laws has the same skeleton $\partial_t\rho+\nabla\!\cdot\mathbf j=0$; only the conserved "stuff" and its flux change:

Written out, the phase-space row is the collisionless Boltzmann/Liouville (Vlasov) equation $\partial_t f+\mathbf v\cdot\nabla_x f+\mathbf a\cdot\nabla_v f=0$. What EM-14 adds is a source: bookkeeping field energy (Poynting) instead of charge turns bare continuity into a balance law

where $-\mathbf J\cdot\mathbf E$ is the rate the field does work on charges (Joule heating) — energy is exchanged between field and matter, so the field-energy flux is not source-free. Probability has no such source because $V$ is real: in P1 the potential terms cancelled exactly, leaving $\partial_t\rho+\nabla\!\cdot\mathbf j=0$ with zero right-hand side. (A complex optical potential $V=V_R-iV_I$ breaks this, adding a sink $-\tfrac{2V_I}{\hbar}|\Psi|^2$ that models absorption — the direct analog of $-\mathbf J\cdot\mathbf E$.) This is the B2 capstone: the QM-04 code verifies the probability row numerically — continuity_residual $\sim10^{-4}$ and total_probability fixed at $1.000000$ — and the residual stays at discretization noise precisely because the law is source-free.

Start from the definition $\langle\phi|\psi\rangle=\sum_i\phi_i^{*}\psi_i$ and swap the arguments:

which is conjugate symmetry. Setting $\phi=\psi$ gives $\langle\psi|\psi\rangle=\sum_i\phi_i^{}\psi_i\big|_{\phi=\psi}=\sum_i|\psi_i|^2$, a sum of non-negative reals, so it is real and $\ge0$. Numerically the bra is conjugated: inner_product([1,1j],[1,2]) $=1^{}\cdot1+(i)^{*}\cdot2=1-2i$, while reversing the slots gives $1+2i$ — the conjugate, exactly as the check reports.

Reality of $\langle Q\rangle=\langle\psi|\hat Q\psi\rangle$ means $\langle Q\rangle=\langle Q\rangle^{*}$. Conjugate symmetry turns the right side into a statement about the adjoint:

using $\langle\hat Q\psi|\psi\rangle=\langle\psi|\hat Q^\dagger\psi\rangle$ from the definition of $\hat Q^\dagger$. So $\langle\psi|(\hat Q-\hat Q^\dagger)\psi\rangle=0$ for every $\psi$; polarization (apply it to $\psi+\phi$ and $\psi+i\phi$) extends this to all matrix elements $\langle\phi|(\hat Q-\hat Q^\dagger)\psi\rangle=0$, forcing $\hat Q=\hat Q^\dagger$. Hence expectation(A, psi).imag $\approx0$ exactly when A is Hermitian, and is nonzero for a generic matrix.

Let $\hat A|n\rangle=\lambda_n|n\rangle$ with $\hat A=\hat A^\dagger$. (1) Real eigenvalues: $\lambda_n\langle n|n\rangle=\langle n|\hat A n\rangle=\langle\hat A n|n\rangle=\lambda_n^{}\langle n|n\rangle$, and $\langle n|n\rangle\neq0$, so $\lambda_n=\lambda_n^{}\in\mathbb R$. (2) Orthogonality: for $\lambda_m\neq\lambda_n$,

so $(\lambda_n-\lambda_m)\langle m|n\rangle=0\Rightarrow\langle m|n\rangle=0$. The eigenvectors thus form a complete orthonormal basis, and inserting $\hat I=\sum_n|n\rangle\langle n|$ gives the spectral form $\hat A=\hat A\hat I=\sum_n\lambda_n|n\rangle\langle n|$. The check confirms w.imag$\approx0$, dagger(V)@V$=I$, and reconstruct_from_spectrum(w,V) rebuilds the matrix.

Antisymmetry is immediate: $[A,B]=AB-BA=-(BA-AB)=-[B,A]$. For the Jacobi identity, expand one term, $[A,[B,C]]=A(BC-CB)-(BC-CB)A=ABC-ACB-BCA+CBA$, and likewise for the two cyclic rotations $A\to B\to C\to A$:

Adding all three, each of the six orderings $ABC,ACB,\dots$ appears once with $+$ and once with $-$, so the sum is $0$. Operators under $[\,\cdot\,,\cdot\,]$ therefore satisfy the Lie-algebra axioms, matching the check that commutator(A,B)$=-$commutator(B,A) and the Jacobi combination $\approx0$.

With $A=A^\dagger$ and $B=B^\dagger$, the adjoint reverses the product (and conjugates):

so $[A,B]$ is anti-Hermitian. Multiplying by $i$ flips the sign that the dagger pulls through:

so $i[A,B]$ is Hermitian — a genuine observable. This is the factor of $i$ in the generalized uncertainty bound. Accordingly is_hermitian(1j*commutator(A,B)) returns True for the random Hermitian pair.

Suppose $[A,B]=0$ and $A|n\rangle=\lambda_n|n\rangle$ with $\lambda_n$ nondegenerate. Then

so $B|n\rangle$ is an eigenvector of $A$ with the same eigenvalue $\lambda_n$; nondegeneracy forces $B|n\rangle\propto|n\rangle$, i.e. $B|n\rangle=\mu_n|n\rangle$. Thus the eigenbasis of $A$ diagonalizes $B$ simultaneously. The check builds $A=U\,\mathrm{diag}(1..5)\,U^\dagger$ and $B=U\,\mathrm{diag}(-2,.5,7,1,3)\,U^\dagger$ from a common $U$, so they commute (np.linalg.norm(commutator(A,B))$\approx1.3\times10^{-14}$) and $A$'s eigenvectors $V_A$ render $V_A^\dagger B\,V_A$ diagonal. Two generic Hermitian matrices have $[A,B]\neq0$ and no shared eigenbasis.

The scalar prefactors multiply to

and the operator part expands using $[a,a]=[a^\dagger,a^\dagger]=0$ and $[a^\dagger,a]=-[a,a^\dagger]$:

Therefore $[\hat x,\hat p]=\tfrac{i\hbar}{2}\cdot2[a,a^\dagger]=i\hbar\,[a,a^\dagger]$, which is $i\hbar$ wherever $[a,a^\dagger]=1$. With $\hbar=1$, canonical_commutator(8)[:7,:7]$\approx i\,I_7$ on the interior, as the check confirms.

Any commutator of finite matrices is traceless, since $\mathrm{tr}(\hat x\hat p)=\mathrm{tr}(\hat p\hat x)$, so $\mathrm{tr}[\hat x,\hat p]=0$; but $\mathrm{tr}(i\hbar I)=i\hbar N\neq0$ — the two cannot be equal, so $[\hat x,\hat p]=i\hbar I$ has no finite-dimensional representation. The deviation is localized: $a^\dagger a=\mathrm{diag}(0,\dots,N-1)$ is exact, but $aa^\dagger$ loses the rung $a^\dagger|N-1\rangle$ (it would land on the discarded $|N\rangle$), giving

So it is $i\hbar$ on the interior but the corner reads $-i\hbar(N-1)$. For $N=8$ the code gives np.trace(C)$\approx0$ and C[7,7]$=-7i=-i\hbar(N-1)$; zeroing that one entry makes C - 1j*np.eye(8)$\approx0$.

Let $\{|n\rangle\}$ be orthonormal, $\langle n|m\rangle=\delta_{nm}$, and complete (it spans $\mathcal H$). Any state expands as $|\psi\rangle=\sum_m c_m|m\rangle$ with $c_m=\langle m|\psi\rangle$. Apply the candidate operator:

Since this reproduces every $|\psi\rangle$, the operator is the identity: $\sum_n|n\rangle\langle n|=\hat I$. The content is completeness — the projectors leave nothing out only if the basis spans the whole space. Numerically, summing the projectors onto the 6 columns of a unitary, sum(projector(V[:,n]) for n in range(6)), returns $\approx I_6$.

Let $\hat A=\hat A^\dagger$ with $\hat A|a_n\rangle=a_n|a_n\rangle$ and $\langle a_n|a_n\rangle=1$. Sandwiching the operator and moving it onto the bra (Hermiticity) gives

so $a_n\in\mathbb R$ — real meter readings. For two eigenvectors, $\langle a_m|\hat A|a_n\rangle$ read rightward is $a_n\langle a_m|a_n\rangle$ and read leftward (Hermiticity, $a_m$ real) is $a_m\langle a_m|a_n\rangle$, so $(a_n-a_m)\langle a_m|a_n\rangle=0$ forces $\langle a_m|a_n\rangle=0$ when $a_n\neq a_m$ (degenerate blocks are orthogonalized by Gram–Schmidt). The spectral theorem then makes the normalized set complete, $\sum_n|a_n\rangle\langle a_n|=\mathbb 1$. This is exactly why eigensystem(A) returns real eigenvalues and V.conj().T @ V = V @ V.conj().T = I.

A measurement can only return one of the eigenvalues $a_n$ of $\hat A$. Projecting the expansion $|\psi\rangle=\sum_n c_n|a_n\rangle$ onto $\langle a_m|$ isolates the amplitude $c_m=\langle a_m|\psi\rangle$, and the generalized Born rule assigns

These form a probability distribution because completeness $\sum_n|a_n\rangle\langle a_n|=\mathbb 1$ and normalization give

Hence outcome_probabilities(psi, A)[1].sum() $=1$ for any state and any Hermitian A.

Insert $|\psi\rangle=\sum_n c_n|a_n\rangle$ into the sandwich and use $\hat A|a_n\rangle=a_n|a_n\rangle$ with $\langle a_m|a_n\rangle=\delta_{mn}$:

So the "average of outcomes" and the "sandwich" are the same number. It is real because every $a_n\in\mathbb R$ and every $|c_n|^2\ge0$ — equivalently $\langle\psi|\hat A|\psi\rangle^*=\langle\psi|\hat A^\dagger|\psi\rangle=\langle\psi|\hat A|\psi\rangle$. expectation(psi, A) evaluates both the spectral and sandwich forms and asserts they agree, returning that real value.

Write the spread as $\sigma_A^2=\langle(\hat A-\langle A\rangle)^2\rangle=\lVert(\hat A-\langle A\rangle)|\psi\rangle\rVert^2\ge0$, which vanishes iff $(\hat A-\langle A\rangle)|\psi\rangle=0$, i.e. $|\psi\rangle$ is an eigenstate ($\hat A\psi=q\psi$ with $q=\langle A\rangle$). For $|0\rangle$: $\sigma_z|0\rangle=+|0\rangle$ is an eigenstate, so $\mathrm{Var}=0$. For $|{+x}\rangle=\tfrac1{\sqrt2}(1,1)$ the diagonal $\sigma_z=\mathrm{diag}(1,-1)$ gives $\langle\sigma_z\rangle=\tfrac12(1-1)=0$, while $\sigma_z^2=\mathbb 1\Rightarrow\langle\sigma_z^2\rangle=1$, so

This matches variance(spin_state(0), sigma_z) $\approx0$ and variance(spin_state(np.pi/2,0), sigma_z) $\approx1$.

Collapse sends $|\psi\rangle\to|\psi'\rangle=\hat P_n|\psi\rangle/\lVert\hat P_n|\psi\rangle\rVert$. The weight of the same outcome on re-measurement is $P'(a_n)=\lVert\hat P_n|\psi'\rangle\rVert^2$. Using idempotence $\hat P_n^2=\hat P_n$,

so $\hat P_n$ leaves $|\psi'\rangle$ untouched and $P'(a_n)=\lVert\psi'\rVert^2=1$. The collapsed state is already an $a_n$-eigenstate, so the second projection changes nothing — this is the reproducibility checked by s = collapse(psi, A, n) then outcome_probabilities(s, A)[1][n] $=1$.

The prepared state is $|{+n}\rangle=(\cos\tfrac\theta2,\,e^{i\phi}\sin\tfrac\theta2)$. (a) With $\sigma_z=\mathrm{diag}(1,-1)$,

the spin-½ Malus law (so $\langle S_z\rangle=\tfrac\hbar2\cos\theta$); at $\theta=60^\circ$ this is $\cos60^\circ=\tfrac12$. (b) For $\theta=90^\circ$, $|{+x}\rangle=\tfrac1{\sqrt2}(1,1)$, and the $\sigma_z$-eigenstates are $|0\rangle,|1\rangle$, so

Confirmed by expectation(spin_state(np.pi/3), sigma_z) $=0.5$ and outcome_probabilities(spin_state(np.pi/2,0), sigma_z) $\to([-1,1],[0.5,0.5])$.

(a) If $[\hat B,\sigma_x]=0$ the two share an eigenbasis, so measuring $\sigma_x=+1$ leaves the state a simultaneous eigenstate and a later $\sigma_x$ still reads $+1$ — no disturbance. (b) Write $|{+x}\rangle=\tfrac1{\sqrt2}(|0\rangle+|1\rangle)$. The intervening $\sigma_z$ measurement collapses it to $|0\rangle$ (or $|1\rangle$). Re-expanding that result in the $\sigma_x$-basis, $|0\rangle=\tfrac1{\sqrt2}(|{+x}\rangle+|{-x}\rangle)$, gives

Because $\sigma_z$ fails to commute, $[\sigma_x,\sigma_z]=-2i\sigma_y\neq0$, it shares no common eigenbasis with $\sigma_x$ and randomizes the once-certain outcome. So outcome_probabilities(spin_state(np.pi/2,0), sigma_x)[1] $=[0,1]$, but after collapse(..., sigma_z, index_of(sigma_z,+1)) the $\sigma_x$ probabilities become $[0.5,0.5]$.

With $\hbar=1$, the density $|\psi|^2\propto e^{-(x-x_0)^2/2\sigma^2}$ is a Gaussian of variance $\sigma^2$, so $\langle x\rangle=x_0$ and $\sigma_x=\sigma$. The boost $e^{ip_0x}$ moves $\langle p\rangle$ to $p_0$ without changing the spread; the Fourier transform $\tilde\psi(p)$ is again Gaussian, of width $\hbar/2\sigma=1/2\sigma$, so $\sigma_p=1/2\sigma$. Hence

for every $\sigma,p_0$. Equality is forced because $\hat p\psi=\bigl[i(x-x_0)/2\sigma^2+p_0\bigr]\psi$, i.e. $(\hat p-\langle p\rangle)\psi=ia(\hat x-\langle x\rangle)\psi$ with $a=1/2\sigma^2$ real — exactly the minimum-uncertainty condition (Schwarz tight, $\operatorname{Re}z$ dropped). uncertainty_product(gaussian_packet(x, sigma=1.3, p0=0.7), x, dx) returns $0.49999\ldots\approx0.5$.

Multiplying the two equal spreads,

with equality only at $n=0$. The ground state $\psi_0\propto H_0(x)e^{-x^2/2}=e^{-x^2/2}$ is exactly the minimum-uncertainty Gaussian of P1 (width $\sigma=1/\sqrt2$, so $\sigma_x=\sigma_p=1/\sqrt2$ and the product is $\tfrac12$). Each excitation adds one full quantum $\hbar=1$ to the product, so every $n\ge1$ state is strictly fuzzier. uncertainty_product(ho_eigenstate(x, n), x, dx) returns $0.5,\,1.5,\,2.5,\,3.5$ for $n=0,1,2,3$ — i.e. $n+\tfrac12$.

Being non-Gaussian, the cat must lie strictly above the bound (P1 uniqueness). Concretely, with unit-width bumps at $x=\pm3$ and $\langle x\rangle=0$, the density has $\langle x^2\rangle\approx\sigma^2+3^2=10$, so $\sigma_x\approx\sqrt{10}\approx3.15$; the coherent superposition prints interference fringes in $\tilde\psi(p)$ that hold $\sigma_p\approx0.47$. Their product

is comfortably consistent with — but far from saturating — the inequality. uncertainty_product(cat, x, dx) returns $1.4935$, confirming $>0.5$.

Feed $[\hat S_x,\hat S_y]=i\hbar\hat S_z$ into the Robertson bound:

For $|{\uparrow_z}\rangle=(1,0)$, a $\hat S_z$ eigenstate with $\langle\hat S_z\rangle=\hbar/2$, the RHS is $\hbar^2/4$. Since $\hat S_x=\tfrac\hbar2\sigma_x$ with $\sigma_x^2=\mathbb 1$, we get $\langle S_x\rangle=0$, $\langle S_x^2\rangle=\hbar^2/4$, so $\sigma_{S_x}=\hbar/2$; identically $\sigma_{S_y}=\hbar/2$, making the LHS $\hbar^2/4$ as well — saturation. In $\hbar=1$ units both sides equal $0.25$, exactly what std(Sx,[1,0])*std(Sy,[1,0]) and generalized_bound(Sx,Sy,[1,0]) return.

The Robertson right-hand side is $\tfrac12\lvert\langle[\hat A,\hat B]\rangle\rvert$. If $[\hat A,\hat B]=0$ it is identically $0$, so $\sigma_A\sigma_B\ge0$ is vacuous — nothing forbids both spreads vanishing at once, and indeed commuting Hermitian operators are simultaneously diagonalizable (a shared complete eigenbasis). Two diagonal matrices commute entrywise:

So commutator(A,B) is the zero matrix and generalized_bound(A,B,[1,1,1]) $=0.0$, while each $\sigma_A,\sigma_B\ge0$ separately. Incompatible (noncommuting) observables, by contrast, carry a nonzero RHS and cannot both be sharp.

Ehrenfest's theorem reads $\dfrac{d\langle x\rangle}{dt}=\dfrac{\langle p\rangle}{m}$ and $\dfrac{d\langle p\rangle}{dt}=-\langle V'(\hat x)\rangle$. With $V=0$ the force $V'=0$, so $\langle p\rangle$ stays fixed at its initial value $p_0$ and

a straight drift. The split-step run confirms it: at the midpoint gradient(xs,t) $=1.500$ equals ps/M $=1.500$, and ps is constant to $\sim5\times10^{-15}$ across all 120 steps. So $d\langle x\rangle/dt=\langle p\rangle/m$ with $d\langle p\rangle/dt=0$, as required.

Here $V=\tfrac12 m\omega^2x^2$ has the linear derivative $V'(x)=m\omega^2x$, so the mean force is exact: $\langle V'(\hat x)\rangle=m\omega^2\langle x\rangle=V'(\langle x\rangle)$. Ehrenfest then closes into a self-contained pair for $(\langle x\rangle,\langle p\rangle)$,

whose solution is the classical SHO $\langle x\rangle(t)=x_0\cos\omega t+(p_0/m\omega)\sin\omega t$. Over one period the code gives max|xs - classical_sho| $\approx6.0\times10^{-5}\sim10^{-5}$, and at $t=T/4$ gradient(ps,t) $\approx$ Fs (both $\approx-0.0024$) — the packet centre rides the classical trajectory exactly.

(Griffiths's Problem A.5, hint and all — the lemma behind P4 and the whole of §2 of the notes.) If $\langle f|f\rangle=0$ then $f=0$ and both sides vanish, so let $\langle f|f\rangle>0$ and define

$g$ minus its component along $f$. First, $h\perp f$:

Next expand $\langle h|h\rangle$ with $c=\langle f|g\rangle/\langle f|f\rangle$ (antilinear first slot, linear second):

By conjugate symmetry $\langle g|f\rangle=\langle f|g\rangle^*$, each of the last three terms equals $|\langle f|g\rangle|^2/\langle f|f\rangle$ (two cancel), so

Equality forces $\langle h|h\rangle=0$, hence $h=0$, hence $g=cf$ — one state a multiple of the other. (In §3 of the notes, demanding this and $\operatorname{Re}\langle f|g\rangle=0$ makes $c$ pure imaginary and turns the condition into the minimum-uncertainty ODE whose solution is the Gaussian.) Numerically, for 300 random 4-component kets the gap and $\langle f|f\rangle\langle h|h\rangle$ agree to $10^{-9}$ with $\langle f|h\rangle=0$; setting $g=(2-0.7i)f$ closes the gap to $10^{-9}$, and adding an orthogonal piece $w$ reopens it by exactly $\langle f|f\rangle\langle w|w\rangle$.

Inside the well $V=0$, so the TISE is $\psi''=-k^2\psi$ with $k=\sqrt{2mE}/\hbar$, giving $\psi=A\sin kx+B\cos kx$. The wall condition $\psi(0)=0$ kills the cosine ($B=0$); $\psi(L)=0$ then forces $\sin kL=0$, i.e. $kL=n\pi$. Therefore

so $E_1:E_2:E_3=1:4:9$, and for $L=1,\ \hbar=m=1$, $E_1=\pi^2/2\approx4.9348$. The finite-difference eigensolver, which only knows $V=0$ and the Dirichlet walls, returns $[4.9348,\,19.7392,\,44.4131]$ — matching infinite_well_energy(n,1.0) and confirming quantization comes purely from confinement.

A bound state has $-V_0<E<0$: inside it oscillates with $l=\sqrt{2m(E+V_0)}/\hbar$, outside it decays with $\kappa=\sqrt{-2mE}/\hbar$. Matching $\psi$ and $\psi'$ at $x=a$ (using parity) gives the transcendental pair $\tan z=\sqrt{z_0^2-z^2}/z$ (even) and $-\cot z=\sqrt{z_0^2-z^2}/z$ (odd), where $z=la$ and $z_0=\frac a\hbar\sqrt{2mV_0}$. The right-hand side is a quarter-circle of radius $z_0$, so there are finitely many intersections — one even per $(n\pi,n\pi+\tfrac\pi2)$ and one odd per $(n\pi+\tfrac\pi2,(n{+}1)\pi)$ — totalling $\lceil 2z_0/\pi\rceil$. With $V_0=15,a=1$: $z_0=\sqrt{30}\approx5.477$, so $2z_0/\pi\approx3.49$ and $\lceil\cdot\rceil=4$. finite_square_well_bound_count(15,1) returns $4$, and the FD solver finds exactly four levels, all with $-15<E<0$.

Away from the origin $V=0$, so a bound state ($E<0$) is the decaying $\psi=Be^{-\kappa|x|}$ with $\kappa=\sqrt{-2mE}/\hbar$. Integrating the TISE across $x=0$ over $[-\varepsilon,\varepsilon]$, the $V$ term gives $-\alpha\psi(0)$ and the kinetic term gives the slope jump, so

Here $\psi'(0^\pm)=\mp\kappa B$, so $-2\kappa B=-\frac{2m\alpha}{\hbar^2}B$, fixing $\kappa=m\alpha/\hbar^2$ — one value, hence one bound state. Then $E=-\hbar^2\kappa^2/2m=-m\alpha^2/2\hbar^2$, i.e. $E=-\alpha^2/2$ in natural units, independent of $\alpha$. For $\alpha=2$, $E=-2$: delta_well_energy(2.0) returns $-2.0$, and the narrow-well limit recovers it.

With $E<V_0$ the wave is evanescent inside the barrier, $\psi\sim e^{\pm\kappa x}$, $\kappa=\sqrt{2m(V_0-E)}/\hbar$. Matching $\psi,\psi'$ at $x=0$ and $x=a$ gives the standard result

nonzero because $\sinh$ never vanishes — the particle tunnels. For $V_0=10,E=5$: $\kappa=\sqrt{2\cdot5}=\sqrt{10}\approx3.162$, and $T(a{=}1)\approx7.14\times10^{-3}$. For a thick barrier $\sinh\kappa a\to\tfrac12 e^{\kappa a}$, so $T\to\frac{16E(V_0-E)}{V_0^2}e^{-2\kappa a}\propto e^{-2\kappa a}$: each unit of width multiplies $T$ by $e^{-2\kappa}=e^{-2\sqrt{10}}\approx1/560$. Indeed transmission_barrier(5,10,a) gives $7.14\times10^{-3},\,1.28\times10^{-5},\, 2.30\times10^{-8}$ for $a=1,2,3$ — a $\sim$560-fold drop per unit width.

For $E>V_0$ the inside wave is oscillatory, $k_2=\sqrt{2m(E-V_0)}/\hbar$, and matching at both faces replaces $\sinh$ by $\sin$:

This dips below $1$ in general but equals $1$ exactly when $\sin(k_2a)=0$, i.e. $k_2a=n\pi$ — the barrier holds an integer number of internal half-wavelengths and becomes transparent. Squaring, $\frac{2m(E-V_0)}{\hbar^2}a^2=n^2\pi^2$, so

For $V_0=10,a=1,n=1$: $E_1=10+\pi^2/2\approx14.9348$, and transmission_barrier(10+np.pi**2/2,10,1) returns $T=1.0$, while off-resonance $T<1$ — the Ramsauer–Townsend resonance.

Write $\psi=Ae^{ik_1x}+Be^{-ik_1x}$ for $x<0$ and $\psi=Fe^{ik_2x}$ for $x>0$, with $k_1=\sqrt{2mE}/\hbar$, $k_2=\sqrt{2m(E-V_0)}/\hbar$. Continuity of $\psi$ and $\psi'$ at $x=0$ gives $A+B=F$ and $k_1(A-B)=k_2F$, whence

(a) Then $R=|B/A|^2=\big(\tfrac{k_1-k_2}{k_1+k_2}\big)^2$ and, current-weighted, $T=\tfrac{k_2}{k_1}|F/A|^2=\tfrac{4k_1k_2}{(k_1+k_2)^2}$; adding, $R+T=\tfrac{(k_1-k_2)^2+4k_1k_2}{(k_1+k_2)^2}=1$. For $E=1.5,V_0=1$: $k_1=\sqrt3,k_2=1$, giving $R\approx0.072$, $T\approx0.928$ — partial reflection despite $E>V_0$. (b) For $E<V_0$, $k_2$ is imaginary, the right side is evanescent (no current), so $R=1,T=0$. This matches step_RT(1.5,1.0)$\approx(0.072,0.928)$ and step_RT(0.5,1.0)$=(1,0)$.

From $\omega(k)=\hbar k^2/2m$ the two velocities are

so $v_g=2v_p$, and $v_g=p/m$ is exactly the classical particle speed (group_velocity(3)$=3=2\times$phase_velocity(3)). Because $\omega\propto k^2$ is nonlinear, the component waves dephase and a minimum-uncertainty Gaussian broadens as $\sigma_x(t)=\sigma_0\sqrt{1+(\hbar t/2m\sigma_0^2)^2}$. This is smallest at $t=0$ (where $\sigma_x=\sigma_0$) and grows symmetrically in $|t|$: with $\sigma_0=1$, gaussian_packet_sigma(0,1)$=1<$gaussian_packet_sigma(2,1)$=\sqrt2\approx1.414$, confirming the packet spreads.

Let $N|n\rangle=n|n\rangle$. From $[N,a^\dagger]=a^\dagger$ and $[N,a]=-a$, the kets $a^\dagger|n\rangle$ and $a|n\rangle$ are again $N$-eigenstates, with eigenvalues $n+1$ and $n-1$ — the ladder. Eigenvalues cannot go negative, since

So lowering must terminate at a floor with $a|0\rangle=0$, giving $N|0\rangle=0$ and $E_0=\hbar\omega(0+\tfrac12)=\tfrac12\hbar\omega$. Climbing with $a^\dagger$ adds one quantum at a time, so $n=0,1,2,\dots$ and $E_n=\hbar\omega(n+\tfrac12)$. The floor $\tfrac12\hbar\omega\neq0$ because $a|0\rangle=0$ is a minimum-uncertainty state, not $x=p=0$. Matches algebraic_spectrum(6)→[0.5 1.5 2.5 3.5 4.5 5.5] and zero_point_energy()→0.5.

Since $a^\dagger|n\rangle$ has $N$-eigenvalue $n+1$, write $a^\dagger|n\rangle=c_n|n+1\rangle$. Fix $|c_n|$ from the norm, using $aa^\dagger=a^\dagger a+[a,a^\dagger]=N+1$:

so with the standard real-positive phase $a^\dagger|n\rangle=\sqrt{n+1}\,|n+1\rangle$. Likewise $\lVert a|n\rangle\rVert^2=\langle n|a^\dagger a|n\rangle=n$ gives $a|n\rangle=\sqrt{n}\,|n-1\rangle$, and at $n=0$ the factor $\sqrt0$ kills the state, $a|0\rangle=0$. Hence creation(6) @ basis_vector(2,6) $=\sqrt3\,|3\rangle$ (the $\sqrt3=1.7321$ entry in slot $3$) and annihilation(6) @ basis_vector(0,6) $=0$.

The trace is cyclic for any finite square matrices, $\operatorname{tr}(AB)=\operatorname{tr}(BA)$, so

But $\operatorname{tr}\mathbb 1_D=D\neq0$, so $[A,B]=\mathbb 1$ cannot hold in finite dimensions — it needs the infinite ladder. Truncating $a,a^\dagger$ to $D$ levels keeps $[a,a^\dagger]=1$ on every interior rung, but $a^\dagger|D-1\rangle=\sqrt{D}\,|D\rangle$ is dropped, so the bottom-right entry must equal $-(D-1)$ to force the trace to zero. For $D=6$ this is $\operatorname{diag}(1,1,1,1,1,-5)$, exactly what np.diag(commutator(annihilation(6),creation(6))).real→[1 1 1 1 1 -5] returns, with np.trace(...)$=0$.

The eigenfunction $\psi_k$ has parity $(-1)^k$, so the product $\psi_m\psi_n$ has parity $(-1)^{m+n}$. When $m+n$ is odd the integrand is an odd function and $\int_{-\infty}^\infty\psi_m\psi_n\,dx=0$ with no calculation — this kills $(0,1),(1,2),\dots$. The diagonal terms are the normalization integrals, set to $1$ by the prefactor $(m\omega/\pi\hbar)^{1/4}/\sqrt{2^n n!}$ together with $\int H_n(\xi)^2e^{-\xi^2}\,d\xi=\sqrt\pi\,2^n n!$. Hence $\int\psi_m\psi_n\,dx=\delta_{mn}$, confirmed by overlap(2,2)≈$1.0$ and overlap(1,2)≈$0.0$.

With $H_2(\xi)=4\xi^2-2=2(2\xi^2-1)$, Eq. 2.85 gives

Under $x\to-x$ we have $\xi\to-\xi$, but $\xi^2$ and $e^{-\xi^2/2}$ are unchanged, so $\psi_2(-x)=\psi_2(x)$ — even. Generally $H_n(-\xi)=(-1)^nH_n(\xi)$, so $\psi_n$ has parity $(-1)^n$; for odd $n$ the function is odd and $\psi_n(0)=0$. Hence psi(2,-x) == psi(2,x) while psi(1,-x) == -psi(1,x).

The algebraic route diagonalizes $H=\hbar\omega(N+\tfrac12)$ with $N=\operatorname{diag}(0,1,2,\dots)$, returning $E_n=\hbar\omega(n+\tfrac12)$ exactly. The analytic route discretizes $-\tfrac12 d^2/dx^2+\tfrac12 x^2$ with a 3-point Laplacian and diagonalizes that matrix — no ladder operators anywhere. Both give the same numbers: algebraic_spectrum(20)[:6]$=[0.5,1.5,2.5,3.5,4.5,5.5]$ and fd_spectrum(k=6)$=[0.5,1.500,2.500,3.500,4.499,5.499]$, the tiny deficits being $O(\Delta x^2)$ grid error that vanishes as the box refines. One spectrum, two independent derivations.

Sample the analytic $\psi_n$ on the grid and apply the finite-difference Hamiltonian $H_{\rm fd}=-\tfrac12 D^2+\tfrac12 x^2$ (3-point stencil). Each $\psi_n$ solves the continuum TISE $H\psi_n=E_n\psi_n$ with $E_n=n+\tfrac12$ exactly, so the only error is the stencil's $D^2\psi_n=\psi_n''+O(\Delta x^2)$. The relative residual is then

Numerically fd_residual(0..4)$=\{3.0,\,8.9,\,19,\,34,\,53\}\times10^{-5}$, all $<2\times10^{-3}$ and shrinking as the grid refines — every Hermite–Gaussian solves the discretized Schrödinger equation.

Write $x=\sqrt{\tfrac{\hbar}{2m\omega}}(a+a^\dagger)$ and $p=i\sqrt{\tfrac{\hbar m\omega}{2}}(a^\dagger-a)$. Each shifts $|n\rangle$ to $|n\pm1\rangle$, orthogonal to $|n\rangle$, so $\langle x\rangle=\langle p\rangle=0$. For the squares only the number-conserving terms survive:

giving $\langle x^2\rangle=\tfrac{\hbar}{2m\omega}(2n+1)=\tfrac{\hbar}{m\omega}(n+\tfrac12)$ and likewise $\langle p^2\rangle=\tfrac{\hbar m\omega}{2}(2n+1)=\hbar m\omega(n+\tfrac12)$. With $\langle x\rangle=\langle p\rangle=0$, $\sigma_x\sigma_p=\sqrt{\langle x^2\rangle\langle p^2\rangle}=\hbar(n+\tfrac12)\ge\tfrac\hbar2$, saturated at $n=0$. In natural units $n=2$ gives $\langle x^2\rangle=\langle p^2\rangle=2.5$, matching (v.conj()@(x@x)@v).real→2.5.

Two operators share a complete eigenbasis only if they commute, so test a single common eigenstate: if $|\psi\rangle$ obeyed $L_x|\psi\rangle=a|\psi\rangle$ and $L_y|\psi\rangle=b|\psi\rangle$, then

forcing $L_z|\psi\rangle=0$; cyclically $L_x|\psi\rangle=L_y|\psi\rangle=0$, so $L^2|\psi\rangle=0$ — only the trivial $l=0$ state. Hence (a) no two components share eigenstates. (b) But $[L^2,L_i]=0$, so $L^2$ and one component ($L_z$) do share a basis, labelled $(l,m)$. The code confirms commutator(Lx(1),Ly(1)) $=i\,$Lz(1) and commutator(L_squared(1),Lz(1)) $=0$.

Order the rungs $m=1,0,-1$. Then $L_z=\hbar\,\mathrm{diag}(1,0,-1)$, with spectrum $\{-1,0,1\}$. The ladder element $\hbar\sqrt{l(l+1)-m(m+1)}=\hbar\sqrt{2-m(m+1)}$ equals $\hbar\sqrt2$ for both $m=0$ and $m=-1$, so

$L_x,L_y$ are not diagonal, but their squares conspire: $L^2=L_x^2+L_y^2+L_z^2=\hbar^2l(l+1)\,\mathbb1=2\hbar^2\mathbb1$. So eigvalsh(Lz(1))$=[-1,0,1]$, L_squared(1)$=2I$, and casimir_eigenvalue(1)$=2$.

The raising coefficient is $c_+(m)=\hbar\sqrt{l(l+1)-m(m+1)}$. At the top rung $m=l$,

The lowering coefficient $c_-(m)=\hbar\sqrt{l(l+1)-m(m-1)}$ at $m=-l$ gives $l(l+1)-(-l)(-l-1)=l(l+1)-l(l+1)=0$, so $L_-|l,-l\rangle=0$. Both ends are sealed, leaving exactly $2l+1$ integer-spaced rungs $m=-l,\dots,l$. For $l=2$ this is ladder_coeff(2,2,+1)$=0$ and ladder_coeff(2,-2,-1)$=0$; correspondingly L_plus(2)[:,0] (the $|2,2\rangle$ column) and L_minus(2)[:,-1] (the $|2,-2\rangle$ column) are all-zero.

For $l=\tfrac12$ the rungs are $m=\pm\tfrac12$ (dim 2), so $L_z=\tfrac\hbar2\,\mathrm{diag}(1,-1)$. The single ladder element is $\hbar\sqrt{\tfrac12\cdot\tfrac32-(-\tfrac12)(\tfrac12)}=\hbar\sqrt{\tfrac34+\tfrac14}=\hbar$, giving

Hence $\tfrac2\hbar L_x,\tfrac2\hbar L_y,\tfrac2\hbar L_z=\sigma_x,\sigma_y,\sigma_z$, i.e. $L_i=\tfrac\hbar2\sigma_i$. With $\hbar=1$ the code returns 2Lx(0.5)$=\sigma_x$, 2Ly(0.5)$=\sigma_y$, 2*Lz(0.5)$=\sigma_z$ — the entire spin-$\tfrac12$ content of ~QM-11.

Numerically sphere_inner_product(2,1,2,1)$=1.0000$ and sphere_inner_product(1,0,2,0)$=-8.7\times10^{-17}\approx0$. No integral is needed to know this: $L^2$ and $L_z$ are Hermitian, and eigenfunctions of a Hermitian operator belonging to distinct eigenvalues are orthogonal. $Y_1^0$ and $Y_2^0$ have different $L^2$ eigenvalues ($2\hbar^2$ vs $6\hbar^2$), and any two $Y_l^m$ with $m\neq m'$ have different $L_z$ eigenvalues $\hbar m$. Either mismatch forces the overlap to vanish, which is exactly the $\delta_{ll'}\delta_{mm'}$ the checks return.

Write $Y_l^m=N\,P_l^m(\cos\theta)\,e^{im\phi}$; only the $e^{im\phi}$ factor depends on $\phi$. Then

So $Y_l^m$ is an $L_z$ eigenfunction with eigenvalue $\hbar m$ — the very number sitting on the diagonal of the matrix $L_z$ in P2; the two are one operator in two representations. For $(l,m)=(3,-2)$ the eigenvalue is $-2\hbar$, and the finite-difference Lz_on_Y(3,-2,0.6,0.9)$=0.1222+0.5238i$ matches $-2\cdot$spherical_harmonic(3,-2,0.6,0.9)$=0.1222+0.5238i$.

The algebra of P1–P3 only needs the ladder to terminate, which it does for every $l=0,\tfrac12,1,\tfrac32,\dots$. But a spatial wavefunction must be single-valued: the azimuthal factor $\Phi(\phi)=e^{im\phi}$ has to satisfy $\Phi(\phi+2\pi)=\Phi(\phi)$, i.e. $e^{2\pi i m}=1$, which forces $m\in\mathbb Z$ and hence $l\in\mathbb Z$. The half-integer rungs the algebra permits therefore carry no $Y_l^m$ — they are intrinsic spin (~QM-11), matrices with no $\mathbf r\times\mathbf p$ realization. Accordingly spherical_harmonic(0.5,0.5,…) raises ValueError, while Lz(0.5)$=\tfrac12\mathrm{diag}(1,-1)$ and L_squared(0.5)$=\tfrac34\,\mathbb1$ exist and satisfy $[L_x,L_y]=i\hbar L_z$.

The length of $\mathbf L$ is fixed by $L^2=\hbar^2l(l+1)$, so $|\mathbf L|=\hbar\sqrt{l(l+1)}$; for $l=2$, $|\mathbf L|=\hbar\sqrt6\approx2.449\hbar$, whereas the largest projection is $L_z^{\max}=\hbar l=2\hbar$. Since $l(l+1)>l^2$ for any $l>0$, always $\sqrt{l(l+1)}>l$. Equality would demand $L_x=L_y=0$ with $L_z$ sharp — impossible, because $[L_x,L_y]=i\hbar L_z\neq0$ (P1) forbids that simultaneous eigenstate. Only as $l\to\infty$ does $\sqrt{l(l+1)}/l\to1$, recovering the classical vector. Hence math.sqrt(casimir_eigenvalue(2))$=2.449$ strictly exceeds max(m_values(2))$=2$.

Write the identity once and again with $i\leftrightarrow j$, using $\delta_{ji}=\delta_{ij}$ and $\varepsilon_{jik}=-\varepsilon_{ijk}$:

Adding cancels the antisymmetric term: $\{\sigma_i,\sigma_j\}=\sigma_i\sigma_j+\sigma_j\sigma_i=2\delta_{ij}\mathbb 1$, so $\sigma_i^2=\mathbb 1$. Subtracting cancels the $\delta$ term: $[\sigma_i,\sigma_j]=\sigma_i\sigma_j-\sigma_j\sigma_i=2i\varepsilon_{ijk}\sigma_k$. For $i=x,j=y$ this gives $\sigma_x\sigma_y=\delta_{xy}\mathbb 1+i\varepsilon_{xyz}\sigma_z=i\sigma_z$, exactly the sigma_x@sigma_y $=i\,$sigma_z checked in test_pauli_product_identity.

With $\mathbf S=\tfrac\hbar2\boldsymbol\sigma$ the spin commutator inherits the Pauli commutator of P1:

and cyclically — the orbital algebra $[L_i,L_j]=i\hbar\varepsilon_{ijk}L_k$ verbatim. Using $\sigma_i^2=\mathbb 1$,

Comparing with $S^2=\hbar^2 s(s+1)\mathbb 1$ gives $s(s+1)=\tfrac34$, i.e. $s=\tfrac12$ (since $\tfrac12\cdot\tfrac32=\tfrac34$). This is ~QM-10's ladder at its half-integer rung, confirmed by test_spin_algebra and test_S_squared_is_three_quarter.

Build $\hat{\mathbf n}\cdot\boldsymbol\sigma=n_x\sigma_x+n_y\sigma_y+n_z\sigma_z$; the off-diagonal entries combine as $n_x\mp in_y=\sin\theta\,e^{\mp i\phi}$, so

Since $(\hat{\mathbf n}\cdot\boldsymbol\sigma)^2=n_in_j\sigma_i\sigma_j=|\hat{\mathbf n}|^2\mathbb 1=\mathbb 1$ (using P1), the eigenvalues are $\pm1$, i.e. $\pm\hbar/2$ for every direction. The top row of $(\hat{\mathbf n}\cdot\mathbf S-\tfrac\hbar2)\chi=0$ reads $(\cos\theta-1)a+\sin\theta\,e^{-i\phi}b=0$, giving $b/a=\tfrac{1-\cos\theta}{\sin\theta}e^{i\phi}=\tan\tfrac\theta2\,e^{i\phi}$, hence $\chi_+=(\cos\tfrac\theta2,\ \sin\tfrac\theta2\,e^{i\phi})^T$. This is the spinor for which spin_operator_along(θ,φ) @ spin_eigenstate(θ,φ) $=0.5\,$spin_eigenstate(θ,φ).

Write $\chi_+=(a,b)$ with $a=\cos\tfrac\theta2$, $b=\sin\tfrac\theta2\,e^{i\phi}$. For the diagonal $\sigma_z$, $\langle\sigma_z\rangle=|a|^2-|b|^2=\cos^2\tfrac\theta2-\sin^2\tfrac\theta2=\cos\theta$, so $\langle S_z\rangle=\tfrac\hbar2\cos\theta$. For the off-diagonal pair, $\chi^\dagger\sigma_x\chi=2\,\mathrm{Re}(a^b)$ and $\chi^\dagger\sigma_y\chi=2\,\mathrm{Im}(a^b)$ with $a^*b=\cos\tfrac\theta2\sin\tfrac\theta2\,e^{i\phi}$:

using $2\cos\tfrac\theta2\sin\tfrac\theta2=\sin\theta$. Thus $\langle\mathbf S\rangle=\tfrac\hbar2\hat{\mathbf n}$ — length $\hbar/2$, pointing along $\hat{\mathbf n}$ (the Bloch sphere). At $\theta=1,\phi=2$ the code returns $\langle\mathbf S\rangle=(-0.1751,\,0.3826,\,0.2702)=0.5\,$n_hat(1.0,2.0), matching test_expectation_is_half_n.

The prepared state is $\chi_+=(\cos\tfrac\theta2,\ \sin\tfrac\theta2\,e^{i\phi})$ from P3, and the $S_z$ "up" eigenstate is $\uparrow_z=(1,0)$. The Born rule gives the amplitude as the inner product

and likewise $P(\downarrow_z)=|\sin\tfrac\theta2\,e^{i\phi}|^2=\sin^2\tfrac\theta2$. The phase $e^{i\phi}$ drops out of the modulus, so the result is $\phi$-independent, and $\cos^2\tfrac\theta2+\sin^2\tfrac\theta2=1$. Only two outcomes occur — the $2s+1=2$ Stern–Gerlach spots. For $\theta=1$, $\cos^2(0.5)=0.7702$, which is what prob_up_z(spin_eigenstate(1.0,2.0)) and born_probabilities(Sz, …) return ($\{+\hbar/2:0.7702,\,-\hbar/2:0.2298\}$).

With $\mathbf B=B_0\hat{\mathbf z}$, $H=-\gamma B_0 S_z=-\tfrac{\gamma B_0\hbar}{2}\sigma_z$, so the propagator is diagonal: $U(t)=e^{-iHt/\hbar}=\mathrm{diag}(e^{i\omega t/2},e^{-i\omega t/2})$ with $\omega\equiv\gamma B_0$. Acting on the tilted initial spinor $\chi(0)=(\cos\tfrac\alpha2,\ \sin\tfrac\alpha2)$,

This is the P3/P5 spinor with $\theta=\alpha$ and azimuth $\phi=-\omega t$. Substituting into $\langle\mathbf S\rangle=\tfrac\hbar2\hat{\mathbf n}$:

So $\langle\mathbf S\rangle$ precesses on a cone of half-angle $\alpha$ at $\omega=\gamma B_0$; for $\gamma=2,B_0=3$, larmor_frequency(2,3) $=6$, matching test_larmor_precession_closed_form.

In the rotating frame $H=\tfrac\hbar2(\Delta\sigma_z+\Omega\sigma_x)$ has eigenvalues $\pm\tfrac\hbar2\Omega_R$ with $\Omega_R=\sqrt{\Omega^2+\Delta^2}$. Evolving $|\!\downarrow\rangle$ and projecting on $|\!\uparrow\rangle$ gives the Rabi formula

On resonance $\Delta=0\Rightarrow\Omega_R=\Omega$, so $P_\uparrow=\sin^2(\Omega t/2)$ reaches $1$ at $t=\pi/\Omega$ (a $\pi$-pulse): rabi_probability(np.pi,1,0) $=1.0$. Off resonance the prefactor caps the flop at $\Omega^2/\Omega_R^2<1$; for $\Omega=1,\Delta=1.5$ the maximum is $1/(1+1.5^2)=1/3.25=0.3077$, matching test_rabi_detuned_amplitude.

Because $(\hat{\mathbf n}\cdot\boldsymbol\sigma)^2=\mathbb 1$, the exponential splits into even and odd powers exactly like $e^{-i\theta/2}$ for a number:

The half-angle is the whole story. At $\theta=2\pi$: $\cos\pi=-1,\ \sin\pi=0$, so $R(2\pi)=-\mathbb 1$; at $\theta=4\pi$: $\cos2\pi=1,\ \sin2\pi=0$, so $R(4\pi)=+\mathbb 1$. A $360^\circ$ rotation flips the spinor's sign — it takes $720^\circ$ to return — because spinors carry a representation of $SU(2)$, the double cover of $SO(3)$ (~MA-18). This is exactly np.allclose(rotation(2np.pi), -I2) and np.allclose(rotation(4np.pi), I2).

The effective radial potential is $V_{\text{eff}}(r)=-\frac1r+\frac{l(l+1)}{2r^2}$. Compare the two pieces as $r\to0$: the centrifugal term $\frac{l(l+1)}{2r^2}\sim r^{-2}$ blows up faster than the Coulomb well $-\frac1r\sim r^{-1}$, so for any $l>0$ it dominates and $V_{\text{eff}}\to+\infty$ — a repulsive wall. Keeping the leading terms in $u''=[l(l+1)/r^2-\dots]u$ near the origin forces $u\sim r^{l+1}$, i.e. $R=u/r\sim r^l$, which is suppressed more steeply the larger $l$. Numerically at $r=0.01$: $l{=}1$ gives $V_{\text{eff}}=9900>0$ (barrier) but $l{=}0$ gives $-100<0$ (pure attraction); at $r=0.1$ the values climb $290>90>-10$ for $l=2,1,0$ — exactly what effective_potential returns.

Setting $u=rR$ turns the radial equation into $u''=[l(l+1)/r^2-2/r-2E]u$. A 3-point stencil for $-\tfrac12 u''$ on $r_i=i\,\Delta r$ with Dirichlet walls $u(0)=u(r_{\max})=0$ makes $H$ symmetric tridiagonal, $H_{ii}=1/\Delta r^2+V_{\text{eff}}(r_i)$, $H_{i,i\pm1}=-1/(2\Delta r^2)$. Diagonalizing for $l=0$ gives eigenvalues $E_k$ with $n=k+1$. The Bohr prediction $E_n=-1/(2n^2)$ is $-0.5,-0.125,-0.0\overline{5}$ Ha, and the solver returns radial_solve(0,3)[0]$=[-0.49999,-0.12500,-0.05556]$. The ground state in eV is $-0.5\times27.2114=-13.61$ eV (radial_solve(0,1)[0][0]*HARTREE_EV$\approx-13.61$) — the $-13.6/n^2$ ladder recovered purely numerically, no closed form used.

Series termination caps $l$ at $n-1$, so the allowed orbitals are $l=0,1,\dots,n-1$ ($n$ of them), each carrying $2l+1$ magnetic substates $m=-l,\dots,+l$. Summing the arithmetic series,

For $n=3$: $l\in\{0,1,2\}$ gives $1+3+5=9=3^2$. The $2l+1$ factor is generic — the rotational $L_z$ degeneracy of any central potential (~QM-10); the extra collapse of different $l$ to one energy is the accidental $1/r$ degeneracy, lifted in ~QM-17. Matches allowed_l(3)$=[0,1,2]$ and count_states(3)$=9=$degeneracy(3).

For $n=1,l=0$ the Laguerre degree is $n-l-1=0$ and $L_0^{(\alpha)}=1$, so Eq. 4.89 collapses to $R_{10}(r)=\sqrt{(2/1)^3\,\frac{0!}{2\cdot1!}}\;e^{-r}=\sqrt{4}\,e^{-r}=2e^{-r}$, giving $R_{10}(0)=2$. Normalize with $\int_0^\infty r^k e^{-ar}\,dr=k!/a^{k+1}$:

So $R_{10}=2e^{-r}$ is correctly normalized — radial_wavefunction(1,0,0.0)$=2.0$ and radial_norm(1,0)$\approx1$.

The radial function is $R_{nl}\propto (2r/n)^l e^{-r/n}\,L_{n-l-1}^{2l+1}(2r/n)$; the prefactor $(2r/n)^l e^{-r/n}$ never changes sign, so every interior zero comes from the Laguerre polynomial, which has degree $n-l-1$ and exactly that many positive real roots. Hence the radial node count is $n-l-1$: for $n=3$, $3s\to3{-}0{-}1=2$, $3p\to1$, $3d\to0$. All three share the energy $E_3=-1/18$ Ha but differ in node structure — the visible fingerprint that distinguishes the degenerate orbitals (the $l>0$ vanishing at $r=0$ is a boundary, not a node). Matches count_radial_nodes(3,0)=2, (3,1)=1, (3,2)=0.

The ground-state radial density is $P(r)=r^2|R_{10}|^2=4r^2e^{-2r}$. Its mode solves $P'(r)=0$: $\frac{d}{dr}(r^2e^{-2r})=2r\,e^{-2r}(1-r)=0\Rightarrow r=1=a_0$, the most probable radius. The mean weights by $r$:

The mean exceeds the mode because $P$ is right-skewed with a long large-$r$ tail. In general $\langle r\rangle_{nl}=\tfrac12[3n^2-l(l+1)]$, e.g. $(2,1)\to\tfrac12[12-2]=5$. Confirmed: expectation_r(1,0)$\approx1.5$, the peak of radial_probability(1,0,r) at $r\approx1.0$, and expectation_r_closed(2,1)$=5$.

The generalized Laguerre polynomial is $L_k^{(\alpha)}(x)=\sum_{j=0}^{k}(-1)^j\binom{k+\alpha}{k-j}\frac{x^j}{j!}$. At $k=0$ only the $j=0$ term survives, $\binom{\alpha}{0}=1$, so $L_0^{(\alpha)}(x)=1$ for every $\alpha$ — therefore each $n=l+1$ state ($1s,2p,3d,\dots$) has the nodeless form $R\propto r^l e^{-r/n}$. Setting $\alpha=0$ reproduces the ordinary Laguerre $L_k=L_k^{(0)}$ that ~MA-12 builds from its recurrence. Hence generalized_laguerre(0,3,2.34)$=1$, and generalized_laguerre(k,0,x) equals MA-12's laguerre(k,x).

The radial operator $H_l=-\tfrac12\frac{d^2}{dr^2}+V_{\text{eff}}(r)$ (acting on $u=rR$ on $L^2(0,\infty)$, equivalently $H$ acting on $R$ under the weight $r^2\,dr$) is Hermitian. For fixed $l$ the $R_{nl}$ are its eigenfunctions with distinct eigenvalues $E_n=-1/2n^2$, and eigenfunctions of a Hermitian operator belonging to different eigenvalues are automatically orthogonal — no integral required, the same argument that orthogonalizes ~QM-10's $Y_l^m$. With the P4 normalization this gives $\int_0^\infty R_{nl}R_{n'l}\,r^2dr=\delta_{nn'}$. Numerically the $(2,0)$–$(3,0)$ overlap is $\approx-2\times10^{-11}\approx0$ while $n=n'$ returns $\approx1$, matching the _trapz(...) check.

The four product states $|\!\uparrow\uparrow\rangle,|\!\uparrow\downarrow\rangle,|\!\downarrow\uparrow\rangle,|\!\downarrow\downarrow\rangle$ carry $M=m_1+m_2=+1,0,0,-1$. A single spin-$J$ multiplet must contain each $M$ from $-J$ to $J$ exactly once in unit steps; here $M=+1$ and $M=-1$ appear once each, so the largest multiplet is $J=1$ with $M=1,0,-1$ — and that uses up only one of the two $M=0$ states. The leftover $M=0$, with no partner at $M=\pm1$, can only be a $J=0$ singlet. Hence

This is exactly decomposition(0.5,0.5)→'1/2 (x) 1/2 = 1 (+) 0' with dimension_check(0.5,0.5)→(4,4).

For spin $\tfrac12$ the lowering operator gives $J_{i-}|\!\uparrow\rangle=\hbar|\!\downarrow\rangle$ (since $\hbar\sqrt{s(s+1)-m(m-1)}=\hbar$ at $s=m=\tfrac12$), so $J_-=J_{1-}+J_{2-}$ acts as

On the coupled side $J_-|1,1\rangle=\hbar\sqrt2\,|1,0\rangle$, so $|1,0\rangle=\tfrac1{\sqrt2}(|\!\uparrow\downarrow\rangle+|\!\downarrow\uparrow\rangle)$; one more lowering gives $|1,-1\rangle=|\!\downarrow\downarrow\rangle$. The singlet is the remaining unit vector of the $M=0$ sector orthogonal to $|1,0\rangle$, $|0,0\rangle=\tfrac1{\sqrt2}(|\!\uparrow\downarrow\rangle-|\!\downarrow\uparrow\rangle)$. In the uncoupled order $|\!\uparrow\uparrow\rangle,|\!\uparrow\downarrow\rangle,|\!\downarrow\uparrow\rangle,|\!\downarrow\downarrow\rangle$ these are the columns coupled_basis(0.5,0.5) returns: $[1,0,0,0]$, $[0,\tfrac1{\sqrt2},\tfrac1{\sqrt2},0]$, $[0,0,0,1]$, $[0,\tfrac1{\sqrt2},-\tfrac1{\sqrt2},0]$.

The series runs from the "parallel" maximum $J=j_1+j_2=\tfrac32+2=\tfrac72$ down to the "antiparallel" minimum $J=|j_1-j_2|=\tfrac12$ in unit steps, each value once:

Coupling only re-labels the same Hilbert space, so the two basis counts must agree:

This is multiplet_content(1.5,2.0)→[3.5,2.5,1.5,0.5] and dimension_check(1.5,2.0)→(20,20).

Expand $J^2=(\mathbf J_1+\mathbf J_2)^2=J_1^2+J_2^2+2J_{1z}J_{2z}+J_{1+}J_{2-}+J_{1-}J_{2+}$. The flip-flop term $J_{1\pm}J_{2\mp}$ sends $|m_1,m_2\rangle\to|m_1{\mp}1,m_2{\pm}1\rangle$ — same $M=m_1+m_2$, but a different product state — so $J^2$ picks up off-diagonal entries inside each equal-$M$ block. For $\tfrac12\otimes\tfrac12$ (order $|\!\uparrow\uparrow\rangle,|\!\uparrow\downarrow\rangle,|\!\downarrow\uparrow\rangle,|\!\downarrow\downarrow\rangle$),

whose central ($M=0$) block is non-diagonal. Thus $[J^2,J_{1z}]\ne0$, while the total $[J^2,J_z]=0$. Numerically J_squared(0.5,0.5) shows the off-diagonal $1$'s, commutator($J^2$,J1z) has largest entry $1\ne0$, and commutator($J^2$,Jz)$=0$.

In $1\otimes\tfrac12$ the $M=\tfrac12$ sector holds exactly two product states, $|1,-\tfrac12\rangle$ and $|0,\tfrac12\rangle$ (both with $m_1+m_2=\tfrac12$); they split one each into the $J=\tfrac32$ and $J=\tfrac12$ multiplets. Reading Table 4.8,

orthogonal unit vectors, the Condon–Shortley sign making the largest-$m_1$ amplitude positive. Thus cg_coefficient(1,1,0.5,-0.5,1.5,0.5)$=\sqrt{1/3}=0.5774$ and the $|0,\tfrac12\rangle$ amplitude of $|\tfrac12,\tfrac12\rangle$, cg_coefficient(1,0,0.5,0.5,0.5,0.5)$=-\sqrt{1/3}=-0.5774$, both entries of cg_table(1,0.5).

Each coupled vector is built inside a single $(j_1,j_2)$ product space, so it is automatically an eigenstate of the Casimirs $J_1^2\to\hbar^2j_1(j_1+1)$ and $J_2^2\to\hbar^2j_2(j_2+1)$. It is assembled within one $M$-sector (eigenstate of $J_z\to\hbar M$) and forced to be an eigenstate of $J^2\to\hbar^2J(J+1)$ by the highest-weight$+$ladder construction. Since $\{J^2,J_z,J_1^2,J_2^2\}$ mutually commute a common eigenbasis exists, and the change of basis diagonalizes the mixed operator: $U^\dagger J^2U=\mathrm{diag}\,\hbar^2J(J+1)$. For $1\otimes\tfrac12$ this returns $\mathrm{diag}(3.75,3.75,3.75,3.75,0.75,0.75)$ — $J(J+1)=3.75$ for the $J=\tfrac32$ quartet and $0.75$ for the $J=\tfrac12$ doublet — exactly what U.conj().T @ J_squared @ U gives.

Both bases are orthonormal, so the change of basis $U$ (entries the CG coefficients) is unitary: $U^\dagger U=UU^\dagger=\mathbb1$. The column condition $\sum_{m_1,m_2}|\langle j_1m_1;j_2m_2|JM\rangle|^2=1$ just says each coupled state $|J,M\rangle$ is normalized. The row condition $\sum_{J,M}|\langle j_1m_1;j_2m_2|JM\rangle|^2=1$ is the probabilistic statement: prepare the product state $|m_1,m_2\rangle$, measure total $(J,M)$, and the possible outcomes occur with probabilities $|\langle j_1m_1;j_2m_2|JM\rangle|^2$ that must sum to one. Hence U.conj().T @ U$=I$ and every row (and column) sum of $|U|^2$ equals $1$ (e.g. $[1,1,1,1]$ for $\tfrac12\otimes\tfrac12$).

The highest-weight$+$ladder recursion that builds $U$ is never graded against itself: every entry of cg_table(1,1) is compared to sympy.physics.wigner.clebsch_gordan (the same Condon–Shortley convention) and agrees to machine precision. For instance the $1\otimes1$ singlet comes out as

each coefficient $\pm1/\sqrt3=\pm0.57735$ matching sympy. The identical coupling reappears as $\mathbf J=\mathbf L+\mathbf S$ in fine structure (~QM-17, where $\mathbf L\!\cdot\!\mathbf S$ is diagonal in the coupled basis) and as the maximally entangled Bell singlet $(|\!\uparrow\downarrow\rangle-|\!\downarrow\uparrow\rangle)/\sqrt2$ in ~QM-21.

With orthonormal $|a\rangle,|b\rangle$ form $\Psi_\pm=\tfrac1{\sqrt2}(|a\rangle|b\rangle\pm|b\rangle|a\rangle)$. The exchange operator swaps the two slots, $P_{12}|a\rangle|b\rangle=|b\rangle|a\rangle$ and $P_{12}|b\rangle|a\rangle=|a\rangle|b\rangle$, so

Orthonormality also fixes the norm, $\langle\Psi_\pm|\Psi_\pm\rangle=\tfrac12(1+1)=1$. Thus the symmetric state is a $+1$ eigenstate (boson) and the antisymmetric a $-1$ eigenstate (fermion): P@symmetrize(a,b) $=$ symmetrize(a,b) and P@antisymmetrize(a,b) $=$ -antisymmetrize(a,b).

Set $\psi_b=\psi_a$ in the two combinations. The antisymmetric one has two identical tensor terms that cancel exactly,

"no wave function at all" — two identical fermions cannot share a state, the Pauli exclusion principle. The symmetric one does not cancel: $\Psi_+=A\cdot2\,|a\rangle|a\rangle$, which normalizes to the plain product $|a\rangle\otimes|a\rangle$, so any number of bosons may pile into one state. Hence norm(antisymmetrize(a,a)) $=0$ while norm(symmetrize(a,a)) $=1$.

Expand the squared norm of the raw combination using $\langle a|a\rangle=\langle b|b\rangle=1$:

(a) Orthonormal ($s=0$): the norm is $\sqrt2$, so $A=1/\sqrt2$. (b) Bosons with $\psi_a=\psi_b$: the symmetric state is $2A\,|a\rangle|a\rangle$ of norm $2|A|=1$, so $A=\tfrac12$. In general $A=1/\sqrt{2(1\pm|s|^2)}$. Numerically the orthonormal raw norm is $\sqrt2=1.4142$, while for $b=(1,1)/\sqrt2$ (so $|s|^2=\tfrac12$) the sum has norm $\sqrt3=1.7321$ and the difference norm $1$.

The completely antisymmetric three-fermion state is the antisymmetrizer applied to the product,

the tensor-space form of the $3\times3$ Slater determinant $\det[\phi_i(\mathbf r_j)]$. Interchanging particles $p$ and $q$ permutes two columns of that determinant, which flips its sign, so $P_{pq}\Psi=-\Psi$ for every pair. With $\phi_i=e_i$ the code confirms apply_pair_swap(Psi,3,3,p,q) $=-\Psi$ for all three pairs $(0,1),(0,2),(1,2)$.

If two orbitals coincide the determinant has two equal rows, and a determinant with two equal rows is zero — Pauli exclusion for $N$ fermions. For the prefactor, expand the raw antisymmetric sum's norm with orthonormal orbitals: cross terms between distinct permutations vanish (each carries a factor $\langle\phi_i|\phi_j\rangle=0$), leaving one surviving term per permutation,

so the raw norm is $\sqrt{N!}$ and dividing by it normalizes the state. Hence norm(slater_determinant([e[0],e[1],e[0]])) $=0$ and the raw norm norm(...,normalize_state=False) $=\sqrt{3!}=\sqrt6=2.449$.

For the well $\langle x\rangle_n=\tfrac L2$ and $\langle x^2\rangle_n=L^2(\tfrac13-\tfrac1{2n^2\pi^2})$, while the off-diagonal element is

The distinguishable value is $\langle x^2\rangle_1+\langle x^2\rangle_2-2\langle x\rangle_1\langle x\rangle_2=L^2(\tfrac16-\tfrac5{8\pi^2})\approx0.1033$, and the exchange term is $2|\langle x\rangle_{12}|^2=\tfrac{512L^2}{81\pi^4}\approx0.0649$. Bosons subtract it (closer, $\approx0.0384$), fermions add it (farther, $\approx0.1682$). This matches r["distinguishable"] $\approx0.1033$, r["exchange_term"] $\approx0.0649$, with r["boson"] < r["distinguishable"] < r["fermion"].

In natural units $\langle x\rangle_n=0$ and $\langle x^2\rangle_n=n+\tfrac12$, so $\langle x^2\rangle_0=\tfrac12$ and $\langle x^2\rangle_1=\tfrac32$. Writing $x=\tfrac1{\sqrt2}(a+a^\dagger)$ gives the one nonzero off-diagonal element

Then distinguishable $=\langle x^2\rangle_0+\langle x^2\rangle_1-0=2$ and the exchange term $2|\langle x\rangle_{01}|^2=1$, so boson $=2-1=1$ and fermion $=2+1=3$. Confirmed by r["distinguishable"] $\approx2$, r["boson"] $\approx1$, r["fermion"] $\approx3$.

The entire identical-particle effect rides on the off-diagonal element $\langle x\rangle_{ab}=\int x\,\psi_a^(x)\psi_b(x)\,dx$. If the orbitals have disjoint support — $\psi_a$ on $[0,1]$, $\psi_b$ on $[2,3]$ — then the product $\psi_a^(x)\psi_b(x)=0$ at every $x$, so the integrand vanishes identically and $\langle x\rangle_{ab}=0$. The exchange term $2|\langle x\rangle_{ab}|^2$ is then $0$ and boson $=$ distinguishable $=$ fermion: spatially separated identical particles behave as distinguishable. The code gives abs(r["exchange_term"]) < 1e-10 (here exactly $0$), with the three $\langle(\Delta x)^2\rangle$ values equal.

The full two-electron state must be antisymmetric. In $1s^2$ both electrons share the same spatial orbital, so the spatial part $1s(1)\,1s(2)$ is symmetric, $P_{12}=+1$. To make the product antisymmetric the spin part must carry $P_{12}=-1$, i.e. the antisymmetric singlet $\tfrac1{\sqrt2}(|\!\uparrow\downarrow\rangle-|\!\downarrow\uparrow\rangle)$; the product of exchange eigenvalues is $(+1)(-1)=-1$. "Both spin up" is the symmetric $|\!\uparrow\uparrow\rangle$, which paired with the symmetric spatial state is overall symmetric — forbidden; equivalently antisymmetrizing two identical $|\!\uparrow\rangle$ spins gives the zero vector. Hence P@s $=-s$, norm(antisymmetrize(spin_up(), spin_up())) $=0$, and each triplet() state satisfies P@t $=+t$.

Insert the series into $H|\psi_n\rangle=E_n|\psi_n\rangle$ and collect $O(\lambda)$: $H^0|\psi_n^{(1)}\rangle+H'|\psi_n^0\rangle=E_n^0|\psi_n^{(1)}\rangle+E_n^{(1)}|\psi_n^0\rangle$. Project onto $\langle\psi_n^0|$; since $H^0$ is Hermitian, $\langle\psi_n^0|H^0|\psi_n^{(1)}\rangle=E_n^0\langle\psi_n^0|\psi_n^{(1)}\rangle$ cancels the matching term on the right, leaving $E_n^{(1)}=\langle\psi_n^0|H'|\psi_n^0\rangle$. With $x=\tfrac1{\sqrt2}(a+a^\dagger)$ the number-conserving terms give $\langle n|x^2|n\rangle=\tfrac12\langle n|(aa^\dagger+a^\dagger a)|n\rangle=\tfrac12(2n+1)=n+\tfrac12$, while $\langle n|x|n\rangle=0$ because $x$ sends $|n\rangle\to|n\pm1\rangle$. At $n=2$ these are $2.5$ and $0$ — matching first_order_energy(...,ho_matrix_power(40,2),2)→2.5 and (...,ho_position(40),2)→0.0.

Projecting the $O(\lambda^2)$ equation onto $\langle\psi_n^0|$ and substituting $|\psi_n^{(1)}\rangle=\sum_{m\neq n}\frac{\langle m|H'|n\rangle}{E_n^0-E_m^0}|m\rangle$ gives $E_n^{(2)}=\sum_{m\neq n}|\langle m|H'|n\rangle|^2/(E_n^0-E_m^0)$. For $H'=x^2=\tfrac12(a+a^\dagger)^2$ the only off-diagonal links are $\langle n{+}2|x^2|n\rangle=\tfrac12\sqrt{(n{+}1)(n{+}2)}$ and $\langle n{-}2|x^2|n\rangle=\tfrac12\sqrt{n(n{-}1)}$, with denominators $\mp2$:

For $n=2$ this is $-1.25$. For the ground state every $E_m^0>E_0^0$ makes each denominator negative while $|\langle m|H'|n\rangle|^2\ge0$, so $E_0^{(2)}\le0$ — the levels repel. Matches second_order_energy(...,ho_matrix_power(60,2),2)→-1.25.

The exact eigenvalue expands as $E_n=E_n^0+\lambda E_n^{(1)}+\lambda^2 E_n^{(2)}+\lambda^3 E_n^{(3)}+\cdots$. Truncating after first order leaves the remainder $\lambda^2 E_n^{(2)}+O(\lambda^3)=O(\lambda^2)$; through second order leaves $\lambda^3 E_n^{(3)}+O(\lambda^4)=O(\lambda^3)$. So halving $\lambda$ should divide the first-order error by $2^2=4$ and the through-second-order error by $2^3=8$. Diagonalizing $\operatorname{diag}(E^0)+\lambda x^2$ exactly at $\lambda=0.04,0.02,0.01$ gives first-order errors $\{3.85,0.98,0.25\}\times10^{-4}$ (ratios $\times3.9,\times4.0$) and through-second errors $\{1.52,0.20,0.025\}\times10^{-5}$ (ratios $\times7.8,\times7.9$) — exactly the $O(\lambda^2)$ and $O(\lambda^3)$ scalings the test__error_is_order_lambda_ checks assert.

Complete the square: $H=\tfrac{p^2}{2}+\tfrac12(x+\lambda)^2-\tfrac12\lambda^2$ is a harmonic oscillator recentred at $x=-\lambda$, so the exact spectrum is $E_n=(n+\tfrac12)-\tfrac12\lambda^2$. Perturbatively, $H'=x$ gives $E_n^{(1)}=\langle n|x|n\rangle=0$ by parity, and only $m=n\pm1$ contribute to second order, with $|\langle n{\pm}1|x|n\rangle|^2= \tfrac12(n{+}1),\tfrac12 n$ and denominators $\mp1$:

The shift $-\tfrac12\lambda^2$ is reproduced already at second order and every higher correction vanishes (the square completes exactly), so PT is exact here. At $\lambda=0.3$, $\tfrac12\lambda^2=0.045$, so perturbed_energy(...,n,0.3,order=2)$=(n+\tfrac12)-0.045$, i.e. $0.455,1.455,2.455$ for $n=0,1,2$.

Ordinary PT puts $E_n^0-E_m^0=0$ in the denominator for the two degenerate states, so it diverges; the cure is that $H'$ itself selects the "good" zeroth-order states. Projecting the $O(\lambda)$ equation into the degenerate subspace turns the first-order problem into the eigenvalue equation $\det(W-E^{(1)}\mathbb 1)=0$ with $W_{ij}=\langle i|H'|j\rangle$. For $W=\begin{pmatrix}0.3&0.4\\0.4&-0.2\end{pmatrix}$,

The off-diagonal $0.4$ pushes the eigenvalues outside the diagonal pair $\{W_{00},W_{11}\}= \{0.3,-0.2\}$ (eigenvalue interlacing), so the naive answer $[-0.2,0.3]$ is wrong. degenerate_first_order_split(Hp,[0,1])→[-0.422, 0.522] matches $(E^{\text{exact}}-E^0)/\lambda$ as $\lambda\to0$.

Expand any normalized trial state in the (unknown) energy eigenbasis, $\psi=\sum_n c_n\psi_n$ with $\sum_n|c_n|^2=1$; then $\langle H\rangle=\sum_n|c_n|^2E_n\ge E_{gs}\sum_n|c_n|^2=E_{gs}$ since every $E_n\ge E_{gs}$, with equality only for $\psi=\psi_{gs}$. For the Gaussian on $V=\tfrac12 x^2$,

The bound is saturated because at $b=\tfrac12$ the trial $\psi_b\propto e^{-x^2/2}$ is the exact ground state $\psi_0$. Hence gaussian_variational_min(lambda t: 0.5*t**2, x)→(0.5, 0.5): minimum energy $0.5=E_0$ at width $b=0.5$.

For the Gaussian, $\langle T\rangle=\tfrac{\hbar^2 b}{2m}=\tfrac{b}{2}$ and $\langle x^4\rangle=3\langle x^2\rangle^2=3/(16b^2)$, so on $V=\tfrac14 x^4$

giving $\langle H\rangle_{\min}\approx0.4293$. Unlike the oscillator the Gaussian is not the true quartic ground state, so this is a strict over-estimate: the finite-difference truth is $E_{gs}\approx0.4208$, about $2\%$ below ($0.4293\ge0.4208$). The inequality in fact holds at every width $b$. So gaussian_variational_min(lambda t: 0.25*t**4, x)[0]$=0.4293\ge$ fd_ground_energy(...)$=0.4208$.

With $p=\sqrt{2(E-\tfrac12 x^2)}$ the turning points are $x_0=\pm\sqrt{2E}$ and the action is the area of a momentum ellipse,

exactly the oscillator spectrum. WKB drops the $A''/A$ term; for the quadratic potential this does not disturb the quantization condition (the result is exact), but for $V=\tfrac14 x^4$ it does, leaving an $O(1/n)$ semiclassical error. Numerically bohr_sommerfeld_energy(lambda t:0.5*t**2, n, -80, 80) returns $0.5,1.5,2.5,3.5$ (exact), while the quartic error falls from $18\%$ at $n=0$ to $1.3\%$ ($n=1$) and $0.6\%$ ($n=2$) — WKB is a large-$n$ method.

In the forbidden region the momentum magnitude is constant, $|p|=\sqrt{2m(V_0-E)}\equiv\hbar\kappa$ with $\kappa=\sqrt{2(V_0-E)}=4$ here, so $\gamma=\kappa a$ and $T_{\text{WKB}}=e^{-2\kappa a}$. QM-08's exact $T=[1+\tfrac{V_0^2\sinh^2\kappa a}{4E(V_0-E)}]^{-1}$ uses $\sinh\kappa a\to\tfrac12 e^{\kappa a}$ for a thick barrier, giving

Both decay with the same exponent $e^{-2\kappa a}$; WKB misses only the constant $O(1)$ prefactor. tunneling_rectangular(2,10,a) and transmission_barrier(2,10,a) give ratio $2.56000$ for all thick $a$ ($a=3,5,8$) — the exponent that powers Gamow's alpha-decay law.

The integral is elementary:

Factor $e^{i\omega_{fi}t}-1=e^{i\omega_{fi}t/2}\,(e^{i\omega_{fi}t/2}-e^{-i\omega_{fi}t/2})=e^{i\omega_{fi}t/2}\,2i\sin(\omega_{fi}t/2)$, so $|e^{i\omega_{fi}t}-1|^2=4\sin^2(\omega_{fi}t/2)$ and

It oscillates in $t$ and is largest for small $\omega_{fi}$ (near-degenerate levels), valid only while $P\ll1$. With $\hbar=1$, first_order_probability(0.013, 0.7, 9.0)$=9.75\times10^{-8}$, matching $0.013^2\cdot4\sin^2(0.7\cdot9/2)/0.7^2$ to the integrator tolerance.

(a) Setting $\delta=0$ in the generalized Rabi formula, the cap $\Omega_R^2/(\Omega_R^2+\delta^2)\to1$ and $\sqrt{\Omega_R^2+\delta^2}\to\Omega_R$, so $P=\sin^2(\Omega_R t/2)$. (b) This first equals $1$ when $\Omega_R t/2=\pi/2$, i.e. $t=\pi/\Omega_R$ — a $\pi$-pulse that completely inverts the population. (c) On resonance $H=\tfrac12\Omega_R\sigma_x$, so

and $|\langle f|\psi(t)\rangle|^2=\sin^2(\Omega_R t/2)$ — identical to (a). Hence rabi_probability(np.pi, 1.0, 0.0)$=1.0$ and two_level_exact_probability(np.pi, 0.0, 1.0)$=1.0$, full inversion at $t=\pi$.

(a) The flopping factor $\sin^2(\cdots)$ ranges over $[0,1]$, so the maximum of $P_{i\to f}$ over $t$ is just the prefactor,

reached whenever $\sqrt{\Omega_R^2+\delta^2}\,t/2=\tfrac\pi2+n\pi$. The detuning steals amplitude — you can never deposit the system in $|f\rangle$ with certainty. (b) The sine's argument advances at the generalized Rabi frequency $\sqrt{\Omega_R^2+\delta^2}$, faster than the bare $\Omega_R$. For $\Omega_R=1,\delta=3$ the cap is $1/(1+9)=1/10$; a fine-grid scan of rabi_probability(t,1,3) peaks at $0.0999999\approx0.1$, and generalized_rabi_frequency(1,3)$=\sqrt{10}=3.162$.

With constant rotating-frame coupling $V_{fi}=\Omega_R/2$ and Bohr frequency $\omega_{fi}=-\delta$, the P1 result gives $P^{(1)}=(\Omega_R/2)^2\,4\sin^2(\delta t/2)/\delta^2=(\Omega_R^2/\delta^2)\sin^2(\delta t/2)$. (a) For $\Omega_R\ll|\delta|$ the exact cap $\Omega_R^2/(\Omega_R^2+\delta^2)\to\Omega_R^2/\delta^2$ and the flopping frequency $\sqrt{\Omega_R^2+\delta^2}\to|\delta|$, so the exact Rabi probability collapses onto $P^{(1)}$. (b) For short times $\sin^2x\to x^2$, and both become

the $\delta$-independent quadratic onset. Numerically first_order_probability(0.005, -0.6, 7.5)$=1.68\times10^{-4}\approx$ rabi_probability(7.5, 0.01, 0.6)$=1.68\times10^{-4}$, and at $\Omega_R=0.5,t=0.02$ both give $(0.5\cdot0.02/2)^2=2.5\times10^{-5}$.

Write $d=\omega_{fi}-\omega$. (a) As $\omega\to\omega_{fi}$, $d\to0$ and $\sin^2(dt/2)/d^2\to(t/2)^2$, so the line peaks at resonance $\omega=\omega_{fi}$ with height

(b) The profile is $\mathrm{sinc}^2(dt/2)$: it vanishes when $dt/2=n\pi$, i.e. first zeros at $\omega=\omega_{fi}\pm2\pi/t$, so the width $\propto1/t$ while the peak $\propto t^2$ — drive longer and the line grows taller and sharper (spectroscopic resolution). For ws=np.linspace(0.7,1.3,4001), ws[np.argmax(sinusoidal_probability(40.0, ws, 1.0, 0.01))]$=1.0000$ and the peak height is $(0.01\cdot40/2)^2=0.04$; going $t=20\to40$ multiplies it by $(40/20)^2=4$.

The honest lab-frame coupling is the full cosine, $\langle f|H'|i\rangle=V_{fi}\cos\omega t'=\tfrac12V_{fi}(e^{i\omega t'}+e^{-i\omega t'})$, so the first-order amplitude keeps both the co- and counter-rotating terms (denominators $\omega_{fi}\mp\omega$),

evolve_tdse integrates $i\partial_t|\psi\rangle=H(t)|\psi\rangle$ exactly — no RWA, no weak-coupling assumption — so for $V_{fi}\ll1$ the population stays $\ll1$ and the neglected $O(V^2)$ orders are tiny, leaving the honest dynamics and first order in agreement. Here evolve_tdse(make_driven_H(0.9, 1.0, 0.01), [1,0], 20.0) gives $|\psi_1|^2=7.14\times10^{-3}$ against first_order_probability(...)$=7.17\times10^{-3}$ — matching to $\sim0.4\%$, the residual being the dropped Bloch–Siegert physics.

(a) The probability carries $|V_{fi}|^2=|\langle f|H'|i\rangle|^2$ with $H'=-qE_0z\cos\omega t$ Hermitian, so $|V_{ba}|=qE_0|\langle b|z|a\rangle|=qE_0|\langle a|z|b\rangle|=|V_{ab}|$. Absorption ($a\to b$) and stimulated emission ($b\to a$) thus ride the same resonance line of §4, with identical peak height $|V_{fi}|^2t^2/4$ — a resonant photon is exactly as likely to drive the atom down (emitting a second, coherent photon) as up. (b) The net upward rate $\propto(N_a-N_b)$: absorption consumes photons while stimulated emission creates them, so amplification (gain $>$ loss) needs $N_b>N_a$, a population inversion. Operationally this is just that $|V_{fi}|^2$ in sinusoidal_probability is invariant under $a\leftrightarrow b$, so the up and down probabilities coincide.

(a) Summing the §4 sinc² over a continuum of density $\rho$, the narrow resonant slice dominates and the rate settles to Fermi's golden rule

(b) Short times: the band-summed first-order probability $P_{\text{band}}(t)=\sum_k|c_k^{(1)}|^2$ grows linearly with slope $\Gamma$ (first order, no saturation); intermediate times: the exact survival $P_i(t)=|\langle i|e^{-iHt}|i\rangle|^2\approx e^{-\Gamma t}$ — the linear law is the short-time face of the exponential. (c) $\Gamma\propto|V|^2$ (doubling $V$ gives $\times4$) and $\propto\rho$ (doubling $\rho$ gives $\times2$). Here golden_rule_rate(0.0309, 50.0)$=2\pi\cdot0.0309^2\cdot50\approx0.3000$, and both the band_first_order_probability slope and the exact band decay reproduce it.

To first order (and in the golden rule) the transition probability is proportional to $|\langle f|H'|i\rangle|^2$; for dipole radiation $H'\propto\mathbf r$, so the rate carries $|\langle f|\mathbf r|i\rangle|^2$. Since $\mathbf r$ is odd under parity, $\langle f|\mathbf r|i\rangle=0$ unless $\psi_f,\psi_i$ have opposite parity ($\Delta\ell$ odd), and the angular (Wigner–Eckart) integral sharpens this to $\Delta\ell=\pm1$, $\Delta m=0,\pm1$. "Forbidden" therefore means literally $\langle f|H'|i\rangle=0$, hence $P_{i\to f}=0$ to this order — no drive strength causes the jump. In code this is the trivial limit first_order_probability(0.0, w, t)$=0.0$ and golden_rule_rate(0.0, rho)$=0.0$.

The combination $\alpha=e^2/4\pi\varepsilon_0\hbar c$ is a pure number: $e^2/4\pi\varepsilon_0$ carries energy$\times$length, $\hbar$ action, $c$ speed, and the units cancel exactly, leaving

The Bohr energy is this $\alpha$ squared below the rest energy, $\mathrm{Ry}=\tfrac12\alpha^2 m_ec^2=\tfrac12(7.297\times10^{-3})^2(511\,\mathrm{keV})=13.606$ eV. Fine structure carries a further $\alpha^2\sim5\times10^{-5}$, so it drops from $13.6$ eV to $\sim10^{-4}$ eV. Matches 1/fine_structure_constant() → $137.036$ and rydberg_energy_eV() → $13.6057$.

With $E^1_{\rm fs}=E_n\dfrac{\alpha^2}{n^2}\!\left[\dfrac{n}{j+\frac12}-\dfrac34\right]$ and $E_2=-3.4$ eV, the bracket is $2/1-\tfrac34=\tfrac54$ for $j=\tfrac12$ and $2/2-\tfrac34=\tfrac14$ for $j=\tfrac32$, so

Dividing by $h$ gives $f=\Delta E/h\approx10.9$ GHz. This is exactly fine_structure_correction_eV(2,1.5)-fine_structure_correction_eV(2,0.5) → $4.528\times10^{-5}$ eV, and $\times e/h/10^9$ → $10.95$ GHz (the lab value is $\sim10.97$).

The grand result $E^1_{\rm fs}(n,j)$ carries no $\ell$. For $2S_{1/2}$ ($\ell=0$) spin-orbit vanishes, so the shift is relativistic $+$ Darwin, $(-1.472+0.906)\times10^{-4}=-5.66\times10^{-5}$ eV; for $2P_{1/2}$ ($\ell=1$) the Darwin term vanishes, so it is relativistic $+$ spin-orbit, $(-0.264-0.302)\times10^{-4}=-5.66\times10^{-5}$ eV — identical, though assembled from different mechanisms. The $\ell$-degeneracy is therefore broken only through $j$: $2P_{3/2}$ sits higher at $-1.13\times10^{-5}$ eV. Hence fine_structure_correction_eV(2,0.5) → $-5.660\times10^{-5}$ (matching both $\ell$-decompositions) $<$ fine_structure_correction_eV(2,1.5) → $-1.132\times10^{-5}$.

From $\langle\mathbf L\!\cdot\!\mathbf S\rangle=\tfrac{\hbar^2}{2}[j(j{+}1)-\ell(\ell{+}1)-s(s{+}1)]$ with $\ell=1,\ s=\tfrac12$ (so $\ell(\ell{+}1)=2,\ s(s{+}1)=\tfrac34$):

Weighting by the multiplicities $2j+1=2,4$: $2(-\hbar^2)+4(\tfrac12\hbar^2)=0$, so spin-orbit only redistributes the doublet about its center of gravity (it adds no net energy). Matches LS_coupling(0.5,1) $=-1$, LS_coupling(1.5,1) $=0.5$, and the weighted sum $=0$ (in units of $\hbar^2$).

With $g_J=1+\dfrac{j(j{+}1)+s(s{+}1)-\ell(\ell{+}1)}{2j(j{+}1)}$ and $s=\tfrac12$:

They bracket $1$ and $2$ because only the projection of $\mathbf S$ on the conserved $\mathbf J$ survives the precession: $s=0$ leaves the orbital value $1$, $\ell=0$ leaves the spin value $2$. Matches lande_g_factor → $2,\ \tfrac23,\ \tfrac43,\ \tfrac65$.

The internal (spin-orbit) scale in $n=2$ is the P2 interval $\Delta E_{\rm fs}\approx4.5\times10^{-5}$ eV; the external Zeeman energy is $\mu_B B$ with $\mu_B=5.79\times10^{-5}$ eV/T. They become comparable when

so the weak-field Landé picture gives way to Paschen–Back only at several tesla. At $B=1$ T both limits happen to return $\mu_B B$ here: zeeman_weak_field_eV(0.5,0,0.5,1.0)$=\mu_B g_J B m_j=\mu_B(2)(\tfrac12)=5.79\times10^{-5}$ eV and zeeman_strong_field_eV(0,0.5,1.0)$=\mu_B B(m_\ell+2m_s)=\mu_B(1)=5.79\times10^{-5}$ eV — both the order of bohr_magneton_eV_per_T()*1.0, confirming the crossover near $1$ T.

Evaluating the closed form with CODATA constants,

(microwave/radio). The gap is the triplet–singlet spin-spin splitting: $\langle\mathbf S_p\!\cdot\!\mathbf S_e\rangle=\tfrac{\hbar^2}{2}[F(F{+}1)-\tfrac32]$ gives $+\tfrac14\hbar^2$ for $F{=}1$ (triplet, raised) and $-\tfrac34\hbar^2$ for $F{=}0$ (singlet, depressed), a gap of $\hbar^2$. Matches hyperfine_splitting_eV()1e6 → $5.877$, hyperfine_frequency()/1e6 → $1421$, hyperfine_wavelength()100 → $21.1$, and spin_spin_coupling(1)-spin_spin_coupling(0) → $1$ ($\hbar^2$).

The three rungs are

The first step is $E^1_{\rm fs}/E_1=\tfrac14\alpha^2\approx1.3\times10^{-5}$ — the $\alpha^2$ of fine structure; the second, $\Delta E_{\rm hf}/E^1_{\rm fs}\approx0.03$, is $\sim m_e/m_p$ ($=5.4\times10^{-4}$) enhanced by $g_p$ and numerical factors. The Lamb shift ($\sim\alpha^3$, QED, ~QM-22) sits between fine and hyperfine. The ladder bohr_energy_eV(1), fine_structure_correction_eV(1,0.5), hyperfine_splitting_eV() returns $-13.606$, $-1.811\times10^{-4}$, $5.877\times10^{-6}$ eV.

The incident plane wave $Ae^{ikz}$ carries probability current $j_{\rm inc}=|A|^2v$ with $v=\hbar k/m$. The scattered wave $Af(\theta)e^{ikr}/r$ has outward radial current $j_r=|A|^2|f|^2v/r^2$, so the rate into solid angle $d\Omega$ is $j_r\,r^2\,d\Omega=|A|^2v\,|f|^2\,d\Omega$. Dividing the scattered rate by the incident flux $|A|^2v$ gives the area

The $1/r^2$ in $j_r$ exactly cancels the $r^2$ in the spherical area element — which is why the $1/r$ amplitude was forced by flux conservation. differential_cross_section returns $|$partial_wave_amplitude$|^2$ identically, and Gauss–Legendre integration of it over $4\pi$ gives $7.860153=$total_cross_section(k,deltas).

(a) Classically every ray with impact parameter $b<a$ is deflected and every $b>a$ misses, so the cross-section is the geometric disc $\sigma_{\rm cl}=\pi a^2$. (b) At low energy only $\ell=0$ survives ($\sin^2\delta_\ell\sim(ka)^{2(2\ell+1)}$ suppresses $\ell\ge1$), and $\delta_0=-ka$, so

The ratio is exactly $4$: a wave of wavelength $\gg a$ diffracts around the entire sphere, not just the head-on disc. total_cross_section(k, phase_shifts_hard_sphere(k,1.0,8))/(pi*1.0**2) climbs $3.98682$ ($ka{=}0.1$) $\to3.99880$ ($ka{=}0.03$) $\to3.99997$ ($ka{=}0.005$), tending to $4$.

Outside the core the $\ell=0$ reduced radial wave is free, $u_0(r)\propto\sin(kr+\delta_0)$. A hard core imposes the Dirichlet condition $u_0(a)=0$, so $ka+\delta_0=0\pmod\pi$, giving

Equivalently $\tan\delta_0=j_0(ka)/n_0(ka)=\dfrac{\sin ka/ka}{-\cos ka/ka}=-\tan ka$, the same $-ka$. The wave function is just the free sine pushed outward by the radius $a$. hard_sphere_phase_shift(0, ka, 1.0) returns exactly $-ka$ ($-0.010000$ at $ka{=}0.01$, $-0.200000$ at $ka{=}0.2$).

Inside the well, regularity at the origin forces $R=j_\ell(k_{\rm in}r)$, $k_{\rm in}=\sqrt{k^2+V_0}$ (with $\hbar=2m=1$); outside, $R=\cos\delta_\ell\,j_\ell(kr)-\sin\delta_\ell\,n_\ell(kr)$. Continuity of the logarithmic derivative $rR'/R$ at $r=a$, with $L=(k_{\rm in}a)\,j_\ell'(k_{\rm in}a)/j_\ell(k_{\rm in}a)$, gives

For $\ell=0$ ($j_0=\sin x/x,\ n_0=-\cos x/x$) this collapses to $\delta_0=-ka+\arctan[(k/k_{\rm in})\tan(k_{\rm in}a)]$; as $V_0\to0$, $k_{\rm in}\to k$ and $\delta_\ell\to0$. square_well_phase_shift(0,k,a,V0) matches the closed form in $\sin^2\delta_0$ ($0.575489$ at $k{=}0.5,V_0{=}2$; $0.882053$ at $k{=}1.5,V_0{=}4$), and $V_0\to0$ gives $\sin^2\delta_\ell\to0$ (no scattering).

Set $\theta=0$ in $f(\theta)=\frac1k\sum_\ell(2\ell+1)e^{i\delta_\ell}\sin\delta_\ell P_\ell(\cos\theta)$; every $P_\ell(1)=1$ and $\operatorname{Im}(e^{i\delta_\ell}\sin\delta_\ell)=\sin^2\delta_\ell$, so

Physically, total scattering depletes the forward beam, and that depletion is the destructive interference of the incident wave with the forward-scattered one — which requires $f(0)$ to carry a positive imaginary part. optical_theorem_sigma(k,deltas) and total_cross_section(k,deltas) agree to machine precision, both $=7.860153$ for the $k{=}1.5,V_0{=}4$ well.

For a spherically symmetric $V$ the Born amplitude reduces to $f(\theta)=-\frac{2m}{\hbar^2q}\int_0^\infty rV(r)\sin(qr)\,dr$, $q=2k\sin(\theta/2)$. With $V=\beta e^{-\mu r}/r$ the $r$ cancels, and $\int_0^\infty e^{-\mu r}\sin(qr)\,dr=q/(\mu^2+q^2)$:

Since $q$ grows monotonically from $0$ to $2k$, $|f|$ falls monotonically — forward-peaked — and $\mu>0$ keeps $f(0)=-2m\beta/\hbar^2\mu^2$ finite. The numerical born_amplitude_radial(theta,k,V) equals the closed-form born_amplitude_yukawa(theta,k,beta,mu) ($-2.062907$ at $\theta{=}0.3$, $-0.598600$ at $\theta{=}1.0$).

Drop the screening ($\mu\to0$) in the Yukawa amplitude: $f=-\frac{2m\beta}{\hbar^2q^2}$, so

with $E=k^2$ ($\hbar=2m=1$). Hence $d\sigma/d\Omega=(\beta/4E)^2/\sin^4(\theta/2)$ — the Rutherford formula, identical to the classical ~CM-08 result, the famous agreement of classical mechanics, the Born approximation and QED for $1/r$. rutherford_cross_section(theta,k,beta)*sin(theta/2)**4 is angle-independent at $\beta^2/(16E^2)=0.0625$ ($\beta{=}E{=}1$), and $|$born_amplitude_yukawa(...,mu=1e-6)$|^2=0.465933$ equals rutherford_cross_section at $\theta{=}1.3$.

Each partial wave contributes $\sigma_\ell=\frac{4\pi}{k^2}(2\ell+1)\sin^2\delta_\ell$. Since $\sin^2\delta_\ell\le1$,

the unitarity bound, reached only when $\delta_\ell=\pi/2$ — a resonance, where that partial wave scatters as strongly as conservation of probability allows. partial_cross_section(l,k,delta_l) stays under $4\pi(2\ell+1)/k^2$ for every $\ell$ (e.g. the near-resonant $\ell=1$ wave of the $k{=}1.8,V_0{=}7$ well gives $\sigma_1=11.44140\le11.63553$).

The kernel is $K_0=\sqrt{\tfrac{m}{2\pi i\hbar t}}\,\exp\!\big[\tfrac{im(x-x')^2}{2\hbar t}\big]$. The exponent is purely imaginary, so $|e^{iS/\hbar}|=1$ and the modulus lives entirely in the prefactor; since $|1/\sqrt{i}|=1$,

independent of $x,x'$ — every endpoint is reached with equal amplitude, all the physics is in the phase. That phase is $\tfrac{m(x-x')^2}{2\hbar t}=S_{\rm cl}/\hbar$, the action of the straight-line path. Finally $1/\sqrt{i}=e^{-i\pi/4}$, so the prefactor adds the Maslov phase $-\pi/4$. Numerically free_propagator(np.array([1.0]),-0.5,0.8) $=0.36279+0.25947i$ with $|K_0|=0.44603=1/\sqrt{2\pi\cdot0.8}$ and $\arg K_0+\pi/4=1.40625=$ classical_action_free(1.0,-0.5,0.8).

$K_0(x,t;x')$ is a complex Gaussian in $x-x'$ of width $\sim\sqrt{\hbar t/m}$ with unit integral, $\int K_0\,dx'=1$. Expanding a smooth $g$ about $x$ inside the convolution,

because the Gaussian's second moment is $\propto t$. As $t\to0$ the kernel concentrates at $x'=x$ and the integral returns $g(x)$ — a nascent $\delta(x-x')$. The leading error is $O(t\,g'')$, linear in $t$, so it only shrinks for decreasing $t$; a fixed or growing $t$ would not converge. The check confirms this: at $t=0.2,0.1,0.05,0.02$ the errors are $\{3.07,\,1.54,\,0.77,\,0.31\}\times10^{-2}$ — halving as $t$ halves, monotone down to $<10^{-2}$.

The $y$-integral is Gaussian. Combining the two exponents, $\tfrac{im}{2\hbar}\big[\tfrac{(x-y)^2}{t_2}+\tfrac{(y-x')^2}{t_1}\big]$, and completing the square in $y$ gives a $y$-independent remainder $\tfrac{im(x-x')^2}{2\hbar(t_1+t_2)}$ plus a Gaussian in $y$. Doing the Fresnel integral, its $\sqrt{2\pi i\hbar\,t_1t_2/m(t_1+t_2)}$ exactly cancels one factor in the product of prefactors:

reproducing $K_0(x,t_1+t_2;x')$ — the Gaussian-convolution form of $e^{-i\hat Ht_2}e^{-i\hat Ht_1}=e^{-i\hat H(t_1+t_2)}$. Numerically with $t_1=0.5,t_2=0.3$ the convolution gives $0.36117+0.26176i$ vs the closed form $K_0(\dots,0.8)=0.36279+0.25947i$, a $2.8\times10^{-3}$ residual that is real-time Fresnel-truncation, not physics.

Propagating a Gaussian by $\Psi(x,t)=\int K_0(x,t;x')\Psi(x',0)\,dx'$ is again a Gaussian integral, so $|\Psi(x,t)|^2$ stays Gaussian. The mean momentum $p_0$ rides through unchanged, so the centroid moves ballistically at the group velocity,

the width growing because the packet's momentum spread $\hbar/2\sigma$ disperses it. In natural units $\hbar=m=1$: at $t=1,3$ with $x_0=-3,p_0=1.5,\sigma=1$, packet_center_width returns centres $-1.5,\,+1.5$ ($=x_0+p_0t$) and widths $1.118,\,1.803$, matching $\sigma\sqrt{1+(t/2)^2}$. And propagate(...) agrees with the independent split-step FFT evolve_free_spectral(...) to $\sim10^{-14}$ — two algorithms, one $\Psi(x,t)$.

Mehler's kernel is $K=\sqrt{\tfrac{m\omega}{2\pi i\hbar\sin\omega t}}\,\exp\!\big\{\tfrac{im\omega}{2\hbar\sin\omega t}[(x^2+x'^2)\cos\omega t-2xx']\big\}$. (a) As $\omega\to0$, $\sin\omega t\to\omega t$ and $\cos\omega t\to1$, so the prefactor $\to\sqrt{\tfrac{m}{2\pi i\hbar t}}$ and the bracket $\to(x^2+x'^2)-2xx'=(x-x')^2$, giving the free phase $\tfrac{im(x-x')^2}{2\hbar t}$ — exactly $K_0$. (b) Since $K$ is built from the oscillator's stationary states it solves $i\hbar\partial_tK=(-\tfrac{\hbar^2}{2m}\partial_x^2+\tfrac12 m\omega^2x^2)K$. The checks bear both out: harmonic_propagator(...,omega=1e-5) differs from free_propagator(...) by $4.6\times10^{-12}$, and the finite-difference TDSE residual is $7.0\times10^{-9}$ — zero up to stencil error.

Slicing $[0,t]$ into $N$ steps and inserting completeness at each gives $K=\int\!\prod dx_k\prod_k k_\varepsilon(x_{k+1};x_k)$ with the free short-step kernel $k_\varepsilon$. Each intermediate integral is Gaussian and, by the semigroup law of P3, collapses exactly back to $K_0$ — so for the free particle the sum is exact for every $N$ (path_integral_realtime_free matches $K_0$ to $0$ at $N=1$, to $3\times10^{-3}$ at $N=2$, the residual being finite-$R$ truncation). The trouble is that the real-time integrand has constant modulus $|e^{iS/\hbar}|=1$ and never damps — convergence rests on fragile stationary-phase cancellation (the sign problem). Wick-rotating $t\to-i\tau$ turns $e^{iS/\hbar}$ into the positive, decaying heat kernel $e^{-S_E/\hbar}$, and the Strang–Trotter product $(D_{V/2}FD_{V/2})^N\to e^{-\tau\hat H}$ with error $O(1/N^2)$. For $V=\tfrac12\omega^2x^2$ the errors at $N=1,4,16,64$ are $\{9.5\times10^{-3},6.7\times10^{-4},4.3\times10^{-5},2.7\times10^{-6}\}$, falling $\sim16\times$ per $4\times$ refinement — the advertised $O(1/N^2)$ — to $<10^{-3}$.

With $\psi_n(x)=\sqrt{2/L}\sin(n\pi x/L)$ and $E_n=n^2\pi^2\hbar^2/2mL^2$, write $K=\sum_n\psi_n(x)\psi_n^(x')e^{-iE_nt/\hbar}$. At $t=0$ every phase is $1$ and the sum is the completeness relation $\sum_n\psi_n(x)\psi_n^(x')=\delta(x-x')$, so $\int K\,g\,dx'\to g(x)$ (sharper as $n_{\max}$ grows). For an eigenstate $\psi_m$, orthonormality $\int\psi_n^*\psi_m\,dx'=\delta_{nm}$ collapses the sum to its single $n=m$ term:

a pure phase, and truncation is exact once $m\le n_{\max}$. The check confirms it: for $m=1,2,3$, box_propagator(0.4,Lx,0.5,60) integrated against box_eigenfunction(m,Lx) equals $e^{-iE_m\cdot0.5}\psi_m(0.4)$ to $\sim10^{-15}$ (e.g. $m=1$: $-1.05073-0.83964i$ either way).

With $G^R=\theta(t)K$, the product rule gives $i\hbar\partial_tG^R=i\hbar\,\delta(t)\,K(x,0;x')+\theta(t)\,i\hbar\partial_tK$. For $t>0$ the kernel solves the homogeneous TDSE, $\theta(i\hbar\partial_t-\hat H)K=0$; and $K(x,0;x')=\delta(x-x')$, so

the retarded Green's function (the $\theta$-step's derivative supplies $\delta(t)$, exactly as MA-14's causal Green's function does). Replacing $e^{-iE_nt/\hbar}\to1/E_n$ in the eigen-sum gives $\sum_n\psi_n(x)\psi_n^*(x')/E_n=\langle x|\hat H^{-1}|x'\rangle$. For the box $\hat H=-\tfrac12\partial_x^2=\tfrac12 L_{\rm MA14}$, so $\hat H^{-1}=2L^{-1}$ and the kernel is $2\,x_<(1-x_>)$. The check matches term-by-term: at $(0.3,0.7)$, box_resolvent_static $=2\cdot$green_series to $\lesssim10^{-16}$ and both $\to0.180000=2\cdot(0.3)(1-0.7)=2\,$green_dirichlet.

With $|+\rangle=\tfrac1{\sqrt2}\binom{1}{1}$, the projector is

It is real-symmetric, so $\rho=\rho^\dagger$, and $\mathrm{Tr}\,\rho=\tfrac12(1+1)=1$. Squaring, $\rho^2=\tfrac14\begin{pmatrix}2&2\\2&2\end{pmatrix}=\rho$, so $\rho$ is an idempotent rank-1 projector with eigenvalues $\{1,0\}$ and $\mathrm{Tr}(\rho^2)=\sum_k\lambda_k^2=1$ — pure. This is exactly density_matrix_pure(plus), for which is_density_matrix(rho) is True, purity(rho) $=1$, and is_pure(rho) is True.

Adding the two projectors,

It is Hermitian with $\mathrm{Tr}\,\rho=1$, but $\mathrm{Tr}(\rho^2)=\mathrm{Tr}(\tfrac14 I)=\tfrac12<1$ — not pure. In general, in $\rho$'s eigenbasis $\mathrm{Tr}(\rho^2)=\sum_k\lambda_k^2$; since the $\lambda_k\ge0$ sum to $1$, $\sum_k\lambda_k^2\le(\sum_k\lambda_k)^2=1$, with equality iff a single $\lambda_k=1$ (all others $0$), i.e. $\rho^2=\rho$. Hence density_matrix_mixed([ket0, ket1], [0.5, 0.5]) equals maximally_mixed(2), with purity(rho) $=0.5$ and is_pure(rho) False.

Insert completeness $\sum_n|n\rangle\langle n|=I$ for the pure case:

and linearity of the trace gives the mixture $\mathrm{Tr}(\rho A)=\sum_i p_i\langle\psi_i|A|\psi_i\rangle$. For $|+\rangle$, $\sigma_x|+\rangle=+|+\rangle$, so $\langle\sigma_x\rangle=+1$ (determinate); for $I/2$, $\langle\sigma_x\rangle=\mathrm{Tr}(\tfrac12 I\,\sigma_x)=\tfrac12\mathrm{Tr}\,\sigma_x=0$. These reproduce expectation(density_matrix_pure(plus), sigma_x) $=1.0$ and expectation(maximally_mixed(2), sigma_x) $=0.0$ — sharp in the pure eigenstate, fully random in the mixture.

Expand $\rho_{AB}=|\Phi^+\rangle\langle\Phi^+|=\tfrac12\big(|00\rangle+|11\rangle\big)\big(\langle00|+\langle11|\big)$ and take the partial trace $\langle a|\rho_A|a'\rangle=\sum_b\langle a,b|\rho_{AB}|a',b\rangle$. The diagonal terms survive; the cross terms $|00\rangle\langle11|,\,|11\rangle\langle00|$ vanish because their $B$-factor is $\langle1|0\rangle=0$:

The global state is pure ($S=0$) but $\rho_A$ is maximally mixed ($S=\ln2$): the mixedness is entanglement, not ignorance — we know $|\Phi^+\rangle$ exactly. Hence is_pure(bell) is True while partial_trace(bell, [2,2], keep=0) equals maximally_mixed(2).

The product state $\rho_{AB}=|0\rangle\langle0|\otimes|+\rangle\langle+|$ factorizes, so $\rho_A=\mathrm{Tr}_B\,\rho_{AB}=|0\rangle\langle0|\,\mathrm{Tr}(|+\rangle\langle+|)=|0\rangle\langle0|$ with eigenvalues $\{1,0\}$, giving

For the Bell state $\rho_A=I/2$ has eigenvalues $\{\tfrac12,\tfrac12\}$, so $S_A=-2\cdot\tfrac12\ln\tfrac12=\ln2\approx0.6931$. The reduced-state entropy is therefore the entanglement measure: $0$ for the unentangled product, maximal $\ln2$ for the Bell pair. This matches the code — the product von_neumann_entropy(...) $=0$ and the Bell version $=\ln 2\approx0.6931$.

Differentiate $\rho(t)=U\rho_0U^\dagger$ using $i\hbar\dot U=HU$ (hence $i\hbar\dot U^\dagger=-U^\dagger H$):

(a) $\rho(t)=U\rho_0U^\dagger$ is unitarily similar to $\rho_0$, so it shares the same spectrum; thus $\mathrm{Tr}(\rho^2)=\sum_k\lambda_k^2$ is conserved and a pure state ($\lambda=\{1,0\}$) stays pure. (b) $\dot\rho=-\tfrac{i}{\hbar}[H,\rho]$ vanishes iff $[H,\rho]=0$, i.e. $\rho$ is diagonal in the energy eigenbasis. Numerically purity(evolve(rho,H,t)) $=1$ for a pure $\rho$, and a state built diagonal in $H$'s eigenbasis gives von_neumann_rhs(rho,H) $=0$ (norm $0$) and unchanged evolve.

Closed evolution $\rho(t)=U\rho_0U^\dagger$ keeps the eigenvalues of $\rho_0$ fixed for all $t$. The start $|+\rangle\langle+|$ has spectrum $\{1,0\}$ ($\mathrm{Tr}(\rho^2)=1$), while $I/2$ has spectrum $\{\tfrac12,\tfrac12\}$ ($\mathrm{Tr}(\rho^2)=\tfrac12$); the two spectra differ, so no unitary maps one to the other. (a) Hence $\rho$ can never become $I/2$ — unitary evolution only rotates the off-diagonal coherences, never erases them. (b) Erasing them needs coupling to an environment: the coherences then decay irreversibly — decoherence — described by a Lindblad master equation (~QO-04). So purity(density_matrix_pure(plus)) $=1$ and stays $1$ under any evolve(...), making purity $=0.5$ ($I/2$) unreachable unitarily.

A product state expands as $(a_0|0\rangle+a_1|1\rangle)\otimes(b_0|0\rangle+b_1|1\rangle)=a_0b_0|00\rangle+a_0b_1|01\rangle+a_1b_0|10\rangle+a_1b_1|11\rangle$. Matching $\alpha|01\rangle+\beta|10\rangle$ needs $a_0b_0=0$ and $a_1b_1=0$, yet $a_0b_1=\alpha\neq0$ forces $a_0,b_1\neq0$ and $a_1b_0=\beta\neq0$ forces $a_1,b_0\neq0$ — then $a_0b_0\neq0$, a contradiction. So no factorization exists. Equivalently the concurrence $C=2|ad-bc|=2|0\cdot0-\alpha\beta|=2|\alpha\beta|$ is nonzero. With $(\alpha,\beta)=(0.6,0.8)$, concurrence(np.array([0,a,b,0]))$=0.96>0$ and is_product_state(...) is False.

With $\rho=|\psi\rangle\langle\psi|$, trace out Bob. The product $|0\rangle\otimes|1\rangle$ keeps its factorization, $\rho_A=|0\rangle\langle0|\,\langle1|1\rangle=|0\rangle\langle0|$, so $\mathrm{Tr}(\rho_A^2)=1$ — pure. For the singlet $|\Psi^-\rangle=\tfrac1{\sqrt2}(|01\rangle-|10\rangle)$,

The maximally mixed $\rho_A$ means Alice's qubit alone reads 50/50 in every basis — the global pure state pins down nothing local; that gap is the entanglement. Matches purity(...)$=1.0$ for the product vs $=0.5$ for the singlet, whose $\rho_A=I/2$.

On the singlet the two-spin operator is isotropic, $\langle\Psi^-|\sigma_i\otimes\sigma_j|\Psi^-\rangle=-\delta_{ij}$, so

Parallel detectors ($\theta=0$) give $E=-1$ — whenever Alice reads $+1$, Bob reads $-1$ (perfect anti-correlation); anti-parallel ($\theta=\pi$) give $E=+1$; perpendicular ($\theta=\pi/2$) give $E=0$. Confirmed by E_singlet(ndir(0),ndir(0))$=-1$, E_singlet(ndir(0),ndir(np.pi))$=+1$, and E_singlet(ndir(0),ndir(np.pi/2))$\approx0$.

Using $E(\mathbf a,\mathbf b)=-\cos\theta_{ab}$ with $\mathbf a\perp\mathbf b$ and $\mathbf c$ at $45°$ to each,

The inequality $|P(\mathbf a,\mathbf b)-P(\mathbf a,\mathbf c)|\le1+P(\mathbf b,\mathbf c)$ then reads

which is false: quantum mechanics violates Griffiths' Bell inequality. bell_3setting(ndir(0),ndir(np.pi/2),ndir(np.pi/4)) returns $(0.707,0.293)$ with LHS $>$ RHS.

(a) For a local deterministic assignment $A,A',B,B'\in\{\pm1\}$,

Since $B,B'=\pm1$, exactly one of $B-B'$, $B+B'$ vanishes and the other is $\pm2$, so $|S|=2$ for every strategy; a general LHV theory is a convex mixture of these vertices and stays $|S|\le2$. (b) At $(0°,90°,45°,135°)$ the singlet gives $E(a,b)=E(a',b)=E(a',b')=-\tfrac1{\sqrt2}$ and $E(a,b')=+\tfrac1{\sqrt2}$, so $S=-\tfrac1{\sqrt2}-\tfrac1{\sqrt2}-\tfrac1{\sqrt2}-\tfrac1{\sqrt2}=-2\sqrt2$. Thus lhv_deterministic_strategies() has max $|S|=2$, while chsh_from_angles(singlet(),0,90,45,135)$=-2.828=-2\sqrt2$, exactly tsirelson_bound()$=2.828$.

For isotropic $\boldsymbol\lambda$ the readings $\mathrm{sgn}(\mathbf a\!\cdot\!\boldsymbol\lambda)$ and $\mathrm{sgn}(\mathbf b\!\cdot\!\boldsymbol\lambda)$ disagree with probability $\eta/\pi$ (the fraction of the sphere between the two flip-planes), $\eta$ being the angle between $\mathbf a,\mathbf b$. With Bob's extra minus sign,

linear in $\eta$ rather than the quantum $-\cos\eta$. Being local it must obey $|S|\le2$; at the optimal angles $\eta=45°,135°$ it gives $S=-2$, sitting right at the classical bound and never reaching $2\sqrt2$. Accordingly the Monte-Carlo test_chsh_physical_lhv_model_below_2 returns $|S|\lesssim2$, well under $2.4$.

Partial trace is cyclic within Bob's factor, so for any unitary $U$

A non-selective measurement inserts $\sum_kP_k=\mathbb I$, so $\rho'=\sum_k(\mathbb I\otimes P_k)\rho(\mathbb I\otimes P_k)$ likewise gives $\mathrm{Tr}_B\rho'=\mathrm{Tr}_B[\rho\sum_kP_k]=\rho_A$. Bob cannot touch Alice's marginal — for a Bell state it stays $\mathbb I/2$ whatever he does, so her data alone are completely random and carry no signal. Hence no_signaling_unitary(singlet(),U) and no_signaling_measurement(singlet(),U) return $\rho_A^{\text{before}}=\rho_A^{\text{after}}$ to machine precision.

Write the singlet antisymmetrically, $|\Psi^-\rangle=\tfrac1{\sqrt2}\sum_{ij}\varepsilon_{ij}|i\rangle|j\rangle$ with $\varepsilon_{01}=-\varepsilon_{10}=1$. Under $U\otimes U$ the coefficient of $|k\rangle|l\rangle$ becomes $\sum_{ij}\varepsilon_{ij}U_{ki}U_{lj}=\varepsilon_{kl}\det U$, so

The singlet is the rotational scalar $J=0$ — the same from every orientation, which is exactly why "up along $\hat{\mathbf n}$" for Alice forces "down along $\hat{\mathbf n}$" for Bob on any axis. The symmetric triplet $|\Phi^+\rangle$ ($J=1$) instead picks up relative phases $e^{\mp i\theta}$ on $|00\rangle,|11\rangle$ under a $z$-rotation. Hence kron(U,U) @ singlet() equals singlet(), while for bell_state("Phi+") it does not.

Acting on $\phi=e^{-ip\cdot x}$, each $\partial_\mu$ pulls down $-ip_\mu$, so $\Box\phi=\partial_\mu\partial^\mu\phi=(-ip_\mu)(-ip^\mu)\phi=-(p\cdot p)\,\phi$ and

which vanishes iff $E^2=\vec p^{\,2}+m^2$. This is quadratic in $E$, so $E=\pm\sqrt{\vec p^{\,2}+m^2}$ — both roots solve KG and the spectrum is unbounded below. The conserved density $\rho=\frac{i}{2m}(\phi^\partial_t\phi-\phi\,\partial_t\phi^)=\frac{E}{m}|\phi|^2$ flips sign with $E$, so the negative-energy root gives $\rho<0$ — not a probability. With $\vec p=(0.6,0,0.8)$, $m=1$, $E=\sqrt2$; kg_operator_fd returns $\approx2.9\times10^{-7}$ (the finite-difference floor) for both the $+1$ and $-1$ signs — i.e. $\approx0$, confirming each root sits on shell.

Set $\mu=\nu=0$: $\{\gamma^0,\gamma^0\}=2(\gamma^0)^2=2g^{00}\mathbb1=+2\mathbb1\Rightarrow(\gamma^0)^2=+\mathbb1$. Set $\mu=\nu=i$: $2(\gamma^i)^2=2g^{ii}\mathbb1=-2\mathbb1\Rightarrow(\gamma^i)^2=-\mathbb1$. For $\mu\ne\nu$, $g^{\mu\nu}=0$ forces

Ordinary numbers always commute ($ab=ba$), so objects that anticommute yet square to the nonzero values $\pm\mathbb1$ cannot be scalars — they must be matrices, and the smallest faithful realization of this Clifford algebra is $4\times4$. The Dirac representation satisfies it for every pair: clifford(0,0)$=+2\mathbb1_4$ (entry $2.0$), clifford(1,1)$=-2\mathbb1_4$ (entry $-2.0$), clifford(0,1)$=0$ — all 16 pairs equal $2g^{\mu\nu}\mathbb1_4$.

Because $p_\mu p_\nu$ is symmetric, only the symmetric (anticommutator) part of $\gamma^\mu\gamma^\nu$ survives:

Hence $(\not{\!p}-m)(\not{\!p}+m)=\not{\!p}^{\,2}-m^2\mathbb1=(p\cdot p-m^2)\mathbb1_4$ — the Klein–Gordon operator. The eigenvalues of $\not{\!p}$ are $\pm\sqrt{p\cdot p}$ (each doubly degenerate), so $\det(\not{\!p}-m)=(\sqrt{p\cdot p}-m)^2(\sqrt{p\cdot p}+m)^2=(p\cdot p-m^2)^2$, vanishing iff $p\cdot p=m^2$, i.e. $E^2=\vec p^{\,2}+m^2$. For $\vec p=(0.3,0.4,0)$, $m=1$, on shell $p\cdot p=1.0=m^2$, so dirac_squared$\approx8\times10^{-17}$ and dirac_determinant$\approx1\times10^{-32}$ (both $\approx0$); off shell the determinant reproduces $(p\cdot p-m^2)^2$ exactly (e.g. $1.1881$).

In the Dirac rep $\not{\!p}=\begin{pmatrix}E\,\mathbb1&-\vec\sigma\cdot\vec p\\\vec\sigma\cdot\vec p&-E\,\mathbb1\end{pmatrix}$. The lower block of $(\not{\!p}-m)u=0$ gives $\chi=\frac{\vec\sigma\cdot\vec p}{E+m}\varphi$, so with $\varphi=\xi^s$ and $N=\sqrt{E+m}$,

With $\bar w=w^\dagger\gamma^0$ and $(\vec\sigma\cdot\vec p)^2=\vec p^{\,2}$, $\bar uu=N^2\big[1-\frac{p^2}{(E+m)^2}\big]=\frac{(E+m)^2-p^2}{E+m}=2m$ (using $(E+m)^2-p^2=2m(E+m)$); swapping upper/lower blocks flips the bracket, so $\bar vv=-2m$. The minus sign marks the negative-energy/antiparticle (positron) sector. Numerically dirac_solutions([.3,.4,0],1.0) gives $\bar uu_0=2.0$, $\bar vv_0=-2.0$, $u_0^\dagger u_0=2.236=2E$, $\bar u_0v_0=0$, and dirac_matrix(p4,1.0)@u0$\approx0$.

Apply $(\vec\sigma\cdot\vec a)(\vec\sigma\cdot\vec b)=\vec a\cdot\vec b+i\vec\sigma\cdot(\vec a\times\vec b)$ with $\vec a=\vec b=\vec\pi$. The kinetic momenta no longer commute: $(\vec\pi\times\vec\pi)_z=[\pi_x,\pi_y]=iqB$ (and cyclic), so

Dividing by $2m$, the Pauli Hamiltonian splits into an orbital piece $\vec\pi^{\,2}/2m$ and a spin piece $H_{\text{spin}}=-\frac{q}{2m}\vec\sigma\cdot\vec B$. Writing $\vec\sigma=2\vec S$ recasts this as $-\frac{q}{2m}g\,\vec S\cdot\vec B$ with $\boxed{g=2}$ — twice the orbital $g=1$, and spin was never inserted by hand. dirac_g_factor() returns $(2.0,\,1.55\times10^{-15})$, the tiny residual confirming the operator identity $(\vec\sigma\cdot\vec\pi)^2-\vec\pi^{\,2}=-q\,\vec\sigma\cdot\vec B$ holds exactly.

Factor out the rest mass and expand the square root for $|\vec p|\ll m$:

The three terms are the rest energy $m$, the Schrödinger kinetic energy $p^2/2m$, and the leading relativistic correction $-p^4/8m^3$. Treated as a perturbation in hydrogen, this last term is one of the two contributions to the fine structure (~QM-17); the other, spin–orbit coupling, likewise drops out of the Dirac reduction. For $\vec p=(0.2,0.1,0)$, $m=1$, energy_expansion(...,order=2)$=1.0246875$ lies within $7.6\times10^{-6}$ of the exact $m+$nonrel_energy$=1.0246951$ — far closer than the order=1 (Schrödinger) value $1.025$, whose error is $3.0\times10^{-4}$.