pith. sign in

arxiv: 2604.14797 · v2 · submitted 2026-04-16 · 🧮 math.NA · cs.NA· physics.comp-ph

High-order kernel regularization of singular and hypersingular Helmholtz boundary integral operators

Luiz M. Faria , Carlos Perez-Arancibia , Svetlana Tlupova This is my paper

Pith reviewed 2026-05-10 10:31 UTC · model grok-4.3

classification 🧮 math.NA cs.NAphysics.comp-ph
keywords boundary integral operatorsHelmholtz equationkernel regularizationhypersingular operatorerror functionsnumerical quadraturescattering problemshigh-order methods
0
0 comments X

The pith

Error-function smoothing with polynomial moment corrections regularizes all four singular and hypersingular Helmholtz boundary integral operators at high order for standard quadrature in three dimensions.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper extends a kernel regularization approach to the complete set of on-surface boundary integral operators arising in the three-dimensional Helmholtz Calderon calculus. Each singular kernel receives a smooth replacement formed from error functions plus a polynomial correction whose coefficients are fixed by matching the first several moments of the original kernel. A unified error analysis couples the regularization width to the mesh size so that the smoothing error and the quadrature error balance, producing explicit convergence rates set by the regularization order and the quadrature rule's degree of exactness. The resulting method reduces every operator evaluation to the numerical integration of a smooth integrand with ordinary surface quadrature, without local solves, singularity-specific rules, or precomputations. Numerical tests confirm the predicted rates and apply the operators to sound-soft and sound-hard scattering problems on smooth obstacles.

Core claim

Each of the four kernels (single-layer, double-layer, adjoint double-layer, hypersingular) is replaced by a smooth modification constructed from the error function and a low-degree polynomial correction. The polynomial coefficients are chosen so that the first several moments of the modified kernel match those of the original singular kernel. When the regularization parameter scales with mesh size, the combined regularization-plus-quadrature error converges at a rate determined by the minimum of the moment-matching order and the quadrature exactness degree. The construction is uniform across operators and applies to both Helmholtz and Laplace kernels.

What carries the argument

The regularizing function: an error-function-based smoothing of the singular part together with a polynomial correction whose coefficients are fixed by moment conditions that enforce agreement with the original kernel through low-order moments.

If this is right

  • The same regularization construction and error bound apply uniformly to the single-layer, double-layer, adjoint double-layer, and hypersingular operators.
  • Overall convergence rates are explicit once the regularization parameter is scaled with mesh size and are governed by the minimum of regularization order and quadrature degree of exactness.
  • Numerical evaluation reduces entirely to standard quadrature of smooth surface integrals, with no need for element-local solves or specialized singularity rules.
  • H-matrix acceleration can be applied in black-box fashion to the resulting dense matrices.
  • The method produces accurate solutions for sound-soft and sound-hard acoustic scattering on smooth obstacles while verifying the predicted rates.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The reduction to ordinary quadrature could be paired with existing H-matrix libraries to accelerate large-scale scattering simulations without rewriting fast multipole or other kernel-specific accelerators.
  • If analogous moment conditions can be derived for other singular kernels, the framework may extend directly to electromagnetic or elastic boundary integral equations.
  • Local adaptation of the regularization width to local mesh size could support adaptive surface discretizations for problems with localized features or corners.

Load-bearing premise

The moment conditions on the polynomial corrections, together with the error-function smoothing, produce the stated high-order accuracy for the hypersingular operator without additional restrictions on surface curvature or solution regularity beyond those stated in the error analysis.

What would settle it

Apply the regularized hypersingular operator to a known smooth solution on a sequence of uniformly refined meshes of a smooth closed surface and check whether the observed error decays at the rate predicted by the chosen regularization order and quadrature degree.

Figures

Figures reproduced from arXiv: 2604.14797 by Carlos Perez-Arancibia, Luiz M. Faria, Svetlana Tlupova.

