How to calculate thin film interference when the incident medium is absorbing?

Question

Edit: I found a text that mentions that when the incident medium is absorbing, the usual equations for interference don't apply, so I have refined my question in the title.

---

I have been working on creating a Blender simulation for thin-film interference using the complex versions of the Fresnel and interference equations.

Part of the implementation of the simulation includes reversing the order of the layers when we see backfacing geometry, such that the substrate becomes the medium and the medium becomes the substrate in the calculations. The reasoning is that light also exits through the (transparent) material. That's where the problem comes in.

These reversed layers do calculate correctly for non-absorbing substrates, but the moment a non-zero extinction coefficient is entered, it starts spitting out values greater than 1 for both reflectance and transmittance.

In the first implementation of the equations, I used

$$ \hat a_x = \hat n_x cos(\theta_x) = \sqrt{ \hat n_x^2 - \hat n_1^2 sin(\theta_1) } \ , \; \hat n = n + ik$$

within the Fresnel equations

$$ \begin{align} \hat r_{s|ij} & = \frac{ \hat a_i - \hat a_j }{ \hat a_i + \hat a_j } \\ \hat r_{p|ij} & = \frac{ \hat n_j^2 \ \hat a_i - \hat n_i^2 \ \hat a_j }{ \hat n_j^2 \ \hat a_i + \hat n_i^2 \ \hat a_j } \\ \hat t_{s|ij} & = \frac{ 2 \hat a_i }{ \hat a_i + \hat a_j } \\ \hat t_{p|ij} & = \frac{ 2 \hat n_i \hat n_j \hat a_i }{ \hat n_j^2 \ \hat a_i + \hat n_i^2 \ \hat a_j } \\ \end{align} $$

and within the interference equations

$$ \begin{align} \hat r & = \frac{\hat r_{12} + \hat r_{23} e^{i 2 \beta}}{1 + \hat r_{12} \hat r_{23} e^{i 2 \beta}} \\ \hat t & = \frac{\hat t_{12} \hat t_{23} e^{i \beta}}{1 + \hat r_{12} \hat r_{23} e^{i 2 \beta}} \\ \beta & = \frac {2 \pi t}{\lambda} \hat a_{film} \end{align} $$

Taking the calculations to Excel, I found that something about these equations were causing values for reflectance and transmittance to be greater than 1 (I even saw values of 14 for reflectance!). These values generally increased with the thickness.

After more investigation, I found that when the thickness is zero (thereby eliminating the film), the flipped layers do calculate correctly. This lead me to believe that the OPD $(\beta)$ used in the interference calculations is the problem.

---

Looking around for possible solutions, for which there seemed to be no content, I found this spectral reflectance calculator which showed that there was a way to calculate it correctly. There's no way to actually enter complex refractive indices into the calculator, so I simply used the provided metals to test the flipped layer idea. I also used it to corroborate the results of my implementation.

Wanting to test another method of interference calculations, I also created an implementation that uses the 2x2 matrix method, but the results were the same.

More searching brought me to Born and Wolf's Principles of Optics (7th Edition) where they derive calculations for absorbing films on a transparent surface. Using their derivation of the Fresnel and interference equations for that specific situation as a guide for my purposes, the following was found:

$$ \hat n_x cos(\theta_x) = u_x + iv_x $$

where

$$ \begin{align} 2u_x^2 & = a_x + \sqrt{ a_x^2 + 4 b_x^2 } \\ 2v_x^2 & = -a_x + \sqrt{ a_x^2 + 4 b_x^2 } \end{align} $$

and

$$ \begin{align} a_x & = n_x^2 - k_x^2 - (n_1^2 - k_1^2) sin(\theta_1) \\ b_x & = n_x k_x - n_1 k_1 sin(\theta_1) \end{align} $$

Using $ u_x + iv_x $ instead of $ \hat a_x $ in the Fresnel equations and the OPD, the reflectances and transmittances do not reach as high as before. In fact, I have yet to find values that reach above 1. Although, when checking for energy conservation by adding the values together $ (R + T) $, the result is greater than 1 at and close to normal incidence.

Comparing the usage of $ u_x + iv_x $ to $ \hat a_x $ showed that the problem lies in how the imaginary part is calculated; switching back and forth between the real parts showed no change.

To confirm whether the OPD was the problem, I used the original Fresnel equations with the new OPD equation, and I also used the old OPD with the new Fresnel equations. This showed exactly the same results as what resulted from the new equations and old equations respectively, suggesting the problem does lie in the OPD. (Edit: it seems that the OPD may not truly be the problem as described by the note I added at the top of this post)

Also, the reflectance calculator seems to suggest that, unless the film is also absorbing, there is no absorbance when light passes out of an absorbing substrate into the medium. I also found another calculator, but it gave exactly the same results as my new calculations.

