1. Introduction

Nanophotonic structures often consist of large collections of repeated meta-atoms as in, for example, many functional photonic crystals or metamaterials [1–11]. Furthermore, the individual meta-atoms often contain fine, subwavelength features that makes simulating large collections of these meta-atoms challenging. In simulating these structures, finding a way to exploit the high degree of repeated structure could provide significant efficiency improvements.

A simple and natural way to exploit the repetition of meta-atoms in such structures is to use a non-overlapping domain decomposition approach, for example a Schur complement domain decomposition [12]. Domain decomposition has been investigated for nanophotonics extensively in the context of finite element methods using sophisticated iterative substructuring techniques, for example in finite element tearing and interconnect (FETI) methods [13, 16]. Additionally, variations focusing on solution of scattering problems using surface integral methods as in [14, 15] have also been studied. However, given the simplicity of finite difference methods relative to finite element methods, investigating simple domain decomposition methods for the former is worthwhile, and we will show that applying a Schur complement domain decomposition to domains consisting of many repeated meta-atoms proves quite advantageous.

In the domain decomposition approach, one solves Maxwell’s equations on a large domain by breaking that domain into smaller subdomains, and solving relatively small problems associated with each. In the Schur complement domain decomposition, one transforms the original problem on the large domain by first solving Maxwell’s equations on a number of smaller subdomains in isolation. One then constructs a secondary, usually smaller, problem on the interfaces between subdomains to account for coupling between the subdomains. One uses the solution to this interfacial problem, combined with the solutions to the subdomain problems, to recover the solution on the complete domain. To take advantage of repeated meta-atoms, one can solve the subproblem associated with a repeated meta-atom once for all meta-atoms of the same kind, and reuse the solution at negligible cost wherever necessary.

The domain decomposition approach thus reduces the original problem to a problem with many fewer degrees of freedom, namely those associated with the interfaces between the subdomains. Such an interfacial problem is cheap to form and small in dimension relative to the original system. Also, the domain decomposition approach is useful in design optimization problems since one can solve an associated subproblem on subdomains that do not change during the optimization, and then perform calculations only on the subdomain involving the optimization region and the interfaces [12, 17–19].

Furthermore, iterative solution of Maxwell’s equations using the standard finite-difference frequency-domain method can suffer from indefiniteness and poor conditioning of the matrix resulting from discretization, even when standard preconditioning methods are used [20]. We show that applying a Schur complement domain decomposition can mitigate such problems, leading to substantial improvements in convergence rate. Additionally, we show that the ratio of source wavelength to subdomain size is an important consideration in achieving such improvement.

The rest of this paper proceeds as follows. Section 2 describes the Schur complement domain decomposition via an example, using it to accelerate the convergence of iterative solution of a finite-difference frequency-domain (FDFD) discretization of Maxwell’s equations. Section 3 provides several numerical experiments showing improvements in convergence when applying the Schur complement domain decomposition method. Section 4 shows how the numerical results in Sec. 3 can be heuristically understood by analyzing the conditioning of the matrices that arise from the domain decomposition method. Section 5 addresses memory and wall-clock time requirements of the domain decomposition method, and briefly comments on the use of subgridding to further enhance the method, particularly in three dimensions. Finally Sec. 6 provides some concluding remarks.

2. Formulation of the Schur complement domain decomposition

For most of this paper, we will consider two-dimensional structures such as the structure shown in Fig. 1. The structure consists of a collection of metallic or dielectric objects distributed on a square lattice with lattice constant a. The objects and the bounding cell will be referred to as meta-atoms. The space between the objects in each meta-atom is air. Either a periodic boundary condition or a perfectly matched layer (PML) boundary condition surrounds the domain.

 

Fig. 1 Two-dimensional simulation domain consisting of an 8 × 8 grid of square meta-atoms. Each meta-atom contains a dielectric inclusion shown as a gray shape. A point source is placed in the center. The dotted lines denote the interfaces between meta-atoms. The lightly shaded region on the borders of the domain represents the perfectly matched layer boundary condition.

