## Abstract

The coordinate transformation method (C method) is a powerful tool for modeling photonic structures with curved boundaries of discontinuities. As a modal method upon the Fourier basis, the C method has superior computational efficiency and rich physical intuitiveness compared to other full-wave numerical methods. But presently the C method is limited to two-dimensional (2D) structures if the boundaries between adjacent *z*-invariant layers are of generally different profiles [with (*x*,*y*,*z*) being the Cartesian coordinate]. Here we report a nontrivial extension of the C method to the general case of three-dimensional (3D) structures with curved boundaries of different profiles between adjacent layers. This extension drastically enlarges the applicability of the C method to the various interesting structures in nanophotonics and plasmonics. The extended 3D-C method adopts a hybrid coordinate transformation which includes not only the *z*-direction coordinate transformation in the classical C method but also the *x*- and *y*-direction matched coordinates adopted in the Fourier modal method (FMM), so as to exactly model the curved boundaries in all the three directions. The method also incorporates the perfectly matched layers (PMLs) for aperiodic structures and the adaptive spatial resolution (ASR) for enhancing the convergence. A modified numerically-stable scattering-matrix algorithm is proposed for solving the equations of boundary condition between adjacent *z*-invariant layers, which are derived via a transformation of the full 3D covariant field-components between the different curvilinear coordinate systems defined by the different-profile top and bottom boundaries of each layer. The validity of the extended 3D-C method is tested with several numerical examples.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

The Fourier modal method (FMM) [1,2], also known as the rigorous coupled wave analysis [3], is nowadays one of the most popular full-wave numerical methods for modeling periodic photonic structures. Compared with the methods that require a discretization of space in all the three *x*-, *y*- and *z*-directions [with (*x*,*y*,*z*) being the Cartesian coordinate], such as the finite-difference time-domain method (FDTD) [4] and the finite element method (FEM) [5], the FMM [6–11] as well as other modal methods [12,13] requires a discretization of space only in *x*- and *y*-directions (upon the Fourier basis for the FMM), and has the merit of field analyticity in *z*-direction in terms of waveguide eigenmodes [14,15] in different *z*-invariant layers that constitute the whole structure. The latter reduces the computational amount and provides more physical intuitiveness. The FMM also benefits from the well-developed theories and algorithms of Fourier series [16].

The FMM has experienced substantial developments in terms of its convergence, accuracy and applicability. The slow convergence of the FMM for high-contrast gratings under transverse-magnetic (TM) polarization is drastically improved [7,8]. This improvement is attributed to the use of the correct Fourier factorization rules [9], which are further applied to crossed gratings [10] and anisotropic medium [11]. For matching the boundary condition between adjacent *z*-invariant layers, the numerical instability of the transfer-matrix algorithm, which arises from the exponentially decaying waveguide modes in central layers, is removed with the use of the various scattering-matrix algorithms (see [17] for a review). The concept of adaptive spatial resolution (ASR) [18,19] is proposed to improve the slow convergence of the FMM arising from the field discontinuity at medium discontinuities. By incorporating perfectly matched layers (PMLs) to build up an artificial periodicity, the algorithm of FMM can be applied as well for aperiodic structures [20]. Matched coordinates [21,22] are proposed to avoid the zigzag approximation of curved boundaries in the transversal *x*- and *y*-directions while ensuring a legitimate use of the correct Fourier factorization rules.

For the structure with curved boundaries between adjacent *z*-invariant layers, classical FMM requires a staircase approximation of the curved boundaries [23,24], which suffers from large computational amount and low accuracy especially for metallic boundaries. For such structures, the coordinate transformation method (C method) [25–32] is proposed to map the curved boundaries to be planar boundaries and solve Maxwell’s equations under the curvilinear coordinate system with the algorithm of the FMM. Thus the C method can avoid the staircase approximation of curved boundaries but at the expense of solving more complex Maxwell’s equations. The C method is applicable to two-dimensional (2D, i.e. invariant along a transversal *y*-direction) [25] or three-dimensional (3D, i.e. varying along both transversal *x*- and *y*-directions) [26,27] structures with identical-profile boundaries between adjacent layers, and is further extended to 2D structures with different-profile boundaries between adjacent layers [28–32].

In this paper, we extend the C method to a very general case, which is the 3D structure with boundaries of different profiles between adjacent *z*-invariant layers (see Fig. 1 for a sketch). This extension is nontrivial compared to the classical C method [25–32] as explained below.

First, for the proposed 3D-C method, it is allowed to exist discontinuities in the transversal *x*- and *y*-directions for the permittivity and permeability tensors and for the derivatives of the profile functions of boundaries. These discontinuities should be treated properly so as to ensure a legitimate use of the correct Fourier factorization rules [9]. For this purpose, besides the *z*-direction coordinate transformation adopted in the classical C method [26,27], the transversal *x*- and *y*-direction transformation of matched coordinates [21,22] is further incorporated to map the curves of transversal discontinuities to be aligned with the new coordinate lines (as explained in Section 2.2). Moreover, additional transversal coordinate transformation is incorporated to perform the PMLs [20] and the ASR [18,19].

Second, for matching the field-continuity boundary conditions at the different-profile boundaries between adjacent *z*-invariant layers, substantial extension of the classical C method for 2D structures [28–32] is required when considering the full 3D field-components under the above hybrid coordinate transformation. This extension involves a transformation of the tangential covariant field-components between the different curvilinear coordinate systems defined by the different-profile top and bottom boundaries of each layer (see Section 2.3), and a numerically-stable scattering-matrix algorithm for solving the equations of boundary condition (Section 2.4).

This paper is organized as follows. The 3D-C method is introduced in Section 2 (with a flow chart at the end as a summary and a guidance for an easy reading). The validity of the method is tested by several numerical examples in Section 3, which also show the possible wide applicability of the method to 3D dielectric/plasmonic and periodic/aperiodic photonic structures. Conclusions are summarized in Section 4.

## 2. Method

#### 2.1 Specification of the structure

As sketched in Fig. 1, the considered 3D structure is composed of *L* central layers sandwiched by a bottom and a top semi-infinite layers, which are assigned numbers *l*=0, 1, …, *L*+1 from bottom to top. The medium is assumed to be generally anisotropic described by a permittivity tensor **ε**(**r**) and a permeability tensor **μ**(**r**) as functions of the Cartesian spatial coordinate **r**=(*x*,*y*,*z*). Each layer is defined to be invariant in the *z*-direction, i.e. **ε**(**r**) and **μ**(**r**) being independent of the *z*-coordinate. All layers are assumed to be periodic with the same periods Λ* _{x}* and Λ

*in the transversal*

_{y}*x*- and

*y*-directions, respectively. For a structure that is aperiodic in

*x*- or

*y*-direction, it can be mapped to be an artificial periodic structure with the use of the PML [20]. The bottom boundary of layer

*l*(

*l*=1, 2, …,

*L*+1) is assumed to be generally a curved surface described by a function

*z*=

*f*

^{(l)}(

*x*,

*y*), and

*f*

^{(l)}(

*x*,

*y*) with different

*l*are assumed to be of generally different profiles, i.e.

*f*

