1. Introduction

Diffraction gratings play an important role in numerous edge optical technologies. Most of applications require optimization procedures which often rely on direct simulation methods. When characteristic features of gratings scale to wavelengths of diffracting waves only rigorous numerical methods of the Maxwell’s equations solution can provide sufficiently accurate results. For large aperture diffractive structures the complexity of diffraction simulation methods defines whether a structure can be analyzed or not. This is why the search of computationally efficient ways to simulate various types of grating is currently of particular importance.

We proposed recently an efficient implementation of the generalized source method (GSM) [1] for grating diffraction simulation with low numerical complexity and memory consumption [2,3]. Similarly to the well-known Fourier modal method (FMM) [4] the GSM implementation implies operation in the Fourier space and slicing, nevertheless, the GSM overwhelms widely the FMM in terms of required computational resources and complexity of possibly solvable structures. For a large variety of dielectric gratings and diffractive optical elements the GSM is very probably the best choice of simulation method but for low-loss metal-dielectric structures it faces the same intrinsic difficulties as the FMM. The Fourier space poorly represents the field in metal-dielectric gratings where it is generally highly concentrated at corrugation surface, this leads to rapid increase of numerical errors.

Simulation of diffraction on metal-dielectric gratings is quite well performed by the C-method [4,5]. The C-method searches for the Fourier decomposition of modes in a curvilinear coordinates related with the grating profile. This approach demonstrates an exceptional convergence compared to the FMM, but exhibits high numerical complexity. Additionally, the C-method relies implicitly on the Rayleigh hypothesis [6], and cannot be applied, for example, to overhanging gratings or closely placed gratings of different shapes.

In this work we aim at extending the GSM applicability domain and to amend the method for simulation of metal-dielectric gratings. We reformulate the method in curvilinear coordinates, which transforms a corrugation profile to a plane interface similarly to the C-method. The new method will be referred to as GSMCC (Generalized Source Method in Curvilinear Coordinates). It will be shown that the GSMCC admits a wider class of allowable coordinate transforms in comparison with the C-method, and, since it does not use the Rayleigh hypothesis in any form, the GSMCC potentially has larger domain of applicability. The proposed method keeps the main advantage of the GSM: the linear time and memory consumption with respect to the problem complexity.

The paper is organized as follows. First, we define the diffraction problem to be solved, discuss properties of curvilinear coordinate systems needed for grating representation and review the coordinate transformation of the Maxwell’s equations. All the derivations in the article are made for the simplest case of collinear diffraction on a 1D grating. We show how to split the Maxwell’s equations in curvilinear coordinates into the “basis” and the “source” parts, and formulate the GSM in the form of an implicit self-consistent integral equation. This integral equation reduces to the system of linear algebraic equations. Numerical solution of such system is pretty the same as for the previously developed GSM [2,3], therefore, we review briefly the algorithm. Finally, we demonstrate the GSMCC numerical implementation and give an illustrative example of the enhanced transmission by a sinusoidally modulated thin metal layer.

2. Diffraction problem

To formulate the new method we start with the following diffraction problem. Consider a monochromatic plane wave with time dependence factor $\exp \left(-i\omega t\right)$ incident on a periodically corrugated interface between two semi-infinite homogeneous media (Fig. 1). The interface is defined in Cartesian coordinates $x_1,x_2,x_3$ by continuously differentiable function $x_3=\vartheta \left(x_1\right)$, with period $\Lambda$ along axis $X_1$, so that $\vartheta \left(x_1+m\Lambda \right)=\vartheta \left(x_1\right)$, $m\in Z$. Diffraction is described by the Maxwell’s equations

 

Fig. 1 Problem under consideration: a plane wave diffraction on a smooth-profiled 1D metal- dielectric grating.

Download Full Size | PDF

3. Curvilinear coordinates

The GSMCC is formulated further in a simple case of the diffraction on a 1D grating. The GSM as it was presented in [2,3] solves the diffraction problem by treating it in one-dimensional Cartesian coordinate space and two-dimensional Fourier reciprocal space. In order to adapt the method to metal-dielectric structures and to get rid of Fourier representation of the grating volume we employ the cornerstone C-method idea and rewrite Eq. (1) in a curvilinear coordinate system $z^1,z^2,z^3$ where the grating interface becomes a plane interface $z^3=0$ [Fig. 2(a)–2(b)]. For a 1D grating it is natural to take $z^2\equiv x_2$. Note that translation invariance of curvilinear coordinates (which is essential for the C-method) is not supposed here since we will not search for eigensolutions below and above plane $z^3=0$. Moreover, we introduce curvilinear coordinates in some bounded semi-infinite region Ω containing the grating $\Omega =\left\{-b\le x_3\le b|b\ge c\right\}$ and not in the whole space (here $\max f=c$, and $\min f=-c$ is assumed, without loss of generality). We assume also that new coordinates $z^1,z^2,z^3$ continuously become the Cartesian coordinates $x_1,x_2,x_3$ at Ω boundaries $\partial \Omega =\left\{x_3=\pm b\right\}$ and coincide with them everywhere outside Ω.

 

