In mathematical physics, the WKB approximation or WKB method is a technique for finding approximate solutions to linear differential equations with spatially varying coefficients. It is typically used for a semiclassical calculation in quantum mechanics in which the wave function is recast as an exponential function, semiclassically expanded, and then either the amplitude or the phase is taken to be changing slowly.

The name is an initialism for Wentzel–Kramers–Brillouin. It is also known as the LG or Liouville–Green method. Other often-used letter combinations include JWKB and WKBJ, where the "J" stands for Jeffreys.

Brief history

This method is named after physicists Gregor Wentzel, Hendrik Anthony Kramers, and Léon Brillouin, who all developed it in 1926. In 1923,

Earlier appearances of essentially equivalent methods are: Francesco Carlini in 1817,

The important contribution of Jeffreys, Wentzel, Kramers, and Brillouin to the method was the inclusion of the treatment of turning points, connecting the evanescent and oscillatory solutions at either side of the turning point. For example, this may occur in the Schrödinger equation, due to a potential energy hill.

Formulation

The WKB method approximates the solution of a differential equation whose highest derivative is multiplied by a small parameter <math>\varepsilon</math>. For a differential equation<math display="block"> \varepsilon \frac{d^ny}{dx^n} + a(x)\frac{d^{n-1}y}{dx^{n-1 + \cdots + k(x)\frac{dy}{dx} + m(x)y= 0,</math>assume a solution of the form of an asymptotic series expansion<math display="block"> y(x) \sim \exp\left[\frac{1}{\delta}\sum_{n=0}^{\infty} \delta^n S_n(x)\right]</math>in the limit <math>\delta \rightarrow 0</math>. The asymptotic scaling of <math>\delta</math> in terms of <math>\varepsilon</math> will be determined by the equation. See the example below.

Substituting the above ansatz into the differential equation and cancelling out the exponential terms allows one to solve for an arbitrary number of terms <math>S_n (x)</math>in the expansion.

WKB theory is a special case of multiple scale analysis.

An example

This example comes from the text of Carl M. Bender and Steven Orszag.<math display="block">n_\max \approx \frac{2}{\varepsilon} \left| \int_{x_0}^{x_\ast} \sqrt{-Q(z)}\,dz \right| , </math><math display="block">\delta^{n_\max}S_{n_\max}(x_0) \approx \sqrt{\frac{2\pi}{n_\max e^{-n_\max}, </math>where <math>x_0</math> is the point at which <math>y(x_0)</math> needs to be evaluated and <math>x_{\ast}</math> is the (complex) turning point where <math>Q(x_{\ast}) = 0</math>, closest to <math>x = x_0</math>.

The number <math>n_\max</math> can be interpreted as the number of oscillations between <math>x_0</math> and the closest turning point.

If <math>\varepsilon^{-1}Q(x)</math> is a slowly changing function,

<math display="block">\varepsilon\left| \frac{dQ}{dx} \right| \ll Q^2 , ^{\text{[might be }Q^{3/2}\text{?]</math>

the number <math>n_\max</math> will be large, and the minimum error of the asymptotic series will be exponentially small.

Application in non-relativistic quantum mechanics

thumb|WKB approximation to the indicated potential. Vertical lines show the turning points

thumb|Probability density for the approximate wave function. Vertical lines show the turning points

The above example may be applied specifically to the one-dimensional, time-independent Schrödinger equation,<math display="block">-\frac{\hbar^2}{2m} \frac{d^2}{dx^2} \Psi(x) + V(x) \Psi(x) = E \Psi(x),</math>which can be rewritten as<math display="block">\frac{d^2}{dx^2} \Psi(x) = \frac{2m}{\hbar^2} \left( V(x) - E \right) \Psi(x).</math>Approximation away from the turning points

The wave function can be rewritten as the exponential of another function (closely related to the action), which could be complex,

<math display="block">\Psi(\mathbf x) = e^{i S(\mathbf{x}) / \hbar}, </math>so that its substitution in Schrödinger's equation gives:<math display="block">i\hbar \nabla^2 S(\mathbf x) - \left(\nabla S(\mathbf x)\right)^2 = 2m \left( V(\mathbf x) - E \right),</math>Next, the semi-classical approximation is used. This means that each function is expanded as a power series in <math>\hbar</math>:<math display="block">S = S_0 + \hbar S_1 + \hbar^2 S_2 + \cdots </math>Substituting in the equation, and only retaining terms up to first order in <math>\hbar</math>, one gets<math display="block">\left(\nabla S_0+\hbar \nabla S_1\right)^2 - i\hbar\left(\nabla^2 S_0\right) = 2m\left(E-V(\mathbf x)\right), </math>which gives the following two relations:<math display="block">\begin{align}

\left(\nabla S_0\right)^2 = 2m \left(E - V(\mathbf x)\right) &= \left(p(\mathbf x)\right)^2 \\[1ex]

2\nabla S_0 \cdot \nabla S_1 - i \nabla^2 S_0 &= 0.

\end{align}</math>These can be solved for one-dimension systems. The first equation can be solved by<math display="block">S_0(x) = \pm \int \sqrt{ 2m \left( E - V(x)\right) } \,dx=\pm\int p(x) \,dx, </math>and the second equation computed for the possible values of the above, is generally expressed as:<math display="block">\Psi(x) \approx C_+ \frac{ e^{+ \frac i \hbar \int p(x)\,dx} }{\sqrt{|p(x)| + C_- \frac{ e^{- \frac i \hbar \int p(x)\,dx} }{\sqrt{|p(x)| </math>Thus, the resulting wave function in first order WKB approximation is presented as,

In the classically allowed region, namely the region where <math>V(x) < E</math> the integrand in the exponent is imaginary and the approximate wave function is oscillatory. In the classically forbidden region <math>V(x) > E</math>, the solutions are growing or decaying. It is evident in the denominator that both of these approximate solutions become singular near the classical turning points, where <math>E = V(x) </math>, and cannot be valid. (The turning points are the points where the classical particle changes direction.)

Hence, when <math>E > V(x)</math>, the wave function can be chosen to be expressed as:<math display="block">\Psi(x') \approx \frac{1}{\sqrt{|p(x)| \left[ C \cos\left(\frac 1 \hbar \int \left|p(x)\right| dx + \alpha\right) + D \sin\left(- \frac 1 \hbar \int \left|p(x)\right| dx +\alpha\right)\right] </math>and for <math>V(x) > E</math>,<math display="block">\Psi(x') \approx \frac{ C_{+} e^{- \frac{1}{\hbar} \int |p(x)|\,dx{\sqrt{|p(x)| + \frac{ C_{-} e^{+ \frac{1}{\hbar} \int |p(x)|\,dx} }{ \sqrt{|p(x)|} } . </math>The integration in this solution is computed between the classical turning point and the arbitrary position <math>x' </math>.

Validity of WKB solutions

From the condition:<math display="block">\left(S_0'(x)\right)^2 - \left(p(x)\right)^2 + \hbar \left(2 S_0'(x)S_1'(x)-iS_0(x)\right) = 0 </math>

It follows that: <math display="inline">\hbar\left| 2 S_0'(x)S_1'(x)\right| + \hbar\left| i S_0(x)\right| \ll \left|(S_0'(x))^2\right| + \left| (p(x))^2\right| </math>

For which the following two inequalities are equivalent since the terms in either side are equivalent, as used in the WKB approximation:

<math display="block">\begin{align}

\hbar \left| S_0(x)\right| &\ll \left|(S_0'(x))^2\right| \\

2\hbar \left| S_0'S_1' \right| &\ll \left|(p'(x))^2\right|

\end{align}</math>

The first inequality can be used to show the following:

<math display="block">\begin{align}

\hbar \left| S_0(x)\right| &\ll \left|p(x)\right|^2 \\[6pt]

\frac{1}{2} \frac{\hbar}{|p(x)|} \left|\frac{dp^2}{dx}\right| &\ll \left|p(x)\right|^2 \\[6pt]

\lambda \left|\frac{dV}{dx}\right| &\ll \frac{\left|p\right|^2}{m}\\

\end{align}</math>

where <math display="inline">|S_0'(x)|= |p(x)| </math> is used and <math display="inline">\lambda(x) </math> is the local de Broglie wavelength of the wave function. The inequality implies that the variation of potential is assumed to be slowly varying. This condition can also be restated as the fractional change of <math display="inline">E-V(x) </math> or that of the momentum <math display="inline">p(x) </math>, over the wavelength <math display="inline">\lambda </math>, being much smaller than <math display="inline">1 </math>.

Similarly it can be shown that <math display="inline">\lambda(x) </math> also has restrictions based on underlying assumptions for the WKB approximation that:<math display="block">\left|\frac{d\lambda}{dx}\right| \ll 1 </math>which implies that the de Broglie wavelength of the particle is slowly varying.<math display="block">\begin{align}

\Psi(x) &= C_A \operatorname{Ai}\left( \sqrt[3]{U_1} \cdot (x - x_1) \right) + C_B \operatorname{Bi}\left( \sqrt[3]{U_1} \cdot (x - x_1) \right) \\[6pt]

&= C_A \operatorname{Ai}\left( u \right) + C_B \operatorname{Bi}\left( u \right).

\end{align}</math>Although for any fixed value of <math>\hbar</math>, the wave function is bounded near the turning points, the wave function will be peaked there, as can be seen in the images above. As <math>\hbar</math> gets smaller, the height of the wave function at the turning points grows. It also follows from this approximation that:<math display="block">\begin{align}

\frac{1}{\hbar}\int p(x) \, dx &= \sqrt{U_1} \int \sqrt{x-a}\, dx \\

&= \frac 2 3 \left[\sqrt[3]{U_1} \left(x-a\right)\right]^{\frac 3 2}

= \frac 2 3 u^{\frac 3 2}.

\end{align}</math>

Connection conditions

It now remains to construct a global (approximate) solution to the Schrödinger equation. For the wave function to be square-integrable, we must take only the exponentially decaying solution in the two classically forbidden regions. These must then "connect" properly through the turning points to the classically allowed region. For most values of , this matching procedure will not work: The function obtained by connecting the solution near <math>+\infty</math> to the classically allowed region will not agree with the function obtained by connecting the solution near <math>-\infty</math> to the classically allowed region. The requirement that the two functions agree imposes a condition on the energy , which will give an approximation to the exact quantum energy levels.[[File:WKB approximation example.svg|thumb|WKB approximation to the indicated potential. Vertical lines show the energy level and its intersection with potential shows the turning points with dotted lines. The problem has two classical turning points with <math>U_1 < 0</math> at <math>x=x_1

</math> and <math>U_1 > 0</math> at <math>x=x_2

</math>.]]The wavefunction's coefficients can be calculated for a simple problem shown in the figure. Let the first turning point, where the potential is decreasing over x, occur at <math>x=x_1

</math> and the second turning point, where potential is increasing over x, occur at <math>x=x_2

</math>. Given that we expect wavefunctions to be of the following form, we can calculate their coefficients by connecting the different regions using Airy and Bairy functions.

<math display="block">\begin{align}

\Psi_{V>E}(x) &\approx u^{-\frac{1}{4 \left[ A \exp\left(\tfrac 2 3 u^\frac{3}{2}\right) + B \exp\left(-\tfrac 2 3 u^\frac{3}{2}\right) \right] \\[6pt]

\Psi_{E>V}(x) &\approx u^{-\frac{1}{4 \left[C \cos\left(\tfrac 2 3 u^\frac{3}{2} - \alpha \right) + D \sin\left(\tfrac 2 3 u^\frac{3}{2} - \alpha\right) \right]\\

\end{align} </math>

First classical turning point

For <math>U_1 < 0</math> ie. decreasing potential condition or <math>x=x_1 </math> in the given example shown by the figure, we require the exponential function to decay for negative values of x so that wavefunction for it to go to zero. Considering Bairy functions to be the required connection formula, we get:

<math display="block">\begin{align}

\operatorname{Bi}(u) &\to -\frac{1}{\sqrt \pi}\frac{1}{\sqrt[4]{u \sin\left(\frac 2 3 |u|^{\frac 3 2} - \frac \pi 4\right) & \text{where} \quad u \to -\infty\\[1ex]

\operatorname{Bi}(u) &\to \frac{1}{\sqrt \pi}\frac{1}{\sqrt[4]{u \exp\left(\frac 2 3 u^{\frac 3 2}\right) & \textrm{where} \quad u \to +\infty

\end{align}</math>

We cannot use Airy function since it gives growing exponential behaviour for negative x. When compared to WKB solutions and matching their behaviours at <math>\pm \infty </math>, we conclude:

<math>A=-D=N </math>,

<math>B=C=0 </math> and <math>\alpha = \frac \pi 4 </math>.

Thus, letting some normalization constant be <math>N </math>, the wavefunction is given for increasing potential (with x) as:

It is possible to show that after piecing together the approximations in the various regions, one obtains a good approximation to the actual eigenfunction. In particular, the Maslov-corrected Bohr–Sommerfeld energies are good approximations to the actual eigenvalues of the Schrödinger operator. Specifically, the error in the energies is small compared to the typical spacing of the quantum energy levels. Thus, although the "old quantum theory" of Bohr and Sommerfeld was ultimately replaced by the Schrödinger equation, some vestige of that theory remains, as an approximation to the eigenvalues of the appropriate Schrödinger operator.

General connection conditions

Thus, from the two cases the connection formula is obtained at a classical turning point, <math>x = a</math>: Since the classical particle's velocity goes to zero at the turning points, it spends more time near the turning points than in other classically allowed regions. This observation accounts for the peak in the wave function (and its probability density) near the turning points.

Applications of the WKB method to Schrödinger equations with a large variety of potentials and comparison with perturbation methods and path integrals are treated in Müller-Kirsten.

Examples in quantum mechanics

Although WKB potential only applies to smoothly varying potentials,

Bound states within 2 rigid wall

The potential of such systems can be given in the form:

<math display="block">V(x) = \begin{cases}

\infty & \text{if } x > x_2 \\

V(x) & \text{if } x_2 \geq x \geq x_1\\

\infty & \text{if } x < x_1 \\

\end{cases} </math>

where <math display="inline">x_1 < x_2 </math>.

For <math display="inline">E \geq V(x) </math> between <math display="inline">x_1 </math> and <math display="inline">x_2 </math> which are thus the classical turning points, by considering approximations far from <math display="inline">x_1 </math> and <math display="inline">x_2 </math> respectively we have two solutions:

<math display="block">\begin{align}

\Psi_{\text{WKB(x) &= \frac{A}{\sqrt{|p(x)| \sin\left(\frac 1 \hbar \int_x^{x_1} |p(x)| dx \right) \\

\Psi_{\text{WKB(x) &= \frac{B}{\sqrt{|p(x)| \sin\left(\frac 1 \hbar \int_x^{x_2} |p(x)| dx \right)

\end{align} </math>

Since wavefunctions must vanish at <math display="inline">x_1 </math> and <math display="inline">x_2 </math>. Here, the phase difference only needs to account for <math>n \pi</math> which allows <math>B = (-1)^n A </math>. Hence the condition becomes:

<math display="block">\int_{x_1}^{x_2} \sqrt{2m \left( E-V(x)\right)}\,dx = n\pi \hbar </math>

where <math display="inline">n = 1,2,3,\cdots </math> but not equal to zero since it makes the wavefunction zero everywhere.

<math display="block">E = {\left(3\left(n-\frac 1 4\right)\pi\right)^{\frac 2 3} \over 2}(mg^2\hbar^2)^{\frac 1 3}. </math>

This result is also consistent with the use of equation from bound state of one rigid wall without needing to consider an alternative potential.

Quantum tunneling

The potential of such systems can be given in the form:

<math display="block">V(x) = \begin{cases}

0 & \text{if } x < x_1 \\

V(x) & \text{if } x_2 \geq x \geq x_1\\

0 & \text{if } x > x_2 \\

\end{cases} </math>

where <math display="inline">x_1 < x_2 </math>.

Its solutions for an incident wave is given as

<math display="block">\psi(x) = \begin{cases}

A e^{i k_0 x} + B e^{-i k_0 x} & \text{if } x < x_1 \\[1ex]

\frac{C}{\sqrt{|p(x)|\exp\left(-\frac 1 \hbar \int_{x_1}^{x} |p(x)| dx \right) & \text{if } x_2 \geq x \geq x_1\\[1ex]

D e^{i k_0 x} & \text{if } x > x_2

\end{cases}</math>

with <math>k_0 = p_0/\hbar</math>, where the wavefunction in the classically forbidden region is the WKB approximation but neglecting the growing exponential. This is a fair assumption for wide potential barriers through which the wavefunction is not expected to grow to high magnitudes.

By the requirement of continuity of wavefunction and its derivatives, the following relation can be shown:<math display="block">\frac {|D|^2} {|A|^2} = \frac{4}{(1+{a_1^2}/{p_0^2} )} \frac{a_1}{a_2}\exp\left(-\frac 2 \hbar \int_{x_1}^{x_2} |p(x')| dx'\right) </math>

where <math>a_1 = |p(x_1)|</math> and <math>a_2 = |p(x_2)| </math>.

Using <math display="inline">\mathbf J(\mathbf x,t) = \frac{i\hbar}{2m} \left(\psi^* \nabla\psi - \psi\nabla\psi^*\right) </math> we express the values without signs as:

<math display="block">\begin{align}

J_{\text{inc. &= \tfrac{\hbar}{2m} \left(\tfrac{2p_0}{\hbar}|A|^2\right) \\

J_{\text{ref. &= \tfrac{\hbar}{2m} \left(\tfrac{2p_0}{\hbar}|B|^2\right) \\

J_{\text{trans. &= \tfrac{\hbar}{2m} \left(\tfrac{2p_0}{\hbar}|D|^2\right)

\end{align} </math>

Thus, the transmission coefficient is found to be:

<math display="block">T = \frac {|D|^2} {|A|^2} = \frac{4}{\left(1 + {a_1^2}/{p_0^2} \right)} \frac{a_1}{a_2}\exp\left(-\frac 2 \hbar \int_{x_1}^{x_2} |p(x')| dx'\right) </math>

where <math>a_1 = |p(x_1)|</math> and <math>a_2 = |p(x_2)| </math>. The result can be stated as <math display="inline">T \sim ~ e^{-2\gamma} </math> where <math display="inline">\gamma = \int_{x_1}^{x_2} |p(x')| dx' </math>.

Further reading

  • (An application of the WKB approximation to the scattering of radio waves from the ionosphere.)