Transfer and Scattering Matrices

Inspect the interface and propagation matrices of the same physical stack. Compare transfer and scattering representations and test a deliberately ill-conditioned transfer product. This coherent, isotropic calculation uses planar electromagnetic boundary conditions.

Two ways to chain layers

For one film you can add up the bounces by hand. For twenty layers you cannot. Matrix methods keep track of two waves in every layer, one travelling forward and one backward, and describe each interface and each layer by a 2×2 matrix. Chaining the matrices chains the layers. This page shows the matrices for a real stack and what happens when the simplest way of chaining them breaks.

Why two waves? At a fixed wavelength and angle, the field in each layer is one wave going forward plus the wave reflected back from further down, so two amplitudes, $E^+$ and $E^-$ , describe it completely. At an interface the tangential electric and magnetic fields must match. The layer's optical admittance $Y$ sets the magnetic field per unit electric field; for TE it equals $q$ , the normal wave number in units of $2\pi/\lambda$ , as defined in the calculation section.

Transfer matrices

A transfer matrix gives the two wave amplitudes on the left of an element from those on its right. At the boundary between layers $j$ and $j+1$ , and across layer $j$ of thickness $d_j$ ,

\begin{gathered}I_{j,j+1}=\frac12\begin{pmatrix}1+\rho & 1-\rho\\ 1-\rho & 1+\rho\end{pmatrix}\\ P_j=\begin{pmatrix}1/p & 0\\ 0 & p\end{pmatrix}\\ \rho=\frac{Y_{j+1}}{Y_j},\qquad p=e^{iq_j2\pi d_j/\lambda}\end{gathered}

Multiplying them in physical order, $M=I_{01}P_1I_{12}\cdots$ , gives one matrix for the whole stack. With nothing arriving from the far side, the stack's amplitudes follow directly:

r=\frac{M_{21}}{M_{11}},\qquad t=\frac{1}{M_{11}}

It is compact and fast.

Why it breaks

Look at $1/p$ in the propagation matrix. In a transparent layer $|p|=1$ and nothing unusual happens. In an absorbing layer, or where the wave is evanescent, $p$ decays exponentially with thickness, so $1/p$ grows exponentially. The product then carries numbers like $e^{+1000}$ alongside numbers like $e^{-1000}$ . Double precision stops near $10^{308}$ , so beyond some thickness the growing entry overflows, or $p$ itself underflows to zero and $1/p$ cannot be formed. In longer stacks, and when internal fields are solved from the product, huge entries must also cancel almost exactly, and rounding destroys that cancellation.

Scattering matrices

A scattering matrix instead relates the waves leaving an element to the waves arriving at it. Each wave is referenced at the side where it enters an element, so a layer contributes its propagation factor $p$ and never $1/p$ . In a passive layer $|p|\le 1$ , so the numbers stay bounded. Combining two scattering matrices takes more algebra than a matrix product, because the amplitudes bouncing between them have to be eliminated. The solver on these pages uses an equivalent reflection recursion built from the same bounded factors; the transfer product is shown for comparison.

Derivation from Maxwell's equations

Start from the single-interface derivation: fields vary as $\exp(-i\omega t)$ ; in layer $j$ the normal wave number is $q_jk_0$ ( $k_0=2\pi/\lambda$ , $\lambda$ the vacuum wavelength), and each partial wave's tangential $H$ is $\pm Y_0Y_j$ times its tangential $E$ —plus for forward waves, minus for backward, times an overall polarization sign that cancels at every boundary. Here $Y_0$ is the vacuum admittance and $Y_j$ is defined in the calculation section. Below, $H$ is in units of $Y_0$ .

Fields in one layer. The layers are stacked along $z$ , increasing away from the incident medium; “forward” means travelling toward $+z$ . In layer $j$ the tangential fields are a forward and a backward wave. Measured from a reference plane $z_j$ , with $\phi=q_jk_0(z-z_j)$ ,

\begin{aligned}E(z)&=E_j^+e^{i\phi}+E_j^-e^{-i\phi}\\ H(z)&=Y_j\left(E_j^+e^{i\phi}-E_j^-e^{-i\phi}\right)\end{aligned}

For transfer matrices each $E_j^\pm$ is evaluated at the plane the matrix connects: at interface $j,j+1$ , layer $j$ 's amplitudes are taken at its exit and layer $j+1$ 's at its entrance.

Interface. At the boundary between layers $j$ and $j+1$ , tangential $E$ and $H$ are continuous. Dividing the magnetic condition by $Y_j$ and writing $\rho=Y_{j+1}/Y_j$ ,

\begin{aligned}E_j^++E_j^-&=E_{j+1}^++E_{j+1}^-\\ E_j^+-E_j^-&=\rho\left(E_{j+1}^+-E_{j+1}^-\right)\end{aligned}

