Curvilinear coordinate Generalized Source Method for gratings with sharp edges

TL;DR

This paper introduces a new curvilinear coordinate Generalized Source Method for gratings with sharp edges, enhancing computational accuracy and efficiency in diffractive and metasurface optics applications.

Contribution

It extends the Generalized Source Method to handle sharp edges in gratings by treating corners as effective media interfaces, improving accuracy and computational efficiency.

Findings

The method achieves higher accuracy for gratings with sharp edges.

It maintains $O(N \, log N)$ complexity, ensuring efficiency.

Corners are effectively modeled as interfaces, simplifying calculations.

Abstract

High-efficient direct numerical methods are currently in demand for optimization procedures in the fields of both conventional diffractive and metasurface optics. With a view of extending the scope of application of the previously proposed Generalized Source Method in the curvilinear coordinates, which has theoretical $O\left(N\log N\right)$ asymptotic numerical complexity, a new method formulation is developed for gratings with sharp edges. It is shown that corrugation corners can be treated as effective medium interfaces within the rationale of the method. Moreover, the given formulation is demonstrated to allow for application of the same derivation as one used in classical electrodynamics to derive the interface conditions. This yields continuous combinations of the fields and metric tensor components, which can be directly Fourier factorized. Together with an efficient algorithm the new formulation is demonstrated to substantially increase the computation accuracy for given computer resources.

Equations

$$\begin{array}[]{c}\xi_{\alpha\beta\gamma}\partial_{\beta}E_{\gamma}\left(\bm{r}\right)=i\omega\mu_{0}H_{\alpha}\left(\bm{r}\right)-M_{\alpha}\left(\bm{r}\right),\\ \xi_{\alpha\beta\gamma}\partial_{\beta}H_{\gamma}\left(\bm{r}\right)=-i\omega\varepsilon\left(\bm{r}\right)E_{\alpha}\left(\bm{r}\right)+J_{\alpha}\left(\bm{r}\right)\end{array}$$

$$\begin{array}[]{l}\tilde{E}_{3}=E_{3}-J_{3}/i\omega\varepsilon_{b}\\ \tilde{H}_{3}=H_{3}-M_{3}/i\omega\mu_{b}\end{array}$$

$${\tilde{E}}_{\alpha }\left(\mathbf{r}\right)=\sum _{m=-\infty }^{\infty }\sum _{\sigma =\pm }\left[{\tilde{a}}_m^{e\sigma }\left(x_3\right){\hat{e}}_{m\alpha }^{e\sigma }+{\tilde{a}}_m^{h\sigma }\left(x_3\right){\hat{e}}_m^{h\sigma }\right]\exp \left(i{k}_{m1}{x}_1+i\sigma k_{m3}{x}_3\right),$$

$${\tilde{H}}_{\alpha }\left(\mathbf{r}\right)=\frac{k_b}{\omega {\mu }_0}\sum _{m=-\infty }^{\infty }\sum _{\sigma =\pm }\left[{\tilde{a}}_m^{h\sigma }\left(x_3\right){\hat{e}}_{m\alpha }^{e\sigma }-{\tilde{a}}_m^{e\sigma }\left(x_3\right){\hat{e}}_{m\alpha }^{h\sigma }\right]\exp \left(i{k}_{m1}{x}_1+i\sigma k_{m3}{x}_3\right),$$

$$\begin{array}[]{c}\hat{\bm{e}}_{m}^{e\pm}=\dfrac{\bm{k}_{m}^{\pm}\times\hat{\bm{e}}_{3}}{\left|\bm{k}_{m}^{\pm}\times\hat{\bm{e}}_{3}\right|},\\ \hat{\bm{e}}_{m}^{h\pm}=\dfrac{\left(\bm{k}_{m}^{\pm}\times\hat{\bm{e}}_{3}\right)\times\bm{k}_{m}^{\pm}}{\left|\left(\bm{k}_{m}^{\pm}\times\hat{\bm{e}}_{3}\right)\times\bm{k}_{m}^{\pm}\right|}.\end{array}$$

$$\left(\!\!\begin{array}[]{c}\tilde{a}_{m}^{e+}\left(x_{3}\right)\\ \tilde{a}_{m}^{h+}\left(x_{3}\right)\end{array}\!\!\right)=\left(\!\!\begin{array}[]{c}\tilde{a}_{m}^{ext,e+}\left(x_{3}\right)\\ \tilde{a}_{m}^{ext,h+}\left(x_{3}\right)\end{array}\!\!\right)-\intop_{-\infty}^{x_{3}}d\zeta\dfrac{\mathcal{G}_{m}\left(x_{3},\zeta\right)}{2k_{m3}}\left(\!\!\begin{array}[]{c}\omega\mu_{b}\hat{\bm{e}}_{m}^{e+}\cdot{\bf J}_{m}\left(\zeta\right)-k_{b}\hat{\bm{e}}_{m}^{h+}\cdot{\bf M}_{m}\left(\zeta\right)\\ \omega\mu_{b}\hat{\bm{e}}_{m}^{h+}\cdot{\bf J}_{m}\left(\zeta\right)+k_{b}\hat{\bm{e}}_{m}^{e+}\cdot{\bf M}_{m}\left(\zeta\right)\end{array}\!\!\right)$$

