# New reformulation of the Fourier modal method with spatial adaptive resolution for multilayered metallic strip grating

H. Elamine, B. Guizal, Meherzi Oueslati, T. Gharbi

► To cite this version:

H. Elamine, B. Guizal, Meherzi Oueslati, T. Gharbi. New reformulation of the Fourier modal method with spatial adaptive resolution for multilayered metallic strip grating. Optics Communications, 2010, 283 (21), pp.4392-4396. 10.1016/j.optcom.2010.06.044. hal-00535894

HAL Id: hal-00535894

https://hal.science/hal-00535894v1

Submitted on 6 May 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# New reformulation of the Fourier modal method with spatial adaptive resolution for multilayered metallic strip grating

H. Elamine<sup>a,c,*</sup>, B. Guizal<sup>b</sup>, M. Oueslati<sup>a</sup>, T. Gharbi<sup>c</sup>

The parametric formulation of the combined boundary conditions method (CBCM) with spatial adaptive resolution is extended to multilayered structures of strip gratings. Furthermore, it is shown that it is not necessary to solve the eigenvalue problem in all the layers of these structures; leading to drastic reduction of the computational load.

## 1. Introduction

Among the methods used to model diffraction gratings, there is a general class dedicated to perfect strip gratings that are assumed to be infinitely thin and perfectly conducting (PEC). This approach is very interesting and very competitive when one is dealing with a problem supporting such assumptions, as is the case for metallic gratings in the microwave and terahertz frequency domains. Thus, the concept has been successfully applied for many applications ranging from photolithographic masks to frequency selective surfaces. Besides, one of the most versatile and efficient methods tailored for perfect strip gratings is the CBCM [1]. The principle of this differential approach is the combination of the continuity relations of the electric and magnetic fields in unified equations that hold over one period of the structure. Thus all the continuity relations can be projected over the Fourier basis permitting straightforward calculations. Since its introduction in [1], this method has been extended to bulk perfectly metallic gratings [2] and successfully adapted to aperiodic gratings [3], circular strip cavities [4] and curved strip gratings [5]. Later on, with the discovery of spatial adaptive resolution (ASR) [6], the CBCM has been substantially improved [7] in the sense that it was rid of the problems inherent to Fourier series (slow convergence and Gibbs phenomenon). In [8], the ASR has been extended to multilayered structures in the context of the Fourier modal method and very

recently improved further [9]. In the present paper, we (i) extend the CBCM equipped with the ASR to multilayered strip gratings and (ii) show that the additional computational steps introduced by the ASR (time consuming eigenlaupe problems) can be drastically reduced.

## 2. Theory

The structure under study is depicted in Fig. 1, it consists of a set of Q homogeneous dielectric layers (with heights $h_q$ and dielectric permittivities $\epsilon_q$, $1 \le q \le Q$) separated by Q + 1 strip gratings sharing the same period $d$ and lying between two homogeneous dielectric media with permittivities $\epsilon_{in}$ and $\epsilon_{out}$. The strip gratings are not necessarily aligned. The structure is invariant along the z direction and the time dependence will be expressed (and suppressed through the paper) through the factor $e^{-i\omega t}$.

The device is illuminated from the lower part by a monochromatic plane wave under the angle of incidence $\theta$, $\lambda$ stands for the vacuum wavelength and $k_0 = 2\pi/\lambda$ for its corresponding wavenumber.

As stated above, we will use the concept of ASR in order to remove convergence problems stemming from the singular behavior at the edges of the strips. In our implementation, we opt for the following coordinate transformation $x(u)$ taken from [8]

$$x_l(u) = a_1 + a_2 u + \frac{a_3}{2\pi} \sin \left( 2\pi \frac{u-u_{l-1}}{u_l-u_{l-1}} \right) \quad (1)$$

$$a_1 = \frac{u_l x_{l-1} - u_{l-1} x_l}{u_l - u_{l-1}} \quad (2)$$

<footnote><sup>a</sup> Unité de Recherche Spectroscopie Raman - Faculté des Sciences de Tunis Campus Universitaire El Manar 1-Tunis-1060, Tunisie</footnote>

