The Transfer Matrix Method

To investigate quantum mechanical tunneling one must extract the transmission coefficient from the solution to the one-dimensional, time-independent Schrödinger equation. The transmission coefficient is the ratio of the flux of particles that penetrate a potential barrier to the flux of particles incident on the barrier. It is related to the probability that tunneling will occur. Consider a rectangular potential barrier of height V0 as shown in Figure 1. The general solution of the Schrödinger equation in each region is

$\displaystyle \psi_1 (x)$ = $\displaystyle A e^{ik_1 x} + B e^{-ik_1x} \quad
(x< 0)$ (2)
$\displaystyle \psi_2 (x)$ = $\displaystyle C e^{ik_2 x} + D e^{-ik_2x} \quad
(0\le x \le a)$ (3)
$\displaystyle \psi_3 (x)$ = $\displaystyle F e^{ik_1 x} + G e^{-ik_1x} \quad
(a < x)$ (4)

where $k_1 = \sqrt{2mE/\hbar^2}$ and $k_2 = \sqrt{2m(E-V_0)/\hbar^2}$. This solution can be rewritten as a vector dot product so, for example, in Region 2 of Figure 1

\begin{displaymath}\psi_2 (x) = \pmatrix{ e^{ik_2x} & e^{-ik_2x} \cr}
\pmatrix{ C \cr
D \cr }
= \pmatrix{ e^{ik_2x} & e^{-ik_2x} \cr}
\phi_2
\end{displaymath} (5)

where the $\phi_i$ are the coefficient vectors representing the wave function in each region. To generate relationships among the coefficients in Equations 2-4 one requires the wave function and its first derivative to be continuous at the boundaries of each region in Figure 1. At x=0 this leads to two expressions.
A+B = $\displaystyle C+D \hfil$ (6)
i k1A - i k1 B = i k2 C - i k2 D (7)

In matrix notation Equations 6-7 can be expressed in the following way.
$\displaystyle \left ( \begin{array}{cc}
1 & 1 \\
ik_1 & -ik_1
\end{array} \right )
\left ( \begin{array}{c}
A\\
B
\end{array} \right )$ = $\displaystyle \left ( \begin{array}{cc}
1 & 1 \\
ik_2 & -ik_2
\end{array} \right )
\left ( \begin{array}{c}
C\\
D
\end{array} \right )$ (8)
$\displaystyle {\bf m} \left ( \begin{array}{c}
A\\
B
\end{array} \right )$ = $\displaystyle {\bf n} \left ( \begin{array}{c}
C\\
D
\end{array} \right )$ (9)

Using the definition of the matrix inverse yields the following result

\begin{displaymath}\begin{array}{c}
\left ( \begin{array}{c}
A\\
B
\end{arr...
...( \begin{array}{c}
C\\
D
\end{array} \right )
\end{array}\end{displaymath} (10)

which can be expressed as

\begin{displaymath}\phi_1 = \left ( \begin{array}{c} A\\ B \end{array} \right ) ...
...c} C \\ D \end{array} \right )
= {\bf d_{12}} \phi_2 \qquad .
\end{displaymath} (11)

