An FFT-Accelerated Boundary Integral Equation Method for Wave Scattering by Smooth Surfaces in Three Dimensions

Abstract: For wave scattering by axisymmetric surfaces, the fast Fourier transform (FFT) method provides an effective tool to accelerate standard boundary integral equation (BIE) solvers. Surface integral equations can be decoupled into a series of curve integral equations on the generating curve, due to the convolution-like integral operators. The Fourier coefficients of the three-dimensional fundamental kernels can be rapidly computed through three-term recurrence relations based on Miller's algorithm. Such well-established techniques break down for nonaxisymmetric surfaces. This paper proposes a novel FFT-accelerated boundary integral method for wave scattering by smooth surfaces of arbitrary shapes. The Fourier coefficients of the singular kernels now satisfy higher-order recurrence relations. Although they can be solved with an optimal linear complexity by the standard Olver's algorithm, it turns out that a singularity swapping approach, that rewrites each kernel as the product of a smooth function and an axisymmetric-related singular factor, is realistically much faster. Consequently, Miller's algorithm together with the standard FFT convolution yields an ${\cal O}(M\log M)$ approach for evaluating the ${\cal O}(M)$ Fourier modes of the kernels, attaining exactly the same order of complexity for axisymmetric surfaces! With such FFT-based efficient procedures, we rewrite the surface integral equations in terms of ${\cal O}(M)$ weakly singular curve integrals, discretize them by panel-based generalized Gaussian quadratures, and obtain highly accurate linear systems to approximate the wavefields. Extensive numerical experiments are carried out to demonstrate the effectiveness of the new approach.

• The paper develops FFT-accelerated Fourier-coefficient algorithms for Helmholtz layer-potential kernels on non-axisymmetric surfaces, restoring O(M log M) evaluation complexity.
• The method enables spectrally accurate Nyström discretizations of single-, double-, and adjoint double-layer operators, achieving errors as low as 10^-12 to 10^-13 for tested analytic geometries.
• The approach handles ellipsoids, wavy and twisted tori, and other closed surfaces with little timing overhead versus axisymmetric cases, while leaving hypersingular operators and fast iterative solves for future work.

This paper develops an FFT-accelerated boundary integral equation (BIE) method for three-dimensional acoustic scattering by smooth surfaces of arbitrary shape, extending a framework previously restricted to axisymmetric geometries (2608.16208). The central contribution is a fast algorithm for evaluating the azimuthal Fourier coefficients of the Helmholtz layer-potential kernels on non-axisymmetric surfaces, which restores the $\mathcal{O}(M\log M)$ complexity achieved on axisymmetric surfaces and enables spectrally accurate Nyström discretizations of the single-layer, double-layer, and adjoint double-layer operators.

Problem setting and motivation

The authors consider interior and exterior Dirichlet and Neumann problems for the Helmholtz equation $(\Delta + k^2)u = 0$ posed on domains separated by a closed smooth surface $\Gamma$ parameterized as $\mathbf{r}(t,\theta)$ over $[0,L]\times[0,2\pi]$, with positive Jacobian everywhere. The four boundary value problems are reformulated via standard layer potentials: single-layer representations yield first-kind equations involving $\mathcal{S}$ for Neumann-type data, while double-layer representations give second-kind Fredholm equations $(\mathcal{K} \mp \mathcal{I})[\psi] = g$ for Dirichlet problems; the adjoint double-layer operator $\mathcal{K}'$ handles Neumann problems. The hypersingular operator $\mathcal{T}$ is explicitly deferred to future work.

The motivating observation is that prior FFT-accelerated schemes for axisymmetric surfaces exploit the fact that $\rho^2$ depends only on $(\Delta + k^2)u = 0$0, containing exactly Fourier modes $(\Delta + k^2)u = 0$1 and $(\Delta + k^2)u = 0$2. This yields a three-term recurrence for the kernel Fourier modes, solvable by Miller's algorithm in optimal complexity. For general surfaces this structure collapses, since $(\Delta + k^2)u = 0$3 acquires an infinite Fourier expansion in $(\Delta + k^2)u = 0$4.

High-order discretization of the single-layer operator

Discretization proceeds by splitting the surface integral into a $(\Delta + k^2)u = 0$5-integral of $(\Delta + k^2)u = 0$6, itself an azimuthal integral. A lemma establishes that for analytic surfaces, $(\Delta + k^2)u = 0$7 is smooth for $(\Delta + k^2)u = 0$8 and at most logarithmically singular at $(\Delta + k^2)u = 0$9, proved via a local Taylor expansion showing $\Gamma$0 near the singular point. Well-separated panels are handled by direct Gauss-Legendre quadrature in $\Gamma$1 combined with the trapezoidal rule in $\Gamma$2; adjacent and self panels use panel-based generalized Gaussian quadrature (GGQ) [Bremer–Gimbutas–Rokhlin], with 16 moving nodes per target point rather than fixed nodes, reducing cost relative to the earlier axisymmetric scheme.

For close interactions, where the trapezoidal rule would require prohibitively large $\Gamma$3, the density is expanded in a truncated Fourier series of length $\Gamma$4, reducing the problem to evaluating $\Gamma$5 Fourier coefficients $\Gamma$6. The paper presents two algorithms:

A linear-complexity recurrence approach: integrating the total derivative of $\Gamma$7 yields a convolution identity $\Gamma$8. Truncation at $\Gamma$9 modes—chosen adaptively so that tail coefficients fall below a threshold such as $\mathbf{r}(t,\theta)$0—produces a $\mathbf{r}(t,\theta)$1-term banded linear recurrence. With the outermost $\mathbf{r}(t,\theta)$2 coefficients precomputed by adaptive Gauss-Legendre quadrature, the remaining unknowns solve a banded Hermitian system; a sufficient condition for strict diagonal dominance is $\mathbf{r}(t,\theta)$3. The axisymmetric case recovers $\mathbf{r}(t,\theta)$4 and the classical three-term recurrence, confirming consistency with [Lai–O'Neil].

A singularity-swapping approach: the nearest complex root $\mathbf{r}(t,\theta)$5 ($\mathbf{r}(t,\theta)$6) of $\mathbf{r}(t,\theta)$7 is located by Newton iteration, and $\mathbf{r}(t,\theta)$8 is factored as $\mathbf{r}(t,\theta)$9, where $[0,L]\times[0,2\pi]$0 is holomorphic at the removed singularity pair and admits a short Fourier expansion of length independent of $[0,L]\times[0,2\pi]$1. The coefficients of $[0,L]\times[0,2\pi]$2 are scaled half-integer Legendre functions $[0,L]\times[0,2\pi]$3, evaluated stably through a difference-recurrence, so all required coefficients follow from a single zero-padded FFT convolution in $[0,L]\times[0,2\pi]$4 operations.

Empirically, the $[0,L]\times[0,2\pi]$5 route dominates even for large $[0,L]\times[0,2\pi]$6: reported CPU times are nearly identical across a sphere, ellipsoid, and twisted torus, meaning loss of axial symmetry incurs no appreciable overhead. Since the final assembly already carries an $[0,L]\times[0,2\pi]$7 factor, the nominally suboptimal method does not degrade overall complexity. Matrix construction splits into a preparation stage ($[0,L]\times[0,2\pi]$8) and a correction stage of complexity $[0,L]\times[0,2\pi]$9, with timings consistent with these rates.

Extension to double-layer operators

For $\mathcal{S}$0 and $\mathcal{S}$1, the kernels decompose into weakly singular terms built from $\mathcal{S}$2, $\mathcal{S}$3, multiplied by bounded analytic functions $\mathcal{S}$4 and $\mathcal{S}$5. The coefficient $\mathcal{S}$6 follows directly from convolutions of $\mathcal{S}$7 with $\mathcal{S}$8. For the more singular $\mathcal{S}$9, the paper derives a forward recurrence for $(\mathcal{K} \mp \mathcal{I})[\psi] = g$0 with an explicit initialization involving complete elliptic integrals of the second kind; a translated-singularity variant stabilizes the case of small $(\mathcal{K} \mp \mathcal{I})[\psi] = g$1 (threshold $(\mathcal{K} \mp \mathcal{I})[\psi] = g$2), following recent work on stabilizing singularity swapping. Cancellation errors in evaluating $(\mathcal{K} \mp \mathcal{I})[\psi] = g$3 are avoided using integral representations along the segment between source and target points, evaluated with few-point Gauss rules.

Numerical validation

Four geometries—an ellipsoid, a non-axisymmetric closed surface, a wiggly torus, and a twisted torus—are tested against manufactured solutions with exact field $(\mathcal{K} \mp \mathcal{I})[\psi] = g$4, wavenumbers up to $(\mathcal{K} \mp \mathcal{I})[\psi] = g$5. Across all examples, discretization errors below $(\mathcal{K} \mp \mathcal{I})[\psi] = g$6 are attained, reaching roughly $(\mathcal{K} \mp \mathcal{I})[\psi] = g$7 to $(\mathcal{K} \mp \mathcal{I})[\psi] = g$8 at low wavenumber; convergence curves show exponential decay in both the number of azimuthal modes $(\mathcal{K} \mp \mathcal{I})[\psi] = g$9 and the panel count $\mathcal{K}'$0, consistent with spectral accuracy for analytic data. The most demanding twisted-torus geometry requires $\mathcal{K}'$1 up to $\mathcal{K}'$2, still yielding errors around $\mathcal{K}'$3.

Limitations and open questions

The paper is candid about scope: the method assumes analytic surface parameterizations, so the exponential convergence guarantees do not extend to merely smooth or nonsmooth boundaries; the solvability condition on $\mathcal{K}'$4 for the banded system is sufficient rather than necessary; the hypersingular operator $\mathcal{K}'$5 remains undiscussed; and the resulting dense linear systems are solved directly, without a fast iterative solver—the existence of low-rank structure among curve operators coupling distinct Fourier indices (which vanishes identically for axisymmetric surfaces) is left as an open question. Electromagnetic and elastic extensions are likewise deferred.

Conclusion

The paper demonstrates that singularity swapping converts the higher-order recurrence structure of non-axisymmetric kernel Fourier coefficients into a form compatible with Miller-type stable evaluation plus one FFT convolution, achieving the same asymptotic per-mode cost as the axisymmetric theory and enabling high-order Nyström BIE solvers on general smooth closed surfaces. The combination of rigorous recurrence derivations, stabilized small-$\mathcal{K}'$6 handling, and extensive numerical verification makes this a substantive advance in fast high-order BIE discretization for 3D acoustic scattering.