Fig. 2 Transformation of a grating corrugation profile. (a) Curvilinear coordinates introduced in a bounded region Ω = {–b ≤□x³ ≤□b | b ≥□c}. (b) Grating profile becomes a plane interface z³ = 0, and everywhere outside Ω the new coordinates coincide with the Cartesian ones.

Download Full Size | PDF

Generally, a curvilinear coordinate system allows for defining two bases. The first basis is composed of tangent vectors to curves $x_{\alpha }=x_{\alpha }\left(z^1,Z^2,Z^3\right)$, $x_{\alpha }=x_{\alpha }\left(Z^1,z^2,Z^3\right)$, and $x_{\alpha }=x_{\alpha }\left(Z^1,Z^2,z^3\right)$, $\alpha =1,2,3$ with $Z^{\beta }$ being some fixed constants. This is the contravariant basis $e_{\alpha }={\displaystyle {\sum }_{\beta }\left(\partial x_{\beta }/\partial z^{\alpha }\right)i_{\beta }}$. Symbols $i_{\beta }$, $\beta =1,2,3$ denote the Cartesian orthonormal basis vectors. The second curvilinear basis is formed by normal vectors to surfaces $z^{\alpha }=const$. It is the covariant basis $e^{\alpha }={\displaystyle {\sum }_{\beta }\left(\partial z^{\alpha }/\partial x_{\beta }\right)i_{\beta }}$. Scalar products $g^{\alpha \beta }=e^{\alpha }\cdot e^{\beta }$ and $g_{\alpha \beta }=e_{\alpha }\cdot e_{\beta }$ form the metric tensors. Covariant and contravariant components of any vector $v$ are identified by lower and upper indices, as ${\text{v}}_{\alpha }$ and ${\text{v}}^{\alpha }$, respectively. With this notation, Eq. (1) are written in the new coordinates (see, for example [4,7],)

After having discussed the coordinate transformations and retrieved Maxwell’s equations in curvilinear coordinates we can proceed to the analysis of Eq. (2) by the GSM and to formulating the new method.

4. GSM implementation in curvilinear coordinates

To transform Eq. (2) to a self-consistent integral equation we exploit the formal similarity between Eqs. (1) and (2). The core of the GSM consists in constructing the solution in a basis structure defined by some functions ${\epsilon }_b\left(z\right)$ and ${\mu }_b\left(z\right)$ (the subscript means “basis”). In principle, such basis structure can be chosen arbitrarily. For example, one can simply take some scalar constants ${\epsilon }_b$ and ${\mu }_b$. The solution of covariant Maxwell Eqs. (2) in such simple structure can be found for any given distribution of electric J(z) and magnetic M(z) sources:

Because of the structure periodicity along $z^1$axis the diffraction problem will be solved in two-dimensional Fourier space (reciprocal to axes $z^1$ and $z^2$) and one-dimensional coordinate space (axis $z^3$). Therefore, the sources are determined in the form of plane currents

At the next step of the GSM Eq. (2) is rewritten in the form (3) introducing the generalized sources:

The last step consists in transition to the Fourier space conjugate to axis $Z^1$. By definition of the discussed coordinate transformation, tensors ${\widehat{\epsilon }}^{\alpha \beta }$ and ${\widehat{\mu }}^{\alpha \beta }$ are periodic functions along axis $z^1$. We assume that in new coordinates the period still equals $\Lambda$. Thus, the Fourier decomposition of fields and matrix V components is used:

Substituting generalized sources (7) and the modified field into (5) yields the integral equation on the unknown vector of the field Fourier components:

5. Numerical method

Reduction of Eq. (12) to a linear algebraic equation system is done in two steps. First the Fourier series are truncated to some finite maximum number $N_O$ (so that$-N_O/2\le n\le N_O/2$). It is important here to notice that in this work we consider continuous functions comprising matrix V (8), i.e. a continuously differentiable grating corrugation profile is assumed and two media below and above the corrugation interface are homogeneous. In this case no problem arises with the Fourier representation of discontinuous products (as in the FMM [10–12] and GSM [2,3] implementations), and truncation of the Fourier series is mathematically justified.

Modified fields $\tilde{E}$ and $\tilde{H}$ can be represented as superposition of plane waves propagating in the basis homogeneous medium with ${\epsilon }_b$ and ${\mu }_b$. Introducing two types of TE and TM diffraction plane waves with amplitudes $a_e{}^{\pm }$ and $a_h{}^{\pm }$, respectively, the relation between the field components and plane wave amplitudes writes (see [3] for the details):

The second discretization step implies replacement of the integral in Eq. (12) by a finite sum over slices along coordinate $z^3$. Denoting slice thicknesses as $\Delta h$ and the number of slices as $N_S=2b/\Delta h$, slice center coordinates become $z_p^3=\Delta h\left(2p+1-N_S\right)/2$, $p=0,1,\dots ,N_S-1$. Exponential factors $\exp \left[\pm i{k}_{3m}\left(z^3-\zeta \right)\right]$ in Eq. (12) due to propagation of plane waves are combined in matrix $R_{mpq}^{\pm }$:

Thus, Eq. (12) is reduced to a set of linear algebraic equations on unknown amplitudes $a={\left(\begin{array}{cccc}{a}_{emp}^{+}& a_{emp}^{-}& a_{hmp}^{+}& a_{hmp}^{-}\end{array}\right)}^{\text{T}}$ of each harmonic (index n) in each slice (index p):

The curvilinear metric is treated in form of electromagnetic sources. Therefore, incident field amplitudes are written in non-transformed space as if the new coordinates coincide with the Cartesian ones. Thus, amplitudes of a plane incident wave are found in each slice by exponential propagation factors $\exp \left(i{k}_{3m}\left|z^3\right|\right)$ similarly to the standard Cartesian GSM formulation [2,3].

Equation (19) allows for calculating $4N_S{N}_O$ diffracted field amplitudes in each slice. In order to find $4N_O$ amplitudes of the diffracted field at upper and lower boundaries of region Ω $a^{out}={\left(\begin{array}{cccc}{a_{en}^{out+}|}_{z^3=b}& {a_{en}^{out-}|}_{z^3=-b}& {a_{hn}^{out+}|}_{z^3=b}& {a_{hn}^{out-}|}_{z^3=-b}\end{array}\right)}^{\text{T}}$ we apply discretized Eq. (12) once again:

Formally, Eqs. (19) and (20) is similar to the corresponding equations derived in [2,3] for diffraction in the Cartesian coordinates. Complexity of matrix inversion in (19) and (20) defines the complexity of the method. An exceptional performance of $O\left(N_O{N}_S\log \left(N_O{N}_S\right)\right)$ for diffraction simulation relies on two powerful numerical techniques: the generalized minimal residual method (GMRES) [13] and the fast Fourier transform (FFT). The GMRES is used to iteratively invert matrices in (19), (20), and requires performing only matrix-vector multiplications. When performing these multiplications we use the fact that product $\text{RPVQ}$ is a product of block-diagonal and block-Toeplitz matrices, while a product of Toeplitz matrix V by a vector can be rewritten as a convolution product by Toeplitz-extended-to-circulant matrix and calculated with the FFT. This allows performing matrix-vector multiplications required by each GMRES iteration in time which linearly grows relative to the matrix size. The memory consumption is also linear: $O\left(N_O{N}_S\right)$. A more detailed explanation of the numerical algorithm can be found in [2,3].

6. Numerical examples

The developed method was applied to calculate the diffraction on a sinusoidal grating $x_3=c\sin \left(K{x}_1\right)$. In this case, the solution can be compared with well recognized reference methods: either with the Rayleigh method [14] or with the C-method [5]. We define the following coordinate transform

 

Fig. 3 (a) Gold sinusoidal grating example and coordinate lines of the introduced curvilinear coordinate transformation. (b) Thin sinusoidal gold film placed in the air.

Download Full Size | PDF

To demonstrate the convergence consider a metallic sinusoidal grating with refractive index n_g = 0.25 + 6.25i placed on a substrate of index n_s = 1.5, and covered by air (n_c = 1) [Fig. 3(a)]. Profile parameters are c = 0.1 μm, b = 0.11 μm, and Λ = 1 μm. We calculate the diffraction of a plane wave of wavelength λ = 0.6328 μm incident under 10° from the substrate side. Figure 4(a) shows the convergence of the calculated diffraction efficiency, and comparison with the Rayleigh method versus increasing number of slices N_S with fixed number of diffraction orders N_O = 64 (this number of harmonics guarantees the accuracy of the Rayleigh method solution to be better than the single floating point precision ${10}^{-8}$). Convergence is calculated as maximum absolute difference between corresponding diffraction amplitudes for two subsequent slice numbers. The GMRES breakpoint was chosen to be ${10}^{-10}$. The graph demonstrates that for a given number of diffraction orders the GSMCC converges directly to the Rayleigh method (or alternatively to the C-method) solution. This is a quite important result since the conventional GSM exhibits poorer convergence even in the case of dielectric sinusoidal grating simulation (see the numerical examples in [2,3]). The next Fig. 4(b) shows the convergence of the method with the increasing number of diffraction orders for large number of slices N_S = 2048. It is seen that the method quite rapidly and almost monotonically converges to better than a single floating point precision. Note also that for the both TE and TM polarizations the convergence is near the same.

 

Fig. 4 Simulation results for the diffraction of a plane wave of wavelength λ = 0.6328 μm on a metal-dielectric sinusoidal grating with refractive index n_g = 0.25 + 6.25i placed onto a substrate of index n_s = 1.5, and covered by air (n_c = 1) with c = 0.1 μm, b = 0.11 μm, and Λ = 1 μm. (a) Convergence of the GSMCC and comparison of the GSMCC with the Rayleigh method with better than 10⁻⁸ accuracy. (b) Convergence of the GSMCC with the increasing number of diffraction orders for large number of slices N_S = 2048.

Download Full Size | PDF

Figure 4(a) demonstrates that solution obtained by the GSMCC and the Rayleigh method are pretty the same. To provide an additional demonstration of the method validity we compared the results for the diffraction on a sinusoidal dielectric grating (period 1 μm, depth 0.1 μm, refractive index 2.5; plane wave of wavelength 0.6328 μm is incident under angle $\theta ={10}^{\circ }$) with the solution obtained by the FMM. Corresponding diffraction efficiencies are shown in Table 1. Both the GSMCC and the FMM solutions are calculated with N_O = 63 and N_S = 1024. It is seen that solutions coincide up to the 5th digit.

Table 1. Comparison of diffraction efficiencies calculated by the GSMCC and the FMM for the diffraction of a plane wave on a sinusoidal corrugation interface separating air and a substrate of refractive index 2.5. Grating period is 1 μm, grating depth is 0.1 μm. A plane wave of wavelength 0.6328 μm is incident under 10° angle from the air side

View Table

In order to demonstrate the method abilities we calculate the enhanced transmission by a thin sinusoidal gold film (n_g = 0.25 + 6.25i) [15] [Fig. 3(b)]. For the same corrugation depth of 0.1 μm and film thickness of 0.02 μm the reflection and transmission curves of a plane wave of wavelength λ = 0.6328 μm incident under 10° angle are shown in Fig. 5 versus the grating period. The enhanced transmission phenomenon is clearly observed for a range of grating periods around 0.53 μm and around 0.76 µm. The resonances are due to the coupling of the gold film surface plasmon in the ( ± 1)^st grating orders. A detailed description of the phenomenon can be found, for example, in [16].

 

Fig. 5 Reflection and transmission curves of a plane wave of wavelength λ = 0.6328 μm incident under 10° angle on a thin sinusoidal gold film (n_g = 0.25 + 6.25i) of thickness 0.02 μm and corrugation depth 0.1 μm versus grating period.

Download Full Size | PDF

7. Conclusion

We present the integral curvilinear coordinate Fourier method for grating diffraction simulation. In this method, which we refer to as the GSMCC, advantages of the both GSM and C-method are combined together. The cornerstone idea of the C-method, namely coordinate transformation which converts a grating corrugation profile to a plane interface, was used to transform the corrugated interface into the continuous volume modification of the material properties. Such transformed diffraction problem is successfully solved by the GSM. Maxwell’s equations are split into the basis part with known analytical solution and the source part containing the curvilinear metric of permittivity and permeability modification. Contrary to the C-method, the curvilinear coordinate transformation is defined in a bounded region, and do not possess the translation invariance property which is essential in the C-method. Such transformation allows for considering closely placed gratings with profiles described by different functions. Additionally, the GSMCC does not rely on the Rayleigh hypothesis and potentially could be applicable to a larger variety of gratings in comparison with the C-method. The GSMCC is shown to have the linear numerical complexity relative to the number of discretization points as opposed to the cubic complexity of both the FMM and the C-method. To our knowledge this is the first work presenting a method with the described rationale, thus, all the derivation was made for a simple case of diffraction on 1D gratings with smooth corrugation profiles. Details of the GSMCC formulation for 2D gratings require a separate thorough explanation and this will be addressed in a future publication.