\(\DeclareMathOperator{\Div}{div}
\DeclareMathOperator{\Rot}{rot}
\newcommand{\parder}[2]{\frac{\partial {#1}}{\partial {#2}}}\)$\setCounter{0}$
This post should be a straightforward walkthrough of the Berreman’s paper1 using the approach we implemented in the previous posts, that is: minimizing the amount of constants by working in the HL system of units and using the geometrical components of the wave-vectors. It should be stated that the original paper is already very clear in this respect. In this post we will follow it up to the point where \(\mathbf{\hat{S}}\) matrix and \(\boldsymbol{\hat{\Delta}}\) matrix are defined – and remove one unnecessary minus.
In following blog post we will treat the interface and discuss some practical questions of wave propagation in quadrochromic thin films and multilayers.
Maxwell equations in 6×6 form
We start from scratch with Maxwell’s equations in the form as seen before:
\(\begin{align}
\begin{aligned}
\Div\mathbf{D} &=0, &\Div\mathbf{B} &=0,
\\
-\Rot\mathbf{E} &=\frac{1}{c}\parder{\mathbf{B}}{t}, &\Rot\mathbf{H} &=\frac{1}{c}\parder{\mathbf{D}}{t} .
\end{aligned}
\end{align}\)
So only the curly part of the electro-magnetic field is of interest, and the corresponding two equations can be written explicitly in a 6×6 matrix form:
\(
\begin{align}
\begin{bmatrix}
0 & 0 & 0 & 0 & -\parder{}{z} & \parder{}{y} \\
0 & 0 & 0 & \parder{}{z} & 0 & -\parder{}{x} \\
0 & 0 & 0 & -\parder{}{y} & \parder{}{x} & 0 \\
0 & \parder{}{z} & -\parder{}{y} & 0 & 0 & 0 \\
-\parder{}{z} & 0 & \parder{}{x} & 0 & 0 & 0 \\
\parder{}{y} & -\parder{}{x} & 0 & 0 & 0 & 0
\end{bmatrix}
\begin{bmatrix}
E_x \\ E_y \\ E_z \\ H_x \\ H_y \\ H_z
\end{bmatrix}
=
\frac{1}{c}\parder{}{t}
\begin{bmatrix}
D_x \\ D_y \\ D_z \\ B_x \\ B_y \\ B_z
\end{bmatrix}.
\end{align}
\)
We will use the same abbreviation as Berreman:
\(\begin{align}\label{eqRGC}
\mathbf{\hat{R}}\mathbf{G}=\frac{1}{c}\parder{}{t}\mathbf{C},
\end{align}\)
where \(\mathbf{\hat{R}}\) is the rotation differential operator, \(\bf{G}\) is the vector of electric and magnetic fields and \(\bf{C}\) is the vector of the electric displacement and magnetic induction. The constitutive relations will have a general form:
\(\begin{align}
\begin{aligned}
\mathbf{D} &=\boldsymbol{\hat{\varepsilon}}\mathbf{E} + \boldsymbol{\hat{\rho}}\mathbf{H},
\\
\mathbf{B} &=\boldsymbol{\hat{\mu}}\mathbf{H}+ \boldsymbol{\hat{\rho}’}\mathbf{E}.
\end{aligned}
\end{align}\)
explicitly as
\(
\begin{align}
\begin{bmatrix}
D_x \\ D_y \\ D_z \\ B_x \\ B_y \\ B_z
\end{bmatrix}
=
\begin{bmatrix}
\varepsilon_{11} & \varepsilon_{12} & \varepsilon_{13} & \rho_{11} & \rho_{12} & \rho_{13} \\
\varepsilon_{21} & \varepsilon_{22} & \varepsilon_{23} & \rho_{21} & \rho_{22} & \rho_{23} \\
\varepsilon_{31} & \varepsilon_{32} & \varepsilon_{33} & \rho_{31} & \rho_{32} & \rho_{33} \\
\rho_{11}’ & \rho_{12}’ & \rho_{13}’ & \mu_{11} & \mu_{12} & \mu_{13} \\
\rho_{21}’ & \rho_{22}’ & \rho_{23}’ & \mu_{21} & \mu_{22} & \mu_{23} \\
\rho_{31}’ & \rho_{32}’ & \rho_{33}’ & \mu_{31} & \mu_{32} & \mu_{33}
\end{bmatrix}
\begin{bmatrix}
E_x \\ E_y \\ E_z \\ H_x \\ H_y \\ H_z
\end{bmatrix},
\end{align}
\)
abbreviated as
\(\begin{align}\label{eqCMG}
\mathbf{C}=\mathbf{\hat{M}}\mathbf{G},
\end{align}\)
with \(\mathbf{\hat{M}}\) being the material’s linear response tensor, written in a coordinate system defined below.
Stationary solution
We will immediatelly assume homogeneous bulk and search for a solution of \eqref{eqRGC} with \eqref{eqCMG} in the form of plane waves
\(\begin{align}\label{eqKaGamma}
\mathbf{G}=\boldsymbol{\Gamma} e^{i(\mathbf{kr}-\omega t)} ,
\end{align}\)
where \(\boldsymbol{\Gamma}\) is the vector of field amplitudes and we will want to choose the coordinate system so the wave \(\mathbf{k}\) propagates in the x-z plane as we are used to. We also introduce the geometrical components:
\(\begin{align}\label{eqKaxiq}
\mathbf{k}=\frac{\omega}{c}\boldsymbol{\kappa} = \frac{\omega}{c}(\xi,0,q),
\end{align}\)
where \(\xi\) is the conserved real x-component of the wave-vector and we will need to find possible \(q\), the z-components, just like we did in the Part 2 of the isotropic problem, only here we will keep Berreman’s notation \(q\) for the z-components. There will be 4 solutions and for each of them we will find the corresponding amplitude/polarization vector.
The following figure illustrates the situation. Althought, at this point, we simply look for the modes possible in the bulk, we will later couple them through interfaces to the ambient. In the transparent ambient we are able to set and measure the real angle of incidence, which fixes the invariant \(\xi\). The angle of refraction loses clear sense in an absorbing medium, as it would need to be expressed as complex quantity. That is due to the fact that the wave vector becomes complex, but its imaginary parts have only z-component – as the x-component \(\xi\) is real and conserved. Moreover, the angle of refraction is hardly ever needed for understanding the ellipsometeric experiment.
The problem we will be solving is an eigen-problem, and the four possible z-components of the wave-vector are eigenvalues of some operator – to which we will arrive at the end of this post. The eigenvectors associated with the eigenvalues describe the polarization states that propagate through the medium without change, apart from evolving phase and decreasing amplitude. For example, linearly polarized light in solution of some sugars will experience rotation of the polarization direction – so it is not the polarization eigen-state of such medium. But there will be a circular polarization that propagates and stays circular, experiencing certain refractive index, while the opposite circular polarization also stays in its state, just experiencing different refractive index. That is called the circular birefringence.
The Berreman’s method is an expansion of the simple isotropic case, where we found that for both linear polarizations, the z-component of the wave-vector was \(q=\pm\sqrt{\varepsilon\mu-\xi^2}\), with the right sign for the down- and up-propagating mode. The material response tensor \(\mathbf{\hat{M}}\) has some physical constrains, but can lead to situations where the four modes (two down and two up) propagating in the medium have distinct properties and carry generally elliptical polarization states. Such general media are called quadrochromic in the litarature.
On the left-hand side of equation \eqref{eqRGC} we apply the differential operator \(\mathbf{\hat{R}}\) on the vector \(\mathbf{G}\) in the form of equation \eqref{eqKaGamma}, with the wave-vector \eqref{eqKaxiq}, obtaining:
\(
\begin{align}
\mathbf{\hat{R}}\mathbf{G}=
i\frac{\omega}{c}
\begin{bmatrix}
0 & 0 & 0 & 0 & -q & 0 \\
0 & 0 & 0 & q & 0 & -\xi \\
0 & 0 & 0 & 0 & \xi & 0 \\
0 & q & 0 & 0 & 0 & 0 \\
-q & 0 & \xi & 0 & 0 & 0 \\
0 & -\xi & 0 & 0 & 0 & 0
\end{bmatrix}
\boldsymbol{\Gamma} e^{i(\mathbf{kr}-\omega t)}.
\end{align}
\)
On the right-hand side of \eqref{eqRGC} we apply the time derivative on the displacement-induction vector \(\mathbf{C}\), while the response tensor \(\mathbf{\hat{M}}\) is constant.
\(\begin{align}
\frac{1}{c}\parder{}{t}\mathbf{C} = \frac{1}{c}\mathbf{\hat{M}}\parder{}{t}\mathbf{G} = -i\frac{\omega}{c}\mathbf{\hat{M}}\boldsymbol{\Gamma} e^{i(\mathbf{kr}-\omega t)}.
\end{align}\)
Now we combine the last two expressions, cancel the constants and the non-zero complex exponential factor. In that we have converted the differential physical problem to a dimensionless, algebraic or geometrical, one.
\(
\begin{align}\label{eqGMG}
\begin{bmatrix}
0 & 0 & 0 & 0 & -q & 0 \\
0 & 0 & 0 & q & 0 & -\xi \\
0 & 0 & 0 & 0 & \xi & 0 \\
0 & q & 0 & 0 & 0 & 0 \\
-q & 0 & \xi & 0 & 0 & 0 \\
0 & -\xi & 0 & 0 & 0 & 0
\end{bmatrix}
\boldsymbol{\Gamma}
=
-\mathbf{\hat{M}}\boldsymbol{\Gamma}.
\end{align}
\)
In this equation, the unknowns are \(q\) – the geometric z-component of the wave-vector – and the elements of \(\boldsymbol{\Gamma}\), the electric and magnetic field amplitudes, which both have the same unit in HLU. It is immediatelly apparent that the fields vector will be found up to a constant prefactor. It describes a polarization state, so multiplying it by constant does not change the state. What follows now is simply commented algebra.
Reduction to 4×4 form
First we extract the 3rd and 6th rows as they do not contain \(q\), and have on the right hand side the \(D_z\) and \(B_z\) components. As before, we will be interested in the state of the in-plane, xy-components of the fields, even we did not yet defined any interface. Taking apart the two rows allows us to separate the \(E_z\) and \(H_z\) components.
So the third and sixth rows yield:
\(\begin{align}
\begin{aligned}
\xi\Gamma_5 &= -\sum_{j=1}^{6} M_{3j}\Gamma_j,
\\
\xi\Gamma_2 &= \sum_{j=1}^{6} M_{6j}\Gamma_j.
\end{aligned}
\end{align}\)
We write this explicitly, separate the \(\Gamma_3\) and \(\Gamma_6\) and define shorthand symbols \(A_3\) and \(A_6\).
\(\begin{align}
\begin{aligned}
\overbrace{ -M_{31}\Gamma_1 – M_{32}\Gamma_2 – M_{34}\Gamma_4 – (M_{35}+ \xi)\Gamma_5 }^{A_3} &= M_{33}\Gamma_3 + M_{36}\Gamma_6 ,
\\
\overbrace{ -M_{61}\Gamma_1 – (M_{62}- \xi)\Gamma_2 – M_{64}\Gamma_4 – M_{65}\Gamma_5 }^{A_6} &= M_{63}\Gamma_3 + M_{66}\Gamma_6 .
\end{aligned}
\end{align}\)
Then
\(\begin{align}\label{eqDMM}
D \equiv M_{33}M_{66} – M_{36}M_{63},
\end{align}\)
and the solution for \(\Gamma_3 \equiv E_z\) and \(\Gamma_6 \equiv H_z\) is
\(\begin{align}
\Gamma_3 = \frac{M_{66}A_{3} – M_{36}A_{6}}{D} , \quad \Gamma_6 = \frac{M_{33}A_{6} – M_{63}A_{3}}{D}.
\end{align}\)
We will rewrite the solutions with new shorthand symbols \(a_{mn}\) as:
\(\begin{align}
\begin{aligned}
\Gamma_3 &= a_{31}\Gamma_{1} + a_{32}\Gamma_{2} + a_{34}\Gamma_{4} + a_{35}\Gamma_{5} ,
\\
\Gamma_6 &= a_{61}\Gamma_{1} + a_{62}\Gamma_{2} + a_{64}\Gamma_{4} + a_{65}\Gamma_{5} ,
\end{aligned}
\label{eqG3G6}
\end{align}\)
where the \(a_{mn}\) are just various combinations of the material response tensor elements, some modified by \(\xi\).
\(\begin{align}
\begin{aligned}
a_{31} &= \frac{1}{D} ( -M_{31}M_{66} + M_{61}M_{36} ), \\
a_{32} &= \frac{1}{D} ( -M_{32}M_{66} + (M_{62}-\xi)M_{36} ), \\
a_{34} &= \frac{1}{D} ( -M_{34}M_{66} + M_{64}M_{36} ), \\
a_{35} &= \frac{1}{D} ( -(M_{35}+\xi)M_{66} + M_{65}M_{36} ), \\
a_{61} &= \frac{1}{D} ( -M_{61}M_{33} + M_{31}M_{63} ), \\
a_{62} &= \frac{1}{D} ( -(M_{62}-\xi)M_{33} + M_{32}M_{63} ), \\
a_{64} &= \frac{1}{D} ( -M_{64}M_{33} + M_{34}M_{63} ), \\
a_{65} &= \frac{1}{D} ( -M_{65}M_{33} + (M_{35}+\xi)M_{63} ).
\end{aligned}
\end{align}\)
Note that some important role is played by the \(D\), defined above in \eqref{eqDMM} as a combination of the zz-elements of the tensors \(\boldsymbol{\hat{\varepsilon}}\), \(\boldsymbol{\hat{\mu}}\), \(\boldsymbol{\hat{\rho}}\) and \(\boldsymbol{\hat{\rho}’}\). We will have to assume that \(D\neq 0\), which – in most common cases – means non-zero \(\varepsilon_{33}\). Curious effects can develop when the \(\varepsilon_{33}\), or \(D\), are near-zero. Those are called Berreman’s effects, but in reference to another Berreman’s cornerstone paper.2 We shall discuss them separately later.
Now we can write the remaining four equations from eq. \eqref{eqGMG} as
\(
\begin{align}
q
\begin{bmatrix}
0 & 0 & 0 & -1 \\
0 & 0 & 1 & 0 \\
0 & 1 & 0 & 0 \\
-1 & 0 & 0 & 0
\end{bmatrix}
\begin{bmatrix}
\Gamma_1 \\ \Gamma_2 \\ \Gamma_4 \\ \Gamma_5
\end{bmatrix}
=-
\begin{bmatrix}
\sum M_{1j}\Gamma_j \\
\sum M_{2j}\Gamma_j – \xi\Gamma_6 \\
\sum M_{4j}\Gamma_j \\
\sum M_{5j}\Gamma_j – \xi\Gamma_3
\end{bmatrix}
\equiv
-\mathbf{\hat{S}}
\begin{bmatrix}
\Gamma_1 \\ \Gamma_2 \\ \Gamma_4 \\ \Gamma_5
\end{bmatrix},
\label{eqqSGamma}
\end{align}\)
where we inserted the \(\Gamma_3\) and \(\Gamma_6\) from \eqref{eqG3G6} and defined the \(\mathbf{\hat{S}}\) matrix identical to Berreman’s manuscript:
\( \small{
\mathbf{\hat{S}}
=
\begin{bmatrix}
M_{11} + M_{13}a_{31} + M_{16}a_{61} & M_{12} + M_{13}a_{32} + M_{16}a_{62} & M_{14} + M_{13}a_{34} + M_{16}a_{64} & M_{15} + M_{13}a_{35} + M_{16}a_{65} \\
M_{21} + M_{23}a_{31} + (M_{26}-\xi)a_{61} & M_{22} + M_{23}a_{32} + (M_{26}-\xi)a_{62} & M_{24} + M_{23}a_{34} + (M_{26}-\xi)a_{64} & M_{25} + M_{23}a_{35} + (M_{26}-\xi)a_{65} \\
M_{41} + M_{43}a_{31} + M_{46}a_{61} & M_{42} + M_{43}a_{32} + M_{46}a_{62} & M_{44} + M_{43}a_{34} + M_{46}a_{64} & M_{45} + M_{43}a_{35} + M_{46}a_{65} \\
M_{51} + (M_{53}+\xi)a_{31} + M_{56}a_{61} & M_{52} + (M_{53}+\xi)a_{32} + M_{56}a_{62} & M_{54} + (M_{53}+\xi)a_{34} + M_{56}a_{64} & M_{55} + (M_{53}+\xi)a_{35} + M_{56}a_{65}
\end{bmatrix}.
} \)
Reorganize to get Delta form
Now we just reorganize the rows and colums, and extinguish minuses. In the field components, the equation \eqref{eqqSGamma} reads
\(
\begin{align}
q
\begin{bmatrix}
-H_y \\ H_x \\ E_y \\ -E_x
\end{bmatrix}
=
-\mathbf{\hat{S}}
\begin{bmatrix}
E_x \\ E_y \\ H_x \\ H_y
\end{bmatrix}.
\end{align}\)
We will interchange the first, second and last row of the whole system, then swap second, third and last column of the matrix on the right hand side. Finally, we switch some minuses obtaining
\(
\begin{align}
q
\begin{bmatrix}
E_x \\ H_y \\ E_y \\ H_x
\end{bmatrix}
=
\overbrace{
\begin{bmatrix}
S_{41} & S_{44} & S_{42} & S_{43} \\
S_{11} & S_{14} & S_{12} & S_{13} \\
-S_{31} & -S_{34} & -S_{32} & -S_{33} \\
-S_{21} & -S_{24} & -S_{22} & -S_{23}
\end{bmatrix}}^{\boldsymbol{\hat{\Delta}}}
\overbrace{
\begin{bmatrix}
E_x \\ H_y \\ E_y \\ H_x
\end{bmatrix}}^{\boldsymbol{\Psi}},
\end{align}\)
which defines the \(\boldsymbol{\hat{\Delta}}\) matrix and fixes the order of polarization state vector \(\boldsymbol{\Psi}\). Note that the order of the components of \(\boldsymbol{\Psi}\) is rather arbitrary, and is chosen so the components of linear p- and s- polarizations appear together in sub-blocks.
Berreman introduced additional minus in front of the \(H_x\) component, which is related to s-polarization and served merely esthetical purposes: In the right-handed media and wave propagating along z-axis, the p-polarization with positive \(E_x\) amplitude will have positive \(H_y\), but s-polarization with positive \(E_y\) will have negative \(H_x\). Placing the minus there will therefore obliterate some minuses when later solving simple problems in semi-infinite media. Schubert3 defined the \(\boldsymbol{\Psi}\) vector without the decorative minus, but in different order of components.
We have now our eigenproblem:
| \(\begin{align} q\boldsymbol{\Psi} = \boldsymbol{\hat{\Delta}}\boldsymbol{\Psi}, \end{align}\) |
with \(\boldsymbol{\hat{\Delta}}\) matrix explicitly:
\( \small{
\boldsymbol{\hat{\Delta}}
=
\begin{bmatrix}
M_{51} +(M_{53}+\xi)a_{31} +M_{56}a_{61} & M_{55} +(M_{53}+\xi)a_{35} +M_{56}a_{65} & M_{52} +(M_{53}+\xi)a_{32} +M_{56}a_{62} & M_{54} +(M_{53}+\xi)a_{34} +M_{56}a_{64} \\
M_{11} +M_{13}a_{31} +M_{16}a_{61} & M_{15} +M_{13}a_{35} +M_{16}a_{65} & M_{12} +M_{13}a_{32} +M_{16}a_{62} & M_{14} +M_{13}a_{34} +M_{16}a_{64} \\
-M_{41} -M_{43}a_{31} -M_{46}a_{61} & -M_{45} -M_{43}a_{35} -M_{46}a_{65} & -M_{42} -M_{43}a_{32} -M_{46}a_{62} & -M_{44} -M_{43}a_{34} -M_{46}a_{64} \\
-M_{21} -M_{23}a_{31} -(M_{26}-\xi)a_{61} & -M_{25} -M_{23}a_{35} -(M_{26}-\xi)a_{65} & -M_{22} -M_{23}a_{32} -(M_{26}-\xi)a_{62} & -M_{24} -M_{23}a_{34} -(M_{26}-\xi)a_{64}
\end{bmatrix}.
} \)
Inserting the response tensor
And since internet is realy huge, we can occupy some more space by expressing everything in the original response tensors:
\(\begin{align}
D = \varepsilon_{33}\mu_{33} – \rho_{33}\rho_{33}’,
\end{align}\)
\(\begin{align}
\begin{aligned}
a_{31} &= \frac{1}{D} ( -\varepsilon_{31}\mu_{33} + \rho_{31}’\rho_{33} ), \\
a_{32} &= \frac{1}{D} ( -\varepsilon_{32}\mu_{33} + (\rho_{32}’-\xi)\rho_{33} ), \\
a_{34} &= \frac{1}{D} ( -\rho_{31}\mu_{33} + \mu_{31}\rho_{33} ), \\
a_{35} &= \frac{1}{D} ( -(\rho_{32}+\xi)\mu_{33} + \mu_{32}\rho_{33} ), \\
a_{61} &= \frac{1}{D} ( -\rho_{31}’\varepsilon_{33} + \varepsilon_{31}\rho_{33}’ ), \\
a_{62} &= \frac{1}{D} ( -(\rho_{32}’-\xi)\varepsilon_{33} + \varepsilon_{32}\rho_{33}’ ), \\
a_{64} &= \frac{1}{D} ( -\mu_{31}\varepsilon_{33} + \rho_{31}\rho_{33}’ ), \\
a_{65} &= \frac{1}{D} ( -\mu_{32}\varepsilon_{33} + (\rho_{32}+\xi)\rho_{33}’ ).
\end{aligned}
\end{align}\)
\( \small{
\boldsymbol{\hat{\Delta}}
=
\begin{bmatrix}
\rho_{21}’ +(\rho_{23}’+\xi)a_{31} +\mu_{23}a_{61} & \mu_{22} +(\rho_{23}’+\xi)a_{35} +\mu_{23}a_{65} & \rho_{22}’ +(\rho_{23}’+\xi)a_{32} +\mu_{23}a_{62} & \mu_{21} +(\rho_{23}’+\xi)a_{34} +\mu_{23}a_{64} \\
\varepsilon_{11} +\varepsilon_{13}a_{31} +\rho_{13}a_{61} & \rho_{12} +\varepsilon_{13}a_{35} +\rho_{13}a_{65} & \varepsilon_{12} +\varepsilon_{13}a_{32} +\rho_{13}a_{62} & \rho_{11} +\varepsilon_{13}a_{34} +\rho_{13}a_{64} \\
-\rho_{11}’ -\rho_{13}’a_{31} -\mu_{13}a_{61} & -\mu_{12} -\rho_{13}’a_{35} -\mu_{13}a_{65} & -\rho_{12}’ -\rho_{13}’a_{32} -\mu_{13}a_{62} & -\mu_{11} -\rho_{13}’a_{34} -\mu_{13}a_{64} \\
-\varepsilon_{21} -\varepsilon_{23}a_{31} -(\rho_{23}-\xi)a_{61} & -\rho_{22} -\varepsilon_{23}a_{35} -(\rho_{23}-\xi)a_{65} & -\varepsilon_{22} -\varepsilon_{23}a_{32} -(\rho_{23}-\xi)a_{62} & -\rho_{21} -\varepsilon_{23}a_{34} -(\rho_{23}-\xi)a_{64}
\end{bmatrix}.
} \)
In some future post, we will continue with stratified media and coupling to the ambient through the interface.
1 D. W. Berreman, Optics in Stratified and Anistropic Media: 4×4-Matrix Formulation, Journal of the Optical Society of America, vol 62, 4 (1972).
2 D. W. Berreman, Infrared Absorption at Longitudinal Optic Frequency in Cubic Crystal Films, Physical Review, vol 130, 6 (1963).
3 M. Schubert, Polarization-dependent optical parameters of arbitrarily anisotropic homogeneous layered systems, Physical Review B, vol 53, 8 (1996).