Accurate near-field calculation in the rigorous coupled-wave analysis method

The study of the interaction of electromagnetic waves with matter has spawned a large variety of methods to analytically or numerically solve the Maxwell equations (MEs). The results obtained from those methods are invaluable for understanding, validating, predicting, and guiding experimental efforts and for the design process of electromagnetic and optical devices. Amongst the numerical methods used in computational electromagnetics, one can fundamentally distinguish between time-domain methods, which directly incorporate the transient behavior of electromagnetic waves, such as the finite-difference time-domain method [1], and frequency-domain methods, like the finite-element frequency-domain method [2], which directly determine time-harmonic solutions to MEs. The former methods are general purpose methods and are capable of simulating virtually any electromagnetic structure comprising metallic or dielectric objects of arbitrary size and shape. However, for certain structures and applications, other methods can be superior in terms of runtime and accuracy. For example, the multiple scattering method [3] is a series expansion method in the frequency domain that is tailored to calculate interaction of light and clusters of spherical particles [4] or cylindrical rods [5], and as such it is superior to other methods when applied to these particular geometries.

An important field of applications, where several specialized numerical algorithms exist [6], is the study of periodic optical structures, including diffraction gratings [7] and periodic metamaterials [8]. Amongst the numerical methods for periodic structures, such as the differential theory of gratings [9], the ${\mathcal{C}}$-method [10], integral methods [6, 11], and the generalized source method [12, 13], the most widely used method is the rigorous coupled-wave analysis (RCWA) [14–16]. It is a modal method in the frequency domain and is based on the decomposition of the periodic structure and the pseudo-periodic solution of MEs in terms of their Fourier series (FS) expansion, hence the periodicity is naturally incorporated into the numerical method. RCWA was initially developed for modeling one-dimensional (1D) diffraction gratings, but with the introduction of fast converging Fourier factorization rules for modeling 1D [17, 18] and two-dimensional (2D) [19–21] structures, its formulation for isotropic and anisotropic materials [22] as well as multilayered [23] and oblique structures [24], RCWA has evolved to describe arbitrary, 2D-periodic structures. The method has been successfully applied to model diffraction gratings, diffractive optical elements, surface coatings, spectroscopic applications, photonic crystals, and periodic metamaterials. It should be noted, that in most cases the functionality of these applications of periodic structures depends on the electromagnetic far-field, i.e. the propagating diffraction orders, and it is known that the RCWA delivers fast converging and accurate far-field results, as reported in many works [14–24].

There is, however, a range of novel applications, which rely on the optical near-field of a periodic structure, especially at their surface, such as surface-enhanced Raman spectroscopy [25], surface second-harmonic generation [26–28], and near-field sensing [29–32]. These applications require a suitably designed near-field distribution, usually optimized for maximum field enhancement within specific spatial domains. Although there have been significant advances in experimental optical near-field measurement techniques [33], these techniques are still in their development state and not readily available to accurately characterize complex photonic nanostructures. These applications and experimental shortcomings lead to a critical demand for numerical methods for periodic structures that can facilitate an accurate calculation of electromagnetic near-fields, and more importantly, the design of gratings with optimized near-field patterns. With very few exceptions [20, 34, 35], a thorough investigation of numerically calculated near-fields in the RCWA has been largely neglected during the development of the method, as it is often merely considered a post-processing step. Moreover, additional reasons for the scarcity of reports on the convergence and the accuracy of the numerically computed electromagnetic near-field in RCWA, are the slow near-field convergence and spurious oscillations displayed by these fields [35].

In this paper we address the issue of inaccurate near-field calculations in the RCWA in several ways. First, we use some generic cases of diffraction gratings to illustrate the slow convergence of RCWA for near-field calculations and reveal the reasons for this behavior. Based on this analysis and the continuity properties of the electromagnetic fields, a new formulation of the field evaluation is proposed, which yields faster convergence of the electromagnetic near-fields and explicitly fulfills the continuity properties of the electromagnetic field at material interfaces. This improved formulation of the numerical evaluation of the near-fields is then benchmarked against the current formulation, with the aim of making RCWA an effective numerical method for modeling modern nanophotonic applications that rely on highly accurate near-field calculations.

The remainder of the paper is organized as follows: section 2 will introduce the reader to the problem of inaccurate/spurious near-fields in RCWA for 1D structures and explain the overall strategy for the accurate field evaluation. Section 3 will extend the ideas gained from the analysis of 1D structures to arbitrary, straight or obliquely etched 2D gratings and will present the underlying mathematical formalism. Then in section 4, computational results for 2D gratings will be presented and discussed. Section 5 will investigate near- and far-field convergence of the modified method for slanted gratings, before final conclusions about the capability of the improved RCWA for modeling electromagnetic near-fields will be drawn in section 6.

Although the last major roadblock that was precluding RCWA from becoming a highly effective method for modeling 1D periodic structures has been removed when the correct Fourier factorization rules for transverse magnetic (TM) polarization were introduced [17, 18], a few topics have continued to attract attention, namely the convergence for slanted 1D-periodic gratings [24] and the numerical instabilities associated to highly conductive gratings [36]. These improvements and refinements of the method are concerned with the accuracy of the far-field calculation and they achieve excellent convergence results for the coefficients of and power carried by the diffraction orders in the cover and substrate regions of a grating, namely the far-field. This picture changes if the near-field is considered. Thus, in this section we investigate the numerically calculated near-field inside the periodic grating region via RCWA, and explain its slow convergence and the origin of the observed high spatial frequency oscillations. Our approach is based on a numerical method introduced in [34], for improved near-field accuracy and reduced numerical artifacts for straight 1D-periodic gratings. This numerical approach is thoroughly investigated here and will form the basis of the general formalism for accurate near-field evaluation in arbitrary, straight or slanted, 1D- and 2D-periodic gratings proposed in the subsequent sections.

The ideas developed in this paper are applicable to any periodic grating structure with sharp boundaries, but in order to conduct a comprehensive assessment of the accuracy of the near-field calculations achievable with the improved RCWA introduced in this work, a number of generic test structures (1D- and 2D-periodic, straight and slanted) have been chosen, as per figure 1. To widen the spectrum of test configurations, three different materials are considered for the grating: silica (SiO₂) as a dielectric with low index of refraction (${n}_{g}={n}_{{\mathrm{SiO}}_{2}}=1.45$), silicon (Si) as a high refractive index dielectric (n_g = n_Si = 3.4) and a metal, gold (Au), with index of refraction n_g = n_Au = 0.97 + 1.87i [37]. In all examples the substrate is SiO₂ (${n}_{s}={n}_{{\mathrm{SiO}}_{2}}$). Normal incidence and TM-polarization is considered in all of the sections, i.e. the incident electric field amplitude is oriented along the x-axis: ${{\bf{E}}}_{\mathrm{inc}}({\bf{x}})={(1,\;0,\;0)}^{T}{E}_{\mathrm{inc}}({\bf{x}}).$

Zoom In Zoom Out Reset image size

The example considered in this section is depicted in figure 1(a). It consists of a binary grating with period, Λ₁ = 1 μm, height, h = 0.25 μm, and filling factor, ρ = 0.5. The grating is illuminated by a normally incident plane wave with wavelength, λ = 0.51 μm.

As the first step of our investigation, we used RCWA to calculate the far-field. In order to quantify the convergence of the RCWA method, the relative error of the far-field is defined as:

$$\begin{eqnarray}&&{\text{}}{e}_{{\rm{F}}}(N)=\displaystyle \frac{\sqrt{{| {T}^{(N)}-{T}^{\mathrm{ref}}| }^{2}+{| {R}^{(N)}-{R}^{\mathrm{ref}}| }^{2}}}{\sqrt{{| {T}^{\mathrm{ref}}| }^{2}+{| {R}^{\mathrm{ref}}| }^{2}}},\end{eqnarray} \tag{ 1 }$$

where T^(N) (R^(N)) denotes the relative transmitted (reflected) power corresponding to a discretization with 2N + 1 complex FS coefficients. The discretization parameter, N, represents the number of harmonics retained for each dimension. Moreover, T^ref(R^ref) is a reference value that is considered to be the exact solution or a sufficiently good approximation of the exact solution. Due to the absence of the exact solution of the diffraction grating problem, the reference values are chosen to be numerical values obtained by high-N simulations; in our case we chose T^ref = T⁽⁹⁰⁵⁾ and R^ref = R⁽⁹⁰⁵⁾.

The far-field relative error, ${\text{}}{e}_{{\rm{F}}}(N),$ for increasing number of harmonics, N = 5, ..., 640, is depicted in figure 2. As this figure illustrates, RCWA converges quickly for all three materials. Specifically, in order to achieve a self-error of ${\text{}}{e}_{{\rm{F}}}(N)\lt 1\%$ (as a generic criterion adopted here for an accurate far-field calculation), for the three materials, silica (red crosses), silicon (blue stars), and gold (green triangles), a relatively small number of harmonics is necessary, namely N > 5, N > 25, and N > 13, respectively.

Zoom In Zoom Out Reset image size

As mentioned in the introductory section, the far-field of optical gratings and periodic structures has been the physical quantity of most interest from experimental point of view, hence the characterization of a numerical method by means of the far-field convergence has usually been the adopted strategy. This approach, however, largely neglects the electromagnetic near-field predicted by a specific method. On the other hand, the near-field is of fundamental importance for modeling plasmonic effects or optical nonlinear phenomena in devices with size comparable to or smaller than the operating wavelength, effects whose description relies on accurate calculations of the electromagnetic near-field.

In order to characterize the numerically obtained near-fields, we define the grating norm, ${\| \cdot \| }_{G}$, of a scalar or vectorial function, f, in the grating region as follows:

$$\begin{eqnarray}&&{\| f\| }_{G}={\left({\displaystyle \int }_{0}^{h}{\displaystyle \int }_{-{{\rm{\Lambda }}}_{1}/2}^{{{\rm{\Lambda }}}_{1}/2}{| f(x,z)| }^{2}\ {\rm{d}}x\ {\rm{d}}z\right)}^{\frac{1}{2}},\end{eqnarray} \tag{ 2 }$$

where the z-integration extends over the bulk of the periodic region. The grating norm is used to define the near-field error, ${\text{}}e\{{E}_{\alpha }^{(N)}\},$ of the scalar field components E^(N)_α, α = x, y, z, of a near-field, ${{\bf{E}}}^{N}$, numerically obtained using N harmonics:

$$\begin{eqnarray}&&{\text{}}e\{{E}_{\alpha }^{(N)}\}={\| {E}_{\alpha }^{(N)}-{E}_{\alpha }^{\mathrm{ref}}\| }_{G}/{\| {E}_{\alpha }^{\mathrm{ref}}\| }_{G}.\end{eqnarray} \tag{ 3 }$$

Here, E_α^(N) denotes the FS reconstruction given by the 2N + 1 central Fourier coefficients, E_αn, which are calculated by RCWA:

$$\begin{eqnarray*}&&{E}_{\alpha }^{(N)}(x)=\displaystyle \sum _{n=-N}^{N}{E}_{\alpha {\rm{n}}}\mathrm{exp}\left({\rm{i}}n\displaystyle \frac{2\pi }{{{\rm{\Lambda }}}_{1}}x\right).\end{eqnarray*}$$

Similarly to the far-field calculations, ${{\bf{E}}}^{\mathrm{ref}}={{\bf{E}}}^{(905)}$ is obtained by a high-resolution RCWA simulation with N = 905. The near-field self-convergence for the tangential component, E_t = E_z, and the normal component, E_n = E_x, of the electric field is depicted in figure 3. For all three materials, the self-convergence of E_t (blue circles) is fast and comparable to the far-field self-convergence (see figure 2). The normal component, E_n (red crosses), however, exhibits much slower convergence and even at the highest numerical resolution of N = 640 a relative error of ${\text{}}e\{{E}_{x}^{(640)}\}\gt 0.9\%$ still remains, whereas the error of the tangential field, ${\text{}}e\{{E}_{z}^{(640)}\}\lt 0.08\%,$ for all materials.

Zoom In Zoom Out Reset image size