Download Full Size | PDF

Figure 1 shows an exemplary structure on an 8a × 8a lattice. For illustration purposes we consider three types of objects: square (with a side length of 0.5a), isosceles right triangle (with a base length of 0.33a), and circle (with a radius of 0.25a). Each such object resides at the center of a meta-atom. We assume that the objects all have the same dielectric constant ϵ. Throughout this paper, we consider two cases, the dielectric case, where ϵ = 12 approximates the dielectric constant of Si and GaAs in the infrared wavelength range, and the metallic case, where ϵ = −3 − 0.3i is the complex dielectric constant for a metal with some loss. The choice of the metallic dielectric constant was made so that the resolution of the grid picked for the dielectric case would be sufficient to resolve the surface waves that would appear on the metal. The structure shown here does not have any apparent optical functionality, and serves only to demonstrate the algorithms. We note, however, that many published works on nanophotonic structures, including photonic crystals, meta-materials, and meta-surfaces, consider similar structures, where a set of dielectric or metallic objects (not necessarily identical) are distributed on a regular lattice [1–11].

We solve for the field distribution inside the domain induced by a point source at the center of the domain. For concreteness, we consider the TE polarization where the magnetic field H has only one component perpendicular to the plane [9], which we designate the z-component. The governing equation is then:

A standard way to solve Eq. (1) is the finite-difference frequency-domain method (FDFD). In this method, one discretizes the H_z field on a regular grid. In our simulation the grid spacing is 0.012a. Approximating Eq. (1) with finite differences on that grid yields a linear system of the form Ax = b. This linear system, when describing a domain consisting of a collection of meta-atoms like that of Fig. 1, can be written:

The vector x₀ corresponds to the field components that lie at the interfaces between different meta-atoms, as indicated by the dashed lines in Fig. 1. x_i for i ≠ 0 corresponds to the field components in the interior of the ith meta-atom. Each diagonal block in A corresponds to a restriction of the approximation to Eq. (1) to a subdomain and terminated with a Dirichlet boundary condition. The block A₀₀ describes the coupling of field components within the union of the boundaries, or interfaces, of all the meta-atoms, which we call the interface. Each block D_ii describes the couplings between field components within the interior of each meta-atom. The blocks B_0i and C_i0 describes the couplings between the field components at the boundary of the interior of the ith meta-atom with the field components at the interface of the ith meta-atom. Finally, b_i contains information about the sources within the ith meta-atom.

To solve Eq. (2) via the domain decomposition approach, we solve for each x_i in terms of x₀:

Then, we substitute every such expression for x_i into the equation obtained from the first row of A in Eq. (2), yielding a reduced system expressed only in terms of x₀:

S is the Schur complement of the block diagonal matrix D = diag(D₁₁, D₂₂, …, D_nn) in A. b′ will denote the modified source term on the right hand side of Eq. (4). This linear system contains only the degrees of freedom on the interfaces shown in Fig. 1. In this paper, we refer to the system in Eq. (4) as the reduced system, and the system in Eq. (2) as the unreduced system. There is a natural correspondence between solving these two systems whenever D is not singular since one can recover the full solution on the total domain after solving Eq. (4) by substituting the solution x₀ into Eq. (3). Numerically, however, the reduced system S is much smaller in dimension than A, often up to an order of magnitude or more, meaning it could be much more efficient to solve Eq. (4) in practice.

To explicitly form S in Eq. (4), one must evaluate the action of $D_{ii}^{-1}$ on the C_0i blocks, a potentially computationally expensive step. But where many of the meta-atoms are identical, and consequently many of the sub matrices D_ii are identical, explicitly forming S requires only evaluating a suitable factorization (to calculate the action of the inverse of D_ii) of every distinct D_ii once and reusing them at negligible or much lower expense to calculate $D_{ii}^{-1}{C}_{i0}{b}_i$ and $B_{0i}{D}_{ii}^{-1}{C}_{i0}$. For general systems with repeated meta-atoms, the interface degrees of freedom in x₀ should be selected to minimize the number of distinct subdomains, and hence minimize the number of distinct factorizations of D_ii. In addition, there is a memory cost since S is a denser matrix than A due to fill-in from $B_{0i}{D}_{ii}^{-1}{C}_{i0}$. In spite of these upfront or additional costs, we will show that in many cases, solving Eq. (4) can be significantly better than solving Eq. (2). Additionally, one can bypass the costs mentioned above of explicitly forming S by instead repeatedly applying the action of S to a vector during iterative solution. Such a method is not atypical of Krylov subspace iterative solvers [23].

