Transfer matrices for piecewise-constant potentials

Set V(x) = Vⱼ on the interval [xⱼ, xⱼ₊₁]. There H|ψ⟩ = E|ψ⟩, written in the position basis as ψ(x) = ⟨x|ψ⟩, reduces to ψ″ = −kⱼ²ψ with kⱼ = √(2m(E − Vⱼ))/ħ constant, whose general solution is Aⱼe^{(ikⱼ(x−xⱼ))} + Bⱼe^{(−ikⱼ(x−xⱼ))}. Two constants, so the state inside that slab is a vector (Aⱼ, Bⱼ) in ℂ² — the two-fold degeneracy of the continuum, written as a basis. Anchoring the exponentials at the slab's own left edge rather than at a global origin is a bookkeeping choice, and it is the one that keeps the matrices free of stray phases. When E < Vⱼ nothing in the algebra changes: kⱼ = iκⱼ with κⱼ = √(2m(Vⱼ − E))/ħ real, and the right-goer becomes the decaying exponential while the left-goer becomes the growing one. For an electron ħ²/2mₑ = 3.810 × 10⁻² eV nm², so k = 5.123 √(E/eV) nm⁻¹: one electronvolt buys about 5.1 inverse nanometres.

At the wall xⱼ₊₁ the potential jumps but stays finite, so ψ and ψ′ are both continuous there. Let a = Aⱼe^(ikⱼdⱼ) and b = Bⱼe^{(−ikⱼdⱼ)} be the amplitudes carried across the slab of width dⱼ = xⱼ₊₁ − xⱼ. Continuity of ψ gives Aⱼ₊₁ + Bⱼ₊₁ = a + b; continuity of ψ′ gives kⱼ₊₁(Aⱼ₊₁ − Bⱼ₊₁) = kⱼ(a − b). Add and subtract, and with s = kⱼ/kⱼ₊₁ you get Aⱼ₊₁ = ½[(1+s)a + (1−s)b] and Bⱼ₊₁ = ½[(1−s)a + (1+s)b]. That is the interface matrix I = ½[[1+s, 1−s], [1−s, 1+s]], with det I = ¼[(1+s)² − (1−s)²] = s. The carrying step is the diagonal Pⱼ = diag(e^(ikⱼdⱼ), e^{(−ikⱼdⱼ)}), determinant 1, which in a classically forbidden slab reads diag(e^{(−κⱼdⱼ)}, e^(+κⱼdⱼ)). Every piecewise-constant problem in one dimension is built from those two matrices and nothing else.

Compose them: M = I_(N−1→N) P_(N−1) ⋯ P₁ I_(0→1), so that (A_R, B_R)ᵀ = M (A_L, B_L)ᵀ. Because a matrix acts on what stands to its right, the leftmost interface is the rightmost factor — the product is written in reverse spatial order, and getting it backwards silently builds the mirror-image stack. Flux forces that mirror image to have the same T and the same R, but a different reflection phase and a completely different interior field, so the mistake hides from the two numbers you are most likely to check. Determinants are the free diagnostic: slab matrices contribute 1, interface matrices contribute kⱼ/kⱼ₊₁, and the chain telescopes to det M = k₀/k_N = k_L/k_R. Code that returns anything else has an interface ratio upside down, a wrong wavenumber, or a lost phase, and it should be caught there rather than in a transmission curve. A reversed product is the one bug this test cannot see, since det(AB) = det(BA).

The probability current in a lead is j = (ħkⱼ/m)(|Aⱼ|² − |Bⱼ|²) = (ħkⱼ/m) v†σ₃v with v = (Aⱼ, Bⱼ) and σ₃ = diag(1, −1). Equating the current in the two leads for every incoming state gives M†σ₃M = (k_L/k_R)σ₃: M preserves an indefinite form, not a norm. Taking determinants of that identity recovers |det M| = k_L/k_R; entry by entry it gives |M₁₁|² − |M₂₁|² = k_L/k_R, |M₂₂|² − |M₁₂|² = k_L/k_R, and M₁₁*M₁₂ = M₂₁*M₂₂. Because the metric is σ₃ and not the identity, the flux-rescaled matrix √(k_R/k_L) M sits in SU(1,1), a group of hyperbolic boosts whose entries grow without bound. Time reversal is a separate symmetry: for real V and propagating leads ψ* also solves the equation, which swaps (A, B) into (B*, A*) and forces σ₁Mσ₁ = M*, so M₂₂ = M₁₁* and M₂₁ = M₁₂*. Then M = [[a, b], [b*, a*]] with |a|² − |b|² = k_L/k_R — unimodular only when the leads match. Give V an absorptive imaginary part and the conjugate pairing goes with it.

Scattering fixes the boundary condition: unit amplitude in from the left, nothing in from the right, so v_L = (1, r) and v_R = (t, 0). The second row of (t, 0)ᵀ = M(1, r)ᵀ gives 0 = M₂₁ + M₂₂ r, hence r = −M₂₁/M₂₂; the first row, tidied with the determinant, gives t = det M/M₂₂ = (k_L/k_R)/M₂₂. Probabilities are current ratios: R = |r|² and T = (k_R/k_L)|t|² = (k_L/k_R)/|M₂₂|², and R + T = 1 then follows from the pseudo-unitarity identity instead of being imposed by hand. Run it on a rectangular barrier with matched leads and the three factors multiply out to M₂₂ = cosh(κL) − (i/2)(k/κ − κ/k) sinh(κL). Squaring, |M₂₂|² = 1 + sinh²(κL)(k² + κ²)²/(4k²κ²) = 1 + V₀² sinh²(κL)/(4E(V₀ − E)), so T = 1/|M₂₂|² is the textbook barrier formula — reached by multiplying three matrices rather than solving four matching equations.

Each forbidden slab contributes P = diag(e^(−κd), e^(+κd)), singular values e^(∓κd) and condition number e^(2κd). Those numbers multiply, so a total forbidden path L costs e^(2κL) whether you cross it in one slab or five hundred: slicing finer never helps, because the exponents merely add. Double precision carries ε = 2.22 × 10⁻¹⁶, so significance is exhausted at 2κL = ln(1/ε) = 36.0, which is 15.65 decimal digits; for an electron 4.00 eV below a barrier top, κ = 10.25 nm⁻¹ and that width is only 1.76 nm. What survives is whatever lives in the growing direction: M₂₂ is dominated by e^(+κL) and keeps full relative accuracy, so T = 1/|M₂₂|² is reliable even at 10⁻¹⁶. What dies is whatever needs the decaying direction — the field inside the barrier, the resonance peaks of a double-barrier stack where terms of size e^(2κL) must cancel down to order 1, and the root M₂₂ = 0 that locates a bound state. The cure is to compose scattering matrices: in flux-normalised amplitudes S is unitary, every entry has modulus at most 1, and the star product t = t₂(1 − r₁′r₂)⁻¹t₁ holds no growing exponential.