University Quantum Mechanics II · Semiclassical Methods: WKB and the Path Integral · 11.4

The WKB wavefunction from an ħ expansion, and its validity

WKB is usually handed over as an inspired guess. Here it is derived one power of ħ at a time, so you know exactly what was kept and what was dropped, and you leave with a number, δ, that says at each point on the x-axis whether to trust it. Every later result in this unit leans on that number.

Open the lab ↓Practice & worked examples ↓

Build the model

Connect the measurement to the mechanism.

Write ψ = exp(iS(x)/ħ) with S complex and nothing has been approximated: the Schrödinger equation becomes (S′)² − iħS″ = p(x)², with p = √(2m(E − V)), and ħ now multiplies the only term with no classical meaning. Expanding S = S₀ + ħS₁ + ħ²S₂ + … turns that into a hierarchy. Order ħ⁰ is the time-independent Hamilton–Jacobi equation, S₀ = ±∫p dx: the phase advances one radian per reduced wavelength ħ/p, the fixed-energy form of the action whose stationary phase picked out the classical path.

Order ħ¹ fixes the amplitude, S₁ = (i/2) ln p, so ψ ≈ (C/√p) exp(±(i/ħ)∫p dx), with growing and decaying real exponentials where V exceeds E. The 1/√p is not decoration: it holds the probability current at |C|²/m, and in a well it makes the locally averaged |ψ|² equal the classical fraction of time spent near x. What was dropped is ħS″ against (S′)², the ratio δ = |d(ħ/p)/dx| = ħm|V′|/p³.

It is local, and it does not care how deep the well is: small where the wavelength is short and slowly changing, infinite at every turning point, where p = 0, whatever the quantum number. Beyond two orders the series is asymptotic rather than convergent. WKB is therefore a leading-order statement in ħ with a computable, local error, and the turning points are where the next three lessons go to work.

Simple definition: The WKB wavefunction is the solution of the Schrödinger equation correct through first order in ħ, ψ ≈ (C/√p) exp(±(i/ħ)∫p dx), trustworthy wherever the reduced wavelength ħ/p changes by a small fraction of itself over a distance ħ/p.
Example: In an oscillator's n = 10 state, δ = u/(21(1 − u²)³⁄²) at x = u xₜ: zero at the centre, 0.037 halfway out, 0.52 at 0.9 xₜ and infinite at the turning point, so one formula is excellent in the middle and useless at the edge.

Exact equation for the phase(S′)² − iħS″ = p(x)², with p(x) = √(2m(E − V(x)))

Exact, derived by substitution: nothing is dropped yet. ħ multiplies the only non-classical term, which makes it the expansion parameter, and a singular one.

ψ = e^(iS/ħ); S complex, in J s; p in kg m s⁻¹; E and V in J; primes are d/dx

Order ħ⁰: the Hamilton–Jacobi equation(S₀′)² = p(x)² ⇒ S₀(x) = ±∫ˣ p(x′) dx′

Derived. The phase gradient is the local momentum, so λ(x) = 2πħ/p(x) is the local de Broglie wavelength. Exact when V is constant.

S₀ in J s is Hamilton's characteristic function W(x, E); the lower limit only sets a constant phase

Order ħ¹: the WKB wavefunctionS₁ = (i/2) ln p ⇒ ψ ≈ (C/√p) exp(±(i/ħ)∫p dx)

Derived, correct to first order in ħ. The 1/√p holds the current of each branch at j = |C|²/m for every x, as a stationary state demands.

|C|² in kg s⁻¹, so |ψ|² is in m⁻¹; where V exceeds E, |p| = √(2m(V − E)) and the exponent is ±(1/ħ)∫|p| dx

Locally averaged density in a well⟨|ψ|²⟩ = 2|C|²/p(x) = 2/(T v(x)), with |C|² = m/T

Derived: the classical probability of being near x. The particle is likeliest where it is slowest, but 1/v diverges at the turning points, where the true |ψ|² stays finite.

⟨ ⟩ is an average over one local wavelength, not an expectation value; T the classical period in s, v = p/m in m s⁻¹

Validity parameterδ(x) = |d(ħ/p)/dx| = ħm|V′(x)|/p(x)³ ≪ 1

Derived as ħ|S₀″|/(S₀′)², the dropped term over the kept one. Local: λ̄ must change by a small fraction of itself over a distance λ̄.

dimensionless; ħ/p = λ̄ is the reduced wavelength in m, |V′| the force magnitude in N, m in kg

Validity next to a turning point x = aδ = ½(ℓ/|x − a|)³⁄², with ℓ = (ħ²/(2m|V′(a)|))¹⁄³

Derived from the validity parameter with V linear; equivalently δ = λ̄/(2|x − a|). Every smooth turning point has a failure zone a few ℓ wide.

