The FLAME-slab method for electromagnetic wave scattering in aperiodic slabs

The proposed numerical method,"FLAME-slab,"solves electromagnetic wave scattering problems for aperiodic slab structures by exploiting short-range regularities in these structures. The computational procedure involves special difference schemes with high accuracy even on coarse grids. These schemes are based on Trefftz approximations, utilizing functions that locally satisfy the governing differential equations, as is done in the Flexible Local Approximation Method (FLAME). Radiation boundary conditions are implemented via Fourier expansions in the air surrounding the slab. When applied to ensembles of slab structures with identical short-range features, such as amorphous or quasicrystalline lattices, the method is significantly more efficient, both in runtime and in memory consumption, than traditional approaches. This efficiency is due to the fact that the Trefftz functions need to be computed only once for the whole ensemble.


I. INTRODUCTION
A. Electromagnetic wave scattering in slab structures Scattering of electromagnetic waves by a patterned photonic slab is an important but challenging class of problems in computational electromagnetics. As shown in Fig. 1, such a structure consists of one or more layers of finite thickness, stacked along the "vertical" axis (z). Along the other two directions (in the "horizontal" x-y plane), the slab is effectively infinite in extent, but is patterned with structural elements such as pillars or holes, making it spatially non-uniform [1,2]. Slab structures form a basis for a wide range of important photonic devices, including photonic crystal slab lasers [3][4][5][6][7][8] and metasurfaces [9,10].
This problem is commonly tackled using standard finite element (FE) or finite difference (FD) discretization, with absorbing conditions or perfectly matched layers (PMLs) to simulate the infinite free-space region outside the slab. There also exist two special computational methods tailored to slabs with piecewise-uniform dielectric functions: Rigorous Coupled Wave Analysis (RCWA) [11,12] and the Vertical Mode Expansion Method (VMEM) [13,14]. In RCWA, 2D Fourier transforms are performed in each horizontal layer, and the solutions are matched along the layer boundaries. In VMEM, a mode expansion is performed on a set of vertical boundaries (e.g., cylinder or hole surfaces), and boundary integral equations are used to match the expansions in the horizontal plane. All of these methods are fairly efficient for periodic structures, where the computational domain can be reduced to a single horizontal lattice cell. (VMEM can also be applied to a standalone feature such as a single cylinder or hole [13].) The computational problem becomes substantially more difficult when there is no horizontal periodicity. In recent years, there has been increasing research interest in photonic structures that are deliberately patterned in an aperiodic manner, including random or amorphous patterns [15][16][17][18][19][20][21], quasicrystalline patterns [22][23][24][25][26][27], and topological defects [28,29]. Such structures have a variety of interesting properties; for instance, amorphous photonic structures have been found to exhibit highly isotropic spectral gaps [16], while slabs patterned with topological defects have been shown to support laser modes with inherent power-flow vorticity [29]. To model these various types of aperiodic slab structures, it is generally necessary to define a horizontal supercell, which in turn requires a large mesh (for FE and FD methods) or a high-order series expansion (for RCWA and VMEM). Such approaches fail to exploit a common feature in most of these examples, which is the existence of short-range regularity. Amorphous and quasicrystalline patterns, for instance, are often generated by taking one or more fixed elements-e.g., pillars or holes with a single fixed radius-and assigning them to aperiodic positions [16,22]. This paper explores a new computational method adapted to wave scattering problems in slab structures where the horizontal patterning has short-range regularity but no long-range periodicity, and hence no simple unit cell. We call this method "FLAME-slab" because its core component is high-order finite difference (FD) schemes generated by the Flexible Local Approximation MEthod (FLAME) [30,31]. FLAME replaces the local polynomial expansions of classical FD with more general, and often much more accurate, approximations by Trefftz functions, which by definition satisfy the governing equation of the problem and the applicable interface boundary conditions (b.c.)-in our case, Maxwell's equations and their associated b.c.
In FLAME-slab, as explained in detail in the following sections, the Trefftz basis functions are generated via small-scale simulations performed by an auxiliary method such as RCWA.
The key point is that the bases and the corresponding difference schemes are purely local, and can thus be computed on geometrically small but physically representative segments of the overall structure; we shall call these segments "Trefftz cells". (Of course, a global FD scheme for the whole structure must still be assembled and solved, but many standard numerical linear algebra techniques can be brought to bear on that part of the computational procedure.) Assuming the slab has the aforementioned short-range regularity, only a small number of basis functions need to be computed. Moreover, the basis functions can be cached and re-used for different slab structures possessing the same local parameters.
To illustrate the method, in Section II we formulate FLAME-slab for a 2D structure with a large number of randomly-positioned rectangular pillars, as depicted in Fig. 1(b).
In this case, the Trefftz basis functions can be generated from RCWA simulations of one pillar with a surrounding section of the substrate. Although this reference problem is 2D, FLAME-slab is intended to be generalizable to fully 3D geometries. In Section III, we show that the FLAME-slab results agree with RCWA calculations on the full structure. When the FLAME-slab solver re-uses basis functions, its runtime is around two orders of magnitude lower than RCWA alone.
We remark in passing that Trefftz basis functions can also be generated using other numerical methods, such as FEM or classical FD, as well as VMEM [13,14]. We have opted for RCWA because of its advantages: it does not require spatial discretization of the vertical direction (so long as the slabs are piecewise-uniform), handles the scattering boundary conditions exactly, without PML approximations, and is available as free software [12].

