Enhanced Transmittance Matrix Method
As addressed in [Li93, MPGG95, PNeviere00],
solving below equation,
\[\begin{split}\begin{align}
\begin{bmatrix}
\sin\psi\ \boldsymbol\delta_{00} \\
\cos\psi\ \cos\theta\ \boldsymbol\delta_{00}
\\
j\sin\psi\ \mathtt n_{\text{I}} \cos\theta\ \boldsymbol\delta_{00} \\
-j\cos\psi\ \mathtt n_{\text{I}}\ \boldsymbol\delta_{00} \\
\end{bmatrix}
+
\begin{bmatrix}
\mathbf I & \mathbf 0 \\
\mathbf 0 & -j\mathbf Z_I \\
-j\mathbf Y_I & \mathbf 0 \\
\mathbf 0 & \mathbf I
\end{bmatrix}
\begin{bmatrix}
\mathbf R_s \\
\mathbf R_p
\end{bmatrix}
% \\\quad\\
% \begin{align*}
=
\prod_{\ell=1}^{L}
\begin{bmatrix}
\mathbb W_\ell & \mathbb {W_\ell X_\ell} \\
\mathbb V_\ell & -\mathbb {V_\ell X_\ell}
\end{bmatrix}
\begin{bmatrix}
\mathbb {W_\ell X_\ell} & \mathbb W_\ell \\
\mathbb {V_\ell X_\ell} & -\mathbb V_\ell
\end{bmatrix}^{-1}
\begin{bmatrix}
\mathbb F_{L+1} \\
\mathbb G_{L+1} \\
\end{bmatrix}
\begin{bmatrix}
\mathbf T_s \\ \mathbf T_p
\end{bmatrix},
% \end{align*}
\end{align}\end{split}\]
may suffer from the numerical instability coming from the inversion of almost singular matrix when
\(\mathbb X_\ell\) has a very small and possibly numerically zero value.
Meent adopted Enhanced Transmittance Matrix Method (ETM) [MPGG95] to overcome this
by avoiding the inversion of \(\mathbb X_\ell\).
The technique is sequentially applied from the last layer to the first layer.
In the equation, the set of modes at the bottom interface of the last layer \((\ell = L)\) is
(1)\[\begin{split}\begin{equation}
\begin{split}
&\begin{bmatrix}
\mathbb W_L & \mathbb{W}_L \mathbb X_L \\
\mathbb V_L & -\mathbb{V}_L \mathbb X_L
\end{bmatrix}
\begin{bmatrix}
\mathbb W_L \mathbb X_L & \mathbb W_L\\
\mathbb V_L \mathbb X_L & -\mathbb V_L
\end{bmatrix}^{-1}
\begin{bmatrix}
\mathbb F_{L+1} \\
\mathbb G_{L+1}
\end{bmatrix}
\begin{bmatrix}
\mathbf T_s \\ \mathbf T_p
\end{bmatrix}
\\
&=
\begin{bmatrix}
\mathbb W_L & \mathbb W_L \mathbb X_L \\
\mathbb V_L & -\mathbb V_L \mathbb X_L
\end{bmatrix}
\begin{bmatrix}
{\mathbb X_L}^{-1} & \mathbb{0} \\
\mathbb{0} & {\mathbb I} \\
\end{bmatrix}
{
\begin{bmatrix}
\mathbb W_L & \mathbb W_L \\
\mathbb V_L & -\mathbb V_L
\end{bmatrix}
}^{-1}
\begin{bmatrix}
\mathbb F_{L+1} \\ \mathbb G_{L+1}
\end{bmatrix}
\begin{bmatrix}
\mathbf T_s \\ \mathbf T_p
\end{bmatrix}.
\\
\end{split}
\end{equation}\end{split}\]
The matrix to be inverted can be decomposed into two matrices by isolating \(\mathbb X_L\),
which is the potential source of the numerical instability.
The right-hand side can be shortened with new variables \(\mathbb A_L, \mathbb B_L\):
\[\begin{split}\begin{equation}
\begin{split}
\begin{bmatrix}
\mathbb A_L \\
\mathbb B_L
\end{bmatrix}
=
\begin{bmatrix}
{\mathbb W_L} & \mathbb{W_L} \\
\mathbb{V_L} & {-\mathbb V_L} \\
\end{bmatrix}^{-1}
\begin{bmatrix}
\mathbb F_{L+1} \\ \mathbb G_{L+1}
\end{bmatrix},
\end{split}
\end{equation}\end{split}\]
then the right-hand side of Equation (1) becomes
(2)\[\begin{split}\begin{equation}
\begin{split}
\begin{bmatrix}
\mathbb W_L & \mathbb W_L \mathbb X_L \\
\mathbb V_L & -\mathbb V_L \mathbb X_L
\end{bmatrix}
\begin{bmatrix}
{\mathbb X_L}^{-1} & \mathbb{0} \\
\mathbb{0} & {\mathbb I} \\
\end{bmatrix}
\begin{bmatrix}
\mathbb A_L \\ \mathbb B_L
\end{bmatrix}
\begin{bmatrix}
\mathbf T_s \\ \mathbf T_p
\end{bmatrix}.
\\
\end{split}
\end{equation}\end{split}\]
We can avoid the inversion of \(\mathbb X_L\) by introducing the substitution
\(\mathbf T_s = {\mathbb A_L}^{-1} \mathbb X_L \mathbf T_{s,L}\) and
\(\mathbf T_p = {\mathbb A_L}^{-1} \mathbb X_L \mathbf T_{p,L}\).
Equation (2) then becomes
\[\begin{split}\begin{equation}
\begin{split}
&\begin{bmatrix}
\mathbb W_L & \mathbb W_L \mathbb X_L \\
\mathbb V_L & -\mathbb V_L \mathbb X_L
\end{bmatrix}
\begin{bmatrix}
{\mathbb X_L}^{-1} & \mathbb{0} \\
\mathbb{0} & {\mathbb I} \\
\end{bmatrix}
\begin{bmatrix}
\mathbb A_L \\ \mathbb B_L
\end{bmatrix}
{\mathbb A_L}^{-1}{\mathbb X_L}
\begin{bmatrix}
\mathbf T_{s,L} \\ \mathbf T_{p,L}
\end{bmatrix}
\\
&=
\begin{bmatrix}
\mathbb W_L & \mathbb W_L \mathbb X_L \\
\mathbb V_L & -\mathbb V_L \mathbb X_L
\end{bmatrix}
\begin{bmatrix}
{\mathbb X_L}^{-1} & \mathbb{0} \\
\mathbb{0} & {\mathbb I} \\
\end{bmatrix}
\begin{bmatrix}
\mathbb X_L \\ \mathbb B_L \mathbb A_L^{-1} \mathbb X_L
\end{bmatrix}
\begin{bmatrix}
\mathbf T_{s,L} \\ \mathbf T_{p,L}
\end{bmatrix}
\\
&=
\begin{bmatrix}
\mathbb W_L & \mathbb W_L \mathbb X_L \\
\mathbb V_L & -\mathbb V_L \mathbb X_L
\end{bmatrix}
\begin{bmatrix}
\mathbb I \\ \mathbb B_L \mathbb A_L^{-1} \mathbb X_L
\end{bmatrix}
\begin{bmatrix}
\mathbf T_{s,L} \\ \mathbf T_{p,L}
\end{bmatrix}
\\
&=
\begin{bmatrix}
\mathbb W_L(\mathbb I+\mathbb X_L \mathbb B_L \mathbb A_L^{-1} \mathbb X) \\
\mathbb V_L(\mathbb I-\mathbb X_L \mathbb B_L \mathbb A_L^{-1} \mathbb X)
\end{bmatrix}
\begin{bmatrix}
\mathbf T_{s,L} \\ \mathbf T_{p,L}
\end{bmatrix}
\\
&=
\begin{bmatrix}
\mathbb F_L \\ \mathbb G_L
\end{bmatrix}
\begin{bmatrix}
\mathbf T_{s,L} \\ \mathbf T_{p,L}
\end{bmatrix}
.
\end{split}
\end{equation}\end{split}\]
These steps can be repeated until the iteration gets to the first layer \((\ell = 1)\), then the form becomes
\[\begin{split}\begin{align}
\begin{bmatrix}
\sin\psi\ \boldsymbol\delta_{00} \\
\cos\psi\ \cos\theta\ \boldsymbol\delta_{00}
\\
j\sin\psi\ n_{\text{I}}\ \cos\theta\ \boldsymbol\delta_{00} \\
-j\cos\psi\ n_{\text{I}}\ \boldsymbol\delta_{00} \\
\end{bmatrix}
+
\begin{bmatrix}
\mathbf I & \mathbf 0 \\
\mathbf 0 & -j\mathbf Z_I \\
-j\mathbf Y_I & \mathbf 0 \\
\mathbf 0 & \mathbf I
\end{bmatrix}
\begin{bmatrix}
\mathbf R_s \\
\mathbf R_p
\end{bmatrix}
=
\begin{bmatrix}
\mathbb F_1 \\ \mathbb G_1
\end{bmatrix}
\begin{bmatrix}
\mathbf T_{s,1} \\ \mathbf T_{p,1}
\end{bmatrix},
\end{align}\end{split}\]
where
\[\begin{split}\begin{bmatrix}
\mathbf T_s \\ \mathbf T_p
\end{bmatrix}
=
\mathbb A_L^{-1} \mathbb X_L \cdots
\mathbb A_\ell^{-1} \mathbb X_\ell \cdots
\mathbb A_1^{-1} \mathbb X_1
\begin{bmatrix}
\mathbf T_{s,1} \\ \mathbf T_{p,1}
\end{bmatrix}.\end{split}\]
[Li93]
Lifeng Li. Multilayer modal method for diffraction gratings of arbitrary profile, depth, and permittivity. JOSA A, 10(12):2581–2591, 1993.
[MPGG95]
(1,2)
M. G. Moharam, Drew A. Pommet, Eric B. Grann, and T. K. Gaylord. Stable implementation of the rigorous coupled-wave analysis for surface-relief gratings: enhanced transmittance matrix approach. J. Opt. Soc. Am. A, 12(5):1077–1086, May 1995.
[PNeviere00]
Evgeni Popov and Michel Nevière. Grating theory: new equations in fourier space leading to fast converging results for tm polarization. J. Opt. Soc. Am. A, 17(10):1773–1784, Oct 2000.