Eigenmodes Identification

Once the permittivity distribution is mapped to the Fourier space, the next step is to apply Maxwell’s equations to identify the eigenmodes of each layer. In this section, we extend the mathematical formulation of the 1D conical incidence case described in [MGPG95] to the 2D grating case as illustrated in Figure 1. To ensure the consistency and clarity, we adopt the same notations and the sign convention of \((+jwt)\).

We consider the normalized excitation wave at the superstrate to take the following form:

Figure 1: Geometry for the stack. Two-dimensional grating layers and incident ray.

(1)\[\begin{align} \mathbf E_{inc} = \mathbf u \cdot e^{-jk_0\mathtt n_{\text{I}}(\sin{\theta} \cdot \cos{\phi}\cdot x + \sin{\theta}\cdot \sin{\phi}\cdot y + \cos{\theta}\cdot z)}, \end{align}\]

where $mathbf u$ is the normalized amplitudes of the wave in each direction:

(2)\[\begin{align} \mathbf{u} = (\cos{\psi}\cdot \cos{\theta}\cdot\cos{\phi} + \sin{\psi}\cdot\sin{\phi}) \hat{x} + (\cos{\psi}\cdot\cos{\theta}\cdot\sin{\phi} + \sin{\psi}\cdot\cos{\phi})\hat{y} + (\cos{\psi}\cdot\sin{\theta})\hat{z}, \end{align}\]

and \(k_0 = 2\pi / \lambda_0\) with \(\lambda_0\) the wavelength of the light in free space, \(\mathtt n_{\text{I}}\) is the refractive index of the superstrate, \(\theta\) is the angle of incidence, \(\phi\) is the rotation (azimuth) angle and \(\psi\) is the angle between the electric field vector and the plane of incidence.

The electric fields in the superstrate and substrate (we will designate these layers by \(\text{I}\) and \(\text{II}\) as in [MGPG95]) can be expressed as a sum of incident, reflected and transmitted waves as the Rayleigh expansion [HSchoberlSZ09, Pet80, Wil07]:

(3)\[\begin{align} \mathbf{E}_{\text{I}} &= \mathbf{E}_{inc} + \sum_{m=-M}^{M} \sum_{n=-N}^{N} \mathbf{R}_{n,m} e^{-j(k_{x,m} x + k_{y,n} y - k_{\text{I}, z,(n,m)}z)}, \end{align}\]

(4)\[\begin{align} \mathbf{E}_{\text{II}} &= \sum_{m=-M}^{M} \sum_{n=-N}^{N} \mathbf{T}_{n,m}e^{-j\{k_{x,m} x + k_{y,n} y + k_{\text{II}, z,(n,m)} (z-d)\}}, \end{align}\]

where \(M\) and \(N\) are the Fourier Truncation Order (FTO) which is related to the number of harmonics in use, and the in-plane components of the wavevector (\(k_{x,m}\) and \(k_{y,n}\)) are determined by the Bloch’s theorem (this has many names and one of them is Floquet condition) [GGFA15, JSGJ08],

\[\begin{split}\begin{align} k_{x,m} &= k_0 \Big(\mathtt n_{\text{I}} \sin{\theta}\cos{\phi} - m\frac{\lambda_0}{\Lambda_x}\Big), \\ k_{y,n} &= k_0 \Big(\mathtt n_{\text{I}} \sin{\theta}\sin{\phi} - n\frac{\lambda_0}{\Lambda_y}\Big), \end{align}\end{split}\]

where \(\Lambda_x\) and \(\Lambda_y\) are the period of the unit cell, and the out-of-plane wavevector is determined from the dispersion relation:

\[\begin{split}\begin{align} k_{\ell,z,(n,m)} = \begin{cases} +\ [(k_0\mathtt n_\ell)^2 - {k_{x,m}}^2 - {k_{y,n}}^2]^{1/2}&, \quad \text{if}\quad ({k_{x,m}}^2 + {k_{y,n}}^2) < (k_0\mathtt n_\ell)^2 \\ -j[{k_{x,m}}^2 + {k_{y,n}}^2 - (k_0\mathtt n_\ell)^2]^{1/2}&, \quad \text{if}\quad ({k_{x,m}}^2 + {k_{y,n}}^2) > (k_0\mathtt n_\ell)^2 \end{cases}, \quad \ell = \text{I}, \text{II}. \end{align}\end{split}\]

Here, \(k_{\ell,z,(n,m)}\) can be categorized into propagation mode and evanescent mode depending on whether it’s real or imaginary. \(\mathbf R_{n,m} \text{ and } \mathbf T_{n,m}\) are the Rayleigh coefficients (also called the reflection and transmission coefficients): \(\mathbf{R}_{n,m}\) is the normalized (3-dimensional) vector of electric field amplitude which is the (\(m^{th}\) in X and \(n^{th}\) in Y) mode of reflected waves in the superstrate and \(\mathbf{T}_{n,m}\) is the normalized (3-dimensional) vector of electric field amplitude which is the (\(m^{th}\) in X and \(n^{th}\) in Y) mode of transmitted waves in the substrate.

Inside the grating layer, the electromagnetic field can be expressed as a superposition of plane waves by the Bloch’s theorem:

(5)\[\begin{split}\begin{align} \mathbf{E}_{g}(x,y,z) &= \sum_{m=-M}^{M} \sum_{n=-N}^{N} \boldsymbol{\mathfrak{S}}_{g,(n,m)} \cdot e^{-j(k_{x,m}x + k_{y,n}y + k_{g,z}z)}, \\ \end{align}\end{split}\]

(6)\[\begin{align} \mathbf{H}_{g}(x,y,z) &= \sum_{m=-M}^{M} \sum_{n=-N}^{N} \boldsymbol{\mathfrak{U}}_{g,(n,m)} \cdot e^{-j(k_{x,m}x + k_{y,n}y + k_{g,z}z)}, \end{align}\]

where \(k_{g,z}\) is the wavevector in Z-direction (this is unique per layer hence the notation g was kept to distinguish) and \(\boldsymbol{\mathfrak{S}}_{g,(n,m)}\) and \(\boldsymbol{\mathfrak{U}}_{g,(n,m)}\) are the vectors of amplitudes in each direction at \((m, n)^{th}\) order:

\[\begin{split}\begin{align} \boldsymbol{\mathfrak{S}}_{g,(n,m)} &= \mathfrak{S}_{g,(n,m), x}\ \hat x + \mathfrak{S}_{g,(n,m), y}\ \hat y + \mathfrak{S}_{g,z}\ \hat z, \\ \boldsymbol{\mathfrak{U}}_{g,(n,m)} &= \mathfrak{U}_{g,(n,m), x}\ \hat x + \mathfrak{U}_{g,(n,m), y}\ \hat y + \mathfrak{U}_{g,z}\ \hat z. \end{align}\end{split}\]

It is also possible to detach wavevector term on \(z\) from exponent and combine with \(\mathbf{\mathfrak{S}}_{g,(n,m)}\) and \(\mathbf{\mathfrak{U}}_{g,(n,m)}\) in Equations (5) and (6) to make \(\mathbf{S}_{g,(n,m)}(z)\) and \(\mathbf{U}_{g,(n,m)}(z)\) which are dependent on \(z\) as shown below:

(7)\[\begin{split}\begin{align} \mathbf{S}_{g,(n,m)}(z) = \boldsymbol{\mathfrak{S}}_{g,(n,m)} \cdot e^{-jk_{g,z}z}, \\ \end{align}\end{split}\]

(8)\[\begin{align} \mathbf{U}_{g,(n,m)}(z) = \boldsymbol{\mathfrak{U}}_{g,(n,m)} \cdot e^{-jk_{g,z}z}, \end{align}\]

then Equations (5) and (6) become

(9)\[\begin{split}\begin{align} \mathbf{E}_{g}(x,y,z) &= \sum_{m=-M}^{M} \sum_{n=-N}^{N} \mathbf{S}_{g,(n,m)}(z) \cdot e^{-j(k_{x,m}x + k_{y,n}y)}, \\ \end{align}\end{split}\]

(10)\[\begin{align} \mathbf{H}_{g}(x,y,z) &= \sum_{m=-M}^{M} \sum_{n=-N}^{N} \mathbf{U}_{g,(n,m)}(z) \cdot e^{-j(k_{x,m}x + k_{y,n}y)}. \end{align}\]

Equations (5) and (6) are used in [LF12, KL23, YR21] and Equations (9) and (10) in [Rum06, MGPG95]. Whichever is used, the result is the same: we will show the development using (\(\boldsymbol{\mathfrak{S}}_{g,(n,m)}\), \(\boldsymbol{\mathfrak{U}}_{g,(n,m)}\)) with the eigendecomposition and then come back to (\(\mathbf{S}_{g,(n,m)}(z)\) and \(\mathbf{U}_{g,(n,m)}(z)\)) with the partial differential equations.

The behavior of the electromagnetic fields can be described by the formulae, called the Maxwell’s equations. Among them, we will use the third and fourth equations,

(11)\[\begin{align} \nabla \times \mathbf E &= -j\omega\mu_0\mathbf H, \end{align}\]

(12)\[\begin{align} \nabla \times \mathbf H &= j\omega\varepsilon_0\varepsilon_r\mathbf E, \end{align}\]

to find the electric and magnetic field inside the grating layer - \(\mathbf E_g\) and \(\mathbf H_g\). Since RCWA is a technique that solves Maxwell’s equations in the Fourier space, curl operator in real space becomes multiplication and multiplication in real space becomes the convolution operator. For this convolution operation, the full set of the modes of the fields and the geometry are required so we introduce a vector notation in the subscript to denote it’s a vector with all the harmonics in use, i.e.,

\[\begin{align} \boldsymbol{{F}}_{g,\vec r} = \begin{bmatrix} {F}_{g,(-N,-M),r} & \cdots & {F}_{g,(-N,M),r} & {F}_{g,(-N+1,-M),r} & \cdots & {F}_{g,(-N+1,M),r} & \cdots & {F}_{g,(N,M),r} \end{bmatrix} ^T, \end{align}\]

where \(\boldsymbol{F} \in \{S, U, \mathfrak{S}, \mathfrak{U}\}\) and \(r \in \{x, y, z\}\). Some variables will be scaled by some factors:

\[\begin{align} \tilde{\mathbf{H}}_g = -j\sqrt{\varepsilon_0/{\mu_0}}\mathbf{H}_g, \quad \tilde k_x = k_x / k_0, \quad \tilde k_y = k_y / k_0, \quad \tilde k_{g,z} = k_{g,z} / k_0, \quad \tilde z = k_0z. \end{align}\]

Substituting Equations (5) and (6) (\(\mathbf E_g\) and \(\tilde{\mathbf{H}}_g\) with \(\boldsymbol{\mathfrak{S}}_g\) and \(\boldsymbol{\mathfrak{U}}_g\)) into Equations (11) and (12) (Maxwell’s equations) and eliminating Z-directional components (\(\mathbf{E}_{g,z}\) and \(\tilde{\mathbf{H}}_{g,z}\)) derive the matrix form of the Maxwell’s equations composed of in-plane components \((\hat x, \hat y)\) in the Fourier space:

(13)\[\begin{split}\begin{align} (-j\tilde{k}_{g,z}) \begin{bmatrix} \boldsymbol{\mathfrak{S}}_{g,\vec x} \cdot e^{-j\tilde{k}_{g,z}\tilde{z}} \\ \boldsymbol{\mathfrak{S}}_{g,\vec y} \cdot e^{-j\tilde{k}_{g,z}\tilde{z}} \end{bmatrix} =\boldsymbol{{\Omega}}_{g,L} \begin{bmatrix} \boldsymbol{\mathfrak{U}}_{g, \vec x} \cdot e^{-j\tilde{k}_{g,z}\tilde{z}} \\ \boldsymbol{\mathfrak{U}}_{g, \vec y} \cdot e^{-j\tilde{k}_{g,z}\tilde{z}} \end{bmatrix} \end{align}\end{split}\]

(14)\[\begin{split}\begin{align} (-j\tilde{k}_{g,z}) \begin{bmatrix} \boldsymbol{\mathfrak{U}}_{g,\vec x} \cdot e^{-j\tilde{k}_{g,z}\tilde{z}} \\ \boldsymbol{\mathfrak{U}}_{g,\vec y} \cdot e^{-j\tilde{k}_{g,z}\tilde{z}} \end{bmatrix} = \boldsymbol{{\Omega}}_{g,R} \begin{bmatrix} \boldsymbol{\mathfrak{S}}_{g,\vec x} \cdot e^{-j\tilde{k}_{g,z}\tilde{z}} \\ \boldsymbol{\mathfrak{S}}_{g,\vec y} \cdot e^{-j\tilde{k}_{g,z}\tilde{z}} \end{bmatrix} \end{align}\end{split}\]