^{(l)}(

*x*,

*y*)−

*f*

^{(m)}(

*x*,

*y*) (with

*l*≠

*m*) being not a constant.

#### 2.2 Modal solution of electromagnetic field in each z-invariant layer

In this Section, we will consider the solution of the electromagnetic field in each *z*-invariant layer in terms of the waveguide eigenmodes [14,15]. The solution will be obtained under a curvilinear coordinate system defined in accordance with the upper or lower boundary of each *z*-invariant layer, so as to achieve an exact matching of the field-continuity boundary condition between adjacent *z*-invariant layers. This is the key merit of the C method [25–32] and can avoid the staircase approximation [23,24] used in the FMM for curved boundaries between adjacent layers.

For this purpose, the upper or lower boundary of each *z*-invariant layer is described by the function *z *= *f*(*x*,*y*) (see Fig. 1), where the superscript of layer number *l* is omitted for simplicity. Then a curvilinear coordinate system (*u*^{1},*u*^{2},*u*^{3})=(*u*,*v*,*w*) is defined as,

Note that in Eq. (1), besides adopting the coordinate transformation *z *= *w *+ *f*(*x*,*y*) in the classical C method [26,27] such that *z *= *f*(*x*,*y*) is mapped to be *w*=0, here we further introduce the matched coordinate transformation *x *= *x*(*u*,*v*) and *y *= *y*(*u*,*v*) [21,22] which we will explain after Eq. (3). The covariant basis vectors can be calculated to be (**e**_{1},**e**_{2},**e**_{3})=(**x**,**y**,**z**)**P**, where (**x**,**y**,**z**) are the Cartesian basis vectors, and the *i*-th row and *j*-th column element of **P** is (**P**)_{i}_{,j}=∂*x _{i}*/∂

*u*with (

^{j}*x*

_{1},

*x*

_{2},

*x*

_{3})=(

*x*,

*y*,

*z*). Then the contravariant basis vectors can be calculated to be (

**e**

^{1},

**e**

^{2},

**e**

^{3})=(

**x**,

**y**,

**z**)

**Q**with

*x*,

_{u}*x*,

_{v}*y*,

_{u}*y*,

_{v}*f*and

_{x}*f*are defined in the same way as

_{y}*x*=∂

_{u}*x*/∂

*u*, and

*V*=

*x*−

_{u}y_{v}*x*. The covariant form of the Maxwell’s equations is [14,33,34],

_{v}y_{u}*E*=

_{i}**e**

*·*

_{i}**E**and

*H*=

_{i}**e**

*·*

_{i}**H**(

*i*=1, 2, 3) are the covariant components of the electric vector

**E**and the magnetic vector

**H**, respectively,

*ε*

^{i}^{,j}=

**e**

*·*

^{i}**ε**·

**e**

*and*

^{j}*μ*

^{i}^{,j}=

**e**

*·*

^{i}**μ**·

**e**

*are the contravariant components of*

^{j}**ε**and

**μ**, respectively, and

*ω*is the angular frequency. Equation (3) has the same form as the Maxwell’s equations under the Cartesian coordinate system, which is the well-known form invariance of Maxwell’s equations under a curvilinear coordinate system [14,33,34], so that the algorithm of FMM for anisotropic medium [11] can be applied as well for solving Eq. (3).

To avoid the zigzag approximation [10] for treating the discontinuities of *Vε ^{i}*

^{,j}and

*Vμ*

^{i}^{,j}in solving Eq. (3) with the FMM, the matched coordinates

*x*(

*u*,

*v*) and

*y*(

*u*,

*v*) [21,22] are defined such that the discontinuities of

*Vε*

^{i}^{,j}and

*Vμ*

^{i}^{,j}are at curves of

*x*=

*x*(

*u*,

*v*) and

*y*=

*y*(

*u*,

*v*) with

*u*or

*v*being a constant [see Fig. 3(c) for an example]. The discontinuity positions of

*Vε*

^{i}^{,j}and

*Vμ*

^{i}^{,j}should include those of

**ε**(

**r**) and

**μ**(

**r**) [the black solid and dashed curves in Fig. 1(c), for example] and of

*f*and

_{x}*f*[the red dashed curves in Fig. 1(c), for example].

_{y}To incorporate the PMLs [20] for aperiodic structures and the ASR [19] for enhancing the convergence, we further introduce a coordinate transformation *u *= *G*_{1}(*u*′) and *v *= *G*_{2}(*v*′), which are maps of *R*→*C* to perform PMLs with *u*′ and *v*′ in the region of PMLs, and are maps of *R*→*R* otherwise (i.e. with *u*′ and *v*′ in the computational window truncated by the PMLs) to perform the ASR [35]. then Eq. (3) with (*u*′,*v*′,*w*) as the new differential variables should be solved after an replacement,

*g*

_{1}(

*u*′)=[

*dG*

_{1}(

*u*′)/

*du*′]

^{−1}and

*g*

_{2}(

*v*′)=[

*dG*

_{2}(

*v*′)/

*dv*′]

^{−1}are set to be continuous functions of

*u*′ and

*v*′, so that the products in the right side of Eq. (4) can be Fourier factorized with the Laurent’s rule [9].

For the electromagnetic field in the *z*-invariant layer *l* (*l*=0, 1, …, *L*+1) under the curvilinear coordinate system of Eq. (1) defined by $z = {f^{(l^{\prime})}}(x,y)$ (*l*′=*l* and *l*+1 for the bottom and top boundaries of the central layer *l*, respectively, *l*′=1 and *L*+1 for the boundaries of the bottom and top semi-infinite layers, respectively), the solution obtained with the algorithm of the FMM in [11] [which uses the correct Fourier factorization rules [9] and with a slight modification here to incorporate Eq. (4)] can be expressed as,

In Eqs. (5) and (6), **ψ**=[*E*_{1},*E*_{2},*E*_{3},*H*_{1},*H*_{2},*H*_{3}]^{T} is a column vector containing all the covariant field components, $\boldsymbol{\mathrm{\psi}}_{l, + }^{(l^{\prime})}$ and $\boldsymbol{\mathrm{\psi}}_{l, - }^{(l^{\prime})}$ represent the up-going and down-going fields, respectively, $\boldsymbol{\mathrm{\psi}}_{l, + ,r}^{(l^{\prime})}$ and $\boldsymbol{\mathrm{\psi}}_{l, - ,r}^{(l^{\prime})}$ represent the *r*-th up-going and down-going waveguide eigenmodes [14,15], respectively, $u_{l,r}^{(l^{\prime})}$ and $d_{l,r}^{(l^{\prime})}$ are the corresponding mode coefficients, and $k_{l, + ,r}^{(l^{\prime})}$ and $k_{l, - ,r}^{(l^{\prime})}$ are the corresponding mode propagation constants [satisfying $\textrm{Re} (k_{l, + ,r}^{(l^{\prime})}) + {\mathop{\rm Im}\nolimits} (k_{l, + ,r}^{(l^{\prime})}) > 0$ and $\textrm{Re} (k_{l, - ,r}^{(l^{\prime})}) + {\mathop{\rm Im}\nolimits} (k_{l, - ,r}^{(l^{\prime})}) < 0$]. The $\boldsymbol{\mathrm{\psi}}_{l, \pm ,r}^{(l^{\prime})}$ and $k_{l, \pm ,r}^{(l^{\prime})}$ are solved with the FMM, and the $u_{l,r}^{(l^{\prime})}$ and $d_{l,r}^{(l^{\prime})}$ are unknowns to be determined by matching the field-continuity boundary condition between adjacent *z*-invariant layers as to be explained in the next subsection. The $\boldsymbol{\mathrm{\psi}}_{l, \pm ,r}^{(l^{\prime})}(u^{\prime},v^{\prime})$ is expressed as a truncated Fourier series,