$$\left(\!\!\begin{array}[]{c}\tilde{a}_{m}^{e-}\left(x_{3}\right)\\ \tilde{a}_{m}^{h-}\left(x_{3}\right)\end{array}\!\!\right)=\left(\!\!\begin{array}[]{c}\tilde{a}_{m}^{ext,e-}\left(x_{3}\right)\\ \tilde{a}_{m}^{ext,h-}\left(x_{3}\right)\end{array}\!\!\right)-\intop_{x_{3}}^{\infty}d\zeta\dfrac{\mathcal{G}_{m}\left(x_{3},\zeta\right)}{2k_{m3}}\left(\!\!\begin{array}[]{c}\omega\mu_{b}\hat{\bm{e}}_{m}^{e-}\cdot{\bf J}_{m}\left(\zeta\right)-k_{b}\hat{\bm{e}}_{m}^{h-}\cdot{\bf M}_{m}\left(\zeta\right)\\ \omega\mu_{b}\hat{\bm{e}}_{m}^{h-}\cdot{\bf J}_{m}\left(\zeta\right)+k_{b}\hat{\bm{e}}_{m}^{e-}\cdot{\bf M}_{m}\left(\zeta\right)\end{array}\!\!\right)$$

$$\begin{array}[]{l}x_{1,2}=z^{1,2},\\ x_{3}=\begin{cases}z^{3}+\phi\left(z^{3}\right)f\left(z^{1}\right),&\left|x_{3}\right|\leq b\\ z^{3},&\left|x_{3}\right|>b\end{cases}\end{array}$$

$$\left\{g^{\alpha\beta}\right\}=\left(\!\!\!\begin{array}[]{ccc}1&0&-\dfrac{\phi f^{\prime}}{1+\phi^{\prime}f}\\ 0&1&0\\ -\dfrac{\phi f^{\prime}}{1+\phi^{\prime}f}&0&\dfrac{1+\phi^{2}f^{\prime 2}}{\left(1+\phi^{\prime}f\right)^{2}}\end{array}\!\!\!\right),$$

$$\begin{array}[]{c}\xi^{\alpha\beta\gamma}\partial_{\beta}E_{\gamma}\left(\bm{r}\right)=i\omega\mu_{0}\sqrt{g}g^{\alpha\beta}H_{\beta}\left(\bm{r}\right)-M^{\alpha}\left(\bm{r}\right)\\ \xi^{\alpha\beta\gamma}\partial_{\beta}H_{\gamma}\left(\bm{r}\right)=-i\omega\varepsilon\left(\bm{r}\right)\sqrt{g}g^{\alpha\beta}E_{\beta}\left(\bm{r}\right)+J^{\alpha}\left(\bm{r}\right)\end{array}$$

$$\begin{array}[]{c}\xi^{\alpha\beta\gamma}\partial_{\beta}E_{\gamma}\left(\bm{r}\right)=i\omega\mu_{b}\delta^{\alpha\beta}H_{\beta}\left(\bm{r}\right)-M_{gen}^{\alpha}\left(\bm{r}\right)-M^{\alpha}\left(\bm{r}\right),\\ \xi^{\alpha\beta\gamma}\partial_{\beta}H_{\gamma}\left(\bm{r}\right)=-i\omega\varepsilon_{b}\delta^{\alpha\beta}E_{\beta}\left(\bm{r}\right)+J_{gen}^{\alpha}\left(\bm{r}\right)+J^{\alpha}\left(\bm{r}\right)\end{array}$$

$$\begin{array}[]{c}J_{gen}^{\alpha}\left(\bm{r}\right)=-i\omega\varepsilon_{b}\left(\eta_{e}\sqrt{g}g^{\alpha\beta}-\delta^{\alpha\beta}\right)E_{\beta}\left(\bm{r}\right),\\ M_{gen}^{\alpha}\left(\bm{r}\right)=-i\omega\mu_{b}\left(\eta_{h}\sqrt{g}g^{\alpha\beta}-\delta^{\alpha\beta}\right)H_{\beta}\left(\bm{r}\right).\end{array}$$

$$\begin{array}[]{c}\partial_{3}E_{2}=M^{1}-i\omega\mu_{b}H_{1},\\ \partial_{1}E_{2}=-M^{3}+i\omega\mu_{b}H_{3},\\ \partial_{3}H_{1}-\partial_{1}H_{3}=-i\omega\varepsilon_{b}E_{2}+J^{2},\end{array}$$

$$\begin{array}[]{c}\partial_{1}\left[\sqrt{g}\left(H_{1}+g^{13}H_{3}\right)\right]+\partial_{3}\left[\sqrt{g}\left(g^{31}H_{1}+g^{33}H_{3}\right)\right]=0,\\ \partial_{2}\left(\sqrt{g}E_{2}\right)=0.\end{array}$$

$$\begin{array}[]{c}\partial_{3}\tilde{H}_{1}-\partial_{1}C_{2}=-i\omega\varepsilon_{b}\tilde{E}_{2}+J^{2},\\ \partial_{1}C_{1}+\partial_{3}\tilde{H}_{3}=0,\end{array}$$

$$\begin{array}[]{c}C_{1}=\dfrac{1}{\sqrt{g}g^{33}}\tilde{H}_{1}+\dfrac{g^{13}}{g^{33}}\dfrac{1}{\eta_{h}}\tilde{H}_{3}\\ C_{2}=-\dfrac{g^{31}}{g^{33}}\tilde{H}_{1}+\dfrac{1}{\sqrt{g}g^{33}}\dfrac{1}{\eta_{h}}\tilde{H}_{3}\end{array}$$

$$\begin{array}[]{c}\left(\phi f^{\prime}\right)C_{1}-C_{2}=-\left(1+\phi^{\prime}f\right)\dfrac{1}{\eta_{h}}\tilde{H}_{3}\\ C_{1}+\left(\phi f^{\prime}\right)C_{2}=\left(1+\phi^{\prime}f\right)\tilde{H}_{1}\end{array}$$

$$\left(\!\!\!\begin{array}[]{c}\left[C_{1}\right]\\ \left[C_{2}\right]\end{array}\!\!\!\right)=\left(\mathcal{I}+\mathcal{B}\mathcal{B}\right)^{-1}\left(\!\!\!\begin{array}[]{cc}\mathcal{I}&-\dfrac{1}{\eta_{h}}\mathcal{B}\\ \mathcal{B}&\dfrac{1}{\eta_{h}}\mathcal{I}\end{array}\!\!\!\right)\mathcal{A}\left(\!\!\!\begin{array}[]{c}\left[\tilde{H}_{1}\right]\\ \left[\tilde{H}_{3}\right]\end{array}\!\!\!\right),$$

$$\left(\!\!\!\begin{array}[]{c}\left[J^{2}\right]\\ \left[M^{1}\right]\\ \left[M^{3}\right]\end{array}\!\!\!\right)=-i\omega\epsilon_{b}\mathrm{M}^{-1}\mathrm{N}\left(\!\!\!\begin{array}[]{c}\left[\tilde{E}_{2}\right]\\ \left[\tilde{H}_{1}\right]\\ \left[\tilde{H}_{3}\right]\end{array}\!\!\!\right).$$

$$\mathrm{M}=\left(\!\!\!\begin{array}[]{ccc}\mathcal{I}&0&0\\ 0&\mathcal{I}+\mathcal{B}\mathcal{B}&0\\ 0&0&\mathcal{I}+\mathcal{B}\mathcal{B}\end{array}\!\!\!\right),$$

$$\mathrm{N}^{e}=\left(\!\!\!\begin{array}[]{ccc}\eta_{e}\mathcal{A}-\mathcal{I}&0&0\\ 0&\eta_{h}\mathcal{A}-\mathcal{I}-\mathcal{BB}&-\mathcal{BA}\\ 0&-\mathcal{BA}&\mathcal{I}+\mathcal{BB}-\dfrac{1}{\eta_{h}}\mathcal{A}\end{array}\!\!\!\right).$$

$${\mathbf{a}}^e={\mathbf{a}}^{e,ext}+\mathrm{G}{\mathrm{P}}^e{\mathrm{M}}^{-1}{\mathrm{N}}^e{ϵ}^e{\mathcal{F}}^e.$$

$${\mathcal{F}}_{mk}^e={\left({\tilde{E}}_{m2}(z_k^3),{\tilde{H}}_{m1}(z_k^3),{\tilde{H}}_{m3}(z_k^3)\right)}^T.$$

$$\mathrm{P}_{klmn}^{e}=\delta_{kl}\delta_{mn}\left(\!\!\!\begin{array}[]{ccc}\dfrac{\omega\mu_{b}}{k_{m3}}&-1&\dfrac{k_{m1}}{k_{m3}}\\ \dfrac{\omega\mu_{b}}{k_{m3}}&1&\dfrac{k_{m1}}{k_{m3}}\end{array}\!\!\!\right),$$

$$\mathrm{G}_{klmn}=-\delta_{mn}\frac{i\omega h_{s}}{2}\left(\!\!\!\begin{array}[]{cc}\exp\left(ik_{m3}\left|z_{k}^{3}-z_{l}^{3}\right|\right)&0\\ 0&\exp\left(ik_{m3}\left|z_{k}^{3}-z_{l}^{3}\right|\right)\end{array}\!\!\!\right).$$

$$\mathrm{Q}_{klmn}^{e}=\delta_{kl}\delta_{mn}\left(\begin{array}[]{cc}-1&-1\\ \dfrac{k_{m3}}{\omega\mu_{b}}&-\dfrac{k_{m3}}{\omega\mu_{b}}\\ -\dfrac{k_{m1}}{\omega\mu_{b}}&-\dfrac{k_{m1}}{\omega\mu_{b}}\end{array}\right),$$

$${\mathcal{F}}^e={\mathcal{F}}^{e,ext}+{\mathrm{Q}}^e\mathrm{G}{\mathrm{P}}^e{\mathrm{M}}^{-1}{\mathrm{N}}^e{ϵ}^e{\mathcal{F}}^e.$$

$${\mathcal{J}}^e={\mathrm{M}}^{-1}{\mathrm{N}}^e{ϵ}^e{\mathcal{F}}^e={\left(\mathrm{M}-{\mathrm{N}}^e{ϵ}^e{\mathrm{Q}}^e\mathrm{G}{\mathrm{P}}^e\right)}^{-1}{\mathrm{N}}^e{ϵ}^e{\mathcal{F}}^{e,ext}.$$

$${\mathbf{a}}_b^{e,out}={\mathbf{a}}_b^{e,ext}+{\mathrm{G}}_b{\mathrm{P}}^e{\mathcal{J}}^e,$$

$$\mathrm{G}_{b,kmn}=-\delta_{mn}\frac{i\omega h_{s}}{2}\left(\!\!\!\begin{array}[]{cc}\exp\left(ik_{m3}\left|b-z_{k}^{3}\right|\right)&0\\ 0&\exp\left(ik_{m3}\left|z_{k}^{3}+b\right|\right)\end{array}\!\!\!\right).$$

Full text

1 Introduction

Electromagnetic planar grating diffraction problem for resonant domain structures composed of dispersive materials can be solved by various numerical methods depending on particular structure properties [1, 2]. While the direct diffraction problem often is not an issue in practice, related inverse problems impose strict requirements on the mentioned methods, as a huge number of direct problem simulation runs is generally required within inversion procedures [3]. This makes the research directed towards an increase of simulation efficiency in terms of computing resource requirements be in demand in the field of engineering of high-efficient optical structures, like resonant gratings [4]. In particular, efficient optimization procedures are of great interest within the modern trend of metasurface optimization, which is aimed at pushing forward the field so as to substitute conventional diffractive optical components with high-index metastructures.

Fourier space methods are good candidates for the purpose of large scale optimization of periodic dielectric structures being versatile in terms of structure geometries and offering a reliable convergence error analysis [5, 6]. There were proposed several fast numerical schemes exhibiting an $O\left(N\log N\right)$ numerical complexity and $O\left(N\right)$ computational memory resort [7, 8, 9, 10, 11], and providing a large room for parallelization on vector processors. This asymptotic behaviour is due to a specific structure of the scattering operator, which factorizes into block-diagonal and block-Toeplitz components, and due to the possibility to multiply Towplitz matrices by vectors by means of the Fast Fourier Transform (FFT). Owing to the remarkable features, the named methods possess a great potential for applications in optimization procedures, though, some additional specific improvements are required to ensure convergence for complex diffractive optical elements and metasurfaces. This includes the scattering vector algorithm [12], which implements the divide-and-conquer strategy for thick structures, and application of curvilinear coordinate transformations, which enhance the numerical behaviour of the Fourier methods. These transformations are of two types. The first type implies in-plane stretching and contracting a grating, and was developed within the context of the Fourier Modal Method (FMM) [13, 14, 15] and the Modal Method (MM) [16, 17] to simplify the correct treatment of the boundary conditions [18, 19, 20] and improve numerical solution conditioning. A similar procedure can be incorporated into the fast methods, though this possibility is not considered in this work. The second type of transformations, which inevitably affects the coordinate orthogonal to the grating plane, includes transformations converting a complex corrugation interface to a plane. An idea of such transformations gave rise to the Chandezon method [21, 22, 23], which is a Fourier modal method in the curvilinear coordinate space [24]. In [25, 26] this idea was expanded to fit the rationale of computationally efficient methods by introducing the generalized metric sources. The corresponding method will be referred to as the Generalized Source Method in Curvilinear Coordinates (GSMCC), and its implementation paved a way to the efficient Fourier space simulation of metallic structures. However, the formulation given in [25, 26] is based on an assumption of the smoothness of functions describing periodic corrugations, thus, having a limited range of validity. In this work a generalization of the latter methods is developed in case of 1D gratings, which profile functions are allowed to have corners at a limited number of points, so that the profiles can be flattened with appropriate coordinate transformations. Note also, that the method described here represents a rigorous explicit numerical procedure of solving Maxwell’s equations with a controlled accuracy, and its validity is not constrained by convergence limitations being inherent to explicit and perturbative techniques like [27], or Rayleigh series based methods [28].

The article structure is the following. The narrative starts with the grating diffraction problem statement. Then the GSMCC rationale is outlined in form of a sequence of steps leading to a required solution form. Details on the corresponding derivations can be found in the previous papers [25, 26, 29]. A consequent analysis of effective metric sources and the Maxwell’s equations in curvilinear metric yields the derivation of correct handling of truncated Fourier vectors. Next, an efficient numerical method based on these rules is described. Numerical examples and the discussion of the results enclose the paper.

2 Grating diffraction

Let us consider a simple case of a periodic grating corrugation separating two different media, which is supposed to be described by a piece-wise smooth periodic function with bounded derivative $x_{3}=f\left(x_{1}\right)$, with period $\Lambda$, such that $f\left(x_{1}+n\Lambda\right)=f\left(x_{1}\right)$, $n\in\mathbb{N}$, in the Cartesian coordinates $x_{\alpha}$, $\alpha=1,2,3$. Without loss of generality one can suppose that $\min f\left(x_{3}\right)=-a$, $\max f\left(x_{3}\right)=a$. Dielectric permittivities of the media below and above the corrugation interface are denoted as $\varepsilon_{1}$ and $\varepsilon_{2}$ respectively. Fig. 1(a) illustrates the setup. Linear electromagnetic optical diffraction is governed by the time-harmonic Maxwell’s equations with implicit time dependence factor $\exp\left(-i\omega t\right)$

[TABLE]

subject to the continuity of the tangential fields at the corrugation interface, and the radiation condition at $x_{3}\rightarrow\pm\infty$ [2]. Here $J_{\alpha}$, $M_{\alpha}$ are the electric and magnetic source terms respectively. The Greek indices vary in range $1,2,3$, and $\bm{r}=\left(x_{1},x_{2},x_{3}\right)^{T}$. The Levi-Civita symbol $\xi_{\alpha\beta\gamma}=1$ for even permutations of $(1,2,3)$, and $\xi_{\alpha\beta\gamma}=-1$ for odd permutations. Summation over repeated coordinate indices is supposed here and further. Permittivity function is $\varepsilon\left(\bm{r}\right)=\varepsilon_{1}$ when $x_{3}\leq f\left(x_{1}\right)$, and $\varepsilon\left(\bm{r}\right)=\varepsilon_{2}$ when $x_{3}>f\left(x_{1}\right)$. A generalization may include a substrate, such that $\varepsilon\left(\bm{r}\right)=\varepsilon_{3}$ for $x_{3}<b$.

To set up a general solution of Eqs. (1), first, the dielectric permittivity and the magnetic permeability inside some region $D_{g}=\left\{\left|x_{3}\right|\leq b:b\geq a\right\}$ incorporating the grating are assumed to be constant everywhere $\varepsilon\left(\bm{r}\right)=\varepsilon_{b}$, $\mu=\mu_{b}$, with the subscript standing for “basis”. This assumption replaces the grating with a homogeneous plane layer. The electromagnetic field for this basis layer can be found for any sources by means of the volume integral equations with the electric and magnetic Green’s tensors $G_{\alpha\beta}^{E,M}\left(\bm{r},\bm{r}^{\prime}\right)$. Attributing the “difference” between the grating and the layer to effective (generalized) sources, one attains a self-consistent linear equation system for the searched solution of the diffraction problem. For the planar grating problem the Green’s tensors are subject to the plane wave decomposition (see, e.g., [30, 29]). In order to get rid of a singular term proportional to $\delta\left(x_{3}-x^{\prime}_{3}\right)$ present in the explicit form of $G_{\alpha\beta}^{E}\left(\bm{r},\bm{r}^{\prime}\right)$ it is convenient to introduce the modified fields as

[TABLE]

