Notebooks · posts/transfer-matrices-and-resonant-tunneling

Tunneling through one and two barriers by transfer matrices

Adam Zettel · September 25, 2026 · written with an agent

An electron meets a barrier it does not have the energy to climb. Chop the potential into flat pieces, and the wave on each side of every edge is related by a 2 by 2 matrix; the product of those matrices gives the chance of getting through, $T$ , and of bouncing back, $R$ . They add to 1 at every energy, to within $1 0^{- 13}$ . Two barriers 4 apart, each of which alone lets through less than 1% at $E = 0.2108$ , together let through 100% at that energy.

(A_{out} B_{out}) = M (1 r), t = \frac{det M}{M _{22}}, T + R = 1

Units: $ℏ = m = 1$ , so a wave of energy $E$ in a region where the potential is $V$ has wavenumber $k = 2 (E - V)$ . Every barrier here is $V_{0} = 3$ high and 1 wide.

The reference, posts/refs/transfer-matrices-and-resonant-tunneling.py, does not multiply transfer matrices. It writes every matching condition of the whole potential as one linear system and solves it with numpy.linalg.solve, and it finds the resonances with scipy.optimize. The single barrier is also checked against its textbook closed form.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

A wave meets a step

Where the potential is flat, the time-independent Schrödinger equation $- \frac{1}{2} ψ^{''} + V ψ = E ψ$ has two solutions, a wave moving right and one moving left:

ψ (x) = A e^{ik (x - x_{0})} + B e^{- ik (x - x_{0})}, k = 2 (E - V) .

When $E < V$ the square root is imaginary, $k = iκ$ , and the two waves become a decaying and a growing exponential. That is tunneling: the wave does not stop at the wall, it fades inside it. The program takes the complex square root and never needs to ask which case it is in.

