Investigation of the Convergence Behavior with Fluctuation Features in the Fourier Modal Analysis of a Metallic Grating
 Author: Kim Hwi, Park Gwanwoo, Kim Changsoon
 Organization: Kim Hwi; Park Gwanwoo; Kim Changsoon
 Publish: Journal of the Optical Society of Korea Volume 16, Issue3, p196~202, 25 Sep 2012

ABSTRACT
We observe that the transmission and reflection efficiencies of a onedimensional metallic grating under transversemagnetic illumination calculated using the Fourier modal method (FMM) with the Fourier factorization rules have peculiar fluctuations, albeit small in magnitude, as the number of field harmonics increases. It is shown that when the number of Fourier terms for the electromagnetic field is increased from that in the conventional FMM, the fluctuations due to nonconvergent highly evanescent eigenmodes can be eliminated. Our examination reveals that the fluctuations originate from the Gibbs phenomenon inherent in the Fourierseries representation of a permittivity function with discontinuities, and from nonconvergence of highly evanescent internal Bloch eigenmodes.

KEYWORD
Fourier modal method , Rigorous coupled wave analysis , Surface plasmons , Gratings

I. INTRODUCTION
The Fourier modal method (FMM) is one of the most popular numerical electromagnetic analysis methods for periodic optical structures [1]. With recent developments of the initial theoretical framework, such as perfect matched layer (PML) [2, 3] and the generalized scatteringmatrix method [4, 5], the range of applications of the FMM was extended to include aperiodic, anisotropic device structures and more generalized photonic network structures [5]. Furthermore, with the performance improvement of modern computers, the FMM has become a crucial tool in the development of various optical and optoelectronic devices, such as antireflection coatings [6], broadband reflectors [7,8], narrowband bandpass filters [8], solar cells [9], waveguides [2, 10], and optomechanical devices [11].
Perhaps the most important development in the history of the FMM is Li’s Fourier factorization rule for improving the convergence of the FMM in analyzing metallic gratings under transversemagnetic (TM) polarization [1214]. Since the FMM with Li’s factorization rule shows good convergence in analyzing such structures, it has become one of the most important analysis tools in the field of plasmonics [15, 16]. However, through FMM analyses of metallic structures under surface plasmon resonance (SPR), we have found that the coefficients of transmission and reflection have peculiar small fluctuations when plotted versus the number of Fourier terms used, even though the FMM has a large number of Fourier terms and is featured with Li’s factorization rule.
This motivates the investigation of the origins of the fluctuations observed in the FMM analysis of a onedimensional lamellar metallic grating with discontinuities in the permittivity profile. An indepth discussion on the internal modal structure of the FMM and the Gibbs phenomenon arising in the Fourier representation of discontinuous metallic structures is presented. To systematically study this problem, we modify the conventional formulation of the FMM to control the number of Fourier harmonics representing the fields independently of that of Fourier terms associated with the permittivity profile, and examine the respective effects of the numbers of the field and structure harmonics on the convergence behavior of the FMM. Based on the analysis, we discuss the origins of the fluctuations in the convergence of the FMM, and conclude with perspective and a future plan for enhancing the accuracy of the FMM.
This paper is organized as follows. In Sec. 2, the numerical framework of the FMM is reviewed. In Sec. 3, the main problem of interest is stated with an example of the FMM analysis. The reason for the problem is explained within the conventional FMM framework by discussing the Gibbs phenomenon and its physical relation to fluctuations. In Sec. 4, we modify the conventional FMM by decoupling the linkage between the numbers of the field and structure harmonics, examine the resulting convergence behaviors of reflection and transmission efficiencies of a metallic grating, and clarify the origins of the fluctuations, followed by concluding remarks in Sec. 5.
II. NUMERICAL FRAMEWORK OF THE FOURIER MODAL METHOD
Before investigating the convergence behaviors of the FMM for electromagnetic simulations of a metallic grating structure, we briefly review the numerical framework of the FMM [1]. For a onedimensional grating structure investigated in the following sections (Fig. 1), the permittivity in the grating layer,
ε (x ), is approximated by a truncated Fourier series,where
G =2/Λ, and Λ is the grating period. In the FMM, the internal electric and magnetic eigenmodes,E ^{(g)} andH ^{(g)}, in the grating structure take the form of a truncated pseudoFourier series:where
e ^{jk(g)z} is the Bloch eigenphase term, (E_{x} ^{(g)})m , (E_{y} ^{(g)})m , (E_{z} ^{(g)})m and (H_{x} ^{(g)})m , (H_{y} ^{(g)})m , (H_{z} ^{(g)})m are the Fourier coefficients of the electric and magnetic field envelope functions, respectively, andk_{x} andk_{y} are determined by the wave vector of the incident plane wave. The Fourier coefficients and Bloch eigenphase terms are calculated by substituting Eqs. (1) and (2) into Maxwell’s equations and solving the resulting algebraic eigenvalue equation [1]. In doing so for TMpolarized (magnetic field parallel to they axis) waves, the Fourier coefficients of the axial (z directional) electric displacement vector,ε (x )E_{z} ^{(g)}, are given by the Laurent’s rule:where？
ε ？_{l,m} are the (l ,m ) matrix element of the Toeplitz matrix？ε ？ generated from the Fourier coefficients ofε (x ),ε_{m} . The Fourier coefficients of the transversal (x directional) electric displacement vector,ε (x )E_{x} ^{(g)}, are obtained in a different fashion, known as Li’s rule [12, 13]:where
are the (
l ,m ) matrix elements of the inverse of the Toeplitz matrix generated from the Fourier coefficients of 1/ε (x ). Li’s factorization rules are critical for fast convergence of the FMM for the case of metallic gratings under TM illumination [13]. In the case of transverseelectric (TE, electric field parallel to they axis) waves, the Fourier coefficients ofε (x )E_{y} ^{(g)} are obtained as in Eq. (3).In the conventional formulation of the FMM, the 4(2
M +1) internal Bloch eigenmodes are obtained numerically, and the total electromagnetic field inside the grating structure is described by the coherent superposition of 4(2M +1) internal Bloch eigenmodes aswhere
C_{g} is the coupling coefficient that represents the degree of coupling between the external excitation and theg^{th} internal electromagnetic Bloch eigenmode, and is matched by electromagnetic boundary conditions.III. FOURIER MODAL ANALYSIS OF A ONEDIMENSIONAL METALLIC GRATING
A grating structure that we analyze to examine the convergence behaviors of the FMM is illustrated in Fig. 1. An optical plane wave from the freespace region is
incident normal to the
xy plane on a substrate with an index of refraction (n ) of 1.75, occupying the lower half space, in which a Ag periodic binary grating, with a period (Λ ) of 250 nm, a fillfactor of 0.25, and a thickness of 20 nm, is embedded. Figure 2(a) shows the spectra of external transmission and reflection efficiencies (T _{TM} andR _{TM}), calculated using the FMM, for a TMpolarized incident plane wave, whose wavelengthλ is varied from 400 nm to 800 nm. The results calculated for a TEpolarized plane wave (T _{TE} andR _{TE}) are shown in Fig. 2(b). The calculations were performed withM = 500, and the complex permittivity of Ag is taken from the literature [17].The
R _{TM} spectrum has the main peak atλ =520 nm and a small shoulder atλ =420 nm, both of which are attributed to SPRs. In contrast, theT _{TE} andR _{TE} spectra are nearly featureless [Fig. 2(b)]. In Figs. 3(a) and 3(b), the TM and TE electric field intensity profiles forλ =520 nm are shown, respectively. For TM polarization, the electric field is tightly concentrated on the sidewalls and the bottom surface of the grating due to SPRs. In contrast, the field enhancements near the grating boundaries are not observed for the case of TE polarization, as shown in Fig. 3(b).Next, we investigate the convergence characteristics of the FMM calculation. In Fig. 4,
T _{TM} [(a)],R _{TM} [(b)],T _{TE} [(c)], andR _{TE} [(d)] are plotted as a function ofM for the three selected wavelengths, 420 nm (blue), 520 nm (red), and 750 nm (black). For the TM case, the results for all three wavelengths seem to be converged for sufficiently largeM (>100) but close examinations reveal that theyhave fluctuations, as shown in Figs. 4(a) and 4(b). In contrast, the results for the TE cases converge much faster, and their variations in
M > 100 are negligible [Figs. 4(c) and 4(d)].Considering the field profiles shown in Fig. 3, the contrast in convergence characteristics between the TM and TE cases suggests that the fluctuations in
T _{TM} andR _{TM} are related to the inherent inaccuracies in the Fourierseries approximation of the permittivity profile near the grating boundaries. Figure 5 shows the exact permittivity profile in the grating layer (black),ε (x ), compared with its Fourierseries approximation,calculated at
λ =520 nm by Eq. (1) forM = 500 (blue) and 1000 (red). Theprofiles have oscillations near the grating boundary, and in each (metal or dielectric) region, the oscillation occurring nearest to the grating boundary has the largest magnitude ？ the characteristic of the Fourierseries representation of a piecewise continuous function with jump discontinuities, known as the Gibbs phenomenon [18]. This observation provides a qualitative and physical explanation that the fluctuations in
T _{TM} andR _{TM} are, at least in part, ascribed to nonnegligible responses of edgeconcentrated electromagnetic fields to the fluctuation innear discontinuous grating boundaries due to the Gibbs phenomenon. Consequently,
T _{TM} andR _{TM} show the fluctuations even withM = 270 and the correct Fourier factorization rules, as shown in Figs. 4(a) and 4(b). Meanwhile, the Gibbs phenomenon does not result in fluctuations inT _{TE} andR _{TE}, since the extent of regions near the grating boundaries whereand
ε (x ) differ is very small compared withΛ (= 250 nm), and the field concentration near the grating boundaries does not occur for TEpolarized illumination [Fig. 3(b)].IV. ORIGINS OF THE FLUCTUATIONS IN THE CONVERGENCE BEHAVIOR OF THE FMM ANALYSIS
We examine the numerical framework of the FMM and focus on its internal modal structure to identify the origin of the fluctuations in
T _{TM} andR _{TM}. In the conventional FMM framework, 4(2M +1) Bloch eigenmodes are numerically obtained with (2M +1) field harmonics, where the coupling between field harmonics is specified by (4M +1)ε _{m} ’s [1]. Here, it is noted that the number of structure harmonics in Eq. (1) (= 4M + 1) and that of field harmonics in Eq. (2) are linked and varied together.Upon increasing
M in the FMM framework, many highly evanescent Bloch eigenmodes are generated, and the degree of evanescence of each mode can be evaluated by the magnitude of the imaginary part of the Bloch eigenvaluek_{z} ^{(g)} [1]. By examining the convergence ofk_{z} ^{(g)}’s with respect toM , we find that for givenM , the calculated 4(2M + 1) internal Bloch eigenmodes can be categorized into (i) lowerorder nonevanescent or weakly evanescent Bloch modes that are converged, in a sense that their eigenvalues and field profiles are identical to those obtained with smallerM ’s, and (ii) a few highly evanescent modes that are not converged yet. The coupling coefficients, in expressing the total field profile as in Eq. (5), of the internal Bloch eigenmodes in categories (i) and (ii) are denoted byC_{g} ^{(i)}’s, andC_{g} ^{(ii)}’s, respectively.The convergence of the eigenmodes in category (i) means that
k_{z} ^{(g)}’s and the corresponding field profiles of these eigenmodes are not affected by the oscillating change innear the grating edges as
M is further increased. Meanwhile, for the eigenmodes in category (ii) to be converged, more higherorder Fourier terms need to be included by increasingM . Then, the values ofk_{z} ^{(g)} for newly generated higherorder terms are different from their converged values; in other words, for anyM , there are always eigenmodes belonging to category (ii). This finding strongly suggests that the fluctuations inT _{TM} andR _{TM} are mainly ascribed to the Bloch eigenmodes in category (ii). More specifically, for TMpolarized illumination,C_{g} ^{(ii)}’s are not negligible. In viewpoint of the numerical framework, the internal mismatch between the number of field harmonics used (= 2M + 1) and that required to accurately represent all 4(2M + 1) Bloch eigenmodes causes parameters, such asT _{TM} andR _{TM}, to exhibit fluctuations. In contrast,C_{g} ^{(ii)}’s in the case of TE polarization are negligible, meaning that the unconverged evanescent modes do not contribute to the total field distribution. Consequently,T _{TE} andR _{TE} converge without fluctuating features despite the harmonic number mismatch.To address this issue more clearly, we perform numerical experiments by modifying the conventional FMM framework, where the number of structure harmonics (= 4
M + 1) representingε (x ) is directly linked to that of field harmonics (= 2M + 1) as seen in Eqs. (1) and (2), to decouple these two numbers. Specifically,ε (x ) is approximated byand the electromagnetic fields are given by Eq. (2), where
M andN can be independently varied. IfM is increased sufficiently beyond the value in the conventional FMM (=N / 2) for fixedN , we expect that the evanescent Bloch eigenmodes with nonnegligible coupling coefficients will be converged, and consequently theT _{TM} andR _{TM}vs.M plots are expected to converge without fluctuations. This procedure may be considered as bandlimitingHowever, we emphasize that it is a means to provide a sufficient number of field harmonics required to accurately represent Bloch eigenmodes involved in diffraction by the approximate grating profile
To show this, we plot
R _{TM} of the onedimensional grating analyzed in Sec. 3 under SPR (λ = 520 nm) as a function ofM withN equal to its value in the convenional FMM (= 2M ) only until a certain value (Fig. 6). Figure 6(a) shows the result when the maximumN (≡N _{max}) is 40, showing that quick convergence ofR _{TM} is obtained with significantly reduced fluctuations, as evident in the inset. This indicates that ifM > ~300,k_{z} ^{(g)} and the corresponding field profiles of Bloch eigenmodes with nonnegligible coupling with the incident wave are completely converged, and this bandlimiteddoes not couple the incident
wave with the highly evanescent eigenmodes that are not converged yet. Figure 6(b) shows the result for
N _{max} = 100, where the height of the envelope within whichR _{TM} oscillates becomes larger. We also note that the values ofR _{TM} in Figs. 6(a) and 6(b) converge to different values, since the permittivity profile is slightly different, near the grating boundaries in particular, in each simulation.These results indicate that calculations of
R _{TM} andT _{TM} as a function ofN for sufficiently largeM are free from errors due to unconvergedk_{z} ^{(g)} of highly evanescent modes, therefore demonstrating fluctuations attributed to the Gibbs phenomenon only. To confirm this, we plot, in Fig. 7,R _{TM} andT _{TM} of the grating under the same illumination condition as in Fig. 6, as functions ofN forM ？600: for eachN , envelopes outlining oscillations inR _{TM} andT _{TM}vs.M plot, such as those shown in Fig. 6, are calculated, and the values of the upper and lower envelopes atM = 600 are marked as red and black dots in Fig. 7. Also shown is the mean of the upper and lower envelopes atM = 600, marked as blue dots. ForN < ~200, the upper and lower envelopes of bothR _{TM} andT _{TM} are indistinguishable, meaning that they are converged with respect toM for eachThis indicates that
M (= 600) is sufficiently large to accurately evaluatek_{z} ^{(g)} of Bloch eigenmodes whose coupling to the incident plane wave is nonnegligible forwith
N < ~200, and that the very small fluctuations forN < ~200 are due to variation inHowever, as
N is increased beyond ~200, the number of field harmonics (= 2M +1？1200) becomes insufficient, resulting in separation between the upper and lower envelopes. This means that both types of errors ？ one in representing trueε (x ) and the other in evaluatingk_{z} ^{(g)} of highly evanescent Bloch eigenmodes with nonnegligible coupling coefficients ？ contribute to the fluctuations. The magnitude of fluctuations at eachN in Fig. 7, i.e., the difference between the upper and the lower envelopes, is much smaller than that of fluctuations aroundM =N /2 in Fig. 4, which indicates that the errors arising from coupling with unconverged higher order evanescent eigenmodes are reduced in the modified FMM. Although the values ofM in our simulations are limited, we infer that ifN approaches to infinity while keepingM ≫N /2, the modified FMM will converge with negligible, or much smaller, fluctuations compared to the expected asymptotic behavior of the conventional FMM.V. CONCLUSION
We have unveiled that two types of errors are present in the FMM analysis of a onedimensional metallic gratings: one in representing true
ε (x ) by a Fourier series with a finite number of Fourier terms, which is related to the Gibbs phenomenon, and the other in evaluating field profiles of highly evanescent Bloch eigenmodes. We have found that the number of field harmonics in the FMM framework, 2M + 1, is always insufficient to accurately represent a few most highly evanescent eigenmodes, and it is an inherent structural limitation of the conventional FMM framework with Li’s Fourier factorization rule. We have shown that the fluctuations arising from inaccurate calculations of these eigenmodes can be eliminated in the modified FMM, by increasingM sufficiently beyond the value in the conventional FMM (=N / 2)Due to highly concentrated fields near the grating boundaries, TMpolarized cases are sensitive to this inaccuracy, resulting in fluctuations in
R _{TM} andT _{TM} even with largeM . We expect that fluctuations similar to what is observed here to be present in dealing with general metallic structures supporting SPRs, whenever the electric fields are highly concentrated in metal around its discontinuous boundaries. Moreover, in analyzing the electric field enhancement in a deepsubwavelength structure, for example, a metallic slit whose width is significantly smaller than the wavelength of interest [19], the coupling between the excitation field and veryhighorder evanescent eigenmodes must be considered. Since the errors investigated in this paper are the inherent characteristics of the conventional FMM, its application to deepsubwavelength structures can be inaccurate. To obtain accurate results in such cases using the FMM with finite computing resources, the method of representing permittivity profiles may need to be modified. A good strategy may be a form of truncated Fourier series with a specific apodization in the spatialfrequency domain. This topic will be studied based on what is established in this paper.

