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 $25\times$ speedup over solving each angle independently. We solve the scattering from a 1000-layer structure with $3\times 10^5$ unknowns to 9-digit accuracy in 2.5 minutes on a desktop workstation.


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 high-powered laser [50,11] or wideband [31] 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 [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  Figure 1: (a) Geometry of scattering problem. Γ i for i = 1, . . . , I are the material interfaces with one spatial period lying between the vertical blue dotted lines. The medium is uniform in the ith layer, which lies above the ith interface and has wavenumber k i . Our algorithm also uses: L i and R i which are the left and right walls of one period Ω i of the ith layer, P i the proxy circle for this layer, and U , D the upper and lower fictitious interfaces (at y = y U and y = y D ) where the radiation condition is applied. (b) Zoom of the top part of the geometry, showing quadrature nodes for Nyström method and collocation (for clarity, less nodes are shown than actually used).
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 I interfaces Γ i , each of which has the same periodicity d in the horizontal (x) direction. The interfaces lie between I + 1 homogeneous material layers, each filling a domain Ω i ⊂ R 2 , i = 1 . . . , I + 1. The ith layer lies between Γ i and Γ i−1 , whilst the top and bottom layers are semi-infinite. The wavenumber will be k i in the ith layer. A plane wave is incident in the uppermost layer, u inc (r) = e ik·r , r ∈ Ω 1 0, otherwise (1) with wavevector k = (k 1 cos θ inc , k 1 sin θ inc ), at angle −π < θ inc < 0. The incident wave is quasi-periodic (periodic up to a phase), meaning u inc (x + d, y) = αu inc (x, y) for all (x, y) ∈ R 2 , 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 [17], the incident wave causes a scattered wave u to be generated, and the physical wave is their total u inc + u. The scattered wave is given by solving the following boundary value problem (BVP). We have the Helmholtz equation in each layer, where we write u i for the scattered wave in the ith layer, and the following interface, boundary, and radiation conditions: • Continuity of the value and derivative the total wave on each interface, i.e.
• Outgoing radiation conditions in u 1 and u I+1 , namely the uniform convergence of Rayleigh-Bloch expansions [13] in the upper and lower half-spaces, u I+1 (x, y) = n∈Z a D n e iκnx e ik D n (−y+y D ) for y ≤ y D , where the horizontal wavenumbers in the modal expansion are κ n := k 1 cos θ inc + 2πn d , and the upper and lower vertical wavenumbers are k U n = k 2 1 − κ 2 n and k D n = k 2 I+1 − κ 2 n respectively, with all signs taken as positive real or imaginary. The complex coefficients a U n and a D n are the Bragg diffraction amplitudes of the reflected and transmitted nth order; only the orders for which k U n > 0 or k D n > 0 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) [17], or the time-harmonic full Maxwell equations in the case of a z-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 e −iωt , Maxwell's equations state that the divergence-free vector fields E and H satisfy Writing E = (0, 0, u) and H = (H x , H y , 0), and eliminating H gives in the ith layer the Helmholtz PDE (3) with wavenumber k i = ω √ ε i µ i , where ε i and µ i are the permittivity and permeability of the layer.
Continuity of transverse E and H then gives the interface conditions (4)- (5). Note that all of the layer wavenumbers k i 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 θ inc 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 [57], which are scattering parameters (θ inc ,ω) for which k n U = 0 or k n D = 0 for some n, making the upper or lower nth 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 θ inc or ω. For more detail see [41,10], and the extensive review in the three-dimensional (3D) case by Shipman [53].
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 [56]. 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 [5]. 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 [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 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 [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 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 i, is the Hankel function of order zero [1]), 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 i is a sum whose slow convergence renders it computationally useless. Thus a large industry has been build around efficient evaluation of G QP 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 i is at a Wood anomaly if either κ n = k i or κ n = −k i (or both) for some n ∈ Z.
The problem then with G QP -based methods is that (12) does not exist (the sum diverges) whenever the ith 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: I. Replace (12) by the quasi-periodic Green's function for the Dirichlet [43] or impedance half-plane problems [28], or their generalization to larger numbers of images [15].
II. Return to the free-space Green's function (11) for the central unit cell plus immediate neighbors, plus a new representation of far-field contributions, imposing the quasi-periodicity condition in the least-squares sense via additional rows in the linear system [9,10,21].
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 [10] 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 [24]. Related representations have been used for some time in the engineering community [25]. 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 [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 finite-thickness layers, and impose outgoing conditions only in the semi-infinite layers 1 and I + 1. 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 blocktridiagonal factorization. The tridiagonal structure arises simply because layer i couples only with layers i − 1 and i + 1. The overall scaling is O (IN 3 ), i.e. linear in the number of layers and cubic in N the number of unknowns per layer. This allows our solver to tackle problems with N I of order 10 6 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 I = 1000 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.

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 Γ i , i = 1, . . . , I 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 d lying between x = −d/2 and x = d/2; of course its horizontal displacement is arbitrary. The blue dashed vertical lines are the left and right boundaries of the layer domains {Ω i } I+1 i=1 lying inside the unit cell. The proxy points for layer i lie on the circle P i (shown by red dotted lines). The magenta dashed lines U and D are fictitious interfaces for radiation conditions located at y = y U and y = y D , touching Ω 1 and Ω I+1 , respectively.

Representation of the scattered wave
Using (11) we define standard potentials for the Helmholtz equation, the single-and double-layer representations [17] lying on a general curve W , at wavenumber k i for the ith layer, where n is the unit normal on the curve W at r , and ds the arclength element. Shortly we will set W to be either Γ i−1 or Γ i , 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 {y i p } P p=1 ∈ R 2 lie uniformly on the circle P i of radius R which is centered on the domain Ω i . As is well known in MFS theory, increasing R allows a higher convergence rate with respect to P [7, Thm. 3]; however, since the proxy points are representing the contributions from periodic interface copies {−∞, . . . , −3, −2} and {2, 3, . . . , ∞}, which thus have singularities at |x| > 3d/2, the proxy charge strengths will turn out to be exponentially large [7,Thm. 7] if R exceeds 3d/2 by much. In practice we choose R ∈ [3d/2, 2d]. In order to make the proxy representation more robust at high wavenumbers we use a "combined field" approach, choosing the proxy basis functions for ith layer, where n p is the outwards-directed unit normal to the circle P i at the pth 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 P i , 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 u i for u in Ω i , By construction, for all layers i = 1, . . . , I + 1, and for all density functions σ i and τ i and proxy unknown vectors c i := {c i p } P p=1 , 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.

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 S i V,W to indicate the single-layer operator at wavenumber k i from a source curve W to target curve V . Similarly we use D i V,W for the double-layer operator, D i, * V,W for the target-normal derivative of the single-layer operator, and T i V,W 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 x ∈ V , When the target curve is the same as the source (V = W ), 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 T operator is hypersingular.
At the first interface Γ 1 , u 1 and u 2 are coupled. The functions u 1 , u 2 , ∂u1 ∂n , and ∂u2 ∂n at Γ 1 can be found by letting r in (16)-(17) approach Γ 1 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 D and D * 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 −τ 1 in the first equation and a +σ 1 in the second; these terms appear for every layer and will cause the BIE to be of Fredholm second kind. On the middle interfaces Γ i , i = 2, . . . , I −1, we similarly match u i and u i+1 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 Γ I , 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 c i = {c i p } P p=1 , and form a vector of right-hand side functions, Then all of the coupled BIEs and functional equations (25)-(30) can be compactly grouped into the matrixtype notation, where A is a I-by-I matrix, each of whose entries A i,j is a 2 × 2 block of operators which maps η j to a pair of functions (i.e. values then normal derivatives) on Γ i . Every block of A is zero apart from the following tridiagonal entries, , i = 1, 2, · · · , I, where I is the identity operator. B is an I-by-(I + 1) matrix, each of whose entries B i,j is a stack of P continous function columns (sometimes called a quasi-matrix) expressing the effect of each proxy point strength c j p on the value and normal derivative functions on Γ i . The only nonzero blocks of B are The term Aη 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.

Imposing the quasi-periodicity conditions
Quasi-periodicity (6) will be enforced in each layer by matching both values and normal derivatives between the left L i and right R i = L i + d 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 L 1 and R 1 walls, and then for similar conditions for all other layers i = 2, . . . , I + 1, gives equations that can be written compactly with a matrix notation as follows: where C is an (I + 1)-by-(I + 1) matrix, each entry of which is a 2 × 2 block of operators mapping interface densities to wall values and normal derivatives. Every block of C is zero apart from the bidiagonal blocks, for i = 1, 2, · · · , I + 1 and i = 2, 3, · · · , I + 1, respectively. Q is an (I + 1)-by-(I + 1) matrix, each entry of which is a stack of P function columns (as with B i,j ), but only the diagonal entries are nonzero,

Imposing the radiation conditions
First we enforce the upward radiation condition (7) at the artificial interface U (with upward-pointing normal), substituting the layer-1 representation (16) to get, Matching values at U is not enough: we also need to match normal (y) derivatives, to ensure that the second-order PDE solution continues smoothly through U , thus, We will truncate the Rayleigh-Bloch expansion to 2K +1 terms, from n = −K to K, since it is exponentially convergent once |κ n | exceeds k 1 (in the upper layer) and k I+1 (lower layer). We also apply the downward radiation condition (8) at D to the representation (18), giving a second set of homogeneous linear conditions. We choose the normals of U and D both to point in the upward sense. As with η and c, we stack all coefficients into a single vector, The resulting conditions can again be written in a simple matrix form: where in which To clarify, in W U and W D , the 2K +1 columns are pairs of Fourier functions evaluated over x ∈ (−d/2, d/2), the x-coordinate extent of the lines U and D.

Discretization of functions and operators
Finally, combining the linear conditions from the previous three subsections, we have the coupled BIE and functional equations, Recall that η contains unknown density functions, while c and a are discrete coefficient vectors. On the right hand side, f involves functions (from the incident wave), and each 0 is a stack of zero functions. Thus A, C, and Z are blocks of operators, while the other six matrix blocks involve quasi-matrices (stacks of function columns). For numerical computation the continuous variables must be discretized, turning each function into a discrete set of values, and each operator block into a matrix. This is simple for the functional conditions in the 2nd and 3rd block rows of (49): we just sample at a discrete set of collocation points. For the 2nd block row, we use M w nodes {x i m } Mw m=1 of a Gauss-Legendre quadrature living on the left wall L i for the ith layer. See Fig. 1(b). Thus each diagonal block of Q (40) is replaced by a 2M w -by-2P matrix Q i with elements Similarly for the 3rd row, we use M equally-spaced (trapezoid rule) nodes {x U m } M m=1 on U , and {x D m } M m=1 on D. The trapezoid rule is appropriate here since functions will be periodic. Inserting these nodes, the formulae for the matrices V U , V D , W U , and W D discretizing (47)- (48) are similar to (50) above.
To discretize the remaining blocks A, B, C and Z, we need to fix a set of quadrature nodes {z i j } Ni j=1 on each interface Γ i . These nodes have associated weights {w i j } Ni j=1 , such that for any smooth d-periodic function f on Γ i , the quadrature rule holds to high accuracy. To choose these nodes and weights, we first consider the case of Γ i a smooth interface (e.g. Γ 1 in Fig. 1(b)). Let one period of the interface be parametrized by a vector function Z(s) for 0 ≤ s < 2π. By changing variable, the periodic trapezoid rule s j = 2π(j − 1/2)/N i , j = 1, . . . , N i in parameter s gives a quadrature rule z i j = Z(s j ) and w i j = (2π/N i )|Z (s j )|. Then for (C ∞ ) smooth d-periodic integrands on Γ i the error in (51) will be superalgebraically convergent [19, (2.9.16)]. On the other hand, if Γ i has corners, it breaks into one or more "segments" (e.g. Γ 2 in Fig. 1(b)). These need not be straight lines, merely smooth. Say a segment is again parametrized by a function Z(s) for 0 ≤ s < 2π. Then we reparametrize it via Z(w(s)) using the corner grading function suggested by Kress [34, (6.9)], where q controls the grading at endpoints. Higher q will cause more nodes to be close to the endpoints; typically we choose q = 6 or higher. Let n i,l be the number of nodes used for the lth segment of Γ i . Then a separate trapezoid rule s j = 2πj/n i,l , j = 1, . . . , n i,l is used on each segment, making N i = l n i,l nodes in total. The formula for the nodes and weights 1 are the same as in the smooth case, with the proviso that Z is replaced by the composed function Z • w. The grading function, since its derivative vanishes to high order at the endpoints, insures that the trapezoid rule achieves high order accuracy (typically order q). We find this more efficient for up to 10-digit accuracy than dyadically-refined panel quadratures [23]. With interface quadratures defined, blocks B, C, and Z are easy to discretize by restriction of the continuous variable to the set of nodes. For example, each block B i,i in (35), describing the interaction of the ith proxy basis with Γ i , is replaced by a 2N i -by-P matrix B i,i with elements The matrix for B i,i+1 is similar. Operator blocks C and Z involve boundary integral representations over interfaces: each integral is replaced by a sum according to (51). For example, the M w -by-N i matrix discretizing the upper-right block of C i,i in (38) . . , M w and j = 1, . . . , N i . Other blocks of C and Z are discretized similarly; for the reader's sanity we refrain from giving all formulae.
Finally we discretize A. The operators in the blocks A i,j for i = j, and for the neighboring terms l = ±1 in the local sums (19)-(20) even when i = j, involve only interactions between differing interfaces, and thus may be replaced simply by substitution of the native quadrature rule (51) for the sources, and evaluation at discrete target nodes, as above. This method of discretizing an integral operator is called Nyström's method [35,Sec. 12.2]. This leaves only the self-interaction terms l = 0 in A i,i , which involve operators that are logarithmically singular. To achieve high-order accuracy for these operators, we use the standard "plain" Nyström matrix entries of the form A(z i m , z i j )w i j , for m, j = 1, . . . , N i (here A symbolizes a generic operator block) for entries far from the diagonal. Near-diagonal entries are adjusted by local Lagrange interpolation of the smooth density from the existing periodic trapezoid nodes onto a set of auxiliary nodes special to the logarithmic singularity, due to Alpert [2]; see [26,Sec. 4] for the full recipe. We use 30 auxiliary nodes per target node, which achieves high-order convergence with error O(N −16 i log N i ). With this Nyström matrix, the non-zero right-hand side terms in f become the samples at the first interface nodes z 1 j , and all the unknown vectors η i become samples of the densities τ i and σ i at the nodes of all interfaces.
The size of the resulting matrix A is N den -by-N den , where the total number of density unknowns is the factor of two coming from the two types of layer potential per interface.

Remark 2
The Alpert correction to the periodic trapezoid rule [2] is designed for closed curves, where the kernel and densities are periodic. However, our interfaces Γ i do not close on each other, moreover the solution density η i is quasi-periodic with Bloch phase α. This means that, when Γ i is smooth and hence has a single periodic quadrature rule, we must modify the Alpert correction entries in the northeast and southwest corners of the matrix, to account for the continuation of the interface into the next unit cell, and the phase factors α and α −1 . 2 Kress' periodic logarithmic correction would also be a slightly more accurate option [34] [26, Sec. 6]; however, we found that it was less convenient to adjust this scheme for quasi-periodic densities and open interface segments. (See Meier et al. [43] for an example of this; extra cut-off functions are required.) We will see that, due to the high-order convergence, the numbers of collocation nodes N i , n i,l , M w , and M can be small, of order a hundred, even for 10-digit accuracies. Note that for smooth interfaces the periodic trapezoid rule is in fact inaccurate for the interactions from neighboring interfaces, e.g. terms like S i Γi±d,Γi , due to "dangling" ends of these interfaces, but that this is handled by the periodizing scheme, retaining high-order convergence.
To summarize, the linear algebraic system which discretizes (49) has identical structure to (49), but with all of its blocks matrices as constructed above. We notate these blocks using standard (non-bold) font. It will now be rearranged in order to solve it in a fast direct fashion, exploiting its tridiagonal structure.

Rearrangement of equations, Schur complement, and fast solver
The discretized BIEs of the previous section require a few hundred to a couple of thousand unknowns per layer, to represent typical geometries as shown in Fig. 1 to accuracies of around 10 digits. The total number of unknowns includes the densities, proxy strengths, and scattered amplitudes, namely N = N den + IP + 2(2K + 1) .
We will find P ∼ 10 2 , thus there is an order of magnitude less proxy unknowns than density unknowns. In any case, a device with dozens or more layers leads to linear systems that are too large for direct O(N 3 ) inversion or solution. On the other hand, due to the proxy point (MFS) representation, the full linear system is exponentially ill-conditioned [7], so iterative solution of the full system is impossible. Therefore, for robustness, we describe in this section a direct solution technique, that will be "fast" (optimally scaling in I the number of layers), by exploiting the algebraic structure. This also will make it easier to solve for the technologically-important case of multiple incident angles θ inc at the same frequency ω.

Rearrangement
The first step is to rearrange the unknowns, in a form amenable to elimination of the proxy and scattering amplitude unknowns. We reorder our vector of all unknowns to be where Now similarly rearranging the (discretized) blocks of (49) we get the full linear system, Because the amplitudes a were separated into a U and a D then merged into c 1 and c I+1 to make x 1 and x I+1 , respectively, the first and last blocks of B, C and Q matrix had to be expanded to

Schur complements
The rearrangement in subsection 3.1 enables us to use I + 1 independent Schur complements to eliminate all the unknown vectors x i ; this corresponds to periodizing all of the Green's functions. For example, consider the equation in the first row, The first two equation rows starting the C block are With the first of these x 1 can be eliminated, and with the second x 2 can, giving where Q † i denotes the pseudo-inverse of the rectangular matrix Q i , for i = 1, . . . , I +1. Filling the matrix Q † i and then using matrix-matrix multiplication with the C blocks to fill matrices in (57) would lose accuracy, due to round-off error with the exponentially large matrix entries of the pseudo-inverses. To retain full accuracy, we do not in fact ever evaluate the pseudo-inverses. Rather, for product matrices such as X = Q † 1 C 1,1 appearing above, we solve the (ill-conditioned) linear system Q 1 X = C 1,1 with a standard backward-stable direct dense solver. 3 Similar computations eliminate x 1 , x 2 , and x 3 from the second equation to give By repeating the same computation for all the equations, all the x j are eliminated and a block tridiagonal system for η is obtained as with new interaction matrices for the density unknowns, Remark 3 The Schur complements (59)-(63) replacing matrix blocks A i,j by A i,j correspond to replacing each layer free-space Green's function G i (11) by an arbitrary quasi-periodic Green's function G i that obeys the layer's correct wall boundary conditions G i (x + d, y) = αG i (x, y), x ∈ L i , y ∈ Ω i . Crucially, since upward and downward radiation conditions are never imposed together in any single layer i, this is not the standard quasi-periodic G i QP of (12), which diverges at that layer's Wood anomalies. The G i selected by the backward-stable solves are always non-divergent; intuitively, since proxy strength vectors of small norm are possible, hence these are selected.

Block tridiagonal solve and evaluation of scattered wave
The block tridiagonal system (58) can be efficiently solved with a block LU decomposition [22,Sec. 4.5.1], at a cost of dense direct inversion of diagonal blocks. Since they derive from a second-kind integral equation, these diagonal blocks are all well conditioned, and no significant rounding error occurs when their inverses are multiplied as matrices. We write f i for the block vectors of the right-hand side of (58). The algorithm is initialized by settingf 1 = f 1 andÃ 1,1 = A 1,1 , then the forward sweep, for i = 2 to I in order, To save RAM, it is possible to haveÃ i,i overwrite A i,i . The backward sweep starts with solvingÃ I,I η I =f I , then for i = I − 1 down to 1 in order, solve for η i iñ If each N i ≈ N , for some constant N , the cost is O(N 3 I) due to two block inversions per layer. This is roughly I 2 times faster than naive inversion of the whole system. Note that the matrix filling time is only O(N 2 I), but dominates for the small N in our settings, due the large prefactor in evaluating special functions and applying Alpert corrections.
Once all of the density vectors η i are known, the proxy and scattering amplitude vectors x i are easily recovered by Here the products of the type Q † i C i,j can be reused from their prior computation in the Schur complements (59)-(63).
Finally, the scattered wave solution u i in each layer can be evaluated from their representations (16)-(18) by applying the interface quadrature rules to the single-and double-layer potentials, using the discrete density vectors η i . The evaluation involves only free-space Green's functions and monopole/dipole sources, which would be compatible with a standard FMM [18]. Above y = y U and below y = y D , the solution is evaluated from the Rayleigh-Bloch expansions (truncated versions of (7)-(8)), using the Bragg amplitudes a U and a D .

Accelerated sweep over multiple incident angles at one frequency
Modeling many real-world devices requires characterizing transmission and reflection over a wide range of incident angles θ inc at one frequency ω. In general, changing θ inc changes both the operators and the righthand side in (49), thus naively an independent solution is required for each angle, making the task very expensive. Here we show how to exploit two independent types of structure that speeds this up by typically an order of magnitude.
Firstly, we exploit the fact that the operators in (49), and hence the entire system matrix, depends on θ inc only through α in (2). Thus we may group together multiple incident angles that share α (as in [21]), and solve them together at a cost that is essentially the same as a single angle, using the precomputed inverse blocks in the tridiagonal solve of Sec. 3.3. Roughly there are k 1 d/π such angles, hence this is the speed-up factor. It grows in proportion to the incident wavenumber. For example, if k 1 = 40 and d = 1, the speedup is around a factor of 6. To set up a sweep of all incident angles, we choose a uniform grid in θ inc with spacing 2π/(ndk 1 ) for some integer n, insuring that only n independent α-value solves are needed.
Secondly, we exploit the fact that filling (rather than solving) the matrix blocks in (56) often accounts for the bulk of the solution time. Consider one of the integral operators used in the BIE matrix (defined in (19)), The integral part is independent of the incident angle or α. Therefore can be precomputed for l ∈ {−1, 0, 1} then used to assemble (D i V,W τ )(r) whenever a new α is given. Exactly the same argument applies to (S i V,W σ)(r), (T i V,W τ )(r), and (D i, * V,W σ)(r). Therefore, the A, B, C, and Q matrices can be assembled for any α simply by adding and subtracting precomputed integral operators; this speeds up the solution at multiple α values at the expense of using extra RAM.

Numerical Results
In all numerical examples, we let the first layer be "vacuum" with ε 1 = 1, set µ i = 1 for all layers, and choose periodicity d = 1. All the computations are conducted using MATLAB 2014a running on a workstation with two Intel Xeon E5-2687W processors (total 16 cores) with 256 GB memory and CentOS 6.5. For filling of the Nyström matrix blocks we use MPSpack [8], which has an interface to Fortran codes for Hankel function evaluations by Vladimir Rokhlin. The parfor command in MATLAB's Parallel Computing Toolbox is also used to fill the A and C matrices, using only 12 threads. In the following we present: numerical tests on convergence, scaling of timing and memory with frequency ω and the number of layers I, solution plots for a 1000-interface case, and accelerated computation of transmission and reflection spectra at multiple angles.

Remark 4
The Bragg coefficients {a U n } and {a D n } will be used as an independent measure of the accuracy of our numerical scheme based on conservation of the flux (energy) [41,53], namely This holds when all the material properties ε i and µ i are real. Therefore, we will define the relative flux error as

Convergence
We study the convergence of the scattered field at a fixed point (0.15, 0.6) in the first layer, relative to its converged value, and convergence of flux errors defined in (70). The frequency ω = 10 corresponds to a period of about 1.6 wavelengths in the first layer, and larger numbers of wavelengths in the other layers.      Thirty sine interfaces with random ε i chosen between 1 and 2 are considered in Fig. 2. The actual geometry is depicted in the inset of Fig. 2(b). All the interfaces have a periodic trapezoid rule with the same number of quadrature points N i = N , i = 1, 2, · · · 30. Fig. 2(a) shows convergence of the scattered field error (blue square) and flux error (red circle) as a function of N : the expected 16th-order rate is observed, and the fact that the two types of error are essentially equivalent. The same convergence test is conducted as a function of the number of proxy points P per layer, in Fig. 2(b). This time super-algebraic convergence is observed. The fact that 10-digit accuracy results with only N = P = 50 is a testament to the extremely rapid convergence of the method. Secondly, in order to study convergence for non-smooth interfaces, every other sine interface is replaced by a rectangular-ridge interface consisting of five line segments (see the inset in Fig. 3(c)), with Kress grading parameter q = 6. We study the convergence of the two types of interfaces independently. Fig. 3(a) shows convergence in N on sine interfaces, where N i = N , i = 1, 3, · · · , 29, while the quadrature points on rectangle is fixed at N i = 110 × 5, i = 2, 4, · · · , 30 (i.e. n i,l = 110 for all l = 1, . . . , 5). It shows the same convergence as before. Then, N i is fixed at 70 on all sine interfaces (i = 1, 3, · · · , 29), and N i = N × 5 on all rectangle interfaces (i = 2, 4, · · · , 30), is increased. Both the scattered field and flux error appear to converge as O(N −6 ) in Fig. 3(b), the order expected from the q = 6 grading. Fig. 3(c) shows super-algebraic convergence in P .
These tests confirm that flux error is a good indicator of pointwise error in u; from now on we quote flux error.

Performance
Tables 1 and 2 present the CPU time, memory usage, and flux error for periodic dielectric structures with I = 1, 3, 10, 30, 100, 300, and 1000 mixed sine, triangle, and rectangle interfaces. The first is at frequency ω = 5 (0.8 wavelengths per period), the second ω = 40 (6.4 wavelengths per period). Triangle interfaces tables, the computation time to solve the matrix (Schur complement and block matrix solve) and memory usage increases linearly as number of interfaces increases as expected in Sec. 3.3. For 1000 interfaces, it took 151 sec and 634 sec for ω = 5 and 40, respectively; note that more quadrature nodes are needed to accurately discretize the more oscillatory functions for higher ω. The sizes of the full matrix (56) for 1000 interfaces are 468260 × 287922 and 832000 × 751762, respectively. At ω = 40 the structure is around 8000λ tall.
In order to present a performance of the numerical method in some extreme cases, three numerical examples are presented. First, we considered 100 interfaces chosen randomly as sine, triangle, or rectangle type, with random heights, phases, layer thicknesses (while preventing collisions), and layer permittivities; see Fig. 4(a). We chose a frequency corresponding to a period of 4.5 wavelengths in the top layer, and an incident angle making the top layer precisely at a Wood anomaly. Due to space limitations, the real part of total field in only the first 10 and last 10 layers (regions enclosed by rectangles in Fig.4(a)) is plotted in Fig. 4(b). Matrix filling time was 19 sec, the Schur complements 9 sec, and the tridiagonal solve 8 sec. We achieve 9-digit accuracy in flux; in general we find that the flux error is no worse at Wood anomalies than at other angles. The second and third examples are for 1000 interfaces. The real part of total field u + u inc in only the top 10 layers is displayed due to figure resolution limitations. Fig. 5(b) is the real part of the total field from the 1000-interface case in the last column of Table 2. The third example is presented to highlight the geometric flexibility: we considered 1000 interfaces consisting of seven complex interface shapes on top of 993 sine interfaces ( Fig. 6(a)). Here ε i , i = 2, . . . , 1000 are chosen randomly between 1 and 3. All other parameters are given in the figure caption. The total field is computed and displayed in Figure 6(b). The CPU time for the solve was about 400 sec with 7 × 10 −9 flux error; evaluation of the solution at 1 million points took a similar time.

Transmission and Reflection Spectrum
We now compute transmission and reflection spectra as a function of incident angle, −π < θ inc < 0, and benchmark the acceleration technique of Sec. 3.4. The 30-interface structure shown in Fig. 7(a) is used, the  Tables 1 and 2. First, ω = 2 (0.3 wavelengths per period) and periodic permittivities ε = {1, 4, 1, 4, · · · , 4, 1} are considered. The spectrum clearly shows Bragg mirror or Fabry-Perot characteristics (there are ranges of incident angles that have total reflection), and symmetry, in Fig. 7(b). Because the wavelength is larger than the geometric features, the interface shape does not play an important role in determining the scattering. However, when ε i are set to random numbers between 1 and 4, the total reflection regime disappears in Fig. 7(c). Since ωd/π < 1, essentially no benefit comes from multiple angles sharing α values, but precomputation of matrix blocks does help. For 200 incident angles, the computation took 50 sec to produce with acceleration but 352 sec without, a 7× speed-up.
The frequency is now increased to ω = 10 (1.6 wavelengths per period), with 641 incident angles. Regardless of the ε distribution, the transmission-reflection spectrum behaves very abruptly in both periodic ( Fig. 7(d)) and random (Fig. 7(e)) ε cases. Also notice that the symmetry of the spectrum is broken because the layered structure has rectangle and triangle shaped interfaces which are now resolved by the wavelength. The computation took 121 sec with acceleration, and 1619 sec without, a speed-up of 13×. Finally, at ω = 20 (3.2 wavelengths per period), using 1279 incident angles, the speed-up was about 25× (we don't show these spectra since they do not show any new phenomena). This is consistent with the acceleration factor growing linearly with ω. Thus for ω = 40 as in Table 2 a speedup of 50× is expected.

Conclusion
We presented a new robust and fast integral equation method for 2D scattering from a periodic dielectric grating with an arbitrary number of layers of general shape. There are three main features: (1) The computational cost of the new method scales optimally (linearly) in the number of layers, allowing 1000layer structures to be solved rapidly. (2) The method is stable for all scattering parameters including Wood anomalies, since it is based on free-space rather than quasi-periodic Green's functions. (3) The periodizing scheme is simple, high-order accurate, and largely supersedes the Sommerfeld integral method of the second author and Greengard in [10]. This solver is expected to be useful in variety of wave applications in engineering and experimental physics, including the high accuracy modeling and optimization of optical, electromagnetic, and acoustic devices, and meta-materials.
There are several natural extensions of this work. Allowing dielectric inclusions in the layers (as in [36]), or material triple-junctions (e.g. incorporating robust representations of Lee-Greengard as in [23]), is simply a matter of bookkeeping, as long as the number of unknowns per layer remains small (e.g. less than 10 4 ). Since we use high-order quadrature schemes, this would allow significantly more complex unit cell shapes than presented here. Beyond this, an iterative FMM solution of the combined system would be appropriate when the number of layers is small [36]; for robustness with many layers a hierarchical fast direct solver such as in [21,23] could be used in each layer, combined with our tridiagonal block solve. Our scheme generalizes to 3D without needing new ideas, given a surface quadrature for the integral operators. However, the number of unknowns per layer could then easily exceed 10 4 , demanding something more elaborate than dense direct linear algebra within each layer.
For a production code, matrix filling and evaluation should be implemented in C or Fortran, and a parallel implementation would allow more simultaneous filling of matrix blocks, as well exploiting a parallel tridiagonal solve. In terms of analysis, an extension of the rigorous framework of free-space integral equations to include the presented periodizing scheme is needed.
The MATLAB codes which implement the methods of this paper and generate some of the figures can be downloaded from http://math.dartmouth.edu/∼mhcho/software These rely on layer-potential quadrature codes in the MPSpack toolbox by the second author, which can be downloaded from http://code.google.com/p/mpspack