Adding and subtracting gives

\begin{pmatrix}E_j^+\\ E_j^-\end{pmatrix}=I_{j,j+1}\begin{pmatrix}E_{j+1}^+\\ E_{j+1}^-\end{pmatrix}

with $I_{j,j+1}$ exactly the interface matrix given above.

Propagation. Across layer $j$ of thickness $d_j$ , the forward wave gains the factor $p=e^{iq_jk_0d_j}$ and the backward wave, travelling the other way, gains the same factor from exit to entrance. So $E^+_\text{exit}=pE^+_\text{entrance}$ and $E^-_\text{entrance}=pE^-_\text{exit}$ , and entrance amplitudes follow from exit amplitudes through $P_j=\operatorname{diag}(1/p,\,p)$ .

The whole stack. Chaining these gives $M=I_{01}P_1I_{12}\cdots$ , which maps the exit-medium amplitudes at the last interface to the incident-medium amplitudes at the first. Send in a unit wave, $E_0^+=1$ , with nothing arriving from the exit side. Then $E_0^-=r$ and the exit medium carries only $t$ :

\begin{gathered}\begin{pmatrix}1\\ r\end{pmatrix}=M\begin{pmatrix}t\\ 0\end{pmatrix}\\ t=\frac{1}{M_{11}},\qquad r=M_{21}t=\frac{M_{21}}{M_{11}}\end{gathered}

Scattering form. The same two continuity equations can be solved for the waves leaving the interface, $E_j^-$ and $E_{j+1}^+$ , in terms of those arriving, $E_j^+$ and $E_{j+1}^-$ . Substitute $E_{j+1}^+$ from the first continuity equation into the second and solve for the outgoing waves:

\begin{gathered}\begin{pmatrix}E_j^-\\ E_{j+1}^+\end{pmatrix}=S_{j,j+1}\begin{pmatrix}E_j^+\\ E_{j+1}^-\end{pmatrix}\\ S_{j,j+1}=\begin{pmatrix}r_{j,j+1} & 1-r_{j,j+1}\\ t_{j,j+1} & -r_{j,j+1}\end{pmatrix}\\ r_{j,j+1}=\frac{Y_j-Y_{j+1}}{Y_j+Y_{j+1}}=\frac{1-\rho}{1+\rho}\\ t_{j,j+1}=1+r_{j,j+1}\end{gathered}

The first row is the wave leaving to the left, the second the wave leaving to the right; the first column is the wave arriving from the left. The off-diagonal $1-r_{j,j+1}$ is the transmission from the right and $-r_{j,j+1}$ the reflection seen from the right. In terms of the transfer matrix $I=I_{j,j+1}$ ,

\begin{gathered}S_{j,j+1}=\frac{1}{I_{11}}\begin{pmatrix}I_{21} & \det I\\ 1 & -I_{12}\end{pmatrix}\\ \det I=\rho\end{gathered}

The two transmissions differ, $(1-r_{j,j+1})/t_{j,j+1}=\rho$ , because the amplitudes are tangential electric fields, not power-normalized; between transparent media the transmitted power is the same in both directions. For a layer, reference each arriving wave at the face where it enters. The propagation relations above then give the layer block $S_j$ :

\begin{gathered}\begin{pmatrix}E^-_\text{entrance}\\ E^+_\text{exit}\end{pmatrix}=S_j\begin{pmatrix}E^+_\text{entrance}\\ E^-_\text{exit}\end{pmatrix}\\ S_j=\begin{pmatrix}0 & p\\ p & 0\end{pmatrix}\end{gathered}

(left = entrance, right = exit.) Only $p$ appears, never $1/p$ . For a passive layer the square root can be taken with $\operatorname{Im}q_j\ge 0$ —the decaying branch the solver uses—so $|p|=e^{-k_0d_j\operatorname{Im}q_j}\le 1$ .

Three things to try

Read the matrices. The Matrices result tab opens by default. Select the film in the editor: its propagation matrix and the interface on its exit side are highlighted. In the default glass film the film matches the glass beneath it, so $I_{12}$ is the identity. Switch between transfer and scattering representations for the same stack.

A coating. Load the quarter-wave antireflection example, designed for normal incidence. The direct product and the stable solution agree; at the page's 30° default the residual reflection is about 0.12% for TE.

Break it on purpose. The ill-conditioned example is a 20 µm film with a complex index of $0.3+6i$ . Across that thickness $p$ underflows to zero, the transfer product cannot be formed and the page reports the overflow. The scattering solution still returns a finite answer: about 97% reflection at the default settings.

UNIFORM LAYERS

Matrix-method editor

