Global rigorous coupled wave analysis for design of multilayer metasurface absorbers

Lin Wang; Lin Wang; Debao Fang; Haibo Jin; Haibo Jin; Jingbo Li; Jingbo Li

doi:10.1364/OE.505139

1. Introduction

Metasurface absorbers have attracted much attention due to their abnormal electromagnetic properties. Their effective permittivity and effective permeability can be easily adjusted by the structure’s parameters. The periodic patterned radar absorbing structures can be modified to significantly improve the absorbing properties and broaden the bandwidth. Up to now, many efforts have been made to achieve a broadband periodic stepped absorbing structure in radar [1–6]. However, all the simulations of those studies at microwave frequencies depend on time-consuming numerical techniques in real space, e.g., the finite difference time domain (FDTD), and the finite element method (FEM).

Rigorous coupled wave analysis (RCWA, or called Fourier modal method, FMM) is a powerful, rapid semi-analytical method to numerically investigate the optical properties of periodic gratings in the infrared region and visible optical regions. After it was developed by Moharam and Gaylord [7] for the analysis of planar diffraction by one-dimensional (1D) binary gratings in 1981, the theory was rapidly extended to the processing of multilayer gratings [8], conical incidence [9], two-dimensional (2D) gratings [10], aperiodic structures [11], and anisotropic gratings [12]. Several approaches like the R matrix [13], the S matrix [14,15], and the enhanced transmittance matrix (ETM) [16] were introduced to settle the numerical stability issue in matching the boundary conditions between layers of this theory. The subsequent developments such as Li’s factorization rules [17] and normal vector fields (NV) [18–20] particularly focused on the convergence problem owing to the wrong treatment of convolution computation of products in truncated Fourier space. In these methods, NV can easily apply the correct Fourier factorization rules over the entire period for arbitrary shape gratings [20] and greatly enhances the convergence of the results.

However, those RCWA research is rarely employed in the study of microwave metasurface absorbers, and the effectiveness and convergence of this method for absorbers composed of dielectric and magnetic materials have not been fully established. Furthermore, considering the vast parameter space involved in practical structural design, it is necessary to propose a method to accelerate the complete search speed of this approach in order to ensure obtaining globally optimal results.

This study proposes an efficient and fast global NV-ETM RCWA method for designing periodic patterned radar absorbing structures in the microwave range. Firstly, a brief introduction to RCWA using NV fields and the ETM approach is presented. Specifically, a new Global NV-ETM RCWA method is proposed to accelerate parameter optimization of complete search. Finally, the effectiveness and convergence of the NV-ETM RCWA method are investigated for dielectric or magnetic metasurface absorbers. In addition, the advantages of computational efficiency and application of Global NV-ETM RCWA method are verified.

2. RCWA with NV field and ETM

The standard RCWA method can be divided into main components: (1) the construction and solution of modal equations for each layer of the periodic grating, and (2) the determination of unknown coefficients in the analytical solutions of these modal equations by enforcing the boundary conditions between layers, thereby obtaining the final reflection and transmission properties. The NV method is utilized during the construction of modal equations to ensure the correct Fourier factorization rules, while the ETM is employed to address the potential numerical instability issue when dealing with specific layer thicknesses in matching the boundary conditions. In this section, the standard transmittance matrix (T-matrix) and ETM method based on the RCWA approach with NV field for multilayer gratings are given. The system depicted in Fig. 1 consists of a trilayer 2D grating sandwiched by reflection region I and transmission region II, for microwave absorbers the region I and II are air and metallic substrate, respectively. Each layer is homogeneous in the z-direction, defined with a thickness h${}_{i}$ and corresponding spatial periods $\Lambda _{x}$ and $\Lambda _{y}$ in the x and y directions respectively with general different filling fractions f${}_{i}$ (see Fig. 1). The incident plane wave is defined by the propagation vector $\mathbf {k}_{inc}$ and the amplitudes in the s and p directions (i.e., perpendicular and parallel to the plane of incidence, respectively). The time harmonic dependence is assumed by exp($+j\omega t$).

Fig. 1. Geometry of general three layer RCWA setup.

Download Full Size | PDF

2.1 Modal equation in periodic region

Referring to the derivation of Eqs. (S1)-(S10), by expressing the periodic distributions of permittivity and permeability in Fourier series and expanding the entire set of Maxwell’s equations in the Fourier domain, we can deduce the modal equations without the correct Fourier factorization rules:

(1)$$\frac{\mathrm{d}}{\mathrm{d} \tilde{z}}\left[\begin{array}{l} \mathbf{u}_{x} \\ \mathbf{u}_{y} \end{array}\right]=\mathbf{Q}\left[\begin{array}{l} \mathbf{s}_{x} \\ \mathbf{s}_{y} \end{array}\right] , \frac{\mathrm{d}}{\mathrm{d} \tilde{z}}\left[\begin{array}{l} \mathbf{s}_{x} \\ \mathbf{s}_{y} \end{array}\right]=\mathbf{P}\left[\begin{array}{l} \mathbf{u}_{x} \\ \mathbf{u}_{y} \end{array}\right].$$

