This method involves representing the wave function as a product of a transfer matrix and a state vector, where the matrix relates the wave function and its derivative at one point to another point. The core advantage of this approach is its ability to elegantly handle discontinuities in potential and to compute transmission and reflection coefficients by propagating the wave function across potential barriers or wells, making it highly effective for studying quantum mechanical wave propagation in layered media.
The idea is to approximate the actual potential V(x) with a series of constant steps, and therefore to approximate the problem to a sum of layer with constant potential (for which the solution is known) with appropriate boundary conditions between the various layers.
So, if there are N layer within the material, there are N+2 total layers considering the one before and after the block of material itself.
Let’s consider an incident wave, and there will be then reflected and transmitted waves and for each layer it is possible to compute the ratio between the various amplitudes in a matrix form and progressively move from one layer to the next.
Each layer j will have a potential energy V_j, a thickness d_m and an effective mass m_{eff{_i}}.
In any layer, if E > V_j then the solution has the form:
\psi_i(x) = A_j e^{ik_j(x-x_{j-1})} + B_j e^{-ik_j(x-x_{j-1})}, \quad k_j \equiv \sqrt{\frac{2m_{eff{_i}}}{\hbar^2}\left(E - V_j\right)}
In any layer, if E < V_j then the solution has the form:
\psi_i(x) = A_j e^{\kappa_j(x-x_{j-1})} + B_j e^{-\kappa_j(x-x_{j-1})}, \quad \kappa_j \equiv \sqrt{\frac{2m_{eff{_i}}}{\hbar^2}\left(V_j - E\right)}
These can be carried explicitly or can be combined in a single expression just defining:
k_j \equiv \sqrt{\frac{2m_{eff{_i}}}{\hbar^2}\left(E - V_j\right)}
and:
\psi_i(x) = A_j e^{ik_j(x-x_{j-1})} + B_j e^{-ik_j(x-x_{j-1})}
As for the second case (E < V_j):
ik_j = i\sqrt{\frac{2m_{eff{_i}}}{\hbar^2}\left(E - V_j\right)} = i^2 \sqrt{\frac{2m_{eff{_i}}}{\hbar^2}\left(V_j - E\right)} = - \sqrt{\frac{2m_{eff{_i}}}{\hbar^2}\left(V_j - E\right)} = -\kappa
So, while there is no mathematical difference, in the practical algorithmic implementation using this formulation is simplify the code as it is using the same formulation for both cases.
The next step it to look at the boundary conditions, first for the continuity of the wavefunction \psi, for the layer j:
A_{j}(L) + B_{j}(L) = A_{j+1}(0) + B_{j+1}(0)
As for the continuity of the derivative:
\frac{\mathrm d\psi}{\mathrm dx} = ik(A-B)\psi
Then:
A_{j}(L) - B_{j}(L) = \frac{k_{j+1}}{k_j}\left(A_{j+1}(0) - B_{j+1}(0)\right) = \Delta_j \left(A_{j+1}(0) - B_{j+1}(0)\right)
This condition can be extended to semiconductors as \frac{1}{m_{eff}}\frac{\mathrm d\psi}{\mathrm dx} if needed.
Combining these two equations:
\begin{aligned} A_{j}(L)&=A_{j+i}(0)\left[\frac{1+\Delta_j}{2}\right] + B_{j+i}(0)\left[\frac{1-\Delta_j}{2}\right] \\ B_{j}(L)&=A_{j+i}(0)\left[\frac{1-\Delta_j}{2}\right] + B_{j+i}(0)\left[\frac{1+\Delta_j}{2}\right] \end{aligned}
That can be written in matrix form:
\begin{bmatrix}A_{j}(L) \\ B_{j}(L) \end{bmatrix} = \begin{bmatrix} \frac{1+\Delta_j}{2} &\frac{1-\Delta_j}{2} \\ \frac{1-\Delta_j}{2} & \frac{1+\Delta_j}{2}\end{bmatrix} \begin{bmatrix}A_{j+i}(0) \\ B_{j+i}(0) \end{bmatrix}
or compactly:
\mathbf v_j(L) = \mathbf D_j \mathbf v_{j+1}(0)
The final step is to add the matrix which is explicitly computing the value at L which is:
\begin{bmatrix}A_{j}(0) \\ B_{j}(0) \end{bmatrix} = \begin{bmatrix} e^{-k_jd_j} & 0 \\ 0 & e^{k_jd_j} \end{bmatrix} \begin{bmatrix}A_{j}(L) \\ B_{j}(L) \end{bmatrix} = \mathbf P_j \mathbf v_{j}(L)
So the full matrix transfer matrix can be written as:
\mathbf v_1 = \mathbf T\, \mathbf v_{N+2}
where:
\mathbf T = D_1 \prod_{i=2}^{N+1} P_iD_i
Once the transfer matrix is computed for a specific energy E, it is possible to compute the fraction of the incident particles; we assume that there is no wave incident from the right and no backward waves, then in this case:
\begin{bmatrix}A \\ B \end{bmatrix} = \begin{bmatrix} T_{11} & T_{12} \\ T_{21} & T_{22} \end{bmatrix} \begin{bmatrix} F \\ 0 \end{bmatrix}
Then:
\begin{aligned} A & = T_{11}F \\ B & = T_{21}F \end{aligned}
then the fraction transmitted is:
\eta = \frac{|A|^2 - |B|^2}{|A|^2} = 1 - \frac{|T_{21}|^2}{|T_{11}|^2}
This method is simple to program and therefore is well suited for numerical calculation and therefore an useful technique for a generic one dimensional potential. Since the coefficients are known, it is possible to compute the wavefunction inside the layers as the sum of the forward and backward wavefunctions:
\psi_i(x) = A_j e^{ik_j(x-x_{j-1})} + B_j e^{-ik_j(x-x_{j-1})}, \quad k_j \equiv \sqrt{\frac{2m_{eff{_i}}}{\hbar^2}\left(E - V_j\right)}
With this framework, it is possible to use the transfer matrix to find eigenvalues of bound states, and are the one for which:
T_{11} = 0
as in this case A_1 = 0, A_{N+2}\ne 0. This condition can be used for a numerical calculation by varying E.
For a quantum barrier to be in resonance, the condition generally involves the matching of the energy of an incoming wave to specific discrete energy levels that allow the wave to transmit through the barrier with high probability. Specifically, for a potential barrier in quantum mechanics, resonance occurs when the phase shift of the wave function across the barrier leads to constructive interference, effectively matching the boundary conditions that maximize the transmission coefficient. In mathematical terms, resonance in quantum tunneling through a barrier can be described as follows: - The transmission coefficient T approaches unity (T \approx 1). - The energy E of the incident particle matches one of the quasi-bound state energies within the barrier. This can be expressed as a condition on the phase \delta of the wave function, such that:
\delta(E) = n\pi
where n is an integer, indicating that the phase increment across the barrier leads to constructive interference. Additionally, for a simple rectangular barrier, the resonance condition can be derived from the matching of the wave vector inside the barrier k_2 to the barrier width L as:
k_2 L = n\pi
Resonance occurs when these conditions lead to a perfect constructive interference of the wave function across the barrier, greatly enhancing the probability of tunneling. That can be demonstrated using the transfer matrix method. Let’s consider two layers of zero potential with a barrier of length d and potential V_0 such as:
k_b d = n\pi
Then the transition matrix takes the simple form:
\mathbf T = \mathbf D_1 \mathbf P_2 \mathbf D_2
The easier to compute it the propagation matrix:
\mathbf P_2 = \begin{bmatrix} e^{-ik_bd} & 0 \\ 0 & e^{ik_b d} \end{bmatrix}
since e^{-i\theta} = \cos(\theta) + i \sin(\theta) with \theta = n\pi:
e^{\pm i n\pi} = \cos(\pm n\pi) + i\sin(\pm n\pi) = \cos(n\pi) + 0 = (-1)^n
then:
\mathbf P_2 = \begin{bmatrix} e^{-in\pi} & 0 \\ 0 & e^{in\pi} \end{bmatrix} = (-1)^n \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} = (-1)^n \mathbf I
as the potential outside the barrier as the same (k), it is possible to define a quantity a = \frac{k}{k_b}, which gives for \mathbf D_2:
\mathbf D_2 = \begin{bmatrix} \frac{1+a}{2} &\frac{1-a}{2} \\ \frac{1-a}{2} & \frac{1+a}{2}\end{bmatrix} = \frac{1}{2}\begin{bmatrix} 1+a & 1-a \\ 1-a & 1+a \end{bmatrix}
and for \mathbf D_1:
\mathbf D_1 = \begin{bmatrix} \frac{1+\frac{1}{a}}{2} & \frac{1-\frac{1}{a}}{2} \\ \frac{1-\frac{1}{a}}{2} & \frac{1+\frac{1}{a}}{2}\end{bmatrix} = \frac{1}{2a}\begin{bmatrix} a+1 & a-1 \\ a-1 & a+1 \end{bmatrix}
Then it is possible to compute the transition matrix:
\begin{aligned} \mathbf T & = \mathbf D_1 \mathbf P_2 \mathbf D_2 = \mathbf D_1 (-1)^n \mathbf I \mathbf D_2 = (-1)^n \mathbf D_1 \mathbf D_2 \\ & = (-1)^n\frac{1}{4a}\begin{bmatrix} a+1 & a-1 \\ a-1 & a+1 \end{bmatrix}\begin{bmatrix} 1+a & 1-a \\ 1-a & 1+a \end{bmatrix} \\ & = (-1)^n\frac{1}{4a}\begin{bmatrix} (a+1)^2 -(a-1)^2 & (a+1)(a-1) +(1-a)(a+1) \\ (1-a)(a+1) + (a+1)(a-1) & -(a-1)^2 + (a+1)^2 \end{bmatrix} \\ & = (-1)^n\frac{1}{4a}\begin{bmatrix} 4a & 0 \\ 0 & 4a \end{bmatrix} = (-1)^n\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} = (-1)^n \mathbf I \end{aligned}
The fraction transmitted is:
\eta = 1 - \frac{|T_{21}|^2}{|T_{11}|^2} = 1
so all the particles are transmitted as expected.
As, example, it is possible to show the results for a double barrier of equal width between two regions of zero potential energy.
The first graph shows \Re\{\psi(x)\}, the real part of the wavefunction, the second the probability density.
Increasing the energy it is possible to notice clearly the resonance is occurring.
This example is quite straight forward and potentially could be solved analytically, but the power of numerical method is that than be easily and quickly adapted for different potentials that be complex and for which an analytical solution does not exist.