Alex H. Barnett

Alex H. Barnett, Gary Marple, Shravan Veerapaneni et al.

We present a spectrally-accurate scheme to turn a boundary integral formulation for an elliptic PDE on a single unit cell geometry into one for the fully periodic problem. Applications include computing the effective permeability of composite media (homogenization), and microfluidic chip design. Our basic idea is to exploit a small least squares solve to apply periodicity without ever handling periodic Green's functions. We exhibit fast solvers for the two-dimensional (2D) doubly-periodic Neumann Laplace problem (flow around insulators), and Stokes non-slip fluid flow problem, that for inclusions with smooth boundaries achieve 12-digit accuracy, and can handle thousands of inclusions per unit cell. We split the infinite sum over the lattice of images into a directly-summed "near" part plus a small number of auxiliary sources which represent the (smooth) remaining "far" contribution. Applying physical boundary conditions on the unit cell walls gives an expanded linear system, which, after a rank-1 or rank-3 correction and a Schur complement, leaves a well-conditioned square system which can be solved iteratively using fast multipole acceleration plus a low-rank term. We are rather explicit about the consistency and nullspaces of both the continuous and discretized problems. The scheme is simple (no lattice sums, Ewald methods, nor particle meshes are required), allows adaptivity, and is essentially dimension- and PDE-independent, so would generalize without fuss to 3D and to other non-oscillatory elliptic problems such as elastostatics. We incorporate recently developed spectral quadratures that accurately handle close-to-touching geometries. We include many numerical examples, and provide a software implementation.

Alex H. Barnett, Leslie Greengard, Tom Hagstrom

We develop a high-order, explicit method for acoustic scattering in three space dimensions based on a combined-field time-domain integral equation. The spatial discretization, of Nyström type, uses Gaussian quadrature on panels combined with a special treatment of the weakly singular kernels arising in near-neighbor interactions. In time, a new class of convolution splines is used in a predictor-corrector algorithm. Experiments on a torus and a perturbed torus are used to explore the stability and accuracy of the proposed scheme. This involved around one thousand solver runs, at up to 8th order and up to around 20,000 spatial unknowns, demonstrating 5-9 digits of accuracy. In addition we show that parameters in the combined field formulation, chosen on the basis of analysis for the sphere and other convex scatterers, work well in these cases.

Yuxiang Liu, Alex H. Barnett

We present a high-order accurate boundary-based solver for three-dimensional (3D) frequency-domain scattering from a doubly-periodic grating of smooth axisymmetric sound-hard or transmission obstacles. We build the one-obstacle solution operator using separation into P azimuthal modes via the FFT, the method of fundamental solutions (with N proxy points lying on a curve), and dense direct least-squares solves; the effort is O(N^3P) with a small constant. Periodizing then combines fast multipole summation of nearest neighbors with an auxiliary global Helmholtz basis expansion to represent the distant contributions, and enforcing quasi-periodicity and radiation conditions on the unit cell walls. Eliminating the auxiliary coefficients, and preconditioning with the one-obstacle solution operator, leaves a well-conditioned square linear system that is solved iteratively. The solution time per incident wave is then O(NP) at fixed frequency. Our scheme avoids singular quadratures, periodic Green's functions, and lattice sums, and its convergence rate is unaffected by resonances within obstacles. We include numerical examples such as scattering from a grating of period 13 λ x 13λ of highly-resonant sound-hard "cups" each needing NP = 64800 surface unknowns, to 10-digit accuracy, in half an hour on a desktop.

Alex H. Barnett, Andrew Hassell

