Finite-Slab Multipole Methods

• Finite-Slab Multipole methods are specialized multipole-expansion techniques that retain finite-edge and finite-thickness effects in periodic acoustic, layered media, and gravitational domains.
• They enable efficient simulation of reaction fields and long-range interactions by leveraging block Toeplitz structures, equivalent polarization sources, and FFT-accelerated translation operators.
• The approaches optimize computational complexity and accuracy through exponential convergence bounds, tailored boundary treatments, and adaptations for acoustics, electromagnetics, and potential problems.

Searching arXiv for recent and foundational papers on finite-slab multipole methods and closely related layered-media multipole/FMM formulations. Finite-slab multipole denotes a family of multipole-expansion techniques specialized to geometries that are finite in a structurally significant sense: a bounded number of unit-cell replicas in periodic acoustic structures, a finite-thickness layer bounded by planar interfaces in layered-media potential theory and electromagnetics, or a bounded computational domain whose exterior field is represented by multipoles on an artificial boundary. Across these settings, the common objective is to preserve finite-edge or finite-thickness effects while accelerating long-range interactions through multipole expansions, translation operators, and structure-exploiting reductions such as block Toeplitz convolution or spherical-harmonic boundary operators (Jelich et al., 2022, Wang et al., 2020, Yuan et al., 24 Jul 2025, Yu et al., 12 Feb 2026).

1. Scope of the term

In the literature considered here, the term appears in three technically distinct but related senses. In acoustic boundary-element analysis, finite-slab multipole refers to finite periodic structures composed of a bounded number of unit-cell replicas arranged along one or more periodic directions, with explicit treatment of edge cells and aperiodic excitation (Jelich et al., 2022). In layered-media analysis for the Laplace and Maxwell equations, it refers to finite-thickness slabs bounded by planar interfaces, where multipole and local expansions are built for reaction fields generated by multiple reflections and transmissions through the slab (Wang et al., 2020, Yuan et al., 24 Jul 2025). In finite-element gravitation, the phrase is used for a bounded spherical buffer layer together with a multipole boundary treatment that eliminates the infinite exterior by enforcing the Neumann data implied by a truncated exterior expansion (Yu et al., 12 Feb 2026).

Setting	Governing field problem	Structural device
Finite periodic acoustics	Helmholtz equation with admittance boundaries	One FMM box per unit cell; block Toeplitz M2L structure
Layered finite slab	Layered Laplace or Maxwell Green’s functions	Equivalent polarization sources and reaction-field Sommerfeld integrals
Truncated gravitational exterior	Poisson equation on a bounded mesh	Spherical-harmonic multipole boundary flux

These variants share a methodological pattern. The geometry is not treated as infinitely periodic or fully unbounded in the computational sense. Instead, finite structural information is retained explicitly, while the long-range part is compressed through expansions whose algebra follows either free-space solid harmonics or spherical harmonics, and whose medium dependence is confined to reaction densities or boundary coefficients. This suggests that finite-slab multipole is better understood as a structural specialization of multipole methodology than as a single fixed algorithm.

2. Finite periodic acoustic structures

For linear acoustics with time-harmonic dependence $e^{-i\omega t}$ in a homogeneous fluid of density $\rho$ and sound speed $c$, the pressure satisfies

$$\nabla^2 p(\mathbf{x}) + k^2 p(\mathbf{x}) = 0,\qquad k=\omega/c,$$

together with an admittance boundary condition on $\Gamma$,

$$\frac{\partial p(\mathbf{x})}{\partial n(\mathbf{x})}=i\omega\rho\, v_\mathrm{f}(\mathbf{x}),\qquad i\omega\rho\, v_\mathrm{f}(\mathbf{x})=ik\,\beta(\mathbf{x})\,p(\mathbf{x}).$$