In Eqs. (7) and (8), [$u_0^{{\prime}}$,$u_0^{{\prime}}$+$\mathrm{\Lambda }_u^{\prime}$] and [$v_0^{{\prime}}$,$v_0^{{\prime}}$+$\mathrm{\Lambda }_v^{\prime}$] are one periods with respect to *u*′ and *v*′, respectively,$\; k_{u,m}^{\prime}$=$k_{u,0}^{\prime}$+*m*2*π*/$\mathrm{\Lambda }_u^{\prime}$ and $k_{v,p}^{\prime}$=$k_{v,0}^{\prime}$+*p*2*π*/$\mathrm{\Lambda }_v^{\prime}$, with *m* an *p* being integers, and $k_{u,0}^{\prime}$ and $k_{v,0}^{\prime}$ being the pseudo-periodic phase shift constants with respect to *u*′ and *v*′, respectively. For the structure that is periodic in *x*- or *y*-direction, there is $k_{u,0}^{\prime}$=*k _{x}*

_{,0}Λ

*/$\mathrm{\Lambda }_u^{\prime}$ or $k_{v,0}^{\prime}$=*

_{x}*k*

_{y}_{,0}Λ

*/$\mathrm{\Lambda }_v^{\prime}$, with*

_{y}*k*

_{x}_{,0}or

*k*

_{y}_{,0}being the pseudo-periodic phase shift constant with respect to

*x*or

*y*(or,

*u*or

*v*), and Λ

*or Λ*

_{x}*being the period with respect to*

_{y}*x*or

*y*(or,

*u*or

*v*[35]), respectively. There is $k_{u,0}^{\prime}$=0 or $k_{v,0}^{\prime}$=0 for the structure that is aperiodic in

*x*- or

*y*-direction, respectively [20]. Instead of using Eq. (7), an improved method is used for calculating discontinuous field components without the Gibbs oscillation [36,37].

#### 2.3 Derivation of equations of the boundary condition between adjacent z-invariant layers

The determination of the waveguide eigenmode coefficients $u_{l,r}^{(l^{\prime})}$ and $d_{l,r}^{(l^{\prime})}$ in Eq. (6) relies on a matching of the field-continuity boundary condition at the different-profile boundaries between adjacent *z*-invariant layers, for which we adopt the hybrid-spectrum method proposed in [30] for 2D structures and further extend it to 3D structures. With this method, the matching of the boundary condition at *z *= *f*^{(l)}(*x*,*y*) requires the use of the fields $\boldsymbol{\mathrm{\psi}}_{l, + }^{(l)}$ and $\boldsymbol{\mathrm{\psi}}_{l, - }^{(l + 1)}$ in layer *l* and $\boldsymbol{\mathrm{\psi}}_{l - 1, + }^{(l - 1)}$ and $\boldsymbol{\mathrm{\psi}}_{l - 1, - }^{(l)}$ in layer *l*−1 [see their definitions in Eq. (5)]. The other two fields $\boldsymbol{\mathrm{\psi}}_{l, + }^{(l + 1)}$ and $\boldsymbol{\mathrm{\psi}}_{l, - }^{(l)}$ in layer *l* and $\boldsymbol{\mathrm{\psi}}_{l - 1, + }^{(l)}$ and $\boldsymbol{\mathrm{\psi}}_{l - 1, - }^{(l - 1)}$ in layer *l*−1 are never used, which holds as well for calculating the fields in layers *l* and *l*−1. We first consider the case that *z *= *f*^{(l)}(*x*,*y*) is between two central layers, i.e. *l*=2, 3, …, *L*, and the case that *z *= *f*^{(l)}(*x*,*y*) is the boundary of the bottom or top semi-infinite layer (i.e. *l*=1 or *L*+1) will be discussed at the end of this subsection.

Due to the generally different profiles of *z *= *f*^{(l−1)}(*x*,*y*) [respectively, *z *= *f*^{(l+1)}(*x*,*y*)] and *z *= *f*^{(l)}(*x*,*y*), a transformation of $\boldsymbol{\mathrm{\psi}}_{l - 1, + }^{(l - 1)}$ [respectively, $\boldsymbol{\mathrm{\psi}}_{l, - }^{(l + 1)}$] from the curvilinear coordinate system of Eq. (1) defined by *z *= *f*^{(l−1)}(*x*,*y*) [respectively, *z *= *f*^{(l+1)}(*x*,*y*)] to that by *z *= *f*^{(l)}(*x*,*y*) is required, as explained hereafter. For the *r*-th up-going mode (*r*=1, 2, …, *N*) contained in $\boldsymbol{\mathrm{\psi}}_{l - 1, + }^{(l - 1)}$, its electric-field vector at *z *= *f*^{(l)}(*x*,*y*) can be expressed as

*i*=1, 2, 3) are elements of the $\boldsymbol{\mathrm{\psi}}_{l - 1, + ,r}^{(l - 1)}$ defined in Eq. (6a), and (

**e**

^{1},

**e**

^{2},

**e**

^{3,(l−1)}) are the (

**e**

^{1},

**e**

^{2},

**e**

^{3}) defined by Eq. (2), with

**e**

^{1}and

**e**

^{2}being shared by all the layer boundaries, and

**e**

^{3,(l−1)}being defined by

*z*=

*f*

^{(l−1)}(

*x*,

*y*). From Eq. (9), one can obtain the tangential covariant components of ${ {{\textbf E}_{l - 1, + ,r}^{(l - 1)}} |_{z = {f^{(l)}}(x,y)}}$ at

*z*=

*f*

^{(l)}(

*x*,

*y*),

In Eq. (11), the superscript ′ is used to denote the quantities that contain the exponentially decaying factor $\varphi _{l - 1, + ,r}^{(l - 1)}$ arising from mode propagation, (${\textbf e}_1^{(l)}$,${\textbf e}_2^{(l)}$) are the (**e**_{1},**e**_{2}) defined before Eq. (2) by *z *= *f*^{(l)}(*x*,*y*), $h_u^{(l - 1)} = f_u^{(l)} - f_u^{(l - 1)}$, $h_v^{(l - 1)} = f_v^{(l)} - f_v^{(l - 1)}$, $f_u^{(l)} = \partial {f^{(l)}}/\partial u$ and $f_v^{(l)} = \partial {f^{(l)}}/\partial v$ with ${f^{(l)}} = {f^{(l)}}(x(u,v),y(u,v))$. With Eq. (11), one can obtain the Fourier expansion coefficients of $E^{{\prime}{(l - 1)}}_{1,l - 1, + ,r}$ and $E^{{\prime}{(l - 1)}}_{2,l - 1, + ,r}$ as,