3. Main numerical results

3.1. Methods

In this section, we demonstrate the advantage of solving Eq. (4) over Eq. (2) in simulating the system described in Sec. 2. To solve these linear systems, we apply an iterative method, which is preferable over direct methods for large systems [23]. Here we use the method of quasi-minimal residual (QMR), which has been demonstrated to work with the FDFD method, particularly for systems with inhomogenous dielectric distributions [21]. We also have tried a few other iterative solvers such as the generalized minimal residual (GMRES) and biconjugate gradient (BICG) method, neither of which converge as fast as QMR for our systems (but produce qualitatively similar results).

3.2. Results for two-dimensional systems

First, we verify that the solution recovered from the reduced system using Eq. (3) is accurate. To that end, since our test case is small enough, we use a direct solver to robustly and accurately produce a solution from the unreduced system (Eq. (2)) for comparison. In Fig. 2 we consider the structure in Fig. 1 surrounded by a periodic boundary condition. Figures 2(a) and 2(b) correspond to the dielectric case, where both the matrix A in Eq. (2) and the matrix S in Eq. (4) are real symmetric. Figures 2(c) and 2(d) correspond to the metallic case where both A and S are complex symmetric. In both cases, we see that the solutions from the direct solver (Figs. 2(a) and 2(c)) agree well with the solutions obtained from the domain decomposition approach with the solution to x₀ acquired with an iterative solver (Figs. 2(b) and 2(d)). The relative norm error between these two solutions in both cases is less than 10⁻¹⁰.

 

Fig. 2 The first column shows the fields (real for dielectric, absolute value for metallic) from the iterative solution (using QMR) of Eq. (4) for a) the dielectric case (ϵ = 12) and c) the metallic case (ϵ = −3 − 0.3i). The second column shows the equivalents obtained with a direct solver for b) the dielectric case and d) the metallic case. Since the fields are complex in the metallic case, we show the absolute value of the fields. The source wavelength in all cases is 1.5a

Download Full Size | PDF

To compare the convergence rates of iteratively solving Eq. (2) versus Eq. (4), we plot the relative residual error ρ = ║Ax − b║/║b║ for a solution x as a function of the number of iterations. For the unreduced system x, is obtained directly by solving Ax = b. For the reduced system, we solve for x₀ in Eq. (4) and then obtain x using Eq. (3).

Figure 3 shows the relative residual error as a function of the number of iterations. For both the dielectric case (Fig. 3(a)) and the metallic case (Fig. 3(b)), we see that the number of iterations to reach a relative residual error |ρ| < 10⁻¹¹ is roughly an order of magnitude lower when solving the reduced system than when solving the unreduced system. We note however that there is a scale difference between the dielectric and metallic cases. The dielectric case in general appears to require significantly more iterations to solve than a comparable metallic system.

 

Fig. 3 Convergence rates for the unreduced and reduced systems on the domain containing an 8 × 8 lattice of meta-atoms from Fig. 1 with a periodic boundary condition for a) the dielectric case and b) the metallic case. The source wavelength in both cases is 1.5a.

Download Full Size | PDF

Figure 4 shows the convergence rates of the same two systems considered in Fig. 3, but with a stretched-coordinate perfectly-matched-layer (PML) boundary condition surrounding the domain. The perfectly matched layer is 15 grid points deep. When applying the domain decomposition method, the perfectly matched layer is split into 36 subdomains as shown in Fig. 1. For both the dielectric system (Fig. 4(a)) and the metallic system (Fig. 4(b)), the unreduced system exhibits some stagnation while the reduced system still converges. Even with the preconditioner from [24], which symmetrizes the PML-terminated FDFD linear system, the problem persists. For both cases considered in Figs. 3 and 4, the convergence improvements from solving the reduced system also translate to significant wall-clock time improvements, an observation further discussed in Sec. 5.

 