ℓ in m, V linearised about a; δ = 0.1 at |x − a| = 2.9ℓ, and δ → ∞ as x → a

Put ħ in the phase, and nothing is approximated yet

Write each travelling branch as ψ(x) = exp(iS(x)/ħ) with S complex: its real part is a phase and its imaginary part a log-amplitude. Then ψ′ = (i/ħ)S′ψ and ψ″ = [(i/ħ)S″ − (S′)²/ħ²]ψ, and −(ħ²/2m)ψ″ + Vψ = Eψ becomes (S′)² − iħS″ = 2m(E − V) ≡ p(x)², exactly. A second-order linear equation has become a first-order nonlinear one for S′, and ħ appears once, in front of the only term with no classical meaning. That makes ħ the expansion parameter, but a singular one: it multiplies the highest derivative, so setting ħ = 0 turns a differential equation into an algebraic one, and the information lost, here the amplitude, has to come back at the next order. Small ħ means large action: the real parameter is ħ divided by ∫p dx over the region of interest, the comparison that made stationary phase work in the path integral.

Order ħ⁰: the classical Hamilton–Jacobi equation

Put S = S₀ + ħS₁ + ħ²S₂ + … into the exact equation and collect powers. At ħ⁰, (S₀′)² = p², so S₀′ = ±p and S₀ = ±∫ˣ p(x′) dx′. Divided by 2m it reads (S₀′)²/2m + V = E, the time-independent Hamilton–Jacobi equation, with S₀ Hamilton's characteristic function W(x, E), quoted from classical mechanics; W − Et is the action whose stationary phase gave the classical path in the last lesson. The phase gradient is the local momentum, so the local wavelength is λ(x) = 2πħ/p(x): an electron with 1.00 eV of kinetic energy has λ̄ = ħ/p = 0.195 nm, and one with 4.00 eV has 0.0976 nm. For constant V the order-zero phase is a plane wave, every higher term vanishes, and the result is exact. Where V exceeds E, p = i|p| with |p| = √(2m(V − E)) and S₀ = ±i∫|p| dx: the waves turn into real exponentials, growing or decaying by one e-fold per ħ/|p|.

Order ħ¹: the amplitude that keeps the current constant

The ħ¹ terms give 2S₀′S₁′ − iS₀″ = 0, so S₁′ = (i/2)S₀″/S₀′ = (i/2)p′/p for either sign, and S₁ = (i/2) ln p. Since exp(iS₁) = p^(−1/2), two orders give ψ ≈ (C/√p) exp(±(i/ħ)∫p dx), and in a forbidden region ψ ≈ (C/√|p|) exp(±(1/ħ)∫|p| dx). The amplitude is conservation at work. For the + branch ψ′ = (ip/ħ − p′/2p)ψ, so j = (ħ/m) Im(ψ*ψ′) = |C|²/m at every x. A stationary state needs dj/dx = 0, and any other power of p would pile probability up or drain it away: where the particle moves slower it must be denser, |ψ|² ∝ 1/v. A single real exponential has ψ*ψ′ real and carries no current at all; current under a barrier needs both exponentials at once, which is how tunnelling is built at the end of this unit. S₀ and S₁ together are the WKB wavefunction, correct to first order in ħ.

Averaged over a wavelength, |ψ|² is the classical dwell time

A bound state carries no net current, so its two branches have equal magnitude and combine into ψ = (2|C|/√p) cos(∫p dx/ħ + θ), with θ a phase the connection formulas will fix. Its density, (4|C|²/p)cos²(…), oscillates, and averaging over one local wavelength replaces cos² by ½, leaving 2|C|²/p = 2|C|²/(mv). A classical particle of period T crosses dx twice per period, taking dx/v each time, so its probability of being in dx is 2 dx/(Tv). The two agree when |C|² = m/T, which is also what normalises the WKB density. For the oscillator the dwell density is 1/(π√(xₜ² − x²)). A check, computed numerically: for n = 30 the probability of |x| ≤ xₜ/2 is 0.3400 from the exact Hermite function, 0.3401 from the WKB standing wave at the exact energy, and 1/3 from the dwell density, which averages over fringes that the cut at xₜ/2 slices through. At the turning points, though, 1/v diverges while the true |ψ|² stays finite.

The validity parameter, measured point by point

