Optics

Schematic

This note documents my implementation of the transfer matrix method (TMM) and the scattering matrix method (SMM) for plane-wave propagation and scattering in one-dimensional multilayer structures. The focus is on a consistent field-amplitude definition and on numerical robustness when dealing with thick layers, evanescent waves, or large refractive-index contrasts. All MATLAB codes discussed here are available in my GitHub repository.

The physical problem considered throughout this note is the propagation of a monochromatic plane wave in a stratified medium, where each layer is homogeneous and isotropic, and interfaces are planar and parallel. Both TE and TM polarizations are supported. The formulation closely follows standard electromagnetic boundary conditions but is written in a way that is convenient for numerical implementation and for reconstructing the field everywhere inside the structure.

Field representation in multilayer structures

In each layer, the electric field is written as a superposition of forward- and backward-propagating plane waves along the direction. With the reference position chosen at the right interface of each layer, the field takes the form

$E(z) = A_m \, e^{i k_{m,z} (z - z_m)} + B_m \, e^{-i k_{m,z} (z - z_m)},$

where and denote the complex amplitudes evaluated at the right boundary of layer . This choice of reference is important, because it fixes how propagation and interface relations are implemented and ensures consistency between the TMM and SMM formulations used in the code.

Transfer matrix formulation

The transfer matrix method relates the field amplitudes on the left side of the structure to those on the right side through a product of interface and propagation matrices. At each interface, the continuity conditions of the tangential electric and magnetic fields can be written in matrix form as

$\begin{pmatrix} A_{m-1} \\ B_{m-1} \end{pmatrix} = D_{m-1}^{-1} D_m \begin{pmatrix} \mathcal{A}_m \\ \mathcal{B}_m \end{pmatrix},$

where is the polarization-dependent interface matrix of layer . The propagation inside layer with thickness is described by

$\begin{pmatrix} \mathcal{A}_m \\ \mathcal{B}_m \end{pmatrix} = P_m \begin{pmatrix} A_m \\ B_m \end{pmatrix}, \qquad P_m = \begin{pmatrix} e^{-i k_{m,z} d_m} & 0 \\ 0 & e^{+i k_{m,z} d_m} \end{pmatrix}.$

By cascading all interfaces and propagation matrices, one obtains the global transfer matrix

$\begin{pmatrix} A_1 \\ B_1 \end{pmatrix} = M_{1\rightarrow N+1} \begin{pmatrix} \mathcal{A}_{N+1} \\ \mathcal{B}_{N+1} \end{pmatrix},$

which directly relates the fields in the input and output semi-infinite regions. With the usual boundary conditions $A1=1\mathcal{B}{N+1}=0B1=1\mathcal{A}{N+1}=0$ for right incidence, reflection and transmission coefficients can be obtained.

While the transfer matrix method is conceptually simple, it is well known that it may suffer from numerical instabilities when exponentially growing and decaying factors coexist. This motivates the use of the scattering matrix method.

Scattering matrix formulation

The scattering matrix method reformulates the problem in terms of incoming and outgoing waves at each side of a layer or a composite structure. For a generic subsystem, the scattering relation is written as

$\begin{pmatrix} B_L \\ A_R \end{pmatrix} = \begin{pmatrix} S_{11} & S_{12} \\ S_{21} & S_{22} \end{pmatrix} \begin{pmatrix} A_L \\ B_R \end{pmatrix},$

where and are the incoming waves from the left and right, and and are the corresponding outgoing waves.

For a single homogeneous layer of thickness $dme^{\pm i k{m,z} d_m}$. Importantly, the phase accumulation inside the layer always appears through the propagation factor and therefore depends explicitly on wavelength and incidence angle, ensuring physical correctness.

When multiple layers are present, the total scattering matrix is obtained by cascading individual scattering matrices. Given two subsystems with scattering matrices and , the combined scattering matrix is given by

$\begin{aligned} S^{(12)}_{11} &= S^{(1)}_{11} + S^{(1)}_{12} \left(I - S^{(2)}_{11} S^{(1)}_{22}\right)^{-1} S^{(2)}_{11} S^{(1)}_{21}, \\ S^{(12)}_{12} &= S^{(1)}_{12} \left(I - S^{(2)}_{11} S^{(1)}_{22}\right)^{-1} S^{(2)}_{12}, \\ S^{(12)}_{21} &= S^{(2)}_{21} \left(I - S^{(1)}_{22} S^{(2)}_{11}\right)^{-1} S^{(1)}_{21}, \\ S^{(12)}_{22} &= S^{(2)}_{22} + S^{(2)}_{21} \left(I - S^{(1)}_{22} S^{(2)}_{11}\right)^{-1} S^{(1)}_{22} S^{(2)}_{12}. \end{aligned}$

By recursively applying this formula, one can build the global scattering matrix of an arbitrary multilayer stack in a numerically stable way.

Field reconstruction from the scattering matrix

Once the total scattering matrix is known, the incoming and outgoing wave amplitudes at the outer boundaries are directly related. For left incidence, one sets and , and solves for and . To reconstruct the field inside the structure, the same scattering relations are applied locally, layer by layer, to back-propagate the known boundary amplitudes and recover the coefficients and defined at each layer’s right interface. These coefficients can then be used to evaluate the exact field distribution everywhere inside the multilayer.

Implementation and examples

The MATLAB functions provided in the GitHub repository implement both TMM and SMM using a unified field-amplitude convention. The function CoeAB_layer_TMM computes the layer-resolved coefficients using the transfer matrix method, while the scattering-matrix-based implementation （CoeAB_layer_SMM）provides improved numerical stability for thick or resonant structures. Additional routines allow direct reconstruction of spatial field profiles.

Several benchmark examples are included, such as Bragg mirrors and standing-wave formation in photoresist layers, where the numerical results agree well with RCWA calculations and classic analytical results reported in the literature.

FieldDistribution

Closing remarks

This note serves both as a theoretical reference and as a practical guide to the accompanying MATLAB code. By presenting the transfer matrix and scattering matrix methods within a consistent framework, it becomes straightforward to switch between the two approaches depending on the numerical requirements of the problem. Readers interested in the implementation details or in extending the code are encouraged to consult the GitHub repository linked in the blog.