Fig. 4 Comparison of convergence rates of the unreduced and reduced systems on the domain containing an 8 × 8 lattice of meta-atoms from Fig. 1 with a perfectly matched layer 15 grid points deep in a) the dielectric case and b) the metallic case. The source wavelength in both cases is 1.5a

Download Full Size | PDF

In Fig. 5 we focus on the dielectric case, and plot the number of iterations to convergence as a function of domain size. For each domain consisting of an N × N lattice of meta-atoms, we randomly place one of the three inclusions discussed above (circle, triangle, square) in each meta-atom. We place a point source at the center of the domain with a free-space wavelength of 1.5a. In Figs. 5(a) and 5(b) we consider the cases with either a perfectly-matched-layer or periodic boundary condition respectively, and compare the performance of iteratively solving the reduced and the unreduced cases. In all cases, the number of iterations to convergence increases as the domain expands. For the systems with a periodic boundary condition, the factor of improvement, defined as the ratio between the number of iterations to convergence for the unreduced and reduced systems, remains approximately constant (Fig. 5(a)). For the systems with a perfectly matched layer boundary condition, the unreduced system fails to converge when the domain size exceeds 8a, whereas the reduced system converges for all domain sizes considered, as shown in Fig. 5(b). We also note that in the simulations with a PML, simply applying the domain decomposition exclusively to the PML region suffices to make the iterative solver converge. In general, we observe that applying the domain decomposition method significantly improves convergence rates for the finite-difference frequency-domain discretizations of Maxwell’s equations discussed here.

 

Fig. 5 Number of iterations to convergence for the unreduced and reduced systems, as a function of increasing domain size. Here, the horizontal axis shows the number of meta-atoms per side in an N × N lattice of meta-atoms. We consider the dielectric case with a) perfectly matched layer boundary condition and b) periodic boundary condition. The source wavelength is 1.5a.

Download Full Size | PDF

3.3. Extension to three-dimensional systems

The domain decomposition method discussed above can be applied in three dimensions as well. In three dimensions, we solve:

To illustrate the domain decomposition method, we consider the structure shown in Fig. 6. The structure is a 3 × 3 × 3 lattice of cubic meta-atoms containing dielectric or metallic cubes of side length a/2 where a is the lattice constant. We use the same dielectric constants for the cubes as in the dielectric and metallic cases for the two-dimensional simulations. A dipole source with a free-space wavelength of 2.5a resides at the interface between the central meta-atom and a diagonal neighbor meta-atom (so the source is not inside the dielectric/metallic inclusion inside the central meta-atom). We discretize Eq. (5) with a Yee’s grid with a grid spacing of 0.053a. A periodic boundary condition terminates the domain.

 

Fig. 6 A lattice of cubic meta-atoms containing dielectric or metallic cubes. The volume between the dielectric or metallic cubes is air.

Download Full Size | PDF

For this structure, we consider the convergence rate of solving the reduced and unreduced systems. In the unreduced case, we directly solve Eq. (5) using the method discussed in [26]. In the reduced case, we solve Eq. (4); the x₀ in Eq. (3) corresponds to the field components at the boundaries between neighboring meta-atoms. In the dielectric case (Fig. 7(a)), iterative solution of the unreduced system stagnates while that of the reduced system converges. In the metallic case (Fig. 7(b)), solving the reduced system converges significantly faster than solving the unreduced system.

 

Fig. 7 Comparison of convergence rates for the unreduced and reduced systems terminated by a periodic boundary condition. a) corresponds to the dielectric case and b) the metallic case. In both cases, the reduced system converges in fewer iterations than the unreduced system. The source wavelength is 2.5a.

Download Full Size | PDF

