# WKB approximation

(Redirected from WKB)

In mathematical physics, the WKB approximation or WKB method is a method 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 wavefunction is recast as an exponential function, semiclassically expanded, and then either the amplitude or the phase is taken to be slowly changing.

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 Wentzel, Kramers, and Brillouin, who all developed it in 1926. In 1923, mathematician Harold Jeffreys had developed a general method of approximating solutions to linear, second-order differential equations, which includes the Schrödinger equation. Even though the Schrödinger equation was developed two years later, Wentzel, Kramers, and Brillouin were apparently unaware of this earlier work, so Jeffreys is often neglected credit. Early texts in quantum mechanics contain any number of combinations of their initials, including WBK, BWK, WKBJ, JWKB and BWKJ. An authoritative discussion and critical survey has been given by R B Dingle.[1]

Earlier references to the method are: Carlini in 1817, Liouville in 1837, Green in 1837, Rayleigh in 1912 and Gans in 1915. Liouville and Green may be said to have founded the method in 1837, and it is also commonly referred to as the Liouville–Green or LG method.[2][3]

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.

## WKB method

Generally, WKB theory is a method for approximating the solution of a differential equation whose highest derivative is multiplied by a small parameter ε. The method of approximation is as follows.

For a differential equation

${\displaystyle \epsilon {\frac {d^{n}y}{dx^{n}}}+a(x){\frac {d^{n-1}y}{dx^{n-1}}}+\cdots +k(x){\frac {dy}{dx}}+m(x)y=0,}$

assume a solution of the form of an asymptotic series expansion

${\displaystyle y(x)\sim \exp \left[{\frac {1}{\delta }}\sum _{n=0}^{\infty }\delta ^{n}S_{n}(x)\right]}$

in the limit δ → 0. The asymptotic scaling of δ in terms of ε 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 Sn(x) in the expansion.

WKB theory is a special case of multiple scale analysis.[4][5][6]

## An example

This example comes from the text of Bender and Orszag.[6] Consider the second-order homogeneous linear differential equation

${\displaystyle \epsilon ^{2}{\frac {d^{2}y}{dx^{2}}}=Q(x)y,}$

where ${\displaystyle Q(x)\neq 0}$. Substituting

${\displaystyle y(x)=\exp \left[{\frac {1}{\delta }}\sum _{n=0}^{\infty }\delta ^{n}S_{n}(x)\right]}$

results in the equation