(15)\[\begin{split}\begin{align} \label{eqn:curlyS-omega_LR} (-j\tilde{k}_{g,z})^2 \begin{bmatrix} \boldsymbol{\mathfrak{S}}_{g,\vec x} \cdot e^{-j\tilde{k}_{g,z}\tilde{z}} \\ \boldsymbol{\mathfrak{S}}_{g,\vec y} \cdot e^{-j\tilde{k}_{g,z}\tilde{z}} \end{bmatrix} = \boldsymbol{{\Omega}}_{g,LR}^2 \begin{bmatrix} \boldsymbol{\mathfrak{S}}_{g,\vec x} \cdot e^{-j\tilde{k}_{g,z}\tilde{z}} \\ \boldsymbol{\mathfrak{S}}_{g,\vec y} \cdot e^{-j\tilde{k}_{g,z}\tilde{z}} \end{bmatrix} \end{align}\end{split}\]

where

(16)\[\begin{split}\begin{align} \boldsymbol{{\Omega}}_{g,L} = \begin{bmatrix} (-\tilde{\mathbf K}_x \left[\!\!\left[\varepsilon_{r,g}\right]\!\!\right] ^{-1}\tilde{\mathbf K}_y) & (\tilde{\mathbf K}_x\left[\!\!\left[\varepsilon_{r,g}\right]\!\!\right]^{-1}\tilde{\mathbf K}_x - \mathbf I) \\ (\mathbf I-\tilde{\mathbf K}_y\left[\!\!\left[\varepsilon_{r,g}\right]\!\!\right]^{-1}\tilde{\mathbf K}_y) & (\tilde{\mathbf K}_y\left[\!\!\left[\varepsilon_{r,g}\right]\!\!\right]^{-1}\tilde{\mathbf K}_x) \end{bmatrix}, \end{align}\end{split}\]

(17)\[\begin{split}\begin{align} \boldsymbol{{\Omega}}_{g,R} = \begin{bmatrix} (-\tilde{\mathbf K}_x\tilde{\mathbf K}_y) & (\tilde{\mathbf K}_x^2 - \left[\!\!\left[\varepsilon_{r,g}\right]\!\!\right]) \\ (\left[\!\!\left[\varepsilon_{r,g}^{-1}\right]\!\!\right] ^{-1} - \tilde{\mathbf K}_y^2) & (\tilde{\mathbf K}_y\tilde{\mathbf K}_x) \end{bmatrix}, \end{align}\end{split}\]

(18)\[\begin{split}\begin{align} \boldsymbol{{\Omega}}_{g,LR}^2 = \begin{bmatrix} {\tilde{\mathbf K}_y}^2 + (\tilde{\mathbf{K}}_x \left[\!\!\left[\varepsilon_{r,g}\right]\!\!\right]^{-1} \tilde{\mathbf{K}}_x - \mathbf{I}) \left[\!\!\left[\varepsilon_{r,g}^{-1}\right]\!\!\right]^{-1} & \tilde{\mathbf K}_x(\left[\!\!\left[\varepsilon_{r,g}\right]\!\!\right]^{-1}\tilde{\mathbf K}_y\left[\!\!\left[\varepsilon_{r,g}\right]\!\!\right] - \tilde{\mathbf K}_y) \\ \tilde{\mathbf K}_y(\left[\!\!\left[\varepsilon_{r,g}\right]\!\!\right]^{-1}\tilde{\mathbf K}_x\left[\!\!\left[\varepsilon_{r,g}^{-1}\right]\!\!\right]^{-1} - \tilde{\mathbf K}_x) & {\tilde{\mathbf K}_x}^2 + (\tilde{\mathbf{K}}_y \left[\!\!\left[\varepsilon_{r,g}\right]\!\!\right]^{-1} \tilde{\mathbf K}_y - \mathbf{I})\left[\!\!\left[\varepsilon_{r,g}\right]\!\!\right] \end{bmatrix}, \end{align}\end{split}\]

and

\[\begin{split}\begin{align} \tilde{\mathbf{K}}_r = \begin{bmatrix} \tilde k_{r,(-N,-M)} & 0 & \cdots & 0 \\ 0 & \tilde k_{r,(-N,-M+1)} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0& \cdots & \tilde k_{r,(N,M)} \end{bmatrix}, \quad r \in \{x, y\}, \end{align}\end{split}\]

and \(\left[\!\!\left[ ~~ \right]\!\!\right]\) is the convolution (a.k.a Toeplitz) matrix: \(\left[\!\!\left[\varepsilon_{r,g}\right]\!\!\right]\) and \(\left[\!\!\left[\varepsilon_{r,g}^{-1}\right]\!\!\right]^{-1}\) are convolution matrices composed of Fourier coefficients of permittivity and one-over-permittivity (by the inverse rule presented in [Li96] and [Li14]).

Equation (15) is a typical form of the eigendecomposition of a matrix.

The vector [\(\boldsymbol{\mathfrak{S}}_{g,\vec x} \cdot e^{-j\tilde{k}_{g,z}\tilde{z}} \quad \boldsymbol{\mathfrak{S}}_{g,\vec y} \cdot e^{-j\tilde{k}_{g,z}\tilde{z}}]^T\) is an eigenvector of \(\boldsymbol{{\Omega}}_{g,LR}^2\) and \(j\tilde k_{g,z}\) is the positive square root of the eigenvalues. This intuitively shows how the eigenvalues are connected to the Z-directional wavevectors.

It is also possible to use \(\mathbf S_{g,\vec x}(\tilde z)\) and \(\mathbf S_{g,\vec y}(\tilde z)\) instead of \(\boldsymbol{\mathfrak{S}}_{g,\vec x}\) and \(\boldsymbol{\mathfrak{U}}_{g,\vec x}\) because they satisfy the following relations:

\[\begin{split}\begin{align} \frac{\partial^2}{\partial(\tilde z)^2} \begin{bmatrix} \mathbf S_{g,\vec x}(\tilde z)\\ \mathbf S_{g,\vec y}(\tilde z) \end{bmatrix} = \frac{\partial^2}{\partial(\tilde z)^2} % (-j\tilde{k}_z)^2 \begin{bmatrix} \boldsymbol{\mathfrak{S}}_{g,\vec x} \cdot e^{-j\tilde{k}_{g,z}\tilde{z}} \\ \boldsymbol{\mathfrak{S}}_{g,\vec y} \cdot e^{-j\tilde{k}_{g,z}\tilde{z}} \end{bmatrix} = (-j\tilde{k}_{g,z})^2 \begin{bmatrix} \boldsymbol{\mathfrak{S}}_{g,\vec x} \cdot e^{-j\tilde{k}_{g,z}\tilde{z}} \\ \boldsymbol{\mathfrak{S}}_{g,\vec y} \cdot e^{-j\tilde{k}_{g,z}\tilde{z}} \end{bmatrix}. \end{align}\end{split}\]

Hence it is just a matter of choice and we will use PDE form (\(\mathbf S_g\) and \(\mathbf U_g\)) for the seamless connection to the 1D conical case in the previous work [MGPG95]. Then Equations (13), (14) and (15) become

(19)\[\begin{split}\begin{align} \frac{\partial}{\partial(\tilde z)} \begin{bmatrix} \mathbf S_{g,\vec x}(\tilde z)\\ \mathbf S_{g,\vec y}(\tilde z) \end{bmatrix} = \boldsymbol{{\Omega}}_{g,L} \begin{bmatrix} \mathbf U_{g,\vec x}(\tilde z) \\ \mathbf U_{g,\vec y}(\tilde z) \end{bmatrix}, \end{align}\end{split}\]

(20)\[\begin{split}\begin{align} \frac{\partial}{\partial(\tilde z)} \begin{bmatrix} \mathbf U_{g,\vec x}(\tilde z)\\ \mathbf U_{g,\vec y}(\tilde z) \end{bmatrix} = \boldsymbol{{\Omega}}_{g,R} \begin{bmatrix} \mathbf S_{g,\vec x}(\tilde z) \\ \mathbf S_{g,\vec y}(\tilde z) \end{bmatrix}, \end{align}\end{split}\]

(21)\[\begin{split}\begin{align} \frac{\partial^2}{\partial(\tilde z)^2} \begin{bmatrix} \mathbf S_{g,\vec x}(\tilde z)\\ \mathbf S_{g,\vec y}(\tilde z) \end{bmatrix} = \boldsymbol{{\Omega}}_{g,LR}^2 \begin{bmatrix} \mathbf S_{g,\vec x}(\tilde z) \\ \mathbf S_{g,\vec y}(\tilde z) \end{bmatrix}, \end{align}\end{split}\]