Unlike the two-dimensional systems considered in the previous section, in three dimensions, for example the metallic case above, the reduction in the number of iterations to convergence does not necessarily translate to a reduction in net computational cost since the matrix for the reduced system can be quite dense. A number of ways to circumvent this problem exist; Sec. 5 provides a more detailed discussion of this issue.

4. Heuristics for understanding convergence improvements

The previous section provided numerical observations that the reduced system typically converges faster than the unreduced system. This section presents intuition regarding such observations and discusses the specific conditions under which the reduced system shows improvement over the unreduced system.

4.1. Mathematical background

As a basis for our heuristic argument, there is an extensive mathematics literature regarding the convergence behavior of iterative solution of symmetric positive definite linear systems Ax = b. For such systems, some insight into convergence behavior can be gained from analyzing the condition number of the matrix A:

The condition number bounds the convergence rate of specific classes of iterative solvers applied to symmetric positive definite linear systems [29, 30]. Though in general one cannot use condition number to understand the convergence behavior of linear systems that are not symmetric positive definite, in practice it often proves a useful guide. In that spirit, here we compare the conditioning of the unreduced system A with the reduced system S.

The Schur complement S of any sub-block D in a positive definite matrix A is at least as well conditioned as the original matrix A [27,28]. This weak bound can be further strengthened under various conditions. For a matrix A in the form of Eq. (2), [27, 31] show that:

In Eq. (7) ${\gamma }_B^2$ is the Cauchy-Bunyakowski-Schwarz (CBS) constant, given by:

4.2. Application to FDFD

The results above apply only to symmetric positive definite matrices. For other matrices, particularly indefinite matrices, the applicability of the results in the previous section is not guaranteed. However, in this section, we prove that an extreme conditioning contrast between A₀₀ and A is possible in FDFD systems, which matches the theoretical requirement that cond(A₀₀) < cond(A) discussed above. Additionally, we know that appropriate FDFD formulations of Maxwell’s equations, such as those used above, can result in systems that are very close to being spectrally definite [26], which may make the arguments in the previous section more applicable in practice. We also note that, in numerical experiments like those above, not using the appropriate formulations of Maxwell’s equations as mentioned above often diminished the reduced system’s convergence improvement over the unreduced system’s, or even made the reduced system’s convergence poorer than the unreduced system’s.

In the domain decomposition approach above, the diagonal block A₀₀ in Eq. (8) corresponds to the degrees of freedom in the domain in Fig. 1 on the interfaces of the meta-atoms. Below, we will show that A₀₀ is far better conditioned than the FDFD matrix A. Therefore, Eq. (8) indicates that the reduced matrix S can be better conditioned than A, which we show is numerically true for our systems even though they are indefinite.

For the matrix A of the unreduced system, its singular values (σ_i(A)) can be estimated as:

A₀₀’s maximum singular value is also approximately determined by Eq. (10), since A₀₀ is the matrix A restricted to the interfacial degrees of freedom, and the grid spacing remains the same on this subdomain. As a result, σ_max(A₀₀) ≈ σ_max(A).

To compare the smallest singular values of A and A₀₀, we first note that the smallest eigenfrequency ω_min,A of the derivative operator is overall determined by the size of the full domain as ${\omega }_{min,A}\approx \sqrt{{\displaystyle {\sum }_{j=1}^d{\left(c_A\pi /L_j\right)}^2}}$, where L_js are side lengths of the full domain and j indexes the spatial dimensions (d) to sum over. Also, the driving wavelength is typically chosen comparable to the lattice constant a in order to maximize the interaction between the wave and the meta-atoms. We have a ≪ L_j because the simulation domain is composed of several meta-atoms in each direction. Therefore, we have

Finally, we can approximate ω²µ₀ since the source wavelength is approximately the size of a:

On the other hand, the interface is very thin, typically only one or a few grid spacings thick, and therefore its thickness is much less than the transverse dimensions of the interface. Therefore, we have ${\omega }_{\min ,A_{00}}=\sqrt{{\displaystyle {\sum }_j{\left(c\pi /L_{j,interface}\right)}^2}}\approx c_{A_{00}}\pi /t$, where t is the thickness of the interface. Since t is much less than the lattice constant a, we have