B. Formulation of the scattering problem
We consider a slab of a finite thickness d slab in free space, extending over 0 ≤ z ≤ d slab .
The slab has magnetic permeability µ = 1 everywhere, and can be characterized by a (possibly complex) dielectric permittivity (ω, r). In this paper, the problem is formulated in the frequency domain, and dependence of parameters on the frequency ω will not henceforth be explicitly indicated. No theoretical limitations on the position dependence of (r) are imposed, although in practice (r) is typically piecewise-constant, representing a modest number of different constituent materials in an actual device.
The present paper demonstrates a proof-of-concept implementation of FLAME-slab for an effectively-2D geometry. As shown in Fig. 1(b), the slab is assumed to be translationally invariant along y. The computational domain is assumed to be a supercell, with quasiperiodic boundary conditions relating the fields on the left (x = 0) and right (x = L x ) sides. These conditions do not necessarily reflect an actual periodicity of the sample, but are imposed by necessity, to mimic the field behavior in an infinite (or very large) structure with as few numerical artifacts as possible.
The electromagnetic field in and around the slab is governed by Maxwell's equations, where k 0 = ω/c is the wavenumber in free space. The exp(−iωt) phasor convention is adopted. We consider the s-mode, with a one-component electric field E = Eŷ and a twocomponent magnetic field H = H xx + H zẑ . Further, we assume that the slab is illuminated by a plane wave incident from z < 0: where θ in is the angle of incidence. The quasi-periodic boundary conditions on the left and right sides of the computational domain are Finally, we introduce the usual splitting of the total fields into the incident and scattered parts. The formal definition of the scattered field is E s is subject to the standard radiation boundary conditions ensuring that E s propagates away from the slab (i.e., toward z → +∞ for z > d slab and toward z → −∞ for z < 0).

C. Rigorous Coupled Wave Analysis (RCWA)
In this section, we briefly review the Rigorous Coupled Wave Analysis (RCWA) algorithm for solving the electromagnetic wave scattering problem in slabs [11,12]. Our implementation of FLAME-slab uses the RCWA solver S 4 (the Stanford Stratified Structure Solver) [12] to calculate Trefftz basis functions, and also to provide comparisons for the numerical results.
The slab (see Fig. 1) is divided into layers in the vertical (z) direction. It is assumed that the dielectric function ε(x, y, z) is periodic in the x-y plane (if necessary, by defining a supercell), and piecewise-uniform in the z direction. We let k denote the projection of the incident wave-vector, k in , onto the plane. In-plane Fourier decomposition of the electromagnetic field in each layer [44] is then carried out. The electric field is expanded in the form where r = (x, y) denotes in-plane coordinates, {G} is the set of reciprocal lattice vectors, E G is a Fourier expansion coefficient, and q z is the Bloch wavenumber in the z direction. The magnetic field is expanded similarly. The in-plane Fourier series are truncated by specifying an upper bound for |G|. The number of terms retained is denoted by N G , so that N G ∼ |G| 2 for a 3D problem (with 2D periodicity), and N G ∼ |G| for 2D (with 1D periodicity).
The Fourier harmonics are then matched across adjacent layers via the transfer matrix method [45]. After repeatedly applying the procedure to each layer, one obtains a transfer matrix relation between the fields at the outermost layers of the slab, and hence the reflection and transmission coefficients. From these, all other quantities of interest can be retrieved, including the fields within the slab as well as the far field. For details, see [11,12].
For our present purposes, the key point to note is that the runtime of the RCWA calculation depends principally on the integer N G , the number of reciprocal lattice vectors included in the in-plane Fourier expansion (5). The method involves running an eigensolver on a full matrix with O(N G ) rows/columns. Hence, the scalings of the runtime T RCWA and the memory usage M RCWA are We have verified that the performance of the S 4 solver [12] is consistent with these scalings.

D. FLAME
In this section, we give a brief general description of FLAME [30,31], and explain how it can be specialized to the slab problem under consideration. This exposition follows previous publications on this subject, especially [41,42].
Consider a differential equation of the form where D(r) is a differential operator, and u(r) is either a scalar or vector field. Mathematically, Eq. (7) should be understood in weak form: in physical terms, it includes proper interface boundary conditions in addition to the differential equation itself.
The general setup for FLAME has two ingredients: (i) a set of n local Trefftz functions solving Eq. (7)  Within each patch, the solution u is approximated locally using a linear combination of Trefftz functions ϕ α (α = 1, 2, . . . , n) : where c ∈ C n is a coefficient vector and ϕ is a vector of basis functions (both generally complex). For each patch (i.e. for each grid "molecule") we seek an FD equation of the where s = (s 1 , s 2 , . . . , s m ) T is a vector of complex coefficients (a "scheme") to be determined. In the simplest version of FLAME, the scheme is required to be exact for any linear combination (8) of basis functions (Least-squares schemes can also be useful [46,47]). After straightforward algebra, it can be shown that [30,31] The construction of the scheme is purely local, which is consistent with the local nature of the differential operator D in Eq. (7); therefore the choices of Trefftz bases and DoF in different patches are independent of one another. Specifically, in this paper, the patches centered near the middle of the slab use different DoF from those in the boundary patches.
In Section II, we consider a 2D geometry with rectangular pillars distributed on a slab substrate, shown in Fig. 1(b). We will define three types of patches; one has DoF corresponding to nodal values of the E-field, while the other two types have DoF consisting of nodal values of the E-field, plus tangential H field components at the boundary nodes. The assignment of these DoF is further elaborated on in Section II A.
Each Trefftz basis only needs to be local, i.e. valid over a given patch. Therefore, it can be generated by solving a set of relatively small auxiliary problems. For instance, for the slab considered in Section II, the local problem can be geometrically limited to "Trefftz cells" consisting of a single pillar with a segment of the surrounding substrate. We obtain basis functions by solving the wave scattering problem in the Trefftz cells, with different illumination conditions. A large part of the computational savings comes from the fact that the local problem is much smaller in size than the global one. This procedure is described in detail in Section II B. In this paper, we use RCWA to solve the local problem, but in principle any suitable method could be used.
Once a Trefftz basis and DoF have been determined for any given patch, Eq. (10) yields the local difference scheme immediately. In our setup, the scheme has nine terms in each patch, both inside the slab and at its boundary. (It is tacitly assumed that the null space of N T is one-dimensional.) The FLAME block of the global FD matrix is then assembled from the local schemes in a standard way, as is done in all FD algorithms; see Section II D.

A. Discretization and Trefftz patches
Our FLAME-slab implementation applies to the 2D slab geometry described in Section I B and depicted in Fig. 1(b). We adopt the convenient discretization shown in Fig. 2 (a). To demonstrate that sufficient accuracy can be achieved in FLAME-slab, we have deliberately chosen only three grid layers (N z = 3) across the thickness of the slab, including the two outer layers. These outer layers of nodes lie in the air regions slightly above/below the structure. In the x direction, the 2D computational domain (or "supercell") of total length L x contains N x grid lines, which are equally-spaced for simplicity. Importantly, the grid need not resolve the fine geometric features of the structure, as the relevant information will be contained in the Trefftz bases. In particular, even though we only use three layers of grid points, the slab could consist of multiple layers of dielectric material, and the grid lines need not coincide with any material interfaces. The total number of unknowns in the system of equations to be solved is N x (N z + 2); in the present setup, with N z = 3, this number is 5N x . Of these unknowns, N x · N z are the values of the E field at all nodes, and the remaining 2N x are the values of H x at the boundary nodes. Correspondingly, N x · N z equations of the system are the Trefftz-FLAME difference schemes described above. The remaining 2N x equations come from the radiation boundary conditions, as explained in Section II C.
In matrix form, the FD-FLAME equations can be written as where A FLAME is a rectangular N x N z × N x (N z + 2) matrix and ψ is a vector containing all N x (N z + 2) DoF (nodal values of the fields).

B. Generation of Trefftz basis functions
As already noted, the computational savings in FLAME-slab come from the fact that Trefftz functions in FD-FLAME only need to be local, and can thus be computed over a relatively small representative segment of the structure. (In practice, they can be precomputed and stored.) For the present 2D slab geometry, we define a "Trefftz cell" consisting  A natural set of illumination conditions comes from plane waves impinging on the segment of the slab at different angles, as schematically shown in Fig. 3. The scattered wave is subject to the standard radiation conditions above and below the slab. One also needs to impose some boundary conditions on the left and right sides of the segment. Extensive computational evidence in photonics, molecular dynamics, and other areas-where very large structures are modeled via a finite representative sample-has shown that the least intrusive conditions are (quasi-)periodic ones. In our case, such conditions read where x is the "horizontal" coordinate relative to the pillar, and θ i is the angle of illumination. Due to these quasi-periodic conditions, the segment of length L can be viewed as the unit cell of a fictitious photonic crystal. The fact that the actual amorphous structure under investigation does not have this periodicity is unimportant because the Trefftz basis approximates the solution only locally, over a small patch, and its long-range behavior is inconsequential. The fictitious period L is much smaller than the computational domain, but its dielectric function matches the actual structure in the vicinity of a given pillar.
In our examples, nine-point FD stencils are used, and eight Trefftz basis functions are needed to produce a valid FLAME scheme. There are several ways to accomplish this. We settled on the use of two different lengths L (12) for the local lattice, L 1 and L 2 , and four angles of illumination {θ 1 , . . . , θ 4 } for each of these lengths. The 2 × 4 = 8 solutions (Trefftz functions) are precomputed with the RCWA solver S 4 [12], and stored. Using these, we can apply the FLAME equation (10) to all grid stencils, as described in Section II A, thereby populating the 3N x rows of the A FLAME matrix (11).
We consider structures with short-range regularity, in the sense that all the pillars are geometrically identical; only their placement is random. Many photonic structures of interest, including amorphous and quasicrystalline structures [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29], have similar forms of shortrange regularity. This feature is exploited here to minimize the number of Trefftz functions that need to be pre-computed; in our examples, this number is just eight for fixed-frequency simulations. Each solution is calculated with a high spatial resolution, cached, and later retrieved during the actual construction of the Trefftz-FLAME schemes. This procedure is even more appealing for the simulation of ensembles of multiple independent realizations of the slab geometry (e.g., different samples of an amorphous lattice with the same elementary structural parameters), which is an important practical case. Then the cache can be reused in all of the calculations. Note also that for certain quasicrystalline designs, it may be advantageous to choose Trefftz cells with more than one pillar; for example, for Fibonacci quasicrystals [26,27], we might use Trefftz cells corresponding to two different two-pillar segments.

C. Radiation conditions
We apply the standard splitting of the electromagnetic field into incident and scattered parts; this is needed because radiation conditions apply to the scattered component: (These equations can be viewed as a formal definition of the scattered field.) Applying this decomposition to Eq. (11) yields where ψ in and ψ s are the decompositions of the vector of field components into incident and scattered parts.
Radiation conditions can be imposed via a Fourier decomposition in the empty semiinfinite strips above and below the slab. This is facilitated by the fact that the outer grid layers have been placed inside the free-space regions rather than in the slab itself (see Section II A). Maxwell's equations (1) relate H s (x, z) and E s (x, z) as where H s refers to the tangential (x) component of the magnetic field, while E s is the electric field of the s-mode under consideration. For simplicity of presentation, we will describe the case of normal incidence, θ in = 0; then the quasi-periodic boundary conditions turn into periodic ones. The Fourier expansion of the scattered field in the free space regions is where n runs over integer values and the horizontal wavenumbers k nx = 2πn/L x .
Combining Eqs. (15) and Eq. (16) gives From the dispersion relation in free space, k nx and k nz are related by The choice of the ± sign depends on which semi-infinite strip the expression applies to.
Referring to the layout of Fig. 2, the plus sign applies for z = z + , and minus for z = z − .
These choices ensure that the scattered waves are outgoing from the slab.
We can now evaluate (16) and (17) along each outer layer of grid points, to produce relations between the electric and magnetic field components: Here, k nz uses the ± sign in Eq. (18) corresponding to z ± . The notation FT m {f n } denotes D. Assembling the finite-difference system Equation (19) contributes N x equations for each outer layer. Combining this with Eq. (14) yields   A FLAME Here, full sub-matrix representing the boundary conditions (19). The vector ψ in contains the field components of the incident wave, evaluated at the grid points. Figure 4 illustrates the structure of the sub-matrices A FLAME and A BC , for the test case where N z = 3. The sub-matrix A FLAME contains nine nonzero entries per row, based on the DoF assignment described in Section II A. The sub-matrix A BC has two full N x × N x subblocks, consisting of the Fourier components in Eq. (19).
Since A BC is full, its explicit calculation and direct solution of Eq. (21) will require an overall runtime of O(N 3 x ). A speedup can be achieved if, instead, an iterative solver such as GMRES [48] is applied to Eq. (21). Each iteration involves a matrix-vector product which could be computed on the fly with Fast Fourier Transforms for the dense blocks. The computational cost will be O(N iter N 2 x log N x ), where N iter is the number of iterations. With a right choice of preconditioners, this could scale better than O(N 3 x ).

A. Test problems
Our main test case is depicted in Fig. 1(b). It consists of 10 dielectric pillars placed randomly on a flat dielectric substrate. Except where stated otherwise, we assume a normally incident input wave whose frequency is f = 0.25. We adopt a system of units [12] where the speed of light (c), vacuum permeability (µ 0 ) and vacuum permittivity (ε 0 ) are all normalized to unity; this implies that f = 1/λ, where λ is the free-space wavelength, i.e. the frequency has units of inverse length. The normalized pillar height and substrate thickness are both d = 0.25, and the pillar width is w = 0.2. The pillars and substrate have the same dielectric The horizontal grid layers are placed at z + = 0, z 0 = 0.26, and z − = 0.5. For the Trefftz cells (Section II B), we take L 1 = 1.73 and L 2 = 4.12, along with the incidence angles θ 1 = 36 • , θ 2 = 72 • , θ 3 = 108 • , and θ 4 = 144 • , measured anti-clockwise from the +x axis.
The accuracy of the FLAME-slab method depends on three adjustable parameters, which we denote by N x , N T , and N G . The first, N x , is the number of nodes in the x direction (see Section II A). To make this physically meaningful, we express it as the ratio N x w/L x , the number of nodes per pillar width. The second parameter, N T , refers to the number of nodes in the x direction for the Trefftz cells. We normalize this parameter as the ratio N T w/L 1 , the number of nodes per pillar width for the smaller of the two Trefftz cells. The third and final parameter, N G , controls the accuracy of the RCWA method. As described in Section I C, this corresponds to the number of reciprocal lattice vectors included in the horizontal Fourier expansion.

B. Electromagnetic field and reflectance calculations
As described in Section II A, each FLAME-slab calculation yields the scattered electric field E y,s along the three layers of grid points, denoted by {E + s , E 0 s , E − s }, as well as the scattered magnetic field H x,s along the two outer layers, denoted by {H + s , H − s }. From these field components, we can extract two other quantities that are of key interest in electromagnetic wave scattering problems: the reflectance R and transmittance T . For s-polarized fields, the z-component of the time-averaged Poynting vector is which involves precisely the field variables produced by FLAME-slab. To obtain R, we average S z over the N x grid points along z = z − , using only the scattered fields. The calculation of T uses the total fields along z = z + , the side of the slab opposite to the incident wave.   R + T = 1 to high precision. At higher frequencies, f 0.5, numerical errors become evident. A similar frequency-dependence of the numerical errors has also been observed in the field components. It is to be expected that higher frequencies require finer grids and other parameter adjustments.

C. Error Analysis
As mentioned in Section III A, there are three main sources of error in FLAME-slab: 1. Even though FLAME schemes are qualitatively more accurate than the classical ones, FD-FLAME is still subject to discretization errors arising from the finite size of the Trefftz basis and DoF. This error depends on the dimensionless parameter N x w/L x .
2. The RCWA solutions are calculated on a grid, and interpolation is used to evaluate the Trefftz basis functions in each local stencil; this is controlled by N T w/L 1 .
3. The basis functions are calculated using RCWA, which itself has finite accuracy; this is controlled by the paremeter N G .
To quantify how the accuracy of FLAME-slab depends on these parameters, we perform a comparison to the reference solution noted above. We define a normalized "solution error", Here, ψ Fs i denotes components of the FLAME-slab solution vector (see Section II), consisting of the E fields on all nodes and H x field components along the boundary nodes; and ψ ref The specific choices of N x , N T , and N G determine which of the three parameters dominate the solution error δ at any one time. For instance, Fig. 6(a) shows that for a fixed value of N x , δ decreases rapidly and then saturates. This saturated value decreases as we increase N G from 51 to 97, but stops falling further as we increase N G to 150. Next, suppose that we fix N T w/L 1 = 160 (well within the region where δ is saturated), and increase N x . As shown in Fig. 6(b), this causes δ to fall rapidly, before saturating again.

D. Runtime and memory performance
We have profiled the runtime and memory performance for FLAME-slab calculations against calculations performed purely using RCWA. For RCWA, we observe that an increase in N G reduces the solution error, while increasing the runtime (as N 3 G ). For FLAME-slab, the total runtime can be competitive with or superior to the runtime of an RCWA calculation at a comparable level of solution error, so long as the parameters {N x , N T , N G } are chosen appropriately. In Fig. 7, for example, the FLAME-slab total runtime is shown to be ≈ 70% of the RCWA runtime with comparable solution error, for the 10-pillar case; for the 20-pillar case, it is ≈ 60%. FLAME-slab acquires a dramatic advantage if the Trefftz basis functions can be saved and re-used. Most of the total runtime of FLAME-slab is spent on the initial calculation of the Trefftz basis functions. As we have previously noted, numerical studies of aperiodic photonic structures are often concerned with ensembles of independent samples with the same regular elements, such as amorphous or quasicrystalline lattices built out of identical pillars. With FLAME-slab, the Trefftz basis functions can be computed just once and reused; when we use RCWA or other traditional methods, by comparison, the calculation must be re-done entirely for each sample.
To illustrate the magnitude of this advantage, Fig. 7 shows also the FLAME-slab runtimes excluding the time taken to pre-compute the Trefftz basis functions. These are around 2 orders of magnitude smaller than RCWA runtimes with comparable solution error.
This numerical data was obtained using a direct solver, which required all matrix elements to be computed explicitly. As discussed in Section II D, it should be possible to achieve an additional speedup by avoiding direct construction of the A BC matrix, and instead using an iterative solver with matrix-vector products evaluated by Fast Fourier Transforms.
The memory consumption of FLAME-slab is significantly lower than that of full-scale RCWA calculations. For N G = 400, the S 4 program uses 198 MB of memory. As seen in Thus, the overall memory usage of FLAME-slab in this example is an order of magnitude lower than that of S 4 .

IV. CONCLUSION
We have devised a numerical method, FLAME-slab, which is suitable for solving the wave scattering problem in photonic slabs lacking in-plane periodicity. FLAME [30] is used to generate finite-difference stencils from a set of basis functions. In this case, the basis functions are computed by running a subroutine-the S 4 RCWA solver [12]-on a set of local problems defined on small-scale "Trefftz cells".
One advantage of FLAME-slab over traditional finite-difference and finite-element methods is that a relatively coarse grid can be used without sacrificing accuracy. In the present implementation, the grid consists of only three horizontal layers (z + , z 0 , and z − ). Another key advantage is that for slab structures with short-range regularities, such as pillars that are placed randomly but have identical geometries and materials, the basis functions only have to be evaluated once. A further improvement, which we have not yet extensively explored, is to use an iterative solver.
In future work, we seek to implement the scheme for fully 3D slab structures. We will also investigate how best to treat the outgoing boundary conditions; it may be possible to use boundary element techniques, perfectly-matched layers, or absorbing boundary conditions instead of FFTs. It will also be interesting to explore methods other than RCWA to generate the basis functions, including the Vertical Mode Expansion Method [13].
V. FUNDING