<footnote><sup>b</sup> Equipe de Nanophotonique, Groupe d'Etude des Semiconducteurs UMR 5650, Bât. 21, CC074 Place Bataillon, 34095 Montpellier Cedex 05, France</footnote>

<footnote><sup>c</sup> FEMTO-ST, UMR CNRS No. 6174, Route de Gray 25030 Besaçnon Cedex, France</footnote>

<footnote>* Corresponding author. FEMTO-ST, UMR CNRS No. 6174, Route de Gray 25030 Besaçnon Cedex, France.</footnote>

<footnote>E-mail addresses: helamine@univ-fcomte.fr (H. Elamine), guizal@ges.univ-montp2.fr (B. Guizal).</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Multilayered strip gratings.</figcaption></figure>

$$a_2 = \frac{x_l - x_{l-1}}{u_l - u_{l-1}} \qquad (3)$$

$$a_3 = G(u_l - u_{l-1}) - (x_l - x_{l-1}) \qquad (4)$$

where $G=f(u_{l-1})=f(u_l)$.

It is important at that stage to emphasize that we manage to use the same transformation for all the layers. In order to solve the problem of diffraction, one has to solve Maxwell's equations in all the homogeneous media and then match the boundary conditions at the interfaces between them. In the present case and due to the z invariance, the problem reduces to the two classical cases of polarization TE (electric field parallel to the strips) and TM (magnetic field parallel to the strips). For homogeneous media, the wave equation has been shown to have the same expression for both of these polarizations [9]:

$$\frac{1}{f} \frac{\partial}{\partial u} \left( \frac{1}{f} \frac{\partial \Psi_q}{\partial u} \right) + \frac{\partial^2 \Psi_q}{\partial y^2} + k_0^2 \varepsilon_q \Psi_q = 0 \qquad (5)$$

where $\Psi_q(u,y)$ stands for $E_{zq}(u,y)$ or $H_{zq}(u,y)$ in the qth layer and $f=\partial x/\partial u$. Written in the Fourier space, and using matrix notations, this yields:

$$\frac{\partial^2 \hat{\Psi}_q}{\partial y^2} = \left( [[f]]^{-1} \alpha [[f]]^{-1} \alpha - k_0^2 \varepsilon_q I_d \right) \hat{\Psi}_q = A_q \hat{\Psi}_q \qquad (6)$$

where $\alpha = (k_0\sqrt{\varepsilon_{in}}\sin(\theta) + 2\pi n / d)$, $n \in \mathbb{Z}$, $I_d$ is the identity matrix and $\hat{\Psi}_q$ is a vector containing the Fourier coefficients of $\Psi_q(u,y)$ and $[[f]]$ is a matrix whose elements are given by $[[f]]_{mn}=f_{m-n}$; the $(m-n)$th Fourier coefficient of function $f$.

Matrix $A_q$ can be written under the form $A_q = P_q D_q^2 P_q^{-1}$ where $D_q^2$ is a diagonal matrix containing the eigenvalues of $A_q$, and $P_q$ is the matrix of the corresponding eigenvectors. Then the general solution can be written:

$$\hat{\Psi}_q(y) = P_q \left( e^{(y-y_{q-1})D_q} a_q + e^{-(y-y_{q-1})D_q} b_q \right) \qquad (7)$$

In order to express the solution in each layer, it is clear, from Eq. (7), that one needs to solve the eigenvalue problem associated to this layer. Thus, for our problem, we have to solve $Q+2$ eigenvalue

problems (Q layers and the in and out semi-infinite media). This represents a non negligible supplementary computational effort compared to the classical CBCM. But as (i) we are dealing with homogeneous layers and (ii) we use the same coordinate transformation for all the layers, it is possible to reduce drastically the computational load as has been shown recently [9]. The key idea is the link between the matrices of the input, the output region and the different layers:

$$\left\{ \begin{aligned} A_{in} &= [[f]]^{-1} \alpha [[f]]^{-1} \alpha - k_0^2 \varepsilon_{in} I_d \\ A_q &= [[f]]^{-1} \alpha [[f]]^{-1} \alpha - k_0^2 \varepsilon_q I_d \\ A_{out} &= [[f]]^{-1} \alpha [[f]]^{-1} \alpha - k_0^2 \varepsilon_{out} I_d \end{aligned} \right. \qquad (8)$$

All the matrices share the term $[[f]]^{-1} \alpha [[f]]^{-1} \alpha$ and one can write:

$$A_{in} + k_0^2 \varepsilon_{in} I_d = A_q + k_0^2 \varepsilon_q I_d = A_{out} + k_0^2 \varepsilon_{out} I_d \qquad (9)$$

These last expressions show that it is sufficient to solve numerically the eigenvalue problem only for a given medium and then deduce the solution of the others from it. Suppose that we have the solution of the eigenvalue of the input region through the eigenvalues $D_{in}^2$ and the associated eigenvectors $P_{in}$, then left multiplying Eq. (9) by $P_{in}^{-1}$ and right multiplying by $P_{in}$ gives:

$$D_{in}^2 + k_0^2 \varepsilon_{in} I_d = P_{in}^{-1} A_q P_{in} + k_0^2 \varepsilon_q I_d = P_{in}^{-1} A_{out} P_{in} + k_0^2 \varepsilon_{out} I_d \qquad (10)$$

and finally:

$$\left\{ \begin{aligned} P_{in}^{-1} A_q P_{in} &= D_{in}^2 + k_0^2 (\varepsilon_{in} - \varepsilon_q) I_d \\ P_{in}^{-1} A_{out} P_{in} &= D_{in}^2 + k_0^2 (\varepsilon_{in} - \varepsilon_{out}) I_d \end{aligned} \right. \qquad (11)$$

These expressions clearly show that all the matrices $A_q$ and $A_{out}$ share the same eigenvectors $P_{in}$ and that their eigenvalues can be, easily, obtained through the simple translations:

$$\left\{ \begin{aligned} D_q^2 &= D_{in}^2 + k_0^2 (\varepsilon_{in} - \varepsilon_q) I_d \\ D_{out}^2 &= D_{in}^2 + k_0^2 (\varepsilon_{in} - \varepsilon_{out}) I_d \end{aligned} \right. \qquad (12)$$

Now that we have derived the solutions in all the regions, we must write the boundary conditions over the interfaces. For this purpose, we are going to treat separately the two cases of polarization TE and TM and formulate the solution trough the S-matrix algorithm.

## 2.1. TE polarization

Under this polarization, the application of the CBCM concept leads to the following continuity equations, that we will express for each interface:

<table><caption>**Table 1**
Reflected and transmitted efficiencies calculated for TE and TM cases for $M=128$.</caption><thead><tr><th></th><th>R<sub>-2</sub></th><th>R<sub>-1</sub></th><th>R<sub>0</sub></th><th>T<sub>0</sub></th></tr></thead><tbody><tr><td>TE</td><td>0.0256</td><td>0.1801</td><td>0.3711</td><td>0.2174</td></tr><tr><td>TM</td><td>0.0256</td><td>0.1801</td><td>0.2175</td><td>0.3711</td></tr></tbody></table>

$c=0.5$, $d=1$, $\theta=26^\circ$, $\lambda=0.7$, $g_{TE}=-10^{-3}$, $g_{TM}=377i$ truncation order $M=128$.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Reflected efficiencies in the orders 0, -1, -2, for the TE case versus the thickness of the layer h, λ = 0.7, θ = 30°.</figcaption></figure>

Interface $y=y_q$:

$$ \forall 0 \le u \le d \left\{ \begin{array}{l} E_q(u, y_q) = E_{q+1}(u, y_q) \\ \chi_q(u)E_q(u, y_q) + g\tilde{\chi}_q(u)f(\partial_y E_{q+1}(u, y_q+1) - \partial_y E_q(u, y_q)) = 0 \end{array} \right. \quad (13) $$

where g is an arbitrary constant, $\chi_q$ the characteristic function of the set of strips lying in the interface given by $y=y_q$ and $\tilde{\chi}_q(u) = 1 - \chi_q(u)$.

$$ \chi_q(u) = \begin{cases} 1 & \text{over the metal} \\ 0 & \text{elsewhere} \end{cases} \quad (14) $$

Expressed in the Fourier space, the continuity equations become:

$$ \begin{cases} P_q(\phi_q a_q + \phi_q^{-1} b_q) = P_{q+1}(a_{q+1} + b_{q+1}) \\ \chi_q P_q(\phi_q a_q + \phi_q^{-1} b_q) + g f \tilde{\chi}_q (P_{q+1} D_{q+1}(a_{q+1} - b_{q+1}) - P_q D_q(\phi_q a_q - \phi_q^{-1} b_q)) = 0 \end{cases} \quad (15) $$

where $\phi_q = \text{diag}(e^{h_q D_q})$ then the S-matrix of this interface is given through:

$$ \begin{pmatrix} b_q \\ a_{q+1} \end{pmatrix} = S_{\text{TE}} \begin{pmatrix} a_q \\ b_{q+1} \end{pmatrix} = \begin{pmatrix} \phi_q & 0 \\ 0 & I_d \end{pmatrix} S_{\text{TE}} \begin{pmatrix} \phi_q & 0 \\ 0 & I_d \end{pmatrix} \begin{pmatrix} a_q \\ b_{q+1} \end{pmatrix} \quad (16) $$

where

$$ S_{q\text{TE}} = \begin{pmatrix} P_q & -P_{q+1} \\ \chi_q P_q + g\chi_q f(P_q D_q) & g\tilde{\chi}_q f P_{q+1} D_{q+1} \end{pmatrix}^{-1} \quad (17) \\ \times \begin{pmatrix} -P_q & P_{q+1} \\ g\chi_q f P_q D_q - \chi_q P_q & g\tilde{\chi}_q f P_{q+1} D_{q+1} \end{pmatrix} $$

<figure><img src="image_4.png" /><figcaption>Fig. 3. The zero-order efficiency on TE mode versus the ratio $\lambda/d$, $\theta = 30°$, the grating dimension is defined in terms of the groove period $d$: $c = d/2$ and $h = d/4$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Computation time versus the number of layers for M = 35.</figcaption></figure>

## 2.2. TM polarization

For this polarization, the CBCM boundary conditions are:

$$ \forall 0 \le u \le d \left\{ \begin{array}{l} \frac{1}{\varepsilon_q} \partial_y H_q(u, y_q) = \frac{1}{\varepsilon_q + 1} \partial_y H_{q+1}(u, y_q) \\ \chi_q \frac{f}{\varepsilon_q} \partial_y H_q(u, y_q) + g \tilde{\chi}_q (H_q(u, y_q + 1) - H_q(u, y_q)) = 0 \end{array} \right. \quad (18) $$

following the same procedure as for the TE case, we obtain:

$$ \begin{pmatrix} b_q \\ a_{q+1} \end{pmatrix} = S_{\text{TM}} \begin{pmatrix} a_q \\ b_{q+1} \end{pmatrix} = \begin{pmatrix} \phi_q & 0 \\ 0 & I_d \end{pmatrix} S_{q\text{TM}} \begin{pmatrix} \phi_q & 0 \\ 0 & I_d \end{pmatrix} \begin{pmatrix} a_q \\ b_{q+1} \end{pmatrix} \quad (19) $$

with

$$ S_{q\text{TM}} = \begin{pmatrix} \frac{1}{\varepsilon_q} P_q D_q & \frac{1}{\varepsilon_q + 1} P_{q+1} D_{q+1} \\ \chi_q \frac{f}{\varepsilon_q} P_q D_q + g \chi_q P_q & -g \chi_q + 1 P_{q+1} \end{pmatrix}^{-1} \quad (20) \\ \times \begin{pmatrix} \frac{1}{\varepsilon_q} P_q D_q & \frac{1}{\varepsilon_q + 1} P_{q+1} D_{q+1} \\ \chi_q \frac{f}{\varepsilon_q} P_q D_q - g \chi_q P_q & g \chi_q P_{q+1} \end{pmatrix} $$

At this stage, in order to implement the method based on the above equations (TE and TM) it's necessary to truncate the Fourier series. In the numerical implementations the series will be truncated to order M which means that only 2M + 1 Fourier Harmonics are kept.

## 3. Numerical results

In this section we provide numerical examples to demonstrate the reliability of the method and the improvements in terms of numerical speed. We constructed numerical codes based on the formalism presented in the last section and checked that they verify the classical criteria of reciprocity and energy conservation. We also checked that we retrieve the results for single strip gratings by using the multilayered strip grating codes. For instance we replicated the original single grating to several stacks and we let the thicknesses between them tend to zero. More concretely, we consider the configuration used in Ref. [2] and recall, in Table 1, the values of the efficiencies $R_{-2,-1,0}$ and $T_0$. Fig. 2 gives the evolution, in a semilogarithmic scale, of $R_{-2,-1,0}$ versus the thickness h between

each couple of strip gratings (cf. inset in the figure). It can be seen that as h is decreased, the results tend to the values of the single strip grating. The presented computations are performed with three stacks of strip gratings but we have verified that it works for larger numbers of stacked gratings.

As a second example, we consider the case of a lamellar PEC grating illuminated under TE polarization. It has been shown in [10] that one efficient way to model diffraction from PEC surface relief gratings, is to cascade several PEC strip gratings following the original shape (staircase approximation). This approach is valid only in the case TE of polarization and is based on the so called filtering effect: the fact that the electric field does not penetrate inside the very thin waveguides formed by two consecutive strips. In Fig. 3, we reproduce the zero order reflected efficiency diffracted from a lamellar PEC grating as published in [11] and later in [12] where the computations were performed by means of a modal method and successfully compared to experiments (the different parameters are given in the figure caption). These results are in excellent agreement with those of Refs. [11,12] and validate our approach indirectly against experiments. It is important to stress that, in this case, as few as 20 layers were necessary to obtain this agreement.

Now, after validating our generalization of the ASR-CBCM to multilayered structures, we turn to the second major issue of this work: the computational time reduction obtained through the solution of only one eigenvalue problem. It is expected that time reduction will scale almost linearly with the number of layers; as the time spent to compute eigenvalues through Eq. (12) is negligible compared to fully solving the corresponding eigenvalue problem. To numerically demonstrate these statements, we consider the case of the lamellar PEC grating treated above by the multilayer approach. In Fig. 4 we report the evolution of the commutation time versus the number of layers used to compute a set of efficiencies for a given wavelength; for instance $\lambda = d$. For these computations, the truncation order M has been fixed to 35. It can be seen, clearly, that the new approach is faster than the old one. Moreover the speedup grows with the numbers of layers. It is important to add that the computation time difference between the new and the old approaches as drawn in Fig. 4 will be more pronounced for a higher truncation order M.

## 4. Conclusion

A generalization of the CBCM equipped with the ASR to multiple stacks of strip gratings has been presented. Moreover it has been shown that it is not necessary to solve the eigenvalue problem, inherent to the ASR concept, in all the layers which leads to substantial computation time saving. The approach has been validated against previously published theoretical and experimental results.
Acknowledgments

This work is realized in the framework of the research project
CNRS-DGRS 08/13-17.

References

[1] F. Montiel, M. Neviere, Opt. Commun. 101 (1993) 151–156.

[2] F. Montiel, M. Neviere, Opt. Commun. 138 (1997) 172–184.

[3] B. Guizal, D. Felbacq, Opt. Commun. 165 (1999) 1–6.

[4] B. Guizal, D. Flebacq, Phys. Rev. E 66 (2002) 026602.

[5] B. Guizal, G. Granet, J. Opt. Soc. Am. A Vol. 24 (No.3) (2007) 669.

[6] G. Granet, J. Opt. Soc. Am. A Vol. 16 (No.10) (1999) 2510–2516.

[7] G. Granet, B. Guizal, Opt. Commun. 255 (2005) 1–11.

[8] T. Vallius, M. Honkanen, Opt. Express Vol. 10 (No.1) (2002) 14.

[9] H. Yala, B. Guizal, D. Felbacq, J. Opt. Soc. Am. A Vol. 26 (No.12) (2009) December.

[10] F. Montiel, M. Neviere, Volume 144, Issues 1–3, 1 December 1997, 82–88.

[11] D. Maystre, M. Neviere, R. Petit, Electromagnetic Theory of Gratings, in: R. Petit (Ed.), Top. Curr. Phys., Vol. 22, Springer-Verlag, Berlin, 1980, p. 159.

[12] Yon-Lin Kok, Neal C. Gallagher, J.Opt.Soc. A Vol. 5 (N 1) (1988) 65–73 January.