We present a new algorithm for numerical computation of large eigenvalues and associated eigenfunctions of the Dirichlet Laplacian in a smooth, star-shaped domain in $\mathbb{R}^d$, $d\ge 2$. Conventional boundary-based methods require a root-search in eigenfrequency $k$, hence take $O(N^3)$ effort per eigenpair found, using dense linear algebra, where $N=O(k^{d-1})$ is the number of unknowns required to discretize the boundary. Our method is O(N) faster, achieved by linearizing with respect to $k$ the spectrum of a weighted interior Neumann-to-Dirichlet (NtD) operator for the Helmholtz equation. Approximations $\hat{k}_j$ to the square-roots $k_j$ of all O(N) eigenvalues lying in $[k - ε, k]$, where $ε=O(1)$, are found with $O(N^3)$ effort. We prove an error estimate $$ |\hat k_j - k_j| \leq C \Big(\frac{ε^2}{k} + ε^3 \Big), $$ with $C$ independent of $k$. We present a higher-order variant with eigenvalue error scaling empirically as $O(ε^5)$ and eigenfunction error as $O(ε^3)$, the former improving upon the 'scaling method' of Vergini--Saraceno. For planar domains ($d=2$), with an assumption of absence of spectral concentration, we also prove rigorous error bounds that are close to those numerically observed. For $d=2$ we compute robustly the spectrum of the NtD operator via potential theory, Nyström discretization, and the Cayley transform. At high frequencies (400 wavelengths across), with eigenfrequency relative error $10^{-10}$, we show that the method is $10^3$ times faster than standard ones based upon a root-search.

Zijian Wang, Jeremy G. Hoskins, Manas Rachh et al.

In this paper, we study the scale-invariant quantity \[\mathcal{G}(Ω)=\frac{\|\partial_n u_1\|_{L^\infty(\partialΩ)}}{λ_1},\]where $u_1$ is the first $L^2$-normalized Dirichlet Laplace eigenfunction of a Euclidean domain $Ω$ and $λ_1$ is its eigenvalue. This is related to the peak boundary heat flux in the long time limit. For convex domains we prove that $\|\partial_n u_1\|_{L^\infty(\partialΩ)}$ is upper-bounded by a (domain-independent) constant multiple of $λ_1$. Using layer potentials, we derive shape-derivative formulae for efficient gradient computations. When combined with high-order Nyström discretization, a fast boundary integral equation solver, and eigenvalue rootfinding, this allows us to numerically optimize $\mathcal{G}$ over a class of rounded polygonal discretized domains. Based on extensive numerical experiments, we then conjecture that, over the set of convex domains, $\mathcal{G}$ is maximized by the semidisk, with the peak flux at the center of the diameter. To lend analytical support to this conjecture, we prove that the semidisk is a critical point of $\mathcal{G}$ under infinitesimal perturbations of its circular arc.

Nour G. Al Hassanieh, Leslie Greengard, Alex H. Barnett

We present a fast algorithm for evaluating the (non-smooth) solution of the free-space two-dimensional (2D) scalar wave equation with many point sources, each with a high-frequency band-limited time signature. Such an algorithm is key to an efficient time-domain scattering solver using spatially-discretized hyperbolic layer potentials. Given $M$ sources/targets and $N_t$ time steps, direct evaluation costs $O(M^2N_t^2)$, due to the history dependence. We develop a quasi-linear scaling algorithm that splits the solution at a given time into (a) a non-smooth time-local part, (b) a (smooth) near history involving sources up to ${\mathcal O}(1)$ domain traversal times into the past, plus (c) a (very smooth) far history comprising all waves emitted before the near history. The local part is computed directly via high-order quadrature. A naive spatial Fourier transform for (b) plus (c) would be both slowly converging and arbitrarily oscillatory as time progresses. Yet in (b) the oscillations are controlled, so we use the recent truncated windowed Fourier projection (TK-WFP) method to give rapid convergence. For (c) -- present due to the weak Huygens' principle -- we exploit a new large-time sum-of-exponentials approximation of the free-space wave kernel. Numerical examples with up to a million sources and targets, a domain of $300\times 300$ wavelengths, and 6-digit accuracy, show an acceleration of five orders of magnitude relative to direct evaluation.

Alex H. Barnett, Charles L. Epstein, Leslie Greengard et al.

