N Layer Fresnel Equations for TM Polarization

The Fresnel reflection for an arbitrary number of layers in the case of TM polarization is

\[\begin{align} \tilde{r}_{j,l,m \ldots n} = \frac{\tilde{r}_{j,l} + \tilde{r}_{l,m \ldots n} \exp\left(2 \mathrm{i} k_{z,l} d_l\right)} {1+\tilde{r}_{j,l} \tilde{r}_{l,m \ldots n} \exp\left(2 \mathrm{i} k_{z,l} d_l\right)} \end{align}\]

where \(d_l\) is the thickness of the \(l\)th metal layer and

\[\begin{align}\tilde{r}_{i,j}&= \left.\left(\frac{k_{z,i}}{\epsilon_i}-\frac{k_{z,j}}{\epsilon_j}\right) \middle/ \left(\frac{k_{z,i}}{\epsilon_i}+\frac{k_{z,j}}{\epsilon_j}\right)\right.\\ &=\frac{\epsilon_j k_{z,i}-\epsilon_i k_{z,j}} {\epsilon_j k_{z,i}+\epsilon_i k_{z,j}}\end{align}\] is the two layer Fresnel relation.

This equation can be iterated to produce the correct reflectivity for any number of layers, as the following Python snippet shows:

# n layer Fresnel
def r(kx,epsilon,d):
    # two layer Fresnel
    rtwo = lambda kx,epsilon: \
        (epsilon[1]*kz(kx,epsilon[0])-epsilon[0]*kz(kx,epsilon[1]))/ \
        (epsilon[1]*kz(kx,epsilon[0])+epsilon[0]*kz(kx,epsilon[1]))

  if len(epsilon) == 2:
    return rtwo(kx,epsilon)
  else:
    return\
    (rtwo(kx,epsilon)+r(kx,delete(epsilon,0),delete(d,0))*exp(2.0j*kz(kx,epsilon[1])*d[1]))/\
    (1+rtwo(kx,epsilon)*r(kx,delete(epsilon,0),delete(d,0))*exp(2.0j*kz(kx,epsilon[1])*d[1]))

or in octave

% returns the complex n layer Fresnel reflectivity for TM polarization
function out = nlayerfresnel(k0,kx,epsilon,d)
    kz = @(kx,epsilon) sqrt(k0.^2.*epsilon-kx.^2);
    % two layer Fresnel
    r12 = @(kx,epsilon) \
                (epsilon(2).*kz(kx,epsilon(1))-epsilon(1).*kz(kx,epsilon(2)))./ \
                (epsilon(2).*kz(kx,epsilon(1))+epsilon(1).*kz(kx,epsilon(2)));

    % if the length is two, return the two layer Fresnel, otherwise
    % recursively go through the layers
    if length(epsilon) == 2
        out = r12(kx,epsilon);
    else
        out = (r12(kx,epsilon)+nlayerfresnel(k0,kx,epsilon(2:end),d(2:end)).*exp(2.0i.*kz(kx,epsilon(2)).*d(2)))./\
                    (1+r12(kx,epsilon).*nlayerfresnel(k0,kx,epsilon(2:end),d(2:end)).*exp(2.0i.*kz(kx,epsilon(2)).*d(2)));
    endif
endfunction

Note here that epsilon and d are ordered arrays holding the ordered list of complex permittivities of the system you’re trying to model.