One intriguing question raised by the data plotted in figures 2 and 3 is why the normal component of the near-field converges much more slowly than the far-field power and the tangential component of the near-field. Or put it the other way around: How can the far-field converge quickly when the near-fields have not converged yet? There are two factors that explain this behavior: (i) the far-field consists of a superposition of a small number of propagating diffraction orders, i.e., plane waves. Hence, the far-field requires a small number of FS components to be reconstructed and does not suffer from the Gibbs phenomenon. (ii) RCWA does not depend on the representation of the discontinuous normal field E_n. Instead, the method relies on the correct Fourier factorization of the continuous normal component of the displacement field, D_n, which can be accurately described by a FS, i.e. without spurious oscillations. The second observation is at the core of the accurate near-field calculation that is introduced in this work. Specifically, instead of reconstructing a discontinuous physical quantity, i.e. E_n$({\bf{x}})$, directly from its FS coefficients, it is more effective to reconstruct a continuous quantity, the normal component of the displacement field, D_n$({\bf{x}})$, and then divide it by a discontinuous quantity, the electric permittivity, ₀$({\bf{x}}),$ which is known in the space-domain where one seeks to solve the diffraction problem.

For 1D-periodic gratings, we define the modified normal component of the electric field:

$$\begin{eqnarray}&&{\tilde{E}}_{{\rm{n}}}^{(N)}({\bf{x}})={\varepsilon }_{0}^{-1}{D}_{{\rm{n}}}^{(N)}({\bf{x}})/\varepsilon ({\bf{x}}),\end{eqnarray} \tag{ 4 }$$

where ₀ is the vacuum permittivity and $({\bf{x}})$ the relative permittivity at position, ${\bf{x}}.$ Note that ${E}_{{\rm{n}}}^{(N)}({\bf{x}})$ and ${\tilde{E}}_{{\rm{n}}}^{(N)}({\bf{x}})$ represent the same physical quantity, namely the normal component of the electric field. However, whereas the former is found by using RCWA to solve directly for the electric field, the latter one is determined by first calculating the displacement field and then the electric field via equation (4).

The error, ${\text{}}e\{{\tilde{E}}_{{\rm{n}}}^{(N)}\},$ of the modified normal component, ${\tilde{E}}_{{\rm{n}}},$ is shown in figure 3 and encoded in green triangle. It is found to converge as fast as the fast-convergent tangential component, E_t, and as fast as the power in the far-field shown in figure 2. Even at the highest considered resolution, N = 640, the conventional formulation of the RCWA yields a self error of ${\text{}}e\{{E}_{{\rm{n}}}^{(640)}\}\gt 9\times {10}^{-3}$ for all three materials. By contrast, the same self error ${\text{}}e\{{\tilde{E}}_{{\rm{n}}}^{(N)}\}\lt 9\times {10}^{-3}$ of the modified normal field ${\tilde{E}}_{{\rm{n}}}$ can be achieved by using as few as N = 70 harmonics.

The spatial profile of the electric field in the grating region illustrates the full benefits of the modified field calculation. Figure 4(a) depicts the conventional normal field component $| {E}_{{\rm{n}}}^{(21)}| =| {E}_{x}^{(21)}|$ in the gold grating for a moderately coarse discretization of N = 21 harmonics. It can be seen that $| {E}_{{\rm{n}}}^{(21)}|$ exhibits unphysical oscillations with a spatial frequency equal to the period of the smallest spatial frequency component in the FS expansion of the solution. This is the well known Gibbs phenomenon, which occurs when describing a discontinuous function with a truncated FS. On the other hand, the modified normal field, ${\tilde{E}}_{{\rm{n}}},$ does not suffer from such spurious oscillations at the interface. In particular, even for a small number of harmonics, N = 21, ${\tilde{E}}_{{\rm{n}}}^{(21)}$ is smooth, as per figure 4(b). At very large number of harmonics, N = 640, the modified normal field is free of any numerical artifacts, as can be seen in figure 4(c).

Zoom In Zoom Out Reset image size

The improved formulation of RCWA exhibits another benefit, namely ${\tilde{E}}_{{\rm{n}}}$ is by construction discontinuous and exactly fulfills the corresponding boundary condition

$$\begin{eqnarray}&&{\varepsilon }^{(\mathrm{in})}{{\bf{E}}}^{(\mathrm{in})}({{\bf{x}}}_{s})\cdot {\bf{n}}={\varepsilon }^{(\mathrm{out})}{{\bf{E}}}^{(\mathrm{out})}({{\bf{x}}}_{s})\cdot {\bf{n}},\end{eqnarray} \tag{ 5 }$$

at surface points of the grating, ${{\bf{x}}}_{s},$ where ${{\bf{E}}}^{\mathrm{in}}$ and ⁽ⁱⁿ⁾ (${{\bf{E}}}^{\mathrm{out}}$ and ^(out)) denote the electric field and permittivity inside (outside) the grating, respectively, and ${\bf{n}}$ is the unit vector normal to the surface. In the conventional formulation of the RCWA, the field E_n^N does not satisfy equation (5), because ${E}_{{\rm{n}}}^{(n)}$ is—as a FS containing a finite number of terms—inherently continuous. The closeup of the interfacial field around x = 0.25Λ in figure 5 emphasizes these ideas. Thus, the normal component of the field calculated using the conventional formulation of RCWA shows spurious oscillations for both small N = 21 (dashed green) and high N = 640 (solid red) number of harmonics, whereas the modified formulation is free of such unphysical oscillations and discontinuous for any number of harmonics, N (see the dotted purple and dashed–dotted blue lines corresponding to N = 21 and N = 640, respectively).

Zoom In Zoom Out Reset image size

It should be clear now that the modified normal-field evaluation in 1D by means of the displacement field represents an improvement of the conventional evaluation of E_n in two ways: first, it exhibits optimal self-convergence in the sense that it is as accurate as the far-field and the continuous tangential field component. Second, it explicitly fulfills the boundary condition (5) at the material interfaces. Equally important, the improved near-field calculation is achieved with minimal additional computational cost and does not alter the mathematical framework of the core RCWA algorithm in 1D.

The ideas of the previous section are straightforward to implement because in 1D-periodic structures there is a trivial and unambiguous distinction between the tangential and normal components of ${\bf{E}}.$ However, this situation becomes more intricate in the case of 2D-periodic diffraction gratings.

