Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group
Journal Home About Issues in Progress Current Issue All Issues Feature Issues

2D rigorous coupled wave analysis with adaptive spatial resolution for a multilayer periodic structure

Lin Wang, Debao Fang, Haibo Jin, and Jingbo Li

Open Access Open Access

Abstract

We utilize the 2D modal solution of rigorous coupled wave analysis with adaptive spatial resolution to enhance the reliability of standard RCWA for multilayer gratings with metal-dielectric structures. The conversion relation in modal solutions between the Cartesian system and the transformed system is for the first time established. Based on the proposed conversion relation, the modal solution of the metal-dielectric structure obtained in the transformed system can match the boundary conditions with modal solutions of other different grating layers in the Cartesian system. Numerical results show that even for metal patches in microwave band, the above method can still achieve good convergence and is perfectly suitable for multi-layer structure calculation.

© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

The optical properties of two-dimensional periodic structures, such as photonic crystals and metamaterials, have attracted extensive attention in recent years [1,2]. Several numerical techniques, such as finite difference time domain (FDTD), finite element method (FEM), and Rigorous coupled wave analysis (RCWA), have been established to calculate the optical properties of these structures. Among these techniques, RCWA has significant advantages because the periodic structure of interest can be calculated without spatial or temporal discretization.

The research on the convergence of RCWA mainly focuses on the following two aspects: First, numerical stability in matching boundary conditions of neighboring layers was successfully resolved by S-matrix (S) [3] and the enhanced transmittance matrix (ETM) [4]. Second, the product of discontinuous function in truncated Fourier space was solved by correct factorization rules [5] and normal vectors (NV) field [6]. However, the Gibbs phenomenon of discontinuous functions still exists, which hinders fast convergence, especially for metal-dielectric structures. This problem was addressed by applying the adaptive spatial resolution (ASR) technique, which introduces a new coordinate transformation to increase the spatial resolution around the discontinuities in the permittivity distribution [7]. This mathematical procedure can substantially speed up the convergence and has already proved useful in analyzing gratings at optical frequencies [8,9], metallic lamellar gratings at microwave frequencies [10] and plasmonic resonances [11].

Unfortunately, the ASR method requires solving eigenvalue problem for all layers including homogeneous regions under the exact same ASR transformation function, which enables matching the boundary conditions in the common transformed system. However, the structure shape in different grating layers may vary independently and similar ASR transformation in each layer is not feasible. Therefore, based on the approach which enable the analysis of 1D multilayer structures with ASR [8], we propose a new method for the analysis of 2D multilayer structures with ASR.

In this work, we briefly review the RCWA modal solution by using the ASR and correct Fourier factorization rules in the transformed system. Next, we deduce the conversion relation in modal solutions between the Cartesian system and the transformed system, and propose a modified ETM and S-matrix method to match boundary conditions in the Cartesian system for multilayer periodic structures. Various numerical examples are provided to attest the applicability and efficiency of the proposed method.

2. RCWA in the transformed system

The system depicted in Fig. 1 consists of trilayer 2D lamellar gratings. Each layer is homogeneous in the z-direction, and is defined with a thickness ${h}_{i}$ and the same spatial periods $\Lambda _{x}$ and $\Lambda _{y}$ in the x- and y-directions respectively with different filling fractions f${}_{i}$ (see Fig. 1). The electromagnetic field is defined by the propagation vector k${}_{inc}$ and the amplitudes in the s and p directions (i.e., perpendicular and parallel to the plane of incidence, respectively).

 figure: Fig. 1.

Fig. 1. Geometry of three layer RCWA setup.

Download Full Size | PDF

2.1 Coordinate transformation

Although the correct application of factorization rules such as NV field tremendously improves the convergence of the RCWA in the Cartesian system, the Gibbs phenomenon still hampers convergence, especially for metal-dielectric structures. This problem can however be solved by introducing the concept of adaptive spatial resolution (ASR). In this technique, a new orthogonal coordinate system (u, v, z) is transformed from the Cartesian system (x, y, z), which increases the spatial resolution around discontinuities in the permittivity function. In this work, the transformation equation is chosen in accordance with that is proposed by Vallius [8]:

$$x(u)=a_{1x} +a_{2x} u+\frac{a_{3x} }{2\pi } \sin \left(2\pi \frac{u-u_{l-1} }{u_{l} -u_{l-1} } \right) .$$
where $a_{1x} =\frac {u_{l} x_{l-1} -u_{l-1} x_{l} }{u_{l} -u_{l-1} }$, $a_{2x} =\frac {x_{l} -x_{l-1} }{u_{l} -u_{l-1} }$, $a_{3x} =G\left (u_{l} -u_{l-1} \right )-\left (x_{l} -x_{l-1} \right )$ , and G is an almost zero constant, here chosen to be $G=0.001$. Although jump points x${}_{l}$ in the original Cartesian system are already known, the transformed jump points u${}_{l}$ in the new coordinate system are not predetermined and are chosen to make the third derivative of x(u) continuous, which can further increase the convergence speed of the Fourier series in the transformed space [12]. The choice of jump points is based on the following formula:
$$\Delta u_{k} =\frac{\sqrt[{3}]{\Delta x_{k} } }{\sqrt[{3}]{\sum _{l}\Delta x_{l} } } .$$
where $\Delta u_{k} =u_{k} -u_{k-1}$, $\Delta x_{k} =x_{k} -x_{k-1}$. Figure 2(a) depicts the plot of the ASR transfer function according to the jump point selection in Eq. (2). The first derivative of ASR function is defined as: $f(u)=\frac {\partial x}{\partial u} =\frac {\partial x(u)}{\partial u}$. As shown in Fig. 2(b), the above jump point selection can ensure that the second derivative of f(u) is continuous.

 figure: Fig. 2.

Fig. 2. (a) Graph of x(u) function, (b) the f(u) function and its first and second derivatives functions.

Download Full Size | PDF

The similar relationship between y and v is defined by substituting y and v for x and u respectively.

After the above coordinate transformation in x and y, the distribution diagram of the dielectric constant is shown in Fig. 3. The Fourier transform of the dielectric constant in the uv space is equal to taking points uniformly in the uv space, which is equivalent to increasing spatial resolution around discontinuities in the xy space, so the influence of Gibbs phenomenon can be ignored, which we will discuss in detail later.

 figure: Fig. 3.

Fig. 3. Grid point and dielectric constant distribution in (a) uv and (b) xy space.

Download Full Size | PDF

2.2 Maxwell’s equations in the transformed system

The general form of Maxwell’s equations in the Cartesian system is

