1. Introduction

Methods of optical modeling of periodic systems (or large systems with periodic boundary conditions) based on the Fourier expansion of spatially arranged material parameters and electromagnetic fields have undergone substantial development during past decades. According to numerical implementation or parameters of interest, these methods are known under various names such as the Fourier modal method (FMM), coupled wave method, plane wave expansion method, etc. In 1990s the Fourier-based approaches received renewed interest owing to high increase in amount and quality of microstructure manufacturing and owing to new concepts such as photonic crystals and devices. Since then, considerable improvements in the numerical efficiency of those methods took place such as stable propagation algorithms or proper Fourier factorization (FF).

One of the crucial problems of poor numerical efficiency was found when the electric field was partially perpendicular to metallic surfaces, i.e., to large permittivity discontinuities [1]. After that several authors suggested improvements by careful choices whether permittivity or impermittivity should be Fourier-factorized via the Laurent convolution rule at particular places of Maxwell’s equations [2–5]. These principles have then been analyzed from the mathematical viewpoint and summarized as the FF rules [6, 7]. Finally these rules have been successfully applied by many authors to more general structures, including an arbitrary shape, dimensionality, anisotropy, or (a)periodicity [8]. However, the FF rules have, to our knowledge, been always treated with respect to permittivity discontinuities.

In this article we generalize this concept to structures with general inhomogeneity of permittivity in the lateral direction, including systems with continuous variations of permittivity, where proper FF rules should also be applied, which we demonstrate on several examples of FMM simulations of diffraction efficiencies for the transverse magnetic (TM) polarization, for which the FF rules are of particular importance. The improved efficiency of the optical modeling of continuous permittivity structures will be useful in various modern complicated systems such as gradient refractive index structures created by holographic lithography [9], photo-induced refractive index changes in chalcogenide glasses [10], nanolayer polymer extrusion [11], or by employing external tuning based on electric, magnetic, acoustic, thermal, or other excitation [12].

2. Fourier modal method for a diffraction grating

The FMM is based on solving modal propagation equations in a periodic medium and boundary conditions on interfaces. From the numerical viewpoint the spatial dependences of fields and permittivity are Fourier-expanded for subsequent manipulation by linear-algebraic methods. Normalized Maxwell’s equations in a periodic medium for the TM polarization can be written

 

Fig. 1 Geometric configuration with planar diffraction mounting.

Download Full Size | PDF

In the frame of the FMM we can expand into the Fourier series either the relative permittivity function, $\epsilon \left(y\right)={\displaystyle {\sum }_{n=-\infty }^{+\infty }{\epsilon }_n\exp \left(iqny\right)}$, or the impermittivity function, $\eta \left(y\right)=1/\epsilon \left(y\right)={\displaystyle {\sum }_{n=-\infty }^{+\infty }{\eta }_n\exp \left(iqny\right)}$, where q = λ/Λ, with Λ denoting the grating’s period in the SI units, and—according to Floquet’s theorem—all the field components into the pseudo-Fourier series, $H_x\left(y,z\right)={\displaystyle {\sum }_{n=-\infty }^{+\infty }{h}_{xn}\exp \left(i{q}_{n}y\right)}$, $E_j\left(y,z\right)={\displaystyle {\sum }_{n=-\infty }^{+\infty }{e}_{jn}\exp \left(i{q}_{n}y\right)}$, with q_n = q₀ + nq, where q₀ = sin ϑ_i is determined by the angle of incidence ϑ_i, assuming vacuum in the superstrate.