where Equation (21) is the second order matrix differential equation which has the general solution of the following form

\[\begin{split}\begin{align} \begin{bmatrix} \mathbf{S}_{g,\vec{x}}(\tilde z) \\ \mathbf{S}_{g,\vec y}(\tilde z) \end{bmatrix} = \boldsymbol{w}_{g,1}&(c_{g,1}^+ e^{-q_{g,1}\tilde z} + c_{g,1}^- e^{+q_{g,1}\tilde z}) % + \boldsymbol{w}_{g,2}(c_{g,2}^+ e^{-q_{g,2}\tilde z} + c_{g,2}^- e^{+q_{g,2}\tilde z}) + \cdots + \boldsymbol{w}_{g,\xi}(c_{g,\xi}^+ e^{-q_{g,\xi} \tilde z} + c_{g,\xi}^- e^{+q_{g,\xi} \tilde z}) \\ &= \sum_{i=1}^{\xi} \boldsymbol{w}_{g,i}(c_{g,i}^+e^{-q_{g,i}\tilde z} + c_{g,i}^-e^{+q_{g,i}\tilde z}), \end{align}\end{split}\]

where \(\xi=(2M+1)(2N+1)\), the total number of harmonics, and \(\boldsymbol{w}_g\) is the eigenvector, \(q_g\) is the positive square root of the corresponding eigenvalue (\(j\tilde{k}_{g,z}\)) and \(c_g^\pm\) are the coefficients (amplitudes) of the mode in each propagating direction (+Z and -Z direction). This can be written in matrix form

(22)\[\begin{split}\begin{align} \begin{bmatrix} \mathbf{S}_{g,\vec{x}}(\tilde z) \\ \mathbf{S}_{g,\vec y}(\tilde z) \end{bmatrix} &= \mathbf{W}_g \mathbf{Q}_g^- \mathbf{c}_g^+ + \mathbf{W}_g \mathbf{Q}_g^+ \mathbf{c}_g^{-} \\ &= \mathbf W_g \begin{bmatrix} \mathbf Q_g^- & \mathbf Q_g^+ \\ \end{bmatrix} \begin{bmatrix} {\mathbf c}_g^+ \\ {\mathbf c}_g^- \\ \end{bmatrix}, \\ \label{eqn:S=WQC} &= \begin{bmatrix} \mathbf W_{g,11} & \mathbf W_{g,12} \\ \mathbf W_{g,21} & \mathbf W_{g,22} \end{bmatrix} \begin{bmatrix} {\mathbf{Q}_{g,1}^-} & 0 & {\mathbf{Q}_{g,1}^+} & 0 \\ 0 & {\mathbf{Q}_{g,2}^-} & 0 & {\mathbf{Q}_{g,2}^+} \end{bmatrix} \begin{bmatrix} \mathbf c_{g,1}^+ \\ \mathbf c_{g,2}^+ \\ \mathbf c_{g,1}^- \\ \mathbf c_{g,2}^- \end{bmatrix}, \end{align}\end{split}\]

where \(\mathbf Q_g^\pm\) are the diagonal matrices with the exponential of eigenvalues

\[\begin{split}\begin{align} \mathbf Q_g^\pm = \begin{bmatrix} e^{\pm{q}_{g,1}} & & 0 \\ & \ddots & \\ 0 & & e^{{\pm q_{g,\xi}}} \end{bmatrix}, \end{align}\end{split}\]

and \(\mathbf W_g\) is the matrix that has the eigenvectors in columns and \(\mathbf c_g^\pm\) are the vectors of the coefficients.

Now we can find the general solution of the magnetic field that shares same \(\mathbf Q_g\) and \(\mathbf c_g^\pm\) with the electric field in corresponding mode. It can be written in a similar form of Equation (22) as

(23)\[\begin{split}\begin{align} \begin{bmatrix} \mathbf{U}_{g,\vec{x}}(\tilde z) \\ \mathbf{U}_{g,\vec y}(\tilde z) \end{bmatrix} &= -\mathbf{V}_g \mathbf{Q}_g^- \mathbf{c}_g^+ + \mathbf{V}_g \mathbf{Q}_g^+ \mathbf{c}_g^{-}. \end{align}\end{split}\]