At an edge where $V$ jumps, $ψ$ and $ψ^{'}$ must be continuous. Those are two linear equations between $(A, B)$ on the left and $(A, B)$ on the right, so the right pair is a 2 by 2 matrix times the left pair. Crossing a flat region of width $d$ multiplies $A$ by $e^{ik d}$ and $B$ by $e^{- ik d}$ , another 2 by 2 matrix. The whole potential is the product of these matrices, one matrix $M$ .

Far to the left the wave is the incoming one plus a reflected one, $(A, B) = (1, r)$ . Far to the right there is only the transmitted wave, $(t, 0)$ . So $M (1, r)^{T} = (t, 0)^{T}$ , and the second row gives $r = - M_{21} / M_{22}$ , the first $t = det M / M_{22}$ . With the same potential on both sides, the transmission probability is $T = ∣ t ∣^{2}$ and the reflection probability $R = ∣ r ∣^{2}$ .

[1]lines 52–127

# every barrier is V0 high
V0 = 3.0
# the single barrier's width, and each wall's width in the double barrier
WALL = 1.0
# the gap between the two walls
GAP = 4.0
# added to 2(E - V) so that k is never exactly 0 when E equals V
K_NUDGE = 1.0e-24
PI = 3.141592653589793

#| k = sqrt(2 (E - V)), imaginary under a barrier.
def wavevector(E, V) {
    return sqrt(complex(2.0 * (E - V) + K_NUDGE, 0.0))
}

#| The matrix taking (A, B) just left of an edge to (A, B) just right of it:
#| continuity of psi and psi', with rho = k_from / k_to.
def interface(kf: Complex, kt: Complex) {
    rho = kf / kt
    D = czeros(2, 2)
    D[0, 0] = 0.5 * (1.0 + rho)
    D[0, 1] = 0.5 * (1.0 - rho)
    D[1, 0] = 0.5 * (1.0 - rho)
    D[1, 1] = 0.5 * (1.0 + rho)
    return D
}

#| Crossing a flat region of width d: A gains the phase exp(ikd), B loses it.
def propagate(k: Complex, d) {
    P = czeros(2, 2)
    P[0, 0] = exp(complex(0.0, 1.0) * k * d)
    P[1, 1] = exp(complex(0.0, -1.0) * k * d)
    return P
}

#| The product of two complex 2 by 2 matrices.
def mul2(X: CArr, Y: CArr) {
    Z = czeros(2, 2)
    for i to 2 {
        for j to 2 {
            Z[i, j] = X[i, 0] * Y[0, j] + X[i, 1] * Y[1, j]
        }
    }
    return Z
}

#| The transfer matrix of flat slabs of potential vs[j] and width ws[j],
#| left to right, with V = 0 on both sides.
def transfer(E, vs: Arr, ws: Arr) {
    k_lead = wavevector(E, 0.0)
    k_prev = k_lead
    M = czeros(2, 2)
    M[0, 0] = complex(1.0, 0.0)
    M[1, 1] = complex(1.0, 0.0)
    for j to len(vs) {
        k = wavevector(E, vs[j])
        M = mul2(interface(k_prev, k), M)
        M = mul2(propagate(k, ws[j]), M)
        k_prev = k
    }
    return mul2(interface(k_prev, k_lead), M)
}

#| T = |det M / M22|^2.
def transmission(E, vs: Arr, ws: Arr) {
    M = transfer(E, vs, ws)
    t = (M[0, 0] * M[1, 1] - M[0, 1] * M[1, 0]) / M[1, 1]
    return abs(t) * abs(t)
}

#| R = |M21 / M22|^2.
def reflection(E, vs: Arr, ws: Arr) {
    M = transfer(E, vs, ws)
    r = M[1, 0] / M[1, 1]
    return abs(r) * abs(r)
}

# every barrier is V0 high
V0 = 3.0
# the single barrier's width, and each wall's width in the double barrier
WALL = 1.0
# the gap between the two walls
GAP = 4.0
# added to 2(E - V) so that k is never exactly 0 when E equals V
K_NUDGE = 1.0e-24
PI = 3.141592653589793

#| k = sqrt(2 (E - V)), imaginary under a barrier.
def wavevector(E, V) {
    return sqrt(complex(2.0 * (E - V) + K_NUDGE, 0.0))
}

#| The matrix taking (A, B) just left of an edge to (A, B) just right of it:
#| continuity of psi and psi', with rho = k_from / k_to.
def interface(kf: Complex, kt: Complex) {
    rho = kf / kt
    D = czeros(2, 2)
    D[0, 0] = 0.5 * (1.0 + rho)
    D[0, 1] = 0.5 * (1.0 - rho)
    D[1, 0] = 0.5 * (1.0 - rho)
    D[1, 1] = 0.5 * (1.0 + rho)
    return D
}

#| Crossing a flat region of width d: A gains the phase exp(ikd), B loses it.
def propagate(k: Complex, d) {
    P = czeros(2, 2)
    P[0, 0] = exp(complex(0.0, 1.0) * k * d)
    P[1, 1] = exp(complex(0.0, -1.0) * k * d)
    return P
}

#| The product of two complex 2 by 2 matrices.
def mul2(X: CArr, Y: CArr) {
    Z = czeros(2, 2)
    for i to 2 {
        for j to 2 {
            Z[i, j] = X[i, 0] * Y[0, j] + X[i, 1] * Y[1, j]
        }
    }
    return Z
}

#| The transfer matrix of flat slabs of potential vs[j] and width ws[j],
#| left to right, with V = 0 on both sides.
def transfer(E, vs: Arr, ws: Arr) {
    k_lead = wavevector(E, 0.0)
    k_prev = k_lead
    M = czeros(2, 2)
    M[0, 0] = complex(1.0, 0.0)
    M[1, 1] = complex(1.0, 0.0)
    for j to len(vs) {
        k = wavevector(E, vs[j])
        M = mul2(interface(k_prev, k), M)
        M = mul2(propagate(k, ws[j]), M)
        k_prev = k
    }
    return mul2(interface(k_prev, k_lead), M)
}

#| T = |det M / M22|^2.
def transmission(E, vs: Arr, ws: Arr) {
    M = transfer(E, vs, ws)
    t = (M[0, 0] * M[1, 1] - M[0, 1] * M[1, 0]) / M[1, 1]
    return abs(t) * abs(t)
}

#| R = |M21 / M22|^2.
def reflection(E, vs: Arr, ws: Arr) {
    M = transfer(E, vs, ws)
    r = M[1, 0] / M[1, 1]
    return abs(r) * abs(r)
}

One barrier, and the check that nothing is lost

For a single rectangular barrier the answer is known in closed form. Below the top, with $κ = 2 (V_{0} - E)$ and width $a$ ,

T = [1 + \frac{V _{0}^{2} sinh ^{2} ( κa )}{4 E ( V _{0} - E )}]^{- 1},

and above it $sinh (κa)$ becomes $sin (q a)$ with $q = 2 (E - V_{0})$ and $V_{0} - E$ becomes $E - V_{0}$ . The cell compares the matrix product with this formula at 500 energies from 0.05 to 10.

The second check needs no formula. The electron is either transmitted or reflected, so $T + R = 1$ . Nothing in the matrix product enforces that; it comes out only if every matching condition is right. Before running: above the barrier top, is $T$ equal to 1? Classically yes. For a wave it is 1 only where a whole number of half wavelengths fit in the barrier, $q a = nπ$ , and below 1 between those energies.

[2]lines 149–211

ONE_V = [V0]
ONE_W = [WALL]
TWO_V = [V0, 0.0, V0]
TWO_W = [WALL, GAP, WALL]
N_CHECK = 500
# the over-barrier resonances checked: n = 1, 2, 3
N_OVER = 3
E_CHECK = linspace(0.05, 10.0, N_CHECK)

#| The closed-form transmission through one barrier of height V0 and width WALL.
def barrier_exact(E) {
    if E < V0 {
        kap = sqrt(2.0 * (V0 - E))
        return 1.0 / (1.0 + V0 * V0 * sinh(kap * WALL) * sinh(kap * WALL) / (4.0 * E * (V0 - E)))
    }
    q = sqrt(2.0 * (E - V0))
    return 1.0 / (1.0 + V0 * V0 * sin(q * WALL) * sin(q * WALL) / (4.0 * E * (E - V0)))
}

T_one = zeros(N_CHECK)
T_formula = zeros(N_CHECK)
err_formula = 0.0
err_flux_one = 0.0
err_flux_two = 0.0
for i to N_CHECK {
    E = E_CHECK[i]
    T_one[i] = transmission(E, ONE_V, ONE_W)
    T_formula[i] = barrier_exact(E)
    err_formula = max(err_formula, abs(T_one[i] - T_formula[i]))
    err_flux_one = max(err_flux_one, abs(T_one[i] + reflection(E, ONE_V, ONE_W) - 1.0))
    err_flux_two = max(err_flux_two, abs(transmission(E, TWO_V, TWO_W) + reflection(E, TWO_V, TWO_W) - 1.0))
}
show err_formula
show err_flux_one
show err_flux_two

test "the matrix product matches the closed form at 500 energies to 1e-13" {
    assert_lt(err_formula, 1.0e-13)
}
test "T + R = 1 to 1e-13 at every energy, for one barrier and for two" {
    assert_lt(err_flux_one, 1.0e-13)
    assert_lt(err_flux_two, 1.0e-13)
}
test "one barrier lets through 6.33% at E = 1, the closed form's 0.0632976" {
    # the closed form at E = 1, in the reference output
    assert_near(transmission(1.0, ONE_V, ONE_W), 0.063297625363163315, 1.0e-14)
}
test "above the top T = 1 exactly where q a = n pi: E = 7.9348, 22.739, 47.413" {
    for n from 1 to N_OVER + 1 {
        E_n = V0 + (n * PI) * (n * PI) / (2.0 * WALL * WALL)
        assert_near(transmission(E_n, ONE_V, ONE_W), 1.0, 1.0e-13)
    }
}

first_top = V0 + PI * PI / (2.0 * WALL * WALL)
figure("One barrier, V0 = 3 and width 1", "energy E", "transmission T", 0.0, 10.0, 0.0, 1.05)
plot(E_CHECK, T_one, "transfer matrices")
plot(E_CHECK, T_formula, "closed form")
style("closed form", "dash s2")
vline("top", V0)
label("toplabel", "barrier top", V0 + 0.9, 0.2)
vline("res1", first_top)
label("reslabel", "$qa = \\pi$", first_top + 0.7, 0.2)

ONE_V = [V0]
ONE_W = [WALL]
TWO_V = [V0, 0.0, V0]
TWO_W = [WALL, GAP, WALL]
N_CHECK = 500
# the over-barrier resonances checked: n = 1, 2, 3
N_OVER = 3
E_CHECK = linspace(0.05, 10.0, N_CHECK)

#| The closed-form transmission through one barrier of height V0 and width WALL.
def barrier_exact(E) {
    if E < V0 {
        kap = sqrt(2.0 * (V0 - E))
        return 1.0 / (1.0 + V0 * V0 * sinh(kap * WALL) * sinh(kap * WALL) / (4.0 * E * (V0 - E)))
    }
    q = sqrt(2.0 * (E - V0))
    return 1.0 / (1.0 + V0 * V0 * sin(q * WALL) * sin(q * WALL) / (4.0 * E * (E - V0)))
}

T_one = zeros(N_CHECK)
T_formula = zeros(N_CHECK)
err_formula = 0.0
err_flux_one = 0.0
err_flux_two = 0.0
for i to N_CHECK {
    E = E_CHECK[i]
    T_one[i] = transmission(E, ONE_V, ONE_W)
    T_formula[i] = barrier_exact(E)
    err_formula = max(err_formula, abs(T_one[i] - T_formula[i]))
    err_flux_one = max(err_flux_one, abs(T_one[i] + reflection(E, ONE_V, ONE_W) - 1.0))
    err_flux_two = max(err_flux_two, abs(transmission(E, TWO_V, TWO_W) + reflection(E, TWO_V, TWO_W) - 1.0))
}
show err_formula
show err_flux_one
show err_flux_two

test "the matrix product matches the closed form at 500 energies to 1e-13" {
    assert_lt(err_formula, 1.0e-13)
}
test "T + R = 1 to 1e-13 at every energy, for one barrier and for two" {
    assert_lt(err_flux_one, 1.0e-13)
    assert_lt(err_flux_two, 1.0e-13)
}
test "one barrier lets through 6.33% at E = 1, the closed form's 0.0632976" {
    # the closed form at E = 1, in the reference output
    assert_near(transmission(1.0, ONE_V, ONE_W), 0.063297625363163315, 1.0e-14)
}
test "above the top T = 1 exactly where q a = n pi: E = 7.9348, 22.739, 47.413" {
    for n from 1 to N_OVER + 1 {
        E_n = V0 + (n * PI) * (n * PI) / (2.0 * WALL * WALL)
        assert_near(transmission(E_n, ONE_V, ONE_W), 1.0, 1.0e-13)
    }
}

first_top = V0 + PI * PI / (2.0 * WALL * WALL)
figure("One barrier, V0 = 3 and width 1", "energy E", "transmission T", 0.0, 10.0, 0.0, 1.05)
plot(E_CHECK, T_one, "transfer matrices")
plot(E_CHECK, T_formula, "closed form")
style("closed form", "dash s2")
vline("top", V0)
label("toplabel", "barrier top", V0 + 0.9, 0.2)
vline("res1", first_top)
label("reslabel", "$qa = \\pi$", first_top + 0.7, 0.2)

Two barriers

Now put a second barrier 4 to the right of the first. Between them is a well, and a particle trapped in a well of width 4 with infinitely high walls could have only the energies $(nπ /4)^{2} /2$ : 0.308, 1.234, 2.776. With walls of finite height the particle can leak out, so these are not exact states, but near such an energy a wave that gets into the well reflects back and forth between the walls and adds up in phase.

Each wall alone lets through about 1% at the lowest of these energies. A naive guess for two walls in a row is the product, $0.009 4^{2} = 8.8 \times 1 0^{- 5}$ . Before running: what is $T$ at the resonance?

The cell scans 3000 energies below the top, takes every local maximum of $T$ , and refines each by golden-section search, a bracket that shrinks by the golden ratio each step, until it is $1 0^{- 12}$ wide. Then it finds the two energies around each peak where $T = 1/2$ by bisection; their distance is the width of the resonance.

[3]lines 234–353

N_SCAN = 3000
E_SCAN = linspace(0.005, V0 - 0.005, N_SCAN)
# golden-section and bisection stop when the bracket is this wide
E_TOL = 1.0e-12
GOLD = 0.6180339887498949

#| The energy of largest T in [lo, hi], by golden-section search.
def golden_peak(lo, hi) {
    a = lo
    b = hi
    c = b - GOLD * (b - a)
    d = a + GOLD * (b - a)
    fc = transmission(c, TWO_V, TWO_W)
    fd = transmission(d, TWO_V, TWO_W)
    while b - a > E_TOL {
        if fc > fd {
            b = d
            d = c
            fd = fc
            c = b - GOLD * (b - a)
            fc = transmission(c, TWO_V, TWO_W)
        } else {
            a = c
            c = d
            fc = fd
            d = a + GOLD * (b - a)
            fd = transmission(d, TWO_V, TWO_W)
        }
    }
    return 0.5 * (a + b)
}

#| Where T crosses 1/2 between e_in (T above) and e_out (T below), by bisection.
def half_point(e_in, e_out) {
    a = e_in
    b = e_out
    while abs(b - a) > E_TOL {
        m = 0.5 * (a + b)
        if transmission(m, TWO_V, TWO_W) > 0.5 {
            a = m
        } else {
            b = m
        }
    }
    return 0.5 * (a + b)
}

#| The distance between the two energies around the peak at e_res where T = 1/2:
#| step outward, doubling the step, until T is below 1/2, then bisect.
def half_width(e_res) {
    step = 1.0e-6
    while transmission(e_res - step, TWO_V, TWO_W) > 0.5 {
        step = 2.0 * step
    }
    lo = half_point(e_res, e_res - step)
    step = 1.0e-6
    while transmission(e_res + step, TWO_V, TWO_W) > 0.5 {
        step = 2.0 * step
    }
    hi = half_point(e_res, e_res + step)
    return hi - lo
}

T_scan = zeros(N_SCAN)
for i to N_SCAN {
    T_scan[i] = transmission(E_SCAN[i], TWO_V, TWO_W)
}
n_res = 0
for i from 1 to N_SCAN - 1 {
    if T_scan[i] > T_scan[i - 1] and T_scan[i] >= T_scan[i + 1] {
        n_res = n_res + 1
    }
}
E_res = zeros(n_res)
T_res = zeros(n_res)
W_res = zeros(n_res)
T_alone = zeros(n_res)
kk = 0
for i from 1 to N_SCAN - 1 {
    if T_scan[i] > T_scan[i - 1] and T_scan[i] >= T_scan[i + 1] {
        E_res[kk] = golden_peak(E_SCAN[i - 1], E_SCAN[i + 1])
        T_res[kk] = transmission(E_res[kk], TWO_V, TWO_W)
        W_res[kk] = half_width(E_res[kk])
        T_alone[kk] = transmission(E_res[kk], ONE_V, ONE_W)
        kk = kk + 1
    }
}
show n_res
show E_res
show T_res
show W_res
show T_alone
T_between = transmission(0.5, TWO_V, TWO_W)
show T_between

test "three resonances below the top, at the reference's 0.21079395, 0.83247284 and 1.82964090" {
    assert_near(n_res, 3.0, 0.0)
    # reference: the zero of dT/dE for the direct linear solve (posts/refs/...out).
    # T is flat at its top, T = 1 - 4 (E - E_res)^2 / W^2 roughly, so a T known
    # to 1e-15 places the peak only to about W * 1e-7.5: 1e-9 at the widest.
    assert_near(E_res[0], 0.21079395104633, 1.0e-9)
    assert_near(E_res[1], 0.83247283693652, 1.0e-9)
    assert_near(E_res[2], 1.82964089984709, 1.0e-9)
}
test "at each resonance T = 1 to 1e-12, though one wall alone passes 0.936%, 4.89% and 16.4%" {
    for j to n_res {
        assert_near(T_res[j], 1.0, 1.0e-12)
    }
    assert_near(T_alone[0], 9.362446e-3, 1.0e-9)
    assert_near(T_alone[1], 4.892834e-2, 1.0e-8)
    assert_near(T_alone[2], 1.642376e-1, 1.0e-7)
}
test "the widths at half height are 1.2572e-3, 1.3063e-2 and 6.6907e-2, as the reference's root finder gives" {
    assert_near(W_res[0], 1.257211710863e-3, 1.0e-10)
    assert_near(W_res[1], 1.306259538942e-2, 1.0e-10)
    assert_near(W_res[2], 6.690700145784e-2, 1.0e-10)
}
test "between the first two resonances, at E = 0.5, two walls pass 1.676e-4" {
    assert_near(T_between, 1.676323351521634e-4, 1.0e-15)
}

N_SCAN = 3000
E_SCAN = linspace(0.005, V0 - 0.005, N_SCAN)
# golden-section and bisection stop when the bracket is this wide
E_TOL = 1.0e-12
GOLD = 0.6180339887498949

#| The energy of largest T in [lo, hi], by golden-section search.
def golden_peak(lo, hi) {
    a = lo
    b = hi
    c = b - GOLD * (b - a)
    d = a + GOLD * (b - a)
    fc = transmission(c, TWO_V, TWO_W)
    fd = transmission(d, TWO_V, TWO_W)
    while b - a > E_TOL {
        if fc > fd {
            b = d
            d = c
            fd = fc
            c = b - GOLD * (b - a)
            fc = transmission(c, TWO_V, TWO_W)
        } else {
            a = c
            c = d
            fc = fd
            d = a + GOLD * (b - a)
            fd = transmission(d, TWO_V, TWO_W)
        }
    }
    return 0.5 * (a + b)
}

#| Where T crosses 1/2 between e_in (T above) and e_out (T below), by bisection.
def half_point(e_in, e_out) {
    a = e_in
    b = e_out
    while abs(b - a) > E_TOL {
        m = 0.5 * (a + b)
        if transmission(m, TWO_V, TWO_W) > 0.5 {
            a = m
        } else {
            b = m
        }
    }
    return 0.5 * (a + b)
}

#| The distance between the two energies around the peak at e_res where T = 1/2:
#| step outward, doubling the step, until T is below 1/2, then bisect.
def half_width(e_res) {
    step = 1.0e-6
    while transmission(e_res - step, TWO_V, TWO_W) > 0.5 {
        step = 2.0 * step
    }
    lo = half_point(e_res, e_res - step)
    step = 1.0e-6
    while transmission(e_res + step, TWO_V, TWO_W) > 0.5 {
        step = 2.0 * step
    }
    hi = half_point(e_res, e_res + step)
    return hi - lo
}

T_scan = zeros(N_SCAN)
for i to N_SCAN {
    T_scan[i] = transmission(E_SCAN[i], TWO_V, TWO_W)
}
n_res = 0
for i from 1 to N_SCAN - 1 {
    if T_scan[i] > T_scan[i - 1] and T_scan[i] >= T_scan[i + 1] {
        n_res = n_res + 1
    }
}
E_res = zeros(n_res)
T_res = zeros(n_res)
W_res = zeros(n_res)
T_alone = zeros(n_res)
kk = 0
for i from 1 to N_SCAN - 1 {
    if T_scan[i] > T_scan[i - 1] and T_scan[i] >= T_scan[i + 1] {
        E_res[kk] = golden_peak(E_SCAN[i - 1], E_SCAN[i + 1])
        T_res[kk] = transmission(E_res[kk], TWO_V, TWO_W)
        W_res[kk] = half_width(E_res[kk])
        T_alone[kk] = transmission(E_res[kk], ONE_V, ONE_W)
        kk = kk + 1
    }
}
show n_res
show E_res
show T_res
show W_res
show T_alone
T_between = transmission(0.5, TWO_V, TWO_W)
show T_between

test "three resonances below the top, at the reference's 0.21079395, 0.83247284 and 1.82964090" {
    assert_near(n_res, 3.0, 0.0)
    # reference: the zero of dT/dE for the direct linear solve (posts/refs/...out).
    # T is flat at its top, T = 1 - 4 (E - E_res)^2 / W^2 roughly, so a T known
    # to 1e-15 places the peak only to about W * 1e-7.5: 1e-9 at the widest.
    assert_near(E_res[0], 0.21079395104633, 1.0e-9)
    assert_near(E_res[1], 0.83247283693652, 1.0e-9)
    assert_near(E_res[2], 1.82964089984709, 1.0e-9)
}
test "at each resonance T = 1 to 1e-12, though one wall alone passes 0.936%, 4.89% and 16.4%" {
    for j to n_res {
        assert_near(T_res[j], 1.0, 1.0e-12)
    }
    assert_near(T_alone[0], 9.362446e-3, 1.0e-9)
    assert_near(T_alone[1], 4.892834e-2, 1.0e-8)
    assert_near(T_alone[2], 1.642376e-1, 1.0e-7)
}
test "the widths at half height are 1.2572e-3, 1.3063e-2 and 6.6907e-2, as the reference's root finder gives" {
    assert_near(W_res[0], 1.257211710863e-3, 1.0e-10)
    assert_near(W_res[1], 1.306259538942e-2, 1.0e-10)
    assert_near(W_res[2], 6.690700145784e-2, 1.0e-10)
}
test "between the first two resonances, at E = 0.5, two walls pass 1.676e-4" {
    assert_near(T_between, 1.676323351521634e-4, 1.0e-15)
}

Every resonance reaches $T = 1$ to within $1 0^{- 12}$ (T_res), while a single wall at the same energy passes 0.94%, 4.9% and 16% (T_alone). The two walls are identical, so the waves reflected by the far wall and leaking back through the near one cancel the near wall's own reflection exactly at resonance. The resonances sit below the infinite-well energies (0.211 against 0.308), because the wave reaches into the walls and the well is effectively wider. The widths grow with the single-wall transmission: a leakier well holds the wave for a shorter time, and a shorter-lived state has a broader resonance.

The plot energies below are 400 evenly spaced points plus 41 points across each peak, since the first peak is narrower than the spacing of the even ones. The single wall is the dashed curve.

[4]lines 370–402

N_PLOT = 400
N_PEAK = 41
# the dense points span this many half-widths on each side of a peak
PEAK_SPAN = 4.0
E_plot = zeros(N_PLOT + n_res * N_PEAK)
E_even = linspace(0.005, 3.5, N_PLOT)
for i to N_PLOT {
    E_plot[i] = E_even[i]
}
for j to n_res {
    E_near = linspace(E_res[j] - PEAK_SPAN * W_res[j], E_res[j] + PEAK_SPAN * W_res[j], N_PEAK)
    for i to N_PEAK {
        E_plot[N_PLOT + j * N_PEAK + i] = E_near[i]
    }
}
sort(E_plot)
T_two_plot = zeros(len(E_plot))
T_one_plot = zeros(len(E_plot))
for i to len(E_plot) {
    T_two_plot[i] = transmission(E_plot[i], TWO_V, TWO_W)
    T_one_plot[i] = transmission(E_plot[i], ONE_V, ONE_W)
}
figure("Two walls 4 apart: T on a log scale", "energy E", "transmission T", 0.0, 3.5, 1.0e-6, 2.0)
yscale("log")
plot(E_plot, T_two_plot, "two walls")
plot(E_plot, T_one_plot, "one wall")
style("one wall", "dash s2")
scatter(E_res, T_res, "resonances found")
style("resonances found", "s3")
vline("top2", V0)

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

N_PLOT = 400
N_PEAK = 41
# the dense points span this many half-widths on each side of a peak
PEAK_SPAN = 4.0
E_plot = zeros(N_PLOT + n_res * N_PEAK)
E_even = linspace(0.005, 3.5, N_PLOT)
for i to N_PLOT {
    E_plot[i] = E_even[i]
}
for j to n_res {
    E_near = linspace(E_res[j] - PEAK_SPAN * W_res[j], E_res[j] + PEAK_SPAN * W_res[j], N_PEAK)
    for i to N_PEAK {
        E_plot[N_PLOT + j * N_PEAK + i] = E_near[i]
    }
}
sort(E_plot)
T_two_plot = zeros(len(E_plot))
T_one_plot = zeros(len(E_plot))
for i to len(E_plot) {
    T_two_plot[i] = transmission(E_plot[i], TWO_V, TWO_W)
    T_one_plot[i] = transmission(E_plot[i], ONE_V, ONE_W)
}
figure("Two walls 4 apart: T on a log scale", "energy E", "transmission T", 0.0, 3.5, 1.0e-6, 2.0)
yscale("log")
plot(E_plot, T_two_plot, "two walls")
plot(E_plot, T_one_plot, "one wall")
style("one wall", "dash s2")
scatter(E_res, T_res, "resonances found")
style("resonances found", "s3")
vline("top2", V0)

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors