# Domain Decomposition Spectral Method Applied to Modal Method: Direct and Inverse Spectral Transforms

Kofi Edee, G. Granet, Francoise Paladian, Pierre Bonnet, Ghida Al Achkar, Lana Damaj, Jean-Pierre Plumey, Maria Cristina Larciprete, Brahim Guizal

► To cite this version:

Kofi Edee, G. Granet, Francoise Paladian, Pierre Bonnet, Ghida Al Achkar, et al., Domain Decomposition Spectral Method Applied to Modal Method: Direct and Inverse Spectral Transforms. Sensors, 2022, 22 (8131), pp.8131. 10.3390/s22218131. [https://hal.science/hal-03826961](https://hal.science/hal-03826961)

HAL Id: hal-03826961

https://hal.science/hal-03826961

Submitted on 29 Nov 2022

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>
Article

# Domain Decomposition Spectral Method Applied to Modal Method: Direct and Inverse Spectral Transforms

Kofi Edee 1,2,* , Gérard Granet 1, Francoise Paladian 1, Pierre Bonnet 1, Ghida Al Achkar 1, Lana Damaj 1, Jean-Pierre Plumey 1, Maria Cristina Larciprete 2,3 and Brahim Guizal 2,4

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

check for
updates

**Citation:** Edee, K.; Granet, G.; Paladian, F.; Bonnet, P.; Al Achkar, G.; Damaj, L.; Plumey, J.-P.; Larciprete, M.C.; Guizal, B. Domain Decomposition Spectral Method Applied to Modal Method: Direct and Inverse Spectral Transforms. *Sensors* **2022**, *22*, 8131. https://doi.org/10.3390/s22218131

Academic Editor: Jean-Baptiste Thomas

Received: 27 September 2022

Accepted: 20 October 2022

Published: 24 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2022 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** We introduce a Domain Decomposition Spectral Method (DDSM) as a solution for Maxwell's equations in the frequency domain. It will be illustrated in the framework of the Aperiodic Fourier Modal Method (AFMM). This method may be applied to compute the electromagnetic field diffracted by a large-scale surface under any kind of incident excitation. In the proposed approach, a large-size surface is decomposed into square sub-cells, and a projector, linking the set of eigenvectors of the large-scale problem to those of the small-size sub-cells, is defined. This projector allows one to associate univocally the spectrum of any electromagnetic field of a problem stated on the large-size domain with its footprint on the small-scale problem eigenfunctions. This approach is suitable for parallel computing, since the spectrum of the electromagnetic field is computed on each sub-cell independently from the others. In order to demonstrate the method's ability, to simulate both near and far fields of a full three-dimensional (3D) structure, we apply it to design large area diffractive metalenses with a conventional personal computer.

**Keywords:** metasurfaces; metalens

## 1. Introduction

Regardless of the numerical method used to solve an electromagnetic (EM) problem, the number of unknowns of the field components increases as well as the electrical size of the device. Consequently, the modeling of extended structures can be time and memory consuming for a conventional computer. This makes the numerical simulation of large-scale problems more challenging for numerical simulations, and it is still an open task. Of course, in practice, the meaning of fast or large depends on the performance of the employed computer resources.

In order to overcome this challenging task and to ensure fast and reliable numerical solutions, two main points could be considered. The first one is the way the code is implemented, and the second one is focused on the numerical method used to solve equations at hand.

For a given electromagnetic problem stated on a large-scale domain, one of the first intuitive methods to overcome the large size issues is the sub-structuring of its computational domain. This technique consists of dividing the whole domain into a sequence of sub-domains. The solution inside these sub-domains may be found using an iterative or parallel strategy. Concerning the iterative process, the solution of the equation is found after solving the boundary conditions equations. This approach is considered within the context of the iterative domain decomposition techniques such as Schwarz's method [1-6]. However, when the number of interfaces increases, the strategy can be time consuming. On the other hand, the parallel computing strategy relies on overlapping [7] or stitching [8] the

<footnote n="1">1 INP, CNRS, Institut Pascal, Université Clermont Auvergne, F-63000 Clermont-Ferrand, France</footnote>

<footnote n="2">2 Kavli Institute of Theoretical Physics (KITP), University of California, Santa Barbara, CA 93106, USA</footnote>

<footnote n="3">3 Dipartimento di Scienze di Base ed Applicate per l'Ingegneria, Sapienza Universita di Roma, Via A. Scarpa 16, I-00161 Roma, Italy</footnote>

<footnote n="4">4 Laboratoire Charles Coulomb (L2C), UMR 5221 CNRS-Université de Montpellier, F-34095 Montpellier, France</footnote>

<footnote>* Correspondence: kofi.edee@uca.fr</footnote>
grid of the sub-domains without applying any boundary conditions. Concerning the numerical method, a large class of electromagnetic field simulation for photonics applications requires solving the partial differential equations (PDE) obtained from Maxwell's equations and a set of ad hoc boundary conditions. A priori, there is no unique and universal method to efficiently solve the PDEs obtained from Maxwell's equations in the general case. Among these methods, some clearly appear to be more used than others, either for their simplicity of implementation, or for their ability to take into account complex-geometry problems, or because they are the most efficient in terms of numerical accuracy. Yet, finite difference methods [9,10], finite element methods [11,12], or spectral modal methods [13–21] are considered from these methods. Furthermore, all these numerical methods involve different amounts of theoretical effort which yield different computational efficiencies. It is clear that improving the efficiency of the numerical implementation and process will often require more pre-computational analytical effort as well as more in-depth understanding and usage of electromagnetic theory. For instance, spectral modal methods aim to approximate an unknown function, as a solution of a differential or an integral equation, by a finite sum of eigenmodes. These eigenmodes are expanded into a so-called basis function. In the case of finite difference and finite elements methods, many studies have been conducted to make them suitable to large-scale problems [7,22,23].

Contrary to these methods, the spectral modal methods, in spite of many previous developments [7,8], still face some major challenges to find numerical solutions for some large-scale problems.

In this paper, we propose a mathematical formalism based on a projector enabling the projection of any eigenvector defined on a domain $\Omega$ into a sequence of eigenfunctions based on smaller-scale boundary value problems defined on subdomains $\tilde{\Omega}_{ij}$. In the proposed domain decomposition method, a large-scale surface is squared into sub-cells. Hence, a projector, i.e., a link between the large problem and the small one, is introduced. This projector, which is defined from the sequence of eigenfunctions of both large and small problems, univocally enables the association of the spectrum of any electromagnetic field component related to a large-size problem, with its restriction on the small-size systems. Consequently, the spectrum of the electromagnetic field component on each sub-cell may be simulated independently. This makes the proposed approach suitable to be integrated in a parallel computing scheme in order to allow for accounting larger size systems.

This paper is organized as follows: After the introduction (Section 1), we introduce in Section 2 the physical system and the proposed domain decomposition spectral method (DDSM). The DDSM is explained and described in detail, and the concept of direct and inverse transform is presented. Moreover, Section 3 is devoted to applications in the framework of the modal method. In this section, we introduce the bases of the modal method applied to a 3D problem. In a first illustration, we demonstrate its ability to accurately describe a field radiated by an electric dipole. In a second example, we prove the capability of the DDSM to simulate the near and far fields of a high-refractive-index-deeply-etched Fresnel plate zone, metalenses made of dielectric cylindrical waveguides (nanorods) with variables cross-section and a dielectric Pachatarnam-Berry phase metalens.

## 2. Methods

The framework of the proposed domain decomposition spectral method (DDSM) is the spectral modal method. For the reader that is not familiar with this method, we highlight in the following first two subsections the main underlying principle both of the modal and the spectral methods.

### 2.1. Maxwell's Equations and Eigenvalue Problem

Let us consider Maxwell's equations [20] in a Cartesian coordinates system ($O, e_x, e_y, e_z$) and the constitutive relations stated on a 2D bounded domain $\Omega = [-0.5D_x, 0.5D_x] \times [-0.5D_y, 0.5D_y] \subset \mathbb{R}^2$. We assume that the physical parameters of the medium, namely the relative permittivity $\epsilon(x,y)$ and the permeability $\mu(x,y)$ functions, do not depend on
the z variable. In this case, any component of the electromagnetic field can be written as $\Phi(x, y)\exp(-ik_0\gamma z)$, where $k_0 = \omega\sqrt{\epsilon_0\mu_0}$ is the wavenumber in the vacuum, $\omega$ is the angular frequency and $k_0\gamma$ is the propagation constant along the z direction. $\epsilon_0$ and $\mu_0$ are the permittivity and the permeability of the vacuum. It can be shown that the transverse components of the electromagnetic field $E_x, E_y, H_x$ and $H_y$ satisfy an eigenvalue equation:

$$ \mathcal{L}_{EH}\mathcal{L}_{HE}|\Phi\rangle = -\gamma^2|\Phi\rangle \quad (1) $$

where $\Phi = [E_x, E_y]^t$, and

$$ \left\{ \begin{array}{l} \mathcal{L}_{EH} = \frac{1}{k_0^2} \begin{bmatrix} -\partial_x \mathcal{E}^{-1} \partial_y & \partial_x \mathcal{E}^{-1} \partial_x + k_0^2 \mu \\ -\partial_y \mathcal{E}^{-1} \partial_y - k_0^2 \mu & \partial_y \partial_x \end{bmatrix} \\ \\ \mathcal{L}_{HE} = \frac{1}{k_0^2} \begin{bmatrix} -\partial_x \mu^{-1} \partial_y & \partial_x \mu^{-1} \partial_x + k_0^2 \mathcal{E} \\ -\partial_y \mu^{-1} \partial_y - k_0^2 \mathcal{E} & \partial_y \partial_x \end{bmatrix} \end{array} \right. , \quad (2) $$

The eigenvalue Equation (1) is solved using the spectral method briefly described in the next subsection.

## 2.2. Spectral Method and Modal Method

In electromagnetism, a spectral method aims to describe a function $|\Phi\rangle$ representing any electromagnetic field component with a set of finite expansion coefficients $\Phi^{nm}$ on a selected basis also known as expansion functions $|e_{nm}\rangle$:

$$ |\Phi\rangle = \Phi^{nm}|e_{nm}\rangle \quad (3) $$

by using Einstein's convention. Expansion or weighted coefficients $\Phi^{nm}$ are the spectrum of the expanded function $\Phi$ with respect to the sequence of expansion vectors $|e_{nm}\rangle$. In a Cartesian coordinates system with axis $(O, x), (O, y), (O, z)$, and in the framework of the modal method, the structure under consideration is often divided into some layers denoted $I_z^{(1)}$ with respect to the propagation direction that is assumed to be $(O, z)$ in our case. For example, in Figure 1, a 2D-diffraction grating is divided into 4 layers $I_z^{(1)}, I_z^{(2)}, I_z^{(3)}$ and $I_z^{(4)}$ with respect to the propagation direction.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Sketch of a 2D grating. In the framework of the modal method, the 2D structure is divided into 4 layers $I_z^{(1)}, I_z^{(2)}, I_z^{(3)}$ and $I_z^{(4)}$ with respect to the propagation direction that is assumed to be $(O, z)$.</figcaption></figure>

In any of the layers $I_z^{(l)}$, the physical parameters, namely the permittivity and permeability, are described as functions that do not depend on the variable z. Consequently, each component of the electromagnetic field may be expanded on a set of eigenfunctions, i.e., solutions $|\Phi_{kl}\rangle$ of an eigenvalue equation similar to Equation (1):

$$ \mathcal{L}|\Phi\rangle = -\gamma^2|\Phi\rangle. \quad (4) $$
The completeness of the set of eigenvectors (|Φ<sub>kl</sub>⟩<sub>kl</sub> allows one to expand any function |f⟩ defined on a domain Ω as a linear combination of the eigenvectors |Φ<sub>kl</sub>⟩:

$$|f\rangle = A^{kl}|\Phi_{kl}\rangle. \qquad (5)$$

Regarding Equation (5), the modal method can be viewed as a spectral method using as basis functions a set of eigenvectors |Φ<sub>kl</sub>⟩ of the operator L. The Fourier modal method (FMM) consists of expanding the eigenvectors |Φ<sub>kl</sub>⟩ in terms of a generalized Fourier basis:

$$|\Phi_{kl}\rangle = \Phi_{kl}^{nm} |e_{nm}\rangle, \qquad (6)$$

where |e<sub>nm</sub>⟩ = |e<sub>n</sub>⟩ ⊗ |e<sub>m</sub>⟩, with

$$e_{nm}(x, y) = e^{-ik_0\alpha_n x} e^{-ik_0\beta_m y}. \qquad (7)$$

Taking into account Equation (6), the vector |f⟩ of Equation (5) may be expressed in terms of the Fourier basis functions as:

$$|f\rangle = A^{kl} \Phi_{kl}^{nm} |e_{nm}\rangle. \qquad (8)$$

From a numerical point of view, only a finite number of (2M<sub>x</sub> + 1) × (2M<sub>y</sub> + 1), ((M<sub>x</sub>, M<sub>y</sub>) ∈ N<sup>2</sup>) Fourier coefficients is kept by selecting a finite sequence of α<sub>n</sub> and β<sub>n</sub> in the expression Equation (7):

$$\left\{ \begin{array}{l} \alpha_n = \alpha_0 + n\lambda/D_x, n \in [-M_x : M_x] \\ \beta_m = \beta_0 + m\lambda/D_y, m \in [-M_y : M_y] \end{array} \right. \qquad (9)$$

where D<sub>x</sub> (resp. D<sub>y</sub>) is the period along the (O, x) (resp. (O, y)) axis and λ is the operating wavelength. The parameters α<sub>0</sub> and β<sub>0</sub> are linked to the incident wave. The choice of the truncation orders, M<sub>x</sub> and M<sub>y</sub>, strongly depends on the values of D<sub>x</sub> and D<sub>y</sub>. A large-scale problem involves large values of D<sub>x</sub> and/or D<sub>y</sub>, yielding large truncation orders M<sub>x</sub> and M<sub>y</sub> as well as prohibitive computation time and memory storage capacity. This is one of the main limitations of the Fourier modal method when a large-scale problem issue is targeted such as the one addressed when designing some classes of photonics devices as metasurfaces. In order to accurately design highly efficient metasurfaces, near field (evanescent modes) are required to account for complex interactions among its sub-wavelength-scale constituents (meta-atoms). These evanescent modes, which are high-frequency harmonic modes, generally are highly sensitive to the truncation orders M<sub>x</sub> and M<sub>y</sub>. The larger the size of the computational domain, the larger the values of the truncation orders that should be taken into account in order to describe the rapid oscillations of the evanescent modes. However, increasing M<sub>x</sub> and M<sub>y</sub> yields high computational cost and memory storage capacity. One of the trivial ways to compute high harmonic modes without increasing the truncation orders M<sub>x</sub> and M<sub>y</sub> is to reduce the size of the computational domain. This issue is addressed in the next subsection.

## 2.3. Domain Decomposition Spectral Method: Direct Transform

Here, we aim to find an approximation |f̃⟩ of |f⟩ of Equation (8) on a domain Ω̃ ⊂ Ω, of width d<sub>x</sub>, d<sub>y</sub> with respect to (O, x) resp. (O, y), assuming that a set of eigenvalues and eigenvectors are known for a boundary value problem stated on Ω̃:

$$\tilde{\mathcal{L}}|\tilde{\Phi}_{pq}\rangle = -\tilde{\gamma}^2|\tilde{\Phi}_{pq}\rangle. \qquad (10)$$

We set |f̃⟩ as

$$|f\rangle = \tilde{A}^{pq} \tilde{\Phi}_{pq}^{\eta\mu} |\tilde{e}_{\eta\mu}\rangle, \qquad (11)$$

where

$$\tilde{e}_{\eta\mu}(x, y) = e^{-ik_0\tilde{\alpha}_{\eta}x} e^{-ik_0\tilde{\beta}_{\mu}y} \qquad (12)$$
with $\tilde{\alpha}_{[\eta]} = \alpha_0 + [-\eta_x : \eta_x]\lambda/d_x$ and $\tilde{\beta}_{[\mu]} = \beta_0 + [-\mu_y : \mu_y]\lambda/d_y$. The weighted coefficients $\tilde{A}^{pq}$ (Equation (11)) are computed using a Galerkin-weighted-residual method. It can be easily shown that:

$$[\tilde{A}^{pq}] = [\mathcal{P}_{kl}^{vw}][A^{kl}], \qquad (13)$$

where $\mathcal{P}$ is equal to:

$$[\mathcal{P}_{kl}^{vw}] = [\tilde{\Phi}_{pq}^{\eta\mu}]^{-1} [\langle \tilde{e}_{vw}^{\dagger}, \tilde{e}_{\eta\mu} \rangle]^{-1} [\langle \tilde{e}_{vw}^{\dagger}, e_{nm} \rangle] [\Phi_{kl}^{nm}]. \qquad (14)$$

Here, $\mathcal{P}$ clearly appears as the projector of the space spanned by $(|\Phi_{kl}\rangle = \Phi_{kl}^{nm}|e_{nm}\rangle)_{kl}$ onto the subspace spanned by $(|\tilde{\Phi}_{pq}\rangle = \tilde{\Phi}_{pq}^{\eta\mu}|\tilde{e}_{\eta\mu}\rangle)_{pq}$. In the current case of orthonormal Fourier basis functions, the set of $|\tilde{e}_{nm}\rangle$ satisfies $\langle \tilde{e}_{vw}^{\dagger}, \tilde{e}_{nm} \rangle = \delta_{v-n}\delta_{w-m}$. Therefore, Equation (14) yields

$$[\mathcal{P}_{kl}^{vw}] = [\tilde{\Phi}_{pq}^{\eta\mu}]^{-1} [\langle \tilde{e}_{vw}^{\dagger}, e_{kl} \rangle] [\Phi_{kl}^{nm}]. \qquad (15)$$

Here, we define the inner product of two periodic or pseudo-periodic functions $g_n$ and $h_m$ defined on an interval $[x_a, x_b]$ as:

$$\langle g_n^\dagger, h_m \rangle = \frac{1}{x_b - x_a} \int_{x_a}^{x_b} g_n^\dagger(x) h_m(x) dx, \qquad (16)$$

where $g_n^\dagger$ denotes the complex conjugate of $g_n$. The inner product of Equation (15) may then be expressed as

$$\left\{ \begin{array}{l} \langle \tilde{e}_v^\dagger, e_k \rangle_{x \in [a, b]} = e^{-i(x_b+x_a)\tau_x} \text{sinc}[(x_b - x_a)\tau_x] \\ \langle \tilde{e}_w^\dagger, e_l \rangle_{y \in [a, b]} = e^{-i(y_b+y_a)\tau_y} \text{sinc}[(y_b - y_a)\tau_y] \end{array} \right. . \qquad (17)$$

where $\tau_x = (\frac{v}{d_x} - \frac{k}{D_x})$, $\tau_y = (\frac{w}{d_y} - \frac{l}{D_y})$. The definition of the inner product can be extended to the case of a two-dimensional (2D) problem by using the Kronecker product:

$$\langle \tilde{g}_{vw}^{\dagger}, h_{kl} \rangle = \langle \tilde{g}_{v}^{\dagger}, h_{k} \rangle \cdot \langle \tilde{g}_{w}^{\dagger}, h_{l} \rangle. \qquad (18)$$

Let us consider now that the bounded plane $\Omega$ is squared into several elementary sub-cells $\Omega = \cup_{i,j} \tilde{\Omega}_{ij}$, with $\tilde{\Omega}_{ij} = [x_i, x_{i+1}] \times [y_j, y_{j+1}]$. Let us assume the known spectral coefficients $\tilde{A}^{pq}(i, j)$, related to each sub-domain $\tilde{\Omega}_{ij}$ such as any restriction of the expanded function $\tilde{f}_{ij}$ describing any component of the electromagnetic field on $\tilde{\Omega}_{ij}$, is written as:

$$|\tilde{f}_{ij}\rangle = \tilde{A}^{pq}(ij)\tilde{\Phi}_{pq}^{\eta\mu}|\tilde{e}_{\eta\mu}\rangle. \qquad (19)$$

The proposed DDSM allows us to obtain access to the near-field (i.e., electric and magnetic field) distribution of the device under study using Equation (19). However, in modeling an electromagnetic device, far-field information may be required in order to capture an accurate picture of the structure's performance. Within the framework of the spectral method presented here, the problem we face consists of reconstructing the spectral coefficients of the electromagnetic field on the whole domain $\Omega$, knowing the coefficients on each sub-domain $\tilde{\Omega}_{ij}$ obtained from the direct spectral decomposition transformation. This inverse transform issue will be detailed in the next subsection.
## 2.4. Domain Decomposition Spectral Method: Inverse Transform

The process we denoted as inverse transform consists of computing the global coefficient $A^{kl}$. These spectral coefficients are required to describe any electromagnetic field component function $f$ in the whole space as:

$$f(x, y, z) = \sum_{kl} A^{kl} e^{-ik_0 \gamma_{kl} z} \sum_{n,m} \Phi_{kl}^{nn} e_{nm}(x, y). \qquad (20)$$

Using once again the Galerking-weighted-residual method, one can demonstrate that the set of $A^{kl}$ is obtained from $\tilde{A}^{pq}(ij)$ throughout the following matrix relation:

$$[A^{kl}] = [ \begin{matrix} Q(1,1) & Q(1,2) & \dots & Q(N,N) \end{matrix} ] \begin{bmatrix} \tilde{A}^{pq}(1,1) \\ \tilde{A}^{pq}(1,2) \\ \vdots \\ \tilde{A}^{pq}(N,N) \end{bmatrix} \qquad (21)$$

with

$$[Q_{vw}^{kl}(i,j)] = [\Phi_{kl}^{nm}]^{-1} [\langle e_{kl}^{\dagger}, \tilde{e}_{vw} \rangle(i,j)] [\tilde{\Phi}_{pq}^{\eta\mu}]. \qquad (22)$$

We summarize in Figure 2 the main steps of DDSM when it is applied to compute near and far fields of a photonics device in the framework of modal methods. For a better understanding, we also report in Figure 3 the flowchart and parallelization structure of the employed method [24]. As a first step, the computational domain $\Omega$ is squared into sub-domains $\tilde{\Omega}_{ij}$. Secondly, using the direct transform, the incident field is projected on the eigenvectors defined on each $\tilde{\Omega}_{ij}$ in the incident medium. Forward simulations can then be independently performed for each sub-domain, yielding an approximation of the near field on each $\tilde{\Omega}_{ij}$. Afterwards, the set of simulated near fields on each sub-cell is stitched together to reproduce the electromagnetic field on the whole surface. Finally, the inverse transform can be used to obtain the spectral coefficients of the field related to the whole domain if needed.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Sketch of the direct and inverse transform of DDSM applied to the near and far fields simulation of a photonics device.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 3. Flowchart and parallelization structure of the direct and inverse transform of DDSM applied to the near and far fields simulation of a photonics device.</figcaption></figure>

Having introduced the mathematical formalism of the proposed approach, we are now ready to go further into some applications.

### 3. Applications

#### 3.1. DDSM Applied to a Dipole Field Expansion in the Framework of Modal Method

The purpose of the present section is merely to show what the DDSM looks like when applied to a simple case. In a Cartesian coordinates system (O, x, y, z), we are interested in the footprint of the field radiated by an electric dipole located at the origin O(0, 0, 0) in a given plane z = -z_s, as shown in Figure 4a. This plane will be the 2D computation domain Ω. This kind of simulation may be useful while performing the so-called adjoint simulation in an optimization problem of metalenses. In the case of the forward simulation, the device under optimization is excited by a planewave, while the adjoint simulation required the use of a dipole source located at the desired foci Figure 4b. In the framework of modal expansion, any component of the electromagnetic field of this dipole source can be fully described in the whole space by a set of known expansion coefficients A<sup>pq</sup>, B<sup>pq</sup> and a set of eigenvectors |Φ<sup>nm</sup><sub>pq</sub>⟩ as:

$$|f(z)\rangle = [A^{pq}\Phi_{pq}^{nm}e^{ik_0\gamma_{pq}z} + B^{pq}\Phi_{pq}^{nm}e^{-ik_0\gamma_{pq}z}]|e_{nm}\rangle. \quad (23)$$

For our illustration, only the first part of the above expression describing the downward wave is kept:

$$|f(z)\rangle = A^{pq}\Phi_{pq}^{nm}e^{ik_0\gamma_{pq}z}|e_{nm}\rangle. \quad (24)$$

The set of eigenvectors |Φ<sup>nm</sup><sub>pq</sub>⟩ and their associated eigenvalues γ<sub>pq</sub> can be computed thanks to Equation (1). Practical details for the numerical computation of the weighted coefficients [A<sup>pq</sup>] and [B<sup>pq</sup>] can be found in [25]. We report in Figure 5a,b, the phase and the real part, respectively, of E<sub>x</sub>(x, y, z_s = -20λ), in vacuum ∀(x, y) ∈ [-0.5D<sub>x</sub>, 0.5D<sub>x</sub>] × [-0.5D<sub>y</sub>, 0.5D<sub>y</sub>]. The truncation orders M<sub>x</sub> and M<sub>y</sub> are set to M<sub>x</sub> = M<sub>y</sub> = 31, D<sub>x</sub> = D<sub>y</sub> = 26.67λ.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Sketch of a dipole field decomposition. The fundamental strategy of the DDSM to compute the footprint of a dipole field on an arbitrary large-area surface consists in dividing the surface into elementary sub-cells. Any component of the dipole field can then be locally simulated, on each sub-cell. Here, the surface $z = -z_s$ is subdivided into 3 × 3 sub-cells. (a) Dipole field decomposition and recomposition. (b) Sketch of simulation of metasurface under a dipole illumination.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Footprints of the phase (a) and the real part (b) of an $E_x$ electric dipole on the plane $z = -20\lambda$. The numerical computation is performed on the whole domain $[-0.5D_x, 0.5D_x] \times [-0.5D_y, 0.5D_y]$. The operating wavelength: $\lambda = 0.64 \mu m$, $D_x = D_y = 26.67\lambda$, truncation orders $M_x = M_y = 31$.</figcaption></figure>

A typical well-defined wavefront of a dipole source can be distinguished. These results will serve as a baseline and will be compared to those obtained with the proposed domain decomposition method. First, the plane $z = -z_s$ is squared into several elementary sub-cells $\Omega = \cup_{i,j} \tilde{\Omega}_{ij}$, with $\tilde{\Omega}_{ij} = [x_i, x_{i+1}] \times [y_j, y_{j+1}]$. Then, the eigenvalue equation Equation (10), obtained from Maxwell's equations, equipped with PMLs (perfectly matched coordinates), are numerically and independently solved in each sub-cell. From a geometrical point of view, in a homogenous medium, the eigenvalue Equation (10) only depends on the size of the computation area. Consequently, for a set of sub-cells $\tilde{\Omega}_{ij}$ of the same size, only one numerical simulation based on Equation (10) is performed to obtain the sequence of eigenvectors $|\tilde{\Phi}_{pq}\rangle$. The truncation orders are set to $m_x = m_y = 10$. The sequence of the spectral coefficients $[\tilde{A}^{pq}(i, j)]$ in each sub-cell $\tilde{\Omega}_{ij} = [x_i, x_{i+1}] \times [y_j, y_{j+1}]$ is then obtained from the weighted coefficients $[A^{pq}]$ using Equation (13). The restriction of any component of the electromagnetic field $|f_{i,j}(z)\rangle$ on a sub-cell $\tilde{\Omega}_{ij}$ is then written as:

$$|\tilde{f}_{(i,j)}(z)\rangle = \tilde{A}^{pq}(i, j) \tilde{\Phi}_{pq}^{nm} e^{ik_0 \tilde{\gamma}_{pq} z} |\tilde{e}_{nm}\rangle. \quad (25)$$

The set of the simulated fields on each sub-cell is then stitched together (geographically, without any post-processing) to reproduce the electromagnetic field on the whole surface $\Omega$.
Figure 6a,b represent the real parts and the phase, respectively, computed using Equation (25) by dividing the computational domain $[-0.5D_x, 0.5D_x] \times [-0.5D_y, 0.5D_y]$ into 3 × 3 sub-domains. In Figure 7a,b, this number is extended to 5 × 5 sub-domains. In all these figures, where no post-processing treatment is applied, the field distribution in the sub-cells areas and the effect of the PMLs areas are clearly shown. Comparing these results with the baseline results of Figure 5a,b, one can also remark that out of the PMLs areas (ratio PMLs-area-wide/Ω<sub>ij</sub>-wide<< 1/10), the footprints of the dipole source are well-described on each sub-cell $\tilde{\Omega}_{ij}$ by the sequence of computed local spectral coefficients [$\tilde{A}^{pq}(i, j)$]. The method then shows its ability to represent the footprint of a given electromagnetic field component. This leads to the next step where we are ready to deploy the proposed concept to compute the electromagnetic field diffracted by a large-scale device. Three potential devices are investigated: a dielectric Fresnel Plate Zone (FPZ), a metalens made of dielectric cylindrical waveguides (nanorodes) with variables cross-sections and a Pachatarnam-Berry phase dielectric metalens.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Footprints of the phase (a) and the real part (b) of an $E_x$ electric dipole in the plane $z = -20\lambda$. The whole domain $[-0.5D_x, 0.5D_x] \times [-0.5D_y, 0.5D_y]$ is divided into a 3 × 3 sub-cell $\Omega_{ij}$. Based on the proposed parallel strategy, the numerical computation is performed on each sub-cell $\Omega_{ij}$. The operating wavelength: $\lambda = 0.64$ µm, $d_x = d_y = 8.89\lambda$, truncation orders $m_x = m_y = 10$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Footprints of the phase (a) and the real part (b) of an $E_x$ electric dipole in the plane $z = -20\lambda$. The whole domain $[-0.5D_x, 0.5D_x] \times [-0.5D_y, 0.5D_y]$ is divided onto a 5 × 5 sub-cell $\Omega_{ij}$. The operating wavelength: $\lambda = 0.64$ µm, $d_x = d_y = 5.33\lambda$, truncation orders $m_x = m_y = 6$.</figcaption></figure>

## 3.2. Modal Method Applied to a Large Binary Fresnel Plate Zone (Binary FPZ)

In this section, the first proposed case will be presented. The task consists of a forward simulation of both the near and far-field components of a 2D dielectric binary FPZ with a height $h$ under a linearly polarized incident plane wave, as shown in Figure 8a. The
FPZ consists of a set of concentric dielectric rings, known as Fresnel zones, operating as a dielectric lens using diffraction to focus light into a given focal point. The considered dielectric material has a refractive index of 1.99 deposited on a SiO₂ substrate (refractive index = 1.45). Its height is set to h = 400 nm for an operating wavelength of λ = 532 nm. To obtain a constructive interference at the focus point, the radius of the Fresnel zones should satisfy the equation:

$$R_n = k_0 \sqrt{n\lambda f + \left(\frac{n\lambda}{2}\right)^2}, \qquad (26)$$

where f is the lens' focal length, λ is the operating wavelength and k₀ = 2π/λ is the wavenumber in the focusing medium. The aperiodic Fourier modal method (AFMM) [26–28], a full-wave analysis based on the Fourier modal method equipped with PMLs (perfectly matched layers) [29–33], is used to compute the electromagnetic field components. It is well known that in the case of the FMM or AFMM, the total computation time grows as O(N⁶) (bi-periodic problem) with respect to the total truncation order N required for all the field components (i.e, the four field components used for the continuity equation). On the other hand, memory consumption time grows as O(N⁴). Therefore, scaling the computation [−0.5Dₓ, 0.5Dₓ] × [−0.5Dᵧ, 0.5Dᵧ] surface into some sub-sections will certainly provide advantages to scaling both the computational time and the required memory capacities as (2mₓ + 1)(2mᵧ + 1)/(2Mₓ + 1)(2Mᵧ + 1). In our illustration, the lens' focal length is set to 70λ, and the number of FPZ rings is 15, yielding a 100λ-width device. To apply the DDSM to the proposed geometry, the 100λ × 100λ × 400 nm structure is divided into 5 × 5 × 400 nm sub-cells of dₓ × dᵧ × h = 6.7277λ × 6.7277λ × 400 nm sub-voxels. Each sub-voxel is treated independently, with very low truncation numbers mₓ and mᵧ, thereby enabling the simulation to be more feasible in terms of memory and time consumption. In the current case (2mₓ + 1) × (2mᵧ + 1) = (2 × 12 + 1) × (2 × 12 + 1), Fourier harmonics are enough to describe the electromagnetic field behavior on each sub-domain. On the contrary, the electromagnetic field expansion considering the whole domain Ω involves at least (2Mₓ + 1) × (2Mᵧ + 1) = (2 × 62 + 1) × (2 × 62 + 1) Fourier harmonics yielding a simulation that is extremely consuming in both time and memory. The insets in Figure 8c,d display the top-view layout of the FPZ lens and the (x, y)-plane phase distribution of the transmitted electric field on the emerging surface of the lens, respectively. Regarding Figure 8d, one can remark the full [−π, π] phase coverage characterizing a binary phase zone plate. Making good use of the inverse transform post-processing, we can obtain access to field distribution in the whole 3D space. In particular, we compute and display in Figure 8e the transmitted field in the 3D space. Regarding this result, one can remark that the incident electric field intensity at the foci considerably increases up to 3500 times. As reported in Table 1, the computational times for the simulation of the near field (i.e., emerged field on the FPZ surface) using the Fourier modal method 2D-AFMM implemented with mₓ = mᵧ Fourier basis functions (so that the size of the eigenvalues equations matrix is 2(2mₓ + 1) × 2(2mₓ + 1)) are 35.07 s for mₓ = mᵧ = 9, 57.29 s for mₓ = mᵧ = 10, 90.77 s for mₓ = mᵧ = 11, and 152.67 s for mₓ = mᵧ = 12 (Mₓ = Mᵧ = 62). These simulations are performed on a classical personal computer DELL PRECISION 3640, with an intel CORE i9 (3.10 Ghz) processor. Note that since the [−0.5Dₓ, 0.5Dₓ] × [−0.5Dᵧ, 0.5Dᵧ] surface is strongly inter-meshed, 2¹² × 2¹² nodes have been used both for the highly pixelated profile of the FPZ lens surface and also for the computation nodes of the near field.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** A forward simulation of a 100λ-width FPZ using the DDSM. The computation domain is subdivided into 5 × 5 sub-domains Ω<sub>ij</sub>. (c) presents the sketch of the diffractive Fresnel lens, (d) shows the phase distribution on the top face of the lens and (e) displays the transmitted field. Numerical parameters: λ = 0.532 µm, z<sub>s</sub> = 70λ, n = 31 zones (15 rings), D<sub>x</sub> = D<sub>y</sub> = 100λ h = 400 nm, refractive index of the lens material 1.99, refractive index of the substrate 1.45, m<sub>x</sub> = m<sub>y</sub> = 12, (M<sub>x</sub> = M<sub>y</sub> = 62).</figcaption></figure>

<table><caption>**Table 1.** Panel of the computation times for the simulation of the FPZ near field on each sub-cell using the aperiodic Fourier modal method (AFMM). Simulation is performed on a classical computer DELL PRECISION 3640 with a processor intel CORE i9 (3.10 Ghz). The whole [-0.5D<sub>x</sub>, 0.5D<sub>x</sub>] × [-0.5D<sub>y</sub>, 0.5D<sub>y</sub>] surface is strongly inter-meshed 2<sup>12</sup> × 2<sup>12</sup> nodes have been used both for the highly pixelated profile the FPZ lens surface and also for the computation nodes of the near field.</caption><tr><td>m<sub>x</sub> = m<sub>y</sub></td><td>9</td><td>10</td><td>11</td><td>12</td></tr><tr><td>computation times (s)</td><td>35.07</td><td>57.29</td><td>90.77</td><td>152.67</td></tr></table>

### 3.3. Modal Method and Large Metalens Consisted of Set of Different Cross-Sections Nanorods

The considered second example is a dielectric metalens consisting of a set of sub-wave-length-scale dielectric Si waveguides with a refractive index of 3.6082, which is deposited on a substrate with a refractive index of 2.4626 in Figure 8b. At a given point (x, y) on the metalens, the radius of each waveguide is chosen so that the spatial distribution of the phase profile θ(x, y) follows the following equation:

$$ \theta(x, y) = k_0[f - \sqrt{x^2 + y^2 + f^2}]. \qquad (27) $$

The metalens height is set to h = 5.57 µm and the operating wavelength is λ = 7 µm. The lens is designed to focus on a linearly polarized incident plane wave at a focus length of f = 10λ.
We perform forward simulation on a $\Omega = 12.8571\lambda \times 12.8571\lambda$-width device. The results are presented in Figure 9. The device surface is squared into 3 × 3 sub-domains $\Omega_{ij}$, as shown in Figure 9a. The truncation orders are set to $m_x = m_y = 17$, which is equivalent to $M_x = M_y = 52$ on $\Omega$. A $[-\pi, \pi]$ coverage is clearly distinguished in Figure 9b, where the phase distribution of the emerged electric field on the top face of the metalens is presented. Applying the inverse transform allows one to obtain the spectral coefficients of the global structure. The emerged near and far fields can then be computed and the focus action of the metalens can be visualized, as shown in Figure 9c.

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Forward simulation of a metalens consisting of sub-wavelength dielectric nanorodes with different cross-section. Power distributions of the focusing spot in the $(x, y = 0, z)$ plane, $\lambda = 7$ µm, $z_s = 10\lambda$, domain $\Omega = 12.8571\lambda \times 12.8571\lambda = 3 \times 3$ sub-cell $\Omega_{ij}$, $h = 5.57$ µm, materials: dielectrics Si 3.6082 deposited on material with refractive index (2.4626), $m_x = m_y = 17$, $M_x = M_y = 52$.</figcaption></figure>

In Figure 10a, and for the same device, we increase the number of the sub-domains from 3 × 3 up to 5 × 5. This allows one to decrease the truncation orders down to $m_x = m_y = 10$ ($M_x = M_y = 52$). As shown in Figure 10b,c, results related to this new partition are still accurate.

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Forward simulation of a metalens consisting of sub-wavelength dielectric nanorodes with different cross-section. Power distributions of the focusing spot in the $(x, y = 0, z)$ plane, $\lambda = 7$ µm, $z_s = 10\lambda$, domain $\Omega = 12.8571\lambda \times 12.8571\lambda = 5 \times 5$ sub-cell $\Omega_{ij}$, $h = 5.57$ µm, materials: dielectrics Si $v = 3.6082$ deposited on material with refractive index 2.4626, $m_x = m_y = 10$, $M_x = M_y = 52$.</figcaption></figure>

### 3.4. Large-Scale Pachatarnam-Berry Phase-Based (PB) Metalens Simulation

In the third example, we are interested in a Pachatarnam-Berry phase metalens. It consists of an arrangement of rotated dielectric meta-atoms having a predefined shape and a refractive index of 1.99 deposited on a SiO₂ substrate (refractive index= 1.45). In our example, a rectangular shape is considered. The dimensions of the rectangular meta-atoms are (length = 300 nm, width = 105 nm) and the metalens height is $h = 400$ nm. These parameters are chosen such that the PB phase-based metasurface converts and focuses a fraction of the left-hand circularly polarized (LCP) incident light injected from the substrate into right-hand circular polarization (RCP). To obtain the desired in plane $(x, y)$ phase
distribution for focusing, each meta-atom located at position (x, y) in the plane is rotated at an angle θ(x, y) with respect to (O, x) axis such as

$$ \theta(x, y) = \frac{1}{2}k_0[f - \sqrt{x^2 + y^2 + f^2}]. \qquad (28) $$

The operating wavelength is set to λ = 0.532 µm. The simulation was run on the same personal computer as the previous case. Results are presented in Figure 11. The computation domain Ω = 23.30λ × 23.30λ is squared into 3 × 3 sub-domains leading to d_x = d_y = 7.76λ-width sub-cells. The phase distribution of the co-polarization E_x^RCP(x, y), E_y^RCP(x, y) components of the emerged RCP electric field at the top face of the metalens, obtained using the direct and the inverse transform, are presented in Figure 11c,d. A [-π, π] coverage is clearly distinguished. Applying the inverse transform yields in Figure 11b, an accurate reproduction of the transmitted power of the metalens.

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Forward simulation of a Ω = 23.3λ × 23.3λ-width Pachatarnam-Berry phase-based metalens. The lens is designed to have a focus length of f = 20λ. The number of sub-cells and truncation orders are set to 3 × 3 and m_x = m_y = 12 (M_x = M_y = 37), respectively. (b) displays the power distribution in the (x, y = 0, z) plane. The phase distribution of the co-polarization E_x^RCP(x, y), E_y^RCP(x, y) components of the emerged RCP electric field at the top face of the metalens, obtained using the direct and the inverse transform, are presented in (c,d). A [-π, π] coverage is clearly distinguished. Numerical parameters: λ = 0.532 µm, focus length z_s = 20λ, D_x = D_y = 23.30λ, h = 400 nm.</figcaption></figure>

Finally, we increase the PB-based metalens up to 30.6λ × 30.6λ width. One can easily admit that increasing the lens size will lead to a more efficient device. The spatial distribution of the transmitted power and phases of the electric field of the cross-polarization are shown in Figure 12. As expected, both the power energy flow and the resolution of the designed lens are increased compared to the previous case of a 23.30λ width metalens (Figure 11). However, increasing the lens' geometrical features requires increasing the
number of sub-cells in order to keep a low value of the truncation numbers $m_x$ and $m_y$. The device is divided into 5 × 5 sub-cells of 6.16λ × 6.16λ, where the truncation orders are set to $m_x = m_y = 10$.

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Forward simulation of a Ω = 30.8λ × 30.8λ-width Pachatarnam-Berry phase-based metalens. The lens is designed to have a focus length of *f* = 20λ. Power distributions of the focusing spot in the (*x*,*y* = 0, *z*) plane. Phase distribution of *E*<sub>*x*</sub>(*x*,*y*) component of the cross-polarization (RCP) on the top of the PB-based metalens. The number sub-cells and truncation orders are set to 5 × 5 and *m*<sub>*x*</sub> = *m*<sub>*y*</sub> = 10, respectively. Numerical parameters: λ = 0.532 µm, focus length *z*<sub>s</sub> = 20λ, *h* = 400 nm.</figcaption></figure>

**4. Conclusions and Outlook**

To summarize, we introduce an approach, namely a domain decomposition spectral
method DDSM, based on spectral methods, which can be optimally and simply integrated in
a parallel computing scheme to provide a solution for large-size systems. On one hand, the
power of a spectral method relies on the use of a few coefficients to describe both near and
far-electromagnetic fields in a three-dimensional space. On the other hand, modal methods
provide a better representation of physical phenomena involving modal interactions such as
resonance phenomena through dielectric and metallic periodic and aperiodic metasurfaces.
While many previous works, mainly based on spatial domain methods such as finite
difference time and frequency domain methods, finite elements methods, etc., address
large-scale problems efficiently, despite some previous developments, the simulation of
large-scale aperiodic devices remains a challenging problem for full spectral modal methods.
We expect that our work will help to tackle this issue. In order to show the efficiency of the
proposed work, we provided numerous examples, such as an electric dipole simulation
and some monochromatic metalenses illuminated under both circular and linear polarized
incident plane waves. In fact, the proposed method is based on two processes: the direct
and inverse transforms, which allow one to compute, efficiently, both near and far-field
spectrums. The DDSM direct transform consists of decomposing the spectral coefficients
of the total field onto some spectral coefficients related to the sub-domains. Meanwhile, the inverse transform consists of reconstructing the spectral coefficients of the total field starting from the coefficients of decomposition of each field obtained from the direct spectral decomposition transformation.

Our domain decomposition spectral method provides a promising pathway for handling large aperiodic, but it is still limited, and further enhancements are needed to extend this method to a very, very large-scale problem. Indeed, the process of decomposition of the incident field requires the resolution of a unique eigenvalue equation stated on a very-very wide domain. However, in an electromagnetism problem, the incident field is not unknown to the problem. This field is generally well known, and its modal characteristic can be efficiently and reliably computed based on interpolating methods.

Therefore, we plan to calculate the modal characteristics of incident fields based on neural network models. This strategy could ultimately and reliably increase the performance of the method in the case of very, very large-scale problems.

**Author Contributions:** Formal analysis, M.C.L. and B.G.; Funding acquisition, P.B.; Validation, G.G. and J.-P.P.; Visualization, F.P., G.A.A. and L.D.; Writing—review & editing, K.E. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Further Computational Details about the Method's Implementation along with a Version of the Code. Available online: [photonicsnum.com](http://photonicsnum.com) (accessed on 26 September 2022).

**Acknowledgments:** The authors acknowledge the KITP program "Emerging Regimes and Implications of Quantum and Thermal Fluctuational Electrodynamics" 2022, where part of this work has been completed. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.

**Conflicts of Interest:** The author declares there are no conflict of interest.

## References

1. Gander, M.J.; Zhang, H. A class of iterative solvers for the Helmholtz equation: Factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized Schwarz methods. *SIAM Rev.* **2019**, *61*, 1–76. [CrossRef]

2. Zepeda-Núñez, L.; Demanet, L. Nested domain decomposition with polarized traces for the 2d Helmholtz equation. *SIAM J. Sci. Comput.* **2018**, *40*, B942–B981.

3. Stolk, C.C. A rapidly converging domain decomposition method for the Helmholtz equation. *J. Comput. Phys.* **2013**, *241*, 240–252. [CrossRef]

4. Chen, Z.; Xiang, X. A source transfer domain decomposition method for Helmholtz equations in unbounded domain. *SIAM J. Numer. Anal.* **2013**, *51*, 2331–2356. [CrossRef]

5. Vouvakis, M.; Zhao, K.; Lee, J.-F. Modeling large almost periodic structures using a non-overlapping domain decomposition method. In Proceedings of the IEEE Antennas and Propagation Society Symposium, Monterey, CA, USA, 20–25 June 2004; pp. 343–346.

6. Dolean, V.; Gander, M.J.; Gerardo-Giorda, L. Optimized schwarz methods for Maxwell's equations. *SIAM J. Sci. Comput.* **2009**, *31*, 2193–2213. [CrossRef]

7. Lin, Z.; Johnson, S.G. Overlapping domains for topology optimization of large-area metasurfaces. *Opt. Express* **2019**, *27*, 32445–32453. [CrossRef] [PubMed]

8. Phan, T.; Sell, D.; Wang, E.W.; Doshay, S.; Edee, K.; Yang, J.; Fan, J.A. High-efficiency, large-area, topology-optimized metasurfaces. *Light. Sci. Appl.* **2019**, *8*, 48. [CrossRef]

9. Taflove, A.; Hagness, S.C. *Computational Electrodynamics*; Artech House Publishers: Norwood, MA, USA, 2000; Volume 28.

10. Rumpf, R.C. Simple implementation of arbitrarily shaped total-field/scattered-field regions in finite-difference frequency-domain. *Prog. Electromagn. Res.* **2012**, *36*, 221–248. [CrossRef]

11. Reddy, J.N. *Introduction to the Finite Element Method*; McGraw-Hill Education: New York, NY, USA, 2019.

12. Babuska, I.; Szabo, B. On the rates of convergence of the finite element method. *Int. J. Numer. Methods Eng.* **1982**, *18*, 323–341. [CrossRef]

13. Knop, K. Rigorous diffraction theory for transmission phase gratings with deep rectangular grooves. *J. Opt. Soc. Am. A* **1978**, *68*, 1206–1210. [CrossRef]

14. Granet, G.; Guizal, B. Efficient implementation of the coupled-wave method for metallic lamellar gratings in TM polarization. *J. Opt. Soc. Am. A* **1996**, *13*, 1019–1023. [CrossRef]
15. Lalanne, P.; Morris, G.M. Highly improved convergence of the coupled-wave method for TM polarization. *J. Opt. Soc. Am. A* **1996**, *13*, 779–784. [CrossRef]

16. Li, L. Use of Fourier series in the analysis of discontinuous periodic structures. *J. Opt. Soc. Am. A* **1996**, *13*, 1870–1876. [CrossRef]

17. Granet, G. Reformulation of the lamellar grating problem through the concept of adaptive spatial resolution. *J. Opt. Soc. Am. A* **1999**, *16*, 2510–2516. [CrossRef]

18. Lalanne, P. Numerical performance of finite-difference modal methods for the electromagnetic analysis of one- dimensional lamellar gratings. *J. Opt. Soc. Am. A* **2000**, *17*, 1033–1042. [CrossRef]

19. Boyd, J. *Chebyshev and Fourier Spectral Methods*, 2nd ed.; Dover Books on Mathematics; Dover Publications: Mineola, NY, USA, 2013; ISBN 9780486141923.

20. Edee, K.; Plumey, J.-P.; Guizal, B. Unified Numerical Formalism of Modal Methods in Computational Electromagnetics and Latest Advances: Applications in Plasmonics. In *Advances in Imaging and Electron Physics*; Elsevier: Amsterdam, The Netherlands, 2016; Volume 197, Chapter 2.

21. Semenikhin, I. Improving accuracy of the numerical solution of Maxwell's equations by processing edge singularities of the electromagnetic field. *J. Comput. Phys.* **2021**, *441*, 110440. [CrossRef]

22. Hammond, A.M.; Oskooi, A.; Chen, M.; Lin, Z.; Johnson, S.G.; Ralph, S.E. High-performance hybrid time/frequency-domain topology optimization for large-scale photonics inverse design. *Opt. Express* **2022**, *30*, 4467–4491. [CrossRef]

23. Lin, Z.; Liu, V.; Pestourie, R.; Johnson, S.G. Topology optimization of free-form large-area metasurfaces. *Opt. Express* **2019**, *27*, 15765–15775. [CrossRef]

24. Further Computational Details about the Method’s Implementation along with a Version of the Code. Available online: *photonicsnum.com* (accessed on 26 September 2022).

25. Plumey, J.-P.; Edee, K.; Granet, G. Modal expansion for the 2D Greens function in a non orthogonal coordinate system. *Prog. Electromagn. Res.* **2006**, *59*, 101–112. [CrossRef]

26. Silberstein, E.; Lalanne, P.; Hugonin, J.-P.; Cao, Q. Use of gratings in integrated optics. *J. Opt. Soc. Am. A* **2001**, *18*, 2865–2875. [CrossRef] [PubMed]

27. Hugonin, J.P.; Lalanne, P. Perfectly matched layers as nonlinear coordinate transforms: A generalized formalization. *J. Opt. Soc. Am. A* **2005**, *22*, 1844–1849. [CrossRef] [PubMed]

28. Edee, K.; Granet, G.; Plumey, J.-P. Complex coordinate implementation in the curvilinear coordinate method: Application to plane-wave diffraction by nonperiodic rough surfaces. *J. Opt. Soc. Am. A* **2007**, *24*, 1097–1102. [CrossRef]

29. Bérenger, J. A perfectly matched layer for the absorption of electromagnetic waves. *J. Comput. Phys.* **1994**, *114*, 185–200. [CrossRef]

30. Sacks, Z.S.; Kingsland, D.M.; Lee, R.; Lee, J.-F. A perfectly matched anisotropic absorber for use as an absorbing boundary condition. *IEEE Trans. Antennas Propag.* **1995**, *43*, 1460–1463. [CrossRef]

31. Chew, W.C.; Weedon, W.H. A 3D perfectly matched medium from modified Maxwell's equations with stretched coordinates. *Microw. Opt. Technol. Lett.* **1994**, *7*, 599–604. [CrossRef]

32. Chew, W.C.; Jin, J.M.; Michielssen, E. Complex coordinate stretching as a generalized absorbing boundary condition. *Microw. Opt. Technol. Lett.* **1997**, *15*, 363–369. [CrossRef]

33. Teixeira, F.L.; Chew, W.C. Differential forms, metrics, and 435 the reflectionless absorption of electromagnetic waves. *J. Electromagn. Waves Appl.* **1999**, *13*, 665–686. [CrossRef]