The term dropped at order zero was ħS″ against (S′)². On S₀′ = p their ratio is δ = ħ|p′|/p² = |d(ħ/p)/dx|, and since pp′ = −mV′, δ = ħm|V′|/p³: dimensionless, since J s × kg × N and (kg m s⁻¹)³ are both kg³ m³ s⁻³. It is the fractional change of λ̄ = ħ/p across a distance λ̄; equivalently the potential may change by only a small part of 2(E − V) over one λ̄. Nothing asks V to be small, and a deep well is fine. For the oscillator's n-th state at x = u xₜ, δ = u/((2n + 1)(1 − u²)³⁄²). The 1/(2n + 1) is ħ in disguise, since ħ/∮p dx = 1/(2π(n + ½)), but the other factor diverges at u = 1. It tracks the exact solution, computed numerically: for n = 10 the WKB fringe heights match the Hermite function's to 0.6% out to the fourth fringe from the centre (δ ≤ 0.07), while the fifth and last, at δ ≈ 0.46, is 12% too high. Near any turning point, δ = ½(ℓ/|x − a|)³⁄² = λ̄/(2|x − a|): you need to be many reduced wavelengths from it, and λ̄ grows as you approach.

The next order, and why the series is only asymptotic

The ħ² terms give 2S₀′S₂′ + (S₁′)² − iS₁″ = 0, so on the + branch S₂′ = 3p′²/(8p³) − p″/(4p²). It is real where E exceeds V, so the second order corrects the phase, not the amplitude, and ħS₂ is of order δ. For a linear ramp it integrates in closed form to ħS₂ = (5/24)δ, which at 3ℓ from the turning point is 0.020 rad against a leading phase of 3.46 rad. In the oscillator at n = 30 the exact fringe peaks sit about δ/4 radian from the WKB ones, as ħS₂ predicts, and the outer peak heights are off by about δ², 0.6% at δ = 0.08. So keeping S₀ and S₁ leaves a local error controlled by δ, the precise sense in which WKB is leading order in ħ. More orders do not converge on the exact answer: at fixed ħ the coefficients eventually grow factorially, so the series is asymptotic (quoted). Each extra term helps only while δ is small, and none repairs a turning point, where every term is more singular than the last. Keep two orders, check δ, and hand the turning points to the linearised potential of the next lesson.

Change one variable at a time

Make the relationship visible.

Interactive model

Quantum number n10

Probe position u = x/xₜ0.80

Leave the probe at u = 0.80 and raise n from 0 to 30: δ there falls as 1/(2n + 1), from 3.70 to 0.061, and the fringes crowd until their average is the dashed curve. The δ curve still spikes at both turning points for every n, and at n = 30 the shaded zones hold 35% of the classical dwell time.

δ AT THE PROBE0.176

ħ / ∮p dx0.0152

δ = 0.1 REACHED AT |x|/xₜ0.716

DWELL TIME WITH δ ≥ 0.149.2 %

Live interpretationδ AT THE PROBE: 0.176. ħ / ∮p dx: 0.0152. δ = 0.1 REACHED AT |x|/xₜ: 0.716. DWELL TIME WITH δ ≥ 0.1: 49.2 %

Catch the common trap

Explain before calculating.

Misconception check

“WKB is the large-n approximation, so a highly excited state is described well everywhere.” The error is uniformity. At a fixed fraction u of the way to the turning point, δ does fall as 1/(2n + 1), but at u = 1 it is infinite for every n, and the zone where δ exceeds 0.1 shrinks slowly: in the oscillator δ ≤ 0.1 reaches |x| = 0.10xₜ at n = 0, 0.85xₜ at n = 30 and 0.93xₜ at n = 100. The classical particle lingers exactly there. At n = 30 it spends 35% of its time where δ exceeds 0.1, 23% at n = 100 and 11% at n = 1000, a fraction falling only as n^(−1/3). Anything that depends on the neighbourhood of a turning point, the wavefunction at the edge, the phase lost on reflection, a tunnelling amplitude, needs the turning-point analysis at every n. Its twin is trusting WKB wherever E exceeds V: δ, not the sign of E − V, is the test, and WKB works well deep inside a forbidden region too.

A particle in a harmonic well is in the state n = 200. Is the WKB wavefunction reliable at x = 0.99xₜ, where xₜ is the classical turning point?

Choose an answer to test the model.

Practice & worked examples

Reason from the model, then test the result.

EasyCaesium fountain clocks toss atoms straight up at about 4 m s⁻¹. Take a caesium-133 atom, m = 2.207 × 10⁻²⁵ kg, launched at 4.0 m s⁻¹. Evaluate the WKB validity parameter at launch, then find the speed and the depth below the apex at which δ reaches 0.1.

Uniform gravity: V = mgz, so |V′| = mg and p = mv. Then δ = ħm|V′|/p³ = ħg/(mv³), which depends on the height only through v.
At launch: δ = (1.055 × 10⁻³⁴ J s)(9.81 m s⁻²)/[(2.207 × 10⁻²⁵ kg)(4.0 m s⁻¹)³] = 7.3 × 10⁻¹¹. WKB is about as good as an approximation gets.
δ = 0.1 when v³ = 10ħg/m = 4.69 × 10⁻⁸ m³ s⁻³, so v = 3.6 × 10⁻³ m s⁻¹.
That speed is reached h = v²/2g = 6.6 × 10⁻⁷ m below the apex. Check with the turning-point form: ℓ = (ħ²/(2m²g))¹⁄³ = 0.227 μm, and 5²⁄³ℓ = 2.92ℓ = 0.66 μm, the same.
The apex is (4.0 m s⁻¹)²/(2 × 9.81 m s⁻²) = 0.82 m up, so WKB fails only in the top 0.66 μm, less than one part in 10⁶ of the climb, and for 0.09% of the flight time.