As it will be shown now, accurate near-field evaluation is strongly related to correct Fourier factorization in 2D-periodic structures. Fourier factorization means in this context the decomposition of a product of periodic functions, f = g · h, into its periodic factors g and h. Depending on the continuity properties of the factors and the product, different rules must be applied in order to obtain fast convergence when increasing the number of FS terms According to these rules, if f and h possess simultaneous discontinuities, but g is continuous at those locations, the product rule yields fast convergence with respect to the number of FS terms, namely one uses f = g · h. Moreover, if g and h are simultaneously discontinuous, but f is continuous, the inverse rule should be used, i.e., f must be factorized as f = (1/g)⁻¹ · h. A rigorous explanation of these rules and the solution to the 1D Fourier factorization problem is given in [18].

Because of the trivial distinction between the continuous (tangential) and discontinuous (normal) components of the electric field in 1D, Fourier factorization is straightforward in that case. But for 2D-periodicity, three different approaches to achieve the correct Fourier factorization have been proposed: (i) approximate the material boundaries by a coordinates-aligned staircase-contour [19]. (ii) Devise a coordinate system in which a given grating geometry is coordinate system aligned and use approach (i) to obtain the correct Fourier factorization [20, 38]. (iii) Construct a normal vector field (NVF) to decompose ${\bf{E}}$ into its normal and tangential components and then apply the corresponding correct factorization rules to them [21].

Since the accurate near-field evaluation relies on the decomposition of the electric field into normal and tangential components, it is natural to use the factorization approach (iii), the NVF approach. It is out of the scope of this paper to derive the full formulation of 2D-RCWA, so that only the crucial and unconventional steps will be given here. Thus, the normal and tangential components of ${\bf{D}}=\varepsilon {\bf{E}}$ have to be decomposed using the product- and inverse rule, respectively. This leads to the following relations for the α-component of the displacement field [21, 24, 39]:

$$\begin{eqnarray}&&[{D}_{\alpha }]={\varepsilon }_{0}\displaystyle \sum _{\beta =1}^{3}\left({\delta }_{\alpha ,\beta }[\kern-2pt[ \varepsilon ]\kern-2pt] -{{\rm{\Delta }}}_{\alpha \beta }\right)[{E}_{\beta }],\end{eqnarray} \tag{ 6 }$$

where δ_αβ is the Kronecker delta. Here, $[f]$ denotes the vector of FS coefficients of a scalar function f, $[\kern-2pt[ g]\kern-2pt]$ is the Toeplitz matrix of FS coefficients of g, and the matrix Δ_αβ is given by ${{\rm{\Delta }}}_{\alpha \beta }=\frac{1}{2}\left({\rm{\Delta }}[\kern-2pt[ {N}_{\alpha }{N}_{\beta }]\kern-2pt] +[\kern-2pt[ {N}_{\alpha }{N}_{\beta }]\kern-2pt] {\rm{\Delta }}\right)$ with ${\rm{\Delta }}=[\kern-2pt[ \varepsilon ]\kern-2pt] -[\kern-2pt[ 1/\varepsilon {]\kern-2pt] }^{-1}$ and N_α is the α-component of the NVF, ${\bf{N}}={({N}_{1},{N}_{2},{N}_{3})}^{T},$ of the material boundary. The matrix ${\delta }_{\alpha ,\beta }[\kern-2pt[ \varepsilon ]\kern-2pt] -{{\rm{\Delta }}}_{\alpha \beta }$ implements the three steps of the Fourier factorization: the decomposition of ${\bf{E}}$ into normal and tangential components (N_β in $[\kern-2pt[ {N}_{\alpha }{N}_{\beta }]\kern-2pt]$), factorization using either the inverse rule ($[\kern-2pt[ 1/\varepsilon {]\kern-2pt] }^{-1}$) or the product rule ($[\kern-2pt[ \varepsilon ]\kern-2pt]$ ), and back-projection to Cartesian coordinates (N_α in $[\kern-2pt[ {N}_{\alpha }{N}_{\beta }]\kern-2pt]$ ).

Note that one can also use asymptotically equivalent definitions of Δ_αβ. However, our investigations have shown that despite the fact that the choices ${{\rm{\Delta }}}_{\alpha \beta }={\rm{\Delta }}[\kern-2pt[ {N}_{\alpha }{N}_{\beta }]\kern-2pt]$ and ${{\rm{\Delta }}}_{\alpha \beta }=[\kern-2pt[ {N}_{\alpha }{N}_{\beta }]\kern-2pt] {\rm{\Delta }}$ produce similar convergence speed, they do not conserve the power for lossless structures, whereas the choice ${{\rm{\Delta }}}_{\alpha \beta }=[\kern-2pt[ {N}_{\alpha }]\kern-2pt] {\rm{\Delta }}[\kern-2pt[ {N}_{\beta }]\kern-2pt]$ yields power conservation but at the price of slower convergence. All three formulations are asymptotically equivalent with respect to the number of terms in the FS expansion, due to the commutation of Toeplitz operators [39].

Normal vector fields can be constructed analytically for a variety of structures and automated algorithms to obtain a NVF for arbitrary grating geometries have been developed [40]. It should be noted that this formulation allows inclined NVFs, i.e. NVFs with simultaneously non-vanishing x, y, and z components. Such a NVF is required to accurately model obliquely etched structures (see figure 1(d)) and hence can be viewed as a generalization of the methods presented in [24] (N_y = 0, figure 1(b)) and [21] (N_z = 0, figure 1(c)).

Given the constitutive relation (6), one can derive the RCWA eigenvalue-problem

$$\begin{eqnarray}&&{\mathsf{M}}\left(\begin{array}{c}[{S}_{x}]\\ [{S}_{y}]\\ [{U}_{x}]\\ [{U}_{y}]\end{array}\right)=\beta \left(\begin{array}{c}[{S}_{x}]\\ [{S}_{y}]\\ [{U}_{x}]\\ [{U}_{y}]\end{array}\right)\end{eqnarray} \tag{ 7 }$$

for the electromagnetic modes described by $[S]$ and $[U]$ and propagation constant, β. The most general formulation of the system-matrix ${\mathsf{M}}$ reads:

$$\begin{eqnarray*}{\mathsf{M}} & = & \left(\begin{array}{cc}{{\mathsf{K}}}_{x}{\mathsf{B}}{{\rm{\Delta }}}_{{zx}} & {{\mathsf{K}}}_{x}{\mathsf{B}}{{\rm{\Delta }}}_{{zy}}\\ {{\mathsf{K}}}_{y}{\mathsf{B}}{{\rm{\Delta }}}_{{zx}} & {{\mathsf{K}}}_{y}{\mathsf{B}}{{\rm{\Delta }}}_{{zy}}\\ {{\rm{\Delta }}}_{{yx}}+{{\rm{\Delta }}}_{{yz}}{\mathsf{B}}{{\rm{\Delta }}}_{{zx}}-{{\mathsf{K}}}_{x}{{\mathsf{K}}}_{y} & {{\mathsf{K}}}_{x}{{\mathsf{K}}}_{x}+{{\rm{\Delta }}}_{{yz}}{\mathsf{B}}{{\rm{\Delta }}}_{{zy}}-{{\mathsf{C}}}_{y}\\ {{\mathsf{C}}}_{x}-{{\rm{\Delta }}}_{{xz}}{\mathsf{B}}{{\rm{\Delta }}}_{{zx}}-{{\mathsf{K}}}_{y}{{\mathsf{K}}}_{y} & {{\mathsf{K}}}_{y}{{\mathsf{K}}}_{x}-{{\rm{\Delta }}}_{{xz}}{\mathsf{B}}{{\rm{\Delta }}}_{{zy}}-{{\rm{\Delta }}}_{{xy}}\end{array}\right.\\ & & \qquad \left.\begin{array}{cc} & {{\mathsf{K}}}_{x}{{\mathsf{BK}}}_{y}{\mathsf{I}}-{{\mathsf{K}}}_{x}{{\mathsf{BK}}}_{x}\\ & {{\mathsf{K}}}_{y}{{\mathsf{BK}}}_{y}-{\mathsf{I}}-{{\mathsf{K}}}_{y}{{\mathsf{K}}}_{x}\\ & {{\rm{\Delta }}}_{{yz}}{{\mathsf{BK}}}_{y}-{{\rm{\Delta }}}_{{yz}}B{{\mathsf{K}}}_{x}\\ & -{{\rm{\Delta }}}_{{xz}}{{\mathsf{BK}}}_{y}{{\rm{\Delta }}}_{{xz}}{{\mathsf{BK}}}_{x}\end{array}\right),\end{eqnarray*}$$

with ${\mathsf{B}}={\left([\kern-2pt[ \varepsilon ]\kern-2pt] -{{\rm{\Delta }}}_{{zz}}\right)}^{-1}$ and ${{\mathsf{C}}}_{\alpha }=[\kern-2pt[ \varepsilon ]\kern-2pt] -{{\rm{\Delta }}}_{\alpha \alpha }.$ The matrices ${{\mathsf{K}}}_{\alpha }={\rm{diag}}\left({k}_{\alpha n}\right)$ are diagonal matrices of the in-plane propagation constants, k_xn for α = x and k_yn for α = y, of the diffraction orders and ${\mathsf{I}}$ is the identity matrix. The size of matrices ${\mathsf{B}},$ ${{\mathsf{C}}}_{\alpha },$ ${{\mathsf{K}}}_{\alpha },$ and ${\mathsf{I}}$ is N₀ × N₀, whereas the size of ${\mathsf{M}}$ is 4N₀ × 4N₀, where N₀ = (2N + 1)² is the total number of FS coefficients in the calculation and the same number of FS terms for the truncation of the FS in the x- and y-direction is assumed.

The mode amplitudes $[S]$ and $[U]$ are determined by matching the tangential components of the electromagnetic field at the top and bottom of the grating region, whereas multilayered structures are described by the staircase approximation in conjunction with the numerically stable ${\mathcal{S}}$-matrix algorithm. This yields the FS coefficients, $[{\bf{E}}],$ of the electric field everywhere in the grating region, in the cover, and the substrate.

Let us now denote by ${\mathcal{R}}\left([f]\right)$ the FS reconstruction of a Fourier coefficient vector $[f],$

$$\begin{eqnarray*}&&{\mathcal{R}}\left([f]\right)(x,y)=\displaystyle \sum _{n=-N}^{N}\displaystyle \sum _{m=-N}^{N}{f}_{{nm}}\mathrm{exp}\left({\rm{i}}m\displaystyle \frac{2\pi }{{{\rm{\Lambda }}}_{1}}x+{\rm{i}}n\displaystyle \frac{2\pi }{{{\rm{\Lambda }}}_{2}}y\right),\end{eqnarray*}$$

where ${f}_{{nm}}={[f]}_{(n+N)(2N+1)+m+N+1}$ are the N₀ FS coefficients of f. Within the conventional RCWA framework, each component of ${\bf{E}}$ is evaluated as:

$$\begin{eqnarray}&&{E}_{\alpha }^{(N)}={\mathcal{R}}\left([{E}_{\alpha }]\right).\end{eqnarray} \tag{ 8 }$$

This relation does not take into account the continuity properties of the different components of ${\bf{E}}$ and hence will lead to spurious oscillations and slow convergence of the near-field as it was seen in section 2.

The modified field evaluation for 2D-periodic diffraction gratings requires one to define the continuous normal component of ${\bf{D}}$ and the tangential component of ${\bf{E}}.$ Their FS coefficient vectors are given by:

$$\begin{eqnarray}&&[{{\bf{D}}}_{{\rm{n}}}]=\displaystyle \frac{1}{2}{\varepsilon }_{0}\left([\kern-2pt[ {{\bf{NN}}}^{T}]\kern-2pt] [\kern-2pt[ 1/\varepsilon {]\kern-2pt] }^{-1}+[\kern-2pt[ 1/\varepsilon {]\kern-2pt] }^{-1}[\kern-2pt[ {{\bf{NN}}}^{T}]\kern-2pt] \right)[{\bf{E}}],\end{eqnarray} \tag{ 9a }$$

$$\begin{eqnarray}&&[{{\bf{E}}}_{{\rm{t}}}]=[\kern-2pt[ {\mathbb{1}}-{{\bf{NN}}}^{T}]\kern-2pt] [{\bf{E}}],\end{eqnarray} \tag{ 9b }$$