$$\left\{\begin{array}{l} \nabla \times \mathbf{H}=j \omega \varepsilon_{0} \varepsilon_{r} \mathbf{E} \\ \nabla \times \mathbf{E}={-}j \omega \mu_{0} \mu_{r} \mathbf{H} \end{array}\right. .$$

By normalizing $\mathbf {H}=j\sqrt {\frac {\varepsilon _{0} }{\mu _{0} } } \mathbf {H}'$, Eq. (3) can be simplified to the following:

$$\left\{\begin{array}{l} \nabla \times \mathbf{H}^{\prime}=k_{0} \varepsilon_{r} \mathbf{E} \\ \nabla \times \mathbf{E}=k_{0} \mu_{r} \mathbf{H}^{\prime} \end{array}\right. ,$$
where $k_{0} {\rm =}\omega \sqrt {\varepsilon _{0} \mu _{0} }$.
$$\left\{\begin{array}{l} \frac{\partial H_{z}^{\prime}}{\partial y}-\frac{\partial H_{y}^{\prime}}{\partial z}=k_{0} \varepsilon_{x y} E_{x} \\ \frac{\partial H_{x}^{\prime}}{\partial z}-\frac{\partial H_{z}^{\prime}}{\partial x}=k_{0} \varepsilon_{x y} E_{y} \\ \frac{\partial H_{y}^{\prime}}{\partial x}-\frac{\partial H_{x}^{\prime}}{\partial y}=k_{0} \varepsilon_{x y} E_{z} \end{array}, \quad \left\{\begin{array}{l} \frac{\partial E_{z}}{\partial y}-\frac{\partial E_{y}}{\partial z}=k_{0} \mu_{x y} H_{x}^{\prime} \\ \frac{\partial E_{x}}{\partial z}-\frac{\partial E_{z}}{\partial x}=k_{0} \mu_{x y} H_{y}^{\prime} \\ \frac{\partial E_{y}}{\partial x}-\frac{\partial E_{x}}{\partial y}=k_{0} \mu_{x y} H_{z}^{\prime} \end{array}\right.\right. .$$

The same vector function expressed in two different coordinate systems is related through the Jacobian matrix as follows:

$$\left[\begin{array}{l} E_{u} \\ E_{v} \end{array}\right]=\left[\begin{array}{ll} \frac{\partial x}{\partial u} & \frac{\partial x}{\partial v} \\ \frac{\partial y}{\partial u} & \frac{\partial y}{\partial v} \end{array}\right]\left[\begin{array}{l} E_{x} \\ E_{y} \end{array}\right]=\left[\begin{array}{cc} f(u) & 0 \\ 0 & g(v) \end{array}\right]\left[\begin{array}{l} E_{x} \\ E_{y} \end{array}\right],$$
where $f(u)=\frac {\partial x}{\partial u} =\frac {\partial x(u)}{\partial u}$, $g(v)=\frac {\partial y}{\partial v} =\frac {\partial y(v)}{\partial v}$.

Thus, we get:

$$E_{x} =\frac{1}{f(u)} E_{u} , E_{y} =\frac{1}{g(v)} E_{v} .$$

In same way, we get:

$$H'_{x} =\frac{1}{f(u)} H'_{u} , H'_{y} =\frac{1}{g(v)} H'_{v} .$$

Substituting Eqs. (7) and (8) into Maxwell’s Eq. (5), we obtain:

$$\left\{\begin{array}{l} \frac{\partial H_{z}^{\prime}}{\partial v}-\frac{\partial H_{v}^{\prime}}{\partial z}=k_{0} \varepsilon^{11} E_{u} \\ \frac{\partial H_{u}^{\prime}}{\partial z}-\frac{\partial H_{z}^{\prime}}{\partial u}=k_{0} \varepsilon^{22} E_{v} \\ \frac{\partial H_{v}^{\prime}}{\partial u}-\frac{\partial H_{u}^{\prime}}{\partial v}=k_{0} \varepsilon^{33} E_{z} \end{array} , \quad\left\{\begin{array}{l} \frac{\partial E_{z}}{\partial v}-\frac{\partial E_{v}}{\partial z}=k_{0} \mu^{11} H_{u}^{\prime} \\ \frac{\partial E_{u}}{\partial z}-\frac{\partial E_{z}}{\partial u}=k_{0} \mu^{22} H_{v}^{\prime} \\ \frac{\partial E_{v}}{\partial u}-\frac{\partial E_{u}}{\partial v}=k_{0} \mu^{33} H_{z}^{\prime} \end{array}\right.\right. ,$$
where $\varepsilon ^{11} =\varepsilon _{uv} \frac {g(v)}{f(u)}$, $\varepsilon ^{22} =\varepsilon _{uv} \frac {f(u)}{g(v)}$, $\varepsilon ^{33} =\varepsilon _{uv} f(u)g(v)$, $\mu ^{11} =\mu _{uv} \frac {g(v)}{f(u)}$, $\mu ^{22} =\mu _{uv} \frac {f(u)}{g(v)}$ and $\mu ^{33} =\mu _{uv} f(u)g(v)$.

We can derive $E_{z}$ and $H'_{z}$ from Eq. (9) and plug them into the remaining Eq. (9), obtaining the system for the tangential components of the fields. Next, we normalize the z coordinate $\tilde {z}=k_{0} z$ and then employ Fourier transforms for the electromagnetic fields ($\mathbf {E}$ and $\mathbf {H}'$) and material parameters ($\varepsilon _{uv}$ and $\mu _{uv}$). Correspondingly, the operator $\frac {\partial }{\partial u}$ and $\frac {\partial }{\partial v}$ are replaced by $j\mathbf {K}_{u}$ and $j\mathbf {K}_{v}$, where ${\mathbf {K}} _{u} {\rm =diag(}k_{u,p} {\rm )}$, $k_{u,p} =k_{0} \sqrt {\varepsilon _{r,} {}_{ref} } \sin \theta \cos \varphi -p\frac {2\pi }{\Lambda _{u} }$; ${\mathbf {K}} _{v} {\rm =diag(}k_{v,q} {\rm )}$, $k_{v,q} =k_{0} \sqrt {\varepsilon _{r,} {}_{ref} } \sin \theta \sin \varphi -q\frac {2\pi }{\Lambda _{v} }$; $\varepsilon ^{11}$, $\varepsilon ^{22}$ , $\varepsilon ^{33}$, $\mu ^{11}$, $\mu ^{22}$ and $\mu ^{33}$ are replaced by corresponding convolution matrices. These replacements allow us to rewrite Maxwell’s equations in Fourier space of this system in the block matrix format:

$$\frac{\mathrm{d}}{\mathrm{d} \tilde{z}}\left[\begin{array}{l} \mathbf{s}_{u} \\ \mathbf{s}_{v} \end{array}\right]=\mathbf{P}_{u v}\left[\begin{array}{l} \mathbf{u}_{u} \\ \mathbf{u}_{v} \end{array}\right], \frac{\mathrm{d}}{\mathrm{d} \tilde{z}}\left[\begin{array}{l} \mathbf{u}_{u} \\ \mathbf{u}_{v} \end{array}\right]=\mathbf{Q}_{u v}\left[\begin{array}{l} \mathbf{s}_{u} \\ \mathbf{s}_{v} \end{array}\right].$$
$$\mathbf{P}_{u v}=\left[\begin{array}{cc} \tilde{\mathbf{K}}_{u} \left[\kern-0.15em\left[ \varepsilon^{33} \right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_{v} & \left[\kern-0.15em\left[ \mu^{22} \right]\kern-0.15em\right]-\tilde{\mathbf{K}}_{u} \left[\kern-0.15em\left[ \varepsilon^{33} \right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_{u} \\ \tilde{\mathbf{K}}_{v} \left[\kern-0.15em\left[ \varepsilon^{33} \right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_{v}-\left[\kern-0.15em\left[ \mu^{11} \right]\kern-0.15em\right] & -\tilde{\mathbf{K}}_{v} \left[\kern-0.15em\left[ \varepsilon^{33} \right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_{u} \end{array}\right],$$
$$\mathbf{Q}_{u v}=\left[\begin{array}{cc} \tilde{\mathbf{K}}_{u} \left[\kern-0.15em\left[ \mu^{33} \right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_{v} & \left[\kern-0.15em\left[ \varepsilon^{22} \right]\kern-0.15em\right]-\tilde{\mathbf{K}}_{u} \left[\kern-0.15em\left[ \mu^{33} \right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_{u} \\ \tilde{\mathbf{K}}_{v} \left[\kern-0.15em\left[ \mu^{33} \right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_{v}-\left[\kern-0.15em\left[ \varepsilon^{11} \right]\kern-0.15em\right] & -\tilde{\mathbf{K}}_{v} \left[\kern-0.15em\left[ \mu^{33} \right]\kern-0.15em\right]^{{-}1} \tilde{\mathbf{K}}_{u} \end{array}\right],$$
where $\tilde {\mathbf {K}}_{u}$ and $\tilde {\mathbf {K}}_{v}$ are normalized wave vectors $\mathbf {K}_{u} / k_{0}$ and $\mathbf {K}_{v} / k_{0}$, $\tilde {z}=k_{0} z$. $\mathbf {s}$ and $\mathbf {u}$ are the 2D Fourier transform coefficients of $\mathbf {E}$ and $\mathbf {H}'$. This equation is consistent with the formula derived from covariant Maxwell’s equations [9].

2.3 2D RCWA with factorization rule in the transformed system

The key to improving the convergence of RCWA in the transformed system is the correct use of Fourier factorization rules in truncated Fourier space. Initially, since $E_{z}$ is continuous in both u and v direction, the displacement ${D}_{z}$ should be factorized by the direct rule in both u and v direction:

$$\left[D_{z}\right]_{12}=\left[\kern-0.15em\left[ \varepsilon^{33} \right]\kern-0.15em\right]_{12}\left[E_{z}\right]_{12}.$$
where the indices behind the brackets indicate that the Fourier expansion has been carried out in u and v direction, respectively. $\left [\kern -0.15em\left [ \varepsilon ^{33} \right ]\kern -0.15em\right ]_{12}$ represents the ’block-Toeplitz Toeplitz-block’ (BTTB) matrix $\left [\kern -0.15em\left [ \varepsilon ^{33} \right ]\kern -0.15em\right ]_{p-p',q-q'}$, because it is a Toeplitz block matrix with its blocks in turn consisting of Toeplitz matrix. $p, q, p'$ and $q'$ are integer. $\left [\kern -0.15em\left [ \varepsilon ^{33} \right ]\kern -0.15em\right ]_{1}$ represents the u-Toeplitz matrix which is made up of elements $\tilde {\varepsilon }_{p-p'}^{33}$ in u direction, $\tilde {\varepsilon }_{p-p'}^{33}$ is the Fourier coefficients of $\varepsilon ^{33}$ in u direction, then the obtained u-Toeplitz matrixes are Fourier-transformed in the v direction, and the v-Fourier-transformed $\left [\kern -0.15em\left [ \varepsilon ^{33} \right ]\kern -0.15em\right ]_{1}$ matrixes are used as elements in the form of a Toeplitz matrix to build up the final BTTB matrix $\left [\kern -0.15em\left [ \varepsilon ^{33} \right ]\kern -0.15em\right ]_{12}$.

From the analysis in Ref. [13], $D_{u} =\varepsilon ^{11} E_{u}$ should be factorized by the inverse rule in the u direction and the direct rule in the v direction, thus:

$$\left[D_{u}\right]_{12}=\left[\kern-0.15em\left[ \left[\kern-0.15em\left[ 1 / \varepsilon^{11}\right]\kern-0.15em\right]_{1}^{{-}1} \right]\kern-0.15em\right]_{2}\left[E_{u}\right]_{12}.$$
where $\left [\kern -0.15em\left [ 1 / \varepsilon ^{11}\right ]\kern -0.15em\right ]_{1}$ represents u-Toeplitz matrix which is made up of Fourier coefficients $\left (1 / \varepsilon ^{11} \right )_{p-p'}$ in u direction, then the inverse of obtained u-Toeplitz matrixes is Fourier-transformed in the v direction, and the v-Fourier-transformed $\left [\kern -0.15em\left [ 1 / \varepsilon ^{11}\right ]\kern -0.15em\right ]_{1}^{-1}$ matrixes are used as elements in the form of a Toeplitz matrix to build up the final BTTB matrix $\left [\kern -0.15em\left [ \left [\kern -0.15em\left [ 1 / \varepsilon ^{11}\right ]\kern -0.15em\right ]_{1}^{-1} \right ]\kern -0.15em\right ]_{2}$.

Next, $D_{v} =\varepsilon ^{22} E_{v}$ must be factorized by the direct rule in the u direction and the inverse rule in the v direction, thus:

$$\left[D_{v}\right]_{12}=\left[\kern-0.15em\left[ \left[\kern-0.15em\left[ \varepsilon^{22}\right]\kern-0.15em\right]_{1}^{{-}1}\right]\kern-0.15em\right]_{2}^{{-}1}\left[E_{v}\right]_{12}.$$
where $\left [\kern -0.15em\left [ \varepsilon ^{22}\right ]\kern -0.15em\right ]_{1}$ represents the u-Toeplitz matrix which is made up of Fourier coefficients $\tilde {\varepsilon }_{p-p'}^{22}$ in u direction, then the inverse of obtained u-Toeplitz matrixes is Fourier-transformed in the v direction, and the v-Fourier-transformed $\left [\kern -0.15em\left [ \varepsilon ^{22}\right ]\kern -0.15em\right ]_{1}^{-1}$ matrixes are used as elements in the form of a Toeplitz matrix to build up the BTTB’ matrix $\left [\kern -0.15em\left [ \left [\kern -0.15em\left [ \varepsilon ^{22}\right ]\kern -0.15em\right ]_{1}^{-1}\right ]\kern -0.15em\right ]_{2}$. The final BTTB matrix is the inverse of the obtained BTTB’ matrix.

A similar relation can be obtained for the magnetic induction B, the permittivity $\mu$, and the magnetic field H. Then the Eq. (11) with factorization rule changes to:

$$\mathbf{P}_{FR}=\left[\begin{array}{cc} \tilde{\mathbf{K}}_{u}\left(\left[\kern-0.15em\left[ \varepsilon^{33} \right]\kern-0.15em\right]_{12}\right)^{{-}1} \tilde{\mathbf{K}}_{v} & \left[\kern-0.15em\left[ \left[\kern-0.15em\left[ \mu^{22} \right]\kern-0.15em\right]_{1}^{{-}1} \right]\kern-0.15em\right]_{2}^{{-}1}-\tilde{\mathbf{K}}_{u}\left(\left[\kern-0.15em\left[ \varepsilon^{33} \right]\kern-0.15em\right]_{12}\right)^{{-}1} \tilde{\mathbf{K}}_{u} \\ \tilde{\mathbf{K}}_{v}\left(\left[\kern-0.15em\left[ \varepsilon^{33} \right]\kern-0.15em\right]_{12}\right)^{{-}1} \tilde{\mathbf{K}}_{v}-\left[\kern-0.15em\left[ \left[\kern-0.15em\left[ 1 / \mu^{11} \right]\kern-0.15em\right]_{1}^{{-}1} \right]\kern-0.15em\right]_{2} & -\tilde{\mathbf{K}}_{v}\left(\left[\kern-0.15em\left[ \varepsilon^{33} \right]\kern-0.15em\right]_{12}\right)^{{-}1} \tilde{\mathbf{K}}_{u} \end{array}\right],$$
$$\mathbf{Q}_{FR}=\left[\begin{array}{cc} \tilde{\mathbf{K}}_{u}\left(\left[\kern-0.15em\left[ \mu^{33} \right]\kern-0.15em\right]_{12}\right)^{{-}1} \tilde{\mathbf{K}}_{v} & \left[\kern-0.15em\left[ \left[\kern-0.15em\left[ \varepsilon^{22} \right]\kern-0.15em\right]_{1}^{{-}1} \right]\kern-0.15em\right]_{2}^{{-}1}-\tilde{\mathbf{K}}_{u}\left(\left[\kern-0.15em\left[ \mu^{33} \right]\kern-0.15em\right]_{12}\right)^{{-}1} \tilde{\mathbf{K}}_{u} \\ \tilde{\mathbf{K}}_{v}\left(\left[\kern-0.15em\left[ \mu^{33} \right]\kern-0.15em\right]_{12}\right)^{{-}1} \tilde{\mathbf{K}}_{v}-\left[\kern-0.15em\left[ \left[\kern-0.15em\left[ 1 / \varepsilon^{11} \right]\kern-0.15em\right]_{1}^{{-}1} \right]\kern-0.15em\right]_{2} & -\tilde{\mathbf{K}}_{v}\left(\left[\kern-0.15em\left[ \mu^{33} \right]\kern-0.15em\right]_{12}\right)^{{-}1} \tilde{\mathbf{K}}_{u} \end{array}\right].$$

2.4 Modal solution in the transformed system

By combining the $\mathbf {P}_{uv}$ and $\mathbf {Q}_{uv}$ matrices, the wave equation is written in terms of S components:

$$\frac{\mathrm{d}^{2}}{\mathrm{~d} \tilde{z}^{2}}\left[\begin{array}{l} \mathbf{s}_{u} \\ \mathbf{s}_{v} \end{array}\right]=\boldsymbol{\Omega}_{u v}^{2}\left[\begin{array}{l} \mathbf{s}_{u} \\ \mathbf{s}_{v} \end{array}\right]=\mathbf{W}_{u v} \boldsymbol{\beta}_{u v}^{2} \mathbf{W}_{u v}^{{-}1}\left[\begin{array}{l} \mathbf{s}_{u} \\ \mathbf{s}_{v} \end{array}\right],$$
where $\boldsymbol {\Omega }_{u v}^{2}=\mathbf {P}_{u v} \mathbf {Q}_{u v}$. $\mathbf {W}_{uv}$ and $\boldsymbol {\beta } _{uv}$ are the eigen-vector matrix and the sqrt eigen-value matrix of the matrix $\mathbf {P}_{u v} \mathbf {Q}_{u v}$, respectively. The analytical solutions of the modal equation for $\mathbf {S}_{uv}$ and $\mathbf {U}_{uv}$ are:
$$\mathbf{S}_{uv}=\left[\begin{array}{c} \mathbf{s}_{u}(z) \\ \mathbf{s}_{v}(z) \end{array}\right]=\mathbf{W}_{u v} e^{-\boldsymbol{\beta}_{u v} k_{0}\left(z-z_{1}\right)} \mathbf{c}_{u v}^{+}+\mathbf{W}_{u v} e^{\boldsymbol{\beta}_{u v} k_{0}\left(z-z_{2}\right)} \mathbf{c}_{u v}^{-},$$
$$\mathbf{U}_{u v}=\left[\begin{array}{c} \mathbf{u}_{u}(z) \\ \mathbf{u}_{v}(z) \end{array}\right]={-}\mathbf{V}_{u v} e^{-\boldsymbol{\beta}_{u v} k_{0}\left(z-z_{1}\right)} \mathbf{c}_{u v}^{+}+\mathbf{V}_{u v} e^{\boldsymbol{\beta}_{u v} k_{0}\left(z-z_{2}\right)} \mathbf{c}_{u v}^{-},$$
where $\mathbf {V}_{u v}=\mathbf {P}_{u v}^{-1} \mathbf {W}_{u v} \boldsymbol {\beta }_{u v}=\mathbf {Q}_{u v} \mathbf {W}_{u v} \boldsymbol {\beta }_{u v}^{-1}$. $\mathbf {c}_{x y}^{+}$ and $\mathbf {c}_{x y}^{-}$ 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, $e^{-\boldsymbol {\beta }_{x y} k_{0}z_{1}}$ and $e^{-\boldsymbol {\beta }_{x y} k_{0}z_{2}}$ are introduced to prevent numerical overflow [14]. Combining Eq. (17) as shown below:
$$\boldsymbol{\Psi}_{uv}(z)=\left[\begin{array}{l} \mathbf{s}_{u}(z) \\ \mathbf{s}_{v}(z) \\ \mathbf{u}_{u}(z) \\ \mathbf{u}_{v}(z) \end{array}\right]=\mathbf{F}_{u v}\left[\begin{array}{cc} e^{-\boldsymbol{\beta}_{uv} k_{0}\left(z-z_{1}\right)} & \mathbf{0} \\ \mathbf{0} & e^{\boldsymbol{\beta}_{uv} k_{0}\left(z-z_{2}\right)} \end{array}\right] \mathbf{c}_{uv},$$
where $\mathbf {F}_{uv}=\left [\begin {array}{cc}\mathbf {W}_{u v} & \mathbf {W}_{u v} \\ -\mathbf {V}_{u v} & \mathbf {V}_{uv}\end {array}\right ]$, $\mathbf {c}_{u v}=\left [\begin {array}{c}\mathbf {c}_{uv}^{+} \\ \mathbf {c}_{uv}^{-}\end {array}\right ]$. The function $\boldsymbol {\Psi }_{uv}(z)$ represents the overall modal solution in the transformed system which is the sum of all the modes at plane $z$. $\mathbf {F}_{uv}$ represents the square matrix whose column vectors describe the modes that exist in the material. $e^{-\boldsymbol {\beta }_{xy} k_{0}z}$ is the diagonal matrix describing how the modes propagate, which includes accumulation of phase and decaying/growing amplitudes. $\mathbf {c}_{uv}$ is the column vector containing the amplitude coefficient of each of the modes.

3. Conversion relation in modal solution between the Cartesian system and the transformed system

In the previous chapters, we have obtained the modal solution Eq. (18) and (S10) in the transformed system and the Cartesian system, respectively. In this chapter, we will derive the conversion relation between these solutions.

In the analysis 2.2, we obtained the conversion relation of the electric field between the Cartesian system and the transformed system in Eq. (7). This means that their Fourier transform coefficients will take the form and the conversion T matrix will be block-diagonal:

$$\left[\begin{array}{l} \mathbf{s}_{x} \\ \mathbf{s}_{y} \end{array}\right]=\mathbf{T}\left[\begin{array}{l} \mathbf{s}_{u} \\ \mathbf{s}_{v} \end{array}\right]=\left[\begin{array}{cc} \mathbf{T}_{x} & \mathbf{0} \\ \mathbf{0} & \mathbf{T}_{y} \end{array}\right]\left[\begin{array}{l} \mathbf{s}_{u} \\ \mathbf{s}_{v} \end{array}\right].$$

Starting from the Fourier decompositions:

$$E_{u}(u, v, z)=\sum_{p} \sum_{q} \mathbf{s}_{u, p q} e^{{-}j\left[k_{x, p} u+k_{y, q} v\right]},$$
$$E_{x}(x, y, z)=\sum_{m} \sum_{n} \mathbf{s}_{x, m n} e^{{-}j\left[k_{x, m} x+k_{y, n} y\right]}.$$

Substituting Eq. (20) into Eq. (7), we easily get:

$$\left[\mathbf{T}_{x}\right]_{p q, m n}=\frac{1}{\Lambda_{x} \Lambda_{y}} \int_{\Lambda_{x}} \int_{\Lambda_{y}} g(v) e^{{-}j\left[k_{x, p} u+k_{y, q} v\right]+j\left[k_{x, m} x(u)+k_{y, n} y(v)\right]} \mathrm{d} u \mathrm{d} v .$$
and it can be efficiently computed as the conjugate result of the Fast Fourier Transform $g(v)e^{-j[k_{x,0} u+k_{y,0} v]{\rm +}j[k_{x,m} x(u)+k_{y,n} y(v)]}$.

In an analogous way, we get:

$$\left[\mathbf{T}_{y}\right]_{p q, m n}=\frac{1}{\Lambda_{x} \Lambda_{y}} \int_{\Lambda_{x}} \int_{\Lambda_{y}} f(u) e^{{-}j\left[k_{x, p} u+k_{y, q} v\right]+j\left[k_{x, m} x(u)+k_{y, n} y(v)\right]} \mathrm{d} u \mathrm{d} v.$$

Substituting the above formula into the modal solution Eq. (17) and (S9), we obtain the conversion relation between E:

$$\mathbf{W}_{x y} e^{-\boldsymbol{\beta}_{x y} k_{0}\left(z-z_{1}\right)} \mathbf{c}_{x y}^{+}+\mathbf{W}_{x y} e^{\boldsymbol{\beta}_{x y} k_{0}\left(z-z_{2}\right)} \mathbf{c}_{x y}^{-}=\mathbf{T}\left(\mathbf{W}_{u v} e^{-\boldsymbol{\beta}_{u v} k_{0}\left(z-z_{1}\right)} \mathbf{c}_{u v}^{+}+\mathbf{W}_{u v} e^{\boldsymbol{\beta}_{u v} k_{0}\left(z-z_{2}\right)} \mathbf{c}_{u v}^{-}\right)$$

In the same way, we obtain the conversion relation between $\mathbf {H}'$:

$$-\mathbf{V}_{x y} e^{-\boldsymbol{\beta}_{x y} k_{0}\left(z-z_{1}\right)} \mathbf{c}_{x y}^{+}+\mathbf{V}_{xy} e^{\boldsymbol{\beta}_{x y} k_{0}\left(z-z_{2}\right)} \mathbf{c}_{x y}^{-}=\mathbf{T}\left(-\mathbf{V}_{u v} e^{-\boldsymbol{\beta}_{u v} k_{0}\left(z-z_{1}\right)} \mathbf{c}_{u v}^{+}+\mathbf{V}_{u v} e^{\boldsymbol{\beta}_{u v} k_{0}\left(z-z_{2}\right)} \mathbf{c}_{u v}^{-}\right)$$

Combining Eqs. (23) and (24) as shown below:

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

The above formula allows us to convert the right-side modal solution obtained in transformed system back to the left-side modal solution in Cartesian system.

4. Application in boundary condition

In order to easily match the boundary conditions, the traditional way is to solve all regions including homogeneous regions and combines boundary conditions in the commmon transformed system, which will lead to a great increase in the amount of calculation and limits the combination of multiple non-uniform layers of different shapes. The study of the above conversion matrix T proposes a simple strategy for combining boundary conditions. This new strategy includes several items: First, the coordinate transformation is only operated for specific layers with huge dielectric constant contrast, such as metamaterial layer for microwave response; and then the modal solution obtained in the transformed system can be transformed back to the original Cartesian system; finally, the solution gets matched to the other layer solution in the Cartesian system. The specific formulas of the commonly used ETM and S algorithms for the above new strategy are given below:

4.1 Standard T matrix and ETM

Assuming that the i-th layer is a metamaterial layer, the modal solution $\mathbf {F}_{uv,i}$ and $\boldsymbol {\beta } _{uv}$ are solved in transformed system. The boundary conditions on the upper and lower sides of the i-th layer are:

$$\mathbf{F}_{x y, i-1}\left[\begin{array}{cc} \mathbf{X}_{x y, i-1} & \mathbf{0} \\ \mathbf{0} & \mathbf{I} \end{array}\right] \mathbf{c}_{x y, i-1}=\left[\begin{array}{cc} \mathbf{T} & \mathbf{0} \\ \mathbf{0} & \mathbf{T} \end{array}\right] \mathbf{F}_{u v, i}\left[\begin{array}{cc} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{X}_{u v, i} \end{array}\right] \mathbf{c}_{uv,i},$$
$$\left[\begin{array}{cc} \mathbf{T} & \mathbf{0} \\ \mathbf{0} & \mathbf{T} \end{array}\right] \mathbf{F}_{u v, i}\left[\begin{array}{cc} \mathbf{X}_{u v, i} & \mathbf{0} \\ \mathbf{0} & \mathbf{I} \end{array}\right] \mathbf{c}_{u v, i}=\mathbf{F}_{x y, i+1}\left[\begin{array}{cc} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{X}_{x y, i+1} \end{array}\right] \mathbf{c}_{x y, i+1},$$
where $\mathbf {X}_{u v, i}=e^{-k_{0} h_{i} \boldsymbol {\beta }_{u v, i}}$, $\mathbf {F}_{u v, i}=\left [\begin{array}{cc}\mathbf {W}_{u v, i} & \mathbf {W}_{u v, i} \\ -\mathbf {V}_{u v, i} & \mathbf {V}_{u v, i}\end {array}\right ]$, $\mathbf {c}_{u v, i}=\left [\begin{array}{l}\mathbf {c}_{u v, i}^{+} \\ \mathbf {c}_{u v, i}^{-}\end {array}\right ]$. In this way, the coordinate transformation is perfectly integrated into the traditional T matrix. Because ETM is an improvement of the standard T, and the formula is similar, it is only necessary to substitute the changes in the T matrix into ETM, so it will not be described in detail here.

4.2 S matrix

Assume that the i-th layer is a metamaterial layer sandwiched between two free-space regions of zero thickness, the Eq. (26) changes to:

$$\left[\begin{array}{cc} \mathbf{W}_{0} & \mathbf{W}_{0} \\ -\mathbf{V}_{0} & \mathbf{V}_{0} \end{array}\right]\left[\begin{array}{l} \mathbf{c}_{x y, i-1}^{+} \\ \mathbf{c}_{x y, i-1}^{-} \end{array}\right]=\left[\begin{array}{cc} \mathbf{T} & \mathbf{0} \\ \mathbf{0} & \mathbf{T} \end{array}\right]\left[\begin{array}{cc} \mathbf{W}_{u v, i} & \mathbf{W}_{u v, i} \\ -\mathbf{V}_{u v, i} & \mathbf{V}_{u v, i} \end{array}\right]\left[\begin{array}{cc} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{X}_{u v, i} \end{array}\right]\left[\begin{array}{l} \mathbf{c}_{u v, i}^{+} \\ \mathbf{c}_{u v, i}^{-} \end{array}\right],$$
$$\left[\begin{array}{cc} \mathbf{T} & \mathbf{0} \\ \mathbf{0} & \mathbf{T} \end{array}\right]\left[\begin{array}{cc} \mathbf{W}_{u v, i} & \mathbf{W}_{u v, i} \\ -\mathbf{V}_{u v, i} & \mathbf{V}_{u v, i} \end{array}\right]\left[\begin{array}{cc} \mathbf{X}_{u v, i} & \mathbf{0} \\ \mathbf{0} & \mathbf{I} \end{array}\right]\left[\begin{array}{c} \mathbf{c}_{u v, i}^{+} \\ \mathbf{c}_{u v, i}^{-} \end{array}\right]=\left[\begin{array}{cc} \mathbf{W}_{0} & \mathbf{W}_{0} \\ -\mathbf{V}_{0} & \mathbf{V}_{0} \end{array}\right]\left[\begin{array}{c} \mathbf{c}_{x y, i+1}^{+} \\ \mathbf{c}_{x y, i+1}^{-} \end{array}\right],$$
where $\mathbf {W}_{0}$ and $\mathbf {V}_{0}$ are eigen-modes for free space in the Cartesian system. The S-matrices of layer i are determined by connecting the coefficients in the left and right media as
$$\left[\begin{array}{l} \mathbf{c}_{x y, i-1}^{-} \\ \mathbf{c}_{x y, i+1}^{+} \end{array}\right]=\left[\begin{array}{ll} \mathbf{S}_{11}^{i} & \mathbf{S}_{12}^{i} \\ \mathbf{S}_{21}^{i} & \mathbf{S}_{22}^{i} \end{array}\right]\left[\begin{array}{l} \mathbf{c}_{x y, i-1}^{+} \\ \mathbf{c}_{x y, i+1}^{-} \end{array}\right].$$

After some algebra, the following expressions can be derived for the S-matrices [15]:

$$\mathbf{S}_{11}^{(i)}=\mathbf{S}_{22}^{(i)}=\left(\mathbf{A}_{i}-\mathbf{X}_{i} \mathbf{B}_{i} \mathbf{A}_{i}^{{-}1} \mathbf{X}_{i} \mathbf{B}_{i}\right)^{{-}1}\left(\mathbf{X}_{i} \mathbf{B}_{i} \mathbf{A}_{i}^{{-}1} \mathbf{X}_{i} \mathbf{A}_{i}-\mathbf{B}_{i}\right),$$
$$\mathbf{S}_{12}^{(i)}=\mathbf{S}_{21}^{(i)}=\left(\mathbf{A}_{i}-\mathbf{X}_{i} \mathbf{B}_{i} \mathbf{A}_{i}^{{-}1} \mathbf{X}_{i} \mathbf{B}_{i}\right)^{{-}1} \mathbf{X}_{i}\left(\mathbf{A}_{i}-\mathbf{B}_{i} \mathbf{A}_{i}^{{-}1} \mathbf{B}_{i}\right),$$
where $\mathbf {A}_{i}=\left (\mathbf {T} \mathbf {W}_{u v, j}\right )^{-1} \mathbf {W}_{0}+\left (\mathbf {T} \mathbf {V}_{u v, j}\right )^{-1} \mathbf {V}_{0}$ and $\mathbf {B}_{i}=\left (\mathbf {T} \mathbf {W}_{u v, j}\right )^{-1} \mathbf {W}_{0}-\left (\mathbf {T} \mathbf {V}_{u v, j}\right )^{-1} \mathbf {V}_{0}$. The combination method with other layers is through the common Redheffer Star Product [15] in the Cartesian system.

5. Numerical validation and application

Since both ETM and S matrix algorithms are correct algorithms to remove the influence of inversion of the ill-conditioned matrix X${}_{i}$, the difference between these results is even lower than 3*10${}^{-15}$ according to our numerical results. Therefore, the following numerical results are calculated based on the ETM algorithm.

As the first numerical example to demonstrate the efficiency and stability of the proposed method, the 7-layer dielectric binary grating in Fig. 4(a) [16] is chosen to compare the computational results of ASR, adaptive spatial resolution with factorization rules (ASR-FR) and High Frequency Structure Simulator (HFSS) simulations. For this dielectric grating, the ASR is only solved within the inhomogeneous layer, and then the boundary condition is matched with the homogeneous layer in the Cartesian system according to the proposed method. Figure 4(b) shows the transmission and reflection of this structure for normal incident, which are calculated by using ASR and ASR-FR with truncation order N=8 and HFSS. The comparison shows good agreement among the three techniques.

 figure: Fig. 4.

Fig. 4. (a) 7-layer dielectric grating slab. ${D}_{x}= {D}_{y} =10\, \rm {mm}$, ${L}_{x} = L_{y} = 7\, \rm {mm}$, $\varepsilon _{rb} = 12$, ${h}_{1} = 2\, \rm {mm}$, $\varepsilon _{r}=2.2$ and ${h}_{2} = 4\, \rm {mm}$. (b) Reflection and Transmission of 7-layer dielectric grating slabs for normal incidence.

Download Full Size | PDF

The reflection and transmission of this structure at 9 GHz is plotted versus the truncation order in Fig. 5(a) to highlight the features of the ASR technique. First, ASR-FR and enhanced transmittance matrix with normal vectors [6] (ETM-NV) with the correct Fourier factorization rule have a considerably faster convergence than ASR and ETM. The curves corresponding to ASR-FR and ETM-NV very quickly converge to stable results, and the points corresponding to N = 3, 4, 5 are basically indistinguishable. Conversely, the results of ASR and ETM converge until N exceeds 8. Next, compared with stable ETM and ETM-NV results obtained in the Cartesian system, the ASR and ASR-FR results show unstable oscillation for truncation order N$\mathrm {>}$22. The incurred instability is further examined in Fig. 5(b), where the condition number of the conversion matrix T increases exponentially with the truncation order N. For both the S algorithm and the ETM algorithm, the inversion of the conversion matrix T is inevitably used in the process of matching the boundary condition. When the truncation order N$\mathrm {>}$22, due to the limited computational precision, the inversion operation of T will bring a certain error, resulting in oscillation. Fortunately, before the conversion matrix reaches instability, the results have already converged, and in practical applications, it takes a long time to compute structures with the truncation order over 22, and the greatest computational advantage of RCWA disappears, so the influence of the inversion of the possibly ill-conditioned T matrix does not need to be considered.

 figure: Fig. 5.

Fig. 5. (a) Convergence curves of transmission and reflection at 9 GHz, (b) The condition number of the conversion matrix T, versus the truncation order N=M, Program language: MATLAB, double precision.

Download Full Size | PDF

Under the application of ETM or S-matrix and correct factorization rule, the accuracy of results in RCWA method only depends on whether the BTTB matrixes formed by the Fourier transform of the permittivity and permeability can correctly represent the actual distribution of permittivity and permeability in real space. For standard RCWA in the Cartesian system, those BTTB matrixes are formed by direct Fourier transform of $\varepsilon _{xy}$ and $\mu _{xy}$ in Eq. (S5) and (S7). For RCWA in the transformed system, those BTTB matrixes are formed by direct Fourier transform of $\varepsilon ^{11}$, $\varepsilon ^{22}$, $\varepsilon ^{33}$, $\mu ^{11}$, $\mu ^{22}$ and $\mu ^{33}$ in Eqs. (11) and (15). It is only necessary to compare the difference between the reconstructed function of the Fourier transform and the actual function in real space to directly demonstrate the accuracy of the BTTB matrix consisting of the Fourier transform coefficient.

As another example, a simple metal square patch in microwave is considered to illustrate the difference of $\varepsilon _{xy}$ in the Cartesian system and $\varepsilon ^{33}$ in the transformed system. The grating parameters are as follows: $\Lambda _{x} =\Lambda _{y} =30\, \rm {mm}$, $f=0.5$, $h=0.01\, \rm {mm}$. The relative 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 [10]: $\varepsilon _{metal} =1-ja\sigma \lambda$, where $a=60\,\rm {S^{-1}}$, there is chosen to be $\varepsilon _{metal} =1-j*10^{6}$. The real parts of permittivity in metal and air are both 1, and the imaginary parts of permittivity in metal and air are 10${}^{6}$ and 0, respectively. So only the reconstruction of the imaginary part needs to be compared with the actual value in real space. Figure 6 shows the imaginary part of the actual $\varepsilon _{xy}$ function and the reconstructed $\varepsilon _{xy}$ function in the Cartesian system. $\varepsilon _{xy}$ is discontinuous in both $x$ and $y$ directions in Fig. 6(a). Significant overshoot around the discontinuity is found in the reconstructed $\varepsilon _{xy}$ function in Fig. 6(b). This overshoot is the Gibbs phenomenon and the increasion of truncation order N does nothing to remove this overshoot but merely moves it closer to the point of discontinuity. The value of overshoot has a fixed value with the jump value between discontinuity, For 1D, this value is about 0.18. For gratings with small dielectric constant contrast, such as first example (1:12), this error is also small, so ETM-NV in the Cartesian system can converge, even in the presence of Gibbs phenomenon. But for the microwave metal square patch with huge dielectric constant contrast ($10^6:0$), the error caused by the Gibbs phenomenon in Fig. 6(c) cannot be ignored, so the ETM-NV method in the Cartesian system cannot converge correctly regardless of the truncation order.

 figure: Fig. 6.

Fig. 6. (a) 3D graphs of the imaginary part of $\varepsilon _{xy}$ and (b) reconstruction of the imaginary part of $\varepsilon _{xy}$ (c) error between (a) actual value and (b) reconstruction in the imaginary part of $\varepsilon _{xy}$, N=M=40.

Download Full Size | PDF

Figure 7 shows the imaginary part of the actual $\varepsilon _{uv}$ function, actual $\varepsilon ^{33}$ function and the reconstructed $\varepsilon ^{33}$ function in the transformed system. $\varepsilon _{uv}$ is discontinuous in both $u$ and $v$ directions in Fig. 7(a), Note, $\varepsilon _{uv}$ is different with $\varepsilon _{xy}$, which is discussed in Fig. 3. From the inset of Fig. 7(b), $\varepsilon ^{33}$ is still a discontinuous function at the discontinuity point, while from the overall view of Fig. 7(b), $\varepsilon ^{33}$ is more like a continuous function. Significant overshoot around the discontinuity is not found in the reconstructed $\varepsilon ^{33}$ function in Fig. 7(c), because Gibbs phenomenon only occurs in discontinuous function and Gibbs phenomenon is greatly suppressed in pseudo-continuous $\varepsilon ^{33}$ function. This is the reason for keeping the second derivative of $f(u)$ continuous in Eq. (2). The error between actual value and reconstruction of $\varepsilon ^{33}$ shown in Fig. 7(d) can be ignored compared with the actual value. This is the essential reason why the ASR method can converge correctly for structures with huge dielectric constant contrast such as metal and air, so other coordinate transformation functions can also be used, as long as the continuity of $\varepsilon ^{33}$ in the transformed system is guaranteed. The same conclusions can be drawn for the analysis of $\varepsilon ^{11}$, $\varepsilon ^{22}$, $\mu ^{11}$, $\mu ^{22}$ and $\mu ^{33}$ in the transformed system.

 figure: Fig. 7.

Fig. 7. (a) 3D graphs of the imaginary part of $\varepsilon _{uv}$, (b) $\varepsilon ^{33}$ (The inset shows a magnified view of $\varepsilon ^{33}$) and (c) reconstruction of the imaginary part of $\varepsilon ^{33}$, (d) error between (b) actual value and (c) reconstruction in the imaginary part of $\varepsilon ^{33}$, N=M=20.

Download Full Size | PDF

The zeroth-order and total reflection and transmission of the metal square patch for normal incidence are displayed in Fig. 8(a). Clearly, the zeroth-order reflection and transmission of ASR-FR and HFSS show consistency, and total reflection and transmission of ASR-FR and ASR also show consistency. The zeroth-order and total reflection and transmission coincide only at frequencies below 10 GHz and differ greatly at frequencies above 10 GHz. We calculate the difference between zeroth-order and total transmission, reflection in Fig. 8(b). 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 for normal incidence, this cutoff condition is $f_{c} =\frac {c}{\Lambda }$, for this example, $f_{c} =10 \, {\rm GHz}$. Therefore, the difference between the zeroth-order and total reflection and transmission is caused by gratings lobes, the zeroth-order reflection and transmission can be considered as total reflection and transmission only at frequencies below $f_{c}$.

 figure: Fig. 8.

Fig. 8. (a) Zeroth-order and total transmission, reflection and (b) difference between zeroth-order and total transmission, reflection of metal square patch at 2-18 GHz for normal incidence, ASR with N=M=20, ASR-FR with N=M=8.

Download Full Size | PDF

We chose the reflection and transmission of this structure at 6 GHz for comparison in Fig. 9 to illustrate the importance of the Fourier factorization rule for convergence. For metamaterials with huge dielectric constant contrast, the application of the Fourier factorization rule has a more prominent effect on the convergence. ASR-FR result shows convergence when N=5, while ASR result slowly converges when N$\mathrm {>}$20.

 figure: Fig. 9.

Fig. 9. Convergence curves for (a) reflection and (b) transmission of metal square patch at 6 GHz under normal incidence.

Download Full Size | PDF

The above computational method for metamaterials are easily integrated into the previously published Global ETM-NV algorithm [17] to optimize broadband absorbing structures. A 3-layer periodic stepped structures with FSS is designed based on a reduced graphene (rGO) and thermoplastic polyurethane (TPU) composite material prepared in our group. The relative permittivity of the rGO/TPU composite material is shown in Fig. 10.

 figure: Fig. 10.

Fig. 10. Relative permittivity of rGO/TPU composite material

Download Full Size | PDF

 figure: Fig. 11.

Fig. 11. (a) 3-layer rGO/TPU grating slab with FSS, transmission region II is a metal substrate. ${L}_{x}=10\, \rm {mm}$, ${L}_{1}=8\, \rm {mm}$, ${L}_{2}=6\, \rm {mm}$, ${L}_{3}=4\, \rm {mm}$, ${h}_{1}=6\, \rm {mm}$, ${h}_{2}=6\, \rm {mm}$, ${h}_{3}=3\, \rm {mm}$. (b) reflection of 3-layer dielectric grating slab with FSS for normal incidence. ASR-FR with N=M=8.

Download Full Size | PDF

Figure 11(b) shows that the ASR-FR result is in good consistence with the HFSS result, so crossed metamaterials can be calculated correctly using RCWA in microwaves and provides guidance for designing absorbing structures.

6. Conclusion

We have utilized the 2D ASR modal solution in the transformed system to enhance the reliability of RCWA for multilayer gratings. The conversion relation in modal solution between the Cartesian system and the transformed system is proposed to match boundary conditions among different grating layers in the Cartesian system. Great convergence rates were achieved and the method proved reliable also for the structures with grating lobes. It must be pointed out that since only ASR coordinate transformation is used, the above conversion relation is only for crossed surface-relief gratings. For gratings of arbitrary shape, additional coordinate transformation such as matched coordinates in Ref. [9] is required.

Funding

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

Acknowledgments

This work was supported by the National Natural Science Foundation of China.

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. Y. Akahane, T. Asano, B. S. Song, and S. Noda, “High-q photonic nanocavity in a two-dimensional photonic crystal,” Nature 425(6961), 944–947 (2003). [CrossRef]  

2. A. Poddubny, I. Iorsh, P. Belov, and Y. Kivshar, “Hyperbolic metamaterials,” Nat. Photonics 7(12), 948–957 (2013). [CrossRef]  

3. L. F. 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]  

4. M. G. Moharam, D. A. Pommet, E. B. Grann, and T. K. Gaylord, “Stable implementation of the rigorous coupled-wave analysis for surface-relief gratings: enhanced transmittance matrix approach,” J. Opt. Soc. Am. A 12(5), 1077–1086 (1995). [CrossRef]  

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

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

7. G. Granet, “Reformulation of the lamellar grating problem through the concept of adaptive spatial resolution,” J. Opt. Soc. Am. A 16(10), 2510–2516 (1999). [CrossRef]  

8. T. Vallius and M. Honkanen, “Reformulation of the fourier modal method with adaptive spatial resolution: application to multilevel profiles,” Opt. Express 10(1), 24–34 (2002). [CrossRef]  

9. T. Weiss, G. Granet, N. A. Gippius, S. G. Tikhodeev, and H. Giessen, “Matched coordinates and adaptive spatial resolution in the fourier modal method,” Opt. Express 17(10), 8051–8061 (2009). [CrossRef]  

10. A. Khavasi and K. Mehrany, “Adaptive spatial resolution in fast, efficient, and stable analysis of metallic lamellar gratings at microwave frequencies,” IEEE Trans. Antennas Propag. 57(4), 1115–1121 (2009). [CrossRef]  

11. T. Weiss, N. A. Gippius, S. G. Tikhodeev, G. Granet, and H. Giessen, “Derivation of plasmonic resonances in the fourier modal method with adaptive spatial resolution and matched coordinates,” J. Opt. Soc. Am. A 28(2), 238–244 (2011). [CrossRef]  

12. A. Khavasi and K. Mehrany, “Regularization of jump points in applying the adaptive spatial resolution technique,” Opt. Commun. 284(13), 3211–3215 (2011). [CrossRef]  

13. L. F. Li, “New formulation of the fourier modal method for crossed surface-relief gratings,” J. Opt. Soc. Am. A 14(10), 2758–2767 (1997). [CrossRef]  

14. 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]  

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. 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]  

17. L. Wang, D. B. Fang, H. B. Jin, and J. B. Li, “A fast modified global rigorous coupled wave analysis for multilayer 2d periodic radar absorber,” IEEE Trans. Antennas Propag. (to be published).

Supplementary Material (1)

NameDescription
Supplement 1       The modal solution of RCWA in the Cartesian system

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.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (11)
Fig. 1.
Fig. 1. Geometry of three layer RCWA setup.

View in Article | Download Full Size | PDF

Fig. 2.
Fig. 2. (a) Graph of x(u) function, (b) the f(u) function and its first and second derivatives functions.

View in Article | Download Full Size | PDF

Fig. 3.
Fig. 3. Grid point and dielectric constant distribution in (a) uv and (b) xy space.

View in Article | Download Full Size | PDF

Fig. 4.
Fig. 4. (a) 7-layer dielectric grating slab. ${D}_{x}= {D}_{y} =10\, \rm {mm}$, ${L}_{x} = L_{y} = 7\, \rm {mm}$, $\varepsilon _{rb} = 12$, ${h}_{1} = 2\, \rm {mm}$, $\varepsilon _{r}=2.2$ and ${h}_{2} = 4\, \rm {mm}$. (b) Reflection and Transmission of 7-layer dielectric grating slabs for normal incidence.

View in Article | Download Full Size | PDF

Fig. 5.
Fig. 5. (a) Convergence curves of transmission and reflection at 9 GHz, (b) The condition number of the conversion matrix T, versus the truncation order N=M, Program language: MATLAB, double precision.

View in Article | Download Full Size | PDF

Fig. 6.
Fig. 6. (a) 3D graphs of the imaginary part of $\varepsilon _{xy}$ and (b) reconstruction of the imaginary part of $\varepsilon _{xy}$ (c) error between (a) actual value and (b) reconstruction in the imaginary part of $\varepsilon _{xy}$, N=M=40.

View in Article | Download Full Size | PDF

Fig. 7.
Fig. 7. (a) 3D graphs of the imaginary part of $\varepsilon _{uv}$, (b) $\varepsilon ^{33}$ (The inset shows a magnified view of $\varepsilon ^{33}$) and (c) reconstruction of the imaginary part of $\varepsilon ^{33}$, (d) error between (b) actual value and (c) reconstruction in the imaginary part of $\varepsilon ^{33}$, N=M=20.

View in Article | Download Full Size | PDF

Fig. 8.
Fig. 8. (a) Zeroth-order and total transmission, reflection and (b) difference between zeroth-order and total transmission, reflection of metal square patch at 2-18 GHz for normal incidence, ASR with N=M=20, ASR-FR with N=M=8.

View in Article | Download Full Size | PDF

Fig. 9.
Fig. 9. Convergence curves for (a) reflection and (b) transmission of metal square patch at 6 GHz under normal incidence.

View in Article | Download Full Size | PDF

Fig. 10.
Fig. 10. Relative permittivity of rGO/TPU composite material

View in Article | Download Full Size | PDF

Fig. 11.
Fig. 11. (a) 3-layer rGO/TPU grating slab with FSS, transmission region II is a metal substrate. ${L}_{x}=10\, \rm {mm}$, ${L}_{1}=8\, \rm {mm}$, ${L}_{2}=6\, \rm {mm}$, ${L}_{3}=4\, \rm {mm}$, ${h}_{1}=6\, \rm {mm}$, ${h}_{2}=6\, \rm {mm}$, ${h}_{3}=3\, \rm {mm}$. (b) reflection of 3-layer dielectric grating slab with FSS for normal incidence. ASR-FR with N=M=8.

View in Article | Download Full Size | PDF

Equations (36)

Equations on this page are rendered with MathJax. Learn more.

x ( u ) = a 1 x + a 2 x u + a 3 x 2 π sin ⁡ ( 2 π u − u l − 1 u l − u l − 1 ) .
Δ u k = Δ x k 3 ∑ l Δ x l 3 .
{ ∇ × H = j ω ε 0 ε r E ∇ × E = − j ω μ 0 μ r H .
{ ∇ × H ′ = k 0 ε r E ∇ × E = k 0 μ r H ′ ,
{ ∂ H z ′ ∂ y − ∂ H y ′ ∂ z = k 0 ε x y E x ∂ H x ′ ∂ z − ∂ H z ′ ∂ x = k 0 ε x y E y ∂ H y ′ ∂ x − ∂ H x ′ ∂ y = k 0 ε x y E z , { ∂ E z ∂ y − ∂ E y ∂ z = k 0 μ x y H x ′ ∂ E x ∂ z − ∂ E z ∂ x = k 0 μ x y H y ′ ∂ E y ∂ x − ∂ E x ∂ y = k 0 μ x y H z ′ .
[ E u E v ] = [ ∂ x ∂ u ∂ x ∂ v ∂ y ∂ u ∂ y ∂ v ] [ E x E y ] = [ f ( u ) 0 0 g ( v ) ] [ E x E y ] ,
E x = 1 f ( u ) E u , E y = 1 g ( v ) E v .
H x ′ = 1 f ( u ) H u ′ , H y ′ = 1 g ( v ) H v ′ .
{ ∂ H z ′ ∂ v − ∂ H v ′ ∂ z = k 0 ε 11 E u ∂ H u ′ ∂ z − ∂ H z ′ ∂ u = k 0 ε 22 E v ∂ H v ′ ∂ u − ∂ H u ′ ∂ v = k 0 ε 33 E z , { ∂ E z ∂ v − ∂ E v ∂ z = k 0 μ 11 H u ′ ∂ E u ∂ z − ∂ E z ∂ u = k 0 μ 22 H v ′ ∂ E v ∂ u − ∂ E u ∂ v = k 0 μ 33 H z ′ ,
d d z ~ [ s u s v ] = P u v [ u u u v ] , d d z ~ [ u u u v ] = Q u v [ s u s v ] .
P u v = [ K ~ u [ [ ε 33 ] ] − 1 K ~ v [ [ μ 22 ] ] − K ~ u [ [ ε 33 ] ] − 1 K ~ u K ~ v [ [ ε 33 ] ] − 1 K ~ v − [ [ μ 11 ] ] − K ~ v [ [ ε 33 ] ] − 1 K ~ u ] ,
Q u v = [ K ~ u [ [ μ 33 ] ] − 1 K ~ v [ [ ε 22 ] ] − K ~ u [ [ μ 33 ] ] − 1 K ~ u K ~ v [ [ μ 33 ] ] − 1 K ~ v − [ [ ε 11 ] ] − K ~ v [ [ μ 33 ] ] − 1 K ~ u ] ,
[ D z ] 12 = [ [ ε 33 ] ] 12 [ E z ] 12 .
[ D u ] 12 = [ [ [ [ 1 / ε 11 ] ] 1 − 1 ] ] 2 [ E u ] 12 .
[ D v ] 12 = [ [ [ [ ε 22 ] ] 1 − 1 ] ] 2 − 1 [ E v ] 12 .
P F R = [ K ~ u ( [ [ ε 33 ] ] 12 ) − 1 K ~ v [ [ [ [ μ 22 ] ] 1 − 1 ] ] 2 − 1 − K ~ u ( [ [ ε 33 ] ] 12 ) − 1 K ~ u K ~ v ( [ [ ε 33 ] ] 12 ) − 1 K ~ v − [ [ [ [ 1 / μ 11 ] ] 1 − 1 ] ] 2 − K ~ v ( [ [ ε 33 ] ] 12 ) − 1 K ~ u ] ,
Q F R = [ K ~ u ( [ [ μ 33 ] ] 12 ) − 1 K ~ v [ [ [ [ ε 22 ] ] 1 − 1 ] ] 2 − 1 − K ~ u ( [ [ μ 33 ] ] 12 ) − 1 K ~ u K ~ v ( [ [ μ 33 ] ] 12 ) − 1 K ~ v − [ [ [ [ 1 / ε 11 ] ] 1 − 1 ] ] 2 − K ~ v ( [ [ μ 33 ] ] 12 ) − 1 K ~ u ] .
d 2   d z ~ 2 [ s u s v ] = Ω u v 2 [ s u s v ] = W u v β u v 2 W u v − 1 [ s u s v ] ,
S u v = [ s u ( z ) s v ( z ) ] = W u v e − β u v k 0 ( z − z 1 ) c u v + + W u v e β u v k 0 ( z − z 2 ) c u v − ,
U u v = [ u u ( z ) u v ( z ) ] = − V u v e − β u v k 0 ( z − z 1 ) c u v + + V u v e β u v k 0 ( z − z 2 ) c u v − ,
Ψ u v ( z ) = [ s u ( z ) s v ( z ) u u ( z ) u v ( z ) ] = F u v [ e − β u v k 0 ( z − z 1 ) 0 0 e β u v k 0 ( z − z 2 ) ] c u v ,
[ s x s y ] = T [ s u s v ] = [ T x 0 0 T y ] [ s u s v ] .
E u ( u , v , z ) = ∑ p ∑ q s u , p q e − j [ k x , p u + k y , q v ] ,
E x ( x , y , z ) = ∑ m ∑ n s x , m n e − j [ k x , m x + k y , n y ] .
[ T x ] p q , m n = 1 Λ x Λ y ∫ Λ x ∫ Λ y g ( v ) e − j [ k x , p u + k y , q v ] + j [ k x , m x ( u ) + k y , n y ( v ) ] d u d v .
[ T y ] p q , m n = 1 Λ x Λ y ∫ Λ x ∫ Λ y f ( u ) e − j [ k x , p u + k y , q v ] + j [ k x , m x ( u ) + k y , n y ( v ) ] d u d v .
W x y e − β x y k 0 ( z − z 1 ) c x y + + W x y e β x y k 0 ( z − z 2 ) c x y − = T ( W u v e − β u v k 0 ( z − z 1 ) c u v + + W u v e β u v k 0 ( z − z 2 ) c u v − )
− V x y e − β x y k 0 ( z − z 1 ) c x y + + V x y e β x y k 0 ( z − z 2 ) c x y − = T ( − V u v e − β u v k 0 ( z − z 1 ) c u v + + V u v e β u v k 0 ( z − z 2 ) c u v − )
F x y [ e − β x y k 0 ( z − z 1 ) 0 0 e β x y k 0 ( z − z 2 ) ] c x y = [ T 0 0 T ] F u v [ e − β u v k 0 ( z − z 1 ) 0 0 e β u v k 0 ( z − z 2 ) ] c u v
F x y , i − 1 [ X x y , i − 1 0 0 I ] c x y , i − 1 = [ T 0 0 T ] F u v , i [ I 0 0 X u v , i ] c u v , i ,
[ T 0 0 T ] F u v , i [ X u v , i 0 0 I ] c u v , i = F x y , i + 1 [ I 0 0 X x y , i + 1 ] c x y , i + 1 ,
[ W 0 W 0 − V 0 V 0 ] [ c x y , i − 1 + c x y , i − 1 − ] = [ T 0 0 T ] [ W u v , i W u v , i − V u v , i V u v , i ] [ I 0 0 X u v , i ] [ c u v , i + c u v , i − ] ,
[ T 0 0 T ] [ W u v , i W u v , i − V u v , i V u v , i ] [ X u v , i 0 0 I ] [ c u v , i + c u v , i − ] = [ W 0 W 0 − V 0 V 0 ] [ c x y , i + 1 + c x y , i + 1 − ] ,
[ c x y , i − 1 − c x y , i + 1 + ] = [ S 11 i S 12 i S 21 i S 22 i ] [ c x y , i − 1 + c x y , i + 1 − ] .
S 11 ( i ) = S 22 ( i ) = ( A i − X i B i A i − 1 X i B i ) − 1 ( X i B i A i − 1 X i A i − B i ) ,
S 12 ( i ) = S 21 ( i ) = ( A i − X i B i A i − 1 X i B i ) − 1 X i ( A i − B i A i − 1 B i ) ,
Select as filters
Select Topics Cancel
Top
       
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.
Privacy | Terms of Use