1. Introduction

1.1. 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–8] and metasurfaces [9, 10].

 

Fig. 1 (a) Schematic of electromagnetic wave scattering from a 3D slab. The slab consists of one or more dielectric layers stacked in the z direction, and is patterned (e.g., with pillars) along the x-y plane. A plane wave is incident from the −z side, with a wave vector k_in. (b) The 2D version of the problem, where the dielectric function does not vary with y, and the wave is incident along the x-z plane.

Download Full Size | PDF

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 applied 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–21], quasicrystalline patterns [22–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.

Previously, Trefftz approximations have been successfully used in classical variational contexts (FE and domain decomposition methods) [32–34], in Discontinuous Galerkin methods [35–38], and in FLAME [30, 31, 39–42]; see also review [43] by Hiptmair et al.

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 2 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 3, 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].

1.2. 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 quasi-periodic 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,

1.3. 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⁴ (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

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 $\mathcal{O}(N_G)$ rows/columns. Hence, the scalings of the runtime T_RCWA and the memory usage M_RCWA are

1.4. 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

The general setup for FLAME has two ingredients: (i) a set of n local Trefftz functions solving Eq. (7) within small local subdomains, or “patches”, of the computational domain (each patch contains a grid “molecule”—a set of geometric entities such as grid nodes—on which the difference scheme is to be formed), and (ii) a set of m degrees of freedom (DoF), also defined for each patch. These DoF are, by definition, linear functionals, l_β(u) (β = 1, 2, …, m), each mapping the field u to a number. The term “DoF” is used here in its mathematical sense, rather than that of kinematics or statistical physics; in the simplest case, if u is a scalar field, we might define l_β(u) ≡ u(r_β), where r₁,…,r_m are a set of grid nodes. Other examples of DoF include the circulation of a field over a given path, the flux over a given surface, the nodal value of any partial derivative of a field component, and so on.

Within each patch, the solution u is approximated locally using a linear combination of Trefftz functions φ_α (α = 1, 2,…,n):

In Section 2, 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 2.1.

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 2, 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 2.2. 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 2.4.

2. Finite-difference schemes for slab geometry

2.1. Discretization and Trefftz patches

Our FLAME-slab implementation applies to the 2D slab geometry described in Section 1.2 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.

 

Fig. 2 (a) A coarse N_x × N_z grid for the slab structure, where N_z = 3. (b) Assignment of DoF for patches centered on column m. Left: patch centered at z⁰, whose DoF are the y-components of the electric field in the nine nodes. Center and right: boundary patches, whose DoF are the electric fields at the six nodes and the x-component of the magnetic field at the boundary nodes.

Download Full Size | PDF

Figure 2(b) shows how the DoF are assigned in each patch. There are three types of patches. The first type contains 3 × 3 = 9 nodes on adjacent grid lines, and we take its DoF to be the nine nodal values of the total E field, $E_{m-1,m,m+1}^{-,0,+}$. (The splitting of the field into “incident” and “scattered” parts is only needed later, to impose the radiation boundary conditions; see Section 2.3.) The second and third types of patch contain 2 × 3 = 6 nodes on the upper or lower boundary; for each, six DoF are the respective values of the E field, and three more DoF are the values of the magnetic field component H_x, tangential to the slab, at the boundary nodes m − 1, m, m + 1.

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 2.3.

In matrix form, the FD-FLAME equations can be written as

2.2. 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 pre-computed and stored.) For the present 2D slab geometry, we define a “Trefftz cell” consisting of a single pillar, with some surrounding substrate of a length L. By solving Maxwell’s equations in the Trefftz cell for different illumination and boundary conditions, we produce the desired Trefftz basis set.

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

 

Fig. 3 Procedure for generating Trefftz basis functions. (a)–(b) Schematic of “Trefftz cells” defined as unit cells of periodic slabs with one pillar per period. Two periods are used, L₁ and L₂; for each, the RCWA solver produces four Trefftz functions corresponding to angles of illumination θ₁,…, θ₄. (c) For a given patch centered at grid coordinate x_m, one uses the eight functions precomputed in (a, b) according to the position x′ of the nearest pillar.

Download Full Size | PDF

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₁ and L₂, and four angles of illumination {θ₁,…,θ₄} for each of these lengths. The 2×4 = 8 solutions (Trefftz functions) are precomputed with the RCWA solver S⁴ [12], and stored. Using these, we can apply the FLAME equation (10) to all grid stencils, as described in Section 2.1, 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–29], have similar forms of short-range 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 re-used 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.

2.3. 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:

Radiation conditions can be imposed via a Fourier decomposition in the empty semi-infinite 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 2.1). Maxwell’s equations (1) relate H_s(x, z) and E_s(x, z) as

Combining Eqs. (15) and Eq. (16) gives

We can now evaluate (16) and (17) along each outer layer of grid points, to produce relations between the electric and magnetic field components:

2.4. Assembling the finite-difference system

Equation (19) contributes N_x equations for each outer layer. Combining this with Eq. (14) yields

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 2.1. The sub-matrix A_BC has two full N_x × N_x subblocks, consisting of the Fourier components in Eq. (19).

 

Fig. 4 Structure of the matrix A for N_z = 3. Non-zero entries are indicated in red.

Download Full Size | PDF

Since A_BC is full, its explicit calculation and direct solution of Eq. (21) will require an overall runtime of $\mathcal{O}(N_x^3)$. 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 $\mathcal{O}(N_{\text{iter}}{N}_x^2\log N_x)$, where N_iter is the number of iterations. With a right choice of preconditioners, this could scale better than $\mathcal{O}(N_x^3)$.

3. Results

3.1. 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 (µ₀) and vacuum permittivity (ε₀) are all normalized to unity; this implies that f = 1/λ, where λ is the free-space wavelength, i.e. 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 constant ε = 12, and are surrounded by air. The pillars are placed at aperiodic positions on the substrate. The entire computational domain has length L_x = 14 in the x direction.

The horizontal grid layers are placed at z₊ = 0, z₀ = 0.26, and z₋ = 0.5. For the Trefftz cells (Section 2.2), we take L₁ = 1.73 and L₂ = 4.12, along with the incidence angles θ₁ = 36°, θ₂ = 72°, θ₃ = 108°, and θ₄ = 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 2.1). To make this physically meaningful, we express it as the ratio N_xw/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_Tw/L₁, 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 1.3, this corresponds to the number of reciprocal lattice vectors included in the horizontal Fourier expansion.

3.2. Electromagnetic field and reflectance calculations

As described in Section 2.1, each FLAME-slab calculation yields the scattered electric field E_y,s along the three layers of grid points, denoted by $\{E_s^{+},E_s^0,E_s^{-}\}$, as well as the scattered magnetic field H_x,s along the two outer layers, denoted by $\{H_s^{+},H_s^0\}$.

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

Figure 5(a)–(b) shows a typical set of electric and magnetic fields calculated for the test structure. Here, we plot $\mathrm{Re}[E_s^0]$ and $\mathrm{Re}[H_s^{-}]$ versus x, comparing the FLAME-slab results to a reference solution obtained by applying RCWA to the entire structure. For the FLAME-slab calculation, we take N_xw/L_x = 4.30, N_Tw/L₁ = 160, and N_G = 150. For the RCWA reference solution, we take an “overkill’ value $N_G^{\text{ref}}=1000$ to ensure high accuracy for the big global structure. From the figure, we observe that the FLAME-slab solution is in excellent agreement with the reference one. The results for the other field components show similar agreement.

 

Fig. 5 (a) Real part of the electric field E_s versus horizontal position x, along the middle layer of grid points at z = z₀. (b) Real part of the magnetic field H_s versus x, along the outer layer of grid points at z = z₋. In both cases, blue dots show FLAME-slab results calculated using N_xw/L_x = 1.44, N_Tw/L₁ = 160, and N_G = 150; red curves show reference solutions calculated using RCWA with $N_G^{\text{ref}}=1000$. The frequency is f = 0.25, and all other parameters are as stated in the main text. (c) Reflectance R versus frequency. Blue dots show FLAME-slab results, obtained using N_xw/L_x = 1.44, N_Tw/L₁ = 160, and N_G = 97; the red curve shows RCWA results calculated with $N_G^{\text{ref}}=1000$.

Download Full Size | PDF

Figure 5(c) plots the reflectance R versus frequency f. For frequencies f ≲ 0.5, the FLAME-slab results are again in close agreement with the reference solution. The transmittance is not shown here, but we verified that it satisfies the flux conservation condition 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.

3.3. Error analysis

As mentioned in Section 3.1, there are three main sources of error in FLAME-slab:

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_Tw/L₁ = 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.

 

Fig. 6 (a) Solution error δ versus the discretization of the local lattice used for computing the Trefftz functions, N_Tw/L₁. (b) Solution error δ versus the discretization of the global lattice, N_xw/L_x. The frequency is f = 0.25, and the reference solution is computed by RCWA using $N_G^{\text{ref}}=1000$.

Download Full Size | PDF

3.4. Runtime and memory performance

We have profiled the runtime and memory performance for FLAME-slab calculations against calculations performed purely using RCWA.

Figure 7 plots the solution error versus runtime for FLAME-slab and RCWA calculations. Two test cases are studied: (i) a slab with 10 pillars and L_x = 14, identical to the structure tested in Sections 3.2–3.3, and (ii) a larger slab with 20 pillars and L_x = 28. In both cases, the pillars are placed randomly, without overlap, along the length of the slab. For the first case, the reference RCWA solution used to determine the solution error was computed with $N_G^{\text{ref}}=1000$; for the second case, the reference RCWA solution was computed with $N_G^{\text{ref}}=1200$, the maximum possible before running out of memory on the desktop computer performing these tests.

 

Fig. 7 Solution error versus runtime for FLAME-slab and RCWA solutions, with two differently-sized test cases: (i) 10 pillars and L_x = 14, and (ii) 20 pillars and L_x = 28. The FLAME-slab solutions were calculated with N_xw/L_x = 4.30 and N_Tw/L₁ = 160; and N_G = 97 for case (i), and N_G = 120 for case (ii). The RCWA solutions were calculated with N_G = 200, 300, 400, 500. The reference solutions were calculated by RCWA using N_G = 1000 for case (i), and N_G = 1200 for case (ii).

Download Full Size | PDF

For RCWA, we observe that an increase in N_G reduces the solution error, while increasing the runtime (as $N_G^3$). 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 re-used; 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 2.4, 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⁴ program uses 198 MB of memory. As seen in Fig. 7, FLAME-slab achieves a similar solution error with the parameters N_xw/L_x = 4.30, N_Tw/L₁ = 160, and N_G = 97. The memory usage for FLAME-slab, with these parameters, consists of: (i) pre-computation of the Trefftz functions using RCWA, which takes ≈ 11.6 MB of memory; (ii) storage of the Trefftz functions, which takes ≈ 1.5 MB; (iii) storage of the A matrix, which occupies ≈ 4.4 MB of memory; and (iv) the memory usage of the direct solver, which we estimated as ≈ 6.6 MB based on the sparsity of the LU decomposition. Thus, the overall memory usage of FLAME-slab in this example is an order of magnitude lower than that of S⁴.

4. 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⁴ 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. Furthermore, the grid lines do not have to conform to the geometric features of the structure, since the relevant physical information is contained in the Trefftz bases rather than in the layout of the grid. In the present implementation, there are only three grid lines across the thickness of the slab. 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].