where ${\mathbb{1}}$ denotes the 3 × 3 identity matrix and ${{\bf{NN}}}^{T}$ is the 3 × 3 projection matrix defined by the NVF at any point in space, ${\bf{x}}.$ Note that by construction, $[{{\bf{D}}}_{{\rm{n}}}]$ and $[{{\bf{E}}}_{{\rm{t}}}]$ are FS coefficient vectors of vector fields that are continuous at material interfaces. Hence their reconstructions, ${\mathcal{R}}\left([{{\bf{D}}}_{{\rm{n}}}]\right)$ and ${\mathcal{R}}\left([{{\bf{E}}}_{{\rm{t}}}]\right),$ do not suffer from the Gibbs phenomenon at the interface. With this observation in mind, the electric field in the improved RCWA method for 2D-periodic structures at a point, ${\bf{x}},$ is given by:

$$\begin{eqnarray}&&{\tilde{{\bf{E}}}}^{(N)}({\bf{x}})={\varepsilon }_{0}^{-1}{\mathcal{R}}\left([{{\bf{D}}}_{{\rm{n}}}]\right)({\bf{x}})/\varepsilon ({\bf{x}})+{\mathcal{R}}\left([{{\bf{E}}}_{{\rm{t}}}]\right)({\bf{x}})\end{eqnarray} \tag{ 10 }$$

and is expected to yield fast near-field convergence and non-oscillatory spatial field profiles. To investigate the validity of these predictions, the accurate field evaluation was implemented in a commercially available RCWA computer software, OmniSim/RCWA [41].

It should be noted that the other electromagnetic fields can be easily calculated with our improved method, too. Specifically, the displacement field, ${\bf{D}},$ can be evaluated using the modified electric field $\tilde{E},$ namely ${{\bf{D}}}^{(N)}({\bf{x}})={\varepsilon }_{0}\varepsilon ({\bf{x}}){\tilde{{\bf{E}}}}^{(N)}({\bf{x}}),$ and hence will have the same convergence properties as $\tilde{{\bf{E}}}.$ The magnetic induction ${\bf{B}}$ and the magnetic field ${\bf{H}}=\mu {\bf{B}}$ do not require special attention, because they are continuous in non-magnetic materials and hence behave similar to the continuous tangential component of the electric field.

In this section, the improved formulation for accurate near-field calculations in 2D-periodic structures is assessed using two test structures under different configurations. To this end, we first extend the definition of the grating norm (2) of a scalar or vector function, f, to 2D-periodic structures in the following straightforward way:

$$\begin{eqnarray}&&{\parallel f\parallel }_{G}={\left({\displaystyle \int }_{0}^{h}{\displaystyle \int }_{-{{\rm{\Lambda }}}_{2}/2}^{{{\rm{\Lambda }}}_{2}/2}{\displaystyle \int }_{-{{\rm{\Lambda }}}_{1}/2}^{{{\rm{\Lambda }}}_{1}/2}{| f(x,y,z)| }^{2}\ {\rm{d}}x\ {\rm{d}}y\ {\rm{d}}z\right)}^{\frac{1}{2}},\end{eqnarray} \tag{ 11 }$$

where the integral is evaluated over the three-dimensional grating region.

4.1. Analysis of a 1D-periodic grating using 2D-RCWA

The first 2D-periodic diffraction grating under consideration is a 1D binary grating, as shown in figure 1(a), rotated by π/4 in the x–y-plane. It should be obvious that it can be modeled as a double-periodic 2D grating with periods Λ₁ = Λ₂ =$\sqrt{2}\tilde{{\rm{\Lambda }}}=\sqrt{2}\mu {\rm{m}},$ where $\tilde{{\rm{\Lambda }}}={{\rm{\Lambda }}}_{1}$ is the period of the grating when it is viewed as a 1D-periodic structure. For clarity, the primary unit cell of the grating is depicted in the inset of figure 6(a). With this choice, the reference quantities R^ref, T^ref, and ${{\bf{E}}}^{\mathrm{ref}}$ can be calculated using 1D simulations, an approach inspired by an example in [21].

Zoom In Zoom Out Reset image size

The error in the calculation of the far-field, ${\text{}}{e}_{{\rm{F}}}$ from equation (1), when 2D simulations are employed is depicted in figure 6(a). The convergence of the calculations for the silicon and gold gratings closely follows the convergence trends observed when 1D simulations are performed and, as expected, it shows somewhat worse, yet still good, agreement for the silica structure. This difference is explained by the fact that the NVF introduced in the 2D-RCWA formulation is discontinuous away from material interfaces and hence it can degrade the convergence rate, especially for the low-index of refraction contrast case.

The conclusions of our analysis of the convergence of the near-field are summarized in figures 6(b)–(d). Thus, for both the silica and gold diffraction gratings, the normal component of the modified electric field, $\tilde{{\bf{E}}},$ exhibits faster convergence than this same field component calculated using the conventional form of RCWA. In both cases, ${\text{}}e\{{\tilde{E}}_{{\rm{n}}}^{(20)}\}$ is one order of magnitude smaller than ${\text{}}e\{{E}_{{\rm{n}}}^{(20)}\}.$ For the silicon grating, only marginal differences between the two formulations can be observed. This is in agreement with the results obtained in the 1D case, as per figure 3(b). For small N, ${\tilde{{\bf{E}}}}^{(N)}$ and ${{\bf{E}}}^{(N)}$ are determined with a comparable degree of accuracy, but for the larger number of harmonics considered in figure 6(c), i.e. N = 30, a higher accuracy of our improved formulation of the RCWA can clearly be observed.

These conclusions are further validated by the profile of the electric field, as presented in figure 7. This figure shows the spatial distribution of the normal component of the electric fields, ${{\bf{E}}}^{(27)}$ and ${\tilde{{\bf{E}}}}^{(27)},$ calculated in the median plane of the grating. A simple examination of these field profiles confirms that the spurious oscillations of the field ${\tilde{{\bf{E}}}}^{(27)}$ near the surface have much smaller amplitude as compared to that of the variations of ${{\bf{E}}}^{(27)}.$ Moreover, a closer inspection of the surface-fields shows that the boundary condition (5) is fulfilled by ${\tilde{{\bf{E}}}}^{(27)}$ only. This first test-case already reveals that the improved near-field evaluation is more accurate in the case of 2D-periodic structures, too.

Zoom In Zoom Out Reset image size

4.2. Near-field calculations for an intrinsically 2D-periodic grating