Loading 3D…

Loading geometry…

Physical scale
Drag to orbit · click a region to inspect

Stable scattering solution · revision 0

Interface and propagation matrices

\begin{pmatrix}E^+\\ E^-\end{pmatrix}_{\text{left}}=M\begin{pmatrix}E^+\\ E^-\end{pmatrix}_{\text{right}},\qquad M=I_{01}\,P_1\,I_{12}\cdots

Stack product

M=\begin{pmatrix}-0.821+1.029\,i & 0.1974+0.2474\,i\\ 0.1974-0.2474\,i & -0.821-1.029\,i\end{pmatrix}

Written in physical order from the incident side. M acts on the exit-side amplitudes, so the last interface is applied to them first. Amplitudes are tangential electric fields; these are wave matrices, not ray-transfer ABCD matrices. Select a matrix to select the region it belongs to; each interface belongs to the region on its incident side.

Direct transfer product: $r = -0.2404,\quad t = -0.4737-0.5938\,i$
Stable scattering solution: $r = -0.2404$

Calculation and evidence

Time dependence is $\exp(-i\omega t)$ . Amplitudes are tangential $E_y$ (TE) or $E_x$ (TM), with

\begin{aligned}q&=\sqrt{\varepsilon-n_0^2\sin^2\theta}\\ Y&=\begin{cases}q & \text{TE}\\ \varepsilon/q & \text{TM}\end{cases}\end{aligned}

Interfaces use Fresnel amplitudes; propagation uses $\exp(iq\,2\pi d/\lambda)$ . Reflection recursion composes scattering responses using only decaying exponentials in passive finite layers. Backward fields are anchored at each layer’s exit.

The exterior media must be transparent and positive-index. Internal layers can be passive and absorbing. All layers are coherent, uniform, isotropic and nonmagnetic. An exact zero longitudinal wave number is reported as a singular amplitude basis; nearby angles remain calculable. No roughness, incoherent thick-layer approximation or beam envelope is included.

Independent fixtures use Steven Byrnes’s tmm 0.2.0 for oblique TE/TM, absorbing films, evanescent layers and multilayers. Powers agree within $10^{-8}$ in the fixtures. Checks also cover boundary fields, antireflection and zero-thickness limits. The Python fixture generator is retained with the source.

History of the method

In 1950 Florin Abelès published a study of the propagation of sinusoidal electromagnetic waves in stratified media in the Annales de Physique. In 1972 Dwight Berreman, building on work by Billard and by Teitler and Henvis, gave a 4×4 matrix formulation that also covers anisotropic layers, where the two polarizations couple and each layer carries four waves.

With evanescent or strongly absorbing layers, growing and decaying exponentials coexist in the transfer product and the calculation becomes unstable. In 1988 D. Y. K. Ko and J. R. Sambles proposed a scattering matrix that treats every such wave as decaying. Scattering-matrix formulations later became popular in rigorous coupled-wave analysis (RCWA) of periodically corrugated multilayers. More recent work still reports double-precision transfer products overflowing for absorbing multilayers.

References

Steven J. Byrnes, Multilayer optical calculations, 2016, revision 5 (2020): field conventions, transfer matrices, flux and absorption. The browser implementation uses scattering recursion for numerical stability; the Python implementation supplies independent fixtures.

Florin Abelès, Recherches sur la propagation des ondes électromagnétiques sinusoïdales dans les milieux stratifiés, Annales de Physique 12, 596–640 (1950), continued in pp. 706–782. Cited by title and venue only; the text was not consulted.

Dwight W. Berreman, Optics in stratified and anisotropic media: 4×4-matrix formulation, Journal of the Optical Society of America 62(4), 502 (1972). Attribution follows Bay et al.; the original text was not consulted.

D. Y. K. Ko and J. R. Sambles, Scattering matrix method for propagation of radiation in stratified media: attenuated total reflection studies of liquid crystals, Journal of the Optical Society of America A 5(11), 1863 (1988). Attribution follows Bay et al.; the original text was not consulted.

Mélanie M. Bay, Silvia Vignolini and Kevin Vynck, PyLlama: a stable and versatile Python toolkit for the electromagnetic modelling of multilayered anisotropic media, Computer Physics Communications 273, 108256 (2022); arXiv version. Its introduction is the source for the Berreman and Ko–Sambles attributions, the evanescent-mode instability and the later use of scattering matrices in RCWA.

Luis Eduardo Puente-Díaz et al., Stable calculation of optical properties of large non-periodic dissipative multilayered systems, Superlattices and Microstructures 145, 106629 (2020); arXiv version: the source for transfer-matrix instability and double-precision overflow with absorption.

Return to a single interface →

Calculation source code · Implementation and validation notes