with the rest components being untouched: $\tilde{E}_{1,2}=E_{1,2}$, $\tilde{H}_{1,2}=H_{1,2}$. The periodicity allows one to fix a zero harmonic wave vector $\bm{k}^{inc}=\left(k_{1}^{inc},0,k_{3}^{inc}\right)$ and to represent the fields and sources in form of Bloch waves. This makes possible a decomposition of the modified electromagnetic field into a set of TE (superscript “$e$”) and TM (superscript “$h$”) polarized plane waves propagating upwards and downwards relative to the vertical coordinate $x_{3}$:

[TABLE]

[TABLE]

where the index $m$ enumerates a discrete set of plane waves having wavevectors $\bm{k}_{m}^{\pm}=\left(k_{m1},0,\pm k_{m3}\right)$ defined by the grating equation $k_{m1}=k_{1}^{inc}+2\pi m/\Lambda$, and the dispersion equation $k_{m1}^{2}+k_{m3}^{2}=k_{b}^{2}$, $\Re k_{m3}+\Im k_{m3}\geq 0$ with $k_{b}=\omega\sqrt{\varepsilon_{b}\mu_{0}}$. Unit vectors $\hat{\bm{e}}_{m}^{e\pm}$ and $\hat{\bm{e}}_{m}^{h\pm}$, which specify the polarization, explicitly write

[TABLE]

Then, the vertical coordinate dependent amplitudes occurring in Eqs. (3), (4) become solutions of the following integral equations (see [29] for a detailed derivation)

[TABLE]

[TABLE]

with the coefficients $\mathcal{G}_{m}\left(\xi,\eta\right)$ coming from the plane wave decomposition of the Green’s tensor. In general case of a planar waveguide they can be found, e.g., in [26, 31], and in particular, when $\varepsilon_{b}=\varepsilon_{1}=\varepsilon_{2}$ they explicitly read $\mathcal{G}_{m}\left(\xi,\eta\right)=\exp\left(ik_{3m}\left|\xi-\eta\right|\right)$ (note, that in this especial case the diffraction problem can still be nontrivial, once $\varepsilon(\bm{r})\neq\varepsilon_{b}$ in the grating region). Here it is supposed that the sources are split in two parts: the first part excites some known external fields $E_{\alpha}^{ext}$, $H_{\alpha}^{ext}$ with known plane wave decomposition amplitudes $\tilde{a}_{m}^{ext,e,h\pm}\left(x_{3}\right)$, while the second part, $J_{\alpha}$, $M_{\alpha}$, present in the integrands corresponds to some local sources to be specified below, which bring an information about the mentioned difference between the grating and the plane homogeneous layer in the region $D_{g}$.

3 Generalized metric sources

Once the general solution equations (6)-(7) are established they can be enclosed on the basis of the curvilinear coordinate transformation idea mentioned in the introduction. Consider curvilinear coordinates $\left(z^{1},z^{2},z^{3}\right)$ such that in $D_{g}=\left\{\left|x_{3}\right|\leq b:b\geq a\right\}$ the coordinate plane $z^{3}=0$ coincides with the corrugation profile, and the coordinates continuously become Cartesian at boundaries $\partial D=\left\{x_{3}=\pm b\right\}$. Such transformation can be defined as

[TABLE]

where $\phi\left(z^{3}\right)=1-\left|z^{3}\right|/b$. An illustrative example is given in Fig. 1(b) . Possible generalizations of the definition (8) concern multilayer structures, and are discussed in [26] with the difference that the corrugation function $f\left(x_{1}\right)$ is not required to be smooth here. Transformation (8) yields the metric tensor

[TABLE]

which components are discontinuous functions of coordinates due to the discontinuity of $f^{\prime}$.

Being written in the curvilinear coordinates the Maxwell’s equations become [32]

[TABLE]

where the lower and upper indices distinguish the covariant and contravariant vector components respectively, and $g$ denotes the determinant of a matrix being inverse to (9): $g=\det\left\{g_{\alpha\beta}\right\}=1/\det\left\{g^{\alpha\beta}\right\}$. An observation that the metric tensor explicitly occurs only behind the field terms in the right-hand parts of Eq. (10) brings the second cornerstone idea of the GSMCC [25]: Eqs. (10) can be rewritten in a form similar to the Maxwell’s equations in the Cartesian coordinates (1) and the rest can be attributed to generalized electromagnetic sources originating from the difference between the curvilinear and the Cartesian formulations. Namely, Eqs. (10) become

[TABLE]

with the generalized metric sources

[TABLE]

Here $\eta_{e}=\varepsilon\left(z^{3}\right)/\varepsilon_{b}$, $\eta_{h}=\mu_{0}/\mu_{b}$. Due to the similarity between the operator parts of Eqs. (1) and Eqs. (11) such decomposition allows directly using solutions (6)-(7) with the sources (12) upon a mere substitution of the coordinates $x_{\alpha}$ with $z^{\alpha}$ and recalling that the modified fields should replace the real ones in Eq. (12) according to Eq. (2). The function $\varepsilon\left(\bm{r}\right)$ in the latter equations depends only on the coordinate $z^{3}$ and takes two constant values below and above the flattened corrugation interface: $\varepsilon\left(z^{3}\right)=\varepsilon_{1}$ for $z^{3}\leq 0$, and $\varepsilon\left(\bm{r}\right)=\varepsilon_{2}$ for $z^{3}>0$, as Eq. (8) states.

Eq. (12) being substituted into (6)-(7) yields a self-consistent system, taking into account Eqs. (2)-(4). At this point in case of a smooth coordinate transformation it is mathematically justified to perform the truncation of infinite Fourier series of the generalized sources [18, 33], so that Eqs. (6)-(7) become a finite linear integral equation system [25]. In case of non-smooth transformations considered here a further effort is required, as the next section explains.

In analogy with the Generalized Source Method [7, 9] and other volume integral methods in the Fourier space, which enclose the volume integral equation with sources of the form ${\bf J}=-i\omega\left[\varepsilon\left(\bm{r}\right)-\varepsilon_{b}\right]{\bf E}$, Eq. (12) demonstrates that the curvilinear metric impact can be treated as inhomogeneous and anisotropic material tensors $\hat{\varepsilon}^{\alpha\beta}=\varepsilon\left(z^{3}\right)\sqrt{g}g^{\alpha\beta}$, and $\hat{\mu}^{\alpha\beta}=\mu_{0}\sqrt{g}g^{\alpha\beta}$. Due to the discontinuity of the metric tensor along the coordinate $z^{1}$ the points of discontinuity appear to effectively act like vertical material interfaces, as Fig. 1(c) illustrates. Thus, special precautions should be taken when calculating Fourier coefficients of the generalized metric sources to be used in Eqs. (6)-(7), in analogy with the Fourier methods in the Cartesian space [18]. In case of the C-method the factorization rules were thoroughly explained in [33] in terms of relations between covariant field components. The same rationale can be applied to the GSMCC, though, here the correct Fourier space matrix-vector relations are demonstrated to appear in a way analogous to the derivation of the electromagnetic interface conditions, as the following section demonstrates.

4 Fourier factorization of discontinuous metric sources

The work [18] demonstrated that justifiability of truncation of infinite Fourier series of function products strongly depends on the continuity of the involved functions. Therefore, in case of discontinuous along $z^{1}$ metric tensor and field components, finite matrix relations between the Fourier components of the generalized sources and the fields should be properly derived. In order to demonstrate the correct truncated Fourier factorization consider the Maxwell’s equations for the TE polarization (due to the presence of the both electric and magnetic sources, derivations for the TM polarization are quite similar):

[TABLE]

where the source terms represent a superposition of real sources exciting incoming diffracting waves, and the local generalized metric sources. In addition to Eqs. (13) let us explicitly recall the Gauss law

[TABLE]

Here the dependency of $\varepsilon\left(\bm{r}\right)$ from the coordinate $z^{3}$ only and $g^{11}=1$ property of Eq. (9) were taken into account. The existence of the right-hand parts of Eqs. (13), (14) requires the existence of corresponding derivatives. Recalling the transition to the modified fields, Eq. (2), and Eqs. (12), the last of Eq. (13), and the first of Eq. (14) become

[TABLE]

where the non-trivial combinations of the field and the metric tensor components

[TABLE]

should be continuous along $z^{1}$ coordinate, and, hence, across the effective vertical interfaces discussed at the end of the previous section. The continuity property follows directly from the same derivations as the ones widely used in university textbooks to attain the electromagnetic interface conditions by integrating Maxwell’s equations in the vicinity of an interface, (see, e.g., [34], Ch. 1.5.). Substitution of the explicit $g^{\alpha\beta}$ components into Eq. (16) yields

[TABLE]

Due to the continuity of $C_{1,2}$ and $(1+\phi^{\prime}f)$ along $z^{1}$ one can directly take the Fourier transform of the latter equations [18] and pass to truncated series to derive relations between the Fourier coefficients $C_{1,2m}$ and $\tilde{H}_{1,3m}$. Denoting the Fourier component vectors with square brackets as $\left[\bullet\right]$, Fourier-Toepltz matrices – with double square brackets as \left\llbracket\bullet\right\rrbracket; expressing explicitly $\left[C_{1,2}\right]$ via $\left[\tilde{H}_{1,3}\right]$:

[TABLE]

with \mathcal{A}=\left\llbracket 1+\phi^{\prime}f\right\rrbracket, \mathcal{B}=\left\llbracket\phi f^{\prime}\right\rrbracket, and identity matrix $\mathcal{I}$; and substituting the modified fields into Eq. (12), one attains the following relations between the Fourier vectors of the generalized metric sources and the fields:

[TABLE]

Here

[TABLE]

[TABLE]

and diagonal matrix $\epsilon_{b}=diag\{\varepsilon_{b},\,\mu_{b},\,\mu_{b}\}$. The terms $\mathcal{BB}$, and $\mathcal{BA}$ imply that vectors should be successively multiplied by the corresponding matrices. Eqs. (18)-(21) hold for each fixed value of $z^{3}$.

5 Numerical method

Definition of the coordinate transformation, Eq. (8), implies that the generalized metric sources (12) are non-zero only when $\left|z^{3}\right|\leq b$. Therefore, the integration limits in Eqs. (6)-(7) do not fall outside the region $D_{g}$. Let us introduce an equidistant mesh (slicing) in $D_{g}$ defined by coordinates $z_{k}^{3}=\left(2k+1-N_{s}\right)b/N_{s}$, $k=0,\dots,N_{s}-1$, with slice thickness $h_{s}=2b/N_{s}$, and evaluate the integrals at the mesh points using the mid-point rule. Also denote a maximum order of truncated infinite Fourier vectors and matrices as $N_{F}$.

Upon substitution of the generalized metric sources, Eqs. (6)-(7) reduce to a set of linear equations on the unknown vector of the TE wave amplitudes in each slice $\bm{a}_{mk}^{e}=\left(\tilde{a}_{m}^{e+}(z_{k}^{3}),\,\tilde{a}_{m}^{e-}(z_{k}^{3})\right)^{T}$:

[TABLE]

Here $\bm{a}^{e,ext}$ is the known amplitude vector of waves coming from the exterior of the layer $D_{g}$, and the unknown field vector $\bm{\mathcal{F}}^{e}$ in case of the TE polarization is

[TABLE]

The Fourier index $m$ runs in range $-N_{F}\leq m\leq N_{F}$ and the spatial mesh index $k$ – in range $0\leq k<N_{s}$. The block-diagonal matrix $\mathrm{P}^{e}$ defines the source-to-amplitude in-slice transformation of Eqs. (6)-(7)

[TABLE]

The matrix $\mathrm{G}$ is defined by operator $\mathcal{G}$ together with factor $-i\omega h_{s}/2$. In case $\varepsilon_{b}=\varepsilon_{1}=\varepsilon_{2}$

[TABLE]

Generally, when the two permittivities $\varepsilon_{1,2}$ present in the regions $x_{3}<-b$ and $x_{3}>b$ are different from $\varepsilon_{b}$, the Green’s tensor becomes more complex, and the latter matrix element should be replaced with another one incorporating multiple reflections at interfaces $x_{3}=\pm b$. For explicit equations in this case see, for example, [26, 31].

Multiplication of Eq. (22) by

[TABLE]

which is composed of components of unit vectors (5), yields the desired linear equation system

[TABLE]

The derivation of Eq. (27) is analogous to what is done in [25, 26], but here one faces with an additional matrix inversion coming from the field-to-source transformation of Eq. (19). Hereof it follows, that when solving this system by an iterative method, each iteration would require performing a numerically expensive operation of matrix inversion. This, in turn, would result in total $O(N_{s}N_{F}^{3})$ numerical complexity instead of the GSMCC breakthrough $O(N_{s}N_{F}\log(N_{s}N_{F}))$ [25], which is based on a decomposition of the corresponding linear system matrix into block-diagonal and block-Toeplitz constituents, and application of the Fast Fourier Transform for speeding-up multiplications. To handle this inversion let us multiply Eq. (27) by $\mathrm{M}^{-1}\mathrm{N}^{e}\epsilon^{e}$ and introduce the new vector

[TABLE]

This system is solved for the unknown vector $\bm{\bm{\mathcal{J}}}^{e}$, it can be substituted into Eq. (27) to attain a vector of diffracted amplitudes at the grating region boundaries:

[TABLE]

Here the output and external amplitude vectors are $\bm{a}_{b,m}^{e,out/ext}=\left(\tilde{a}_{m}^{out/ext,e+}(b),\,\tilde{a}_{m}^{out/ext,e-}(-b)\right)^{T}$, and the matrix operator $\mathrm{G}_{b}$ also comes from tensor $\mathcal{G}$, and “collects” the waves diffracted in each slice at $D_{g}$ boundaries. In case $\varepsilon_{b}=\varepsilon_{1}=\varepsilon_{2}$

[TABLE]

Since the curvilinear metric continuously transforms to the Cartesian one at $\left|z^{3}\right|=b$, and the modified field coincides with the real one outside the grating region, the output amplitudes immediately define the diffracted field outside $D_{g}$, and no other transformation is required. For the TM polarization resulting equations are the same as Eqs. (28), (29), but with slightly different matrices $\mathrm{N}$, $\epsilon$, $\mathrm{P}$, and $\mathrm{Q}$. They are listed in Appendix A.

The provided formulation of Eqs. (28), (29) theoretically preserves the fast and memory efficient conception of the previous papers [25, 26], and allows one to perform the grating diffraction calculation in $O\left(N\log N\right)$ time and $O\left(N\right)$ memory resort with $N=N_{s}N_{F}$. Namely, this comes from the fact that multiplications by block-diagonal matrices $\mathrm{P}^{e}$, $\mathrm{Q}^{e}$, $\mathrm{G}$, and $\mathrm{G}_{b}$ are of linear asymptotic complexity, and multiplications by $\mathrm{M}$, and $\mathrm{N}^{e}$ are of $O\left(N\log N\right)$ asymptotic complexity due to the block-Toeplitz structure.

In order to demonstrate the validity of the developed method and codes as well as to study the asymptotic rate of convergence in dependence of $N_{F}$ let us consider an example of a non-symmetric triangular corrugation profile with the period-to-wavelength ratio $\Lambda/\lambda=1.5$ and $d_{1}/d_{2}=2/3$ (see Fig. 1(a)). Matrices (20), (21) are found analytically in case of an arbitrary piecewise-linear finction $f$. For convergence calculation there was chosen the starting value $N_{F,0}=16$, and a set of method runs for subsequent values of $N_{F}$ was performed. With the increase of $N_{F}$ the size of output amplitude vectors $\bm{a}^{out}$ also increases. Therefore, only central amplitudes with indices $-N_{F,0}\leq m\leq N_{F,0}$ were picked up for comparison. The convergence was plotted as the dependence of the maximum absolute value $\max|a^{out}_{m}(N_{F,k+1})-a^{out}_{m}(N_{F,k})|$, $N_{F,k+1}=N_{F,k}+1$, from the inverse of $N_{F,k}$. Since the validity of the former GSMCC formulation was established previously [25], and generally incorrect and correct formulations of the Fourier methods are known to converge to the same results [6, 18], no other reference method for calculation of diffraction amplitudes is used here. Other polygonal grating shapes can be simulated by means of a sample Matlab code downloadable from [35].

Examples of the method convergence in dependence of $N_{F}$ are given in Figs. 2,3 for the mentioned triangular corrugation profile analysed with the formulation of [25], which is incorrect for grating with corners, and with the correct formulation provided here. Grating depth-to-wavelength ratio, and slice thickness were taken to be $2a/\lambda=1.0$, $h_{s}=0.001\lambda$ for a dielectric grating with $\varepsilon_{1}=1$, $\varepsilon_{2}=\left(2.5\right)^{2}$, and $2a/\lambda=0.5$, $h_{s}=0.0005\lambda$ for a metallic grating with the same $\varepsilon_{1}$ and $\varepsilon_{2}=\left(0.2+3.2i\right)^{2}$. The angle of incidence was $10^{\circ}$. The upper subplots in each figure show changes in absolute values of zero order diffraction efficiencies, and the lower subplots demonstrate the convergence defined above. The trendlines $x^{2}$ and $x^{3}$, reveal, that the new formulation not only allows obtaining substantially more accurate results for a given $N_{F}$ than the old one, but also that it exhibits a superior rate of convergence – close to cubic but in the case of the metallic grating for the TM polarization, which might be related to a poor Fourier space field representation in the vicinity of metallic corners. It can be also noticed that the convergence of the incorrect formulation in case of the TM polarization is rather slow.

