Robust fast direct integral equation solver for quasi-periodic scattering problems with a large number of layers
We present a new boundary integral formulation for time-harmonic wave diffraction from two-dimensional structures with many layers of arbitrary periodic shape, such as multilayer dielectric gratings in TM polarization. Our scheme is robust at all scattering parameters, unlike the conventional quasi-periodic Green’s function method which fails whenever any of the layers approaches a Wood anomaly. We achieve this by a decomposition into near- and far-field contributions. The former uses the free-space Green’s function in a second-kind integral equation on one period of the material interfaces and their immediate left and right neighbors; the latter uses proxy point sources and small least-squares solves (Schur complements) to represent the remaining contribution from distant copies. By using high-order discretization on interfaces (including those with corners), the number of unknowns per layer is kept small. We achieve overall linear complexity in the number of layers, by direct solution of the resulting block tridiagonal system. For device characterization we present an efficient method to sweep over multiple incident angles, and show a speedup over solving each angle independently. We solve the scattering from a 1000-layer structure with unknowns to 9-digit accuracy in 2.5 minutes on a desktop workstation.
Periodic geometries (such as diffraction gratings and antennae) and multilayered media (such as dielectric mirrors) are both essential for the manipulation of waves in modern optical and electromagnetic devices. In an increasing number of applications these features occur in tandem: for instance, performance of a grating with transverse periodicity will be enhanced by using many layers of different indices. Dielectric gratings for high-powered laser [50, 11] or wideband  applications rely on structures with up to 50 dielectric layers. Efficient thin-film photovoltaic cells exploit multiple layers of silicon, transparent conductors, and dielectrics, which are patterned to enhance absorption [4, 32]. For solar thermal power, efficient visible-light absorbers which reflect in the infrared require patterning of several layers . Related wave scattering problems appear in photonic crystals , process control in semiconductor lithography , in the electromagnetic characterization of increasingly multilayered integrated circuits, and in models for underwater acoustic wave propagation. In most of these settings, it is common to solve for the scattering for a large number of incident angles and/or wavelengths, then repeat this inside a design optimization loop. Thus, a robust and efficient solver is crucial. We present such a solver, which scales optimally with respect to the number of layers.
Let us describe the geometry of the problem (Fig. 1(a)). Consider interfaces , each of which has the same periodicity in the horizontal () direction. The interfaces lie between homogeneous material layers, each filling a domain , . The th layer lies between and , whilst the top and bottom layers are semi-infinite. The wavenumber will be in the th layer. A plane wave is incident in the uppermost layer,
with wavevector , at angle . The incident wave is quasi-periodic (periodic up to a phase), meaning for all , where the Bloch phase (phase factor associated with translation by one unit cell) is
Note that is controlled by the period and indicent wave alone. We will seek a solution sharing this quasi-periodic symmetry.
As is standard for scattering theory , the incident wave causes a scattered wave to be generated, and the physical wave is their total . The scattered wave is given by solving the following boundary value problem (BVP). We have the Helmholtz equation in each layer,
where we write for the scattered wave in the th layer, and the following interface, boundary, and radiation conditions:
Continuity of the value and derivative the total wave on each interface, i.e.
Quasi-periodicity in all layers, i.e. for all ,
Outgoing radiation conditions in and , namely the uniform convergence of Rayleigh–Bloch expansions  in the upper and lower half-spaces,
where the horizontal wavenumbers in the modal expansion are
and the upper and lower vertical wavenumbers are and respectively, with all signs taken as positive real or imaginary. The complex coefficients and are the Bragg diffraction amplitudes of the reflected and transmitted th order; only the orders for which or carry power away from the system, and they characterize the far-field diffraction pattern.
The above BVP describes time-harmonic acoustics with layers of different sound speeds (but the same densities) , or the time-harmonic full Maxwell equations in the case of a -invariant multilayer dielectric in TM polarization, as we now briefly review; see [54, 29]. (The modification to the interface conditions for differing densities or TE polarization are simple, and we will not present them.) Letting the time dependence be , Maxwell’s equations state that the divergence-free vector fields and satisfy
Writing and , and eliminating gives in the th layer the Helmholtz PDE (3) with wavenumber , where and are the permittivity and permeability of the layer. Continuity of transverse and then gives the interface conditions (4)-(5). Note that all of the layer wavenumbers are scaled linearly by the overall frequency .
The BVP (3)–(8) has a solution for all parameters [13, Thm. 9.2]. The solution is unique at all but a discrete set of frequencies when is fixed [13, Thm. 9.4]; these frequencies correspond to guided modes of the dielectric structure, where resonance makes the physical problem ill-posed. They are distinct from (but in the literature sometimes confused with) Wood anomalies , which are scattering parameters (,) for which or for some , making the upper or lower th Rayleigh–Bloch mode a horizontally traveling plane wave. A Wood anomaly does not prevent the solution from being unique, although it does become arbitrarily sensitive to changes in or . For more detail see [41, 10], and the extensive review in the three-dimensional (3D) case by Shipman .
There are many low-order numerical methods used to solve multilayer scattering problems, which in test problems may only agree to 1 digit of accuracy . Finite difference time-domain (FDTD) [55, 27] is easy to code but has dispersion errors, and requires artificial absorbing boundary conditions and arbitrarily long settling times near resonances. Direct discretization in the frequency-domain, as in finite difference (FD) or finite element methods, is also possible [6, 20, 44], although they require a large number of unknowns, and, as the frequency grows, “pollution” error means that an increasing number of degrees of freedom per wavelength is needed . They also demand artificial absorbing boundary conditions (perfectly matched layers). The rigorous-coupled wave analysis (RCWA) or Fourier Modal Method is specially designed for multilayer gratings ; to overcome slow convergence a Fourier factorization method is needed [39, 38]. However, RCWA is not easy to apply for arbitrary shapes (relying on an intrinsically low-order “staircase” approximation of layer shapes), nor to generalize to 3D. Other methods include volume integral equations [31, 37], for which it is hard to exceed low-order convergence. In general, when the layers are strictly planar, the problem becomes 1D  and the (unstable) transfer matrix and (stable) scattering matrix approaches are standard . However, we are concerned with interfaces of arbitrary shape.
Since the medium is piecewise constant, boundary integral equations (BIE) formulated on the interfaces are natural and mathematically rigorous [17, 47, 43]. By exploiting the reduced dimensionality, the number of unknowns is much reduced, and high-order quadratures exist, in 2D [34, 26] but also in 3D. Combined with fast algorithms for handling the resulting linear systems, this leads to much higher efficiency and accuracy than FD methods, and has started to be used effectively in periodic problems [43, 3, 48, 49, 14, 21, 23]. In the setting of quasi-periodic scattering, the interfaces are unbounded and the usual approach is to replace the 2D free-space kernel (Green’s function) for waves in layer ,
(where is the Hankel function of order zero ), by the quasi-periodic one obeying (6), in which case the problem may be formulated on a single period of the interface. The usual quasi-periodic Green’s function for layer is
a sum whose slow convergence renders it computationally useless. Thus a large industry has been build around efficient evaluation of using convergence acceleration, Ewald’s method, or lattice sums . It is convenient to expand slightly the definition of Wood anomaly, as follows.
We say that layer is at a Wood anomaly if either or (or both) for some .
The problem then with -based methods is that (12) does not exist (the sum diverges) whenever the th layer is at a Wood anomaly. As the number of different materials increases in a structure, the chances of some layer hitting (or being close to) a Wood anomaly, and thus of failure, increases.
Two classes of solutions to this non-robustness problem have recently been introduced:
In the multilayer dielectric setting robustness using class I would require the impedance Green’s function, which is difficult to evaluate; while the class II contour-integral approach of Barnett–Greengard  does not generalize well to multiple layers.
In this paper we introduce a simpler class II BIE method which combines the free-space Green’s function for the unit cell and neighbors, with a ring of proxy point sources (i.e. the method of fundamental solutions, or MFS [12, 7]) to represent the far-field contributions which “periodize” the field. This combines ideas in [9, Sec. 3.2] and the fast direct solver community [42, Sec. 5], and has been independently proposed recently for Laplace problems by Gumerov–Duraiswami . Related representations have been used for some time in the engineering community . A modern interpretation of the key idea is that the far-field contribution is smooth—the interaction between distant periodic copies and the central unit cell is low rank; eg see —and hence only a small number of proxy points is needed, at least if the frequency is not too high.
Conveniently, in our new formulation we can take Schur complements to eliminate the proxy strength unknowns for each layer without recreating (12) and its associated Wood anomaly problem, as would happen in prior methods [9, 10] (see [9, Remark 8]). The difference is that in  both upward and downward radiation conditions (of the type (7) and (8)) are imposed on the Green’s function, making it equivalent to (12), whereas we do not impose any radiation condition in the finite-thickness layers, and impose outgoing conditions only in the semi-infinite layers 1 and . Since non-divergent Green’s functions do exist which satisfy these minimal radiation conditions, they are selected by the least-squares linear algebra in the Schur complement (see Sec. 3.2).
We present our new representation and its discretization in Sec. 2, then combine it in Sec. 3 with a direct solver which has two steps: Schur complements to eliminate the proxy unknowns, followed by direct block-tridiagonal factorization. The tridiagonal structure arises simply because layer couples only with layers and . The overall scaling is , i.e. linear in the number of layers and cubic in the number of unknowns per layer. This allows our solver to tackle problems with of order unknowns in only a few minutes. Since the solution is direct, as explained in Sec. 3.4 we can solve new incident waves that share the same without extra effort, and, by reusing matrix blocks, handle other values efficiently. We test the solver’s error and speed performance with a variety of interface shapes, with up to layers, with random or periodic permittivities, and multiple incident angles including a Wood anomaly, in Sec. 4. We conclude with a summary and implications for future work in Sec. 5.
2 Boundary integral formulation, periodizing scheme, and its discretization
We now reformulate the BVP as a system of linear second-kind integral equations on the interfaces , which lie in a single unit cell, coupled with linear conditions on fictitious unit cell walls; the complete system will be summarized by (49) below. A little extra geometry notation is needed, as shown in Fig. 1. Let us define the (central) unit cell as the vertical strip of width lying between and ; of course its horizontal displacement is arbitrary. The blue dashed vertical lines and are the left and right boundaries of the layer domains lying inside the unit cell. The proxy points for layer lie on the circle (shown by red dotted lines). The magenta dashed lines and are fictitious interfaces for radiation conditions located at and , touching and , respectively.
2.1 Representation of the scattered wave
where is the unit normal on the curve at , and the arclength element. Shortly we will set to be either or , with the normals pointing down (into the layer below the interface). Integral representations which include phased contributions from the nearest neighbors are indicated with a tilde,
Let the proxy points lie uniformly on the circle of radius which is centered on the domain . As is well known in MFS theory, increasing allows a higher convergence rate with respect to [7, Thm. 3]; however, since the proxy points are representing the contributions from periodic interface copies and , which thus have singularities at , the proxy charge strengths will turn out to be exponentially large [7, Thm. 7] if exceeds by much. In practice we choose . In order to make the proxy representation more robust at high wavenumbers we use a “combined field” approach, choosing the proxy basis functions for th layer,
where is the outwards-directed unit normal to the circle at the th proxy point. This results in smaller coefficients than if monopoles or dipoles alone were used (which can be justified by treating the proxy points as a discrete approximation to a layer potential on , and considering arguments in [58, Sec. 7.1]).
Combining the near-field single- and double-layer potentials and proxy representations in each layer we have, recalling the notation for in ,
By construction, for all layers , and for all density functions and and proxy unknown vectors , this representation satisfies the Helmholtz equations (3). In the following subsections, we describe in turn how to enforce the interface matching, quasi-periodicity, and radiation conditions. Each of these three conditions will comprise a block row of the final linear system (49) that enables us to solve for the densities and proxy unknowns.
2.2 Matching conditions at material interfaces
In the indirect approach, boundary integral operators arise from the restriction of representations (13) to curves . Following  we use notation to indicate the single-layer operator at wavenumber from a source curve to target curve . Similarly we use for the double-layer operator, for the target-normal derivative of the single-layer operator, and for the target-normal derivative of the double-layer operator. As before, we use a tilde to indicate summation over the source curve and its phased nearest neighbors, thus, for a target point ,
When the target curve is the same as the source (), we note that the single-layer operator is a weakly singular integral operator, that the action of the double-layer and its transpose must be taken in their principal value sense, and that the operator is hypersingular.
At the first interface , and are coupled. The functions , , , and at can be found by letting in (16)–(17) approach from the respective side, and using the standard jump relations [17, Thm. 3.1 and p.66] which introduce terms of one half times the density to each and term, giving
On this interface only, the matching conditions (4) include the indicent wave, and enforcing them using the above gives the inhomogeneous coupled BIE and functional equations,
Note that the half density terms added to give a in the first equation and a in the second; these terms appear for every layer and will cause the BIE to be of Fredholm second kind.
On the middle interfaces , , we similarly match and and their normal derivatives, noting that now there is coupling to both the above and below interfaces, but no effect of the incident wave, to get
On the bottom interface , the only change is the absence of coupling from any lower interface, so,
We wish to write these in a more compact form, hence we pair up double- and single-layer densities, then stack them into a single column vector,
Similarly we stack the proxy strength vectors , and form a vector of right-hand side functions,
where is a -by- matrix, each of whose entries is a block of operators which maps to a pair of functions (i.e. values then normal derivatives) on . Every block of is zero apart from the following tridiagonal entries,
where is the identity operator. is an -by- matrix, each of whose entries is a stack of continous function columns (sometimes called a quasi-matrix) expressing the effect of each proxy point strength on the value and normal derivative functions on . The only nonzero blocks of are
The term in (33) is precisely (barring the summation over neighbors) the Müller–Rokhlin formulation [46, 51] for multiple material interfaces. This is of Fredholm second kind since the off-diagonal blocks in (34) have continuous kernels, and cancellation of the leading singularities occurs in the pairs in parentheses in (34), leaving the diagonal operators at most weakly singular, hence compact.
2.3 Imposing the quasi-periodicity conditions
Quasi-periodicity (6) will be enforced in each layer by matching both values and normal derivatives between the left and right walls. Since the PDE is second-order, matching two functions (values and normal derivatives) is sufficient Cauchy data to guarantee extension as a quasi-periodic solution.
We evaluate the first layer representation (16) on the walls, and exploit the following simplification due to translational symmetry (as in [10, 9]) which cancels six terms (three from each near-field sum) down to two,
For quasi-periodicity we wish this function to vanish, so we make it the first operator block row of a homogeneous linear system. Doing the same for the normal derivatives on the and walls, and then for similar conditions for all other layers , gives equations that can be written compactly with a matrix notation as follows:
where is an -by- matrix, each entry of which is a block of operators mapping interface densities to wall values and normal derivatives. Every block of is zero apart from the bidiagonal blocks,
for and respectively. is an -by- matrix, each entry of which is a stack of function columns (as with ), but only the diagonal entries are nonzero,
2.4 Imposing the radiation conditions
Matching values at is not enough: we also need to match normal () derivatives, to ensure that the second-order PDE solution continues smoothly through , thus,
We will truncate the Rayleigh–Bloch expansion to terms, from to , since it is exponentially convergent once exceeds (in the upper layer) and (lower layer). We also apply the downward radiation condition (8) at to the representation (18), giving a second set of homogeneous linear conditions. We choose the normals of and both to point in the upward sense. As with and , we stack all coefficients into a single vector,