In Eq. (12), the superscript ∼ is used *throughout the paper* to denote the quantities of the Fourier expansion coefficients. $\tilde{{\textbf E^{\prime}}}_{1,l - 1, + ,r}^{(l - 1)}$ is a column vector with its (*m*,*p*)-th element being

*F*

_{(u′,v′)}{·}

_{(m,p)}follows the definition in Eq. (8). The quantities $\tilde{{\textbf E^{\prime}}}_{2,l - 1, + ,r}^{(l - 1)}$ and $\tilde{{\textbf E}}_{i,l - 1, + ,r}^{(l - 1)}$ (

*i*=1, 2, 3) are defined in the same way as $\tilde{{\textbf E^{\prime}}}_{1,l - 1, + ,r}^{(l - 1)}$. $\tilde{{{{\boldsymbol{\mathrm{\varphi}}}} }}_{l - 1, + ,r}^{(l - 1)}$ is a matrix with its (

*m*,

*p*)-th row and (

*n*,

*q*)-th column element being ${(\tilde{{{{\boldsymbol{\mathrm{\varphi}}}} }}_{l - 1, + ,r}^{(l - 1)})_{(m,p),(n,q)}} = {F_{(u^{\prime},v^{\prime})}}{\{ \varphi _{l - 1, + ,r}^{(l - 1)}\} _{(m - n,p - q)}}.$ $\tilde{{\textbf h}}_u^{(l - 1)}$ and $\tilde{{\textbf h}}_v^{(l - 1)}$ are defined in the same way as $\tilde{{{{\boldsymbol{\mathrm{\varphi}}}} }}_{l - 1, + ,r}^{(l - 1)}$, i.e. ${(\tilde{{\textbf h}}_\alpha ^{(l - 1)})_{(m,p),(n,q)}} = {F_{(u^{\prime},v^{\prime})}}{\{ h_\alpha ^{(l - 1)}\} _{(m - n,p - q)}}$ (

*α*=

*u*,

*v*). In Eq. (12), the $\tilde{{\textbf E}}_{i,l - 1, + ,r}^{(l - 1)}$ (

*i*=1, 2, 3) have been solved with the FMM [see Eq. (7)], and then the $\tilde{{\textbf E^{\prime}}}_{i,l - 1, + ,r}^{(l - 1)}$ (

*i*=1, 2) can be calculated with Eq. (12).

Equation (12) is derived from a Fourier factorization of Eq. (11) by applying the Laurent’s rule [9]. This is allowed in view that $\varphi _{l - 1, + ,r}^{(l - 1)}$ and $E_{3,l - 1, + ,r}^{(l - 1)}$ are both continuous with respect to *u*′ and *v*′. A direct numerical implementation of Eq. (12) could be time consuming due to its loop for all the up-going waveguide modes (indexed by *r*=1, 2, …, *N*). However, since the four matrix-vector products in Eq. (12) take the form of discrete convolution, their numerical implementation can be accelerated with the use of the fast Fourier transform algorithm [16].

With the electric-field components in Eqs. (11)-(13) replaced by the corresponding magnetic-field components, one can obtain as well the $\tilde{{\textbf H^{\prime}}}_{i,l - 1, + ,r}^{(l - 1)}$ (*i*=1, 2) that represent the counterpart of $\tilde{{\textbf E^{\prime}}}_{i,l - 1, + ,r}^{(l - 1)}$ for the magnetic field. Then according to Eq. (6a) (with replacement *l*→*l*−1 and *l*′→*l*−1), the Fourier expansion coefficients of the tangential covariant components of the electromagnetic field $\boldsymbol{\mathrm{\psi}}_{l - 1, + }^{(l - 1)}$ at *z *= *f*^{(l)}(*x*,*y*) can be expressed as,

*r*-th column of ${\textbf W^{\prime}}_{l - 1, + }^{(l - 1)}$ is

*r*-th element being $u_{l - 1,r}^{(l - 1)}$ in Eq. (6a). In Eq. (15), the semicolon represents a concatenation of several matrices along their row dimension.

Similar to Eqs. (14) and (15), the Fourier expansion coefficients of the tangential covariant components of field $\boldsymbol{\mathrm{\psi}}_{l, - }^{(l + 1)}$ at *z *= *f*^{(l)}(*x*,*y*) can be expressed as,

*l*and under the curvilinear coordinate system of Eq. (1) defined by

*z*=

*f*

^{(l+1)}(

*x*,

*y*), and ${\textbf d}_l^{(l + 1)}$ is a column vector with its

*r*-th element being the $d_{l,r}^{(l + 1)}$ in Eq. (6b). Note that the elements of ${\textbf W^{\prime}}_{l, - }^{(l + 1)}$ contain the exponentially decaying factors that arise from the propagation of down-going waveguide modes.

According to Eq. (6b) (with replacement *l*→*l*−1 and *l*′→*l*), the Fourier expansion coefficients of the tangential covariant components of field $\boldsymbol{\mathrm{\psi}}_{l - 1, - }^{(l)}$ at *z *= *f*^{(l)}(*x*,*y*) can be expressed as,

*r*-th column of ${\textbf W}_{l - 1, - }^{(l)}$ is

*l*→

*l*−1). The $\tilde{{\textbf E}}_{i,l - 1, - ,r}^{(l)}$ and $\tilde{{\textbf H}}_{i,l - 1, - ,r}^{(l)}$ (

*i*=1, 2) in Eq. (18) follow the definition of Eq. (13) but for the field components $E_{i,l - 1, - ,r}^{(l)}$ and $H_{i,l - 1, - ,r}^{(l)}$ in Eq. (6b), respectively.

Similar to Eqs. (17) and (18), the Fourier expansion coefficients of the tangential covariant components of field $\boldsymbol{\mathrm{\psi}}_{l, + }^{(l)}$ at *z *= *f*^{(l)}(*x*,*y*) can be expressed as,

*r*-th column of ${\textbf W}_{l, + }^{(l)}$ is

*l*→

*l*+1). The $\tilde{{\textbf E}}_{i,l, + ,r}^{(l)}$ and $\tilde{{\textbf H}}_{i,l, + ,r}^{(l)}$ (

*i*=1, 2) in Eq. (20) follow the definition of Eq. (13) but for the field components $E_{i,l, + ,r}^{(l)}$ and $H_{i,l, + ,r}^{(l)}$ in Eq. (6a), respectively.

Note that different from the ${\textbf W^{\prime}}_{l - 1, + }^{(l - 1)}$ and ${\textbf W^{\prime}}_{l, - }^{(l + 1)}$, the ${\textbf W}_{l, + }^{(l)}$ and ${\textbf W}_{l - 1, - }^{(l)}$ can be obtained directly from the modal solution of Eqs. (6a) and (6b) without further calculations (as shown in the flow chart of Fig. 2).

The continuity boundary condition of tangential eletromagnetic field components at *z *= *f*^{(l)}(*x*,*y*) can be expressed in terms of the Fourier expansion coefficients as,

With Eqs. (14), (16), (17) and (19) inserted into, Eq. (21) becomes,

Equation (22) is the equations of boundary condition at *z *= *f*^{(l)}(*x*,*y*) between two central layers (with *l*=2, 3, …, *L*, as sketched in Fig. 1). At the boundary of the bottom or top semi-infinite layer [i.e. *z *= *f*^{(l)}(*x*,*y*) with *l*=1 or *L*+1], the equations of boundary condition are simply Eq. (22) with the replacement

#### 2.4 Solution of equations of the boundary condition

In Eq. (22), the ${\textbf u}_0^{(0)}$ and ${\textbf d}_{L + 1}^{(L + 2)}$ [with replacement (23)] are known quantities as the excitation condition, and all the other ${\textbf u}_l^{(l)}$ and ${\textbf d}_l^{(l + 1)}$ are the unknowns to be solved. To seek the solution, we adopt the hybrid-spectrum method proposed in [30] for 2D structures and further extend it to the present 3D structures (as explained at the beginning of Section 2.3). Besides, we make some modification of the method to improve the computational efficiency as explained below.

In the framework of [17], the hybrid-spectrum method in [30] can be treated as a scattering-matrix algorithm of **W**→**s**→**S** variant, where **W**, **s** and **S** represent the eigenvector matrix [such as the 1×2 block matrices in Eq. (22)], the layer scattering matrix [which relates the waveguide mode coefficients of adjacent layers, see Eq. (35) in [30]], and the global scattering matrix [which relates the waveguide mode coefficients of each layer and the bottom layer, see Eq. (24)], respectively.

Differently, here we apply the scattering-matrix algorithm of **W**→**S** variant in [6] with some modification for solving Eq. (22). Compared with the **W**→**s**→**S** algorithm in [30], the **W**→**S** algorithm in [6] does not need the calculation of the **s** matrix and thus is more computational efficient. Concerning the modification, the exponentially decaying factors of mode propagation appear in an explicit form in [6] [see the $\phi _ + ^{(p)}$ and $\phi _ - ^{(p + 1) - 1}$ in Eq. (13.138) therein], while they are contained in the matrices ${\textbf W^{\prime}}_{l - 1, + }^{(l - 1)}$ and ${\textbf W^{\prime}}_{l, - }^{(l + 1)}$ in Eq. (22) in an implicit form. This difference arises from the fact that the bottom and top boundaries of each layer are of identical planar profiles for the former and are of different curved profiles for the latter.

To perform the above modification, the global scattering matrix **S*** _{l}* is defined to satisfy,

*l*=0, 1, …,

*L*+1. The iteration formulas for calculating

**S**

*are,*

_{l}For Eqs. (24)-(26), the replacement (23) should be performed as well. The above modified **W**→**S** algorithm is numerically stable in view that for the matrix **Z*** _{l}* that is the only matrix to be inverted in the algorithm, the stability of its numerical inversion can be ensured by the stability of the numerical inversion of the matrix $[{\textbf W}_{l, + }^{(l)}, - {\textbf W}_{l - 1, - }^{(l)}]$ even if some elements of ${\textbf W^{\prime}}_{l - 1, + }^{(l - 1)}{\textbf R}_{l - 1}^{ud}$ decay exponentially to be numerically null.

The iteration of Eqs. (25) and (26) can be started by setting **S**_{0}=**I** to be an identity matrix. With **S**_{L}_{+1} obtained, one can determine the ${\textbf d}_0^{(1)}$ and ${\textbf u}_{L + 1}^{(L + 1)}$ from the given ${\textbf u}_0^{(0)}$ and ${\textbf d}_{L + 1}^{(L + 2)}$ [with replacement (23)] according to Eq. (24), and then obtain all the electromagnetic quantities in the bottom and top semi-infinite layers (such as diffraction efficiencies of plane waves in homogeneous medium). The mode coefficients ${\textbf u}_l^{(l)}$ and ${\textbf d}_l^{(l + 1)}$ in central layers (*l*=1, 2, …, *L*) can be determined with the iteration formulas,

*l*=

*L*+1,

*L*, …, 2 so as to determine ${\textbf d}_{l - 1}^{(l)}$ from ${\textbf d}_l^{(l + 1)}$, which with Eq. (27b) then determines ${\textbf u}_{l - 1}^{(l - 1)}$.

As a summary and a guidance for an easy reading of this section, a flow chart for carrying out the proposed 3D-C method is provided in Fig. 2.

## 3. Numerical examples

In this section, four numerical examples are provided to test the validity of the proposed 3D-C method, which is implemented with an in-house software [38]. For all the examples, the structure is illuminated by a normally-incident *x*-polarized plane wave at a wavelength *λ*=0.633 μm [see Fig. 3(a)]. The relevant results are summarized in Figs. 4–7, respectively.

As sketched in Fig. 3(a), the first numerical example is a periodic array of dielectric Al_{2}O_{3} nano-spheres on a glass substrate in air. The structural parameters are specified in the caption of Fig. 4. As illustrated in Fig. 3(b), we instead consider a more general structure of a dielectric nano-cylinder (with a central-cylinder length *L _{c}*) with two ends of semi-spheres (with a radius

*R*). Thus the nano-sphere is the special case of

*L*=0. The nano-cylinder array on the substrate can be modeled with four

_{c}*z*-invariant layers [numbered by

*l*=0, 1, 2, 3 in Fig. 3(b) with boundaries shown by the red lines].

Note that if the whole surfaces of the two semi-spheres are treated as the top and bottom boundaries of the central layer 2, the derivatives *f _{x}* and

*f*of the profile function

_{y}*z*=

*f*(

*x*,

*y*) will take infinitely large values at

*z*=0 and

*L*[i.e. at

_{c}*P*

_{7},

*P*

_{8}, ${P^{\prime}_7}$ and ${P^{\prime}_8}$ in Fig. 3(b), the coordinate origin

*O*defined at the center of the bottom semi-sphere]. This will result in numerical instability in the calculation of Eq. (2). To solve this problem, we instead consider two approximate boundaries of layer 2 along with a medium perturbation of the nano-cylinder. For the approximate bottom boundary, we consider a conical surface that is tangential to the bottom sphere surface at

*P*

_{3}and

*P*

_{4}(at a small angle

*θ*) and that intersects the extended cylindrical surface at

*P*

_{2}and

*P*

_{5}, as illustrated in Fig. 3(b). Then the approximate bottom boundary of layer 2 is taken to be composed of a planar surface (

*P*

_{1}

*P*

_{2}and

*P*

_{5}

*P*

_{6}), the conical surface (

*P*

_{2}

*P*

_{3}and

*P*

_{4}

*P*

_{5}) and the sphere surface (

*P*

_{3}

*P*

_{4}). The medium in the tiny region surrounded by

*P*

_{2}

*P*

_{3}

*P*

_{7}and

*P*

_{4}

*P*

_{5}

*P*

_{8}is changed from the surrounding medium of air to be the dielectric medium of the nano-cylinder. The approximate top boundary of layer 2 and the corresponding medium perturbation of the top semi-sphere are set in a similar way, as illustrated in Fig. 3(b) with points ${P^{\prime}_1}$-${P^{\prime}_8}$. Thus, for a small value of