Answerδ = 7.3 × 10⁻¹¹ at launch; δ reaches 0.1 at v = 3.6 mm s⁻¹, 0.66 μm below an apex 0.82 m up.

MediumA ⁴⁰Ca⁺ ion, m = 6.636 × 10⁻²⁶ kg, in a trap of frequency ω/2π = 1.00 MHz is in its n = 30 motional state, with Quantum Mechanics I's exact energy (n + ½)ħω. Find the turning point xₜ, δ at 0.5xₜ and 0.9xₜ, where δ = 0.1, and the fraction of the classical dwell time spent beyond that point.

x₀ = √(ħ/mω) = √(1.055 × 10⁻³⁴ J s/(6.636 × 10⁻²⁶ kg × 6.283 × 10⁶ s⁻¹)) = 15.9 nm, and ½mω²xₜ² = (n + ½)ħω gives xₜ = √61 x₀ = 124 nm.
With V′ = mω²x and p = mω√(xₜ² − x²), δ = ħm|V′|/p³ = u/((2n + 1)(1 − u²)³⁄²) at x = u xₜ. So δ(0.5) = 0.5/(61 × 0.650) = 0.0126 and δ(0.9) = 0.9/(61 × 0.0828) = 0.178.
δ = 0.1 means u = 6.1(1 − u²)³⁄². Put w = 1 − u² and square: 37.21w³ + w − 1 = 0, whose one positive root is w = 0.2697, so u = √0.7303 = 0.8546 and |x| = 106 nm.
Classically the fraction of time with |x| ≤ u xₜ is (2/π) arcsin u = (2/π)(1.0248) = 0.652, so the ion spends 0.348 of its time beyond |x| = 106 nm, in the zone where δ exceeds 0.1.

Answerxₜ = 124 nm; δ = 0.013 at xₜ/2 and 0.18 at 0.9xₜ; δ = 0.1 at |x| = 106 nm, or 0.855xₜ, and 35% of the classical dwell time lies beyond it.

HardAn electron decelerated by a uniform electric field of magnitude 1.00 × 10⁶ V m⁻¹ turns back at x = a. Derive the order-ħ² phase ħS₂ for a linear potential, evaluate it and δ at d = 10.0 nm before the turning point, and find the d at which ħS₂ would equal the leading phase measured from the turning point.

The force is F = e × 1.00 × 10⁶ V m⁻¹ = 1.602 × 10⁻¹³ N, so ℓ = (ħ²/(2mₑF))¹⁄³ = (1.112 × 10⁻⁶⁸ J² s²/2.919 × 10⁻⁴³ kg N)¹⁄³ = 3.36 nm.
With d = a − x, p² = 2mFd, so p′ = −mF/p and p″ = −m²F²/p³. Then S₂′ = 3p′²/(8p³) − p″/(4p²) = (3/8 + 1/4)m²F²/p⁵ = (5/8)m²F²/p⁵.
Integrate, choosing the constant so the correction vanishes far away: S₂ = (5/12)m²F²/((2mF)⁵⁄²d³⁄²). With ℓ³⁄² = ħ/√(2mF) this is ħS₂ = (5/48)(ℓ/d)³⁄² = (5/24)δ, using δ = ½(ℓ/d)³⁄².
At d = 10.0 nm, ℓ/d = 0.3365: δ = 0.0976 and ħS₂ = 0.0203 rad, against the leading phase (1/ħ)∫p dx = (2/3)(d/ℓ)³⁄² = 3.42 rad from the turning point, a 0.60% correction.
The two are equal when (5/48)(ℓ/d)³⁄² = (2/3)(d/ℓ)³⁄², so (d/ℓ)³ = 5/32 and d = 0.539ℓ = 1.81 nm: inside that the expansion has no ordering at all.
Check: with ζ = (2/3)(d/ℓ)³⁄², (5/48)(ℓ/d)³⁄² = (5/72)/ζ, and 5/72 is the first coefficient, quoted, of the asymptotic series of the Airy functions, which solve this potential exactly.

Answerℓ = 3.36 nm; at d = 10.0 nm, δ = 0.098 and ħS₂ = 0.020 rad against a leading phase of 3.42 rad; ħS₂ would overtake it at d = 1.81 nm.