In order to thoroughly test the near-field evaluation for 2D-periodic structures using the improved RCWA presented in this article, in what follows we consider the challenging test structure depicted in figure 1(c). The grating region consists of a coordinate system aligned parallelepiped with the length of the sides aligned to the x-, y-, and z-axis being a = 0.5Λ₁, 2a, and h, respectively, placed adjacently to a semicircular cylinder with radius a and height h, with Λ = Λ₁ = Λ₂ = 0.25 μm. The structure is illuminated normally by a x-polarized plane wave with wavelength λ = 0.5 μm.

Since it is generally computationally time consuming to obtain high-accuracy solutions in the case of 2D-periodic structures and due to the fact that the higher the ratio Λ/λ and the refractive index $| n| ,$ the more harmonics are necessary to achieve convergence [35], a relatively small period-to-wavelength ratio of Λ/λ = 0.25 μm/0.5 μm = 0.5 is chosen for this example.

As reference values in the definition of the far-field error, ${\text{}}{e}_{{\rm{F}}}$ from equation (1), T^ref = T³¹ and R^ref = R³¹, namely results obtained from simulations with N = 31 are chosen. The convergence of the far-field, ${\text{}}{e}_{{\rm{F}}}(N),$ is shown in figure 8(a), where the far-field physical quantity considered are the transmission and reflection coefficients, T and R, respectively. As expected, the fastest convergence can be observed for the silica grating, because it has a low refractive index. The numerical error obtained for the gratings made of gold and silicon are one and two orders of magnitude larger than in the case of the silica grating, respectively. This behavior is similar to that seen in the 1D case for a small number of harmonics, N < 30 (see figure 2).

Zoom In Zoom Out Reset image size

Figures 8(b)–(d) contain the dependence of the near-field self-error determined for the three material configurations, silica, silicon, and gold, respectively. The three physical quantities plotted in each case are the z-component of ${\bf{E}}$ and the in-plane components, ${{\bf{E}}}_{{xy}}:= {({E}_{x},{E}_{y},0)}^{T}$ and ${\tilde{{\bf{E}}}}_{{xy}}:= {({\tilde{E}}_{x},{\tilde{E}}_{y},0)}^{T},$ obtained by using the conventional and improved versions of RCWA, respectively. Note that the in-plane component contains the discontinuous normal component, ${\bf{N}}\cdot {\bf{E}}$ (${\bf{N}}\cdot \tilde{{\bf{E}}}$), and the continuous tangential component, $({\mathbb{1}}-{{\bf{NN}}}^{T}){{\bf{E}}}_{{xy}}$ ($({\mathbb{1}}-{{\bf{NN}}}^{T}){\tilde{{\bf{E}}}}_{{xy}}$).

It can be seen that in all cases, the z-component of ${\bf{E}},$ which is continuous at vertical surfaces inside the grating region, converges much faster than the in-plane component, in both the conventional and modified formulations. In addition, the modified formulation leads to a somewhat smaller error than the conventional formulation. Specifically, it was found that ${\text{}}e\{{\tilde{{\bf{E}}}}_{{xy}}^{(N)}\}\approx 0.5{\text{}}e\{{{\bf{E}}}_{{xy}}^{(N)}\}$ for the silica grating and ${\text{}}e\{{\tilde{{\bf{E}}}}_{{xy}}^{(N)}\}\approx 0.9{\text{}}e\{{{\bf{E}}}_{{xy}}^{(N)}\}$ for the silicon and gold gratings. However, the corresponding convergence speed, i.e. the slope of ${\text{}}e\{{{\bf{E}}}_{{xy}}^{(N)}\}$ and ${\text{}}e\{{\tilde{{\bf{E}}}}_{{xy}}^{(N)}\},$ is the same. Moreover, in the modified formulation of the RCWA the convergence speed of the tangential component of the near-field is larger than that of the in-plane component. Three factors contribute to this behavior: (i) the decomposition of the near-field in a normal and tangential component in 2D-periodic structures relies on the specific definition of the NVF, ${\bf{N}}$ and it is not directly performed in Cartesian coordinates as in the 1D case. Hence, the inexact field decomposition by ${\bf{N}}$ introduces an additional error. (ii) The field of normal vectors characterizing the structure is only uniquely defined at the interfaces defining the grating, except at the corners, and, more importantly, away from the grating surface. This ambiguity can lead to choices of NVFs which are not optimal for the convergence and accurate calculation of the near-field. (iii) The NVF itself has discontinuities, which can cause additional oscillations in the spatial profile of the electromagnetic field.

It is also worthwhile to investigate the spatial profile of the near-field. The dominant x-component of the electric field in a horizontal cross-section through the grating region at z = h/2 is depicted in figure 9. As in the 1D case, these maps show that the field ${\tilde{E}}_{x}$ exhibits spatial oscillations with smaller amplitude as compared to the variations of the field E_x, especially at y-aligned interfaces (outlined with blue dashed lines in figure 9).

Zoom In Zoom Out Reset image size

Oblique diffraction gratings as those shown in figures 1(b) and (d) are modeled in the RCWA within the staircase approximation along a direction perpendicular onto the grating plane. Each computational slice is assumed to be a z-independent structure, for which the eigenmodes can be found numerically by solving equation (7). The field in each computational layer is then found by using the boundary conditions at the top and bottom interfaces of the grating, employing a numerically stable ${\mathcal{S}}$-matrix formulation. The validity of this staircase approximation for 1D-periodic gratings under TE-polarization has been proven in [35]. In the context of the NVF formulation of RCWA for oblique 1D-periodic structures, it has been found that the use of an out-of-plane NVF, i.e. ${\bf{N}}={\left({N}_{x},0,{N}_{z}\right)}^{T},$ is beneficial for the far-field convergence speed [24, 35]. In this section we will study the relation between the accurate field formulation (10) and the near-field convergence speed for slanted 1D-periodic structures.

The situation for oblique 2D-periodic structures has not yet been explored in the context of RCWA. However, based on the results presented so far, one can conjecture that both the near- and far-field convergence of RCWA can be improved by using an out-of-plane 3D NVF, ${\bf{N}}={\left({N}_{x},{N}_{y},{N}_{z}\right)}^{T}\ne 0,$ and that the near-field evaluation would be more accurate as well. The validity of this supposition will be explored in the second part of this section.

5.1. Analysis of slanted 1D-periodic binary diffraction gratings