*θ*(

*θ*=3° for instance), the actual nano-cylinder can be well approximated by the perturbed nano-cylinder with the approximate boundaries of layer 2, for which the derivatives

*f*and

_{x}*f*of the profile function

_{y}*z*=

*f*(

*x*,

*y*) will take a finite value everywhere.

As illustrated in Fig. 3(c), the matched coordinate transformation *x *= *x*(*u*,*v*) and *y *= *y*(*u*,*v*) is defined so that *x*=±Λ* _{x}*/2,

*y*=±Λ

*/2 and the two circular boundaries of the above defined conical surface are at curves of*

_{y}*u*or

*v*being constants (shown by circles). There is (

*x*,

*y*)=(

*u*,

*v*) at the boundaries of one period, i.e. at

*x*=±Λ

*/2 or*

_{x}*y*=±Λ

*/2 [35]. The curves of*

_{y}*u*= constant and

*v*= constant are shown with vertically-oriented red curves and horizontally-oriented blue curves, respectively, and are determined with a linear interpolation between the curves of circles (more details can be found in [21,22]). The ASR for enhancing the spatial resolution and resultantly the convergence [19,35] is applied near the outer circular boundary of the conical surface, where

*f*and

_{x}*f*take large values and are discontinuous along with a discontinuity of permittivity.

_{y}For the first numerical example, its propagative orders of plane waves include the (0,0)-th order reflected plane wave, and the (0,0), (±1,0) and (0,±1)-th order transmitted plane waves, whose diffraction efficiencies are denoted by *R*_{0}, *T*_{0}, *T*_{10} and *T*_{01}, respectively. Then the absorptance can be defined as *A*=1−*R*_{0}−*T*_{0}−2*T*_{10}−2*T*_{01}, which has a theoretical value of zero for the lossless medium of the first numerical example, and thus can be used as a criteria for evaluating the accuracy of the calculated diffraction efficiencies. As shown Figs. 4(a1)-(a5), the calculated diffraction efficiencies exhibit a fast convergence with the increase of the truncated harmonic order *M _{x}*=

*M*=

_{y}*M*[see their definitions in Eq. (7)], and

*A*converges fast to the theoretical value of zero (shown by the horizontal red dashed line).

To see the performance of the 3D-C method in calculating the electromagnetic field, an enhancement factor of magnetic field is defined as EF=|*H _{y}*|

^{2}/|

**H**

_{inc}|

^{2}at (

*x*,

*y*,

*z*)=(0,0,−

*R*−

*g*) [at the gap center and the substrate surface, see the coordinate in Fig. 3(b)]. Here

**E**

_{inc}and

**H**

_{inc}denote the electric and magnetic vectors of the incident plane wave, respectively. A fast convergence of EF can be observed in Fig. 4(a6). The field distributions obtained for

*M*=30 are provided in Fig. 4(b).

The second numerical example is a periodic array of Al_{2}O_{3} nano-cylinders on a glass substrate in air. All the structural parameters are the same as those of the first example except the central-cylinder length *L _{c}*>0 [see Fig. 3(b)]. Figure 5(a) plots the afore defined

*R*

_{0},

*T*

_{0},

*T*

_{10},

*T*

_{01},

*A*and EF as functions of

*L*for

_{c}*M*=30. The calculated values of

*A*for all the values of

*L*are close to the theoretical value of zero (shown by the horizontal red dashed line). This confirms the numerical stability of the modified

_{c}**W**→

**S**algorithm proposed in Section 2.4 even for large values of the central-layer thickness [which means a strong exponential decaying of the evanescent waveguide modes, and the resultant numerically null elements in the ${\textbf W^{\prime}}_{l - 1, + }^{(l - 1)}$ in Eq. (26b)]. The field distributions obtained for

*L*=1.4 μm are provided in Fig. 5(b), showing Fabry-Perot-like standing waves along

_{c}*z*-direction. The latter is a signature of the Fabry-Perot-like interference that results in the oscillating behaviors of the curves in Figs. 5(a1)-(a4) and (a6).

The third numerical example is a periodic array of gold nano-spheres on a gold substrate in air, which is a typical nanoparticle-on-mirror plasmonic structure. Such structures have attracted great research interests in recent years due to their ease to form metallic nano-gaps down to nanometer scales where a drastic enhancement of field can be achieved [39]. Here an enhancement factor of electric field is defined as EF=|**E**|^{2}/|**E**_{inc}|^{2} at (*x*,*y*,*z*)=(35nm,0,−*R *− *g*) in the air region [see the coordinate in Fig. 3(b)]. The refractive index of gold at the illumination wavelength *λ*=0.633 μm is 0.1807 + 2.9970*i* [40]. Figure 6(a) shows that both *R*_{0} and EF converge with the increase of the truncated harmonic order *M _{x}*=

*M*=

_{y}*M*. Their converged values agree well with those obtained with the FEM implemented with the COMSOL Multiphysics software (shown by the horizontal red dashed lines). The distributions of electric field are provided in Fig. 6(b), showing a drastic enhancement of the

*E*component in the metallic nano-gap, and in good agreement with the FEM results shown in Fig. 6(c). This good agreement also confirms that the error introduced by the medium perturbation as described in Fig. 3(b) is negligible, since there is no such perturbation in the FEM calculation.

_{z}The last numerical example is a single gold nano-sphere on a gold substrate in air, which is an aperiodic structure. As stated in Section 2.2 [see Eq. (4)], PMLs in *x*- and *y*-directions are introduced as nonlinear complex coordinate transforms (maps of *R*→*C*) to build up an artificial periodicity [20,35], so as to expand the electromagnetic field upon the Fourier basis. In the computational window truncated by the PMLs, the matched coordinate transformation *x *= *x*(*u*,*v*) and *y *= *y*(*u*,*v*) [21,22] and the ASR [19] are performed in the same way as those of the previous numerical examples of periodic structures [35]. Figure 7(a) shows that both *R*_{0} and EF converge with the increase of the truncated harmonic order *M _{x}*=

*M*=

_{y}*M*. The converged value of

*R*

_{0}agrees well with the specular reflectance of the normally incident plane wave in the absence of the nano-sphere (shown by the horizontal red dashed line), which can be understood in view that the scattering effect of the finite-size nano-sphere simply vanishes in the reflectance of the infinite-extent plane wave. The converged value of EF also coincides with the value obtained with the FEM (shown by the horizontal red dashed line). As shown in Fig. 7(b), there exists a strong enhancement of electric field in the metallic nano-gap, and the results are in good agreement with the FEM results shown in Fig. 7(c).

Compared with other solvers such as FDTD [4] and FEM [5], the unique advantage of the proposed 3D-C method lies in its nature as a modal method. For instance, for the calculation in Fig. 5(a), the 3D-C method does not require a repetitive calculation of the waveguide eigenmodes in each layer [the $\boldsymbol{\mathrm{\psi}}_{l, \pm ,r}^{(l^{\prime})}$ in Eq. (6)] for different values of *L _{c}*. Besides, the computational amount does not increase with the increase of