To ensure uniqueness, the acoustic finite-periodic formulation uses the combined-field Burton–Miller boundary integral equation with coupling parameter $\alpha$ satisfying $\mathrm{Im}(\alpha)\neq 0$, with optimal $\alpha=-i/k$ in the reported setting. After collocation BEM discretization with geometric quadrilateral surface elements and discontinuous Lagrange polynomials for $p$, the discrete system is

$\rho$0

Here $\rho$1 and $\rho$2 are dense, $\rho$3 is block diagonal, and the nodal surface pressures are the unknowns (Jelich et al., 2022).

The finite-slab acceleration aligns fast multipole boxes with unit cells. Well-separated unit cells interact through a truncated Helmholtz multipole series, while near-field interactions are computed directly. In single-level form, the decomposition is

$\rho$4

where $\rho$5 contains near-field interactions, $\rho$6 is particle-to-moment, $\rho$7 collects moment-to-local operators, and $\rho$8 is local-to-particle. The admissibility criterion is

$\rho$9

with $c$0 the box size. In multilevel form, the usual upward and downward shifts add moment-to-moment and local-to-local operators. Single-level complexity scales as $c$1 for a uniform grid of $c$2 boundary unknowns, while multilevel achieves $c$3 or $c$4 depending on tolerance (Jelich et al., 2022).

The defining finite-slab feature is the use of translational repetition without imposing field periodicity. When one box is chosen per unit cell, inter-cell operators depend only on lattice offset. This produces block Toeplitz structure in the far-field and, for regular near-field bands, also in $c$5. In one periodic direction, a Toeplitz matvec is equivalent to a circular convolution,

$c$6

and the multidimensional version follows by multi-indexing. Toeplitz–circulant embedding with size $c$7 in each periodic direction permits FFT-accelerated evaluation, reducing the cost of interaction families from $c$8 to $c$9 and reducing storage to the unique Toeplitz blocks (Jelich et al., 2022).

A common misconception is that Toeplitz structure implies periodic boundary conditions on the acoustic field. In the reported formulation, the Toeplitz representation maps offsets between boxes, not periodic boundary conditions on fields. Edge boxes simply contribute fewer neighbors in near-field bands, and edge effects are treated natively. For half-space problems with a ground plane, translation invariance is preserved when the periodic axes lie parallel to the plane; when an axis is perpendicular, the mirrored contribution is block Hankel, and a permutation converts the relevant operator to Toeplitz form so that the same FFT acceleration can be applied (Jelich et al., 2022).

3. Layered-media finite slabs for the Laplace equation

For the three-dimensional Laplace equation in layered media, the finite-slab configuration consists of two planar interfaces at $$\nabla^2 p(\mathbf{x}) + k^2 p(\mathbf{x}) = 0,\qquad k=\omega/c,$$0 and $$\nabla^2 p(\mathbf{x}) + k^2 p(\mathbf{x}) = 0,\qquad k=\omega/c,$$1 with $$\nabla^2 p(\mathbf{x}) + k^2 p(\mathbf{x}) = 0,\qquad k=\omega/c,$$2, creating a top half-space, a finite slab, and a bottom half-space. Each layer $$\nabla^2 p(\mathbf{x}) + k^2 p(\mathbf{x}) = 0,\qquad k=\omega/c,$$3 has isotropic permittivity $$\nabla^2 p(\mathbf{x}) + k^2 p(\mathbf{x}) = 0,\qquad k=\omega/c,$$4, and the layered Green’s function $$\nabla^2 p(\mathbf{x}) + k^2 p(\mathbf{x}) = 0,\qquad k=\omega/c,$$5 solves

$$\nabla^2 p(\mathbf{x}) + k^2 p(\mathbf{x}) = 0,\qquad k=\omega/c,$$6

subject to transmission conditions across the interfaces and decay as $$\nabla^2 p(\mathbf{x}) + k^2 p(\mathbf{x}) = 0,\qquad k=\omega/c,$$7. For $$\nabla^2 p(\mathbf{x}) + k^2 p(\mathbf{x}) = 0,\qquad k=\omega/c,$$8, the Green’s function decomposes as