${\displaystyle \epsilon ^{2}\left[{\frac {1}{\delta ^{2}}}\left(\sum _{n=0}^{\infty }\delta ^{n}S_{n}'\right)^{2}+{\frac {1}{\delta }}\sum _{n=0}^{\infty }\delta ^{n}S_{n}''\right]=Q(x).}$

To leading order (assuming, for the moment, the series will be asymptotically consistent), the above can be approximated as

${\displaystyle {\frac {\epsilon ^{2}}{\delta ^{2}}}S_{0}'^{2}+{\frac {2\epsilon ^{2}}{\delta }}S_{0}'S_{1}'+{\frac {\epsilon ^{2}}{\delta }}S_{0}''=Q(x).}$

In the limit δ → 0, the dominant balance is given by

${\displaystyle {\frac {\epsilon ^{2}}{\delta ^{2}}}S_{0}'^{2}\sim Q(x).}$

So δ is proportional to ε. Setting them equal and comparing powers yields

${\displaystyle \epsilon ^{0}:\quad S_{0}'^{2}=Q(x),}$

which can be recognized as the Eikonal equation, with solution

${\displaystyle S_{0}(x)=\pm \int _{x_{0}}^{x}{\sqrt {Q(t)}}\,dt.}$

Considering first-order powers of ε fixes

${\displaystyle \epsilon ^{1}:\quad 2S_{0}'S_{1}'+S_{0}''=0.}$

This is the unidimensional transport equation, having the solution

${\displaystyle S_{1}(x)=-{\frac {1}{4}}\ln Q(x)+k_{1},}$

where k1 is an arbitrary constant.

We now have a pair of approximations to the system (a pair, because S0 can take two signs); the first-order WKB-approximation will be a linear combination of the two:

${\displaystyle y(x)\approx c_{1}Q^{-{\frac {1}{4}}}(x)\exp \left[{\frac {1}{\epsilon }}\int _{x_{0}}^{x}{\sqrt {Q(t)}}\,dt\right]+c_{2}Q^{-{\frac {1}{4}}}(x)\exp \left[-{\frac {1}{\epsilon }}\int _{x_{0}}^{x}{\sqrt {Q(t)}}\,dt\right].}$

Higher-order terms can be obtained by looking at equations for higher powers of δ. Explicitly,

${\displaystyle 2S_{0}'S_{n}'+S''_{n-1}+\sum _{j=1}^{n-1}S'_{j}S'_{n-j}=0}$

for n ≥ 2.

### Precision of the asymptotic series

The asymptotic series for y(x) is usually a divergent series, whose general term δn Sn(x) starts to increase after a certain value n=nmax. Therefore, the smallest error achieved by the WKB method is at best of the order of the last included term.

For the equation

${\displaystyle \epsilon ^{2}{\frac {d^{2}y}{dx^{2}}}=Q(x)y,}$

with Q(x) <0 an analytic function, the value ${\displaystyle n_{\max }}$ and the magnitude of the last term can be estimated as follows:[7]

${\displaystyle n_{\max }\approx 2\epsilon ^{-1}\left|\int _{x_{0}}^{x_{\ast }}{\sqrt {-Q(z)}}\,dz\right|,}$
${\displaystyle \delta ^{n_{\max }}S_{n_{\max }}(x_{0})\approx {\sqrt {\frac {2\pi }{n_{\max }}}}\exp[-n_{\max }],}$

where ${\displaystyle x_{0}}$ is the point at which ${\displaystyle y(x_{0})}$ needs to be evaluated and ${\displaystyle x_{\ast }}$ is the (complex) turning point where ${\displaystyle Q(x_{\ast })=0}$, closest to ${\displaystyle x=x_{0}}$.

The number nmax can be interpreted as the number of oscillations between ${\displaystyle x_{0}}$ and the closest turning point.

If ${\displaystyle \epsilon ^{-1}Q(x)}$ is a slowly changing function,

${\displaystyle \epsilon \left|{\frac {dQ}{dx}}\right|\ll Q^{2},}$

the number nmax will be large, and the minimum error of the asymptotic series will be exponentially small.

## Application to the Schrödinger equation

The above example may be applied specifically to the one-dimensional, time-independent Schrödinger equation,

${\displaystyle -{\frac {\hbar ^{2}}{2m}}{\frac {d^{2}}{dx^{2}}}\Psi (x)+V(x)\Psi (x)=E\Psi (x),}$

which can be rewritten as

${\displaystyle {\frac {d^{2}}{dx^{2}}}\Psi (x)={\frac {2m}{\hbar ^{2}}}\left(V(x)-E\right)\Psi (x).}$

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

${\displaystyle \Psi (x)=e^{\Phi (x)},}$

so that

${\displaystyle \Phi ''(x)+\left[\Phi '(x)\right]^{2}={\frac {2m}{\hbar ^{2}}}\left(V(x)-E\right),}$

where Φ ' indicates the derivative of Φ with respect to x. This derivative Φ ' can be separated into real and imaginary parts by introducing the real functions A and B,

${\displaystyle \Phi '(x)=A(x)+iB(x).}$

The amplitude of the wavefunction is then

${\displaystyle \exp \left[\int _{x_{0}}^{x}A(x')\,dx'\right],}$

while the phase is

${\displaystyle \int _{x_{0}}^{x}B(x')\,dx'.}$

The real and imaginary parts of the Schrödinger equation then become

${\displaystyle A'(x)+A(x)^{2}-B(x)^{2}={\frac {2m}{\hbar ^{2}}}\left(V(x)-E\right),}$
${\displaystyle B'(x)+2A(x)B(x)=0.}$

Next, the semiclassical approximation is used. This means that each function is expanded as a power series in ħ. From the above equations, it can be seen that the power series must start with at least an order of 1/ħ to satisfy the real part of the equation. In order to achieve a good classical limit, it is necessary to start with as high a power of Planck's constant ħ as possible:

${\displaystyle A(x)={\frac {1}{\hbar }}\sum _{n=0}^{\infty }\hbar ^{n}A_{n}(x),}$
${\displaystyle B(x)={\frac {1}{\hbar }}\sum _{n=0}^{\infty }\hbar ^{n}B_{n}(x).}$

To the zeroth order in this expansion, the conditions on A and B can be written,

${\displaystyle A_{0}(x)^{2}-B_{0}(x)^{2}=2m\left(V(x)-E\right),}$
${\displaystyle A_{0}(x)B_{0}(x)=0\;.}$

The first derivatives A'(x) and B'(x) were discarded, because they include factors of order 1/ħ, higher than the dominant ħ−2.

Then, if the amplitude varies sufficiently slowly as compared to the phase (${\displaystyle A_{0}(x)=0}$), it follows that

${\displaystyle B_{0}(x)=\pm {\sqrt {2m\left(E-V(x)\right)}},}$

which is only valid when the total energy is greater than the potential energy, as is always the case in classical motion.

After the same procedure on the next order of the expansion, it follows that

${\displaystyle \Psi (x)\approx C_{0}{\frac {e^{\theta +i\int \hbar ^{-1}{\sqrt {2m\left(E-V(x)\right)}}\,dx}}{\hbar ^{-1/2}{\sqrt[{4}]{2m\left(E-V(x)\right)}}}}.}$

On the other hand, if it is the phase that varies slowly (as compared to the amplitude), (${\displaystyle B_{0}(x)=0}$) then

${\displaystyle A_{0}(x)=\pm {\sqrt {2m\left(V(x)-E\right)}},}$

which is only valid when the potential energy is greater than the total energy (the regime in which quantum tunneling occurs).

Finding the next order of the expansion yields, as in the example of the previous section,

 ${\displaystyle \Psi (x)\approx {\frac {C_{+}e^{\int \hbar ^{-1}{\sqrt {2m\left(V(x)-E\right)}}\,dx}+C_{-}e^{-\int \hbar ^{-1}{\sqrt {2m\left(V(x)-E\right)}}\,dx}}{\hbar ^{-1/2}{\sqrt[{4}]{2m\left(V(x)-E\right)}}}}.}$

It is evident in the denominator that both of these approximate solutions become singular near the classical turning points, where E = V(x), and cannot be valid. These are the approximate solutions away from the potential hill and beneath the potential hill.

Away from the potential hill, E>V(x), the particle acts similarly to a free wave—the wave-function is oscillating. Beneath the potential hill, E<V(x), the particle undergoes exponential changes in amplitude: barrier penetration.

To complete the derivation, the approximate solutions must be found everywhere and their coefficients matched to comprise a global approximate solution. The approximate solution near the classical turning points E = V(x) is yet to be found.

For a classical turning point x1 and close to E = V(x1), the term ${\displaystyle {\frac {2m}{\hbar ^{2}}}\left(V(x)-E\right)}$ can be expanded in a power series,

${\displaystyle {\frac {2m}{\hbar ^{2}}}\left(V(x)-E\right)=U_{1}\cdot (x-x_{1})+U_{2}\cdot (x-x_{1})^{2}+\cdots \;.}$

To first order, one finds

${\displaystyle {\frac {d^{2}}{dx^{2}}}\Psi (x)=U_{1}\cdot (x-x_{1})\cdot \Psi (x).}$

This differential equation is known as the Airy equation, and the solution may be written in terms of Airy functions,

${\displaystyle \Psi (x)=C_{A}{\textrm {Ai}}\left({\sqrt[{3}]{U_{1}}}\cdot (x-x_{1})\right)+C_{B}{\textrm {Bi}}\left({\sqrt[{3}]{U_{1}}}\cdot (x-x_{1})\right).}$

This solution should connect the far away and beneath solutions. Given the two coefficients on one side of the classical turning point, the 2 coefficients on the other side of the classical turning point can be determined by using this local solution to connect them. Thus, a relationship between ${\displaystyle C_{0},\theta }$ and ${\displaystyle C_{+},C_{-}}$ can be found.

Fortunately, the Airy functions will asymptote into sine, cosine and exponential functions in the proper limits. The relationship can be found to be as follows (often referred to as the"connection formulas"):

${\displaystyle C_{+}=+{\frac {1}{2}}C_{0}\cos {\left(\theta -{\frac {\pi }{4}}\right)},}$
${\displaystyle C_{-}=-{\frac {1}{2}}C_{0}\sin {\left(\theta -{\frac {\pi }{4}}\right)}.}$

Now the global (approximate) solutions can be constructed. The same can be done at the other turning points; assume there is just another one, x2. The expression there, however, will appear different than the one determined above at x1 by a difference in the argument of these trigonometric functions.

Thus, single-valuedness of the resulting wave-function dictates the use of EBK quantization, an effective Bohr-Sommerfeld quantization phase space constraint, modified by a Maslov correction, here 1/2, which thus serves to specify the energy eigenvalues, E,

${\displaystyle \int _{x_{1}}^{x_{2}}{\sqrt {{\frac {2m}{\hbar ^{2}}}\left(E-V(x)\right)}}\,dx=(n+1/2)\pi \hbar ,}$

where ${\displaystyle x_{1},x_{2}}$ are the turning points of the potential discussed, where the integrand vanishes.

(For an estimate of the errors in this approximation, cf. Ch 15.6 of Hall.)

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.[8]

## References

1. ^ R.B. Dingle, Asymptotic Expansions: Their Derivation and Interpretation (Academic Press, 1973).
2. ^ Adrian E. Gill (1982). Atmosphere-ocean dynamics. Academic Press. p. 297. ISBN 978-0-12-283522-3.
3. ^ Renato Spigler & Marco Vianello (1998). "A Survey on the Liouville–Green (WKB) approximation for linear difference equations of the second order". In Saber Elaydi; I. Győri & G. E. Ladas. Advances in difference equations: proceedings of the Second International Conference on Difference Equations : Veszprém, Hungary, August 7–11, 1995. CRC Press. p. 567. ISBN 978-90-5699-521-8.
4. ^ Filippi, Paul (1999). Acoustics: basic physics, theory and methods. Academic Press. p. 171. ISBN 978-0-12-256190-0.
5. ^ Kevorkian, J.; Cole, J. D. (1996). Multiple scale and singular perturbation methods. Springer. ISBN 0-387-94202-5.
6. ^ a b Bender, C.M.; Orszag, S.A. (1999). Advanced mathematical methods for scientists and engineers. Springer. pp. 549–568. ISBN 0-387-98931-5.
7. ^ Winitzki, S. (2005). "Cosmological particle production and the precision of the WKB approximation". Phys. Rev. D. 72: 104011, 14 pp. arXiv:. Bibcode:2005PhRvD..72j4011W. doi:10.1103/PhysRevD.72.104011.
8. ^ Harald J.W. Müller-Kirsten, Introduction to Quantum Mechanics: Schrödinger Equation and Path Integral, 2nd ed. (World Scientific, 2012).