[FIG. 1.] Schematic of a onedimensional metallic grating embedded in a substrate occupying the lower half space.

[FIG. 2.] (a) Transmission (TTM) and reflection (RTM) efficiencies for TMpolarized illumination. (b) Transmission (TTE) and reflection (RTE) efficiencies for TEpolarized illumination.

[FIG. 3.] The electric field intensity profiles for TM (a) and TE (b) polarization when λ = 520 nm.

[FIG. 4.] Convergence behavior of transmission and reflection efficiencies for TM and TEpolarized plane waves with λ = 420 nm (blue), 520 nm (red), and 750 nm (black). (a) TTM vs. M for 270≤M≤300 (b) RTM vs. M for 270≤M≤300. (c) TTE vs. M for 100≤M≤300. (d) RTE vs. M for 100≤M≤300.

[FIG. 5.] The exact permittivity profile (black), ε(x), compared with its Fourierseries approximation, ε(x), with M = 500 (blue) and 1000 (red). Both real (left) and imaginary (right) parts are shown.

[FIG. 6.] Reflection efficiencies (red) versus M when the number of structure harmonics (2N+1) is limited with (a) N = 40 and (b) N = 100 .

[FIG. 7.] (a) Reflection and (b) transmission efficiencies versus N, when M？600. For each N envelopes outlining oscillations in RTM  and TTM vs.M plots, such as those shown in Fig. 6, are calculated, and the values of the upper and lower envelopes at M = 600 are marked as red and black dots. The blue dots represent the mean of the upper and lower envelopes.