The constitutive relations (formulas for the scaled components of the displacement field) D_j = ε(y)E_j present in Eqs. (2) and (3) can then be treated either by Laurent’s rule, $d_{jn}={\displaystyle {\sum }_{k=-\infty }^{+\infty }{\epsilon }_{n-k}{e}_{jk}}$, where d_jn are the pseudo-Fourier coefficients of the displacement field (analogous to e_jn) or the inverse rule, $e_{jn}={\displaystyle {\sum }_{k=-\infty }^{+\infty }{\eta }_{n-k}{d}_{jk}}$. Both can be rewritten into matrix forms, ${\mathbf{d}}_j〚\epsilon 〛{\mathbf{e}}_j$ or ${\mathbf{d}}_j={〚\eta 〛}^{-1}{\mathbf{e}}_j$, respectively, where ${〚\epsilon 〛}_{nk}={\epsilon }_{n-k}$ denotes an element of the Toeplitz matrix of the Fourier components of the function ε(y) (and analogously for $〚\eta 〛$) and where e_j and d_j are column vectors of the values e_jn and d_jn, respectively. Analogously we define the column vector h_x of the pseudo-Fourier coefficients h_xn. Furthermore, the operation of the ∂_y differentiation can be replaced by [∂_y] = iq, where q is a diagonal matrix, q_mn = q_nδ_mn. For the sake of computer implementation, all the above infinite summations are replaced by finite summations, ${\displaystyle \sum { }_{n=-n_{\max }}^{+n_{\max }}}$, so that the matrix forms of the constitutive relations have the dimension 2n_max + 1 (the order of the matrices 〚ε〛, 〚η〛, and q and the length of the column vectors e_j, d_j, and h_x).

The FMM algorithm consists of two steps. First, to calculate the propagation modes in the periodic medium we assume the exp(isz) dependence of each field component, where s is the mode propagation number. Then we can rewrite Eq. (4) into the form (1 denotes the unit matrix)

Second, to calculate the scattering problem, we need a relation between the field components E_y and H_x, which are tangential to the interfaces between the periodic and ambient media. For the TM polarization we assume the principal component H_x and the derived component E_y determined by Eq. (2), whose matrix version is

3. Numerical investigation of continuous profiles

According to the FF rules derived for discontinuous permittivity profiles (discontinuities assumed along the y axis), we have to put ε_I = 〚η〛⁻¹ and ε_II = 〚ε〛, because D_y, E_z are continuous and E_y, D_z discontinuous [7]. Since we cannot simply generalize these rules to the case of continuous permittivity, we will perform numerical experiments with all the four combinations according to Table 1, where Model A corresponds to the correct FF rules derived for discontinuous profiles and Models B, C, and D denote the three remaining alternatives.

Table 1. Four models according to the used FF rules

View Table | View all tables in this article

Below we compare the numerical performance for one discontinuous (rectangular) and five continuous examples of permittivity profiles, a piecewise linear profile (with the discontinuous first derivative), a piecewise quadratic profile (with the discontinuous second derivative), an analytic profile (with all derivatives continuous and with infinite Fourier series), a profile with sinusoidal permittivity (with all derivatives continuous and with only three terms of the Fourier series), and a profile with sinusoidal impermittivity (with all derivatives continuous and with three terms of impermittivity’s Fourier series). The corresponding formulas are written in Table 2, and their graphical dependences are displayed in Fig. 2.

Table 2. Modeled examples of periodic permittivity profiles (in SI units)

View Table | View all tables in this article

 

Fig. 2 Simulated examples of permittivity profiles, one discontinuous (rectangular) and five continuous. The top row displays the permittivity (ε) profiles, the bottom row the impermittivity (η) profiles. The horizontal axis is scales to unit periodicity.

Download Full Size | PDF

In the permittivity definitions we set ε_min = 4, ε_max = 8, and σ = 5. We also assume permittivity of the substrate ε_sub = 2.25, the depth of the grating medium d = 300 nm, the grating period Λ = 700 nm, the wavelength λ = 700 nm, and the angle of incidence ϑ_i = 70°.

The Fourier coefficients of permittivity (ε_n) and impermittivity (η_n) are calculated numerically. Each function is discretized into N = 2²⁴ equidistant values ε(y_k), η(y_k) with y_k = (k − 0:5)Λ/N (k = 1, 2, …; N), which are used in the fast Fourier transform formula, ${\epsilon }_n=\left(1/N\right){\displaystyle {\sum }_{k=1}^N\epsilon (y_k)\exp \left[-2\pi i\left(n-1\right)\left(k-1\right)/N\right]}$, of which we only use the values with n ∈ {−n_max, …, n_max}, and analogously for η_n.