How do I go about performing interference calculations for thin-films with transparent, slightly absorbing substrates when light is exiting the material?

---

The following are some of the results from the simulation.

Results of the equations using $ \hat a $: (Left: $R+T$ Color of the material, Right: Full material)

$$ \hat n_{sub} = 1.5 + i 0.1, \ \hat n_{film} = 2 + 0i, \ t = 500nm $$

Results of the equations using $ u + i v $: (Left: $R+T$ Color of the material, Right: Full material)

$$ \hat n_{sub} = 1.5 + i 0.1, \ \hat n_{film} = 2 + 0i, \ t = 500nm $$

The following show whether $R+T$ is greater than one. If true, it returns 1 for whichever of the R, G, B values is greater than 1 to make a color. Here's a link to animations of these values changing as thickness increases. (Left: Equations using $ \hat a $, Right: Equations using $ u + i v $)

---

Edit: I realized that $\hat a$ can give the same results as using $u + iv$ because $u + iv$ is simply the other solution to the square root.

I also was able to implement another version of the interference equation where each Fresnel coefficient has been converted to it's exponential form:

$$ \hat r_{ij} = |\hat r_{ij}| e^{i(arg(\hat r_{ij}))} = \rho_{ij} e^{i \phi_{ij}}$$ $$ \hat t_{ij} = |\hat t_{ij}| e^{i(arg(\hat t_{ij}))} = \tau_{ij} e^{i \chi_{ij}}$$

such that the new equations become:

$$ \begin{align} \hat r & = \frac{\rho_{12} e^{i \phi_{12}} + \rho_{23} e^{i (2 \beta + \phi_{23})}}{1 + \rho_{12} \rho_{23} e^{i (2 \beta + \phi_{12} + \phi_{23})}} \\ \hat t & = \frac{\tau_{12} \tau_{23} e^{i (\beta + \chi_{12} + \chi_{23})}}{1 + \rho_{12} \rho_{23} e^{i (2 \beta + \phi_{12} + \phi_{23})}} \end{align} $$

These equations unfortunately give the same results as in the second pair of images above. I guess it's to be expected when it's simply just another form of the same thing.

Answer 1

After more research, it seems I was able to find an answer to my question.

The following were adapted from "Principles of Optics" (7th Edition) by Born and Wolf:

Let's define a wave vector:

$$\kappa = \hat n_t \omega (k_x, \ k_y, \ k_z)$$ $$k_x=sin \theta_t, \ k_y=0, \ \text {and} \ k_z=cos \theta_t$$

According to Snell's law:

\begin{align} sin \theta_t&=\frac{\hat n_i}{\hat n_t} sin \theta_i \\ cos \theta_t&=\sqrt{1-\frac{\hat n_i}{\hat n_t} sin \theta_i} \end{align}

Substituting these into the wave vector:

$$\kappa=\omega \left(\hat n_isin\theta_i, \ 0, \ \hat n_t\sqrt{1-\frac{\hat n_i}{\hat n_t} sin \theta_i} \right)$$

Let:

$$u_t+iv_t=\hat n_t \ cos \theta_t=\hat n_t \sqrt{1-\frac{\hat n_i}{\hat n_t} sin \theta_i}$$

Squaring both sides and equating real and imaginary parts, where $\hat n=n+ik$, results in:

\begin{align} u^2_t-v^2_t&=n_t^2-k_t^2-(n_i^2-k_i^2)sin\theta_i=b_t \\ u_tv_t&=n_tk_t-n_ik_isin\theta_i=a_t \end{align}

From which:

\begin{align} 2u_t^2 & = a_t + \sqrt{ a_t^2 + 4 b_t^2 } \\ 2v_t^2 & = -a_t + \sqrt{ a_t^2 + 4 b_t^2 } \end{align}

Now:

$$\kappa=\omega \left((n_i+ik_i)sin\theta_i, \ 0, \ u_t+iv_t \right)$$

Using trigonometry, we can now find the angles that the real and imaginary phases make with the boundary, respectively:

\begin{align} cos\theta_{re}&=\frac {u_t}{\sqrt{n_i^2sin^2\theta_i+u_t^2}} \\ cos\theta_{im}&=\frac {v_t}{\sqrt{k_i^2sin^2\theta_i+v_t^2}} \end{align}

---

Moving to equations adapted from "Thin-Film Optical Filters" (4th Edition) by H. Angus Macleod:

The optical admittance of a material for $s$ and $p$ polarization are respectively given as:

\begin{align} \eta_s&=ncos\theta_{re}+ikcos\theta_{im} \\ \eta_s&=\frac {\hat n^2}{\eta_s} \end{align}

Really, the optical admittance of free space should be included, but because it would later cancel out, it is left out for simplicity.

Derived from the boundary equations, we can describe a single layer's matrix as:

$$ M_x= \begin{bmatrix} cos\delta_x&-isin\delta_x/\eta_x\\ -i\eta_xsin\delta_x&cos\delta_x \end{bmatrix} $$

Where $\delta=(2\pi d/\lambda) (ncos\theta_{re}+ikcos\theta_{im})$, $d=$ layer thickness, and $\lambda=$ wavelength

The characteristic matrix of a stack with $x$ amount of layers and substrate $m$ is:

$$ \begin{bmatrix} B \\ C \end{bmatrix} = M_1M_2 \cdots M_x \begin{bmatrix} 1 \\ \eta_m \end{bmatrix} $$

From here, we can finally calculate the reflectance and transmittance:

\begin{align} R_{s,p}&=\left| \frac {B\eta_0-C}{B\eta_0^*+C}\right|^2 \\ T_{s,p}&=\frac {4Re(\eta_0)Re(\eta_m)}{\left|B\eta_0^*+C\right|^2} \end{align}

For unpolarized light:

\begin{align} R&=\frac{1}{2}(R_s+R_p) \\ T&=\frac{1}{2}(T_s+T_p) \end{align}

---

The two biggest things that helped solve the energy conservation problem for absorbing incident media were finding the new angles, $\theta_{re}$ and $\theta_{im}$, and the appearance of the complex conjugate of $\eta_0$ within the denominator of the formulas for reflectance and transmittance.

In "Principles of Optics", $\theta_{re}$ was used as the transmitted angle in the Fresnel equations to describe conductors. Naturally, I used that value within the optical admittance $\eta_s$ (and by extension $\eta_p$) which was described by $\eta_s=\hat n cos\theta$ in "Thin-Film Optical Filters". While $R+T=1$ even in absorbing media, the materials I had created before the change were altered vastly despite not being calculated with absorbing incident media both before and after the change.

After discussing the change, it was noted that, of course,

if you use a more general law for a case the simpler law was valid for, you must get the same outcome.

After rereading the article, Ray tracing in absorbing media, I noticed that it had a separate angle for the complex part. That lead to the derivation of $\eta_s=ncos\theta_{re}+ikcos\theta_{im}$. This correction lead to the reversion of the changed materials.

Although, this did not yet allow $R+T=1$. Through Macleod's article, The Mixed Poynting Vector, I found another definition of transmission:

$$T=\frac {4 Re(Y)\left|\eta_0\right|^2}{Re(\eta_0)\left[(Re(\eta_0)+Re(Y))^2+(Im(\eta_0)-Im(Y))^2\right]}$$

Substituting $Y=C/B$ (where $Y$ is the optical admittance of the stack) a relation described in "Thin-Film Optical Filters" produces:

$$T=\frac {4 \left|\eta_0\right|^2 Re(BC^*)}{Re(\eta_0)\left|B\eta_0^* + C\right|^2}$$

$Re(BC^*)=Re(\eta_m)$ according to the book which would mean the differences with the original transmittance formula are $\left|\eta_0\right|^2/Re(\eta_0)$ and $\eta_0^*$ (the formula for transmittance in the book is $(4\eta_0Re(\eta_m))/|B\eta_0+C|^2$ where the incident medium is assummed to be real).

With the differences taken into account, reflectance and transmittance still did not equal 1. In fact, due to $\left|\eta_0\right|^2$, $T$ would keep growing as the extinction coefficient increased. Although, $\eta_0^*$ did remind me of another article I came across previously, Rigorous expressions for the Fresnel equations at interfaces between absorbing media. The article derives the following for $s$-polarized light:

$$r_s=\frac {k_{z|i}-k_{z|t}}{k_{z|i}^*+k_{z|t}}$$

Because the optical admittance is linearly dependent to $k_z$, I could use it to derive the following:

$$r_s=\frac {\eta_{s|0}-Y_s}{\eta_{s|0}^*+Y_s}$$

Because we will be taking the modulus of $r_s$ above to derive transmittance, the same formula can be used for $p$-polarization (as a result of the way that polarization is handled):

$$T=1-|r|^2=\frac {4Re(\eta_0)Re(BC^*)}{\left|B\eta_0^*+C\right|^2}$$

Again, $Re(BC^*)=Re(\eta_m)$ which means that the equation is nearly the same as the formula described in "Thin-Film Optical Filters" (and reflectance also takes on the conjugate value).

Finally, $R+T=1$ in absorbing media.

While I don't know if the conjugate value was intended, because $|a+ib|=|a-ib|$ when $a$ and $b$ are real but not all the time when they are complex (and there didn't seem to be any indication of the knowledge), it still really helped solve the problem.

If you would like to check out the results of the implementation and maybe play with it yourself, check out this Blender Artists forum thread!

Here's a sample of some materials that can be made with it: