Rigorous Coupled-Wave Analysis Algorithm for Stratified Two-Dimensional Gratings with Unconditionally Stable H-Matrix Methods

Abstract

Nanostructures with the two-dimensional (2D) periodicity are attracting increasing attention due to their promising applications in planar optical devices and their potential for scalable industrial production. While Rigorous Coupled-Wave Analysis (RCWA) has proven to be an efficient electromagnetic solver for simulating the diffraction of large-scale periodic nanostructures, it has been largely applied in nanostructures with one-dimensional (1D) periodicity and suffers from potentially low computational stability. In this study, we present a step-by-step formulation of the RCWA algorithm for 2D stratified grating structures. Through dimensionality reduction, we show that the boundary conditions in 2D gratings can be transformed into forms analogous to those of 1D gratings. Additionally, we implement a hybrid matrix algorithm to enhance the computational stability of the RCWA. The stability and accuracy of the hybrid matrix-enhanced RCWA algorithm are compared with other recursive methods. An exemplary application in metalens demonstrates the effectiveness of our algorithms.

1. Introduction

Being a type of two-dimensional (2D) periodic structure, metasurfaces have emerged as a versatile tool for controlling light at the nanoscale [1,2,3]. These periodic subwavelength nanostructures enable the unprecedented manipulation of the phase [4], amplitude [5], and polarization [6] of electromagnetic waves, paving the way for compact and high-performance optical devices [7,8,9,10]. Two-dimensional gratings are one kind of metasurfaces, which pattern the nanostructures periodically in the two directions of the in-plane. They have garnered increasing attention from researchers due to the development of nanofabrication technologies [11]. For the analysis of such 2D periodic structures, several computational approaches are commonly employed, such as Finite-Difference Time Domain [12], boundary element method (BEM) [13,14], the Finite Element Method [15], and Rigorous Coupled-Wave Analysis (RCWA) [16,17,18,19,20]. As a semi-analytical electromagnetic solver, RCWA has the advantages over other approaches of both accuracy and efficiency if a large area (cm² or beyond) of nanostructures is required [11]. For example, the cutting-edge augmented-reality glasses based on diffractive optical elements [21] and the inversely designed [22,23] optical devices usually require numerous computational resources for full-wave simulations, making other numerical approaches almost impossible.

Tamir, Wang, and Oliver [24] appear to have first applied the RCWA to grating problems in the mid-1960s. From 1981 to 1986, Moharam and Gaylord used the RCWA to analyze slanted gratings [25], conical gratings [26], and metallic gratings [27]. In 1996, Lalanne and Morris proposed a novel formulation that significantly enhanced the convergence rates for TM polarization in 1D gratings [17]. Following this work, Lalanne introduced a new formula applicable to 2D gratings [18]. In 1996, Li established a theoretical framework [19] that showed that the reason for the success of Lalanne’s formulation [17] was that it uniformly preserved the continuity of appropriate field components across the discontinuities of the permittivity function. In 1997, Li presented a new formulation [20] of the RCWA that applies correct rules [19] of Fourier factorization for crossed surface-relief gratings. RCWA with such a Fourier factorization converged much faster than traditional approaches. In 2003, Li introduced formulas for crossed anisotropic gratings featuring arbitrary permittivity and permeability tensors [28]. In 2008, Liu and Lalanne developed an RCWA algorithm to discuss 2D metallic gratings, focusing on extraordinary optical transmission (EOT) through subwavelength nanohole arrays [29]. In 2005, Logofatu developed a methodological framework that integrated analytical formulations with computational implementation for 2D gratings [30]. In recent years, several ideas for RCWA algorithms without solving eigenvalues have been progressively proposed [31,32]. While RCWA for single-layer gratings is well-established, its numerical instability in matching boundary conditions across multilayered grating interfaces has become a focal point of the research.

To tackle the computational instability and inefficiency for large-depth gratings, Moharam et al. [33] proposed the enhanced transmittance matrix (T-matrix) method for the stratified structures of 1D gratings. Compared to the standard T-matrix algorithm, the enhanced algorithm effectively reduces computational instability. Apart from the enhanced T-matrix method widely used in optics, other methods, including the reaction matrix (R-matrix) [34] and scattering matrix (S-matrix) [35], were put forward to enhance the computational stability of RCWA for stratified 1D gratings. In 2006, Tan presented enhanced R-matrix algorithms, which delivered superior stability compared to conventional R-matrix algorithms [36]. In 2011, Rumpf proposed an enhanced S-matrix algorithm by sandwiching a grating layer between two virtual thicknessless vacuum layers, accelerating the computational time by potential parallel processes [37]. Tan [38] pointed out that the T-matrix and R-matrix methods might exhibit sever computational instabilities if the thickness of the grating layer becomes excessively large or small, respectively. He proposed a so-called hybrid matrix (H-matrix) method for grating structures with an arbitrary thickness but with 1D periodicity. The H-matrix method demonstrated superior stability over the T-matrix and R-matrix methods [38].

In comparison with the 1D case, the computational stability problem of RCWA for structures with 2D periodicity is more severe considering the increased electromagnetic mode coupling complexity, but has received much less attention. In this work, the RCWA algorithms for the structures with 2D periodicity are derived step by step. Through systematic numerical benchmarking, we demonstrate the superior performance of the enhanced H-matrix relative to widely adopted large-scale computation matrices. The results confirm its optimal stability and computational efficiency as a boundary matching algorithm for reflection analysis. We used the design of metalens, for example, to demonstrate the application of such algorithms.

2. Formulation

2.1. Fourier Expansion of the Permittivity Within the Grating Region

A plane wave with wavelength λ is incident with an incidence angle θ, an azimuth angle ϕ, and a polarization angle ψ. The definition of the angles and the geometry of grating are shown in Figure 1a. The grating is bound by two semi-infinitive homogeneous media, Region I and Region II, with relative permittivity values of εI and εII, respectively. The relative permeability for the three regions is μ = 1. The 2D grating exhibits a spatial period of Λ_x along the x-axis and Λ_y along the y-axis. As Figure 1b shows, the grating has alternating regions of dielectric constants, ε_rd (ridge) and ε_gr (groove); the dielectric permittivity can be represented as [30]:

(1)$\epsilon \left(x,y\right)=\left\{\begin{array}{cc}{\epsilon }_{rd}& f_{x1}{\mathsf{\Lambda }}_x<x<f_{x2}{\mathsf{\Lambda }}_x \mathit{and} f_{y1}{\mathsf{\Lambda }}_y<y<f_{x2}{\mathsf{\Lambda }}_y\\ {\epsilon }_{gr}& else\end{array}\right.,$

where f_x = f_x2 − f_x1, f_y = f_y2 − f_y1. f_x1, f_x2 represent the dimensionless fractions of the two boundaries of a grating ridge with respect to a grating period along the x-axis. f_y1, and f_y2 have the same indications, but along the y-axis. The dielectric function in the grating region is first expanded using a Fourier series, which can be expressed as [30]:

(2)$\epsilon \left(x,y\right)=\sum _{m,n}{\epsilon }_{mn}\mathit{exp}\left[i\left(m{K}_{x}x+n{K}_{y}y\right)\right],$

where ε_mn is the (m, n)th Fourier coefficient of the relative permittivity. K_x and K_y are defined as K_x=2π/Λ_x, K_y=2π/Λ_y. The coefficients of the Fourier expansion (2) can be calculated by:

(3)${\epsilon }_{mn}={\displaystyle \frac{1}{{\mathsf{\Lambda }}_x{\mathsf{\Lambda }}_y}}{\int }_0^{{\mathsf{\Lambda }}_y} {\int }_0^{{\mathsf{\Lambda }}_x} \epsilon \left(x,y\right)\mathit{exp}\left[-i\left(m{K}_{x}x+n{K}_{y}y\right)\right]dxdy.$

According to Li [20], the algorithm demonstrated significantly accelerated convergence when utilizing the following Fourier decomposition.

(4)${⌈⌊\epsilon ⌋⌉}_{mn,ph}={⌈\{\lfloor 1/\epsilon {\rfloor }^{-1}{\}}_{nh}⌉}_{mp}={\displaystyle \frac{1}{{\mathsf{\Lambda }}_x}}{\int }_0^{{\mathsf{\Lambda }}_x} \{\lfloor 1/\epsilon {\rfloor }^{-1}{\}}_{nh}(x\left)exp\right[-i(m-p)K_{x}x]dx,$

(5)${⌊⌈\epsilon ⌉⌋}_{mn,ph}=\lfloor \{{⌈1/\epsilon ⌉}^{-1}{\}}_{mp}{\rfloor }_{nh}={\displaystyle \frac{1}{{\mathsf{\Lambda }}_y}}{\int }_0^{{\mathsf{\Lambda }}_y} \{{⌈1/\epsilon ⌉}^{-1}{\}}_{mp}\left(y\right)exp[-i(n-h\left)K_{y}y\right]dy,$

where:

(6)${⌈1/\epsilon ⌉}_{mn}={\displaystyle \frac{1}{{\mathsf{\Lambda }}_x}}{\int }_0^{{\mathsf{\Lambda }}_x} {\displaystyle \frac{1}{\epsilon (x,y)}}\mathit{exp}\left[-i\left(m-n\right)K_{x}x\right]dx,$

(7)${⌊1/\epsilon ⌋}_{mn}={\displaystyle \frac{1}{{\mathsf{\Lambda }}_y}}{\int }_0^{{\mathsf{\Lambda }}_y} {\displaystyle \frac{1}{\epsilon (x,y)}}\mathit{exp}\left[-i\left(m-n\right)K_{y}y\right]dy.$

2.2. Electromagnetic Fields at Region I and Region II

The incident electric field can be expressed as [28]:

(8)${\overrightarrow{E}}_{inc,\sigma }=I_{\sigma }\overrightarrow{\sigma }\mathrm{e}\mathrm{x}\mathrm{p}\left[i\left(k_x^{0}x+k_y^{0}y-k_z^{00}z\right)\right],$

where σ = x, y, or z, representing the electric component along the x-, y-, or z- direction. $I_{\sigma }$ is the amplitude of the incident electric field. The $k_x^0$, $k_y^0$, $k_z^{00}$ can be defined as [30]:

(9)$\begin{array}{ccc}{k}_x^0=k_{\mathrm{I}}sin\theta cos\varphi ,& k_y^0=k_{\mathrm{I}}sin\theta sin\varphi ,& k_z^{00}=k_{\mathrm{I}}cos\theta \end{array}.$

where $k_{\mathrm{I}}=2\pi \sqrt{{\epsilon }_{\mathrm{I}}\mu }/\lambda$ is defined as the wave number in Region I. $k_{\mathrm{I}\mathrm{I}}=2\pi \sqrt{{\epsilon }_{\mathrm{I}\mathrm{I}}\mu }/\lambda$ is defined as the wave number in Region II. Then, electric fields in Region I and Region II are expressed with the following Rayleigh expansions [20]:

(10)${\overrightarrow{E}}_{\sigma ,\mathrm{I}}={\overrightarrow{E}}_{inc,\sigma }+\sum _{m,n}{R}_{\sigma }^{m,n}\overrightarrow{\sigma }\exp \left[i\left(k_x^{m}x+k_y^{n}y+k_{z,\mathrm{I}}^{mn}z\right)\right],$

(11)${\overrightarrow{E}}_{\sigma ,\mathrm{I}\mathrm{I}}=\sum _{m,n}{T}_{\sigma }^{m,n}\overrightarrow{\sigma }\exp \left[i\left(k_x^{m}x+k_y^{n}y-k_{z,\mathrm{I}\mathrm{I}}^{mn}z\right)\right].$

where σ = x or y. $R_{\sigma }^{m,n}$ and $T_{\sigma }^{m,n}$ are the amplitudes of the (m, n)th-order reflected and transmitted electric fields, respectively. $k_x^m$ and $k_y^n$ are defined as:

(12)$\begin{array}{cc}{k}_x^m=k_x^0+m{K}_x,& k_y^n=k_y^0+n{K}_y\end{array}.$

$k_{z,a}^{mn}$ is to be solved from:

(13)${k_x^m}^2+{k_y^n}^2+{k_{z,a}^{mn}}^2=k_a^2,$

where a = I or II, which represents the wave number in Region $\mathrm{I}$ or Region $\mathrm{I}\mathrm{I}$. The sign of $k_{z,a}^{mn}$ is chosen according to [28]:

$Re\left[k_{z,a}^{mn}\right]+Im\left[k_{z,a}^{mn}\right]>0.$

The Maxwell’s equations, when utilizing Gaussian units, are given by the following equations [28]:

(14)$\begin{array}{cc}\nabla \times \overrightarrow{\mathrm{E}}=i{k}_0\mu \overrightarrow{H},& \nabla \times \overrightarrow{\mathrm{H}}=-i{k}_0\epsilon \overrightarrow{\mathrm{E}}\end{array},$

where k₀ = 2π/λ is the wave number in the vacuum. In Region I and Region II, by substituting Equations (10) and (11) into Maxwell’s Equation (14), the amplitudes of electromagnetic fields in Region II and Region I can be written as [28]:

(15)$\left[\begin{array}{c}{E}_{x,a}^{mn}\\ E_{y,a}^{mn}\\ \begin{array}{c}{H}_{x,a}^{mn}\\ H_{y,a}^{mn}\end{array}\end{array}\right]=\left[\begin{array}{cccc}I& 0& I& 0\\ 0& I& 0& I\\ -A_a^{mn}& -B_a^{mn}& A_a^{mn}& B_a^{mn}\\ C_a^{mn}& A_a^{mn}& {-C}_a^{mn}& -A_a^{mn}\end{array}\right]\left[\begin{array}{cccc}{X}_a^{mn+}& 0& 0& 0\\ 0& X_a^{mn+}& 0& 0\\ 0& 0& X_a^{mn-}& 0\\ 0& 0& 0& X_a^{mn-}\end{array}\right]\left[\begin{array}{c}{c_{x,a}^{mn}}^{+}\\ {c_{y,a}^{mn}}^{+}\\ {c_{x,a}^{mn}}^{-}\\ {c_{y,a}^{mn}}^{-}\end{array}\right],$

where $A_a^{m,n}$, $B_a^{m,n}$, $C_a^{m,n}$ and $X_a^{mn}$ are defined as [28]:

(16)$\begin{array}{ccc}{A}_a^{m,n}={\displaystyle \frac{{k_y^{n}k}_x^m}{k_{z,a}^{mn}\mu k_0}},& B_a^{m,n}={\displaystyle \frac{{-{\left(k_x^m\right)}^2+k}_a^2}{k_{z,a}^{mn}\mu k_0}},& C_a^{m,n}={\displaystyle \frac{{-{\left(k_y^n\right)}^2+k}_a^2}{k_{z,a}^{mn}\mu k_0}},\end{array}$

(17)$\begin{array}{cc}{X_a^{mn}}^{-}=exp\left(-i{k}_{z,a}^{mn}z\right),{X_a^{mn}}^{+}=exp\left(i{k}_{z,a}^{mn}z\right)& ,\end{array}$

(18)$\begin{array}{cccc}{c_{\sigma ,\mathrm{I}}^{mn}}^{+}=R_{\sigma }^{mn}& {c_{\sigma ,\mathrm{I}\mathrm{I}}^{mn}}^{+}=0& {c_{\sigma ,\mathrm{I}}^{mn}}^{-}=I_{\sigma }& {c_{\sigma ,\mathrm{I}\mathrm{I}}^{mn}}^{-}=T_{\sigma }^{mn}\end{array},$

(19)$\begin{array}{cc}{W}_a=\left[\begin{array}{cc}I& 0\\ 0& I\end{array}\right],& V_a=\left[\begin{array}{cc}-A_a^{mn}& -B_a^{mn}\\ C_a^{mn}& A_a^{mn}\end{array}\right]\end{array}.$

The reason for writing the electromagnetic fields in the form of Equation (15) is that the field components share a similar configuration for the three regions, as will be discussed in the subsequent sections in this paper.

2.3. Electromagnetic Fields at the Grating Region

The electromagnetic fields can be expressed using a Fourier expansion in the grating region, which can be written as [30]:

(20)${\overrightarrow{E}}_g=\sum _{m,n}\left[E_x^{mn}\left(z\right)\overrightarrow{x}+E_y^{mn}\left(z\right)\overrightarrow{y}+E_z^{mn}\left(z\right)\overrightarrow{z}\right]exp\left(i{k}_x^{m}x+i{k}_y^{n}y\right),$

(21)${\overrightarrow{H}}_g=\sum _{m,n}\left[H_x^{mn}\left(z\right)\overrightarrow{x}+H_y^{mn}\left(z\right)\overrightarrow{y}+H_z^{mn}\left(z\right)\overrightarrow{z}\right]exp\left(i{k}_x^{m}x+i{k}_y^{n}y\right),$

where $E_x^{mn}$, $H_x^{mn}$ $E_y^{mn}$, $H_y^{mn},$ $E_z^{mn}$, and $H_z^{mn}$ are the complex amplitudes of the electromagnetic field. By substituting Equations (20) and (21) into Maxwell’s Equation (14), the spatial harmonic components within the grating region are derived [30]:

(22)${\displaystyle \frac{k_0}{i}}{\displaystyle \frac{\partial E_x^{m,n}}{\partial z}}=k_0^2\mu H_y^{m,n}-k_x^m\sum _{p,q}{\displaystyle \frac{1}{{\epsilon }_{m-p,n-h}}}\left(k_x^p{H}_y^{p,h}-k_y^h{H}_x^{p,h}\right),$

(23)${\displaystyle \frac{k_0}{i}}{\displaystyle \frac{\partial E_y^{m,n}}{\partial z}}={-k}_0^2\mu H_x^{m,n}-k_y^n\sum _{p,q}{\displaystyle \frac{1}{{\epsilon }_{m-p,n-h}}}\left(k_x^p{H}_y^{p,h}-k_y^h{H}_x^{p,h}\right),$

(24)${\displaystyle \frac{{\mu k}_0}{i}}{\displaystyle \frac{\partial H_y^{m,n}}{\partial z}}=\mu k_0^2\sum _{p,q}{\epsilon }_{m-p,n-h}{E}_x^{p,h}{+k}_y^n\left(k_x^m{E}_y^{m,n}-k_y^n{E}_x^{m,n}\right)$

(25)${\displaystyle \frac{{\mu k}_0}{i}}{\displaystyle \frac{\partial H_x^{m,n}}{\partial z}}=-\mu k_0^2\sum _{p,q}{\epsilon }_{m-p,n-h}{E}_y^{p,h}+k_x^m\left(k_x^m{E}_y^{m,n}-k_y^n{E}_x^{m,n}\right).$

To unify the dimensionality of all the variables, the indices for the 2D grating are redefined as:

(26)$\begin{array}{cc}u=m{s}_1+n,& v=p{s}_2+h\end{array},$

where s₁ = 2n_max + 1 and s₂ = 2h_max + 1. Here, m ranges from −m_max to m_max, n ranges from–n_max to n_max, p ranges from −p_max to p_max, and h ranges from −h_max to h_max, with m_max = p_max and n_max = h_max. The parameters m_max and p_max represent the maximum order of Fourier harmonics in the x-axis. Similarly, n_max and h_max represent the maximum order of Fourier harmonics in the y-axis. Thus, the index u ranges from −m_max (2n_max + 1) − n_max to m_max (2n_max + 1) +n_max, and v ranges from −p_max (2h_max + 1) − h_max to p_max (2h_max + 1) + h_max. The total number of possible values for u is (2m_max + 1) (2n_max + 1), and for v, it is (2p_max + 1) (2h_max + 1). When m_max = p_max = 0(n_max = h_max = 0) and f_x = 1(f_y = 1), the 2D grating turns into a 1D grating with the periodicity in the y- (x-) axis. The dimensionality reduction from 4D to 2D is achieved by Ξ^uv = ε_m-p,n-h, A^uv = 1/ε_m-p,n-h, ${\Xi }_{XY}^{uv}$ =${⌊⌈\epsilon ⌉⌋}_{mn,ph}$, ${\Xi }_{YX}^{uv}$=${⌈\lfloor \epsilon \rfloor ⌉}_{mn,ph}$. The wave vectors are converted to 2D vectors via the following diagonal matrix construction: $K_x^{uv}=k_x^m{\delta }_{uv}/k_0$, $K_y^{uv}=k_y^n{\delta }_{uv}/k_0$. δ_uv is the Kronecker delta function, which is zero everywhere, except u = v.

After unifying the dimensionality of all the variables, Equations (22)–(25) can be reduced to:

(27)${\displaystyle \frac{k_0}{i}}{\partial }_z\left(\begin{array}{c}{E}_x\\ E_y\end{array}\right)=F\left(\begin{array}{c}{H}_x\\ H_y\end{array}\right), {\displaystyle \frac{k_0}{i}}{\partial }_z\left(\begin{array}{c}{H}_x\\ H_y\end{array}\right)=G\left(\begin{array}{c}{E}_x\\ E_y\end{array}\right).$

According to Li [20], the convergence is significantly improved when $F$ and $G$ are defined as:

(28)$F=\left[\begin{array}{cc}{K}_x^{uv}{A}^{uv}{K}_y^{uv}& \mu k_0^2-K_x^{uv}{A}^{uv}{K}_x^{uv}\\ -\mu k_0^2+K_y^{uv}{A}^{uv}{K}_y^{uv}& {-K}_y^{uv}{A}^{uv}{K}_y^{uv}\end{array}\right],$

(29)$G=\left[\begin{array}{cc}-K_x^{uv}{K}_y^{uv}& {-\mu k_0^2{\Xi }_{YX}^{uv}+K}_x^{uv}{K}_x^{uv}\\ -K_y^{uv}{K}_y^{uv}+\mu k_0^2{\Xi }_{XY}^{uv}& K_x^{uv}{K}_y^{uv}\end{array}\right].$

Equation (27) can be further reduced to:

(30)${-k}_0^2{\partial }_z^2\left(\begin{array}{c}{E}_x\\ E_y\end{array}\right)=M\left(\begin{array}{c}{E}_x\\ E_y\end{array}\right)$

where M = F*G. The coupled-wave Equation (30) can be addressed by determining the eigenvectors and the eigenvalues associated with the matrix M. Following this, the amplitudes of the electric and magnetic fields in the grating region can be expressed as [28]:

(31)$E_{g,\sigma }\left(z\right)=\sum _{u,v,j}{W}_{\sigma uvj}\left[c_{\sigma uj}^{-}exp\left(-i{q}_{uvj}z\right)+c_{\sigma uj}^{+}exp\left(i{q}_{uvj}z\right)+\right]\mathit{exp}\left[i\left(K_X^{uv}x+K_Y^{uv}y\right)\right],$

(32)$H_{g,\sigma }\left(z\right)=\sum _{u,v,j}{V}_{\sigma uvj}\left[-c_{\sigma uj}^{-}exp\left(-i{q}_{uvj}z\right)+c_{\sigma uj}^{+}exp\left(i{q}_{uvj}z\right)\right]\mathit{exp}\left[i\left(K_X^{uv}x+K_Y^{uv}y\right)\right].$

Q and W are the eigenvalue and eigenvector matrices of the matrix M, respectively. q_uvj and W_σuvj are the elements of matrices $\sqrt{Q/k_0^2}$ and W, respectively. The sign of the q_uvj should be satisfied to [20]:

(33)$Re\left[q_{uvi}\right]+Im\left[q_{uvi}\right]>0,$

V = GW/(k₀Q). $c_{\sigma uj}^{+}$ and $c_{\sigma uj}^{-}$ can be determined by matching boundary conditions. Equations (31) and (32) can be expressed as a matrix form:

(34)$\left[\begin{array}{c}{E}_{g,xu}\\ E_{g,yu}\\ H_{g,xu}\\ H_{g,yu}\end{array}\right]=\left[\begin{array}{cc}{W}_{xuvj}& W_{xuvj}\\ W_{yuvj}& W_{yuvj}\\ V_{xuvj}& {-V}_{xuvj}\\ V_{yuvj}& -V_{yuvj}\end{array}\right]\left[\begin{array}{cccc}{X}_{uvj+}^{+}& 0& 0& 0\\ 0& X_{uvj-}^{+}& 0& 0\\ 0& 0& X_{uvj+}^{-}& 0\\ 0& 0& 0& X_{uvj-}^{-}\end{array}\right]\left[\begin{array}{c}{c}_{xuj}^{+}\\ c_{yuj}^{+}\\ c_{xuj}^{-}\\ c_{yuj}^{-}\end{array}\right].$

Equation (34) can be further simplified as:

(35)$\left[\begin{array}{c}{E}_g\\ H_g\end{array}\right]=\left[\begin{array}{cc}W& W\\ V& -V\end{array}\right]\left[\begin{array}{cc}{X}^{+}& 0\\ 0& X^{-}\end{array}\right]\left[\begin{array}{c}{c}^{+}\\ c^{-}\end{array}\right],$

where

$\begin{array}{cccc}{E}_g=\left[\begin{array}{c}{E}_{g,xu}\\ E_{g,yu}\end{array}\right],& H_g=\left[\begin{array}{c}{H}_{g,xu}\\ H_{g,yu}\end{array}\right],& W=\left[\begin{array}{c}{W}_{xuvj}\\ W_{yuvj}\end{array}\right],& V=\end{array}\left[\begin{array}{c}{V}_{xuvj}\\ V_{yuvj}\end{array}\right],$

$\begin{array}{cccc}{X}^{-}=\left[\begin{array}{cc}{X}_{uvj+}^{-}& 0\\ 0& X_{uvj-}^{-}\end{array}\right],& X^{+}=\left[\begin{array}{cc}{X}_{uvj+}^{+}& 0\\ 0& X_{uvj-}^{+}\end{array}\right],& c^{-}=\left[\begin{array}{c}{c}_{xuj}^{-}\\ c_{yuj}^{-}\end{array}\right],& c^{+}=\left[\begin{array}{c}{c}_{xuj}^{+}\\ c_{yuj}^{+}\end{array}\right]\end{array}.$

Employing an identical transformation, Equation (15) can be recast in an equivalent structure. Equation (35) is analogous to the RCWA formulas for 1D gratings, where W and V are the eigenvector matrix blocks. X is the propagation matrix block. c⁺ and c⁻ represent unknown constants in different regions, which are determined by applying boundary condition matching.

2.4. Matching the Boundary Conditions in Single-Layer 2D Gratings

To determine the unknown constants, we first considered matching the boundary conditions for the single layer of gratings. The relationship of electromagnetic fields between z = 0 and z = d can be expressed in an H-matrix form [38]:

(36)$\left[\begin{array}{c}{E}_g\left(0\right)\\ H_g\left(d\right)\end{array}\right]=\left[\begin{array}{cc}{H}_{11}^{\left[1\right]}& H_{12}^{\left[1\right]}\\ H_{21}^{\left[1\right]}& H_{22}^{\left[1\right]}\end{array}\right]\left[\begin{array}{c}{H}_g\left(0\right)\\ E_g\left(d\right)\end{array}\right],$

where matrix blocks H^[1] can be obtained from:

(37)$\left[\begin{array}{cc}{H}_{11}^{\left[1\right]}& H_{12}^{\left[1\right]}\\ H_{21}^{\left[1\right]}& H_{22}^{\left[1\right]}\end{array}\right]=\left[\begin{array}{cc}W& W{X}^{+}\left(d\right)\\ V{X}^{+}\left(d\right)& -V\end{array}\right]{\left[\begin{array}{cc}V& -V{X}^{+}\left(d\right)\\ W{X}^{+}\left(d\right)& W\end{array}\right]}^{-1}.$

W and V are the eigenvector matrix blocks in the single grating layer. X is the propagation matrix in the single grating layer. Now, we match the tangential electric and magnetic fields through [38]:

(38)$\begin{array}{cc}{E}_{\sigma ,\mathrm{I}\mathrm{I}}\left(d\right)=E_{\sigma ,g}\left(d\right)=W_{\mathrm{I}\mathrm{I}}T,& H_{\sigma ,\mathrm{I}\mathrm{I}}\left(d\right)=H_{\sigma ,g}\left(d\right)=V_{\mathrm{I}\mathrm{I}}T,\\ E_{\sigma ,\mathrm{I}}\left(0\right)=E_{\sigma ,g}\left(0\right)=W_{\mathrm{I}}{I}_{\sigma }+W_{\mathrm{I}}R,& H_{\sigma ,\mathrm{I}}\left(0\right)=H_{\sigma ,g}\left(0\right)=V_{\mathrm{I}}{I}_{\sigma }-V_{\mathrm{I}}R.\end{array}$

For most grating problems, transmission and reflection amplitudes are expected, which can be obtained by [38]:

(39)$\left[\begin{array}{c}R\\ T\end{array}\right]={\left[\begin{array}{cc}-W_{\mathrm{I}}{-H}_{11}^{\left[1\right]}{V}_{\mathrm{I}}& H_{12}^{\left[1\right]}{W}_{\mathrm{I}\mathrm{I}}\\ -H_{21}^{\left[1\right]}{V}_{\mathrm{I}}& {-V_{\mathrm{I}\mathrm{I}}+H}_{22}^{\left[1\right]}{W}_{\mathrm{I}\mathrm{I}}\end{array}\right]}^{-1}\left[\begin{array}{c}{-H_{11}^{\left[1\right]}{V}_{\mathrm{I}}+W}_{\mathrm{I}}\\ -H_{21}^{\left[1\right]}{V}_{\mathrm{I}}\end{array}\right]\left[I_{\sigma }\right].$

Once the diffraction amplitudes are obtained from Equation(39), the diffraction efficiencies can be determined by using [28]:

(40)$D{E}_{ri}={B_{\mathrm{I}}^{uv}\left|R_x^u\right|}^2+{A_{\mathrm{I}}^{uv}\left|R_y^u\right|}^2+C_{\mathrm{I}}^{uv}\left(R_y^u\overline{R_x^u}+R_x^u\overline{R_y^u}\right),$

(41)$D{E}_{ti}=B_{\mathrm{I}\mathrm{I}}^{uv}{\left|T_x^u\right|}^2+A_{\mathrm{I}\mathrm{I}}^{uv}{\left|T_y^u\right|}^2+C_{\mathrm{I}\mathrm{I}}^{uv}\left(T_y^u\overline{T_x^u}+T_x^u\overline{T_y^u}\right),$

where $A_a^{uv}$, $B_a^{uv}$ and $C_a^{uv}$ are can be derived by through Equations (16) and (26). $R_x^u$, $R_y^u$, $T_x^u$ and $T_y^u$ are defined as:

$\begin{array}{cc}\left[\begin{array}{c}{R}_x^u\\ R_y^u\end{array}\right]=R,& \left[\begin{array}{c}{T}_x^u\\ T_y^u\end{array}\right]=T.\end{array}$

2.5. Matching the Boundary Conditions in Stratified 2D Gratings

For stratified 2D gratings, as shown in Figure 2, the grating region is further divided into L layers. The refraction index of the lth layer is n_l (l = 1, 2, …, L) and the thickness of the lth layer is d_l. Each layer may have distinct a grating period or filling factor. There are connections between the electromagnetic field at the boundaries z = Z_l⁻ and z = Z_L⁺, which can be described using a stacked H^[l] matrix:

(42)$\left[\begin{array}{c}{E}_l\left(Z_l^{-}\right)\\ H_L\left(Z_L^{+}\right)\end{array}\right]=\left[\begin{array}{cc}{H}_{11}^{\left[l\right]}& H_{12}^{\left[l\right]}\\ H_{21}^{\left[l\right]}& H_{22}^{\left[l\right]}\end{array}\right]\left[\begin{array}{c}{H}_l\left(Z_l^{-}\right)\\ E_L\left(Z_L^{+}\right)\end{array}\right].$

The relationship between stack H^[l] and stack H^[l+1] is [38]:

(43)$\begin{array}{c}\left[\begin{array}{cc}{H}_{11}^{\left[l\right]}& H_{12}^{\left[l\right]}\\ H_{21}^{\left[l\right]}& H_{22}^{\left[l\right]}-H_{22}^{\left[l+1\right]}\end{array}\right]=\left[\begin{array}{cc}{W}_l& W_l{X}_l^{+}\left(d_l\right)\\ H_{21}^{\left[l+1\right]}{V}_l{X}_l^{+}\left(d_l\right)& {-H}_{21}^{\left[l+1\right]}{V}_l\end{array}\right]\left[\begin{array}{c}{V}_l\\ \left[W_l-H_{11}^{\left[l+1\right]}{V}_l\right]X_l^{+}\left(d_l\right)\end{array}\right.\\ {\left.\begin{array}{c}{-V}_l{X}_l^{+}\left(d_l\right)\\ W_l+H_{11}^{\left[l+1\right]}{V}_l\end{array}\right]}^{-1}\left[\begin{array}{cc}I& 0\\ 0& H_{12}^{\left[l+1\right]}\end{array}\right].\end{array}$

H^[L+1] is the H-matrix of the output layer, which is:

(44)$\left[\begin{array}{cc}{H}_{11}^{\left[L+1\right]}& H_{12}^{\left[L+1\right]}\\ H_{21}^{\left[L+1\right]}& H_{22}^{[L+1]}\end{array}\right]=\left[\begin{array}{cc}0& I\\ I& 0\end{array}\right].$

By substituting Equation (44) into Equation (43), the H^[L] matrix at layer L can be obtained from the H^[L+1] matrix at layer L+1. Substituting H^[L] into Equation (43) yields H^[L−1]. By repeating this procedure, H^[1] is ultimately acquired. H^[1] provides the connection of electromagnetic fields between the first layer and L layer. After obtaining H^[1], for stratified grating structures, the transmission and reflection amplitudes, as well as the diffraction amplitudes, can be calculated using Equations (39)–(41).

For problems solely concerning reflectance, the enhanced H matrix is expressed by [38]:

(45)$H_{11}^{\left[l\right]}=\left[W_l+W_l{\mathsf{\Gamma }}_l\right]{\left[V_l-V_l{\mathsf{\Gamma }}_l\right]}^{-1},$

where:

(46)${\mathsf{\Gamma }}_l=-X_l^{-}\left(-d_l\right){\left[W_l+H_{11}^{\left[l+1\right]}{V}_l\right]}^{-1}\left[W_l-H_{11}^{\left[l+1\right]}{V}_l\right]X_l^{+}\left(d_l\right).$

We assume the $L+1$ layer is infinite, so $H_{11}^{\left[L+1\right]}$ is determined by:

(47)$H_{11}^{\left[\mathrm{L}+1\right]}={W_{\mathrm{L}+1}\left(V_{\mathrm{L}+1}\right)}^{-1}.$

The reflection efficiency is given by [38]:

(48)$R={\left[H_{11}^{\left[1\right]}V+W\right]}^{-1}\left[{-W+H}_{11}^{\left[1\right]}V\right]I_{\sigma }{\delta }_{m0}{\delta }_{n0}.$

The enhanced H matrix exhibits greater stability than other algorithms. A detailed discussion is provided in the next section.

3. Numerical Result

For numerical validations, we first verify the stability of the enhanced H-matrix method using a symmetric single-layer 2D grating. The incident wavelength λ is 1 μm. The incident angle, azimuth angle, and polarization angle are set as θ = ϕ = ψ = 0°. The key parameters of the grating are listed in Table 1.

Figure 3 compares computational stability as a function of grating thickness d between five algorithms, the T-matrix (circles), S-matrix (dashed line), R-matrix (square), H-matrix (triangular), and enhanced H-matrix (solid line), in a logarithmic coordinate system. The condition number of a matrix is the ratio between the largest singular value and the smallest singular value, implying the sensitivity of output values to the input data. A larger condition number indicates greater matrix instability and a higher likelihood of numerical errors during matrix inversion. We calculated the condition numbers of the matrices to be inverted in these methods. For instance, the condition number of the matrix to be inverted in H^[1] corresponds to the condition number of the second matrix on the right-hand side of Equation (43). We calculated the condition numbers in T^[1] (Equation (16) in Ref. [38]), R^[1] (Equation (17) in Ref. [36]), S^[1] (Equation (19) in Ref. [37]), H^[1] (Equation (32)), and enhanced H^[1] (Equation (45)) versus the layer thickness relative to the wavelength.

As Figure 3 shows, the condition numbers in R[1] decrease as the grating thickness increases in the bilogarithmic diagram, indicating the R-matrix method is suitable for the thick gratings. By contrast, the condition numbers in T[1] increase almost exponentially as the grating thickness increases, meaning that the T-matrix method is suitable for the thin gratings. However, the H-matrix and S-matrix combine the advantages of the R-matrix method and T-matrix method. The condition numbers of H^[1] and S^[1] to be inverted initiate from a much lower level than R and remain quite stable for a large thickness range. Notably, condition numbers in enhanced H^[1] to be inverted demonstrate the smallest value compared to the other methods, which demonstrates unconditional computational stability, even for 2D periodic structures.

To evaluate computational efficiency, an asymmetric single-layer 2D grating was used for comparative analysis. Table 2 lists the details of the testing platform. It can be noted that all the benchmarks are evident on the same platform. The incident wavelength λ is 1 μm. The angles of incidence, azimuth, and polarization are all set to θ = ϕ = ψ = 0°. Parameters are detailed in Table 3. Harmonic orders—representing the maximum Fourier expansion terms in RCWA simulations—directly influence solution accuracy. While higher harmonics improve precision, they increase computational demands. We employed both the enhanced H-matrix and the S-matrix method to compute reflection coefficients for this structure, with the computation time serving as the efficiency metric. The runtime measurement method uses MATLAB R2022a’s tic and toc functions, while memory usage is measured with MATLAB’s memory function. These functions collect the memory usage and time consumption of the S-matrix and H-matrix code blocks, while the memory usage and time consumption of other code blocks, such as eigenvector calculations, are not included. As Figure 4a shows, the enhanced H-matrix method consistently outperforms the S-matrix approach, reducing the computation time by a range of 16–30%. According to Figure 4b, the enhanced H-matrix outperforms the S-matrix, with a memory consumption value reduced by approximately 20%. This indicates that the enhanced H-matrix shows advantages in optimizing memory usage as well as in computational efficiency. It should be noted that during the computation of the S-matrix, transmission-related variables were intentionally excluded—only reflection variables were computed—ensuring a valid comparison.

We next verify the accuracy of the H-matrix using three-layer grating slabs, as shown in Figure 5a. The light is incident from Region II. The incident wavelength is 1.55 µm, with incident angles of θ = 0°, ϕ = 0°, and ψ = 0°. The grating parameters are listed in Table 4. Figure 5b shows the simulation results using the S-matrix method (cube and triangle symbols) and H-matrix method (cross and dashed line). The x-axis represents the first layer thickness relative to the wavelength, and the y-axis shows reflection and transmission efficiencies. When the normalized parameter d₁/λ is relatively small, variations in d₁ are minimal; consequently, diffraction efficiency shows only minor changes. Conversely, when d₁/λ is sufficiently large, variations in d₁ become substantial, resulting in significant alterations in the diffraction efficiency. As can be seen, there is a good agreement between the H-matrix method and S-matrix method, validating the accuracy of our proposed H-matrix algorithm. In Figure 5b, the solid line represents the results of the enhanced H-matrix algorithm. The enhanced H-matrix method provides almost the same results as those of the regular H-matrix method and S-matrix method, but the required computation memory is reduced by almost half [38].

We finally demonstrate the practical application of such an H-matrix-enhanced 2D RCWA algorithm. To validate the algorithm, a metalens consisting of 2D gratings was designed. The incident light with a wavelength of 532 nm is normally incident at θ = ϕ = ψ = 0° on the metalens. The unit cell of the metalens comprises a substrate layer with a thickness of 1 μm and a refractive index of 1.5. The nanopillar within the unit cell of the metalens has a refractive index of 2.04 and a height of 1.2 μm. The periodicity of the unit structure is 0.45 μm. As shown in Figure 6a, a parametric database was established utilizing the 2D RCWA H-matrix algorithm. The database established a map between P_w and the phase of transmission. P_w denotes the cross-sectional width of the cube-shaped nanopillar. According Capasso [39], the metalens unit cells are spatially arranged according to hyperbolic phase requirements. The metalens is designed with a focal length of 100 μm and a radius of 11 μm. An approximate method was employed to compute electromagnetic fields in the far-field beyond the metalens [3]. The distribution of electric field intensity in the x-z plane is shown in Figure 6b. Figure 6c shows the electric field amplitude along the x-direction in the focal plane. The full-width at half maximum of the focus is 2.35 μm. The results demonstrate that the metalens designed using the 2D RCWA H-matrix algorithm achieves excellent phase matching, thereby validating the accuracy of our algorithm.

4. Conclusions

In this paper, we derive the RCWA algorithms for 2D gratings step by step. We also develop a stable enhanced H-matrix algorithm for 2D stratified grating structures. The simulation results of specific grating structures demonstrate that the enhanced H-matrix algorithm offers better stability compared to the T-matrix, R-matrix and S-matrix algorithms. The integration of the enhanced H-matrix with 2D grating structures provides a more stable and efficient approach. Notably, the dimensionality reduction achieved through the enhanced matrix methodology offers a new solution to address the prolonged runtime and excessive memory consumption in the inverse designs of metasurfaces, which will be an interesting research topic in our future study.

Author Contributions

K.S.: Conceptualization, Data Curation, Formal Analysis, Investigation, Methodology, Software, Validation, Visualization, Writing—Original Draft Preparation, and Writing—Review and Editing; J.W.: Conceptualization, Data Curation, Formal Analysis, Investigation, Methodology, Resources, Validation, Writing—Original Draft Preparation, Writing—Review and Editing, Funding Acquisition, and Supervision; G.W.: Methodology, Validation, Writing—Review and Editing, and Supervision. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 (a) Schematics of single 2D-layer grating structures. (b) Schematic diagram of unit 2D grating structures.

Figure 2 Schematics of a stratified unit grating cell.

Figure 3 Comparison of the condition numbers between five recursive methods.

Figure 4 (a) Comparison of the computational speed between the S-matrix and enhanced H-matrix. (b) Comparison of the memory usage between the S-matrix and enhanced H-matrix.

Figure 5 (a) Schematic unit cell of a three-layer grating slabs. (b) Validation of H-matrix computational accuracy against the S-matrix method.

Figure 6 The metalens is designed using the 2D RCWA H-matrix algorithm. (a) The database establishing a map between P_w and the phase of transmission. (b) Simulated electric field in the x-z cross-section. (c) Electric field intensity in the focal plane.

Table 1

The parameter of the symmetric single-layer grating.

Parameter	${\mathsf{\Lambda }}_{\mathit{x}}$	${\mathsf{\Lambda }}_{\mathit{y}}$	${\mathit{f}}_{\mathit{x}1}$	${\mathit{f}}_{\mathit{x}2}$	${\mathit{f}}_{\mathit{y}1}$	${\mathit{f}}_{\mathit{y}2}$	${\mathit{\epsilon }}_{\mathit{I}}$	${\mathit{\epsilon }}_{\mathit{I}\mathit{I}}$	${\mathit{\epsilon }}_{\mathit{r}\mathit{d}}$	${\mathit{\epsilon }}_{\mathit{g}\mathit{r}}$	$\mathit{m}$	$\mathit{n}$
Value	$1.2 \mathsf{\mu }\mathrm{m}$	$1.2 \mathsf{\mu }\mathrm{m}$	0.25	0.75	0.25	0.75	1	2.25	2.25	1	12	12

Table 2

The detailed implementation and testing setup information.

Category	Specification
Processor	HexaCore AMD Ryzen 5 4600H, 3921 MHz (3.0 GHz base, 4.0 GHz turbo) 6 cores, 12 threads
Cache	384 KB L1 cache, 3 MB L2 cache, 8 MB L3 cache
Memory	16GB (2 × 8 GB) Ramaxel RMSA3320MJ78HAF-3200 DDR4
Memory Transfer Speed	Read: [45.68 GB/s], Write: [45.80 GB/s]
Operating System	Windows 10 64-bit
Software	MATLAB R2022a, default numerical libraries

Table 3

The parameter of the asymmetric single-layer grating.

Parameter	${\mathsf{\Lambda }}_{\mathit{x}}$	${\mathsf{\Lambda }}_{\mathit{y}}$	${\mathit{f}}_{\mathit{x}1}$	${\mathit{f}}_{\mathit{x}2}$	${\mathit{f}}_{\mathit{y}1}$	${\mathit{f}}_{\mathit{y}2}$	${\mathit{\epsilon }}_{\mathit{I}}$	${\mathit{\epsilon }}_{\mathit{I}\mathit{I}}$	${\mathit{\epsilon }}_{\mathit{r}\mathit{d}}$	${\mathit{\epsilon }}_{\mathit{g}\mathit{r}}$	$\mathit{d}$
Value	$1.2 \mathsf{\mu }\mathrm{m}$	$1.2 \mathsf{\mu }\mathrm{m}$	0.416	0.833	0.416	0.833	1	2.25	2.25	1	$1 \mathsf{\mu }\mathrm{m}$

Table 4

Three-layer grating slabs.

Parameter	Value	Parameter	Value
d2	0.5 μm	${\epsilon }_{II}$	1
d3	0.5 μm	${\epsilon }_{layer1,rd}$	2.16
${\mathsf{\Lambda }}_x$	0.6 μm	${\epsilon }_{layer1,gr}$	1
${\mathsf{\Lambda }}_y$	0.6 μm	${\epsilon }_{layer2}$	12.06
${\epsilon }_I$	1	${\epsilon }_{layer3}$	2.16
$m_{max}$	10	$n_{max}$	12

References

1. Yu, N.; Capasso, F. Flat optics with designer metasurfaces. Nat. Mater.; 2014; 13, pp. 139-150. [DOI: https://dx.doi.org/10.1038/nmat3839]

2. Yu, N.; Genevet, P.; Kats, M.A.; Aieta, F.; Tetienne, J.-P.; Capasso, F.; Gaburro, Z. Light propagation with phase discontinuities: Generalized laws of reflection and refraction. Science; 2011; 334, pp. 333-337. [DOI: https://dx.doi.org/10.1126/science.1210713] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21885733]

3. Pestourie, R.; Pérez-Arancibia, C.; Lin, Z.; Shin, W.; Capasso, F.; Johnson, S.G. Inverse design of large-area metasurfaces. Opt. Express; 2018; 26, pp. 33732-33747. [DOI: https://dx.doi.org/10.1364/OE.26.033732]

4. Aieta, F.; Kats, M.A.; Genevet, P.; Capasso, F. Multiwavelength achromatic metasurfaces by dispersive phase compensation. Science; 2015; 347, pp. 1342-1345. [DOI: https://dx.doi.org/10.1126/science.aaa2494]

5. Lee, G.-Y.; Yoon, G.; Lee, S.-Y.; Yun, H.; Cho, J.; Lee, K.; Kim, H.; Rho, J.; Lee, B. Complete amplitude and phase control of light using broadband holographic metasurfaces. Nanoscale; 2018; 10, pp. 4237-4245. [DOI: https://dx.doi.org/10.1039/C7NR07154J]

6. Zhao, Y.; Alù, A. Manipulating light polarization with ultrathin plasmonic metasurfaces. Phys. Rev. B; 2011; 84, 205428. [DOI: https://dx.doi.org/10.1103/PhysRevB.84.205428]

7. Khorasaninejad, M.; Chen, W.T.; Devlin, R.C.; Oh, J.; Zhu, A.Y.; Capasso, F. Metalenses at visible wavelengths: Diffraction-limited focusing and subwavelength resolution imaging. Science; 2016; 352, pp. 1190-1194. [DOI: https://dx.doi.org/10.1126/science.aaf6644] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27257251]

8. Wang, J.; Coillet, A.; Demichel, O.; Wang, Z.; Rego, D.; Bouhelier, A.; Grelu, P.; Cluzel, B. Saturable plasmonic metasurfaces for laser mode locking. Light Sci. Appl.; 2020; 9, 50. [DOI: https://dx.doi.org/10.1038/s41377-020-0291-2]

9. Zhang, L.; Sun, X.; Yu, H.; Deng, N.; Qiu, F.; Wang, J.; Qiu, M. Plasmonic metafibers electro-optic modulators. Light Sci. Appl.; 2023; 12, 198. [DOI: https://dx.doi.org/10.1038/s41377-023-01255-7] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37607948]

10. Zhang, C.; Zhang, L.; Zhang, H.; Fu, B.; Wang, J.; Qiu, M. Pulsed polarized vortex beam enabled by metafiber lasers. PhotoniX; 2024; 5, 36. [DOI: https://dx.doi.org/10.1186/s43074-024-00151-6]

11. Dainese, P.; Marra, L.; Cassara, D.; Portes, A.; Oh, J.; Yang, J.; Palmieri, A.; Rodrigues, J.R.; Dorrah, A.H.; Capasso, F. Shape optimization for high efficiency metasurfaces: Theory and implementation. Light Sci. Appl.; 2024; 13, 300. [DOI: https://dx.doi.org/10.1038/s41377-024-01629-5] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/39468011]

12. Yee, K. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media. IEEE Trans. Antennas Propag.; 1966; 14, pp. 302-307. [DOI: https://dx.doi.org/10.1109/TAP.1966.1138693]

13. García de Abajo, F.J.; Howie, A. Retarded field calculation of electron energy loss in inhomogeneous dielectrics. Phys. Rev. B; 2002; 65, 115418. [DOI: https://dx.doi.org/10.1103/PhysRevB.65.115418]

14. Hohenester, U.; Trügler, A. MNPBEM—A Matlab toolbox for the simulation of plasmonic nanoparticles. Comput. Phys. Commun.; 2012; 183, pp. 370-381. [DOI: https://dx.doi.org/10.1016/j.cpc.2011.09.009]

15. Jin, J.-M. The Finite Element Method in Electromagnetics; 3rd ed. John Wiley & Sons: Hoboken, NJ, USA, 2015.

16. Moharam, M.G.; Grann, E.B.; Pommet, D.A.; Gaylord, T.K. Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings. J. Opt. Soc. Am. A; 1995; 12, pp. 1068-1076. [DOI: https://dx.doi.org/10.1364/JOSAA.12.001068]

17. Lalanne, P.; Morris, G.M. Highly improved convergence of the coupled-wave method for TM polarization. J. Opt. Soc. Am. A; 1996; 13, pp. 779-784. [DOI: https://dx.doi.org/10.1364/JOSAA.13.000779]

18. Lalanne, P. Improved formulation of the coupled-wave method for two-dimensional gratings. J. Opt. Soc. Am. A; 1997; 14, pp. 1592-1598. [DOI: https://dx.doi.org/10.1364/JOSAA.14.001592]

19. Li, L. Use of Fourier series in the analysis of discontinuous periodic structures. J. Opt. Soc. Am. A; 1996; 13, pp. 1870-1876. [DOI: https://dx.doi.org/10.1364/JOSAA.13.001870]

20. Li, L. New formulation of the Fourier modal method for crossed surface-relief gratings. J. Opt. Soc. Am. A; 1997; 14, pp. 2758-2767. [DOI: https://dx.doi.org/10.1364/JOSAA.14.002758]

21. Zhao, J.; Tian, X.; Wang, J. Conical diffractions of multilayered gratings modeled by Cartesian rigorous coupled-wave analysis. J. Opt. Soc. Am. A; 2023; 40, pp. 1940-1946. [DOI: https://dx.doi.org/10.1364/JOSAA.499890]

22. Liu, V.; Fan, S. S4: A free electromagnetic solver for layered periodic structures. Comput. Phys. Commun.; 2012; 183, pp. 2233-2244. [DOI: https://dx.doi.org/10.1016/j.cpc.2012.04.026]

23. Sapra, N.V.; Vercruysse, D.; Su, L.; Yang, K.Y.; Skarda, J.; Piggott, A.Y.; Vučković, J. Inverse design and demonstration of broadband grating couplers. IEEE J. Sel. Top. Quantum Electron.; 2019; 25, 6100207. [DOI: https://dx.doi.org/10.1109/JSTQE.2019.2891402]

24. Tamir, T.; Wang, H.C.; Oliner, A.A. Wave propagation in sinusoidally stratified dielectric media. IEEE Trans. Microw. Theory Tech.; 1964; 12, pp. 323-335. [DOI: https://dx.doi.org/10.1109/TMTT.1964.1125815]

25. Moharam, M.G.; Gaylord, T.K. Rigorous coupled-wave analysis of planar-grating diffraction. J. Opt. Soc. Am.; 1981; 71, pp. 811-818. [DOI: https://dx.doi.org/10.1364/JOSA.71.000811]

26. Moharam, M.G.; Gaylord, T.K. Three-dimensional vector coupled-wave analysis of planar-grating diffraction. J. Opt. Soc. Am.; 1983; 73, pp. 1105-1112. [DOI: https://dx.doi.org/10.1364/JOSA.73.001105]

27. Moharam, M.G.; Gaylord, T.K. Rigorous coupled-wave analysis of metallic surface-relief gratings. J. Opt. Soc. Am. A; 1986; 3, pp. 1780-1787. [DOI: https://dx.doi.org/10.1364/JOSAA.3.001780]

28. Li, L. Fourier modal method for crossed anisotropic gratings with arbitrary permittivity and permeability tensors. J. Opt. A Pure Appl. Opt.; 2003; 5, 345. [DOI: https://dx.doi.org/10.1088/1464-4258/5/4/307]

29. Liu, H.; Lalanne, P. Microscopic theory of the extraordinary optical transmission. Nature; 2008; 452, pp. 728-731. [DOI: https://dx.doi.org/10.1038/nature06762]

30. Logofatu, P.C. Rigorous coupled-wave analysis for two-dimensional gratings. Proc. SPIE; 2005; 5972, 59720Q. [DOI: https://dx.doi.org/10.1117/12.639735]

31. Li, J.; Shi, L.; Ma, Y.; Ran, Y.; Liu, Y.; Wang, J. Efficient and stable implementation of RCWA for ultrathin multilayer gratings: T-matrix approach without solving eigenvalues. IEEE Antennas Wirel. Propag. Lett.; 2020; 20, pp. 83-87. [DOI: https://dx.doi.org/10.1109/LAWP.2020.3041299]

32. Xu, J.; Charlton, M.D.B. Highly efficient Rigorous Coupled-Wave Analysis implementation without eigensystem calculation. IEEE Access; 2024; 12, pp. 127966-127975. [DOI: https://dx.doi.org/10.1109/ACCESS.2024.3456798]

33. Moharam, M.G.; Pommet, D.A.; Grann, E.B.; Gaylord, T.K. Stable implementation of the rigorous coupled-wave analysis for surface-relief gratings: Enhanced transmittance matrix approach. J. Opt. Soc. Am. A; 1995; 12, pp. 1077-1086. [DOI: https://dx.doi.org/10.1364/JOSAA.12.001077]

34. Li, L. Multilayer modal method for diffraction gratings of arbitrary profile, depth, and permittivity. J. Opt. Soc. Am. A; 1993; 10, pp. 2581-2591. [DOI: https://dx.doi.org/10.1364/JOSAA.10.002581]

35. Li, L. Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings. J. Opt. Soc. Am. A; 1996; 13, pp. 1024-1035. [DOI: https://dx.doi.org/10.1364/JOSAA.13.001024]

36. Tan, E.L. Enhanced R-matrix algorithms for multilayered diffraction gratings. Appl. Opt.; 2006; 45, pp. 4803-4809. [DOI: https://dx.doi.org/10.1364/AO.45.004803]

37. Rumpf, R.C. Improved formulation of scattering matrices for semi-analytical methods that is consistent with convention. Prog. Electromagn. Res. B; 2011; 35, pp. 241-261. [DOI: https://dx.doi.org/10.2528/PIERB11083107]

38. Tan, E.L. Hybrid-matrix algorithm for rigorous coupled-wave analysis of multilayered diffraction gratings. J. Mod. Opt.; 2006; 53, pp. 417-428. [DOI: https://dx.doi.org/10.1080/09500340500407701]

39. Chen, W.T.; Zhu, A.Y.; Sanjeev, V.; Khorasaninejad, M.; Shi, Z.; Lee, E.; Capasso, F. A broadband achromatic metalens for focusing and imaging in the visible. Nat. Nanotechnol.; 2018; 13, pp. 220-226. [DOI: https://dx.doi.org/10.1038/s41565-017-0034-6] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29292382]