$$\nabla^2 p(\mathbf{x}) + k^2 p(\mathbf{x}) = 0,\qquad k=\omega/c,$$9

and the reaction field is expressed as a sum of four Sommerfeld-type components $\Gamma$0, $\Gamma$1 (Wang et al., 2020).

The reaction components are written as

$\Gamma$2

where $\Gamma$3 and the four mappings $\Gamma$4 encode the interface geometry. For the three-layer slab, the slab thickness appears through

$\Gamma$5

inside the transfer-matrix construction of the reaction densities $\Gamma$6. A central theoretical result is that these densities are analytic and bounded in the right half complex plane $\Gamma$7 when the layer parameters are positive. That analyticity underlies the convergence theory for the slab-adapted multipole expansions (Wang et al., 2020).

The slab formulation introduces equivalent polarization source points $\Gamma$8 produced by reflecting the source depth across the interfaces nearest the target layer. With $\Gamma$9, all reaction components are rewritten as Sommerfeld integrals in Euclidean differences between targets and polarization sources. This leads to multipole expansions and local expansions with the same solid-harmonic algebra as free space:

$$\frac{\partial p(\mathbf{x})}{\partial n(\mathbf{x})}=i\omega\rho\, v_\mathrm{f}(\mathbf{x}),\qquad i\omega\rho\, v_\mathrm{f}(\mathbf{x})=ik\,\beta(\mathbf{x})\,p(\mathbf{x}).$$0

and

$$\frac{\partial p(\mathbf{x})}{\partial n(\mathbf{x})}=i\omega\rho\, v_\mathrm{f}(\mathbf{x}),\qquad i\omega\rho\, v_\mathrm{f}(\mathbf{x})=ik\,\beta(\mathbf{x})\,p(\mathbf{x}).$$1

The only modification relative to free space is the presence of the slab-dependent densities $$\frac{\partial p(\mathbf{x})}{\partial n(\mathbf{x})}=i\omega\rho\, v_\mathrm{f}(\mathbf{x}),\qquad i\omega\rho\, v_\mathrm{f}(\mathbf{x})=ik\,\beta(\mathbf{x})\,p(\mathbf{x}).$$2 inside the Sommerfeld integrals and the vertical reflection $$\frac{\partial p(\mathbf{x})}{\partial n(\mathbf{x})}=i\omega\rho\, v_\mathrm{f}(\mathbf{x}),\qquad i\omega\rho\, v_\mathrm{f}(\mathbf{x})=ik\,\beta(\mathbf{x})\,p(\mathbf{x}).$$3 when $$\frac{\partial p(\mathbf{x})}{\partial n(\mathbf{x})}=i\omega\rho\, v_\mathrm{f}(\mathbf{x}),\qquad i\omega\rho\, v_\mathrm{f}(\mathbf{x})=ik\,\beta(\mathbf{x})\,p(\mathbf{x}).$$4 (Wang et al., 2020).

Because the coefficients retain free-space solid-harmonic form, M2M and L2L translation formulas apply unchanged. M2L remains medium dependent through translation kernels $$\frac{\partial p(\mathbf{x})}{\partial n(\mathbf{x})}=i\omega\rho\, v_\mathrm{f}(\mathbf{x}),\qquad i\omega\rho\, v_\mathrm{f}(\mathbf{x})=ik\,\beta(\mathbf{x})\,p(\mathbf{x}).$$5 containing $$\frac{\partial p(\mathbf{x})}{\partial n(\mathbf{x})}=i\omega\rho\, v_\mathrm{f}(\mathbf{x}),\qquad i\omega\rho\, v_\mathrm{f}(\mathbf{x})=ik\,\beta(\mathbf{x})\,p(\mathbf{x}).$$6 and the interface geometry. The convergence theorem gives explicit remainder bounds of the form

