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

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.


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][2][3][4][5][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

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.

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 Ω = [−0.5D x , 0.5D x ] × [−0.5D y , 0.5D y ] ⊂ R 2 . We assume that the physical parameters of the medium, namely the relative permittivity ε(x, y) and the permeability µ(x, y) functions, do not depend on the z variable. In this case, any component of the electromagnetic field can be written as Φ(x, y)exp(−ik 0 γz), where k 0 = ω √ ε 0 µ 0 is the wavenumber in the vacuum, ω is the angular frequency and k 0 γ is the propagation constant along the z direction. ε 0 and µ 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: The eigenvalue Equation (1) is solved using the spectral method briefly described in the next subsection.

Spectral Method and Modal Method
In electromagnetism, a spectral method aims to describe a function |Φ representing any electromagnetic field component with a set of finite expansion coefficients Φ nm on a selected basis also known as expansion functions |e nm : by using Einstein's convention. Expansion or weighted coefficients Φ nm are the spectrum of the expanded function Φ with respect to the sequence of expansion vectors |e nm . 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 z 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 , I z and I z , I z and I (4) z with respect to the propagation direction that is assumed to be (O, z).
In any of the layers I (l) z , 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 |Φ kl of an eigenvalue equation similar to Equation (1): The completeness of the set of eigenvectors (|Φ kl ) kl allows one to expand any function | f defined on a domain Ω as a linear combination of the eigenvectors |Φ kl : Regarding Equation (5), the modal method can be viewed as a spectral method using as basis functions a set of eigenvectors |Φ kl of the operator L. The Fourier modal method (FMM) consists of expanding the eigenvectors |Φ kl in terms of a generalized Fourier basis: where |e nm = |e n ⊗ |e m , with Taking into account Equation (6), the vector | f of Equation (5) may be expressed in terms of the Fourier basis functions as: From a numerical point of view, only a finite number of (2M x + 1) × (2M y + 1), ((M x , M y ) ∈ N 2 ) Fourier coefficients is kept by selecting a finite sequence of α n and β n in the expression Equation (7): where D x (resp. D y ) is the period along the (O, x) (resp. (O, y)) axis and λ is the operating wavelength. The parameters α 0 and β 0 are linked to the incident wave. The choice of the truncation orders, M x and M y , strongly depends on the values of D x and D y . A largescale problem involves large values of D x and/or D y , yielding large truncation orders M x and M y 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 subwavelength-scale constituents (meta-atoms). These evanescent modes, which are highfrequency harmonic modes, generally are highly sensitive to the truncation orders M x and M y . 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 x and M y yields high computational cost and memory storage capacity. One of the trivial ways to compute high harmonic modes without increasing the truncation orders M x and M y is to reduce the size of the computational domain. This issue is addressed in the next subsection.

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 x , d y 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 Ω: We set | f as where e ηµ (x, y) = e −ik 0 α η x e −ik 0 β µ y with The weighted coefficients A pq (Equation (11)) are computed using a Galerkin-weighted-residual method. It can be easily shown that: where P is equal to: Here, P clearly appears as the projector of the space spanned by |Φ kl = Φ nm kl |e nm kl onto the subspace spanned by | Φ pq = Φ ηµ pq |ẽ ηµ pq . In the current case of orthonormal Fourier basis functions, the set of |ẽ nm satisfies ẽ † vw ,ẽ nm = δ v−n δ w−m . Therefore, Equation (14) yields 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: where g † n denotes the complex conjugate of g n . The inner product of Equation (15) may then be expressed as where The definition of the inner product can be extended to the case of a two-dimensional (2D) problem by using the Kronecker product: Let us consider now that the bounded plane Ω is squared into several elementary sub-cells ]. Let us assume the known spectral coefficients A pq (i, j), related to each sub-domain Ω ij such as any restriction of the expanded function f ij describing any component of the electromagnetic field on Ω ij , is written as: 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 Ω, knowing the coefficients on each sub-domain Ω ij obtained from the direct spectral decomposition transformation. This inverse transform issue will be detailed in the next subsection.

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: Using once again the Galerking-weighted-residual method, one can demonstrate that the set of A kl is obtained from A pq (ij) throughout the following matrix relation: 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 Ω is squared into sub-domains Ω ij . Secondly, using the direct transform, the incident field is projected on the eigenvectors defined on each Ω 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 Ω 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.  Having introduced the mathematical formalism of the proposed approach, we are now ready to go further into some applications.

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 pq , B pq and a set of eigenvectors |Φ nm pq as: For our illustration, only the first part of the above expression describing the downward wave is kept: The set of eigenvectors |Φ nm pq and their associated eigenvalues γ pq can be computed thanks to Equation (1). Practical details for the numerical computation of the weighted coefficients [A pq ] and [B pq ] can be found in [25]. We report in Figure 5a,b, the phase and the real part, re-  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 ] 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) on a sub-cell Ω ij is then written as: 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 Ω. 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 ] × [−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/Ω ij -wide<< 1/10), the footprints of the dipole source are well-described on each sub-cell Ω ij by the sequence of computed local spectral coefficients [ 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.

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 2 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: where f is the lens' focal length, λ is the operating wavelength and k 0 = 2π/λ is the wavenumber in the focusing medium. The aperiodic Fourier modal method (AFMM) [26][27][28], a full-wave analysis based on the Fourier modal method equipped with PMLs (perfectly matched layers) [29][30][31][32][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 6 ) (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 4 ). Therefore, scaling the computation [−0.5D x , 0.5D x ] × [−0.5D y , 0.5D y ] surface into some sub-sections will certainly provide advantages to scaling both the computational time and the required memory capacities as (2m x + 1)(2m y + 1)/(2M x + 1)(2M y + 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 x × d y × h = 6.7277λ × 6.7277λ × 400 nm sub-voxels. Each sub-voxel is treated independently, with very low truncation numbers m x and m y , thereby enabling the simulation to be more feasible in terms of memory and time consumption. In the current case (2m x + 1) × (2m y + 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 x + 1) × (2M y + 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 x = m y Fourier basis functions (so that the size of the eigenvalues equations matrix is

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 subwavelength-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: 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 Ω = 12.8571λ × 12.8571λ-width device. The results are presented in Figure 9. The device surface is squared into 3 × 3 sub-domains Ω 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 Ω. A [−π, π] 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.  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.

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 2 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 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 RCP x (x, y), E RCP y (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.  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.

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.