Anisotropic thin-film optics simulation (Berreman 4×4)

Why 4×4

In an isotropic film the s and p polarisations never mix, and a 2×2 calculation is enough. In a birefringent layer with a tilted optic axis, a biaxial crystal or a layer with an off-diagonal dielectric tensor, the two polarisations exchange energy inside the layer and the polarisation state of the reflected and transmitted light changes. Treating that exactly requires the 4×4 (Berreman) formalism, which carries all four tangential field components together.

Phazic solves the 4×4 problem with a block S-matrix and the Redheffer star product, not with transfer matrices — for the same reason as the isotropic engine: thick layers, absorbing layers and evanescent regions stay finite. Uniaxial layers use a closed-form solution; biaxial and tensor layers solve the eigenproblem and separate the modes with invariant-subspace projectors, so degenerate cases such as normal incidence on a c-cut crystal or Δn = 0 do not break the calculation.

Which anisotropy

• Uniaxial — n_o, n_e and the optic-axis direction (θ_c, φ_c). Calcite, sapphire, quartz, LiNbO₃, liquid crystals.
• Biaxial — three principal indices and the rotation of the principal axes (θ_c, φ_c, ψ_c). Nonlinear crystals such as KTP and LBO, oriented organic films.
• Dielectric tensor — off-diagonal components read from a file, including absorbing tensors whose real and imaginary parts have different principal axes.
• Magneto-optic (gyrotropic) — Faraday rotation in non-reciprocal media.
• Absorbing and negative-permittivity media — graphite, hBN, and hyperbolic media in the infrared reststrahlen bands (sapphire, α-MoO₃).
• Crystalline incident and exit media — a coating on a sapphire substrate, a calcite prism. A sample-azimuth control rotates every optic axis together.
• Incoherent anisotropic slabs — a 2×2 coherency matrix that averages the round-trip phase while keeping the phase between modes.

How the anisotropic axes are placed — turn them yourself

The material axes of a layer (x′, y′, z′) are placed in the laboratory frame (x, y, z) by three angles. In this frame z is the stacking direction (incident medium → stack) and x–z is the plane of incidence, so x is the p and y the s direction. The dielectric tensor is carried into the laboratory frame by

so that the optic axis z′ points along (sin θ_c cos φ_c, sin θ_c sin φ_c, cos θ_c). The three angles act on axes that start aligned with the laboratory (x′ = x, y′ = y, z′ = z) and are read from right to left in R:

• ψ_c — turns x′ and y′ about z′ while z′ stays on the normal.
• θ_c (tilt) — tilts about y, so that z′ leaves the normal by θ_c towards +x. 0° is the c-plane, 90° the a- or m-plane.
• φ_c (azimuth) — turns all three axes together about the normal: the angle by which the projection of the optic axis on the surface (dashed) is turned from the plane of incidence towards y, equivalent to rotating the sample in its own plane. The sample-azimuth field in the program adds its angle to φ_c of every anisotropic layer.
• θ (angle of incidence) — between the incident wave and the normal (yellow arc, drawn at 35°). θ and θ_c share the normal as their reference but are independent: θ belongs to the light, θ_c to the anisotropic layer.

Uniaxial layers

One optic axis z′. A field across z′ meets n_o, a field along it n_e. A turn about the optic axis changes nothing, so θ_c and φ_c suffice. The presets are the c-cut (θ_c = 0°), the a- or m-cut (90°) and the r-cut of sapphire (57.6°).

Figure 1 — turn it with the sliders and presets. The numbers below are the third column of R: the optic axis in laboratory components.

Biaxial layers

Three principal indices n_x′, n_y′, n_z′ along three perpendicular axes, and all three angles matter. With the indices sorted n_α < n_β < n_γ, the two optic axes (pink dashed) lie in the plane of the n_α and n_γ axes at ±V from the n_γ axis,

so the optic axes move with the indices and hence with wavelength. These are the wave-normal optic axes (binormals), along which both plane waves see the same index; the ray axes (biradials) differ slightly. Neither needs to coincide with z′: in LBO the angle 2V about the n_γ axis exceeds 90°, and the acute bisectrix is the n_α axis. Typing indices into the fields shows how the two axes close and open.

Figure 2 — KTP and LBO at 550 nm from refractiveindex.info. Making two indices equal merges the optic axes into one: the uniaxial limit.

Both figures use the rotation and optic-axis formulas of the Phazic solver itself, so the directions they show are the directions the calculation uses.

In every module

• Plane wave — R/T/A, polarisation conversion (s↔p), angle and wavelength scans, per-layer absorption, the internal |E(z)|² and absorption density.
• Ellipsometry — generalised-ellipsometry Ψ, Δ (r_pp/r_ss) and fits of thickness and principal indices. Crystal fits are staged from several starting points so that they do not settle in a wrong valley, and an axis index the data cannot determine is reported as such.
• OLED — dipole emission in stacks with anisotropic layers and anisotropic emitters: uniaxial stacks with the optic axis along the normal use the scalar CPS solver directly; other crystals use a dedicated 4×4 dipole engine. Outcoupling efficiency, Purcell factor, channel decomposition and angular distribution are all available.
• Solar cell — J_sc, EQE and depth-resolved generation with anisotropic absorbers and window layers.
• Optimisation — thickness optimisation of stacks that contain anisotropic layers (their structure is preserved).

Stacks without anisotropic layers go through the existing isotropic engine unchanged, so their results are bit-identical — adding this capability changed no number in an existing project.

Validation

• Agreement within about 10⁻¹⁴ with an independently written Berreman 4×4 code for uniaxial, biaxial and tensor stacks.
• Within 5×10⁻¹³ on eleven cases built from measured refractiveindex.info data (sapphire, calcite, quartz, rutile, LiNbO₃, KTP, LBO, graphite, hBN, and the infrared reststrahlen bands of sapphire and α-MoO₃). Wave plates are also compared with the closed-form result.
• The anisotropic OLED path is checked against a Lorentz-reciprocity calculation that does not use the dipole formulae at all.
• Isotropic input reproduces the existing engine within 10⁻¹²; a dipole in an unbounded crystal reproduces the closed-form rates.

Applications

• Wave plates and retarders, liquid-crystal cells and polarisation optics — design and tolerance
• AR coatings on sapphire c-axis windows and crystalline substrates — how substrate birefringence changes coating performance
• OLED outcoupling and angular colour with birefringent organic layers
• Reflection and absorption of 2D materials (hBN, MoS₂, graphite) and hyperbolic media
• Faraday rotation in magneto-optic films; form-birefringent multilayers

Try it

Anisotropy is included in the Designer Pro edition and above (Designer Pro, Display, Energy, Suite). Calculations on isotropic stacks run in every edition as before. Phazic runs on Windows 10 and 11 and is available from the Microsoft Store.

Related: Mixed layers (EMA) · Thickness-tolerance Monte Carlo · OLED dipole emission