NAME

Physics::Ellipsometry::VASE::TMM - Transfer Matrix Method for multilayer thin film optics

SYNOPSIS

use PDL;
use Physics::Ellipsometry::VASE::TMM qw(tmm_reflect psi_delta
                                      tmm_graded tmm_interface);
use Physics::Ellipsometry::VASE::Dispersion qw(cauchy_nk);

my $lambda = sequence(500) * 2 + 300;   # 300–1298 nm
my $theta  = ones($lambda) * 70;        # 70° incidence

# --- Single-layer example: SiO2 on Si ---

# Ambient (air)
my $N_air = ones($lambda) + i() * zeroes($lambda);

# SiO2 film (Cauchy model, transparent)
my ($n_sio2, $k_sio2) = cauchy_nk($lambda, 1.46, 0.003, 0.0);
my $N_sio2 = $n_sio2 + i() * $k_sio2;

# Si substrate (tabulated values, simplified here)
my $N_si = 3.87 * ones($lambda) + i() * 0.02 * ones($lambda);

my ($psi, $delta) = psi_delta(
    $lambda, $theta,
    [ $N_air, $N_sio2, $N_si ],   # media: ambient, film, substrate
    [ 100.0 ],                     # thickness: 100 nm SiO2
);

print "Psi range:   ", $psi->min,  " – ", $psi->max,  " deg\n";
print "Delta range: ", $delta->min, " – ", $delta->max, " deg\n";

# --- Raw reflection coefficients ---
my ($rp, $rs) = tmm_reflect(
    $lambda, $theta,
    [ $N_air, $N_sio2, $N_si ],
    [ 100.0 ],
);
my $reflectance_s = abs($rs)**2;
my $reflectance_p = abs($rp)**2;

DESCRIPTION

Implements the 2×2 transfer matrix method (TMM) for computing the optical response of planar multilayer thin-film stacks. This is the standard algorithm used in spectroscopic ellipsometry to predict Psi (Ψ) and Delta (Δ) from a physical model of the sample.

Physical background

When polarised light reflects from a surface, the p-polarised (parallel to the plane of incidence) and s-polarised (perpendicular) components experience different amplitude and phase changes described by the complex Fresnel reflection coefficients rp and rs.

Ellipsometry measures the ratio

ρ = rp / rs = tan(Ψ) · exp(iΔ)

For a simple ambient/substrate interface the Fresnel equations give rp and rs directly. When one or more thin films are present, multiple reflections within each layer create interference effects that modify the overall reflection. The transfer matrix method handles this systematically by representing each layer as a 2×2 matrix and multiplying them together.

Algorithm

For a stack of N media (ambient = 0, layers 1..N-2, substrate = N-1):

Snell's law determines the complex propagation angle in each layer:
```
N_0 · sin(θ_0) = N_j · sin(θ_j)    for all j
```
Fresnel coefficients are calculated at every interface for both polarisations.
Starting from the deepest interface (last layer → substrate), the Airy formula recursively combines each interface reflection with the phase accumulated traversing the layer:
```
r = (r_ij + r_below · e^{+2iβ}) / (1 + r_ij · r_below · e^{+2iβ})
```
where β = (2π/λ) · N_j · d_j · cos(θ_j) is the single-pass phase thickness of layer j.
The final rp and rs at the ambient surface give Ψ and Δ.

CONVENTIONS

Phase convention

Physics sign convention: e^{−iωt} time dependence, producing e^{+2iβ} phase accumulation for forward propagation through a layer.

Fresnel rp (Verdet convention)

rp = (N_f · cosθ_i − N_i · cosθ_f) / (N_f · cosθ_i + N_i · cosθ_f)

This convention yields Δ ≈ 180° for bare dielectrics at moderate angles, consistent with WVASE® and most modern ellipsometry software.

Fresnel rs

rs = (N_i · cosθ_i − N_f · cosθ_f) / (N_i · cosθ_i + N_f · cosθ_f)

Delta mapping

Δ = −arg(ρ)   mapped to [0°, 360°)

matching the convention used in WVASE® and RefEllips.

FUNCTIONS

tmm_reflect

my ($rp, $rs) = tmm_reflect($lambda_nm, $theta_deg,
                            \@N_layers, \@d_nm);

Calculates the complex Fresnel reflection coefficients for a complete multilayer stack.

Arguments:

$lambda_nm: PDL piddle of wavelengths in nanometres (npts elements).
$theta_deg: PDL piddle of angles of incidence in degrees (same size as $lambda_nm).
\@N_layers: Arrayref of complex refractive index piddles [N_0, N_1, ..., N_s], one per medium. N_0 is the ambient (usually air, n=1) and N_s is the substrate. Each piddle must have the same length as $lambda_nm.
\@d_nm: Arrayref of layer thicknesses in nanometres [d_1, d_2, ..., d_{s-1}]. There is no thickness entry for the ambient or the substrate (semi-infinite half-spaces), so this array has two fewer elements than \@N_layers.

Returns a list ($rp, $rs) of complex PDL piddles.

# Two-layer stack: SiO2 (100 nm) + Ta2O5 (50 nm) on Si
my ($rp, $rs) = tmm_reflect($lambda, $theta,
    [ $N_air, $N_ta2o5, $N_sio2, $N_si ],
    [ 50.0, 100.0 ],
);

psi_delta

my ($psi, $delta) = psi_delta($lambda_nm, $theta_deg,
                              \@N_layers, \@d_nm,
                              delta_ref => $ref_delta);

Convenience wrapper around "tmm_reflect" that returns the ellipsometric angles Ψ and Δ in degrees.

Ψ = atan( |ρ| )          in [0°, 90°]
Δ = −arg(ρ)              in [0°, 360°)

The optional delta_ref parameter provides a reference Δ piddle; the calculated Δ is shifted by multiples of 360° to minimise the distance to the reference, avoiding discontinuities at the 0°/360° boundary during fitting.

Returns a list ($psi_deg, $delta_deg) of real PDL piddles.

# Multi-angle measurement
my @angles = (65, 70, 75);
for my $ang (@angles) {
    my $theta = ones($lambda) * $ang;
    my ($psi, $delta) = psi_delta(
        $lambda, $theta,
        [ $N_air, $N_film, $N_sub ],
        [ $thickness ],
    );
    printf "Angle %d: Psi(500nm)=%.2f  Delta(500nm)=%.2f\n",
           $ang, $psi->at(100), $delta->at(100);
}

tmm_graded

my ($rp, $rs) = tmm_graded($lambda_nm, $theta_deg,
                            $N_ambient, $N_substrate,
                            $n_func, $d_total,
                            n_slices => 20);

Models a graded layer whose optical constants vary continuously with depth. The layer is discretized into n_slices sublayers, each with uniform optical constants evaluated at the sublayer midpoint.

This approach handles:

- Composition gradients (e.g., SiO2 transitioning to Si3N4)
- Thermal diffusion profiles
- Ion-implanted layers with depth-dependent damage

Parameters:

$n_func - coderef sub($lambda, $z_frac) returning complex N at fractional depth z (0 = top, 1 = bottom)
$d_total - total graded layer thickness [nm]
n_slices - number of sublayers (default 20; increase for accuracy)

Example:

use Physics::Ellipsometry::VASE::TMM qw(tmm_graded);
use Physics::Ellipsometry::VASE::Dispersion qw(cauchy_nk);

# Linear gradient from SiO2 (n=1.46) to Si3N4 (n=2.0)
my $grad_func = sub {
    my ($lambda, $z) = @_;
    my $n = 1.46 + 0.54 * $z;  # linear gradient
    return pdl($n) + i() * pdl(0);
};

my $lambda = sequence(100) * 10 + 300;
my ($rp, $rs) = tmm_graded($lambda, pdl(70),
    pdl(1.0)+i()*pdl(0),       # air ambient
    pdl(3.94)+i()*pdl(0.02),   # Si substrate
    $grad_func, 200,           # 200 nm graded layer
    n_slices => 30,
);

tmm_interface

my ($rp, $rs) = tmm_interface($lambda_nm, $theta_deg,
                               $N_ambient, $N_substrate,
                               $N_a, $N_b, $d_interface,
                               n_slices => 10,
                               ema_model => 'bruggeman');

Models a thin interface layer (interlayer) between two materials using an EMA-graded profile. The composition varies linearly from 100% material A at the top to 100% material B at the bottom.

Physical motivation: Real interfaces are rarely atomically sharp. Interdiffusion, roughness, and reaction products create a transition region that is better modeled as a graded interlayer than as an abrupt boundary.