We estimate numerical precision values as the differences $\mathrm{\Delta }{\rho }_{n_{\max }}={\rho }_{n_{\max }}-{\rho }_{n_{\max }^{{ }^{\prime }}}$, where ${\rho }_{n_{\max }}$ is the efficiency of a zeroth-order diffracted wave for which we retain the Fourier harmonic components up to the order n_max. The value of ${\rho }_{n_{\max }^{{ }^{\prime }}}$ is the last calculated value in the convergence calculation $(n_{\max }=1,2,\dots ,n_{\max }^{{ }^{\prime }})$. We carry out calculations with $n_{\max }^{{ }^{\prime }}=150,100,70,60,14$, and 14, corresponding to the rectangular, piecewise linear, piecewise quadratic, analytic, ε-sinusoidal, and η-sinusoidal profiles, whose convergence properties are displayed in Figs. 3(a)–3(f), respectively, with shorter n_max displayed ranges.

 

Fig. 3 Convergence properties of the four models simulated for a discontinuous (a) and five continuous (b–f) permittivity profiles.

Download Full Size | PDF

Figure 3(a) displays the convergence corresponding to the discontinuous, piecewise constant permittivity profile, frequently discussed in previous papers. The blue curve (Model A) yields the fastest convergence, followed by the curves corresponding to Models D (only η expanded), B (only ε expanded), and C (FF opposite to Model A). Moreover, for high values of n_max the precision of Model A becomes two-orders higher (with smaller $\mathrm{\Delta }{\rho }_{n_{\max }}$) than that of the remaining models. Similar convergence properties are visible in Fig. 3(b), corresponding to piecewise linear (triangular) permittivity profile, which is continuous but whose first derivative is discontinuous. Here for high values of n_max Model A yields almost two-orders higher precision than Model D, which yields one-order higher precision than Models B and C.

On the other hand, the permittivity profile defined as periodic quadratic splines (with the first derivative continuous but the second discontinuous) does not exhibit so apparent convergence differences among the four models [Fig. 3(c)]. Most of the precision values corresponding to Model A are again highest, followed by Models D, B, and C. Similar properties are visible in Fig. 3(d), corresponding to an analytical permittivity profile (with all derivatives continuous). Deviations from the convergence trends are probably due to the range of retained harmonics from −n_max to n_max, which is not entirely suitable for oblique incidence.

Figures 3(e) and 3(f), corresponding to other analytic profiles (ε-sinusoidal and η-sinusoidal) exhibit convergence properties different than those above. Here both Models A and B yield similar convergence properties, followed by Models C and D, whose precision is nearly one-order lower. In the case of the ε-sinusoidal profile, Model B converges even faster for n_max ≤ 6, which can be explained by only three nonzero terms in the Fourier series of the expanded permittivity function.

In general, the importance of the correct FF rules decreases when the number of continuous derivatives increases. It can be explained by the rule of the Fourier theory that the coefficients of the Fourier series of a periodic function f converge as fast as the sequence 1/n^p+1, where p is the number of derivatives after which the function f ^(p) becomes discontinuous [13]. Analytic profiles (with p = ∞) thus yield so fast convergence for all the four models that the correct FF rules significantly loose their advantage. Nevertheless, is can be stated that the rules applied for Model A are correct for any inhomogeneous, not only discontinuous profile.

4. Conclusion

Analysis of the convergence properties of the diffraction efficiencies of several diffraction gratings demonstrated that optical modeling methods based on Fourier expansion (such as the FMM used here) require proper FF rules not only on structures with sharp permittivity discontinuities (such as patterned metallic nanostructures) but also on structures with continuous permittivity (such as gradient refractive index structures), although they do not yield significant advantages on structures with analytic profiles.