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.