User:Gareth Owen/WKB approximation

In physics, the WKB (Wentzel-Kramers-Brillouin) approximation, also known as WKBJ approximation, is the most familiar example of 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.

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 linear, second-order differential equations, which includes the Schrödinger equation. But since the Schrödinger equation was developed two years later, and Wentzel, Kramers, and Brillouin were apparently unaware of this earlier work, Jeffreys is often neglected credit. Early texts in quantum mechanics contain any number of combinations of their initials, including WBK, BWK, WKBJ and BWKJ.

Derivation
We beginning with a one dimensional, time-independent wave equation in which the local wavenumber varies. Such an equation can typically be written as
 * $$\epsilon^2\frac{d^2}{dx^2}\Psi(x) + K^2(x) \Psi(x) = 0$$,

where K is O(1) and $$\epsilon$$ is small.

We recast the wavefunction as the exponential of another function &Phi; (which is closely related to the action):
 * $$\Psi(x) = \displaystyle{e^{\Phi(x)}}$$

The function &Phi; must then satisfy
 * $$\epsilon^2\left[\Phi''(x) + \left[\Phi'(x)\right]^2\right] + K^2(x) = 0$$

where &Phi;' indicates the derivative of &Phi; with respect to x. Now let us separate $$\Phi'(x)$$ into real and imaginary parts by introducing the real functions A and B:
 * $$\displaystyle\Phi'(x) = A(x) + i B(x)$$

The amplitude of the wavefunction is then $$e^{A(x)}$$ while its phase is $$B(x)$$. The governing equation implies that these functions must satisfy:


 * $$\displaystyle \epsilon^2\left(A'(x) + A^2(x) - B^2(x)\right) + K^2(x) = 0$$

and since the right hand side of the differential equation for &Phi; is real,


 * $$\displaystyle B'(x) + 2 A(x) B(x) = 0.$$

Next we want to find an asymptotic approximation to solve this. That means we expand each function as a power series in $$\epsilon$$. From the equations we can already see that the power series must start with at least an order of $$\epsilon^{-1}$$ to satisfy the real part of the equation.


 * $$A(x) = \sum_{n=-1}^\infty \epsilon^n A_n(x)$$


 * $$B(x) = \sum_{n=-1}^\infty \epsilon^n B_n(x)$$

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


 * $$\displaystyle A_{-1}^2(x) - B^2_{-1}(x) + K^2(x) = 0,$$
 * $$\displaystyle A_{-1}(x) B_{-1}(x) = 0$$

Clearly then, the second of these equations tells us that either $$A_{-1}(x)$$ or $$B_{-1}(x)$$ must be identically zero. Since both functions are real, examination of the first equation tells us that if $$K^2 > 0$$ then $$A_{-1}\equiv0$$ and if $$K^2(x) <0$$, then $$B_{-1} \equiv0$$

Oscillatory Regime
In the former case, $$\displaystyle K^2>0$$ the leading order term


 * $$\displaystyle\frac1\epsilon B_{-1}(x) = \pm \frac1\epsilon K(x),$$

constitutes a rapid variation in phase. Thus, the solutions to the equation are predominately oscillatory in nature. For this case, we can calculate the next order correction. Taking the real and imaginary parts of term of order $$\epsilon$$, and noting that $$A_{-1}=A'_{-1}=0$$ we have
 * $$\displaystyle B_0(x)B_{-1}(x) = 0$$

and
 * $$\displaystyle B'_{-1}(x) + 2 A_0(x)B_{-1}(x) = 0,$$

giving us
 * $$\displaystyle B_0(x) = 0\quad\mathrm{and}\quad A_0(x) = -\frac12\frac{B'_{-1}(x)}{B_{-1}(x)} = -\frac12\frac{K'(x)}{K(x)}$$

Thus, to the first two orders
 * $$\Phi'(x) = \frac{i}\epsilon B_{-1}(x) + A_0(x) = \pm\frac{i}{\epsilon}{K(x)}

- \frac12\frac{K'(x)}{K(x)}$$ Integrating directly, we have
 * $$\Psi(x)=e^{\Phi(x)}=\exp\left(\frac{i}{\epsilon}\int^x K(s)ds - \frac12 \log (K(x))\right) = \frac{\exp\left(\pm\frac{i}\epsilon \int^x{K(s)}ds\right)}{\sqrt{K(s)}}$$

Exponential Regime
When $$K^2(x) < 0$$ (corresponding to an imaginary wavenumber, K) the leading order term is given by
 * $$\frac1\epsilon A_{-1}(x) = \pm\frac1\epsilon|K(x)|,$$

leading to either exponential growth or decay. At the next order, noting that $$B_{-1}(x)=0$$, we have
 * $$\displaystyle A'_{-1}(x) + 2 A_0(x)A_{-1}(x) = 0$$

and
 * $$\displaystyle A_{-1}(x)B_0(x) = 0.$$

Therefore, $$B_0(x)=0$$ and, as in the oscillatory regime,
 * $$A_0(x) = - \frac12 \frac{A_{-1}'(x)}{A_{-1}(x)}.$$

Combining these two terms we obtain
 * $$\Psi(x) = \frac{\exp\left(\pm\frac1\epsilon\int^x|K(s)|ds\right)}{|K(x)|^{1/2}}.$$

It is clear that in neither regime is the approximation valid near $$K(x)=0$$, where the denominator becomes singular, as the assumption that $$K(x)$$ is $$O(1)$$ breaks down. In wave dynamics, particularly optics, the location where this happens is known as a caustic.

Across the Caustic
It is apparent from the denominator, that both of these approximate solutions 'blow up' when the local wavenumber $$|K(x)|$$ passes through zero, and cannot be valid. The approximate solutions that we have found are accurate away from this zero, but inaccurate near to it. We can find accurate approximate solutions near to this zero by approximating $$|K(x)|$$ by its Taylor series

Let's label the zero of $$|K(x)|$$ by $$x_0$$. Now, if $$x$$ is near $$x_0$$, we can write
 * $$\frac{K^2(x)}{\epsilon^2} \simeq \frac{K^2(x_0)}{\epsilon^2} + U (x - x_0)  + \cdots = U (x - x_0)  + O((x-x_0)^2)$$,

where $$U$$ is the derivative of $$K^2(x)/\epsilon^2$$ at $$x_0$$

To first order, one finds
 * $$\frac{d^2}{dx^2} \Psi(x) + U(x - x_0) \Psi(x) = 0.$$

This differential equation is Airy equation, and the solution may be written in terms of Airy functions. (Alternatively, with some trickery, it may be transformed into a Bessel equation of fractional order.) The exact form of the solution depends on the sign of $$U$$. The case when $$U<0$$ corresponds to the oscillatory regime being to the left of $$x_0$$ and the exponential regime to the right. In this case, the solution is given by,
 * $$\displaystyle\Psi(x) = C_1 \mathrm{Ai}(|U|^{1/3}(x-x_0)) + C_2 \mathrm{Bi}(|U|^{1/3}(x-x_0))$$

When $$U>0$$, the locations of the oscillatory and exponential regimes are reversed, and the solution is:
 * $$\displaystyle\Psi(x) = C_1 \mathrm{Ai}(-|U|^{1/3}(x-x_0)) + C_2 \mathrm{Bi}(-|U|^{1/3}(x-x_0))$$

This solution is accurate near the zero, and should connect with the solutions. Thus, we should be able to determine the 2 coefficients $$C_1$$ and $$C_2$$ so that the solutions are identical in the region of overlap, where the are both accurate.