The matrix, $\bf d_{12}$, is known as the discontinuity matrix and `connects' the coefficient vectors $\phi_1$ and $\phi_2$ in Region 1 and Region 2.

The wave function and its derivative must also be continuous at x=a. At this point, consider a new coordinate system such that the transition from Region 2 to Region 3 takes place at $x^\prime =0$. By analogy with Equation 11 one can show that,

\begin{displaymath}\phi^\prime_2 = {\bf d_{21}} \phi^\prime_3
\end{displaymath} (12)

where the $\phi^\prime_i$ are the coefficient vectors in the new coordinate system in the equivalent regions of the potential energy curve and $\bf d_{21}$ has the same form as $\bf d_{12}$ in Equation 3 except for the interchange of the indices. To exploit this result, one must relate the coefficient vectors in the primed coordinate system to the ones in the original coordinates. The original coordinate system is transformed such that $x^\prime = x - a$. The new wave function in Region 2 must satisfy

\begin{displaymath}\psi_2(x) = \psi^\prime_2(x-a)
\end{displaymath} (13)

which can be written in matrix form and rearranged to yield the following result (recall Equation 5).

\begin{displaymath}\psi_2(x) =\left ( \begin{array}{cc} e^{ik_2(x-a)} & e^{-ik_2...
...\prime e^{-ik_2a} \\
D^\prime e^{ik_2a} \end{array} \right )
\end{displaymath} (14)

The row vector on the right hand side of Equation 14 is the same as the row vector in Equation 5 so the coefficient vectors representing $\phi_2(x)$ and $\phi^\prime_2(x)$ are related by

\begin{displaymath}\phi_2 = \left ( \begin{array}{c}C\\ D \end{array} \right )
...
...ime \\ D^\prime \end{array} \right )
= {\bf p_2}\phi_2^\prime
\end{displaymath} (15)

where $\bf p_2$ is the propagation matrix in Region 2. Similarly, one can show that

\begin{displaymath}\phi^\prime_3
=
\pmatrix{e^{ik_1a} & 0 \cr
0 &e^{-ik_1a} \cr}
\pmatrix{F\cr G\cr}
=
{\bf p_{-1}} \phi_3
\end{displaymath} (16)

where $ {\bf p_{-1}} \phi_3$ shifts the wave function back to the original coordinate system. Combining Equations 11,12, 15, and 16 and setting G=0 since there are no incoming waves in Region 3 one obtains

\begin{displaymath}\phi_1 = \pmatrix{A\cr B\cr}
= {\bf d_{12} p_2 d_{21} p_{-1}} \phi_3
= {\bf t} \phi_3
= \pmatrix{t_{11}F\cr t_{21}F\cr}
\end{displaymath} (17)

where $\bf t$ is the transfer matrix relating the coefficient vectors in Regions 1 and 3. The transmission coefficient is then

\begin{displaymath}T = {\vert Fe^{ik_1 x}\vert^2 \over \vert Ae^{ik_1 x}\vert^2}
= {1 \over \vert t_{11}\vert^2}
\end{displaymath} (18)

This treatment of the rectangular potential barrier problem can be extended to potential barriers of arbitrary shape. Consider the radioactive $\alpha $-decay of $\rm ^{212}Po$. The potential barrier is shown as the solid curve in Figure 2.

 
Figure: 50pt 50pt Potential energy curve for an $\alpha $ particle in the force field of $\rm ^{208}Pb$, the daughter nucleus of the $\rm ^{212}Po$ decay.
\begin{figure}
\begin{center}
\par\vskip -0.3in
\epsfxsize=3.0in
\leavevmode
\epsfbox{gilfoyle_f2.eps}
\par
\par\end{center}\end{figure}

The central portion of the curve is the Coulomb potential V(x) = Z1 Z2 e2/x where x is the distance between the nuclear centers and the product of the charges is Z1 Z2 e2. It is divided into a sequence of adjacent barriers (the dot-dashed lines)lying between the nuclear radius, x0, and the classical turning point, xmax, for an $\alpha $ particle of total energy, E. The potential energy is taken to be zero outside these limits in the manner of Condon and Gurney.12 One can now use the propagation and discontinuity matrices to relate the wave function inside the barrier to the wave function outside (x>xmax). For the configuration shown in Figure 2 one chooses the origin at the nuclear radius and the two wave functions are related by

\begin{displaymath}\phi_{inside} = {\bf d_{01} p_1 d_{12} p_2 d_{23} p_3 d_{30} p_{-0}}
\phi_{outside}
\end{displaymath} (19)

where $\bf d_{01}$ is the discontinuity matrix between the region where V(x) = 0 and V(x) = V1, $\bf p_1$ is the propagation matrix where V(x) = V1, and so on. The propagation matrix $\bf p_{-0}$ returns the wave function to the appropriate coordinate system. The last matrix $\bf p_{-0}$ is unnecessary for calculating the transmission coefficient since it changes the coordinates, but does not change the ratio of the coefficients. The adequacy of treating the potential energy curve in Figure 2 as a sequence of adjacent rectangular barriers will improve as the number of barriers increases and should converge to some limiting value. The transmission coefficient will be extracted from the transfer matrix using Equation 18.




1998-09-14