$$\frac{\partial p(\mathbf{x})}{\partial n(\mathbf{x})}=i\omega\rho\, v_\mathrm{f}(\mathbf{x}),\qquad i\omega\rho\, v_\mathrm{f}(\mathbf{x})=ik\,\beta(\mathbf{x})\,p(\mathbf{x}).$$7

with geometric ratios determined by box radii, target distances, and well-separation. The finite-slab interpretation is precise: slab thickness modifies the densities $$\frac{\partial p(\mathbf{x})}{\partial n(\mathbf{x})}=i\omega\rho\, v_\mathrm{f}(\mathbf{x}),\qquad i\omega\rho\, v_\mathrm{f}(\mathbf{x})=ik\,\beta(\mathbf{x})\,p(\mathbf{x}).$$8 but does not change the Euclidean geometry of multipole and local expansion convergence once equivalent polarization sources are introduced. As $$\frac{\partial p(\mathbf{x})}{\partial n(\mathbf{x})}=i\omega\rho\, v_\mathrm{f}(\mathbf{x}),\qquad i\omega\rho\, v_\mathrm{f}(\mathbf{x})=ik\,\beta(\mathbf{x})\,p(\mathbf{x}).$$9, $\alpha$0 and the slab reduces to a half-space (Wang et al., 2020).

4. Maxwell finite-thickness layered media

In the electromagnetic setting, finite-slab multipole methods are formulated for Maxwell’s equations in three-dimensional layered media under time dependence $\alpha$1. Using the magnetic vector potential $\alpha$2 under the Lorenz gauge,

$\alpha$3

the fields satisfy

$\alpha$4

and $\alpha$5 solves a vector Helmholtz equation. The layered dyadic Green’s function is represented using three scalar Helmholtz layered Green’s functions, and the interface-induced reaction field is written through unified Sommerfeld integrals (Yuan et al., 24 Jul 2025).

The spectral-domain electric dyadic is decomposed into a free-space part and reaction components indexed by whether the field leaves the source upward or downward and arrives upward or downward. Each reaction component has a unified form

$\alpha$6

where the scalar functions $\alpha$7 and $\alpha$8 are generalized reflection/transmission coefficients obtained from classic $\alpha$9 recurrences for multilayer stacks. Inverse Fourier transform yields eight Sommerfeld-type integrals grouped by five angular harmonics $\mathrm{Im}(\alpha)\neq 0$0 with $\mathrm{Im}(\alpha)\neq 0$1 (Yuan et al., 24 Jul 2025).

A central construction is the use of equivalent polarization coordinates for sources together with effective locations for targets. These are defined by interface reflections $\mathrm{Im}(\alpha)\neq 0$2 and are chosen so that the asymptotic exponential factor matches the minimal vertical transmission distance for each reaction component. The layered far-field criterion is then built from the actual transmission distance

$\mathrm{Im}(\alpha)\neq 0$3

This distance is larger than or equal to prior criteria based solely on image sources. The reported implication is algorithmic: interaction lists built from $\mathrm{Im}(\alpha)\neq 0$4 improve conditioning of multipole and local expansions and reduce near-field workload along interfaces (Yuan et al., 24 Jul 2025).

The reaction-field multipole expansion and local expansion retain the same spherical-harmonic structure as free-space Helmholtz and Maxwell FMMs. M2M and L2L are therefore identical to free space, while M2L is slab dependent and assembled from master Sommerfeld integrals. To reduce the $\mathrm{Im}(\alpha)\neq 0$5 cost of tabulating M2L integrals, the method expands products of analytically continued associated Legendre functions in Chebyshev polynomials. The resulting M2L tables require only $\mathrm{Im}(\alpha)\neq 0$6 one-dimensional Sommerfeld-type integrals per table instead of $\mathrm{Im}(\alpha)\neq 0$7. The full layered-media FMM has overall cost $\mathrm{Im}(\alpha)\neq 0$8 for low-to-moderate frequencies and modest numbers of layers, with memory bounded by a finite number of geometry cases per level (Yuan et al., 24 Jul 2025).