In order to analyze oblique 1D-periodic structures, we consider the slanted binary grating depicted in figure 1(b). The period of the grating is Λ = 1 μm, the filling factor, ρ = 0.5, the height, h = 0.25 μm, and the slanting angle is θ = π/4. Only the gold grating is considered in this section as this would be the most challenging case. If the unit cell is assumed to extend from x = −Λ/2 to x = Λ/2 and the center of the binary grating is set to be x = 0, z = h/2, a suitable out-of-plane NVF is given by $\tilde{{\bf{N}}}(x,y,z)={\rm{sign}}\left(x-z\right)\sqrt{2}{\left(\mathrm{1,0,1}\right)}^{T}.$

The two formulations of the RCWA compared in this section are the in-plane NVF, which is used in conventional RCWA together with the conventional field evaluation (8), and the out-of-plane NVF, $\tilde{{\bf{N}}},$ combined with the improved field evaluation formulation (10). In contrast to the results presented in the previous sections, the two formulations yield different results for both the near- and far-field quantities. Moreover, for the sake of the clarity of the presentation, all physical quantities corresponding to the out-of-plane NVF formulation are denoted with a tilde symbol.

Numerical results for increasing number of harmonics, N = 2, ..., 320, and number of computational layers, M = 2, ..., 256, are presented in figure 10. It can be inferred from this figure that the in-plane formulation requires both a high number of FS coefficients, N, and layers, M, to achieve convergence to a result of R^ref = 0.287 88, whereas the out-of-plane formulation yields fast convergence to ${\tilde{R}}^{\mathrm{ref}}={R}^{(\mathrm{320,256})}=0.288\;37$ with respect to N, as per figures 10(a) and (b), respectively. This behavior is in agreement with the findings reported in [24]. The convergence of the calculated near-field, illustrated in figures 10(c) and (d), exhibits similar features. Specifically, the in-plane NVF formulation requires both high N and M to achieve a small self-error of ${\text{}}e\{{{\bf{E}}}^{(\mathrm{226,256})}\}=4.7\times {10}^{-2},$ whereas this self-error can already be achieved with N = 10, M = 256 in the out-of-plane formulation. This clearly demonstrates a drastically improved efficiency to the calculation of the near-field of oblique diffraction gratings of the approach based on the combination of out-of-plane NVF and the accurate near-field formulation. The highly improved near-field profile is illustrated in figure 11(b), which exhibits no unphysical oscillations near the gold–vacuum interface. This is in sharp contrast to the conventional field evaluation of the in-plane formulation (see figure 10(a)), which clearly suffers from spurious oscillations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.2. Analysis of slanted 2D-periodic cylindrical diffraction gratings

In this section we investigate the efficiency of using the accurate field evaluation and the out-of-plane NVF formulation to model a challenging, slanted 2D-periodic diffraction grating. The grating with periods, Λ₁ = Λ₂ = Λ = 1 μm, is schematically depicted in figure 1(d) and consists of a cylindrical rod with radius r = 0.3 μm and height h = 0.125 μm, which is slanted by θ = π/4 along the x-axis. Again, only the gold grating is considered in this section. Finally, the incident plane wave is impinging onto the grating along the normal direction, is polarized along the x-axis, and has a wavelength of λ = 2 μm.

The reflection coefficient calculated for N = 3, ..., 19 harmonics and M = 2, ..., 32 layers is shown in figures 12(a) and (b) and was determined by using the conventional in-plane NVF and the out-of-plane NVF formulation, respectively. It can be seen that both approaches converge rather slowly, neither one achieves convergence even for the highest considered values of M = 32 and N = 19. Moreover, in order to characterize the error of the near-field calculations in the two formulations, the self-error with respect to the reference solutions obtained with N = 19 and M = 32 is presented in figures 12(c) and (d). The in-plane formulation achieves a relative self-error of ${\text{}}e\{{{\bf{E}}}^{(\mathrm{13,32})}\}=0.625,$ whereas for the same values of N and M, the modified field evaluation in conjunction with the out-of-plane NVF achieves a substantially lower self-error of ${\text{}}e\{{{\bf{E}}}^{(\mathrm{13,32})}\}=0.261.$ It has to be stressed, that the necessary accuracy for full convergence could not be achieved in our simulations. The evolution of the computational results for the shown values of $N\leqslant 19$ and $M\leqslant 32$ however can be interpreted in favor of the out-of-plane NVF formulation, due to the lower near-field self-error ${\text{}}e\{{{\bf{E}}}^{(\mathrm{13,32})}\}\lt {\text{}}e\{{{\bf{E}}}^{(\mathrm{13,32})}\}.$ It can be supposed that future simulations with finer discretization will reveal the practical benefit of the out-of-plane NVF formulation in conjunction with the modified field formulation for 2D-periodic slanted structures, similar to the case of 1D-periodic slanted diffraction gratings.

Zoom In Zoom Out Reset image size

To summarize, we have analyzed the numerical near-field calculated using the RCWA and identified the Gibbs phenomenon as the main reason for their slow convergence and the spurious oscillations from which they suffer. As a solution to these deficiencies of RCWA, we proposed an improvement of this method that can be applied for modeling arbitrary diffraction gratings and, more generally, periodic optical structures. The modified formulation significantly improves the accuracy of 1D-RCWA calculations, for both straight and slanted gratings, where it speeds up convergence and removes the numerical artifacts from the calculated near-fields. The accuracy of 2D-periodic grating simulations can be enhanced, however, to a lesser extent than in the 1D case. The reduced performance in 2D can be attributed to the discontinuity and non-exactness of the numerical NVF, which is at the core of the modified formulation. Therefore, it might be fruitful to investigate more elaborate NVF-formulations and their suitability for near-field calculations, such as a complex valued NVF [42], which unlike the NVF used here is continuous everywhere in the grating region.

We expect that the proposed modification of the RCWA method will greatly advance its computational capabilities, especially for 1D periodic optical structures. In particular, this improved method could prove instrumental to accurate modeling of periodic plasmonic structures, diffraction gratings, and surface-nonlinear devices, namely to simulation of physical systems whose functionality rely on the electromagnetic near-field at interfaces.

The work of M Weismann was supported through a UCL Impact Award graduate studentship funded by UCL and Photon Design Ltd. N C P acknowledges support from European Research Council/ERC Grant Agreement no. ERC-2014-CoG-648328.