The negative sign in the first term was given to adjust the direction of the curl operation, \(E \times H\), to be in accordance with the wave propagation direction, \(\tilde k_{g,z}\). By substituting Equations (22) and (23) into Equation (20), we can get

\[\begin{align} \mathbf V_g = \boldsymbol{{\Omega}}_{g,R} \mathbf W_g \mathbf q_g^{-1}, \end{align}\]

where \(\mathbf q_g\) is the diagonal matrix with the eigenvalues. This can be written in matrix form

\[\begin{split}\begin{equation} \begin{split} \mathbf{V}_g = \begin{bmatrix} \mathbf V_{g,11} & \mathbf V_{g,12} \\ \mathbf V_{g,21} & \mathbf V_{g,22} \end{bmatrix} = \begin{bmatrix} -\tilde{\mathbf K}_x \tilde{\mathbf K}_y & \tilde{{\mathbf K}}_x^2-\left[\!\!\left[\varepsilon_{r,g}\right]\!\!\right] \\ \left[\!\!\left[\varepsilon_{r,g}^{-1}\right]\!\!\right]^{-1} - \tilde{{\mathbf K}}_y^2 & \tilde{\mathbf K}_y \tilde{\mathbf K}_x \end{bmatrix} \begin{bmatrix} \mathbf W_{g,11} & \mathbf W_{g,12} \\ \mathbf W_{g,21} & \mathbf W_{g,22} \end{bmatrix} \begin{bmatrix} \mathbf{q}_{g,1} & 0 \\ 0 & \mathbf{q}_{g,2} \end{bmatrix}^{-1} \end{split}. \end{equation}\end{split}\]

---

[GGFA15]

Pablo Gómez García and José-Paulino Fernández-Álvarez. Floquet-bloch theory and its application to the dispersion curves of nonperiodic layered systems. Mathematical Problems in Engineering, 2015:475364, 2015. URL: https://doi.org/10.1155/2015/475364.

[HSchoberlSZ09]

M. Huber, J. Schöberl, A. Sinwel, and S. Zaglmayr. Simulation of diffraction in periodic media with a coupled finite element and plane wave approach. SIAM Journal on Scientific Computing, 31(2):1500–1517, 2009.

[JSGJ08]

John D. Joannopoulos and Robert D. Meade Steven G. Johnson, Joshua N. Winn. Photonic Crystals: Molding the Flow of Light - Second Edition. Princeton University Press, 2nd edition, 2008. ISBN 0691124566, 978-0691124568. URL: http://ab-initio.mit.edu/book/photonic-crystals-book.pdf.

[KL23]

Changhyun Kim and Byoungho Lee. Torcwa: gpu-accelerated fourier modal method and gradient-based optimization for metasurface design. Computer Physics Communications, 282:108552, 2023.

[Li96]

Lifeng Li. Use of fourier series in the analysis of discontinuous periodic structures. J. Opt. Soc. Am. A, 13(9):1870–1876, Sep 1996.

[Li14]

Lifeng Li. Fourier Modal Method. In ed. E. Popov, editor, Gratings: Theory and Numeric Applications, Second Revisited Edition, pages 13.1–13.40. (Aix Marseille Université, CNRS, Centrale Marseille, Institut Fresnel, April 2014. URL: https://hal.science/hal-00985928.

[LF12]

Victor Liu and Shanhui Fan. S4: a free electromagnetic solver for layered periodic structures. Computer Physics Communications, 183(10):2233–2244, 2012.

[MGPG95] (1,2,3,4)

M. G. Moharam, Eric B. Grann, Drew A. Pommet, and T. K. Gaylord. Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings. J. Opt. Soc. Am. A, 12(5):1068–1076, May 1995.

[Pet80]

R. Petit. Electromagnetic theory of gratings. Springer-Verlag, 1980. ISBN 0387101934.

[Rum06]

Raymond Rumpf. Design and optimization of nano-optical elements by coupling fabrication to optical behavior. Electronic Theses and Dissertations, 2004-2019, 2006. URL: https://stars.library.ucf.edu/etd/1081.

[Wil07]

Strutt John William. On the dynamical theory of gratings. Proc. R. Soc. Lond. A, 79:399–416, 1907.

[YR21]

Gwanho Yoon and Junsuk Rho. Maxim: metasurfaces-oriented electromagnetic wave simulation software with intuitive graphical user interfaces. Computer Physics Communications, 264:107846, 2021.