We study the stability properties of explicit marching schemes for second-kind Volterra integral equations that arise when solving boundary value problems for the heat equation by means of potential theory. It is well known that explicit finite difference or finite element schemes for the heat equation are stable only if the time step $Δt$ is of the order $O(Δx^2)$, where $Δx$ is the finest spatial grid spacing. In contrast, for the Dirichlet and Neumann problems on the unit ball in all dimensions $d\ge 1$, we show that the simplest Volterra marching scheme, i.e., the forward Euler scheme, is unconditionally stable. Our proof is based on an explicit spectral radius bound of the marching matrix, leading to an estimate that an $L^2$-norm of the solution to the integral equation is bounded by $c_dT^{d/2}$ times the norm of the right hand side. For the Robin problem on the half space in any dimension, with constant Robin (heat transfer) coefficient $κ$, we exhibit a constant $C$ such that the forward Euler scheme is stable if $Δt < C/κ^2$, independent of any spatial discretization. This relies on new lower bounds on the spectrum of real symmetric Toeplitz matrices defined by convex sequences. Finally, we show that the forward Euler scheme is unconditionally stable for the Dirichlet problem on any smooth convex domain in any dimension, in $L^\infty$-norm.

Alex H. Barnett, Jeremy F. Magland, Leslie F. Greengard

We describe a suite of validation metrics that assess the credibility of a given automatic spike sorting algorithm applied to a given electrophysiological recording, when ground-truth is unavailable. By rerunning the spike sorter two or more times, the metrics measure stability under various perturbations consistent with variations in the data itself, making no assumptions about the noise model, nor about the internal workings of the sorting algorithm. Such stability is a prerequisite for reproducibility of results. We illustrate the metrics on standard sorting algorithms for both in vivo and ex vivo recordings. We believe that such metrics could reduce the significant human labor currently spent on validation, and should form an essential part of large-scale automated spike sorting and systematic benchmarking of algorithms.

Min Hyung Cho, Alex H. Barnett

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.

Alex H. Barnett, Bowei Wu, Shravan K. Veerapaneni

Dense particulate flow simulations using integral equation methods demand accurate evaluation of Stokes layer potentials on arbitrarily close interfaces. In this paper, we generalize techniques for close evaluation of Laplace double-layer potentials in J. Helsing and R. Ojala, J. Comput. Phys. 227 (2008) 2899-2921. We create a "globally compensated" trapezoid rule quadrature for the Laplace single-layer potential on the interior and exterior of smooth curves. This exploits a complex representation, a product quadrature (in the style of Kress) for the sawtooth function, careful attention to branch cuts, and second-kind barycentric-type formulae for Cauchy integrals and their derivatives. Upon this we build accurate single- and double-layer Stokes potential evaluators by expressing them in terms of Laplace potentials. We test their convergence for vesicle-vesicle interactions, for an extensive set of Laplace and Stokes problems, and when applying the system matrix in a boundary value problem solver in the exterior of multiple close-to-touching ellipses. We achieve typically 12 digits of accuracy using very small numbers of discretization nodes per curve. We provide documented codes for other researchers to use.

Alex H. Barnett, Leslie Greengard

In this paper, we consider band-structure calculations governed by the Helmholtz or Maxwell equations in piecewise homogeneous periodic materials. Methods based on boundary integral equations are natural in this context, since they discretize the interface alone and can achieve high order accuracy in complicated geometries. In order to handle the quasi-periodic conditions which are imposed on the unit cell, the free-space Green's function is typically replaced by its quasi-periodic cousin. Unfortunately, the quasi-periodic Green's function diverges for families of parameter values that correspond to resonances of the empty unit cell. Here, we bypass this problem by means of a new integral representation that relies on the free-space Green's function alone, adding auxiliary layer potentials on the boundary of the unit cell itself. An important aspect of our method is that by carefully including a few neighboring images, the densities may be kept smooth and convergence rapid. This framework results in an integral equation of the second kind, avoids spurious resonances, and achieves spectral accuracy. Because of our image structure, inclusions which intersect the unit cell walls may be handled easily and automatically. Our approach is compatible with fast-multipole acceleration, generalizes easily to three dimensions, and avoids the complication of divergent lattice sums.