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.

 

Fig. 1. Sketch of the studied 3D multilayer structure with boundaries of different profiles between adjacent z-invariant layers. (a)-(c) Views from the side and the sections A-A and B-B, respectively. Different colors represent different regions where permittivity and permeability tensors are of constant values (or more generally, continuous in x- and y-directions and invariant in z-direction). In (a)-(b), the curved boundaries between adjacent layers are shown with the red-solid curves. In (b)-(c), the projections of the discontinuity positions of the medium (or the derivatives of the layer profile functions) are shown by the black solid and dashed (or red dashed) curves. In the figure, the boundaries of blue, yellow, brown and pink regions [from top to bottom in (a)] are assumed to be revolutionary surfaces as examples. The quantities in the figure are defined in the main text [see Section 2.1 and Eq. (22)].

Download Full Size | PDF

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 Λ_y in the transversal 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¹,u²,u³)=(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₁,e₂,e₃)=(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^j with (x₁,x₂,x₃)=(x,y,z). Then the contravariant basis vectors can be calculated to be (e¹,e²,e³)=(x,y,z)Q with

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_x and f_y [the red dashed curves in Fig. 1(c), for example].

To incorporate the PMLs [20] for aperiodic structures and the ASR [19] for enhancing the convergence, we further introduce a coordinate transformation u = G₁(u′) and v = G₂(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,

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₁,E₂,E₃,H₁,H₂,H₃]^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}$+m2π/$\mathrm{\Lambda }_u^{\prime}$ and $k_{v,p}^{\prime}$=$k_{v,0}^{\prime}$+p2π/$\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Λ_x/$\mathrm{\Lambda }_u^{\prime}$ or $k_{v,0}^{\prime}$=k_y,0Λ_y/$\mathrm{\Lambda }_v^{\prime}$, with 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 Λ_x or Λ_y being the period with respect to 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

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₁,e₂) 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

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,

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,

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,

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,

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

 

Fig. 2. A flow chart for carrying out the proposed 3D-C method.

Download Full Size | PDF

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,

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₀=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,

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.

 

Fig. 3. (a) Structure of the numerical examples, which is a bi-periodic array of spheres on a flat substrate in air. (b) Definition of layers l=0-3 for the structure in (a) [L_c=0 for (a)]. (c) Curves of u = constants (red) or v = constants (blue) for the matched coordinate transformation x = x(u,v) and y = y(u,v) as explained in the main text [Λ_x=Λ_y=6μm, R=2μm and θ=60° for (c)].

Download Full Size | PDF

 

Fig. 4. Results for the numerical example of a periodic array of Al₂O₃ nano-spheres on a glass substrate in air. The refractive indices of Al₂O₃ and glass at the illumination wavelength λ=0.633μm take values 1.7662 and 1.4570, respectively. Other parameters are Λ_x=Λ_y=600 nm, R=150 nm, θ=3°, g=10 nm (see Fig. 3 for their definitions). (a1)-(a6) Convergence curves of R₀, T₀, T₁₀, T₀₁, A and EF (defined in the main text) plotted as functions of truncated harmonic order M_x=M_y=M. (b1)-(b3) Distributions of |E_x| at y=0, |H_y| at y=0 and |H_y| at z = −R − g obtained for M=30, with |E_inc| and |H_inc| normalized to be 1. The superimposed black-white lines in (b1)-(b2) show the discontinuous boundaries of the sphere and the substrate.

Download Full Size | PDF

 

Fig. 5. Results for the numerical example of a periodic array of Al₂O₃ nano-cylinders on a glass substrate in air. All the structural parameters are the same as those in Fig. 4 except the central-cylinder length L_c>0. (a1)-(a6) R₀, T₀, T₁₀, T₀₁, A and EF plotted as functions of L_c for M=30. (b1)-(b3) Distributions of |E_x| at y=0, |H_y| at y=0 and |H_y| at z = −R − g obtained for M=30 and L_c=1.4 μm, with |E_inc| and |H_inc| normalized to be 1.

Download Full Size | PDF

 

Fig. 6. Results for the numerical example of a periodic array of gold nano-spheres on a gold substrate in air. The structural parameters are Λ_x=Λ_y=550 nm, R=75 nm, θ=3°, g=10 nm. (a1)-(a2) Convergence curves (blue) of R₀ and EF defined in the main text. The horizontal red dashed lines show the results obtained with the FEM. (b1)-(b3) Distributions of |E_x| at y=0, |E_z| at y=0 and |E_z| at z = −R − g (in the air region) obtained for M=45 and |E_inc|=1 with the proposed 3D-C method. (c1)-(c3) The same as (b1)-(b3) but obtained with the FEM.

Download Full Size | PDF

 

Fig. 7. Results for the numerical example of a single gold nano-sphere on a gold substrate in air. The values of R, θ, g and λ are the same as those in Fig. 6. The width of the computational window truncated by the PMLs is 2R+λ in both x- and y-directions. (a1)-(a2) Convergence curves (blue) of R₀ and EF (following their definitions in Fig. 6). The horizontal red dashed lines in (a1) and (a2) show the theoretical value of R₀ and the value of EF obtained with the FEM, respectively. (b1)-(b3) Distributions of |E_x| at y=0, |E_z| at y=0 and |E_z| at z = −R − g (in the air region) obtained for M=50 and |E_inc|=1 with the proposed 3D-C method. (c1)-(c3) The same as (b1)-(b3) but obtained with the FEM.

Download Full Size | PDF

As sketched in Fig. 3(a), the first numerical example is a periodic array of dielectric Al₂O₃ 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_c=0. The nano-cylinder array on the substrate can be modeled with four 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_y of the profile function z = f(x,y) will take infinitely large values at z=0 and L_c [i.e. at P₇, P₈, ${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₃ and P₄ (at a small angle θ) and that intersects the extended cylindrical surface at P₂ and P₅, as illustrated in Fig. 3(b). Then the approximate bottom boundary of layer 2 is taken to be composed of a planar surface (P₁P₂ and P₅P₆), the conical surface (P₂P₃ and P₄P₅) and the sphere surface (P₃P₄). The medium in the tiny region surrounded by P₂P₃P₇ and P₄P₅P₈ 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_x and f_y of the profile function 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=±Λ_y/2 and the two circular boundaries of the above defined conical surface are at curves of u or v being constants (shown by circles). There is (x,y)=(u,v) at the boundaries of one period, i.e. at x=±Λ_x/2 or y=±Λ_y/2 [35]. The curves of 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_x and f_y take large values and are discontinuous along with a discontinuity of permittivity.

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₀, T₀, T₁₀ and T₀₁, respectively. Then the absorptance can be defined as A=1−R₀−T₀−2T₁₀−2T₀₁, 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|²/|H_inc|² 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₂O₃ 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₀, T₀, T₁₀, T₀₁, A and EF as functions of L_c for M=30. The calculated values of A for all the values of L_c are close to the theoretical value of zero (shown by the horizontal red dashed line). This confirms the numerical stability of the modified 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_c=1.4 μm are provided in Fig. 5(b), showing Fabry-Perot-like standing waves along 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|²/|E_inc|² 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.9970i [40]. Figure 6(a) shows that both R₀ 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_z 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.

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₀ and EF converge with the increase of the truncated harmonic order M_x=M_y=M. The converged value of R₀ 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_c due to the waveguide mode analyticity along the 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_c, 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 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 + 2R + 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].