Combining the results so far, we conclude

Combining the arguments in the preceding paragraphs for both the maximum and minimum singular values, we therefore have cond(A₀₀) ≪ cond(A). Now that we have proved that for FDFD, we can achieve large conditioning contrast between A₀₀ and A, we now demonstrate numerically that cond(S) ≪ cond(A). In Fig. 8, we consider the same system as in Fig. 5 for the dielectric case. We plot the condition numbers of A₀₀, S, and A as a function of domain size, with the domain surrounded by either a periodic boundary condition (Fig. 8a), or a perfectly matched layer boundary condition (Fig. 8b). In both Figs. 8a and 8b, the condition number of A₀₀ is largely independent of domain size, and is four orders of magnitude or more below the condition number of A. The condition number of S is around two orders of magnitude lower than the condition number of A. Such a large contrast in conditioning between A₀₀ and A, combined with the experimental observation that cond(S) ≪ cond(A), suggests that we can apply some of the insights from Sec. 4.1 to our system. In particular, the size of the conditioning improvement between A and A₀₀ is important because, were A symmetric positive definite, this means that conditioning improvement of S is still possible even if ${\gamma }_B^2\approx 1$, as long as $1/\left|1-{\gamma }_B^2\right|$ is less than the magnitude of the ratio of the conditioning of A over A₀₀.

 

Fig. 8 Condition number of the unreduced matrix (A), the reduced matrix (S), and the interface sub-block of A (A₀₀), as a function of domain size, for the physical system described in Fig. 5 where the dielectric inclusions have ϵ = 12 with (a) periodic and (b) perfectly matched layer boundary conditions.

Download Full Size | PDF

This conditioning analysis provides an intuitive explanation for why the use of domain decomposition techniques can significantly improve convergence in FDFD, as observed in Sec. 4. This analysis also indicates that in applying this domain decomposition technique, if possible it is better to choose the interfaces between subdomains to be in the region with lower dielectric constant, since in this case the condition number of A₀₀ should be smaller due to a higher value of σ_min(A₀₀).

5. Memory and wall-clock time considerations

The time required for iterative solution of a linear system Ax = b is roughly T_AI_A, where T_A denotes the time to evaluate the matrix-vector product Ax, and I_A represents the number of iterations to convergence. The preceding analysis focused on I_A. This section discusses T_A in the situation where we explicitly form S, and highlights directions for future work, especially in three dimensions.

When solving a linear system iteratively, the number of nonzero elements in the matrix roughly determines both the memory necessary to store the system, as well as the time for each matrix-vector product (matvec). To simplify our discussion, we consider a domain consisting of a d-dimensional cubic lattice of cubic meta-atoms. Additionally, we consider the thinnest possible interface we can construct. For example, with d = 2, such an interface consists of a layer one grid-spacing thick between each meta-atom, as illustrated by the dotted lines in Fig. 1. We define P_t as the number of nodes per side length of each meta-atom and N the number of meta-atoms along each dimension (for a total of N^d meta-atoms in the grid for a spatial dimension d). C(d) denotes the number of E-field components where d is the spatial dimension of the system. When d = 3, C(3) = 3 (see Eq. (5)) and when d = 2, C(2) = 1 (see Eq. (1)). Finally, Σ(d) denotes the size of the stencil for each node, which takes on a value of Σ(3) = 13 in three dimensions for the standard Yee’s grid and Σ(2) = 5 for a typical compact stencil in two dimensions. Hence, the unreduced system A contains approximately Σ(d)C(d)(NP_t)^d nonzero entries. Then for A the time cost for a matrix-vector product scales approximately as:

To discuss the reduced system, we let P_e denote the number of nodes per side length of the interface surrounding each meta-atom. In general P_e ≠ P_t unless the grid spacing on the interface is the same as that on the interior. The total number of nonzero elements in the reduced matrix S is then approximately $N^d{\left(dC(d)P_e^{d-1}\right)}^2$. To get this expression, we observe that for each of the N^d meta-atoms, we get a dense sub-block of the form $B_{0i}{D}_{ii}^{-1}{C}_{i0}$. We count the number of degrees of freedom on the interface of each sub-block, which is $dC(d)P_e^{d-1}$. Since the block is fully dense, we must then square this term. We also note that for d = 3, C(3) = 3 is not quite accurate as the thinnest possible interface on the Yee grid that decouples two adjacent meta-atoms only requires two components. However, to maintain consistent notation, we use C(3) = 3 as a strong upper bound. Finally, the same cost for S is roughly:

Now taking the ratio of Eq. (15) over Eq. (16):

First, consider a uniform grid such that P_e = P_t. In two dimensions, evaluating Eq. (17) gives T_A/T_S ≈ 5/4, independent of P_t or P_e. A ratio less than one indicates how much more costly evaluating a matvec for S is versus A. Conversely, a ratio greater than one means matvecs with S are cheaper to evaluate than A. Therefore, given a factor of 10 improvement in the convergence rates as demonstrated in Sec. 3, we see that in the two-dimensional case, significant wall-clock speed-up is possible. In three dimensions, however, the ratio in Eq. (17) acquires a factor of 1/P_t. Thus, both the memory cost for storing the matrix and the time required for matrix-vector products makes solving the reduced system less attractive than solving the unreduced system unless the improvement in I_A is substantially better than the increase in T_A. In the three-dimensional case considered in Sec. 5, for example, where P_t = 19, the storage for the reduced system is roughly a factor of 40 greater than that for the unreduced system. However, one can simply evaluate the action of S as discussed in Sec. 2, rather than forming and applying S explicitly, to bypass this issue.

For three-dimensional systems, in order to achieve significant reduction in computational cost via the domain decomposition approach described above, one can consider subgridding [33, 34]. Specifically, one might discretize the interfaces between meta-atoms more coarsely than the interiors of meta-atoms (αP_e = P_t in Eq. (13), where α is a positive integer greater than 1). Doing so could reduce the memory cost and matrix-vector product evaluation time for the reduced system relative to the unreduced system. Additionally, subgridding is natural for the types of systems considered in this paper. The need to capture all complex geometry, as well as all relatively rapid field variation within large-dielectric-constant materials, forces a fine discretization scale. In physical systems like these, all such complex geometry and large-dielectric-constant material typically appears within the meta-atoms’ interiors rather than on their interfaces. Thus, in principle the interfaces between meta-atoms need not be discretized at the same resolution as the interiors of meta-atoms.

6. Conclusions and final remarks

In this paper, we consider iterative solution of finite-difference frequency-domain discretizations of Maxwell’s equations, and show that a Schur complement domain decomposition method can significantly improve convergence rates when considering commonly used nanophotonic structures consisting of many repeated meta-atoms. The improvement in convergence rates can be heuristically understood in terms of the improvement in condition number of the relevant matrices.

The success of the approach here exploits the fact that the size of the meta-atoms is comparable to or smaller than the free space wavelength. The ability of the domain decomposition method to reuse solutions on subdomains containing repeated meta-atoms should facilitate optimization of photonic structures based on repeated meta-atoms. Further improvement of the method can be achieved via subgridding, as discussed in Sec. 5.

Aside from subgridding, one could precondition the reduced system S during iterative solution [32], which is important as the reduced system’s iterations to convergence and condition number increases with increasing domain size as exhibited in Fig. 5 and Fig. 8 respectively. Additionally, work has been done to find sparse approximations of the Schur complement [35], which can also mitigate storage costs mentioned in Sec. 5. Moreover, there are many preconditioning schemes on the unreduced linear system Eq. (2) that exploit domain decomposition and also avoid directly solving or iterating with S. Particular examples include algebraic recursive multi-level schemes (ARMS), which are described in [36]. The development here points to the possibility of using these methods to accelerate electromagnetic simulations.