Quantum
Quest

Transfer matrix method implementation

To complete the exploration of the transfer matrix method (here a conceptual explanation of the method), I developed a Python implementation.

Introduction to the method

The transfer matrix method is used to analyze wave propagation across various media layers, typically in the fields of quantum mechanics and electronics. By representing each layer with matrices, this method allows for the calculation of transmission and reflection coefficients, essential for understanding wave behavior in layered structures.

Implementation details

Initialization

The simulation begins by defining the wave vectors and matrix variables for each layer. The propagation matrices and boundary condition matrices (Dm for boundaries and Pm for propagation through each layer) are crucial.

The vector k_m, which is defined as

k_m \equiv \sqrt{\frac{2m_{eff{_m}}}{\hbar^2}\left(E - V_m\right)}

is computed, with a factor chosen to ensure the scale are normalized to a range suitable for display.

factor = (2 * m_elec * scale) / (hbar ** 2)
km = np.sqrt(np.vectorize(complex)(factor * meff) * (energy - V))

First boundary

I start from the left, handling the first boundary specially because it sets the stage for the propagation calculations.

delta1 = (km[1] / km[0]) * (meff[0] / meff[1])
Dm[0] = calculate_Dm_matrix(delta1)
Pm[0] = np.eye(2, dtype=complex)
Tm = Dm[0].copy()  # Initialize the total transfer matrix

Two matrices are necessary for each step:

1. calculate_Dm_matrix: Calculates the matrix ensuring compatible boundary conditions between adjacent layers, based on their relative properties.
2. calculate_Pm_matrix: Computes the propagation matrix within each layer, incorporating the exponential terms derived from the wave vector and the layer’s thickness.

On the first layer, only the boundary condition matrix is required, so for consistency the propagation matrix is set to be the identity matrix.

Loop through each layer

I loop through each layer, computing two matrices for each: one for propagation through the layer and another for transitioning between layers.

for n in range(1, len(meff) - 1):
    deltan = (km[n + 1] / km[n]) * (meff[n] / meff[n + 1])
    Dm[n] = calculate_Dm_matrix(deltan)
    Pm[n] = calculate_Pm_matrix(km[n], dist[n])
    Tm = np.dot(Tm, np.dot(Pm[n], Dm[n]))  # Update total transfer matrix sequentially

Compute transmission coefficient

Finally, the transmission coefficient is calculated.

eta = 1 - np.abs(Tm[1, 0])**2 / np.abs(Tm[0, 0])**2

Compute the wavefunction coefficients

Once the matrices are computed it is possible to recompute the coefficients for the wave functions across layers:

Am = np.zeros(len(meff), dtype=complex)
Bm = np.zeros(len(meff), dtype=complex)
Am[len(An) - 1] = A / Tm[0, 0]
Am[0] = A
Bm[0] = Tm[1, 0] * Am[len(An) - 1]
for n in range(len(meff) - 2, 0, -1):
    layer_matrix_n = np.dot(Pm[n], Dm[n])
    Am[n], Bm[n] = np.dot(layer_matrix_n, np.array([Am[n + 1], Bm[n + 1]]))

I presume that the incident wavelength amplitude is known, and since the assumption is that there is no incident from the right and no backward waves in the last region Bm_\text{last}=G=0, that condition is sufficient for computing Am_\text{last}=F=Tm_{11}Am[0]=Tm_{11}A and B = Bm[0] = Tm_{21}F.

The intermediate coefficients are computed starting from the right (it would be also possible equivalently to start from the left as A and B are known), and for each layer the wavefunction is:

\psi_m(x) = A_m e^{ik_m(x-x_{m-1})} + B_m e^{-ik_m(x-x_{m-1})}

and can be plot as piecewise functions.

Conclusion

For further exploration and visual demonstrations of this method applied to specific cases like a double potential barrier, see the visualizations here.