The finite-slab aspect is not an explicit image summation. For finite-thickness slabs, multiple reflections are embedded in the recurrence-defined coefficients $\mathrm{Im}(\alpha)\neq 0$9 and $\alpha=-i/k$0, while the vertical exponential factors encode only the minimal transmission path. This distinguishes the method from direct image-series constructions and from free-space FMMs that do not account for interface-induced reaction fields (Yuan et al., 24 Jul 2025).

5. Multipole boundary treatment on truncated domains

A different finite-slab usage appears in finite-element gravitation, where the Poisson equation for the gravitational potential is posed on $\alpha=-i/k$1,

$\alpha=-i/k$2

but computations are performed on a bounded spherical domain $\alpha=-i/k$3 containing the body $\alpha=-i/k$4. The weak form on $\alpha=-i/k$5 is

$\alpha=-i/k$6

The exterior field is then eliminated either through a Dirichlet-to-Neumann map or through a multipole expansion on the artificial boundary (Yu et al., 12 Feb 2026).

For $\alpha=-i/k$7, the harmonic exterior potential admits the spherical-harmonic expansion

$\alpha=-i/k$8

with moments

$\alpha=-i/k$9

Using the relation $p$0, the boundary flux on $p$1 is written directly in terms of interior density moments. Substitution into the weak form yields the symmetric density-to-boundary-flux bilinear form

$p$2

The special case $p$3 recovers the global mass correction (Yu et al., 12 Feb 2026).

This multipole boundary treatment differs structurally from the acoustic and layered-media FMM formulations. It is not a hierarchical many-body accelerator; rather, it is an exact or truncated representation of the infinite exterior on a spherical truncation boundary. In MFEM, the multipole operator factors as $p$4, where $p$5 maps interior density degrees of freedom to multipole coefficients and $p$6 maps boundary test functions to spherical coefficients. For static problems, the operator is applied once to construct a corrected right-hand side, after which the pure Neumann diffusion system is solved with AMG preconditioning. For linearized perturbations with $p$7, the same strategy produces a vector-to-boundary-flux operator involving $p$8 (Yu et al., 12 Feb 2026).

The same work also treats the exact DtN map on a spherical boundary, where each mode satisfies

$p$9

The DtN operator is diagonal in spherical harmonics and is implemented matrix-free as $\rho$00 on boundary degrees of freedom. The comparison clarifies the role of multipole truncation in this setting: the multipole method moves work to a one-time right-hand-side build, whereas DtN applies a low-rank boundary operator at every matrix-vector product (Yu et al., 12 Feb 2026).

6. Numerical behavior, advantages, and limitations

The acoustic finite-periodic formulation is validated on a $\rho$01 array of sound-hard spheres of radius $\rho$02 and spacing $\rho$03 in air with $\rho$04 and $\rho$05. Each sphere uses $\rho$06 quadrilateral elements, for $\rho$07 degrees of freedom. At $\rho$08, FMM with $\rho$09, FMPBEM with $\rho$10, and FMPBEM2 with $\rho$11 for inter-cell and $\rho$12 for intra-cell expansions achieve $\rho$13 matvec error $\rho$14 versus PBEM. Scaling in the number of unit cells shows $\rho$15 assembly and matvec behavior for FMPBEM, while memory is asymptotically $\rho$16; the largest reported refinement, with up to $\rho$17 million total degrees of freedom, assembled in approximately $\rho$18 and performed a matvec in approximately $\rho$19 on a 6-core desktop with $\rho$20 RAM (Jelich et al., 2022).

