Recursive Computation of the Multipole Expansions of Layer Potential Integrals over Simplices for Efficient Fast Multipole Accelerated Boundary Elements Nail A. Gumerov 1 , Shoken Kaneko and Ramani Duraiswami Department of Computer Science and Institute for Advanced Computer Studies University of Maryland, College Park, MD 20742 Abstract In boundary element methods (BEM) in R 3 , matrix elements and right hand sides are typically computed via analytical or numerical quadrature of the layer potential multiplied by some function over line, triangle and tetrahedral volume elements. When the problem size gets large, the resulting linear systems are often solved iteratively via Krylov subspace methods, with fast multipole methods (FMM) used to accelerate the matrix vector products needed. When FMM acceleration is used, most entries of the matrix never need be computed explicitly - they are only needed in terms of their contribution to the multipole expansion coefficients. We propose a new fast method - Quadrature to Expansion (Q2X) - for the analytical generation of the multipole expansion coefficients produced by the integral expressions for single and double layers on surface triangles; charge distributions over line segments and over tetrahedra in the volume; so that the overall method is well integrated into the FMM, with controlled error. The method is based on the O(1) per moment cost recursive computation of the moments. The method is developed for boundary element methods involving the Laplace Green’s function in R 3 . The derived recursions are first compared against classical quadrature algorithms, and then integrated into FMM accelerated boundary element and vortex element methods. Numerical tests are presented and discussed. Keywords: Fast Multipole Method, Boundary Element Method, Vortex Element Method, Multipole Expansions, Laplace and Poisson Equations, Special Function Recursions, Quadrature to Expansion 1. Introduction The fast multipole method (FMM) [5] is often used to accelerate the matrix vector products needed in the Krylov subspace based iterative methods for solution of the large linear systems that result from boundary element methods (BEM) (e.g., for a review see [16]). The surface is usually discretized in patches, which in the simplest case are plane triangles. For the lowest-order consistent approximation, the functions of interest are approximated via low order polynomials over these triangles, often constant or linear. For constant panel methods, the BEM matrix elements can be represented as integrals of the Green’s function and its derivatives over the triangles (“layer potentials”). Techniques for evaluating these layer potentials have been extensively developed including analytical methods, e.g., [15, 14, 20, 25, 13, 11]. In [13], the authors have developed a method for layer potential integral evaluation via dimensionality reduction and recursions for the Laplace and Helmholtz kernel for higher order elements. The conventional BEM solution proceeds by evaluating all matrix elements via quadrature over the surface elements and solving the linear system either via matrix decomposition or iteratively. When used 1 This author passed away during the submission of this paper. Preprint submitted to Elsevier April 6, 2023 arXiv:2107.10942v3 [math.NA] 5 Apr 2023