Robust fast direct integral equation solver for quasiperiodic scattering problems with a large number of layers
Abstract
We present a new boundary integral formulation for timeharmonic wave diffraction from twodimensional 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 quasiperiodic Green’s function method which fails whenever any of the layers approaches a Wood anomaly. We achieve this by a decomposition into near and farfield contributions. The former uses the freespace Green’s function in a secondkind integral equation on one period of the material interfaces and their immediate left and right neighbors; the latter uses proxy point sources and small leastsquares solves (Schur complements) to represent the remaining contribution from distant copies. By using highorder 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 1000layer structure with unknowns to 9digit accuracy in 2.5 minutes on a desktop workstation.
1 Introduction
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 highpowered laser [50, 11] or wideband [31] applications rely on structures with up to 50 dielectric layers. Efficient thinfilm photovoltaic cells exploit multiple layers of silicon, transparent conductors, and dielectrics, which are patterned to enhance absorption [4, 32]. For solar thermal power, efficient visiblelight absorbers which reflect in the infrared require patterning of several layers [52]. Related wave scattering problems appear in photonic crystals [30], process control in semiconductor lithography [44], 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 semiinfinite. The wavenumber will be in the th layer. A plane wave is incident in the uppermost layer,
(1) 
with wavevector , at angle . The incident wave is quasiperiodic (periodic up to a phase), meaning for all , where the Bloch phase (phase factor associated with translation by one unit cell) is
(2) 
Note that is controlled by the period and indicent wave alone. We will seek a solution sharing this quasiperiodic symmetry.
As is standard for scattering theory [17], 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,
(3) 
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.
(4) (5) 
Quasiperiodicity in all layers, i.e. for all ,
(6) 
Outgoing radiation conditions in and , namely the uniform convergence of Rayleigh–Bloch expansions [13] in the upper and lower halfspaces,
(7) (8) 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 farfield diffraction pattern.
The above BVP describes timeharmonic acoustics with layers of different sound speeds (but the same densities) [17], or the timeharmonic 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 divergencefree vector fields and satisfy
(9)  
(10) 
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 illposed. They are distinct from (but in the literature sometimes confused with) Wood anomalies [57], 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 threedimensional (3D) case by Shipman [53].
There are many loworder numerical methods used to solve multilayer scattering problems, which in test problems may only agree to 1 digit of accuracy [56]. Finite difference timedomain (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 frequencydomain, 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 [5]. They also demand artificial absorbing boundary conditions (perfectly matched layers). The rigorouscoupled wave analysis (RCWA) or Fourier Modal Method is specially designed for multilayer gratings [45]; 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 loworder “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 loworder convergence. In general, when the layers are strictly planar, the problem becomes 1D [16] and the (unstable) transfer matrix and (stable) scattering matrix approaches are standard [33]. 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 highorder 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 quasiperiodic scattering, the interfaces are unbounded and the usual approach is to replace the 2D freespace kernel (Green’s function) for waves in layer ,
(11) 
(where is the Hankel function of order zero [1]), by the quasiperiodic one obeying (6), in which case the problem may be formulated on a single period of the interface. The usual quasiperiodic Green’s function for layer is
(12) 
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 [40]. It is convenient to expand slightly the definition of Wood anomaly, as follows.
Definition 1
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 nonrobustness 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 contourintegral approach of Barnett–Greengard [10] does not generalize well to multiple layers.
In this paper we introduce a simpler class II BIE method which combines the freespace 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 farfield 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 [24]. Related representations have been used for some time in the engineering community [25]. A modern interpretation of the key idea is that the farfield contribution is smooth—the interaction between distant periodic copies and the central unit cell is low rank; eg see [21]—and hence only a small number of proxy points is needed, at least if the frequency is not too high.
Remark 1
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 [10] 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 finitethickness layers, and impose outgoing conditions only in the semiinfinite layers 1 and . Since nondivergent Green’s functions do exist which satisfy these minimal radiation conditions, they are selected by the leastsquares 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 blocktridiagonal 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 secondkind 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
Using (11) we define standard potentials for the Helmholtz equation, the single and doublelayer representations [17] lying on a general curve , at wavenumber for the th layer,
(13) 
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,
(14) 
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,
(15) 
where is the outwardsdirected 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 nearfield single and doublelayer potentials and proxy representations in each layer we have, recalling the notation for in ,
(16)  
(17)  
(18) 
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, quasiperiodicity, 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 this subsection, matching conditions (4) and (5) will be enforced at all material interfaces in a standard Müller–Rokhlin [46, 51] scheme.
In the indirect approach, boundary integral operators arise from the restriction of representations (13) to curves [17]. Following [10] we use notation to indicate the singlelayer operator at wavenumber from a source curve to target curve . Similarly we use for the doublelayer operator, for the targetnormal derivative of the singlelayer operator, and for the targetnormal derivative of the doublelayer operator. As before, we use a tilde to indicate summation over the source curve and its phased nearest neighbors, thus, for a target point ,
(19)  
(20) 
When the target curve is the same as the source (), we note that the singlelayer operator is a weakly singular integral operator, that the action of the doublelayer 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
(21)  
(22)  
(23)  
(24) 
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,
(25)  
(26) 
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
(27)  
(28) 
On the bottom interface , the only change is the absence of coupling from any lower interface, so,
(29)  
(30) 
We wish to write these in a more compact form, hence we pair up double and singlelayer densities, then stack them into a single column vector,
(31) 
Similarly we stack the proxy strength vectors , and form a vector of righthand side functions,
(32) 
Then all of the coupled BIEs and functional equations (25)(30) can be compactly grouped into the matrixtype notation,
(33) 
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,
(34) 
where is the identity operator. is an by matrix, each of whose entries is a stack of continous function columns (sometimes called a quasimatrix) expressing the effect of each proxy point strength on the value and normal derivative functions on . The only nonzero blocks of are
(35) 
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 offdiagonal 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 quasiperiodicity conditions
Quasiperiodicity (6) will be enforced in each layer by matching both values and normal derivatives between the left and right walls. Since the PDE is secondorder, matching two functions (values and normal derivatives) is sufficient Cauchy data to guarantee extension as a quasiperiodic 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 nearfield sum) down to two,
(36) 
For quasiperiodicity 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:
(37) 
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,
(38)  
(39) 
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,
(40) 
2.4 Imposing the radiation conditions
First we enforce the upward radiation condition (7) at the artificial interface (with upwardpointing normal), substituting the layer1 representation (16) to get,
(41) 
Matching values at is not enough: we also need to match normal () derivatives, to ensure that the secondorder PDE solution continues smoothly through , thus,
(42) 
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,