*L*due to the waveguide mode analyticity along the

_{c}*z*-direction [14,15]. All these can ensure a higher computational efficiency compared with FDTD or FEM. Specifically, we check the computation on a workstation computer with 32GB RAM and two Intel Xeon id X5450 CPUs of 3.00 GHz and 2.99 GHz. The total computation time of Fig. 5(a) is 12.11 hours for 251 scanned points of

*L*, resulting in 173.63 s per point, which is much shorter than the usual computation time of FDTD or FEM for such 3D structures. For instance, the computation time of FEM for the result in Fig. 6(c) is 1.79 hours, with designated maximum element sizes of

_{c}*R*/5,

*R*/7 and

*λ*/13 within the gold nano-sphere, a surrounding spherical shell of thickness

*λ*/10 and the other region, respectively (all other element parameters taking default values in the software). The distance between the substrate surface and the top (or bottom) port is

*λ*/2 + 2

*R*+

*g*(or

*λ*/3). However, this comparison does not mean that the present method must be faster than FDTD or FEM for a single computation without scanning of parameters. For instance, the computation times of the present method are 0.61, 1.53, 3.07 and 6.15 hours for the result in Fig. 6(a) for

*M*=30, 35, 40 and 45 at convergence, respectively, which are comparable to the computation time of FEM. Here note that the above comparison is only a very crude one because the computation time depends not only on the algorithm itself but also on the efficiency of the code that carries out the algorithm (the 3D-C method being carried out with the code of MATLAB). For the numerical examples, the mirror symmetries of the problem about

*x*=0 and

*y*=0 have been employed to reduce the number of the unknown Fourier expansion coefficients [41]. Concerning the consumed memory, this computer with 32GB RAM can support the calculation in Fig. 6(a) for

*M*≤55 and is out of memory for

*M*≥60, and the consumed memory for the FEM calculation of Fig. 6(c) is 25.73GB.

Besides, the 3D-C method, as well as other modal methods [12,13], can solve the full scattering matrix [the **S**_{L}_{+1} in Eq. (24)] of the system via a single computation (see Section 2.4). The different columns of the full scattering matrix correspond to excitations by different modes, and all columns form a numerically complete set of basis for various excitations. For FDTD or FEM, a single computation can only solve the electromagnetic field for a specific excitation, which corresponds to one column or a specific linear combination of all columns of the full scattering matrix (for single- or multi-mode excitation, respectively). With the use of the full scattering matrix, it is possible to achieve a rigorous modeling of some difficult problems such as semi-infinite *z*-periodic structures [42] and the near-field radiative transfer between two nanostructures [43], and an efficient modeling of *z*-periodic structures with many periods [32]. The 3D-C method can also benefit from the well-developed theories of waveguide modes [14,15]. For instance, with the use of the Lorentz reciprocity theorem between waveguide modes and electric or magnetic current sources [14] to obtain the source-excited mode coefficients, and the use of the full scattering matrix, it is possible to avoid repetitive calculations in solving the Green’s dyadic of the system when scanning the position or polarization of the Dirac sources.

## 4. Conclusions

We report a nontrivial extension of the classical C method to the general 3D periodic/aperiodic photonic structures with curved boundaries. As a modal method, the extended 3D-C method not only has superior computational efficiency due to its analyticity of field in terms of waveguide eigenmodes along the invariant *z*-direction, but also can provide more physical intuitiveness compared with other full-wave numerical methods [44]. The extended 3D-C method drastically enlarges the applicability of the classical C method and may find important applications in modeling the various interesting problems in nanophotonics and plasmonics [32,39,42,43].

## Funding

National Natural Science Foundation of China (62075104, 61775105).

## Disclosures

The author declares no conflicts of interest.

## References

**1. **T. Tamir, H. C. Wang, and A. A. Oliner, “Wave propagation in sinusoidally stratified dielectric media,” IEEE Trans. Microwave Theory Tech. **12**(3), 323–335 (1964). [CrossRef]

**2. **C. B. Burcardt, “Diffraction of a plane wave at a sinusoidally stratified dielectric grating,” J. Opt. Soc. Am. **56**(11), 1502–1509 (1966). [CrossRef]

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

**4. **F. I. Baida and A. Belkhir, “Finite Difference Time Domain Method for Grating Structures,” in * Gratings: Theory and Numeric Applications*, E. Popov, ed., 2nd rev. ed. (Institut Fresnel, AMU2014), Chap. 9.

**5. **G. Demésy, F. Zolla, A. Nicolet, and B. Vial, “Finite Element Method,” in * Gratings: Theory and Numeric Applications*, E. Popov, ed., 2nd rev. ed. (Institut Fresnel, AMU2014), Chap. 5.

**6. **L. Li, “Fourier Modal Method,” in * Gratings: Theory and Numeric Applications*, E. Popov, ed., 2nd rev. ed. (Institut Fresnel, AMU2014), Chap. 13.

**7. **P. Lalanne and G. M. Morris, “Highly improved convergence of the coupled-wave method for TM polarization,” J. Opt. Soc. Am. A **13**(4), 779–784 (1996). [CrossRef]

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

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

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

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

**12. **P. Lalanne and J. P. Hugonin, “Numerical performance of finite-difference modal methods for the electromagnetic analysis of one-dimensional lamellar gratings,” J. Opt. Soc. Am. A **17**(6), 1033–1042 (2000). [CrossRef]

**13. **G. Granet, M. H. Randriamihaja, and K. Raniriharinosy, “Polynomial modal analysis of slanted lamellar gratings,” J. Opt. Soc. Am. A **34**(6), 975–982 (2017). [CrossRef]

**14. **C. Vassallo, * Optical Waveguide Concepts* (Elsevier, 1991), pp. 9–26, 256-257.

**15. **Here the term of waveguide eigenmode/mode is used with a general sense, i.e., the eigenmode with a *z*-dependence of exp(*ikz*) for any *z*-invariant structure (including the homogeneous free space).

**16. **D. Sundararajan, * The Discrete Fourier Transform: Theory, Algorithms and Applications* (World Scientific, 2001).

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

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

**19. **B. Guizal, H. Yala, and D. Felbacq, “Reformulation of the eigenvalue problem in the Fourier modal method with spatial adaptive resolution,” Opt. Lett. **34**(18), 2790–2792 (2009). [CrossRef]

**20. **J. P. Hugonin and P. Lalanne, “Perfectly matched layers as nonlinear coordinate transforms: a generalized formalization,” J. Opt. Soc. Am. A **22**(9), 1844–1849 (2005). [CrossRef]

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

**22. **J. Küchenmeister, T. Zebrowski, and K. Busch, “A construction guide to analytically generated meshes for the Fourier modal method,” Opt. Express **20**(16), 17319–17347 (2012). [CrossRef]

**23. **E. Popov, M. Nevière, B. Gralak, and G. Tayeb, “Staircase approximation validity for arbitrary-shaped gratings,” J. Opt. Soc. Am. A **19**(1), 33–42 (2002). [CrossRef]