(2a)$$\mathbf{Q}=\left[\begin{array}{cc} \tilde{\mathbf{K}}_{x} \left[\kern-0.15em\left[{ \mu_{\mathrm{r}} }\right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_{y} & \left[\kern-0.15em\left[{ \varepsilon_{\mathrm{r}} }\right]\kern-0.15em\right]-\tilde{\mathbf{K}}_{x} \left[\kern-0.15em\left[{ \mu_{\mathrm{r}} }\right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_{x} \\ \tilde{\mathbf{K}}_{y} \left[\kern-0.15em\left[{ \mu_{\mathrm{r}} }\right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_{y}-\left[\kern-0.15em\left[{ \varepsilon_{\mathrm{r}} }\right]\kern-0.15em\right] & -\tilde{\mathbf{K}}_{y} \left[\kern-0.15em\left[{ \mu_{\mathrm{r}} }\right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_{x} \end{array}\right],$$

(2b)$$\mathbf{P}=\left[\begin{array}{cc} \tilde{\mathbf{K}}_{x} \left[\kern-0.15em\left[{ \varepsilon_{\mathrm{r}} }\right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_{y} & \left[\kern-0.15em\left[{ \mu_{\mathrm{r}} }\right]\kern-0.15em\right]-\tilde{\mathbf{K}}_{x} \left[\kern-0.15em\left[{ \varepsilon_{\mathrm{r}} }\right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_{x} \\ \tilde{\mathbf{K}}_{y} \left[\kern-0.15em\left[{ \varepsilon_{\mathrm{r}} }\right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_{y}-\left[\kern-0.15em\left[{ \mu_{\mathrm{r}} }\right]\kern-0.15em\right] & -\tilde{\mathbf{K}}_{y} \left[\kern-0.15em\left[{ \varepsilon_{\mathrm{r}} }\right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_{x} \end{array}\right].$$

where $\tilde {z}=k_{0} z$, in order to perform computations on a computer, it is necessary to truncate the infinite expansion. The truncation order $N$ corresponds to the values of the Fourier components, which are selected as $-N,\ldots, -1, 0, 1,\ldots, N$. $\tilde {\mathbf {K}}_{x}$ and $\tilde {\mathbf {K}}_{y}$ are normalized wave vectors $\mathbf {K}_{x} /k_{0}$ and $\mathbf {K}_{y} /k_{0}$, $\mathbf {K} _{x}={\rm diag}(k_{x}(m,n))$, $\mathbf {K}_{y}= {\rm diag}(k_{y}(m,n))$, $k_{x}(m,n)=k_{0} \sqrt {\varepsilon _{\rm {r,inc}} \mu _{\rm {r,inc}}} \sin \theta \cos \varphi -m\frac {2\pi }{\Lambda _{x} }$, $k_{y} (m,n)=k_{0} \sqrt {\varepsilon _{\rm {r,inc}} \mu _{\rm {r,inc}}} \sin \theta \sin \varphi -n\frac {2\pi }{\Lambda _{{\rm y}} }$, $\left[\kern-0.15em\left[{ \varepsilon _{r} }\right]\kern-0.15em\right]$ and $\left[\kern-0.15em\left[{ \mu _{r} }\right]\kern-0.15em\right]$ correspond to "Block-Toeplitz Toeplitz-Block" (BTTB) matrixes during the convolution calculation, whose elements are $\varepsilon _{m-p,n-q}$ and $\mu _{m-p,n-q}$, respectively. $\mathbf {s}_{x}$, $\mathbf {s}_{y}$, $\mathbf {s}_{\tilde {z}}$, $\mathbf {u}_{x}$, $\mathbf {u}_{y}$ and $\mathbf {u}_{\tilde {z}}$ are the sets of $\left \{S_{x,mn} (z)\right \}$, $\left \{S_{y,mn} (z)\right \}$, $\left \{S_{z,mn} (z)\right \}$, $\left \{U_{x,mn} (z)\right \}$, $\left \{U_{y,mn} (z)\right \}$ and $\left \{U_{z,mn} (z)\right \}$, respectively. This is consistent with Eqs. (7) and (8) in [15].

The product of two discontinuous functions ( $\varepsilon$ and $\mathbf {E}$, $\mu$ and $\mathbf {H}'$) in the Fourier domain must use the correct Fourier factorization rules [17]. A vector function ("normal vector" field) was first suggested by Popov et al [18] for the differential method and adapted to the RCWA by Schuster et al [19]. The normal vector is given by ${\mathbf {N}} =\left (N_{x} ,N_{y} \right )$, $N_{x} =-{\rm sin}\alpha,N_{y} ={\rm cos} \alpha$. $\alpha$ is the angle between the x coordinate axis and the NV vector. For rectangular and circular structures, NV fields are shown in Fig. 2, the actual NV fields in Section 4 correspond to Fig. 2(a).

Fig. 2. NV fields distribution of (a) rectangular and (b) circular structure grating.

Download Full Size | PDF

Substituting Eqs. (S13) and (S14) with correct factorization rules into Maxwell’s equation in Eq. (S9), the Eq. (2a) with NV field changes to:

(3a)$$\mathbf{Q}=\left[\begin{array}{cc} \tilde{\mathbf{K}}_x \left[\kern-0.15em\left[{ \mu_{r} }\right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_y-\left\{\Delta_1 \left[\kern-0.15em\left[{ N_x N_y }\right]\kern-0.15em\right]\right\}& \left\{\left[\kern-0.15em\left[{ \varepsilon_{r} }\right]\kern-0.15em\right]-\Delta_1 \left[\kern-0.15em\left[{ N_y^2 }\right]\kern-0.15em\right]\right\}-\tilde{\mathbf{K}}_x \left[\kern-0.15em\left[{ \mu_{r} }\right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_x \\ \tilde{\mathbf{K}}_y \left[\kern-0.15em\left[{ \mu_{r} }\right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_y-\left\{\left[\kern-0.15em\left[{ \varepsilon_{r} }\right]\kern-0.15em\right]-\Delta_1 \left[\kern-0.15em\left[{ N_x^2 }\right]\kern-0.15em\right]\right\}& \left\{\Delta_1 \left[\kern-0.15em\left[{ N_x N_y }\right]\kern-0.15em\right]\right\}-\tilde{\mathbf{K}}_y \left[\kern-0.15em\left[{ \mu_{r} }\right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_x \end{array}\right] ,$$

where $\Delta _{1}=\left[\kern-0.15em\left[{ \varepsilon _{\mathrm {r}} }\right]\kern-0.15em\right] -\left[\kern-0.15em\left[{ 1 / \varepsilon _{\mathrm {r}} }\right]\kern-0.15em\right] ^{-1}$.This is consistent with Eq. (8) in [19]. Equation (2b) remains unchanged for pure dielectric materials, however, for magnetic materials, we should do the same analysis for the product of $\mu$ and $\mathbf {H}'$, then Eq. (2b) with NV field changes to:

(3b)$$\mathbf{P}=\left[\begin{array}{cc} \tilde{\mathbf{K}}_x \left[\kern-0.15em\left[{ \varepsilon_{r} }\right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_y-\left\{\Delta_2 \left[\kern-0.15em\left[{ N_x N_y }\right]\kern-0.15em\right]\right\}& \left\{\left[\kern-0.15em\left[{ \mu_{r} }\right]\kern-0.15em\right]-\Delta_2 \left[\kern-0.15em\left[{ N_y^2 }\right]\kern-0.15em\right]\right\}-\tilde{\mathbf{K}}_x \left[\kern-0.15em\left[{ \varepsilon_{r} }\right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_x \\ \tilde{\mathbf{K}}_y \left[\kern-0.15em\left[{ \varepsilon_{r} }\right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_y-\left\{\left[\kern-0.15em\left[{ \mu_{r} }\right]\kern-0.15em\right]-\Delta_2 \left[\kern-0.15em\left[{ N_x^2 }\right]\kern-0.15em\right]\right\}& \left\{\Delta_2 \left[\kern-0.15em\left[{ N_x N_y }\right]\kern-0.15em\right]\right\}-\tilde{\mathbf{K}}_y \left[\kern-0.15em\left[{ \varepsilon_{r} }\right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_x \end{array}\right] ,$$

where $\Delta _{2}=\left[\kern-0.15em\left[{ \mu _{\mathrm {r}} }\right]\kern-0.15em\right] -\left[\kern-0.15em\left[{ 1 / \mu _{\mathrm {r}} }\right]\kern-0.15em\right] ^{-1}$.

By combining the $\textbf {P}$ and Q matrices in Eq. (1), the modal equation is written to eliminate the $\mathbf {H}'$:

(4)$$\frac{\mathrm{d}^{2}}{\mathrm{d} \tilde{z}^{2}}\left[\begin{array}{l} \mathbf{s}_{x} \\ \mathbf{s}_{y} \end{array}\right]=\mathbf{\Omega}^{2}\left[\begin{array}{l} \mathbf{s}_{x} \\ \mathbf{s}_{y} \end{array}\right]=\mathbf{W} \boldsymbol{\beta}^{2} \mathbf{W}^{{-}1}\left[\begin{array}{l} \mathbf{s}_{x} \\ \mathbf{s}_{y} \end{array}\right],$$

where $\mathbf{\Omega }^{2}=\mathbf {P Q}$. W and $\boldsymbol{\beta }$ are, respectively, the eigenvectors and positive square root eigenvalues matrices of the matrix PQ. The analytical solutions of the modal equation for E and $\mathbf {H}'$ are:

(5a)$$\left[\begin{array}{l} \mathbf{s}_{x}(z) \\ \mathbf{s}_{y}(z) \end{array}\right]=\mathbf{W} e^{-\boldsymbol{\beta} k_{0}\left(z-z_{1}\right)} \mathbf{c}^+{+}\mathbf{W} e^{\boldsymbol{\beta} k_{0}\left(z-z_{2}\right)} \mathbf{c}^{-},$$

(5b)$$\left[\begin{array}{l} \mathbf{u}_{x}(z) \\ \mathbf{u}_{y}(z) \end{array}\right]={-}\mathbf{V} e^{-\boldsymbol{\beta} k_{0}\left(z-z_{1}\right)} \mathbf{c}^+{+}\mathbf{V} e^{\boldsymbol\beta} k_{0}\left(z-z_{2}\right) \mathbf{c}^{-},$$

where $\mathbf {V}{\rm =}\mathbf {P}^{-1} \mathbf {W}\boldsymbol{\beta }$. $\mathbf {c}^+$ and $\mathbf {c}^{-}$ are the coefficient vectors of the propagating waves in positive and negative directions, respectively. The layer starts at $z ={z}_{1}$ and ends at $z =z_{2}$. The phase factors, exp($-\boldsymbol{\beta }k_{0}z_{1}$) and exp($-\boldsymbol{\beta }k_{0}z_{2}$), are introduced to prevent numerical overflow [21]. Combining Eqs. (5a) and (5b), the overall field solution is as follows:

(6)$$\boldsymbol{\psi}(z)=\left[\begin{array}{c} \mathbf{s}_{x}(z) \\ \mathbf{s}_{y}(z) \\ \mathbf{u}_{x}(z) \\ \mathbf{u}_{y}(z) \end{array}\right]=\mathbf{F}\left[\begin{array}{cc} e^{-\boldsymbol{\beta} k_{0}\left(z-z_{1}\right)} & \mathbf{0} \\ \mathbf{0} & e^{\boldsymbol{\beta} k_{0}\left(z-z_{2}\right)} \end{array}\right] \mathbf{c},$$

where $\mathbf {F}=\left [\begin {array}{cc} \mathbf {W} & \mathbf {W} \\ \mathbf {-V} & \mathbf {V} \end {array}\right ]$, $\mathbf {c} =\left [\begin {array}{c} {\mathbf {c}^+ } \\ {\mathbf {c}^{-} } \end {array}\right ]$. The function $\boldsymbol {\psi }$ represents the overall modal solution of E and $\mathbf {H}'$ in specific layer, which is the sum of all the modes at plane $z$. F represents the square matrix whose column vectors describe the modes existing in the material. $e^{\boldsymbol{\beta } k_{0}z}$ is the diagonal matrix describing how the modes propagate, which includes accumulation of phase and decaying/growing amplitudes. c is the column vector containing the amplitude coefficient of each of the modes.

When a layer is homogeneous, the solution to the eigen-value problem in Eq. (1) is:

(7)$$\mathbf{W}=\left[\begin{array}{cc} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{I} \end{array}\right], \boldsymbol{\beta}=\left[\begin{array}{cc} j \tilde{\mathbf{K}}_{z} & \mathbf{0} \\ \mathbf{0} & j \tilde{\mathbf{K}}_{z} \end{array}\right],$$

where $\tilde {\mathbf {K}}_{z}=\left (\sqrt {\mu _{\mathrm {r}} \varepsilon _{\mathrm {r}} \mathbf {I}-\tilde {\mathbf {K}}_{x}^{2}-\tilde {\mathbf {K}}_{y}^{2}}\right )^{*}$, $\varepsilon _{{\rm r}}$ and $\mu _{{\rm r}}$ are permittivity and permeability of the uniform layer, the superscript * represents complex conjugation.

The eigen–modes for the magnetic fields are simply:

(8)$$\mathbf{V}=\mathbf{Q} \boldsymbol{\beta}^{{-}1}.$$

Both reflection region and transmission region are homogeneous, but they are composed of vacuum and perfect metal respectively for microwave absorbers. For reflection region, W${}_{\rm ref}$ and V${}_{\rm ref}$ can be directly calculated using $\varepsilon _{{\rm r}} {\rm =}1$. For transmission region, the complex permittivity of the metal in microwave frequency can be written in terms of the metal conductivity $\sigma$, and the free space wavelength $\lambda$ in the SI system [22]:

(9)$$\varepsilon_{\text{metal }}=1-j 60 \sigma \lambda.$$

For simplicity, the dielectric constant of the transmission metal region is chosen to be $1-10^{30} j$ to determine W${}_{\rm trn}$ and V${}_ {\rm trn}$.

2.2 Boundary conditions

By matching the tangential electric-field and magnetic-field components at boundaries between different layers, the standard transfer matrix(T-matrix) method for the whole structure solves the unknown components r and t can be obtained:

(10)$$\mathbf{s}+\mathbf{A} \mathbf{r}=\left(\prod_{i=1}^{N} \mathbf{F}_{i}\left[\begin{array}{cc} \mathbf{X}_{i}^{{-}1} & \mathbf{0} \\ \mathbf{0} & \mathbf{X}_{i} \end{array}\right] \mathbf{F}_{i}^{{-}1}\right) \mathbf{B t},$$

where $\mathbf {A}=\left [\begin {array}{c}\mathbf {W}_{\text {ref }} \\ \mathbf {V}_{\text {ref }}\end {array}\right ]$, $\mathbf {B}=\left [\begin {array}{c}\mathbf {W}_{\mathrm {trn}} \\ -\mathbf {V}_{\mathrm {trn}}\end {array}\right ]$, $\mathbf {F}_{i}=\left [\begin {array}{cc}\mathbf {W}_{i} & \mathbf {W}_{i} \\ -\mathbf {V}_{i} & \mathbf {V}_{i}\end {array}\right ]$, $p_{x}$ and $p_{y}$ are electric field polarization vector in x and y direction, $\mathbf {s}={{\left [ \begin {matrix} {{p}_{x}}{{\boldsymbol{\delta }}_{00,pq}} & {{p}_{y}}{{\boldsymbol{\delta }}_{00,pq}} & {j{{k}_{z,\text { inc }}}{{p}_{y}}}{{\boldsymbol{\delta }}_{00,pq}} & {-j{{k}_{z,\text { inc }}}{{p}_{x}}}{{\boldsymbol{\delta }}_{00,pq}} \\ \end {matrix} \right ]}^{T}}$, $\mathbf {X}_{i}=e^{-k_{0} h_{i} \boldsymbol{\beta } _{i}}$. $\boldsymbol{\delta }$ represents Dirac delta function, $k_{z,\mathrm {inc}}=k_{0}\sqrt {\varepsilon _{\mathrm {r,inc}}\mu _{\mathrm {r,inc}}}\cos \theta$.

However, the T-matrix has been proved to be numerically unstable when the structure along the propagation direction is thick. The inversion of ill-conditioned X${}_{i}$ matrices cannot be represented with sufficient accuracy due to the finite precision of computers. To solve this problem, enhanced transmittance matrix (ETM), S matrix, R matrix are studied. Among these methods, the ETM method is more suitable for the analysis of microwave metasurface absorbers. The enhanced transmittance matrix method to prevent the numerical instability associated with the matrix inversion $\textbf {X}_{i}^{-1}$ is as follows [16]:

Consider boundary at the last N layer interface:

(11)$$\left[\begin{array}{l} \mathbf{a}_{N} \\ \mathbf{b}_{N} \end{array}\right]=\mathbf{F}_{N}^{{-}1} \mathbf{B}.$$

The boundary condition at the N-1 layer interface is

(12)$$\left[\begin{array}{l} \mathbf{a}_{N-1} \\ \mathbf{b}_{N-1} \end{array}\right]=\mathbf{F}_{N-1}^{{-}1} \mathbf{F}_{N}\left[\begin{array}{cc} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{X}_{N} \end{array}\right]\left[\begin{array}{c} \mathbf{I} \\ \mathbf{b}_{N} \mathbf{a}_{N}^{{-}1} \mathbf{X}_{N} \end{array}\right].$$

This process repeats for all the layer(s). After working backward through all interfaces, we obtain an equation of the form:

(13)$$\mathbf{s}+\mathbf{A} \mathbf{r}=\mathbf{F}_{1}\left[\begin{array}{cc} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{X}_{1} \end{array}\right]\left[\begin{array}{c} \mathbf{I} \\ \mathbf{b}_{1} \mathbf{a}_{1}^{{-}1} \mathbf{X}_{1} \end{array}\right] \mathbf{t}_{1}.$$

Equation (13) is easily solved for r and t${}_{1}$ without using $\textbf {X}_{i}^{-1}$ for any number of layers.

(14)$$\left[\begin{array}{l} \mathbf{r} \\ \mathbf{t}_{1} \end{array}\right]=\left[\begin{array}{ll} -\mathbf{A} & \mathbf{B}^{\prime} \end{array}\right]^{{-}1} \mathbf{s},$$

where $\mathbf {B}^{\prime }=\mathbf {F}_{1}\left [\begin {array}{cc}\mathbf {I} & \mathbf {0} \\ \mathbf {0} & \mathbf {X}_{1}\end {array}\right ]\left [\begin {array}{c}\mathbf {I} \\ \mathbf {b}_{\mathbf {1}} \mathbf {a}_{1}^{-1} \mathbf {X}_{1}\end {array}\right ]$.

According to t${}_{1}$, t can be solved by stepping forward through the layers.

(15)$$\left\{\begin{array}{c} \mathbf{t}_{2}=\mathbf{a}_{1}^{{-}1} \mathbf{X}_{1} \mathbf{t}_{1} \\ \mathbf{t}_{3}=\mathbf{a}_{2}^{{-}1} \mathbf{X}_{2} \mathbf{t}_{2} \\ \mathbf{t}_{N}=\mathbf{a}_{N-1}^{{-}1} \mathbf{X}_{N-1} \mathbf{t}_{N-1} \\ \mathbf{t}=\mathbf{a}_{N}^{{-}1} \mathbf{X}_{N} \mathbf{t}_{N} \end{array}\right..$$

First, the longitudinal field components can be found:

(16)$$\begin{aligned} &\mathbf{r}_{z}={-}\tilde{\mathbf{K}}_{z , \text{ ref }}^{{-}1}\left(\tilde{\mathbf{K}}_{x} \mathbf{r}_{x}+\tilde{\mathbf{K}}_{y} \mathbf{r}_{y}\right) ,\\ &\mathbf{t}_{z}={-}\tilde{\mathbf{K}}_{z, \text{ ref }}^{{-}1}\left(\tilde{\mathbf{K}}_{x} \mathbf{t}_{x}+\tilde{\mathbf{K}}_{y} \mathbf{t}_{y}\right). \end{aligned}$$

Second, the diffraction efficiencies of different diffraction orders are:

(17)$$\mathbf{R}=\frac{\operatorname{Re}\left[-\tilde{\mathbf{K}}_{z, \text{ref}}\right]}{\operatorname{Re}\left[k_{z, \text{inc}}\right]}|{\mathbf{r}}|^{2}, \mathbf{T}=\frac{\operatorname{Re}\left[\tilde{\mathbf{K}}_{z,\text{trn}} / \mu_{\mathrm{r}, \text{trn}}\right]}{\operatorname{Re}\left[k_{z, \text{inc}} / \mu_{\mathrm{r}, \text{inc}}\right]}|{\mathbf{t}}|^{2},$$

where $|{\mathbf {r}}|^{2}=\left |\mathbf {r}_{x}\right |^{2}+\left |\mathbf {r}_{y}\right |^{2}+\left |\mathbf {r}_{z}\right |^{2}$, $|{\mathbf {t}}|^{2}=\left |\mathbf {t}_{x}\right |^{2}+\left |\mathbf {t}_{y}\right |^{2}+\left |\mathbf {t}_{z}\right |^{2}.$

The zeroth-order reflection $R_{0}$ and zeroth-order transmission $T_{0}$ are the median of vector $\mathbf {R}$ and vector $\mathbf {T}$, respectively. And the total reflection $R$ and total transmission $T$ are the sum of all diffraction orders (the sum of vector $\mathbf {R}$ and vector $\mathbf {T}$, respectively).

Fortunately, Eq. (15) do not need to calculate for microwave absorbers due to the zero transmission for metal substrate.

2.3 Limitations and advantages

Focusing on periodically structured gratings, the above rigorous coupled wave analysis solves the coupled system of second order partial differential Eq. (4) in Fourier domain, so the main limitation of the RCWA method depends on whether the BTTB matrix formed by the Fourier transform of the permittivity in Eq. (S8) can correctly describe the distribution of permittivity in real space. This limitation is specifically divided into the following two aspects. First, RCWA is not suitable for simulating complex-shaped structures because a high truncation order is required for the BTTB matrix to correctly describe the permittivity in real space, so complex periodic structures are more suitable for discrete calculations in real space with FEM or FDTD. Second, although the correct Fourier factorization rule is used, the Gibbs phenomenon of the Fourier transform of discontinuous functions still exists in Eqs. (S13) and (S14), so for gratings with huge dielectric constant contrast, such as metal and air, the existence of the Gibbs phenomenon will seriously hinder RCWA convergence. However this problem was addressed for the cross-shaped [23] and arbitrary-shaped [24] microwave metal metamaterials, and those methods can be easily combined into the following Global NV-ETM algorithm. Other constraints such as uniform grating in the z-direction and periodic structures can be broken by multilayer staircase approximation and PML boundaries [11].

Periodically structured gratings can be essentially regarded as frequency selective surfaces (FSS), FSS will diffract an applied wave into discrete directions (or called gratings lobes) if the frequency is higher than $f_{c}$. Assuming the FSS is operated in air, this cutoff condition (onset of grating lobes) is $f_{c} =\frac {c}{\Lambda }$. When $f>f_{c}$, the occurrence of grating lobes makes the zeroth order reflection and transmission unable to represent the total reflection and transmission of the structure, while RCWA can easily calculate the total reflection and transmission, and distinguishs the contribution of the grating lobe to the total reflection and transmission. We will discuss it later.

3. Global NV-ETM algorithm

The above-mentioned process is aimed at the calculation of reflection and transmission of a specific system with a fixed period $\Lambda$, filling ratio f and layer thickness h at different frequencies. However, in the actual application, RCWA is usually used to search broadband high or low reflection and transmission system. The most effective method is to simulate all the reflection and transmission spectra with all possible $\Lambda$, f and h combination one by one, however such direct complete search may take months to simulate, although NV-ETM is fast compared with FEM and FDTD. To accelerate the global complete search, we propose the Global NV-ETM algorithm. In this algorithm, all repeated calculation can be avoided and the calculation sequence of the formula is modified to suit for computer matrix acceleration. The proposed algorithm is shown below.

Take the three-layer grating shown in Fig. 1 to illustrate the optimization parameters of complete search. We assume that in the grating regions three layers are the same dielectric material, so optimization parameters are the period length ($\Lambda _{x} =\Lambda _{y} =\Lambda$), the filling fraction (f${}_{1}$, f${}_{2}$) and the respective thickness of the three layer (h${}_{1}$,h${}_{2}$,h${}_{3}$).

The period $\Lambda$ is first determined according to the wavelength range of the study. The period can be fine-tuned after calculating the reflection and transmission spectra of different filling ratios and thicknesses in a specific period, because systems with different periods and same filling ratio and thickness show approximate reflection and transmission spectra. The following calculations can be divided into two steps:

In the step 1, the analytical modal solutions for the first step of standard RCWA are independent of thickness and same for the structures have the same filling ratio, which is one of the motivations behind the Global NV-ETM method. The most time-consuming step in the entire RCWA calculation process is the calculation of eigenvalues of Eq. (4) which is only related to the filling ratio and period, so we first calculate the convolution BTTB matrix of the dielectric constant and magnetic permeability with different filling ratios and selected period, and then calculate the corresponding NV field based on the structure shape. Next, we construct the corresponding P and Q matrix with NV field in Eq. (3), and solve the eigenvalues of the P*Q matrix. Finally, we save W, V and k${}_{0}$*$\boldsymbol{\beta }$ with different filling ratios locally, and load them directly in subsequent step without repeated calculation.

In the step 2, we first construct different three-layer filling ratio combinations according to the actual situation that the filling fraction of the bottom layer is bigger than that of the adjacent layer above. For each individual f combination, load corresponding W${}_{i}$, V${}_{i}$ and k${}_{0}$*$\boldsymbol{\beta }$${}_{i}$. Use the following equation to construct F.

(18)$$\mathbf{F}_{i}=\left[\begin{array}{cc} \mathbf{W}_{i} & \mathbf{W}_{i} \\ -\mathbf{V}_{i} & \mathbf{V}_{i} \end{array}\right], i=1,2,3.$$

Inversion is much slower than matrix multiplication, so the number of inverse operations should be reduced as much as possible. $\mathbf {F}_{23}$ and $\mathbf {F}_{12}$ is defined as

(19)$$\mathbf{F}_{23}=\mathbf{F}_{2}^{{-}1} \mathbf{F}_{3},\mathbf{F}_{12}=\mathbf{F}_{1}^{{-}1} \mathbf{F}_{2}.$$

$\mathbf {F}_{3}$, $\mathbf {F}_{23}$, $\mathbf {F}_{12}$, $\mathbf {F}_{1}$ and k${}_{0}$*$\boldsymbol{\beta }$${}_{3}$, k${}_{0}$*$\boldsymbol{\beta }$${}_{2}$, k${}_{0}$*$\boldsymbol{\beta }$${}_{1}$ are used in the next calculation.

(20)$$\left[\begin{array}{l} \mathbf{a}_{3} \\ \mathbf{b}_{3} \end{array}\right]=\mathbf{F}_{3}^{{-}1} \mathbf{B}.$$

The above calculations are conducted only one cycle. Meanwhile, $\mathbf {d}_{3}$ is defined as $\mathbf {d}_{3}=\mathbf {b}_{3} \mathbf {a}_{3}^{-1}$ to reduce number of inverse operations in the next calculations regarding thickness.

Assuming that each thickness h${}_{i}$ has m different values. For each h${}_{1}$ in h${}_{1}$ in Fig. 1, $\mathbf {X}_{3}=e^{-k_{0} \boldsymbol{\beta }_{3} h_{1}}$, calculate a${}_{2}$, b${}_{2}$:

(21)$$\left[\begin{array}{c} \mathbf{a}_{2} \\ \mathbf{b}_{2} \end{array}\right]=\mathbf{F}_{23}\left[\begin{array}{cc} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{X}_{3} \end{array}\right]\left[\begin{array}{c} \mathbf{I} \\ \mathbf{d}_{3} \mathbf{X}_{3} \end{array}\right].$$

$\mathbf {d}_{2}$ is defined as $\mathbf {d}_{2}=\mathbf {b}_{2} \mathbf {a}_{2}^{-1}$. The number of calculations for above part is $m$ times.

For each h${}_{2}$ in h${}_{2}$, calculate $\mathbf {X}_{2}=e^{-k_{0} \boldsymbol{\beta }_{2} h_{2}}$, calculate a${}_{1}$, b${}_{1}$:

(22)$$\left[\begin{array}{l} \mathbf{a}_{1} \\ \mathbf{b}_{1} \end{array}\right]=\mathbf{F}_{12}\left[\begin{array}{cc} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{X}_{2} \end{array}\right]\left[\begin{array}{c} \mathbf{I} \\ \mathbf{d}_{2} \mathbf{X}_{2} \end{array}\right].$$

$\mathbf {d}_{1}$ is defined as $\mathbf {d}_{1}=\mathbf {b}_{1} \mathbf {a}_{1}^{-1}$. The number of calculations for above part is ${m}^{2}$ times.

For each h${}_{3}$ in h${}_{3}$, calculate $\mathbf {X}_{1}=e^{-k_{0} \boldsymbol{\beta }_{1} h_{3}}$,

(23)$$\mathbf{B}^{\prime}=\mathbf{F}_{1}\left[\begin{array}{cc} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{X}_{1} \end{array}\right]\left[\begin{array}{c} \mathbf{I} \\ \mathbf{d}_{1} \mathbf{X}_{1} \end{array}\right],$$

(24)$$\left[\begin{array}{l} \mathbf{r} \\ \mathbf{t}_{1} \end{array}\right]=\left[\begin{array}{ll} -\mathbf{A} & \mathbf{B}^{\prime} \end{array}\right]^{{-}1} \mathbf{s}.$$

The reflection and transmission are calculated based on Eqs. (15)-(17). The number of calculations is m${}^{3}$ times.

Through the above derivation, the global complete search greatly shortens the computing time for solving reflection and transmission. From step 1, it saves at least m${}^{3}$ (h combination) times of time. From step 2, we assume that the calculation time in different parts is same. It saves (m${}^{3}$+ m${}^{3}$+ m${}^{3}$+ m${}^{3}$)/(1+ m+ m${}^{2}$+ m${}^{3}$)$\mathrm {\approx }$ 4 times of time. In fact, the time saved is much greater than this value due to reduction of the number of inverse operations.

4. Numerical validation and application

To verify the effectiveness and consistency of the proposed method in microwave frequency band, one example, a 7-layer dielectric binary grating [25] in Fig. 3(a) is selected to compare the calculation results of O-ETM, NV-ETM method and High Frequency Structure Simulator (HFSS, frequency domain finite element method, assign the periodic lattice pair boundaries and Floquet port). Figure 3(b) shows the transmission and reflection of this system for the case of normally incident plane wave. The results of Fig. 3(b) are calculated by using O-ETM and NV-ETM with truncation order N=10 and HFSS. The comparison shows good agreement among the three techniques. In order to emphasize the advantages of the NV field technique, we start by calculating the reflection and transmission at 9 GHz displayed in Fig. 3(c). The results are the same between the O-ETM and O-S matrix algorithms, as well as NV-ETM and NV-S matrix, but the S matrix algorithm [15] requires longer calculation time than the ETM under the same truncation order $N$, as shown in Fig. 3(d). The NV-ETM method has a considerably faster convergence than the O-ETM. The curve corresponding to NV-ETM quickly converges to a stable result, and the points corresponding to N = 3, 4, 5 are basically indistinguishable. In contrast, O-ETM converges until N exceeds 8.

Fig. 3. (a) 7-layer dielectric grating slab consisting of four 2D grating slabs supported by three homogenous dielectric slabs. D${}_{x}$ = D${}_{y}$ =10 mm, L${}_{x}$ = L${}_{y}$ = 7 mm, $\varepsilon {}_{rb}$ = 12, h${}_{1}$ = 2 mm, $\varepsilon {}_{r}$=2.2 and h${}_{2}$ = 4 mm. (b) Reflection and Transmission coefficients of 7-layer dielectric grating slabs for normal incidence. Truncation order N=10. (c) Convergence curves of transmission and reflection at 9 GHz. (d) The calculation time of O-S matrix (O for old Eq. (2) without correct Fourier factorization rules), NV-S matrix (NV for Eq. (3)), O-ETM, NV-ETM varies with the truncation order $N$ (Processor: 12-core CPU 5900X with 4.0 GHz and 32 G RAM, Programming language MATLAB).

Download Full Size | PDF

Magnetic materials are seldom utilized in infrared and visible applications. Therefore, a 2-layer magnetic structure composed of $\alpha$ Fe reinforced epoxy resin composite is chosen to validate the effectiveness of NV-ETM method for magnetic metamaterials in the microwave frequency band. The dielectric constant and permeability of the composite and periodic patterned radar absorbing structure can be found in Ref. [5]. Figure 4(a) shows the reflection of the 2-layer periodic magnetic grating slab computed by O-ETM, NV-ETM with N=10 and HFSS. The results yielded by the O-ETM, NV-ETM and HFSS simulations are in good agreement. Notably, the NV-ETM method exhibits significantly faster convergence than the O-ETM. The curve corresponding to NV-ETM converges quickly to a stable result, i.e., converges as N $\mathrm {\ge }$ 3. In contrast, O-ETM converges until $N$ exceeds 10.

Fig. 4. (a) Reflection of 2-layer periodic magnetic grating slab, O-ETM (O for old Eq. (2)), NV-ETM (Eq. (3)), $N$=10; (b) Convergence curves of reflection at 10 GHz.

Download Full Size | PDF

Through the above results of dielectric and magnetic models, it can be concluded that the NV field can effectively accelerate the convergence of the calculation for both dielectric and magnetic metamaterials, and reduce the discrepancy with the HFSS calculation. Comparing the calculation time of reflection of 2-layer periodic magnetic grating slab with NV-ETM, O-ETM and HFSS in Table 1, it is apparent that NV-ETM has a significant efficiency advantage over HFSS. Moreover, the truncation order N = 3 is large enough to assure a good convergence for rectangular shape radar absorbing structure using NV-ETM, so N=3-5 is selected for subsequent calculations.

Table 1. Calculation time of 2-layer periodic magnetic grating slab

View Table | View all tables in this article

We present an additional example to illustrate the advantages of the Global NV-ETM algorithm in terms of computation time. Specifically, 3- and 4-layer periodic stepped structures are designed based on a composite material prepared in our group consisting of reduced graphene (rGO) and thermoplastic polyurethane (TPU). Figure 5 shows the permittivity of the rGO/TPU composite material.

Fig. 5. The permittivity of rGO/TPU composite material.

Download Full Size | PDF

The goal of the optimization is to find an optimal structure which can realize the reflection loss less than -10 dB in the whole range of 2-18 GHz. For the 3-layer square stepped absorbers in Fig. 6(a), the system parameters of complete search are summarized as below:

Fig. 6. (a) 3-layer square stepped absorber, (b) reflection of 3-layer rGO/TPU square stepped absorbers for normal incidence calculated by NV-ETM($N$=3) and HFSS, the number is consistent with Table 3.

Download Full Size | PDF

Frequency: 2-18 GHz, 201 point with corresponding permittivity $\varepsilon {}_{r}$.

Period: $\Lambda _{x} =\Lambda _{y} =\Lambda {\rm =\ }$ 10 mm.

Filling fraction: f${}_{1}$-f${}_{2}$= [0.1,0.2,…,0.9].

Layer thickness: h${}_{1}$-h${}_{3}$= [0,1,2,…,10] mm.

There are 36 combinations for f and 1331 combinations for $h$, so we calculate a total of 47916 kinds of 3-layer stepped structure absorbers. The computation time with truncation order is shown in Table 2. The results indicate that the Global NV-ETM algorithm improves computational efficiency by a factor of 38.604 times for $N$ = 3 compared to the direct NV-ETM algorithm. For $N$ = 4, the improvement is 14.886 times, and for $N$ = 5, the improvement is 10.083 times.

Table 2. 3-layer calculation time with truncation order N

View Table | View all tables in this article

According to the calculation result of Global NV-ETM, there are 654 kinds of 3-layer square stepped absorbers which can achieve effective absorbing from 2 GHz to 18 GHz. We select several relatively thinner absorbers whose parameters are listed in Table 3. The thinnest thickness of 3-layer square stepped absorbers which can achieve effective absorption from 2 GHz to 18 GHz is 15mm, as indicated by the overall thickness h in Table 3. Figure 6(b) displays the reflection of the absorbers in the Table 3 calculated by NV-ETM and HFSS for normal incidence. All absorbers satisfying the optimization goal, and results of NV-ETM are in good consistence with results of HFSS.

Table 3. 3-layer square stepped absorber of 10 mm period

View Table | View all tables in this article

Additionally, the parameters of the 654 structural types are not close. Therefore, for the microwave metasurface absorbers, multiple local optimal solutions exist. While using gradient-based or stochastic methods may yield a convergent local optimal solution, it does not guarantee that the convergent solution is the global optimal one. Moreover, based on the overall results of the three-layer structures, the overall solution in the parameter space is continuous, without any discontinuities. With this information, there will be a new strategy to obtain the global optimal solution more quickly. First, the complete information about the entire parameter space will be acquired through a complete search in a large parameter interval. Then, the global optimal result will be obtained through complete searches in a smaller parameter interval near the solutions that meet or approach the objectives in the large interval solution space.

For the 4-layer square stepped absorbers in Fig. 7(a), the parameters of complete search are shown below:

Fig. 7. (a) 4-layer square stepped absorber, (b) reflection of 4-layer rGO/TPU square stepped absorbers for normal incidence calculated by NV-ETM($N$=3) and HFSS, the number is consistent with Table 5

Download Full Size | PDF

Frequency: 2-18 GHz, 201point with corresponding permittivity $\varepsilon {}_{r}$. period: $\Lambda _{x} =\Lambda _{y} =\Lambda {\rm =}$ 10 mm.

Filling fraction: f${}_{1}$-f${}_{3}$=[0.1,0.2,…,0.9].

Layer thickness: $h_{1}$-$h_{4}$=[2,3…,7]mm.

There are 84 combinations for f and 1296 combinations for h, so we calculate a total of 108864 kinds of 4-layer square stepped absorbers. The calculation time with truncation order is shown in Table 4. The results show that the Global NV-ETM algorithm improves the calculation efficiency significantly, achieving a speedup of 38.784 times, 15.987 times, and 11.443 times for $N$ equals to 3, 4, and 5, respectively.

Table 4. 4-layer calculation time with truncation order

View Table | View all tables in this article

According to the result of Global NV-ETM, there are 5557 types of 4-layer square stepped absorbers which can achieve effective absorption from 2 GHz to 18 GHz, so multilayer for square stepped absorbers is conducive for achieving broadband absorption. We select some relatively thinner 4-layer absorbers whose parameters are listed in Table 5. From the overall thickness $h$, the thinnest thickness of 4-layer square stepped absorbers which can achieve effective absorption from 2 GHz to 18 GHz is 15mm, indicating that multilayer does not reduce absorber thickness. Figure 7(b) illustrates the reflection of the absorbers in Table 5 calculated by NV-ETM and HFSS for normal incidence. All absorbers satisfy the optimization goal, and results of NV-ETM are in good consistence with results of HFSS.Finally, the influence of periods on 3-layer square stepped absorbers is studied. The parameters of the 3-layer square-stepped absorber of different periods that can achieve effective absorbing from 2 GHz to 18 GHz are listed in Table 6. Figure 8 shows the total and zeroth-order reflection of the absorbers in Table 6 calculated by NV-ETM and HFSS for normal incidence. Obviously, the zeroth-order reflection of NV-ETM and HFSS results show consistency, and the total and zeroth-order reflection calculated by NV-ETM for the absorbers with 10 mm period are also coincident, while the total and zeroth-order reflection for the absorbers with 20 mm and 30 mm periods differ greatly at frequencies above 15 GHz and 10 GHz, respectively. For 3-layer square-stepped absorbers of 10 mm, 20 mm and 30 mm periods, the onset of grating lobes ${f}_{c}$ is 30 GHz, 15 GHz and 10 GHz, respectively. When grating lobes occur, the zeroth-order reflection and transmission are not equal to the total reflection and transmission, so the difference between zeroth-order and total reflection in Fig. 8(d) is the contribution of grating lobes to the total reflection. From the total reflection of NV-ETM in Fig. 8(a), the total reflections near the ${f}_{c}$ increase suddenly and significantly when grating lobes occur, indicating that the grating lobe is unfavorable for the absorber.

Fig. 8. Total reflection calculated by (a) NV-ETM ($N$=3), zeroth-order reflection calculated by (b) NV-ETM ($N$=3) and (c) HFSS for different periods 3-layer rGO/TPU square stepped absorbers at normal incidence; (d) the difference between total reflection (a) and zeroth-order reflection (b) calculated by NV-ETM ($N$=3), the number is consistent with Table 6.

Download Full Size | PDF

Table 5. 4-layer square stepped absorber of 10 mm period

View Table | View all tables in this article

Table 6. 3-layer rGO/TPU square stepped absorbers of different period

View Table | View all tables in this article

5. Conclusions

The rigorous coupled-wave analysis method based on normal field and enhanced transmittance matrix approach is adapted to the simulation of the reflection of periodic stepped absorbing structure in radar. The modification formula of NV field on RCWA for magnetic material is proposed and verifies that it highly improves the convergence. The NV-ETM RCWA calculations show a substantial efficiency advantage over HFSS. Furthermore, the Global NV-ETM RCWA algorithms are proposed to improve the calculation efficiency of complete search. Using this method, the optimal 3- and 4-layer square stepped absorbers of 15 mm thickness are obtained, which achieve the 2-18 GHz entire band absorption (RL<10 dB). The optimization results confirm the Global NV-ETM RCWA method significantly increases the computational efficiency, achieving 38 times the efficiency of the direct NV-ETM RCWA calculations in complete search for $N$=3. The proposed method is widely applicable to designing metasurfaces at different frequencies, making it a powerful tool for investigating metasurfaces with various desired functionalities.

Funding

National Natural Science Foundation of China (51772029, 51972029).

Acknowledgments

This work was supported by the National Natural Science Foundation of China under Grants 51972029, 51772029.

Disclosures

The authors declare no conflicts of interest.

Data availability

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.

Supplemental document

See Supplement 1 for supporting content.

References

1. Q. Zhou, T. Shi, B. Xue, et al., “Multi-scale integrated design and fabrication of ultra-broadband electromagnetic absorption utilizing multi-walled carbon nanotubes-based hierarchical metamaterial,” Compos. Sci. Technol. 232, 109877 (2023). [CrossRef]

2. J. Wang, Z. Wu, Y. Xing, et al., “Multi-scale design of ultra-broadband microwave metamaterial absorber based on hollow carbon/mxene/mo2c microtube,” Small 19(14), 2207051 (2023). [CrossRef]

3. Y. Li, Y. Duan, and X. Kang, “Multi-scale integrated design and fabrication of ultrathin broadband microwave absorption utilizing carbon fiber/prussian blue/fe3o4-based lossy lattice metamaterial,” J. Mater. Chem. C 9(19), 6316–6323 (2021). [CrossRef]

4. K. Zhang, J. Zhang, Z. Hou, et al., “Multifunctional broadband microwave absorption of flexible graphene composites,” Carbon 141, 608–617 (2019). [CrossRef]

5. Q. Zhou, X. Yin, F. Ye, et al., “A novel two-layer periodic stepped structure for effective broadband radar electromagnetic absorption,” Mater. Des. 123, 46–53 (2017). [CrossRef]

6. W. Li, T. Wu, W. Wang, et al., “Broadband patterned magnetic microwave absorber,” J. Appl. Phys. 116(4), 044110 (2014). [CrossRef]

7. M. Moharam and T. Gaylord, “Rigorous coupled-wave analysis of planar-grating diffraction,” J. Opt. Soc. Am. 71(7), 811–818 (1981). [CrossRef]

8. L. Li, “Multilayer modal method for diffraction gratings of arbitrary profile, depth, and permittivity,” J. Opt. Soc. Am. A 10(12), 2581–2591 (1993). [CrossRef]

9. M. Moharam, E. B. Grann, D. A. Pommet, et al., “Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings,” J. Opt. Soc. Am. A 12(5), 1068–1076 (1995). [CrossRef]

10. R. Bräuer and O. Bryngdahl, “Electromagnetic diffraction analysis of two-dimensional gratings,” Opt. Commun. 100(1-4), 1–5 (1993). [CrossRef]

11. E. Silberstein, P. Lalanne, J.-P. Hugonin, et al., “Use of grating theories in integrated optics,” J. Opt. Soc. Am. A 18(11), 2865–2875 (2001). [CrossRef]

12. L. Li, “Fourier modal method for crossed anisotropic gratings with arbitrary permittivity and permeability tensors,” J. Opt. A: Pure Appl. Opt. 5(4), 345–355 (2003). [CrossRef]

13. L. Li, “Bremmer series, r-matrix propagation algorithm, and numerical modeling of diffraction gratings,” J. Opt. Soc. Am. A 11(11), 2829–2836 (1994). [CrossRef]

14. L. Li, “Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings,” J. Opt. Soc. Am. A 13(5), 1024–1035 (1996). [CrossRef]

15. R. C. Rumpf, “Improved formulation of scattering matrices for semi-analytical methods that is consistent with convention,” Prog. Electromagn. Res. B 35, 241–261 (2011). [CrossRef]

16. M. Moharam, D. A. Pommet, E. B. Grann, et al., “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 (1995). [CrossRef]

17. L. Li, “Use of fourier series in the analysis of discontinuous periodic structures,” J. Opt. Soc. Am. A 13(9), 1870–1876 (1996). [CrossRef]

18. E. Popov, B. Chernov, M. Nevière, et al., “Differential theory: application to highly conducting gratings,” J. Opt. Soc. Am. A 21(2), 199–206 (2004). [CrossRef]

19. T. Schuster, J. Ruoff, N. Kerwien, et al., “Normal vector method for convergence improvement using the rcwa for crossed gratings,” J. Opt. Soc. Am. A 24(9), 2880–2890 (2007). [CrossRef]

20. P. Götz, T. Schuster, K. Frenner, et al., “Normal vector method for the rcwa with automated vector field generation,” Opt. Express 16(22), 17295–17301 (2008). [CrossRef]

21. G. Granet and B. Guizal, “Efficient implementation of the coupled-wave method for metallic lamellar gratings in tm polarization,” J. Opt. Soc. Am. A 13(5), 1019–1023 (1996). [CrossRef]

22. S. Peng and C. Shiao, “Scattering of plane waves by metallic gratings,” 1994 IEEE MTT-S International Microwave Symposium Digest (Cat. No. 94CH3389-4), (IEEE, 1994), pp. 879–882.

23. L. Wang, D. Fang, H. Jin, et al., “2d rigorous coupled wave analysis with adaptive spatial resolution for a multilayer periodic structure,” Opt. Express 30(12), 21295–21308 (2022). [CrossRef]

24. L. Wang, D. Fang, H. Jin, et al., “Rigorous coupled wave analysis upgraded with combined boundary conditions method for 2d metamaterials in microwave,” Opt. Express 30(17), 29856–29867 (2022). [CrossRef]

25. A. M. Attiya and A. A. Kishk, “Modal analysis of a two-dimensional dielectric grating slab excited by an obliquely incident plane wave,” Prog. Electromagn. Res. 60, 221–243 (2006). [CrossRef]

N	Direct(s)	Step1(s)	Step2(s)
3	125492	16.780	3250
4	228750	40.071	15339
5	504842	85.975	49984

number	$f_{1}$	$f_{2}$	$h_{1}$ - $h_{3}$ (mm)	$h$ (mm)
1	0.8	0.6	6, 6, 3	15
2	0.8	0.7	6, 4, 5	15
3	0.9	0.7	5, 4, 6	15

N	Direct (s)	Step1(s)	Step2(s)
3	311677	16.780	8036
4	615843	40.071	38521
5	1450612	85.975	126669

number	$f_{1}$	$f_{2}$	$f_{3}$	$h_{1}$ - $h_{4}$ (mm)	$h$ (mm)
1	0.9	0.8	0.7	5, 3, 3, 4	15
2	0.9	0.8	0.7	5, 3, 4, 3	15
3	0.8	0.7	0.4	6, 6, 2, 2	16

number	Period(mm)	$f_{1}$	$f_{2}$	$h_{1}$ - $h_{3}$ (mm)	$h$ (mm)
1	10	0.8	0.7	6, 4, 5	15
2	10	0.8	0.7	6, 4.5, 4.5	15
3	20	0.8	0.6	4, 7.5, 4	15.5
4	20	0.8	0.6	4, 8, 3.5	15.5
5	30	0.8	0.6	2.5, 7.5, 6	16
6	30	0.8	0.6	3, 7, 6	16

Global rigorous coupled wave analysis for design of multilayer metasurface absorbers

Abstract

1. Introduction

2. RCWA with NV field and ETM

2.1 Modal equation in periodic region

2.2 Boundary conditions

2.3 Limitations and advantages

3. Global NV-ETM algorithm

4. Numerical validation and application

5. Conclusions

Funding

Acknowledgments

Disclosures

Data availability

Supplemental document

References

Supplementary Material (1)

Data availability

Cited By

Figures (8)

Tables (6)

Equations (27)

Optics Express