Despite the theoretical complexity of the method is close to linear, in practice it appears to get worse due to an impact of a linear equation iterative solver. The presented examples were calculated by means of the Generalized Minimal Residual method (GMRes). Fig. 4 shows the dependence of the method run time from $N_{F}$ for the four considered cases of the dielectric and metal gratings under the TE and TM polarizations. The run time appears to grow linearly for small $N_{F}$, and faster for larger numbers of the Fourier harmonics, which is due to the increase of a number of iterations required for the GMRes to converge down to a given tolerance (taken to be $10^{-7}$ here). In case of the dielectric grating the iteration number grew from about 100 to 350. In case of metallic gratings it was found that solution of Eq. (28) for the whole grating exhibits quite slow GMRes convergence, or even stagnation. To solve this issue, the divide-and-conquer strategy described in Section 6 of [26] was applied (see also [12]). A particular case of this strategy implemented in this work is outlined in Appendix B.

For practical applications it is also important to track a dependence of the run time from grating depth-to-period ratio. Such dependence is shown in the right-hand part of Fig. 4 for dielectric gratings of the same geometry as considered above. For this set of simulations with increasing depth the vertical spatial resolution remained constant, $h_{s}=10^{-3}\lambda$, hence the number of slices $N_{s}$ increased linearly with the ratio $2a/\Lambda$. The dependence appears to be nonlinear as the number of GMRes iterations, shown in the right vertical axis, also increased. A corresponding time-depth dependence in case of metal gratings is quite similar.

6 Conclusion

To summarize this work extends the applicability of the GSMCC method by providing an explicit formulation of the method in case of gratings having corrugations profiles with corners. Eqs. (28), (29) reveal that the previously demonstrated low asymptotic numerical complexity and memory consumption can be preserved in the new formulation, and the numerical examples show its importance in attaining accurate diffraction efficiencies. The presented theory can be also applied in the case of conical diffraction by simply extending the involved matrices [7]. Calculation of diffraction by crossed gratings would require a separate revision of matrix-vector products, and will be reported elsewhere.

Concerning "overhanging" gratings, which profiles cannot be defined by single-valued functions, a parametric transformation can be used leading to a metric tensor being very similar to Eq. (9). The rest formulation would remain the same. An evidence for the GSMCC to be capable to handle such cases was reported in [36], where a curvilinear Fourier modal method was demonstrated to work for "overhanging" structures. This means that, in particular, the method should be valid for vertical wall structures. The main issue is either to find an appropriate transformation allowing to analytically calculate necessary Fourier matrices, or to develop a suitable method for numerical calculation of these matrices. In any case such demonstration lies outside the scope of this paper, and might be a subject of a future publication.

Despite the demonstrated practical calculation time growth with the increase of $N_{F}$ is faster than linear, it is still less than the cubic dependence inherent to the modal approaches. Absolute values of run times are strongly dependent on a programming implementation and hardware, and can be substantially decreased by using the graphical processing units. Therefore, the attained results can be further utilized in the field of grating structure optimization, where computationally efficient and reliable rigorous solvers of the Maxwells equations are of great importance.

Funding

Russian Science Foundation (RSF) (17-79-20345).

Appendix A

In case of the TM polarization solution of the diffraction problem is still given by Eqs. (28), (29), but with slightly different matrices $\mathrm{N}^{h}$, $\epsilon^{h}$, $\mathrm{P}^{h}$, and $\mathrm{Q}^{h}$. Namely, $\mathrm{N}^{h}$ is the same as Eq. (21), but with $\eta_{e}$ and $\eta_{h}$ interchanged. The other matrices are

[TABLE]

[TABLE]

and

[TABLE]

Appendix B

As noted in Section 5, application of Eq. (28) to a whole metallic grating results in a poor convergence or even stagnation of the GMRes depending of grating depth. A remedy to this issue is the use of a divide-and-conquer strategy described in Section 6 of [26]. To briefly outline it, denote the diffraction operator explicitly given by Eqs. (27), (29) as $\mathcal{S}$:

[TABLE]

The dependence on the lower and the upper bounds of the grating layer $\pm b$ emphasizes the fact that the operator is applied to the whole grating layer. When the layer is too thick for the linear iterative solver to converge in reasonable time, this layer can be divided in two as Fig. 1d illustrates. Let us introduce intermediate self-consistent amplitude vectors between the layers $\bm{a}^{\pm}(z^{3}=0)$. Then, defining the diffraction operators of the lower and the upper half-layers as $\mathcal{S}(-b,0)$, and $\mathcal{S}(0,b)$ allows to write out equations similar to Eq. (34) for the halves including the intermediate amplitude vectors. Rearrangement of blocks of these equations yields

[TABLE]

This equation system can be in turn solved by the GMRes. Therefore, solutions of the linear equation systems become nested: while solving Eq. (35), each iteration requires two solutions of Eqs. (27), (29) for each half-layer. Once vectors $\bm{a}^{\pm}(z^{3}=0)$ are calculated, the required output amplitudes come from

[TABLE]

This approach can be generalized to any number of sub-layers of the initial grating layer. Such implicit self-consistent method allows to dramatically decrease calculation time for metal gratings providing that all plane material interfaces in the curvilinear coordinates coincide with some sub-layer boundaries, and constant $\varepsilon_{b}$ is chosen to be different for each sub-layer, being equal to an averaged sub-layer permittivity.