**24. **I. Gushchin and A. V. Tishchenko, “Fourier modal method for relief gratings with oblique boundary conditions,” J. Opt. Soc. Am. A **27**(7), 1575–1583 (2010). [CrossRef]

**25. **J. Chandezon, M. T. Dupuis, G. Cornet, and D. Maystre, “Multicoated gratings: a differential formalism applicable in the entire optical region,” J. Opt. Soc. Am. **72**(7), 839–846 (1982). [CrossRef]

**26. **G. Granet, “Analysis of diffraction by surface-relief crossed gratings with use of the Chandezon method: application to multilayer crossed gratings,” J. Opt. Soc. Am. A **15**(5), 1121–1131 (1998). [CrossRef]

**27. **J. B. Harris, T. W. Preist, J. R. Sambles, R. N. Thorpe, and R. A. Watts, “Optical response of bigratings,” J. Opt. Soc. Am. A **13**(10), 2041–2049 (1996). [CrossRef]

**28. **G. Granet, J. P. Plumey, and J. Chandezon, “Scattering by a periodically corrugated dielectric layer with non-identical faces,” Pure Appl. Opt. **4**(1), 1–5 (1995). [CrossRef]

**29. **T. W. Preist, N. P. K. Cotter, and J. R. Sambles, “Periodic multilayer gratings of arbitrary shape,” J. Opt. Soc. Am. A **12**(8), 1740–1748 (1995). [CrossRef]

**30. **L. Li, G. Granet, J. P. Plumey, and J. Chandezon, “Some topics in extending the C method to multilayer gratings of different profiles,” Pure Appl. Opt. **5**(2), 141–156 (1996). [CrossRef]

**31. **L. Li, “Oblique-coordinate-system-based Chandezon method for modeling one-dimensionally periodic, multilayer, inhomogeneous, anisotropic gratings,” J. Opt. Soc. Am. A **16**(10), 2521–2531 (1999). [CrossRef]

**32. **T. Vallius and M. Kuittinen, “Novel electromagnetic approach to photonic crystals with use of the C method,” J. Opt. Soc. Am. A **20**(1), 85–91 (2003). [CrossRef]

**33. **C. T. Tai, * Generalized Vector and Dyadic Analysis: Applied Mathematics in Field Theory*, 2nd ed. (IEEE, 1997), pp. 58–98.

**34. **M. Farhat, P. Y. Chen, S. Guenneau, and S. Enoch, eds., * Transformation Wave Physics: Electromagnetics, Elastodynamics, and Thermodynamics* (Pan Stanford, 2016).

**35. **For more details about the coordinate transformation, at *u *= *u _{c}* that may represent the discontinuity positions of

*Vε*

^{i}^{,j}and

*Vμ*

^{i}^{,j}, the ASR is performed with

*u*=

*G*

_{1}(

*u*′) that maps a small change of

*u*to be a big change of

*u*′, and the

*du*′/

*du*=

*g*

_{1}(

*u*′) that characterizes the spatial resolution increases monotonously (from the order of 1 to the order of 100 for the numerical examples) with the decrease of |

*u*−

*u*| (from 0.1

_{c}*λ*to 0 with

*λ*being the wavelength). The detailed form of

*u*=

*G*

_{1}(

*u*′) for the ASR is slightly different from that in [18,19] but follows its key features. As explained by Eq. (4) and the following statement, the way of performing the ASR follows that in [19] which may exhibit better numerical stability than that in [18]. For performing the PML,

*u*=

*G*

_{1}(

*u*′) adopts the complex coordinate transform in [20] with the parameters of Fig. 3(d) therein, for which the

*g*

_{1}(

*u*′) is a continuous function. The width

*q*of the PML region (i.e. the

_{u}*q*in [20]) takes one wavelength. Concerning the connection of

*u*=

*G*

_{1}(

*u*′) between the PML region and the ASR region, let

*u*′=$u_0^{\prime}$+

*q*/2 and $u_0^{\prime}$+$\mathrm{\Lambda }_u^{\prime}$−

_{u}*q*/2 denote the boundary points between the PML region and the ASR region (with

_{u}*u*′=$u_0^{\prime}$ and $u_0^{\prime}$+$\mathrm{\Lambda }_u^{\prime}$ denoting the starting and end coordinates of one period with respect to

*u*′), and then the connection fulfils

*G*

_{1}(

*u*′)=

*u*′ at

*u*′=$u_0^{\prime}$+

*q*/2 and

_{u}*g*

_{1}(

*u*′) = 1 at

*u*′=$u_0^{\prime}$+

*q*/2 and $u_0^{\prime}$+$\mathrm{\Lambda }_u^{\prime}$−

_{u}*q*/2. The

_{u}*v*=

*G*

_{2}(

*v*′) is performed in the same way as

*u*=

*G*

_{1}(

*u*′) but in the

*y*-direction. For the matched coordinate transformation

*x*=

*x*(

*u*,

*v*) and

*y*=

*y*(

*u*,

*v*), there is (

*x*,

*y*)=(

*u*,

*v*) out of the ASR region (which is the computational window truncated by the PMLs). For the structure that is periodic in

*x*- or

*y*-direction, above specifications still hold simply by setting

*q*=0 or

_{u}*q*=0 (

_{v}*q*being the counterpart of

_{v}*q*in

_{u}*y*-direction), and there is (

*x*,

*y*)=(

*u*,

*v*) at the boundaries of one period in

*x*- or

*y*-direction, respectively.

**36. **P. Lalanne and M. P. Jurek, “Computation of the near-field pattern with the coupled-wave method for transverse magnetic polarization,” J. Mod. Opt. **45**(7), 1357–1374 (1998). [CrossRef]

**37. **The method is an extension of the method in [36] for one-dimensionally periodic gratings to the present two-dimensionally periodic structures. Its details are out of the main scope of the present work and may be provided in future publications.

**38. **H. Liu, * DIF CODE for Modeling Light Diffraction in Nanostructures* (Nankai University, 2010).

**39. **G. C. Li, Q. Zhang, S. A. Maier, and D. Lei, “Plasmonic particle-on-film nanocavities: a versatile platform for plasmon-enhanced spectroscopy and photochemistry,” Nanophotonics **7**(12), 1865–1889 (2018). [CrossRef]

**40. **E. D. Palik, * Handbook of Optical Constants of Solids*, Part II (Academic, (1985), pp. 290–295.

**41. **B. Bai and L. Li, “Group-theoretic approach to enhancing the Fourier modal method for crossed gratings with one or two reflection symmetries,” J. Opt. A: Pure Appl. Opt. **7**(7), 271–278 (2005). [CrossRef]

**42. **G. Lecamp, J. P. Hugonin, and P. Lalanne, “Theoretical and computational concepts for periodic optical waveguides,” Opt. Express **15**(18), 11042–11060 (2007). [CrossRef]

**43. **Y. Yang and L. Wang, “Spectrally enhancing near-field radiative transfer between metallic gratings by exciting magnetic polaritons in nanometric vacuum gaps,” Phys. Rev. Lett. **117**(4), 044301 (2016). [CrossRef]

**44. **E. Popov, ed., 2nd rev. ed., * Gratings: Theory and Numeric Applications* (Institut Fresnel, AMU 2014).