The same acoustic study includes a half-space sound-barrier comparison. A wall barrier periodic model gives insertion loss closely matching a full-scale PBEM model despite omission of top and side surfaces in the periodic model; the reported insertion loss stays above $\rho$21 up to $\rho$22, with peak $\rho$23 at $\rho$24 and minimum $\rho$25 at $\rho$26. For a sonic crystal barrier of cylinders, insertion loss rises from approximately $\rho$27 below $\rho$28 to a peak $\rho$29 near $\rho$30, close to the first Bragg condition. For a barrier of c-shaped Helmholtz resonators, a 2D FEM estimate gives a Helmholtz resonance near $\rho$31, while FMPBEM reports insertion loss $\rho$32 at $\rho$33 and a second band of insertion loss at least $\rho$34 between $\rho$35 and $\rho$36. These results show that finite-slab formulations can capture both edge-sensitive barrier performance and resonance or Bragg phenomena without imposing infinite periodicity (Jelich et al., 2022).

For the Laplace layered-slab problem, the key result is theoretical rather than benchmark centered: the reaction densities are analytic and bounded in $\rho$37, and the resulting multipole, local, and translation expansions satisfy explicit exponential convergence bounds. The error constants depend on $\rho$38, source or target box radii, and geometric margins, while the convergence rate is set by ratios such as $\rho$39, $\rho$40, and $\rho$41. A recurring misconception is that slab thickness changes the algebra of the multipole method. In the reported analysis, slab thickness enters only through the reaction densities via $\rho$42; the solid-harmonic translation algebra remains that of free space (Wang et al., 2020).

For Maxwell layered media, numerical experiments on two- and three-layer stacks show spectral convergence of both free-space and reaction components as $\rho$43 increases. The overall algorithm exhibits $\rho$44 complexity, and reaction-field evaluations are significantly faster than free-space parts because actual transmission distance produces enhanced separations across interfaces and very few near-field pairs. One reported example gives, for two layers and $\rho$45, reaction fields for $\rho$46 particles per layer in approximately $\rho$47 per layer versus approximately $\rho$48–$\rho$49 for the free-space parts on a single CPU core. The stated regime of robustness is low to moderate frequency; at higher frequencies, $\rho$50 must increase and alternative acceleration such as plane-wave forms may be needed (Yuan et al., 24 Jul 2025).

For gravitational FEM on truncated spherical domains, naive Dirichlet truncation can still be acceptable if a very large buffer domain is affordable and coarsening is used in the exterior mesh. The offset-sphere benchmark requires $\rho$51 to reach $\rho$52, with reference solver time $\rho$53. On the same class of problems, both DtN and multipole boundary treatments at $\rho$54 produce rapid error decay with $\rho$55; by $\rho$56, DtN errors are dominated by finite-element discretization, while multipole and DtN have very similar accuracy for a given $\rho$57. Solver times are nearly independent of $\rho$58, but multipole assembly grows rapidly with $\rho$59 and can exceed diffusion assembly by approximately $\rho$60 at $\rho$61. This comparison places the finite-slab multipole boundary treatment in a broader numerical context: it is advantageous when boundary-condition accuracy dominates, but its cost structure differs from both DtN and hierarchical FMM variants (Yu et al., 12 Feb 2026).

Across these formulations, the main distinctions from infinite-periodic, purely free-space, or naive truncation methods are consistent. Finite-slab multipole methods retain finite edges, finite thickness, or finite truncation explicitly; they exploit repeated or separable structure where it exists; and they localize the effect of the surrounding medium or exterior either in reaction densities, translation kernels, or boundary spherical-harmonic operators. Their limitations are equally setting dependent: linear acoustics and homogeneous backgrounds in the acoustic formulation, positive layered parameters and planar interfaces in the Laplace theory, low-to-moderate frequencies in the Maxwell layered FMM, and spherical truncation boundaries for diagonal DtN structure in the gravitation setting (Jelich et al., 2022, Wang et al., 2020, Yuan et al., 24 Jul 2025, Yu et al., 12 Feb 2026).