Suppose we have quantum mechanical plane waves of energy $E$ incident upon a one-dimensional potential barrier $V(x)$ with sloping sides.
One can compare the WKB solutions in the three relevant regions to asymptotics of the Airy function solutions to the linearized Schrödinger equation near the turning points to obtain matching formulae, relating $A, B$ (the coefficients of the right- and left-propagating waves, respectively, left of the barrier) to $F, G$ (the coefficients of the right- and left-propagating waves, respectively, right of the barrier). One can then construct the relevant scattering matrix $S$, defined by:
$$\left[\begin{array}{c}
B\\
F
\end{array}\right]=S\left[\begin{array}{c}
A\\
G
\end{array}\right] \, .$$
I have carried out the above computation and found that
$$S=\left[\begin{array}{cc}
i\left(\dfrac{1-4e^{2\gamma}}{1+4e^{2\gamma}}\right) & \dfrac{4e^{\gamma}}{1+4e^{2\gamma}}\\
\dfrac{4e^{\gamma}}{1+4e^{2\gamma}} & i\left(\dfrac{1-4e^{2\gamma}}{1+4e^{2\gamma}}\right)
\end{array}\right]$$
where
$$\gamma=\frac{1}{\hbar}\int_{x_{1}}^{x_{2}}\sqrt{2m\left[V\left(x\right)-E\right]}dx$$
measures the "strength" of the reflecting, classically-forbidden region.
This scattering matrix $S$ is in fact unitary, and hence this scattering process conserves the probability current across the barrier. But why should this be true?
I know that a solution of the Schrödinger equation will necessarily conserve probability current, but here we have both (1) assumed approximate WKB-type solutions and (2) carried out a complicated asymptotic matching procedure. It is not clear to me why, after this approximation and matching procedure, the resulting solution should still conserve probability current.
Said another way: is there an a priori reason one should expect the WKB-matched solutions to conserve probability current (in spite the fact that they are inexact), or is this simply a happy coincidence that falls out of the computation?
Best Answer
I) Recall that the probability current 1D QM is
$$\tag{1} J(x)~:=~\frac{i\hbar}{2m}W(\psi,\psi^{\ast}) (x), $$
where
$$\tag{2} W(\psi,\psi^{\ast}) (x)~:=~\psi(x) \psi^{\prime}(x)^{\ast}- \psi^{\prime}(x)\psi(x)^{\ast}$$
is the Wronskian. Unitarity of the $S$-matrix is equivalent to the statement that
$$ \lim_{x\to -\infty} W(\psi,\psi^{\ast}) (x)~=~\lim_{x\to +\infty} W(\psi,\psi^{\ast}) (x), \tag{3}$$
cf. e.g. eq. (11) in my Phys.SE answer here. A sufficient condition for unitarity is if the Wronskian (2) does not depend on the $x$-position.
II) Firstly, there exists many semiclassical WKB approximations in the literature, e.g. to what order in $\hbar$ are we talking? If we consider some truncated WKB approximation scheme, there is no a priori reason why the unitarity condition (3) should be satisfied.
Secondly, the one-directional nature of the connection formulas often makes it impossible to determine the exponentially growing branch in a classically forbidden side from knowledge about the classically allowed side of a turning point, in particular in situations with many turning points, cf. e.g. Refs. 1-3.
OP is presumably interested in the case of a single potential barrier with two turning points. This case has been complete resolved in the literature using e.g. uniform approximation techniques, and the resulting WKB formula preserves unitarity, see e.g. Ref. 3. However, we shall not consider these methods here.
For the rest of this answer, we will consider the following WKB formula
$$ \psi(x) ~=~A(x)\sum_{\sigma=\pm 1} C_{\sigma}(x;x_0)~ \exp\left(\frac{i}{\hbar}\sigma S(x;x_0)\right)$$ $$ ~=~A(x)\sum_{\sigma=\pm 1} B_{\sigma}(x;x_0)~ \cos\left(\frac{1}{\hbar} S(x;x_0)+\sigma\frac{\pi}{4}\right),\tag{4}$$
written either in terms of an exponential function or a cosine function, where
$$A(x)~:=~\frac{1}{\sqrt{p(x)}} , \qquad S(x;x_0)~:=~ \int_{x_0}^x\!dx^{\prime}~p(x^{\prime}), $$ $$ p(x)~:=~ \sqrt{2m(E-V(x))},\qquad E,V(x)~\in~\mathbb{R}; \tag{5} $$
$x_0\in \mathbb{R}$ is fixed reference point; $\sigma\in\{\pm 1\}$ is a sign; and $B_{\sigma}(x;x_0),C_{\sigma}(x;x_0)\in\mathbb{C}$ are complex piecewise $x$-independent constants. As we shall see in Section VI, the piecewise constants $B_{\sigma}(x;x_0)$ and $C_{\sigma}(x;x_0)$ may jump at turning points because of the metaplectic correction/Maslov index, cf. e.g. this Phys.SE post and references therein. It is implicitly implied that the two square roots in formula (5) are appropriately chosen branches, not double-valued functions.
III) We have no a priori argument why the Wronskian of the WKB approximation (4) should be $x$-independent to all orders in $\hbar$. Nevertheless, it is interesting and instructive to analyze its form. The Wronskian of the WKB approximation (4) becomes
$$ W(\psi,\psi^{\ast}) ~\stackrel{(2)+(4)}{=}~\sum_{\sigma=\pm 1}\sum_{\sigma^{\prime}=\pm 1} \exp\left[\frac{i}{\hbar}\left(\sigma S -\sigma^{\prime}S^{\ast}\right)\right]$$ $$\left[|A|^2(C_{\sigma} C^{\prime\ast}_{\sigma^{\prime}} - C^{\prime}_{\sigma} C^{\ast}_{\sigma^{\prime}}) + C_{\sigma} C^{\ast}_{\sigma^{\prime}}\left\{ (AA^{\prime \ast} - A^{\prime}A^{\ast}) - \frac{i}{\hbar}|A|^2(\sigma p+\sigma^{\prime} p^{\ast} ) \right\}\right] .\tag{6} $$
Firstly, one may argue (using arguments$^1$ similar to what is done below) that the term
$$AA^{\prime \ast} - A^{\prime}A^{\ast}~=~0\tag{7} $$
vanishes in eq. (6). Secondly, the term $ C_{\sigma} C^{\prime\ast}_{\sigma^{\prime}} - C^{\prime}_{\sigma} C^{\ast}_{\sigma^{\prime}} $ is proportional to delta functions with supports at the turning points. Hence, away from the turning points, the Wronskian (6) simplifies to
$$ W^{\neq}(\psi,\psi^{\ast}) ~\stackrel{(6)}{=}~- \frac{i}{\hbar}|A|^2 \sum_{\sigma=\pm 1}\sum_{\sigma^{\prime}=\pm 1} C_{\sigma} C^{\ast}_{\sigma^{\prime}}\left(\sigma p+\sigma^{\prime} p^{\ast} \right)\exp\left[\frac{i}{\hbar}\left(\sigma S -\sigma^{\prime}S^{\ast}\right)\right] $$ $$~=~\frac{1}{\hbar}|A|^2 \sum_{\sigma=\pm 1}\sum_{\sigma^{\prime}=\pm 1} B_{\sigma} B^{\ast}_{\sigma^{\prime}} \left[p\sin\left(\frac{1}{\hbar} S+\sigma\frac{\pi}{4}\right)\cos\left(\frac{1}{\hbar} S^{\ast}+\sigma^{\prime}\frac{\pi}{4}\right)\right. $$ $$\left. -p^{\ast}\cos\left(\frac{1}{\hbar} S+\sigma\frac{\pi}{4}\right)\sin\left(\frac{1}{\hbar} S^{\ast}+\sigma^{\prime}\frac{\pi}{4}\right)\right]. \tag{8} $$
IV) Classically allowed region $E>V(x)$. Then $p$ and $S$ are both real, and we must choose same signs $\sigma=\sigma^{\prime}$ in the exponential version of eq. (8) to get non-zero contributions. The exponential factor is then just $1$. The Wronskian (8) is $x$-independent
$$\tag{9} W^>(\psi,\psi^{\ast}) ~\stackrel{(8)}{=}~- \frac{2i}{\hbar}|A|^2p\sum_{\sigma=\pm 1}\sigma |C_{\sigma}|^2 ~=~ \frac{1}{\hbar}|A|^2p\sum_{\sigma=\pm 1}\sigma B_{\sigma}B^{\ast}_{-\sigma}, $$
because $|A|^2p$ is $x$-independent, cf. eq. (5). The last expression of eq. (9) follows from the trigonometric version of eq. (8) with the help of a trigonometric addition formula.
V) Classically forbidden region $E<V(x)$. Then $p$ and $S$ are both imaginary, and we must choose opposite signs $\sigma=-\sigma^{\prime}$ in the exponential version of eq. (8) to get non-zero contributions. The exponential factor is then just $1$. The Wronskian (8) is $x$-independent
$$\tag{10} W^<(\psi,\psi^{\ast}) ~\stackrel{(8)}{=}~- \frac{2i}{\hbar}|A|^2p\sum_{\sigma=\pm 1}\sigma C_{\sigma}C^{\ast}_{-\sigma} ~=~ \frac{1}{\hbar}|A|^2p\sum_{\sigma=\pm 1}\sigma |B_{\sigma}|^2. $$
(The superscripts $>$ and $<$ refer to the classically allowed and forbidden regions, respectively.)
VI) A turning point $V(x_0)=E \Leftrightarrow p(x_0)=0$. Let us for notational simplicity assume $x_0=0$. Assume that in a neighborhood around $x_0=0$, the potential can be approximated with a linear$^2$ potential
$$\tag{11} V(x)-E~\propto ~x,$$
i.e. one side is classically allowed and the other side is classically forbidden. Then eq. (5) implies
$$\tag{12} p(x)~\propto ~x^{\frac{1}{2}}, \qquad S(x)~\propto ~x^{\frac{3}{2}},\qquad A(x)~\propto ~x^{-\frac{1}{4}}. $$
The TISE of the linear potential (11) becomes the Airy equation. From the asymptotic expansion of the Airy functions, one derives the Airy connection formulas
$$\tag{13} \exp\left(\sigma\frac{2}{3}|x|^{\frac{3}{2}}\right) \quad\longleftrightarrow\quad\frac{3-\sigma}{2}\cos\left( \frac{2}{3}|x|^{\frac{3}{2}}+\sigma\frac{\pi}{4}\right) $$
between the positive and negative $x$-axis, where $\sigma\in\{\pm 1\}$. The Airy connection formulas (13) lead to the WKB connection formulas
$$\tag{14} C^<_{+1}~=~e^{\frac{i\pi}{4}}B^>_{+1}, \qquad 2C^<_{-1}~=~e^{\frac{i\pi}{4}}B^>_{-1}, $$
for the complex coefficients $B_{\sigma},C_{\sigma}\in\mathbb{C}$, if we adapt the following sign conventions
$$\tag{15} p^<~=~ip^>\qquad A^<~=~e^{-\frac{i\pi}{4}}A^>$$
for the square roots in eq. (5). In matrix-form we have
$$\tag{16} C^<~=~\begin{pmatrix} 1 & i \cr \frac{i}{2} &\frac{1}{2} \end{pmatrix} C^>\qquad\text{and}\qquad C^>~=~\begin{pmatrix} \frac{1}{2} & -i \cr -\frac{i}{2} &1 \end{pmatrix} C^<.$$
Note that the connection formulas (14) and (15) imply that the Wronskian
$$\tag{17} W^>(\psi,\psi^{\ast}) ~=~W^<(\psi,\psi^{\ast}) $$
of the WKB approximation (4) does not jump at the turning point, cf. eqs. (9) and (10).
The Airy connection formulas (13) may superficially look like double-directional connection formulas, but that's an illusion. Because of its asymptotic nature, it involves a (possibly illegitimate) implicit choice, which we have made for simplicity.
It is possibly to use more general complex methods involving the method of steepest descent, which leads to the principle of exponential dominance:
"A multiplier $C_+$, or $C_-$ can only change on a good path when crossing a Stokes line on which its exponential is subdominant."
One may show that this principle obeys unitarity, see e.g. Ref. 3 for details.
TL;DR: We have argued using various simplifying assumptions that the WKB formula (4) preserves unitarity.
--
References:
D. Griffiths, Intro to QM, Chapter 8.
A. Galindo & P. Pascual, QM2, Chapter 9.
M.V. Berry & K.E. Mount, Semiclassical approximations in wave mechanics, Rep. Prog. Phys 35 (1972) 315; Chapter 3 & 4.
$^1$ E.g. at a turning point $x_0$, the transformation $x\to -x$ induces a sign flip on the left-hand side of eq. (7). But minus zero is still zero!
$^2$ We leave it to the reader to analyze the cases of higher-order turning-points and the case of an infinite hard wall.