Rigorous Coupled-Wave Analysis Solver (RCWA)

RCWA expands periodic media and periodic fields into Fourier series, transforming Maxwell's equations in continuous space into a set of linear differential equations in discrete Fourier space. Since the periodic structure is discrete in the plane, the fields couple between different diffraction orders, hence this method is also called "Rigorous Coupled-Wave Analysis". The main physical principles are as follows:

1. Fourier Expansion of Periodic Media

In two-dimensional periodic gratings, the permittivity and permeability satisfy

$$\epsilon_r(x,y)=\epsilon_r(x+\Lambda_x,y)=\epsilon_r(x,y+\Lambda_y),\\ \mu_r(x,y)=\mu_r(x+\Lambda_x,y)=\mu_r(x,y+\Lambda_y)$$

Thus, they can be expanded into Fourier series:

$$\epsilon_r(x,y)=\sum_{m,n} \epsilon_{m,n}\,e^{-j(mG_x x+nG_y y)},\\ \mu_r(x,y)=\sum_{m,n} \mu_{m,n}\,e^{-j(mG_x x+nG_y y)}$$

here, $G_x=2\pi/\Lambda_x$ and $G_y=2\pi/\Lambda_y$ are the components of the reciprocal lattice vectors.

2. Expansion of Electromagnetic Fields into Spatial Harmonics

Let the wave vector of the incident wave in the $x$-$y$ plane be $(k_{x0},k_{y0})$, then the wave vector components corresponding to the $(m,n)$ diffraction order are

$$k_{x,m}=k_{x0}+mG_x,\quad k_{y,n}=k_{y0}+nG_y,\quad k_{z,mn}=\sqrt{k_0^2-k_{x,m}^2-k_{y,n}^2}$$

here, $k_0=2\pi/\lambda$ is the free-space wave number. Therefore, the electric and magnetic fields can be expressed as a superposition of spatial harmonics:

$$\mathbf{E}(x,y,z)=\sum_{m,n}\mathbf{S}_{m,n}(z)e^{-j(k_{x,m}x+k_{y,n}y)}$$

$$\mathbf{H}(x,y,z)=\sum_{m,n}\mathbf{U}_{m,n}(z)e^{-j(k_{x,m}x+k_{y,n}y)}$$

here, $\mathbf{S}_{m,n}(z)$ and $\mathbf{U}_{m,n}(z)$ are the amplitude functions of each Fourier mode at the given $z$ position.

3. Substituting into Maxwell's equations and converting them into an eigenvalue problem

Substituting the above expansions into the frequency-domain Maxwell's equations

$$\nabla\times\mathbf{E}=j\omega\mu_0\mu_r\mathbf{H},\quad \nabla\times\mathbf{H}=-j\omega\epsilon_0\epsilon_r\mathbf{E}$$

Since the permittivity and permeability are periodic functions, multiplication becomes convolution in Fourier space, so different diffraction orders can couple with each other. For each $(m,n)$ order, one obtains a set of first-order differential equations for $\mathbf{S}_{m,n}(z)$ and $\mathbf{U}_{m,n}(z)$, whose general form can be written as

$$\frac{d}{dz}\mathbf{\Psi}(z)=\mathbf{\Omega^2}\mathbf{\Psi}(z)$$

where

$$\mathbf{\Psi}(z)= \begin{bmatrix} \mathbf{S}_{m,n}(z) \\ \mathbf{U}_{m,n}(z) \end{bmatrix}$$

The elements of the matrix $\mathbf{\Omega^2}$ are determined by the wavevector components $k_{x,m},k_{y,n}$, the material parameters $\epsilon_{m,n},\mu_{m,n}$, and the frequency $\omega$. In this way, the matrix essentially describes the coupling between different diffraction orders. For each homogeneous layer, the constant-coefficient matrix $\mathbf{\Omega^2}$ can be further used to solve the eigenvalue problem:

$$\mathbf{\Omega^2}\mathbf{v}_\ell=\eta_\ell\mathbf{v}_\ell$$

where $\eta_\ell$ is the $\ell$-th propagation eigenvalue, and $\mathbf{v}_\ell$ is the corresponding eigenmode. In other words, the field inside a given layer is not a single plane wave, but a linear superposition of a set of propagation modes that satisfy the same boundary conditions. The field components in any layer can therefore be written as

$$\mathbf{\Psi}(z)=\sum_{\ell} c_\ell\mathbf{v}_\ell e^{j\eta_\ell z}$$

where $c_\ell$ are the unknown expansion coefficients.

1. Boundary conditions and matrix cascading

The interface boundary conditions are used to “connect” these modes from adjacent layers at the interface. If the interface is located at $z=z_0$, the tangential electric and magnetic fields must be continuous based on Maxwell’s boundary conditions, leading to a set of matching equations for the coefficients $c_\ell^{(l)}$ and $c_\ell^{(l+1)}$:

$$\sum_{\ell} c_\ell^{(l)}\mathbf{f}_\ell^{(l)}e^{j\eta_\ell^{(l)}z_0} = \sum_{\ell} c_\ell^{(l+1)}\mathbf{f}_\ell^{(l+1)}e^{j\eta_\ell^{(l+1)}z_0}$$

where $\mathbf{f}_\ell^{(l)}$ denotes the tangential field vector of the $\ell$-th propagation mode at the interface. Rewriting this in matrix form gives the linear system

$$\mathbf{M}^{(l)}\mathbf{c}^{(l)}=\mathbf{M}^{(l+1)}\mathbf{c}^{(l+1)}$$

By solving the linear system formed by this interface matching, one can obtain the relationships between the eigenmode coefficients of adjacent layers, and thereby determine the reflection and transmission coefficients of each layer, as well as the amplitude of each diffraction order. For multilayer structures, one can first obtain the local scattering matrix $S$ of each layer, and then use the transmission matrix or the Redheffer star product to combine the layers and obtain the total reflection and transmission properties of the entire device.

In summary, the essence of RCWA is to decompose a complex periodic structure into a set of coupled discrete Fourier diffraction orders. An incident wave excites several diffraction orders inside the grating, and these orders couple to one another through the Fourier coefficients of the periodic medium. By solving the matrix eigenvalue problem, one can accurately obtain the propagation, reflection, and transmission characteristics of each order. Therefore, RCWA is particularly well suited for analyzing periodic gratings, metasurfaces, Bragg reflectors, and other optical devices with spatial periodicity.