The transfer matrix: a stack of thin films

Between two interfaces a wave travels back and forth, picking up phase on each pass, and every reflection spawns another. Summing those multiple reflections term by term is possible for one layer (it gives the Airy formula of the Fabry–Pérot section) and impractical for twenty. The transfer matrix does the bookkeeping instead: in each layer the field is a forward and a backward wave, and a $2\times 2$ matrix carries either the two amplitudes or the two tangential fields from one side of the layer to the other. The Foundations chapter explains why such a matrix must exist for any stratified geometry; this section writes it out for films, the way Abelès did in 1950, and shows what it is good for.

The characteristic matrix of a layer

Take s polarization (p differs only in the admittance) and a layer $j$ of index $n_j$ and thickness $d_j$ at the in-plane wavenumber $k_x = k_0 n_0\sin\theta_0$ fixed by the incidence. Inside the layer the field is $E_y(z) = a\,e^{ik_z z} + b\,e^{-ik_z z}$ with $k_z = k_0 n_j\cos\theta_j$, and the tangential magnetic field, in units that make the admittance of free space 1, is $H_x(z) = \eta_j\,(a\,e^{ik_z z} - b\,e^{-ik_z z})$ with the tilted admittance $\eta_j = n_j\cos\theta_j$. Eliminate $a$ and $b$ between the values at the two faces, a distance $d_j$ apart, and the result is a matrix relation between the tangential fields themselves:

Abelès' characteristic matrix $M_j$, with the phase thickness $\delta_j$ and, for p polarization, $\eta_j = n_j/\cos\theta_j$. Its determinant is 1 (it is a product of two matrices of unit Wronskian), its trace is $2\cos\delta_j$, and it depends on the layer through $\delta_j$ and $\eta_j$ only: two numbers per layer per wavelength, which is why a stack of forty layers is no harder than one.

A stack and its admittance

The tangential fields are continuous at every interface, so what leaves the back of layer $j$ enters the front of layer $j+1$ unchanged, and the matrices of successive layers simply multiply. In the substrate there is only a transmitted wave, whose fields are in the ratio $H/E = \eta_s$, so for $q$ layers on a substrate

$Y$ is the admittance of the whole stack: seen from the incident medium of admittance $\eta_0$, the coated substrate behaves exactly like a single surface of admittance $Y$, and the single-interface formulas of the Fresnel section apply to it:

This is the formulation of Macleod's book and of coating design in general, and its power is that $Y$ can be followed layer by layer. Macleod's admittance locus plots $Y$ in the complex plane as each layer is added: a lossless layer moves $Y$ along a circle centred on the real axis, through an angle $2\delta$, and a design is a path from $\eta_s$ to the target $Y = \eta_0$ (zero reflectance) or as far from it as possible (a mirror). A designer reads off from the diagram which layer to add next and how thick, which is how the classic antireflection and filter designs below were found long before computers.

The amplitude form: the field inside

The same physics can be written for the wave amplitudes instead of the fields, which is what Byrnes's code and PyStratify's planar engine do. Crossing an interface between media $i$ and $j$ maps the forward and backward amplitudes by an interface matrix built from the Fresnel coefficients, and crossing a layer by a diagonal propagation matrix:

The product $\tilde M = I_{01} P_1 I_{12} P_2 \cdots I_{q s}$ relates the incident-medium amplitudes $(1, r)$ to the substrate's $(t, 0)$, so $r = \tilde M_{21}/\tilde M_{11}$ and $t = 1/\tilde M_{11}$. The two forms are equivalent (a change of basis in each layer), but the amplitude form hands over the forward and backward amplitudes in every layer as a by-product. From them the field at any depth,

is a one-line evaluation, and the power absorbed in layer $j$ is the difference of the Poynting flux through its two faces, the per-layer absorptance of the solar-cell section. The calculator's field profile and layer absorptances come from exactly these amplitudes.

Quarter waves and half waves

Two thicknesses make the characteristic matrix trivial at the design wavelength $\lambda_0$, at normal incidence:

• a half-wave layer ($n d = \lambda_0/2$, $\delta = \pi$) has $M = -I$: it is absent, and the stack reflects as if it were not there (which is why the spacer of a bandpass filter is a half wave);
• a quarter-wave layer ($n d = \lambda_0/4$, $\delta = \pi/2$) has $M = \begin{pmatrix} 0 & i/\eta \\ i\eta & 0 \end{pmatrix}$ and turns an admittance $Y$ into $\eta^2/Y$: an admittance transformer, the optical quarter-wave line.