Figure 1
Figure 1. Figure 1: Illustration of the decomposition of the surface Γ into annular regions Γ [PITH_FULL_IMAGE:figures/full_fig_p019_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: Convergence order o ⋆ as given in (72) for each integral operator considered in this work—S, K, K ⊤, H, and W—corresponding to (p, s) = (0, 0), (1, 0), (1, 0), (2, 1), and (2, 0), respectively, for various quadrature degrees of exactness q and regularization orders m. The following theorem summarizes the theoretical results of this section: Theorem 3.1 (Quadrature and discretization error). Let Γ be a smoo… view at source ↗
Figure 3
Figure 3. Figure 3: Normalized error e (B) nrm for B ∈ {S,K,K ⊤,T} defined in (82) as a function of the regularization parameter δ for the single layer (S), double layer (K), adjoint double layer (K ⊤), and hypersingular (T) operators, for regularization orders m ∈ {3, 5, 7, 9} (left to right, top to bottom). The dashed black line shows the theoretical O(δ m) convergence rate. Computations are performed on S 2 with wavenumber… view at source ↗
Figure 4
Figure 4. Figure 4: Convergence of the relative error e (B) for B ∈ {S,K,K ⊤,T} as a function of the mesh size h, with fixed regularization parameter δ = 0.05, wavenumber k = π, and regularization order m = 5, for quadrature orders q ∈ {2, 4, 5}. Dashed lines indicate the expected O(h q+1) convergence rates. haviour as h → 0. This is because the combined error bound satisfies e (B) (δ, h, k) ≲  c0 + c1|I (B) p,m(κ)|  h o ⋆ … view at source ↗
Figure 5
Figure 5. Figure 5: Overall relative error e (B) for B ∈ {S,K,K ⊤,T} defined in (81) as a function of the mesh size h for the single layer (S), double layer (K), adjoint double layer (K ⊤), and hypersingular (T) operators, for regularization orders m ∈ {3, 5, 7, 9} (left to right, top to bottom) and for quadrature rules of degree of exactness q ∈ {2, 4, 5}. The dashed black line shows the theoretical O(h o ⋆ ) convergence rat… view at source ↗
Figure 6
Figure 6. Figure 6: Normalized relative error ee (B) nrm for B ∈ {S,K,K ⊤,T} defined in (84) as a function of the mesh size h for the single layer (S), double layer (K), adjoint double layer (K ⊤), and hypersingular (T) operators, for regularization orders m ∈ {3, 5, 7, 9} (left to right, top to bottom) and for quadrature rules of degree of exactness q ∈ {2, 4, 5}. The dashed black line shows the theoretical O(h o ⋆ ) converg… view at source ↗
Figure 7
Figure 7. Figure 7: displays the far-field errors (89) for CFIE-D (87a) and CFIE-N (87b), discretized via the Nystr¨om method, as a function of mesh size h for the torus (major radius 1, minor radius 1/2) and the bean￾shaped surface from [7], at wavenumber k = π. The incident point source is located at x0 = (1, 1, 0) for the torus and x0 = (0.1, 0.1, 0.1) for the bean. Three parameter combinations (m, q) ∈ {(3, 2),(5, 4),(7, … view at source ↗
Figure 8
Figure 8. Figure 8: Real part of the total fields u tot = u + u inc and v tot = v + u inc for plane-wave scattering u inc(x) = eikx·d at wavenumber k = 5π, computed via the Nystr¨om method with regularization order m = 5 and quadrature degree q = 5. Top row: torus (major radius 1, minor radius 1/2); bottom row: bean-shaped surface [7]. Left column: sound-soft (Dirichlet) problem (85a); right column: sound-hard (Neumann) probl… view at source ↗
read the original abstract

This paper extends and analyzes the high-order kernel regularization framework of Beale & Tlupova (arXiv:2510.13639) to all four on-surface boundary integral operators of the Helmholtz Calderon calculus in three dimensions: the single-layer, double-layer, adjoint double-layer, and hypersingular operators. To the best of our knowledge, this work provides the first high-order kernel regularization of the hypersingular operator for both the Helmholtz and Laplace equations in three dimensions. The regularization replaces each singular kernel with a smooth modification constructed from error functions together with a polynomial correction whose coefficients are determined through moment conditions. Alongside the derivation of the regularizing functions, the paper provides a unified error analysis of the combined regularization and quadrature discretization procedure. By coupling the regularization parameter to the mesh size, the two error contributions can be balanced, leading to explicit overall convergence rates that depend jointly on the order of the regularization and the degree of exactness of the surface quadrature rule. A key practical feature of the method is its implementation simplicity. Once the regularizing functions are determined, the numerical task reduces entirely to the evaluation of smooth surface integrals using standard quadrature, without the need for element-local solves, singularity-specific precomputations, or specialized quadrature rules. Although the modified kernel is generally incompatible with kernel-specific fast methods, this limitation is addressed through H-matrix acceleration, applicable in a black-box manner. Numerical examples -- including verification of the predicted convergence rates and solution of sound-soft and sound-hard scattering problems by smooth obstacles -- demonstrate the accuracy and practicality of the proposed methodology.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

1 major / 2 minor

Summary. The manuscript extends the high-order kernel regularization framework of Beale & Tlupova to the single-layer, double-layer, adjoint double-layer, and hypersingular operators for the 3D Helmholtz Calderón calculus. Regularizers are constructed from error-function smoothing plus polynomial corrections fixed by moment conditions; a unified error analysis couples the regularization width to mesh size to balance regularization and quadrature errors, yielding explicit convergence rates depending on regularization order and quadrature exactness. The approach reduces to standard smooth quadrature plus H-matrix acceleration, with numerical verification of rates and application to sound-soft and sound-hard scattering.

Significance. If the error analysis for the hypersingular operator holds under the stated assumptions, the work provides a practical, unified regularization method for all four operators that avoids element-local solves or specialized quadrature rules. Strengths include the deterministic construction of the regularizers, the explicit balancing argument, and the numerical confirmation of predicted rates for both Laplace and Helmholtz cases. The implementation simplicity and black-box H-matrix treatment are clear advantages for large-scale problems.

major comments (1)
  1. The unified error analysis asserts that the moment conditions derived from the error-function smoothing suffice to cancel singular terms up to the target order for the hypersingular operator under the same surface and solution regularity used for the single- and double-layer cases. However, the hypersingular kernel involves two normal derivatives, which introduce geometry-dependent curvature terms (second fundamental form) not necessarily present in the flat-space moment-matching derivation. It is unclear whether these residuals are controlled by the stated assumptions or whether they produce an O(h^k) error that prevents the claimed balancing when the regularization parameter is tied to mesh size. A explicit expansion of the remainder for a curved patch or a numerical test with known curvature would resolve this.
minor comments (2)
  1. The abstract states that the method is 'incompatible with kernel-specific fast methods' but is addressed via H-matrix acceleration; the numerical section should specify the H-matrix parameters (admissibility, maximum rank, tolerance) used in the scattering examples to allow reproduction.
  2. Table or figure captions for the convergence studies should explicitly list the regularization width scaling (e.g., ε = C h^α) and the quadrature degree of exactness for each operator and mesh family.

Simulated Author's Rebuttal

1 responses · 0 unresolved

We thank the referee for the careful review and constructive feedback on our manuscript extending the high-order kernel regularization to the Helmholtz Calderón operators. We address the major comment below and plan to revise the manuscript to clarify the analysis for the hypersingular operator.

read point-by-point responses

  1. Referee: The unified error analysis asserts that the moment conditions derived from the error-function smoothing suffice to cancel singular terms up to the target order for the hypersingular operator under the same surface and solution regularity used for the single- and double-layer cases. However, the hypersingular kernel involves two normal derivatives, which introduce geometry-dependent curvature terms (second fundamental form) not necessarily present in the flat-space moment-matching derivation. It is unclear whether these residuals are controlled by the stated assumptions or whether they produce an O(h^k) error that prevents the claimed balancing when the regularization parameter is tied to mesh size. A explicit expansion of the remainder for a curved patch or a numerical test with known curvature would resolve this.

    Authors: We thank the referee for pointing out this subtlety in the hypersingular case. The derivation in the manuscript proceeds by expressing the kernels in local tangential coordinates on the surface, where the singular parts (including those from the two normal derivatives) are expanded using the fundamental solution and its derivatives. The moment conditions are set to match and cancel these singular contributions up to the desired order, with the error function providing the smoothing. Curvature effects enter through the metric and normal variations, but under the assumed surface regularity (sufficiently smooth for the Taylor expansions to hold), these appear as regular terms that are absorbed into the O(ε^m) regularization error, where m is the regularization order. Consequently, the balancing ε ~ h^{q/m} with quadrature error of order h^q remains valid without additional O(h^k) terms. To strengthen the presentation, we will add an explicit remainder expansion for a curved surface patch in the revised manuscript, along with numerical verification on a sphere demonstrating the convergence rates for the hypersingular operator. revision: yes

Circularity Check

0 steps flagged

Minor self-citation to prior regularization framework; derivation and error analysis otherwise independent

full rationale

The paper extends the kernel regularization approach from Beale & Tlupova (arXiv:2510.13639) with overlapping authorship (Tlupova), but the central contributions—the explicit construction of regularizing functions for the hypersingular operator via error functions and moment conditions, plus the unified error analysis balancing regularization parameter against quadrature degree—are developed directly in this work using deterministic matching and standard quadrature theory. No prediction reduces to a fitted input by construction, no ansatz is smuggled, and no uniqueness theorem from self-citation is invoked to force the result. The chain remains self-contained.

Axiom & Free-Parameter Ledger

0 free parameters · 1 axioms · 0 invented entities

The regularizing functions are built from error functions and polynomials whose coefficients are fixed by moment conditions; no parameters are fitted to data. The error analysis invokes classical assumptions on surface smoothness and quadrature exactness.

axioms (1)
  • domain assumption The surface is sufficiently smooth and the solution possesses the regularity needed for the quadrature error estimates to hold at the stated orders.
    Invoked to derive the balanced convergence rates when the regularization parameter is coupled to mesh size.

pith-pipeline@v0.9.0 · 5593 in / 1453 out tokens · 27581 ms · 2026-05-10T10:31:06.031519+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

29 extracted references · 29 canonical work pages

  1. [1]

    Abramowitz and I

    M. Abramowitz and I. A. Stegun.Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables, volume 55. Courier Corporation, 1965

    work page 1965

  2. [2]

    T. G. Anderson, M. Bonnet, L. M. Faria, and C. P´ erez-Arancibia. Fast, high-order numerical evalu- ation of volume potentials via polynomial density interpolation.Journal of Computational Physics, 511:113091, 2024

    work page 2024

  3. [3]

    J. T. Beale and S. Tlupova. Extrapolated regularization of nearly singular integrals on surfaces. Advances in Computational Mathematics, 50(4):61, 2024

    work page 2024

  4. [4]

    J. T. Beale and S. Tlupova. High order regularization of nearly singular surface integrals.arXiv preprint arXiv:2510.13639, 2025

    work page arXiv 2025

  5. [5]

    J. T. Beale, W. Ying, and J. R. Wilson. A simple method for computing singular or nearly singular integrals on closed surfaces.Communications in Computational Physics, 20(3):733–753, 2016

    work page 2016

  6. [6]

    O. P. Bruno and E. Garza. A Chebyshev-based rectangular-polar integral solver for scattering by geometries described by non-overlapping patches.Journal of Computational Physics, 421:109740, 2020

    work page 2020

  7. [7]

    O. P. Bruno and L. A. Kunyansky. A fast, high-order algorithm for the solution of surface scattering problems: basic implementation, tests, and applications.Journal of Computational Physics, 2001

    work page 2001

  8. [8]

    L. M. Faria, C. P´ erez-Arancibia, and M. Bonnet. General-purpose kernel regularization of boundary integral equations via density interpolation.Computer Methods in Applied Mechanics and Engineering, 378:113703, 2021

    work page 2021

  9. [9]

    Ganesh and I

    M. Ganesh and I. G. Graham. A high-order algorithm for obstacle scattering in three dimensions. Journal of Computational Physics, 198(1):211–242, 2004

    work page 2004

  10. [10]

    Ganesh and S

    M. Ganesh and S. C. Hawkins. A spectrally accurate algorithm for electromagnetic scattering in three dimensions.Numerical Algorithms, 43(1):25–60, 2006. 39

    work page 2006

  11. [11]

    Geuzaine and J.-F

    C. Geuzaine and J.-F. Remacle. Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities.International journal for numerical methods in engineering, 79(11):1309– 1331, 2009

    work page 2009

  12. [12]

    Greengard and V

    L. Greengard and V. Rokhlin. A new version of the fast multipole method for the Laplace equation in three dimensions.Acta Numerica, 6:229–269, 1997

    work page 1997

  13. [13]

    N. A. Gumerov and R. Duraiswami.Fast multipole methods for the Helmholtz equation in three dimensions. Elsevier, 2005

    work page 2005

  14. [14]

    Hackbusch.Hierarchical Matrices: Algorithms and Analysis

    W. Hackbusch.Hierarchical Matrices: Algorithms and Analysis. Springer, 2015

    work page 2015

  15. [15]

    G. C. Hsiao and W. L. Wendland.Boundary Integral Equations. Springer, 2021

    work page 2021

  16. [16]

    J. Hu, E. Garza, and C. Sideris. A chebyshev-based high-order-accurate integral equation solver for maxwell’s equations.IEEE Transactions on Antennas and Propagation, 69(9):5790–5800, 2021

    work page 2021

  17. [17]

    HMatrices.jl, 2025

    IntegralEquations contributors. HMatrices.jl, 2025

    work page 2025

  18. [18]

    Kl¨ ockner, A

    A. Kl¨ ockner, A. Barnett, L. Greengard, and M. O’Neil. Quadrature by expansion: A new method for the evaluation of layer potentials.Journal of Computational Physics, 252:332–349, 2013

    work page 2013

  19. [19]

    S. G. Krantz and H. R. Parks.The Implicit Function Theorem: History, Theory, and Applications. Springer Science & Business Media, 2002

    work page 2002

  20. [20]

    R. Kress. On the numerical solution of a hypersingular integral equation in scattering theory.Journal of Computational and Applied Mathematics, 61(3):345–360, Aug. 1995

    work page 1995

  21. [21]

    Maltez-Faria, C

    L. Maltez-Faria, C. P´ erez-Arancibia, and C. Turc. Combined field-only boundary integral equations for PEC electromagnetic scattering problem in spherical geometries.SIAM Journal on Applied Math- ematics, 84(1):19–38, 2024

    work page 2024

  22. [22]

    P´ erez-Arancibia, L

    C. P´ erez-Arancibia, L. M. Faria, and C. Turc. Harmonic density interpolation methods for high-order evaluation of Laplace layer potentials in 2D and 3D.Journal of Computational Physics, 376:411–434, 2019

    work page 2019

  23. [23]

    P´ erez-Arancibia, C

    C. P´ erez-Arancibia, C. Turc, and L. Faria. Planewave density interpolation methods for 3D Helmholtz boundary integral equations.SIAM Journal on Scientific Computing, 41(4):A2088–A2116, 2019

    work page 2019

  24. [24]

    Saad and M

    Y. Saad and M. H. Schultz. GMRES: A generalized minimal residual algorithm for solving nonsym- metric linear systems.SIAM Journal on Scientific and Statistical Computing, 7(3):856–869, 1986

    work page 1986

  25. [25]

    S. A. Sauter and C. Schwab.Boundary Element Methods. Springer, 2010

    work page 2010

  26. [26]

    Vioreanu and V

    B. Vioreanu and V. Rokhlin. Spectra of multiplication operators as a numerical tool.SIAM Journal on Scientific Computing, 36(1):A267–A288, 2014

    work page 2014

  27. [27]

    H. Whitney. Differentiable even functions. InHassler Whitney Collected Papers, pages 309–310. Springer, 1943

    work page 1943

  28. [28]

    Wu and P.-G

    B. Wu and P.-G. Martinsson. Corrected trapezoidal rules for boundary integral equations in three dimensions.Numerische Mathematik, 149:1025–1071, 2021. 40

    work page 2021

  29. [29]

    Wu and P.-G

    B. Wu and P.-G. Martinsson. A unified trapezoidal quadrature method for singular and hypersingular boundary integral operators on curved surfaces.SIAM Journal on Numerical Analysis, 61(5):2182– 2208, 2023. 41

    work page 2023