A single quarter-wave layer of index $n_1$ on a substrate $n_s$ therefore presents $Y = n_1^2/n_s$ and gives

which vanishes for $n_1 = \sqrt{n_0 n_s}$, 1.23 for glass in air. No durable material has so low an index; magnesium fluoride (1.38) comes closest, and leaves about 1.3 % instead of 4 %.

import numpy as np
from pystratify import planar

n = lambda wl: [1.0, 1.38, 1.52]  # incident medium, layers, substrate
d = [np.inf, 0.109, np.inf]  # thicknesses, µm
lam = (400 + 2 * np.arange(201)) * 1e-3  # wavelength, µm
th = np.radians(0)
R, T = {}, {}
for pol in ('s', 'p'):
    rows = [planar.coh_tmm(pol, n(wl), d, th, wl) for wl in lam]
    R[pol], T[pol] = np.array([r['R'] for r in rows]), np.array([r['T'] for r in rows])
R_avg, T_avg = (R['s'] + R['p']) / 2, (T['s'] + T['p']) / 2  # unpolarized

The minimum sits at $\lambda_0 = 4 n_1 d$ and is broad, because $\delta$ changes slowly with $\lambda$. Two or three layers (the V-coat and the quarter–half–quarter design) broaden it further, and a graded index takes the idea to its limit in the moth-eye antireflection section.

Herpin's equivalent layer. Any symmetric combination of layers, $p\,q\,p$ say, has a characteristic matrix of the same form as a single layer, with an equivalent index and phase thickness that depend on wavelength (Herpin 1947; Macleod ch. 6). A combination of two available materials can therefore stand in for an index that no material offers, including the 1.23 above, and a periodic stack is a periodic repetition of one equivalent layer, which is the shortest route to the stop band of the Bragg mirror. The equivalent index becomes imaginary exactly inside the stop band.

Oblique incidence and absorbing layers

Nothing above assumes normal incidence or real indices. At an angle, $\cos\theta_j$ shortens each layer's phase thickness, so every feature of a spectrum moves to shorter wavelengths as the coating is tilted, and the admittances of s and p separate, so the two polarizations respond differently (the p admittance $n/\cos\theta$ grows with angle, the s admittance shrinks). In an absorbing layer $\delta_j$ is complex, $\cos\delta_j$ and $\sin\delta_j$ become hyperbolic, and the wave decays across the layer; the formulas stay valid, and $A = 1 - R - T$ is now positive.

Numerics: where the matrix overflows, and what to do

For a thick absorbing layer, or an evanescent one beyond the critical angle, $e^{i\delta_j}$ contains a real exponential $e^{\mathrm{Im}\,\delta_j}$ that grows without bound: 10 µm of a metal at visible wavelengths gives $e^{+700}$, beyond double precision. The characteristic matrix carries both $e^{+}$ and $e^{-}$ and overflows, and even before that the reflected amplitude is the small difference of two huge numbers. Three remedies are in use:

1. Propagate only decaying waves. Byrnes's formulation orders the exponentials so that the growing one never appears in a product that is later subtracted, and clamps a layer beyond which nothing returns as opaque; this is what pystratify.planar does, and it is why the calculator can hold a millimetre of metal.
2. The scattering matrix. Instead of transfer matrices, combine the reflection and transmission matrices of sub-stacks with Redheffer's star product. Every quantity is then a reflection or transmission coefficient, bounded by one for a passive stack, and no exponential ever grows (Li 1996, who proved the stability). PyStratify's Verification page compares the planar engine with an independent scattering-matrix solver on lossy, oblique and opaque stacks to $10^{-10}$.
3. Logarithmic amplitudes. Carry $\log a_j$ and $\log b_j$ rather than $a_j$ and $b_j$, so that a factor $e^{700}$ is the number 700; this is how PyStratify's sphere and cylinder solvers stay stable to any thickness and order, as the multilayered sphere section describes.

The second trap is subtler: the branch of $\cos\theta_j$ in an absorbing or evanescent layer. Taking the wrong sign makes the "forward" wave grow into the layer, and the error does not announce itself; it shows up as $R > 1$ or as $T$ failing to vanish beyond the critical angle. The forward-decaying branch of the Fresnel section must be chosen in every layer, and the energy-conservation checks ($R + T = 1$ below and beyond the critical angle for lossless stacks) are the tests that catch a wrong branch.