Nearly singular, close evaluation, singularity swap quadrature, error control, integral equation, body of revolution

1 Introduction↩︎

We consider the numerical evaluation of layer potentials of the form \[u(\xx) = \int_S\frac{k(\xx,\yy)~\sigma(\yy)}{\left\|\yy-\xx\right\|^{2p}}\dS(\yy),\quad p=1/2,~3/2,~5/2, \label{eq:generic95layer95potential}\tag{1}\] where \(\dS\) is the surface element measure of a smooth axisymmetric surface \(S\subset\mathbb{R}^3\), \(\sigma\) is a smooth scalar- or vector-valued function defined on \(S\), and \(k(\xx,\yy)\) is a smooth function originating from a fundamental solution of an elliptic partial differential equation (PDE). The evaluation (target) point \(\xx\in\mathbb{R}^3\) may lie far from or close to the surface, but not on it.

Layer potentials of the form 1 arise in boundary integral formulations of elliptic PDEs and can be viewed as a convolution of the PDE’s Green’s function with the unknown “density” \(\sigma\) over the boundary. Enforcing boundary conditions leads to an integral equation for \(\sigma\), and once this equation is solved, the solution anywhere in the domain is obtained by evaluating the associated layer potential.

Such settings appear in many applications. At the microscale, for example, Stokes flow governs the motion of fluid–particle systems. These arise in the rheology of fiber and polymer suspensions [1], suspensions of spheres [2], the design of new materials [3], [4], and the self-assembly of biological or synthetic particles [5]. Axisymmetric geometries such as spheroids, rods, rings, and general bodies of revolution occur naturally in these contexts. Axisymmetric geometries are also important in acoustics and wave propagation governed by the Helmholtz equation, particularly in the design of absorbing surfaces and high-power loudspeakers. Boundary integral solvers have proven effective in both Stokes [6]–[10] and Helmholtz [11]–[13] settings. Nevertheless, accurate and error-controlled evaluation of layer potentials close to the surface remains a challenging and active area of research.

The difficulty is well known. As the target approaches the surface, the Green’s function develops a sharp peak and the integrand varies rapidly. A straightforward remedy is to upsample the density, that is, interpolate \(\sigma\) to a finer surface grid and apply a fixed standard quadrature based on the surface discretization. However, the refinement needed to maintain accuracy grows dramatically as the target-surface distance shrinks, making this approach inefficient. This creates a need for special quadrature methods designed specifically for nearly singular integrals.

A large body of work addresses this need. In two dimensions, where boundary integrals reduce to line integrals, the problem is largely considered “solved”, but in three dimensions it remains an active research area. Quadrature by expansion (QBX) [8], [14], [15] was first introduced in two dimensions and later extended to three dimensions. Since QBX integrates naturally with spherical-harmonic-based fast multipole methods, an integration of the two was envisioned. Yet achieving high efficiency proved difficult, and substantial effort was required before robust software for Laplace and Helmholtz problems became available [16]. Important advances include target-specific expansions [17], which significantly reduce near-field cost and enable faster evaluations [18]. Another promising direction, introduced in [19] and improved in [20], transforms a surface integral into a collection of line integrals via Stokes’ theorem. These nearly singular line integrals are then evaluated using the singularity swap quadrature (SSQ) method, first proposed in [21] and subsequently extended in [22], [23], together with the recent stabilization in [24]. SSQ “swaps” the original nearly singular factor for a simpler one with matching singularities in the complex parameter plane. The remaining smooth factor is expanded in a suitable basis, and each basis function is then integrated analytically against the simplified near-singular term. These analytic integrals are typically computed via recurrence relations, whose stability may require additional care. Although the current Stokes’-theorem formulation is limited to Laplace layer potentials, the approach shows considerable potential.

Like QBX, the hedgehog or line-extrapolation method [25], [26] exploits the fact that the layer-potential field remains smooth even when the integrand is sharply peaked. The potential is evaluated by standard quadrature at off-surface points farther away along a line from the surface, and is then extrapolated toward it, optionally including the surface value itself when available. Although appealing in its simplicity, the achievable accuracy is strongly dependent on the extrapolation distance, and the method’s non-trivial parameter selection makes it challenging to use optimally.

Regularization methods replace the singular kernel by a smooth approximation, and recent work offers high-order approximations in the regularization parameter [27]. One remaining difficulty lies in balancing the regularization error against the quadrature error, which increases as the regularization error decreases. For spherical geometries, vectorial spherical-harmonic-based quadrature has been explored [28] and was recently extended to oblate and prolate spheroids [29].

Axisymmetric geometries have also motivated significant work on specialized integral equation solvers and quadrature schemes. For Stokes flows, one strategy for making QBX more efficient is to precompute target-specific quadrature weights [7], [8]. This hides the QBX cost for on-surface evaluations. For off-surface, close evaluations, however, the expansion still must be recomputed because the target locations are not known in advance. This remains a computational bottleneck, and local panel-based forms of QBX could potentially mitigate this cost. For Helmholtz problems, the approach of [11] applies a discrete Fourier transform in the azimuthal direction, reducing the problem to Fourier integrals along the generating curve that are evaluated by analytic recursions. Helsing and Karlsson [13] improved this by introducing analytic product integration in the axial direction based on [30], and by addressing instabilities in the Fourier recurrences. These methods rely on the axisymmetric nature of the entire problem. While the quadrature ideas might extend to problems where only the geometry is axisymmetric, this was not explored.

Because special quadrature schemes are substantially more expensive than standard quadratures, it is desirable in any boundary integral method to activate them only when strictly necessary. This in turn requires a reliable mechanism for determining, for each target point, whether the standard quadrature is sufficiently accurate. In this work, that role is played by the error predictors introduced in [31] and later refined for axisymmetric geometries in [32], both of which build on the asymptotic error analysis developed over a sequence of earlier studies (reviewed in [31]).

These predictors approximate the quadrature error by locating the relevant complex singularities of the integrand and inserting their positions into closed-form asymptotic formulas derived from the complex-variable framework of [33]–[35]. Although they are not rigorous upper bounds, they are designed to be computationally inexpensive and practically effective at flagging when the standard quadrature is inadequate. In practice, they have proven remarkably accurate, making them well suited as a decision mechanism for when to invoke special quadrature.

1.1 Contributions and outline↩︎

The main contribution of this paper is an adaptive semi-analytic quadrature method for evaluating nearly singular layer potentials of the form 1 on smooth axisymmetric surfaces. The method combines singularity swap quadrature (SSQ) in the azimuthal direction, and when necessary, an analogous SSQ-based treatment in the polar direction, combined with adaptive refinement. Assuming only that the underlying surface discretization resolves the geometry and layer density, the method automatically selects all internal parameters. The only required user input is the desired accuracy tolerance, although the order of the Gauss–Legendre rule used in the polar direction may be specified for efficiency. Because the resulting scheme uses SSQ in both surface-parameter directions, we refer to it as singularity swap surface quadrature (S3Q). Figure 1 illustrates the basic workflow.

The contributions of this paper, each of independent interest, are as follows:

• SSQ for near-singular line integrals on closed curves in three dimensions. We extend the SSQ method [21]–[23] to near-singular line integrals over smooth closed curves embedded in three dimensions, discretized by the trapezoidal rule. This forms the one-dimensional building block used later in the azimuthal direction of the surface algorithm.

• Quadrature and interpolation error predictors. We derive error predictors for adaptive refinement in the polar direction, including a new interpolation predictor obtained by combining the Hermite interpolation formula [36] with the complex-variable framework of [31]. This framework, previously used only for quadrature error prediction, is shown to yield accurate and practical interpolation error predictors for functions with near-singular branch-point behavior. These predictors are used to control the adaptive refinement in the polar direction.

• The S3Q algorithm for axisymmetric surface integrals. By combining SSQ in both the azimuthal and polar directions, using the azimuthal SSQ as a building block and adaptive, predictor-based refinement in the polar direction, we obtain a fully automated and error-controlled scheme for evaluating single- and double-layer potentials for the Laplace, Helmholtz, and Stokes equations for target points close to smooth axisymmetric surfaces.

The outline of this paper is as follows. Section 2 introduces notation and basic formulas. Section 3 develops SSQ for near-singular line integrals on closed curves in three dimensions, providing the one-dimensional building block used in the azimuthal direction of S3Q. Section 4 lifts this machinery to axisymmetric surface integrals and combines it with adaptive refinement in the polar direction to obtain the full S3Q algorithm. Section 5 develops the practical quadrature and interpolation error predictors that drive this adaptive refinement. Section 6 discusses the computation of complex roots of distance functions used in the predictors. Section 7 analyzes complexity, and Section 8 presents numerical results. We conclude in Section 9. Appendices 10–12 contain auxiliary proofs, derivations of analytical root formulas, and details on stabilizing the recurrence relations used in the SSQ weights.

2 Preliminaries↩︎

2.1 Geometry and parameterization↩︎

Let \(S\) be a smooth surface in \(\mathbb{R}^3\) parameterized by \(\ggamma:E\rightarrow\mathbb{R}^3\), where \(E=\{E_1\times E_2\}\subset\mathbb{R}^2\). The layer potential 1 can then be written in parameter space as \[u(\xx) = \iint_E \frac{k(\xx,\ggamma(\theta,\varphi))~\sigma(\ggamma(\theta,\varphi))}{\left\|\ggamma(\theta,\varphi)-\xx\right\|^{2p}}~\left\|\frac{\partial\ggamma}{\partial \theta} \times \frac{\partial\ggamma}{\partial\varphi}\right\|\dphi\dtheta = \iint_E\frac{f(\theta,\varphi)}{\left\|\ggamma(\theta,\varphi)-\xx\right\|^{2p}}\dphi\dtheta, \label{eq:parameterized95layer95potential}\tag{2}\] where the function \(f(\theta,\varphi)\) contains all the smooth components, and implicitly depends on the target point \(\xx\).

Throughout this paper we consider axisymmetric surfaces. We take \(E=\{(\theta,\varphi)\in[0,\pi]\times[0,2\pi)\}\) and use the axisymmetric parametrization \[\ggamma(\theta,\varphi) = \left(a(\theta)\sin(\theta)\cos(\varphi),~a(\theta)\sin(\theta)\sin(\varphi),~b(\theta)\cos(\theta)\right), \label{eq:gamma95axi}\tag{3}\] where \(a(\theta)>0\) and \(b(\theta)>0\) are smooth analytic functions.

To simplify later notation, it is convenient to introduce a linear mapping between \(t\in[-1,1]\) and a subinterval \([\tha,\thb]\subseteq[0,\pi]\) \[\theta(t,\tha,\thb) = \frac{\tha+\thb}{2} + \thsc t,\quad \thsc = \frac{\thb-\tha}{2}. \label{eq:t2theta95map}\tag{4}\] with inverse \[t(\theta,\tha,\thb) = \frac{1}{\thsc}\left(\theta-\frac{\tha+\thb}{2}\right). \label{eq:theta2t95map}\tag{5}\] When \([\tha,\thb]=[0,\pi]\), we write simply \(\theta(t)\) and \(t(\theta)\).

Next, we define the squared-distance function for the surface \(\ggamma(\theta,\varphi)=(\gamma_1(\theta,\varphi),\gamma_2(\theta,\varphi)),\gamma_3(\theta,\varphi))\) and a target point \(\xx=(x,y,z)\), \[R^2(\theta,\varphi,\xx) \mathrel{\vcenter{:}}= \|\ggamma(\theta,\varphi)-\xx\|^2 = \left(\gamma_1(\theta,\varphi)-x\right)^2 + \left(\gamma_2(\theta,\varphi)-y\right)^2 + \left(\gamma_3(\theta,\varphi)-z\right)^2. \label{eq:R2}\tag{6}\] Using 4 , we write the layer potential compactly as \[\I[\Theta_p](\xx) = \I_t\I_\varphi[\Theta_p](\xx) = \iint_E \Theta_p(\theta(t,\tha,\thb),\varphi,\xx)\thsc\dphi\dt, \quad \Theta_p(\theta,\varphi,\xx) = \frac{f(\theta,\varphi)}{\left(R^2(\theta,\varphi,\xx)\right)^{p}}. \label{eq:Itheta}\tag{7}\] Here, the notation \(\I_t\) and \(\I_\varphi\) indicate integration in the \(t\)- and \(\varphi\)-direction, respectively. We will later use similar subscripts to distinguish in which direction an operator is applied.

2.2 Discretization↩︎

A natural approximation of 7 is the tensor-product rule formed from an \(\nt\)-point Gauss–Legendre quadrature rule on \(t\in[-1,1]\) nodes \(\{t_k\}_{k=1}^\nt\) and weights \(\{\wt\}_{k=1}^\nt\), and an \(\nphi\)-point trapezoidal rule on \(\varphi\in[0,2\pi)\) with nodes \(\{\varphi_\ell\}_{\ell=1}^\nphi\) and weights \(\{\wphi\}_{\ell=1}^\nphi\). The standard quadrature becomes \[\Q_{\nt,\nphi}[\Theta_p](\xx)=\Q_{t,\nt}\Q_{\varphi,\nphi}[\Theta_p](\xx) = \sum_{k=1}^\nt\sum_{\ell=1}^\nphi \frac{f(\theta(t_k),\varphi_\ell)\wphi\wt\thsc}{\left(R^2(\theta(t_k),\varphi_\ell,\xx)\right)^{p}}. \label{eq:reg95quad}\tag{8}\] with corresponding error \[\left|\E_{\nt,\nphi}[\Theta_p](\xx)\right| = \left|\I[\Theta_p](\xx) - \Q_{\nt,\nphi}[\Theta_p](\xx)\right|. \label{eq:reg95quad95err}\tag{9}\]

2.3 Complex roots of↩︎

Although the integrand in 7 is smooth for all real \((\theta,\varphi)\) with target \(x\notin S\), it develops sharp peaks when the target lies close to the surface. For a fixed value of one parametrization variable the function becomes singular, i.e., \(R^2(\theta,\varphi,\xx)\) has roots at certain points in the complexified plane of the other variable.³ These singularities bound the region of analyticity. Several complex-conjugate root pairs may exist, but their influence decays exponentially with distance from the real interval; only the closest pair is relevant for our purposes.⁴

For fixed \(\theta\), we denote by \(\{\phiroot(\theta,\xx),\overline{\phiroot(\theta,\xx)}\}\) the pair closest to \(E_2\), i.e. the pair that governs the near-singular behavior. They satisfy, \[R^2\left(\theta,\phiroot(\theta,\xx),\xx\right) = R^2\left(\theta,\overline{\phiroot(\theta,\xx)},\xx\right) = 0.\]

Similarly, fixing \(\varphi\) and treating \(\theta\) (or equivalently \(t\)) as complex, the relevant roots are
\(\{\theta(t_0(\varphi,\xx)),\theta(\overline{t_0(\varphi,\xx)})\}\), where \(\troot(\varphi,\xx)\) solves \(R^2(\theta(t),\varphi,\xx)=0\) and lies closest to the interval \(E_1\).

These roots can sometimes be expressed analytically, but in general they must be found by a one-dimensional numerical root search. Their computation is discussed in Section 6. From this point forward, we write simply \(\phiroot\) and \(t_0\), leaving their dependencies implicit.

3 Singularity swap quadrature for near-singular line integrals on closed curves in↩︎

In the axisymmetric setting, evaluation of the layer potential 1 naturally leads to a family of one-dimensional integrals obtained by fixing the polar parameter \(\theta\) and integrating in the azimuthal direction. Each such slice is an integral over a circle in \(\mathbb{R}^3\), and its accurate evaluation is essential for resolving nearly singular behavior when the target point lies close to the surface.

Motivated by this, we develop singularity swap quadrature (SSQ) for the more general task of evaluating integrals of the form 1 over an arbitrary smooth analytic closed curve \(\Gamma\), following [21]. This extension is of independent interest, but in the present paper its main role is to provide the one-dimensional machinery used in the azimuthal direction of the surface quadrature method developed in Section 4.

Let \(\boldsymbol{\xi}:[0,2\pi)\rightarrow\mathbb{R}^3\), where \(\boldsymbol{\xi}(\phi)=(\xi_1(\phi),\xi_2(\phi),\xi_3(\phi))\), be an analytic, \(2\pi\)-periodic parametrization of \(\Gamma\). To avoid confusion with the azimuthal parameter \(\varphi\) used for axisymmetric surfaces, we reserve \(\phi\) for the curve parameter.⁵

Define the squared-distance function \[R(\phi)^2 \mathrel{\vcenter{:}}= \|\boldsymbol{\xi}(\phi)-\xx\|^2 = \left(\xi_1(\phi)-x\right)^2+\left(\xi_2(\phi)-y\right)^2+\left(\xi_3(\phi)-z\right)^2\] Let \(\phi_0\), \(\overline{\phi_0}\) be the complex-conjugate roots of \(R(\phi)^2\) closest to the interval \([0,2\pi)\), assuming they lie within the region of analyticity of \(\boldsymbol{\xi}\). In practice, \(\phi_0\) can be found efficiently through a one-dimensional numerical root search, see [22] for details.

The parametrized integral takes the form \[I_p = I_p(\xx) = \int_{0}^{2\pi} \frac{g(\phi)}{(R(\phi)^2)^p}\dpphi,\quad p=1/2,~3/2,~5/2, \label{eq:Ip}\tag{10}\] with \(g(\phi)\mathrel{\vcenter{:}}= k(\xx,\boldsymbol{\xi}(\phi))\sigma_\Gamma(\phi)\|\boldsymbol{\xi}'(\phi)\|\), where \(\sigma_\Gamma(\phi)\) denotes the restriction of the density \(\sigma\) to the curve \(\Gamma\).

The central idea of [21] is to transform the integral 10 over a curve in \(\mathbb{R}^3\) into an integral with a known singularity structure on a simple domain in \(\mathbb{C}\). Once the singular points—namely the complex roots of the squared-distance function—are identified via analytic continuation, the integrand can be rewritten in terms of a simpler factor that exhibits the same local behavior at those points. In the present setting the relevant singularities are \(\phi_0\), \(\overline{\phi_0}\), and the elementary function that vanishes at the same location is \(|e^{i\phi}-e^{i\phi_0}|^2\).

“Swapping” the singularity in 10 using this factor yields \[I_p = \int_{0}^{2\pi} \frac{G(\phi)}{|e^{i\phi}-e^{i\phi_0}|^{2p}}\dpphi,\] where the “SSQ numerator” \(G(\phi)\), defined by \[G(\phi)\mathrel{\vcenter{:}}= g(\phi)\frac{|e^{i\phi}-e^{i\phi_0}|^{2p}}{(R(\phi)^2)^p}, \label{eq:ssq95G}\tag{11}\] is a smooth and \(2\pi\)-periodic function on \([0,2\pi)\) (in fact analytic in a typically large open neighborhood of \([0,2\pi)\)).

Suppose \(G\) is sampled at \(n_\phi\) equispaced nodes, with \(n_\phi\) even. Expanding \(G\) in a truncated Fourier series and integrating termwise gives \[I_p \approx \sum_{k=-n_\phi/2}^{n_\phi/2-1} \widehat{G}_k(\phi_0)\underbrace{\int_{0}^{2\pi} \frac{e^{ik\phi}}{|e^{i\phi}-e^{i\phi_0}|^{2p}}\dpphi}_{\eqqcolon S_k^p(\phi_0)}, \label{eq:Ip95approx95ssq}\tag{12}\] where \(\widehat{G}_k\) are the Fourier coefficients of \(G\) obtained via a fast Fourier transform (FFT) at \(\mathcal{O}(n_\phi\log n_\phi)\) cost.

Thus \(I_p\) is approximated by an interpolatory quadrature rule on the unit circle, with weights \(S_k^p\) equal to the Fourier coefficients of \(|e^{i\phi}-e^{i\phi_0}|^{-2p}\). A key advantage is that all near-singularity is confined to the Fourier basis integrals \(S_k^p\), which are independent of the density \(\sigma\). These integrals can be evaluated analytically, at cost \(\mathcal{O}(n_\phi)\), using the recurrence relations stated in the following lemma. Its proof is given in Appendix 10.

In the lemma, and throughout the azimuthal SSQ construction, we take \(\varphi_0\) to be the root with positive imaginary part. Choosing the conjugate root gives an equivalent SSQ formulation, differing only by normalization factors in intermediate expressions.

With this one-dimensional SSQ machinery in place, we now return to the surface problem and use it as the core building block of S3Q in both the azimuthal and polar directions.

4 The S3Q method↩︎

In this section, we extend the one-dimensional construction of Section 3 from curves to smooth axisymmetric surfaces. The resulting method, which we call singularity swap surface quadrature (S3Q), gives a fully automated scheme for evaluating layer potentials of the form 1 to prescribed accuracy. We begin with an overview before giving the details.

Using notation from Section 2.1, the parameterized layer potential is \[\mathcal{I}(\xx) = \int_0^\pi\int_0^{2\pi} \frac{f(\theta,\varphi)}{(R^2(\theta,\varphi,\xx))^p}\dphi\dtheta. \label{eq:s3q95Ithetaphi}\tag{14}\] For each fixed \(\theta\), the inner integral in \(\varphi\) is taken over the circle \(\Gamma_\theta = \{\ggamma(\theta,\varphi) : \varphi\in[0,2\pi)\}\). This is precisely the type of one-dimensional integral treated by the SSQ machinery of Section 3.

We therefore use SSQ to evaluate the azimuthal integral. The result of this step is a function of \(\theta\) to be integrated. Because the computation of the inner integral relies on recurrence formulas, simple closed-form expressions for this function are not available. Although its dominant singular behavior is significantly weakened after integrating in \(\varphi\), it can still be sharply peaked when the target lies close to the surface. For this reason, the outer polar integral requires an adaptive strategy equipped with reliable quadrature and interpolation error predictors, together with an SSQ treatment of the residual singularity in the cases where regular Gauss–Legendre quadrature is insufficient.

The S3Q method can be summarized as follows:

0. Target classification. Given a prescribed tolerance, use the quadrature error predictor of [32] to determine whether regular tensor-product quadrature is expected to be insufficient. If so, proceed with the remaining steps.

1. Adaptive polar subdivision. Use the quadrature and interpolation error predictors of Section 5 to adaptively refine the \(\theta\)-interval.

2. Azimuthal SSQ. For each \(\theta\)-panel, apply SSQ to the closed curve \(\Gamma_\theta\) whenever the standard trapezoidal rule does not meet the required accuracy. This yields an approximation of \[\mathcal{I}_\varphi(\theta) = \int_0^{2\pi} \frac{f(\theta,\varphi)}{(R^2(\theta,\varphi,\xx))^p}\dphi. \label{eq:s3q95Iphi}\tag{15}\]

3. Polar SSQ: For each \(\theta\)-panel, apply SSQ to the generating curve whenever the standard Gauss–Legendre rule does not meet the required accuracy. The output is an approximation of \[\mathcal{I}(\xx) = \int_0^\pi\mathcal{I}_\varphi(\theta)\dtheta. \label{eq:s3q95Itheta}\tag{16}\]

This structure enables S3Q to resolve near-singular behavior efficiently, activating special quadrature only when needed and relying on inexpensive standard quadrature elsewhere. The following subsections describe each component in detail.

4.1 Singularity swap quadrature in the azimuthal variable↩︎

The goal of this subsection is to apply the SSQ machinery from Section 3 to evaluate the nearly singular azimuthal integral 15 for each fixed value of the polar value \(\theta\). Our main analytical result is stated in Theorem 1, which follows from two preparatory lemmas, and the SSQ formulas of Section 3.

Let \(\phiroot=\phiroot(\theta)\) denote the complex-valued root of \(\Rth\) with respect to \(\varphi\). For axisymmetric surfaces, \(\varphi_0\) is available in closed form via the formulas of [32], restated in Lemma 2. A key result in the axisymmetric setting is that the factor \(|e^{i\varphi}-e^{i\varphi_0}|^{2}/\Rth\) appearing in the “SSQ numerator”, similar to 11 , is independent of \(\varphi\). This property allows the azimuthal integral 15 to collapse to an expression involving only the Fourier coefficients of \(f(\theta,\cdot)\) and the Fourier basis integrals \(\mu_k^p\).

In what follows, \(\phiroot\) denotes the root in ?? with positive imaginary part, consistent with the convention of Section 3.

With these preparations we now present the main analytical reduction of the layer potential to a one-dimensional integral in \(\theta\). The theorem below shows that, after azimuthal integration and Fourier truncation, the remaining integrand always contains a square-root-type singularity, and for \(p>1/2\) a weakened near-singularity compared to the original.

In the next subsection, we analyze this structure in detail to determine how ?? is best integrated in the polar direction.

4.2 Examination of integral in the polar direction↩︎

Following the azimuthal SSQ reduction in Section 4.1, the layer potential is approximated by a one-dimensional integral in the polar parameter \(\theta\). The purpose of this section is to analyze the structure of the reduced integrand in ?? as a function of \(\theta\), and to identify the features that determine the accuracy of its numerical evaluation.

We write ?? in the form \[u(\xx) = \int_0^\pi \frac{1}{\sqrt{a(\theta)}} \Lsq(\theta,\xx)\Fcal(\theta,p,\xx) \left(\Rlambda(\theta,\xx)\right)^{-(p-1/2)} \dtheta, \label{eq:layer95potential95abc}\tag{20}\] where \[\Lsq(\theta,\xx) = \left(\lamb+\sqrt{\lamb^2-\rho^2\sin^2(\theta)}\right)^{-1/2}, \label{eq:Lsq}\tag{21}\] with \(\lambda\) defined in ?? . The function \(\Fcal\) is written with three arguments to emphasize its dependencies.

Among the factors in the integrand, the inverse \(\Rlambda\) term is the most difficult to resolve numerically, as it exhibits a nearly singular behavior for close target points. Its singularity order is, however, reduced compared to that of the original square-distance function \(R^2\), and for \(p=1/2\) this factor is identically equal to one. The function \(\Fcal(\theta,p,\xx)\), defined via recurrence relations, may have a logarithmic near singularity⁶. The factor \(\Lsq(\theta,\xx)\) has a nearly singular derivative and, while remaining bounded in magnitude, is nevertheless challenging to approximate accurately by global polynomial expansions. The remaining factor \(1/\sqrt{a(\theta)}\) is benign.

Figure 2 illustrates the typical behavior of the functions appearing in 20 for close evaluation of the layer potential 2 on a spheroidal surface. Each factor is represented by its global Chebyshev expansion, and the decay of the corresponding coefficient magnitudes is shown. The results indicate that \(\Lsq(\theta,\xx)\) and \(\Fcal(\theta,p,\xx)\) require comparable polynomial degrees to be resolved to a given accuracy, whereas for \(p>1/2\) the factor \((\Rlambda(\theta,\xx))^{-(p-1/2)}\) exhibits significantly slower coefficient decay and therefore requires much higher polynomial degree.

For \(p>1/2\), the function \((\Rlambda(\theta,\xx))^{-(p-1/2)}\) has a singularity at \(\thetarootlambda\in\mathbb{C}\), a complex root of \(\Rlambda\). This singularity can be regularized by multiplication with \(|\theta-\thetarootlambda|^{2p-1}\). As shown in Figure [fig:cheb95coeffs95p3], this regularization dramatically reduces the number of Chebyshev modes required to resolve the factor, confirming that the slow decay observed in the unregularized factor is directly associated with the complex root structure of \(\Rlambda\).

The functions \(\Lsq(\theta,\xx)\) and \(\Fcal(\theta,p,\xx)\) possess branch-point and logarithmic singularities, respectively, which limit the decay of global polynomial expansions, as observed in Figures [fig:cheb95coeffs95p1] and [fig:cheb95coeffs95p3]. These singularities are integrable and localized near the same complex locations that govern the behavior of \(\Rlambda\). While they affect the global approximation rate of the integrand, they do not introduce stronger singular behavior than that already present in the \(\Rlambda\) factor.

For this reason, the numerical treatment focuses on \(\Rlambda\). We employ an adaptive subdivision of the polar integration interval, and for \(p>1/2\), an SSQ-based approach is used to handle the \(\Rlambda\) factor, while for \(p=1/2\) standard Gauss–Legendre quadrature is applied on each subinterval. Details of the adaptive discretization are given in Section 4.5. We now proceed to describe the special quadrature construction for \(p>1/2\).

4.3 Singularity swap quadrature in the polar variable↩︎

We now describe the application of singularity swap quadrature in the polar variable for the evaluation of the reduced integral ?? when \(p>1/2\). The treatment is applied locally on a subinterval \([\tha,\thb]\subseteq[0,\pi]\) arising from the adaptive subdivision described later.

The procedure mirrors the steps used in the azimuthal direction in Section 4.1. It consists of identifying the complex roots of the relevant squared-distance function, swapping the singularity to a simpler function, expanding the remaining factors in a polynomial basis, and evaluating the resulting nearly singular basis integrals analytically.

Since \(\Rlambda(\theta,\xx)\) is real for real \(\theta\), its roots come in complex conjugate pairs. For spheroidal surfaces these roots are available analytically, while for general axisymmetric surfaces they are obtained via a one-dimensional complex root-finding procedure; see Section 6. When the generating curve is not excessively curved and the target point \(\xx\) is close to it, there is typically only a single pair of roots near the interval under consideration [21]. We denote this pair by \(\{\thetarootlambda(\xx),\overline{\thetarootlambda(\xx)}\}=\{\theta(\trootlambda(\xx),\tha,\thb),\theta(\overline{\trootlambda(\xx)},\tha,\thb)\}\) where \(\theta(t,\tha,\thb)\) is the affine map from \([-1,1]\) to \([\theta_a,\theta_b]\) defined in 4 . For brevity, the dependence of \(\trootlambda\) on \(\xx\) is suppressed in the remainder of this section.

Following af Klinteberg & Barnett [21], we define \[H(\theta) \mathrel{\vcenter{:}}= \frac{1}{\sqrt{a(\theta)}}\Lsq(\theta,\xx)\Fcal(\theta,p,\xx)\frac{|\theta-\thetarootlambda|^{2p-1}}{\left(\Rlambda(\theta,\xx)\right)^{p-1/2}} \label{eq:H}\tag{22}\] which allows the integral 20 over \([\tha,\thb]\) to be written as \[u(\xx) = \thsc^{2(1-p)}\int_{-1}^{1} \frac{H(\theta(t,\tha,\thb))}{|t-t_0^\lambda|^{2p-1}}\dt, \label{eq:Hintegral}\tag{23}\] where \(|\theta-\thetarootlambda|=\thsc|t-\trootlambda|\). The singularity swap has transformed the integral on the curve in the polar direction to a denominator that corresponds to that of a straight line.

The function \(H(\theta(t,\tha,\thb))\) is expanded in a real monomial basis on \([-1,1]\), \[H(\theta(t,\tha,\thb)) \approx \sum_{k=1}^\nt c_kt^{k-1}, \label{eq:H95poly}\tag{24}\] where the coefficients \(c_k\) are obtained by solving a Vandermonde system using the sampled values
\(H(\theta(t_i,\tha,\thb))\) with \(t_i\) being the Gauss–Legendre quadrature nodes. Substituting this expansion into the integral yields \[u(\xx) \approx \thsc^{2(1-p)}\sum_{k=1}^\nt c_k\nu_k^p(\trootlambda) \label{eq:layer95potential95ssq}\tag{25}\] where \[\nu_k^p(\trootlambda) = \int_{-1}^1\frac{t^{k-1}}{|t-\trootlambda|^{2p-1}}\dt,\quad k=1,\dots,\nt.\] These basis integrals are evaluated efficiently using the recurrence relations stated in Lemma 4, with proof and initial values given in Appendix 10.

The accuracy of this quadrature depends on how well \(H\) is resolved by polynomial interpolation on the current subinterval. While the singularity swap removes the dominant algebraic singularity associated with \((\Rlambda)^{-(p-1/2)}\), the analytic continuation of \(H\) remains limited by branch-point and logarithmic singularities inherited from \(\Lsq\) and \(\Fcal\). Consequently, a single global polynomial interpolant of \(H\) over \([0,\pi]\) is generally inefficient and may converge slowly when the complex singularities lie close to the real interval.

Rather than increasing the polynomial degree, we subdivide the interval \([0,\pi]\) in the polar variable. This subdivision is introduced for efficiency: on each subinterval, the nearest complex roots, when mapped to \([-1,1]\), are sufficiently far from the interval so that \(H\) can locally be well approximated by a low-degree polynomial. This enlarges the effective domain of analyticity on each panel and enables accurate polynomial interpolation using a fixed, modest number of Gauss–Legendre nodes (typically 16 or 32). The same subdivision strategy is employed for \(p=1/2\), where standard Gauss–Legendre quadrature is used instead of SSQ. Section 4.5 describes how this discretization refinement is realized adaptively.

4.4 Adjoint method for target-specific quadrature weights↩︎

The evaluation of the azimuthal integral in 15 and the polar integral in 16 using SSQ relies on the corresponding basis coefficients, namely Fourier coefficients in the azimuthal direction and monomial coefficients in the polar direction. For a fixed target point, these coefficients are determined from samples of a smooth numerator function that includes the kernel \(k(\xx,\yy)\) and the layer density \(\sigma(\yy)\). In many applications, however, the same target and kernel must be evaluated for several different densities, for instance in iterative boundary integral solvers. In such cases, it is advantageous to precompute target-specific quadrature weights that act directly on the sampled density. We refer to this as the adjoint method. Its practical construction in the monomial and Fourier settings is described in [21] and [22], respectively.

4.5 Local refinement scheme↩︎

To efficiently represent the function \(H\) in 22 using a monomial basis, we adopt a local refinement strategy in the polar direction based on panel subdivision. The refinement is driven by how well the function \(H\) can be resolved by fixed-order polynomial interpolation, recognizing that its limited regularity is inherited from the branch-point and logarithmic singularities of \(\Lsq\) and \(\Fcal\), and that these difficulties are localized in the polar variable, as illustrated in Figure 2.

Let \(P\in[-1,1]\) denote the base Gauss–Legendre panel on which both the geometry and layer density are well resolved, and let \(\xx\) be a target point with corresponding complex root \(\trootlambda\) of \(\Rlambda\). We subdivide \(P\) into \(\npan\) subpanels \([t_a^i,t_b^i]\), \(i=1,\dots,\npan\), each equipped with \(\nGL\)-point Gauss–Legendre nodes. The subdivision is organized as a binary tree: each panel may be bisected into two child panels of equal length, and refinement proceeds recursively.

The refinement is designed to satisfy \[\E(\nGL,t_a^i,t_b^i,p,\sigma,\xx) \leq \eps_i,\qquad i=1,\dots,\npan, \label{eq:dt95criteria95i}\tag{26}\] where \(\E(\nGL, t_\alpha, t_\beta, p, \sigma, \xx)\) denotes the error predictor for a general \(\nGL\)-point Gauss–Legendre panel on \([t_a, t_b]\), to be discussed in Section 5, and where the local tolerances satisfy \(\sum_{i=1}^{\npan} \eps_i = \eps\), with \(\eps\) being the prescribed global error tolerance. Panels for which the error predictor exceeds the assigned tolerance are bisected, and the process continues until the criterion is satisfied on all leaf panels.

This binary-tree structure has important practical advantages. Since all panels at a given refinement level are obtained by repeated bisection of the reference interval \([-1,1]\), interpolation matrices for transferring surface quantities from parent panels to child panels can be precomputed up to a fixed depth. For deeper levels, refinement is handled dynamically using precomputed local interpolation matrices for left and right subpanels. This approach avoids repeated construction of interpolation operators and leads to an efficient and robust implementation.

The output of the algorithm is a subdivision of \([-1,1]\) into \(\npan\) \(\nGL\)-point Gauss–Legendre panels that allow the reduced polar integral 20 to be evaluated with accuracy \(\eps\). The layer density is interpolated from the quadrature nodes on the base panel \(P\) to the nodes on each subpanel using barycentric Lagrange interpolation [39]. On each subpanel, the contribution is then evaluated by Gauss–Legendre quadrature when this is sufficient, and otherwise by the polar SSQ formula 23 . The subpanel contributions are finally summed.

To complete the description of the refinement strategy it remains to concretize the error predictor in 26 . This is the subject of the following section.

5 Error analysis and practical predictors↩︎

The goal of this section is to develop computable quadrature and interpolation error predictors that are sufficiently reliable to drive adaptive panel refinement in the polar direction, rather than to derive sharp or asymptotically exact bounds. Our starting point is a complex-analytic framework, reviewed in Section 5.1, in which the quadrature and interpolation errors admit exact contour-integral representations. Building on this framework and prior analyses [31], [32], we obtain practical predictor formulas by introducing asymptotic simplifications and retaining only the nearest complex singularities. These predictors, derived in Sections 5.2–5.3, form the basis of the adaptive panel subdivision strategy described in Section 4.5.

5.1 Quadrature and interpolation error formulas using complex analysis↩︎

Let \(g\) be a function defined on a real interval \(E\subset\mathbb{R}\), and consider the integral \[\I[g] = \int_E g(t)\dt, \label{eq:Ig}\tag{27}\] approximated by an \(n\)-point quadrature rule with quadrature nodes \(\{t_\ell\}_{\ell=1}^n\) and corresponding weights \(\{w_\ell\}_{\ell=1}^n\), \[\Q_n[g] = \sum_{\ell=1}^n g(t_\ell)w_\ell. \label{eq:Qng}\tag{28}\] The quadrature error is then given by \[\E_n[g] = \I[g] - \Q_n[g]. \label{eq:E}\tag{29}\] The rate at which \(\E_n[g]\) decay as \(n\) increases depends on the smoothness of \(g\) and the chosen quadrature rule.

Classical estimates of \(\E_n[g]\) for Gauss–Legendre quadrature typically involve high-order derivatives of \(g\). Such estimates are therefore only practical for very smooth integrands and tend to overestimate the error when \(g\) has singularities near \(E\). In nearly singular settings, they may even incorrectly suggest an increase in error as \(n\) grows, despite the actual error decreasing [40].

A more informative approach in this regime is based on complex analysis, as described by Donaldson & Elliot [33]. They express the error as a contour integral in the complex plane, \[\E_n[g] = \frac{1}{2\pi i} \int_C g(t)k_n(t)\dt, \label{eq:Ecomplex}\tag{30}\] where the contour \(C\) in Figure 4, containing the integration interval \(E\), is chosen so that the complex continuation of \(g\) is analytic on and inside \(C\). The error (or remainder) function \(k_n(t)\) depends on the quadrature rule and is detailed in Section 5.2.2 for the Gauss–Legendre quadrature rule.

A central technique in deriving asymptotic estimates of 30 , introduced in [34], involves deforming the contour \(C\) away from \(E\) while avoiding all singularities. Assuming the integrand tends to zero faster than \(|t|^{-1}\) as \(|t|\rightarrow\infty\), contributions from parts of \(C\) far from \(E\) vanish. Thus, \(C\) can be extended to infinity, leaving only the contribution from the singularities and their branch cuts. The particular choice of branch cuts in [34] is essential for obtaining accurate, simplified analytical asymptotic estimates. For functions with multiple pairs of complex conjugate singularities, such as 1 , all pairs contribute in principle, but typically only the pair closest to \(E\) has a significant contribution, and the others can safely be discarded [21], [31]. ⁷ The functions considered in this work exhibit branch-point singularities, meaning the main error contribution comes from neighborhoods of these branch points and their branch cuts, illustrated in Figure 4.

The same complex analysis framework applies also to polynomial interpolation. Given \(n\) distinct sample points \(\{t_k\}_{k=1}^n\subset E\), then there exists a unique polynomial \(p_{n-1}\) of degree at most \(n-1\) such that \(p(t_k)=g(t_k)\), \(k=1,\dots,n\). This polynomial can be represented in various ways, but due to its superior stability and efficiency, we construct it using the barycentric Lagrange interpolation formula [36], [39] \[p_{n-1}(t) = \elln(t)\sum_{k=1}^n\frac{w_kg(t_k)}{t-t_k},\] where the barycentric weights are \[w_k = \frac{1}{\elln^{\prime}(t_k)} = \left(\prod_{j=1,~j\neq k}^n(t_k-t_j)\right)^{-1},\quad k=1,\dots,n \label{eq:barycentric95weights}\tag{31}\] and the node polynomial is \[\elln(t) = \prod_{k=1}^n (t-t_k). \label{eq:elln}\tag{32}\] The polynomial interpolation error at \(\tau\in E\), \[\EI_n[g](\tau) = g(\tau)-p_{n-1}(\tau),\] admits, by the Hermite interpolation formula [36], the contour representation \[\EI_n[g](\tau) = \frac{1}{2\pi i} \int_C \frac{\elln(\tau)}{\elln(t)}\frac{g(t)}{(t-\tau)}\dt. \label{eq:Einterp95complex}\tag{33}\] This formula holds under the same assumptions as stated below 30 , namely that \(g\) is analytic on and inside \(C\), and that \(E\) is contained within \(C\). We will be concerned with the maximum absolute interpolation error over \(E\), \[\left|\EI_n[g]\right| = \max_{\tau\in E}~\left|\EI_n[g](\tau)\right|. \label{eq:Einterp95complex95max}\tag{34}\] Comparing 33 with 30 , we see that, the error function \(k_n(t)\) is replaced in the interpolation case by the factor \(\elln(\tau)/(\elln(t)(t-\tau))\), which depends on the node distribution and evaluation point \(\tau\). Crucially, both factors share the property that their magnitude decays rapidly away from the real interval \(E\).

5.2 Quadrature error predictors for adaptivity in polar direction↩︎

We now construct a quadrature error predictor for the function \(\Lsq\), defined in 21 , repeated here for convenience, \[\Lsq(\theta,\xx) = \left(\lamb+\sqrt{\lamb^2-\rho^2\sin^2(\theta)}\right)^{-1/2}, \label{eq:Lsq95est}\tag{35}\] with \(\lambda\) defined in ?? . This predictor is used to guide adaptive subdivision of the polar interval in the case \(p=1/2\), where no singularity-swap formulation is applied.

We begin by deriving a predictor applicable to a general \(n\)-point quadrature rule, expressed in terms of its error function \(k_n(t)\), and then specialize the construction to the Gauss–Legendre rule. The resulting predictor is summarized in Error predictor 1.

5.2.1 General results↩︎

Consider the integral \[\int_{\tha}^{\thb} \Lsq(\theta,\xx)\dtheta = \int_E \Lsq(\theta(t,\tha,\thb),\xx)\thsc\dt, \label{eq:Bintegral}\tag{36}\] where \(E=[-1,1]\), \(\theta(t,\tha,\thb)\) is the linear map defined in 4 , and \(\thsc\) is the associated scaling factor. We denote the quadrature error of this integral, evaluated with an \(n\)-point rule by \(\E[\Lsq](\xx,n,\tha,\thb)\). Let \(\{\trootlambda,\overline{\trootlambda}\}\) be the pair of branch points of the integrand closest to \(E\), defined as the roots of \(\Rlambda(\theta(t,\tha,\thb),\xx)\) in ?? . For notational convenience, we write \(\theta(t,\tha,\thb)=\thetat\) and denote these roots by \(\{\troot,\overline{\troot}\}\) for the remainder of Section 5.

From the contour representation reviewed in Section 5.1, the quadrature error admits an exact contour integral expression. By deforming the contour away from the real interval and retaining only the contributions associated with the nearest branch points, the error is modeled by the sum of two branch-cut integrals, \[\E[\Lsq](\xx,n,\tha,\thb) \approx \frac{1}{2\pi i}\int_{C_1}\Lsq(\thetat,\xx)k_n(t)\thsc\dt + \frac{1}{2\pi i}\int_{C_2}\Lsq(\thetat,\xx)k_n(t)\thsc\dt\eqqcolon \E_1+\E_2,\] where \(C_1\) and \(C_2\) are the curves going around \(\troot\) and \(\overline{\troot}\), respectively (see Figure 4).

We focus on the contribution \(\E_1\) from the curve \(C_1\). The integrand contains square-root branch-point singularities at \(\troot\) and \(\overline{\troot}\). By isolating these factors, \(\E_1\) can be written as \[\E_1 = \frac{\thsc}{2\pi i}\int_{C_1} \frac{k_n(t)}{\left(\lambda(\thetat)+\Gbar(\thetat)^{1/2}\thsc(t-\troot)^{1/2}(t-\overline{\troot})^{1/2}\right)^{1/2}}\dt, \label{eq:Lsq95quad95E1}\tag{37}\] where \[\Gbar(\theta) \mathrel{\vcenter{:}}= \frac{\lambda(\theta)^2-\rho^2\sin^2(\theta)}{(\theta-\thetaroot)(\theta-\overline{\thetaroot})}. \label{eq:Gbar}\tag{38}\] Let \(C_1^+\) denote the side of \(C_1\) going outwards, from \(\troot\) to \(\infty\), and let \(C_1^-\) denote the other side going in the opposite direction such that \(C_1=C_1^-\cup C_1^+\), illustrated in Figure 4. Defining the jump across the branch cut as \[(t-\troot)\Big.\Big\vert_{C_1^+} = (t-\troot)\Big.\Big\vert_{C_1^-}e^{-2\pi i},\] we obtain \[(t-\troot)^{1/2}\Big.\Big\vert_{C_1^+} = -(t-\troot)^{1/2}\Big.\Big\vert_{C_1^-}.\] This allows \(\E_1\) to be written as the difference of two integrals: an integral from \(t_0\) to \(\infty\) with \((t-\troot)^{1/2}\Big.\Big\vert_{C_1^+}\), and another from \(\infty\) to \(t_0\) with \((t-\troot)^{1/2}\Big.\Big\vert_{C_1^-}\). With the definition \[J_{\pm}(\troot,n) \mathrel{\vcenter{:}}= \int_{\troot}^\infty \frac{k_n(t)}{\left(\lambda(\thetat)\pm\Gbar(\thetat)^{1/2}\thsc(t-\troot)^{1/2}(t-\overline{\troot})^{1/2}\right)^{1/2}}\dt, \label{eq:Jpm}\tag{39}\] we obtain \[\E_1 = \frac{\thsc}{2\pi i}\left(J_+(\troot,n) - J_-(\troot,n)\right) \eqqcolon \frac{\thsc}{2\pi i}J(\troot,n) \label{eq:J}\tag{40}\] The contribution \(\E_2\) from \(C_2\) for the branch point \(\trootbar\) is computed analogously and satisfies \(\E_2=\overline{\E_1}\). Combining the two contributions, discarding the real part, and taking the absolute value, yields the following general quadrature error predictor \[\left|\E[\Lsq](\xx,n,\tha,\thb)\right| \approx \frac{\thsc}{\pi}\left|J(\troot,n)\right|. \label{eq:general95quad95est95B}\tag{41}\] To render this predictor computationally useful, an efficient and sufficiently accurate way to evaluate the integrals defining \(J_\pm\) is needed. We now consider that problem in the case of the Gauss–Legendre rule.

5.2.2 Gauss–Legendre rule↩︎

For Gauss–Legendre quadrature on \(E=[-1,1]\), the error function \(k_n\) does not have a closed form expression, but was shown in [33], [34] to satisfy the asymptotic formula \[k_n(z) \simeq \frac{2\pi}{\xi(z)^{2n+1}} \label{eq:GL95knz}\tag{42}\] as \(n\rightarrow\infty\) with \(\xi(z) = z + \sqrt{z^2-1}\). Here, and for the remainder of this paper, we shall use the symbol \(\simeq\) to mean “asymptotically equal to”. That is, \(a(n)\simeq b(n)\), for large \(n\), if \(\lim_{n\rightarrow\infty}a(n)/b(n)=1\). Moreover, we will define \(\sqrt{z^2-1}\) as \(\sqrt{z+1}\sqrt{z-1}\) with \(-\pi<\arg(z\pm1)\leq\pi\).

The integrands defining 39 contain a branch point at \(\troot\). Following [34], we define a branch cut extending from \(\troot\) to \(\infty\) that does not intersect \([-1,1]\), \[B(\troot) = \left\{\branchparam(s)\in\mathbb{C} : \branchparam(s) = \frac{1}{2}\left(\zeta(s)+\frac{1}{\zeta(s)}\right),\quad 1\leq s<\infty\right\}, \label{eq:branch95cut}\tag{43}\] where \[\zeta(s) = \zeta_0s,\quad \zeta_0 = \troot + \sqrt{\troot^2-1}. \label{eq:branch95cut95extra}\tag{44}\] One may verify that, as intended, \(\branchparam(1)=\troot\). This particular parametrization leads to the following helpful simplification of \(k_n\), \[k_n(\branchparam(s))\simeq \frac{2\pi}{\zeta(s)^{2n+1}}, \label{eq:kn95branch95cut}\tag{45}\] where \[\bernsteinr(\troot)=|\zeta_0| \label{eq:bernstein95radius}\tag{46}\] is the so-called Bernstein radius of \(\troot\).

Substituting this asymptotic form into the expressions for \(J_\pm\), we obtain the following model for the branch-cut contributions \[J_\pm(\troot,n) \simeq \frac{2\pi}{\zeta_0^{2n+1}}\int_1^\infty \frac{1}{s^{2n+1}} \frac{\branchparamprime(s)\ds}{\left(\lambda(\theta(\branchparam(s))\pm\Gbar(\theta(\branchparam(s)))^{1/2}\thsc(\branchparam(s)-\troot)^{1/2}(\branchparam(s)-\overline{\troot})^{1/2}\right)^{1/2}}. \label{eq:Jbranch}\tag{47}\] Inserting this expression into the general predictor derived in the previous subsection leads directly to the Gauss–Legendre quadrature error predictor stated below.

5.3 Interpolation error predictor for adaptivity in polar direction↩︎

We now construct an interpolation error predictor for polynomial interpolation of the function \(\Lsq\) defined in 35 . As in the quadrature case, the purpose is not to obtain a sharp upper bound, but a reliable efficiently computed quantity suitable for driving adaptive panel refinement. The interpolation is performed on the base interval \(E=[-1,1]\), with the linear map \(\theta(t,\tha,\thb)\) defined in 4 , using the \(n\) Gauss–Legendre nodes.

Starting from the complex contour representation of the interpolation error reviewed in Section 5.1, and using the same notation as in Section 5.2, the pointwise interpolation error at \(\tau\in E\) admits the representation \[\left|\EI_n[\Lsq](\xx,\tau)\right| \approx \frac{1}{\pi}\left|\elln(\tau)\left(\JI_+(\troot,n)-\JI_-(\troot,n)\right)\right|, \label{eq:EI}\tag{48}\] where the branch-cut integrals \(\JI_\pm\) are defined as \[\JI_\pm(\troot,n) \approx \int_1^\infty\frac{1}{(\branchparam(s)-\tau)\elln(\branchparam(s))}\frac{\branchparamprime(s)\ds}{\left(\lambda(\theta(\branchparam(s))+\Gbar(\theta(\branchparam(s)))^{1/2}\thsc(\branchparam(s)-\troot)^{1/2}(\branchparam(s)-\overline{\troot})^{1/2}\right)^{1/2}}. \label{eq:JIbranch}\tag{49}\] The main difference between 47 and 49 is the factor \(1/\elln(\branchparam(s))\), which encodes the dependence on the interpolation nodes. For Gauss–Legendre nodes, this factor can be characterized explicitly.

Along the branch cut \(\branchparam\) 43 associated with the singularity \(t_0\), the right-most term in ?? simplifies to \((\zeta_0s)^{n+1/2}\), \(1\leq s<\infty\). Thus the interpolation error inherits an exponential factor \(|\zeta_0|=\varrho(t_0)^{-(n+1/2)}\), recalling 46 , in contrast to the quadrature case, where the decay is governed by \(\varrho(t_0)^{-(2n+1)}\).

By combining Lemma 5 with 49 and maximizing over \(\tau\in\{-1,0,1\}\), we obtain the final interpolation error predictor in ?? . This maximization is supported by numerical observations.

5.4 Numerical evaluation of branch-cut integrals↩︎

Both the quadrature and interpolation error predictors reduce to branch-cut integrals \(\betaQ\) and \(\betaI\) in ?? and ?? , respectively. These integrals take the form \[\int_1^\infty s^{-m}F(s)\ds,\] with \(m=2n+1\) in the quadrature case and \(m=n+1/2\) in the interpolation case. The factor \(s^{-m}\) introduce rapid algebraic decay also for moderate \(n\), localizing the integral near \(s=1\).

We therefore truncate the interval \([1,\infty)\) to \([1,L]\), where \(L\) is chosen such that \(s^{-m}=\kappa\) at \(s=L\). In all experiments, we set \(\kappa=10^{-10}\). The remaining factor \(F(s)\) is found to be surprisingly smooth for all target point locations. The truncated integral is therefore efficiently evaluating using Gauss–Legendre quadrature on \([1,L]\), with eighth nodes typically being sufficient for both predictors.

The resulting computational cost is negligible compared to the surface quadrature itself, ensuring that the adaptive refinement remains efficient.

5.5 Performance of error predictors for adaptivity↩︎

We now assess the performance of the derived quadrature and interpolation error predictors for the function \(\Lsq\) in 35 . The aim is to verify that the predictors correctly capture the dependence on target point location, reproduce the expected decay with increasing \(n\), and remain reliable in the close-evaluation regime.

Figure 5 compares measured and predicted errors for a fixed value of \(n\), for both a convex spheroid and a non-convex “smooth star” geometry. The results show that the predictors perform well for both geometries across all target point locations, including those close to the symmetry axis.

We next examine the decay of the error as \(n\) increases. Since Gauss–Legendre error function 42 is constant along the level sets of the Bernstein radius function \(\bernsteinr\), the quadrature error is approximately constant on these level sets. We therefore construct target points by determining roots \(\trootlambda\) such that \(\varrho(\trootlambda)=P\), for a fixed constant \(P\), and then form the corresponding target points \(\xx(\trootlambda)\) following [31]. These target points are shown in black in Figure 6 (a).

Figure 6 (b) presents the maximum absolute error over these targets as a function of \(n\). Both quadrature and interpolation errors exhibit exponential decay, with rates consistent with their theoretical dependence on the Bernstein radius. The predictors closely follow this decay, capturing both the rate and magnitude of the true error.

To assess robustness in the close-evaluation regime, additional target points are placed along three lines normal to the surface, as illustrated in 6 (a). The lower panels of Figure 5, with \(n=28\) fixed, show the error as a function of target location and distance to the surface. The predictors remain stable across this range and provide reliable approximate upper bounds for both distant and near-surface target points.

6 Root finding↩︎

The error predictors in Section 5 depend on the location of the nearest complex singularities of the relevant integrands. In our setting, these singularities are characterized by complex roots of squared-distance functions.

For a target point \(\xx\), three types of roots arise:

1. The azimuthal root \(\varphi_0(\theta,\xx)\) of the full squared-distance function \(R(\theta,\varphi,\xx)^2\) for fixed \(\theta\). This root is used in singularity swap quadrature (SSQ) in the azimuthal variable and is also required for evaluating the quadrature error predictor of the layer potential 2 . For axisymmetric surfaces, \(\varphi_0\) is available in closed form and is given in Lemma 2.

2. The polar root \(\theta_0(\varphi,\xx)\) of \(R(\theta,\varphi,\xx)^2\) for fixed \(\varphi\). This root is required to evaluate the quadrature error predictor of the layer potential 2 , which is used to determine whether special quadrature is needed to achieve the desired accuracy.

3. The polar root \(\theta_0^\lambda(\xx)\) of the reduced squared-distance function \(R_\lambda(\theta,\xx)^2\) in ?? . This root governs the error predictors for \(\Lsq\) used to drive adaptivity in the polar direction and is also used in SSQ in the polar variable.

In each case, the relevant root is the one closest to the corresponding integration interval. For the azimuthal root this interval is \([0,2\pi)\), whereas for the polar roots it is either the full interval \([0,\pi]\) or, after adaptive refinement, the current polar subpanel mapped to \([-1,1]\). This closest-root principle reflects the fact that the contribution of a singularity to the quadrature and interpolation error decays exponentially with its Bernstein radius, which is the mechanism underlying the predictors in Section 5.

For spheroidal geometries, \(\theta_0^\lambda(\xx)\) and \(\theta_0(\varphi,\xx)\) can be obtained analytically. We begin by presenting these derivations, followed by a discussion of the general axisymmetric case.

6.1 Analytical roots for spheroidal geometries↩︎

We consider the spheroid parametrized by \[\ggamma(\theta,\varphi) = \left(a\sin(\theta)\cos(\varphi),~a\sin(\theta)\sin(\varphi),~b\cos(\theta)\right),\quad 0\leq\theta\leq\pi,~0\leq\varphi<2\pi, \label{eq:gamma95spheroid}\tag{50}\] with constants \(a,b>0\). The derivation of the formula in Lemma 6 follows from reducing the problem to a planar circle and ellipse using a Joukowsky transform. For completeness, these auxiliary results are provided in Appendix 11.

The following result extends the derivation in [32], where the analysis was carried out for the special case of a sphere.

6.2 General axisymmetric geometries↩︎

The closed-form formulas above apply only to spheroids. For a general axisymmetric surface, the azimuthal root \(\varphi_0(\theta,\xx)\) remains available analytically (Lemma 2), while the polar roots \(\theta_0^\lambda(\xx)\) and \(\theta_0(\varphi,\xx)\) are computed numerically.

For instance, to the \(\theta_0(\varphi,\xx)\) of \(R(\theta,\varphi,\xx)^2=0\) at fixed \(\varphi\), we use a complex Newton iteration, following the approach outlined in [32]. The initial guess is chosen from the polar angle \(\theta^*\) corresponding to the discretization point on the surface closest to the target \(\xx\). We then set the initial guess as \(\theta^*+i\beta\). With \(\beta=0.1\), convergence to a strict tolerance is usually achieved within 10 iterations.

The same procedure is used to compute \(\theta_0^\lambda(\xx)\), replacing \(R(\theta,\varphi,\xx)^2\) by \(R_\lambda(\theta,\xx)^2\). In the rare cases where convergence is not achieved, we restart the iteration and use over-relaxation.

7 Cost model for close evaluation↩︎

Here we briefly discuss the per-target cost of the presented algorithm. The goal is not to provide a sharp complexity theorem with explicit theoretical guarantees, but rather to identify the dominant operations in the current implementation and how they scale with the surface discretization and the adaptive refinement in the polar direction. The current MATLAB implementation is not well suited for a detailed throughput study, which is left for future work.

We assume a base surface discretization consisting of an \(\nphi\)-point trapezoidal rule in the azimuthal direction and an \(\nt\)-point Gauss–Legendre panel in the polar direction, with the layer density known at these nodes. If several base Gauss–Legendre panels are used in the polar direction, the discussion below applies panelwise. The cost of constructing the base discretization, solving for the layer density, and performing one-time precomputations is not included.

After adaptive refinement in the polar direction, the target-dependent refined grid contains \(\npan\nGL\) polar nodes for each of the \(\nphi\) azimuthal values, giving a total of \(\nphi\npan\nGL\) refined surface points. In the current implementation, a substantial part of the cost comes from interpolating the layer density from the base polar nodes to this refined grid. For each fixed azimuthal node, this is a one-dimensional interpolation from \(\nt\) source nodes to \(\npan\nGL\) target nodes, resulting in a total interpolation cost of \[\mathcal{O}(\nphi\npan\nGL\nt).\]

Let \(\npanSQphi\) denote the number of refined panels on which the singularity swap quadrature (SSQ) method is used in the azimuthal direction, and let \(\npanTZ\) denote the number of panels on which the trapezoidal rule is sufficient, so that \(\npan=\npanSQphi+\npanTZ\). For each of the \(\npanSQphi\nGL\) polar nodes where azimuthal SSQ is applied, a one-dimensional \(\nphi\)-point FFT is computed, giving a total cost of \[\mathcal{O}\left(\nphi\log(\nphi)\npanSQphi\nGL\right).\] After the FFTs have been computed on the panels requiring azimuthal SSQ, the azimuthal integral must still be evaluated at every refined polar node. This costs \(\mathcal{O}(\nphi)\) per node, both for SSQ weights and for the trapezoidal rule, and therefore contributes \[\mathcal{O}(\nphi\npan\nGL).\] The remaining operations in the polar direction entail evaluation of one-dimensional integrals, resulting in computational cost that is subdominant relative to these steps.

Thus, for a fixed target point \(\xx\), the dominant per-target cost in the current implementation is described by \[\mathcal{O}(\nphi\npan\nGL\nt) + \mathcal{O}\left(\nphi\log(\nphi)\npanSQphi\nGL\right) + \mathcal{O}\left(\nphi\npan\nGL\right). \label{eq:full95cost}\tag{52}\] The number of \(\nGL\)-point Gauss–Legendre panels, \(\npan\), determined by the adaptive algorithm for a given value of \(\nGL\), depends primarily on the target location \(\xx\) and the prescribed error tolerance \(\eps\). The last term in 52 has the same panel-dependent factor as the interpolation term, but lacks the additional factor \(\nt\), and is therefore absorbed into the first term. Accordingly, the dominant-cost model can be summarized as \[C_1(\eps,\xx,\nGL)\nt\nphi + C_2(\eps,\xx,\nGL)\nphi\log(\nphi), \label{eq:simplified95cost}\tag{53}\] where \(C_1,~C_2>0\) capture the target-dependent refinement and the distribution of panels requiring SSQ in the azimuthal direction.

8 Numerical experiments↩︎

We now numerically demonstrate the accuracy, adaptivity, and error-control behavior of the proposed singularity swap surface quadrature (S3Q) method in Section 4. We consider several layer potentials over different smooth axisymmetric geometries and seek to evaluate them to a prescribed error tolerance \(\epsilon\). Both single-panel and multi-panel base discretizations in the polar direction are used. Unless stated otherwise, a single-panel discretization is employed. In the multi-panel case, the geometry is partitioned into segments in the polar direction, which we refer to as patches.

All code for this paper is written and run in MATLAB. This implementation is intended to demonstrate accuracy and adaptive behavior rather than optimized runtime performance.

8.1 Experimental setup↩︎

We consider layer potentials arising from Laplace, Helmholtz, and Stokes problems, all of which fit into the generic form 1 for appropriate choices of the kernel factor \(k(\xx,\yy)\) and singularity power \(p\). In all numerical experiments, the surface parametrization is assumed to be analytically known. The discretization parameters \(\nt\) and \(\nphi\) are chosen so that the geometry and sampled surface values are well resolved on the underlying tensor-product grid, which we refer to as the base discretization. The layer density \(\sigma\) is prescribed analytically so that it can be sampled on this grid, evaluated at complex roots for the predictors, and used to compute reference solutions. Within S3Q, the sampled density values on the base discretization are then interpolated onto the adaptively refined polar grid when needed.

Given a target point \(\xx\) with the minimum distance \(d\) from the surface \(S\), we use the quadrature error predictor of [32] to determine whether the base discretization is sufficient to evaluate \(u(\xx)\) to the prescribed tolerance \(\epsilon\). If the predicted error exceeds \(\epsilon\), S3Q is activated and new quadrature weights are computed through the adjoint formulation described in Section 4.4. In particular, S3Q combines standard quadrature and singularity swap quadrature (SSQ) in the azimuthal and polar directions, while also allowing adaptive refinement in the polar direction using \(\nGL\)-point Gauss–Legendre panels. The number of such panels is denoted by \(\npan\).

The error in evaluating \(u(\xx)\) is measured against a reference solution. For scalar-valued layer potentials, we report the absolute error between the numerical and reference solution. For Stokes potentials, which are vector-valued, we report the maximum componentwise absolute error between the numerical and reference solution. In both settings, we refer to this as the absolute error. The reference solution is computed using an adaptive quadrature routine.

For the Laplace experiments, the target points are sampled on a uniform \(M\times M\) grid in the \(xz\)-plane with \(M=200\) and fixed \(y=0\), with points inside the geometry discarded. In the Helmholtz experiment, we consider two similar target planes, with \(M=120\). In the Stokes experiment, we use the same type of target grid as in the Laplace experiments, with the only difference being that we consider targets both inside and outside the geometry.

In Section 8.4, we also show that S3Q can be applied patchwise and coupled to a fast summation method, namely the fast multipole method (FMM) [42]–[46], to accelerate the standard quadrature. Specifically, we use the point-based FMM from the FMM3D package.⁸

8.2 Harmonic single layer potential: Type S3↩︎

As a first example, we consider the harmonic single layer potential, corresponding to \(p=1/2\), \[u(\xx) = \int_S \frac{\sigma(\yy)}{\|\yy-\xx\|}\dS(\yy) \label{eq:laplace95single95layer}\tag{54}\] on a type S3 spheroid. Here a type \(SX\) spheroid denotes a spheroid with aspect ratio \(1:X\), with semi-axes \(a(\theta)=1\) and \(b(\theta)=X\). The layer density function is chosen as \[\sigma(\ggamma(\theta,\varphi))=\sin(5\theta)e^{-\cos^2(\varphi)}+1.03, \label{eq:sigma95spheroid95s3}\tag{55}\] and the potential is evaluated for the three error tolerances \(\eps=10^{-4},~10^{-6},~10^{-8}\).

Figure 7 presents the results with \(\nt=\nphi=40\) and \(\nGL=32\). The figure shows that the standard quadrature error predictor accurately identifies the target points where S3Q must be used, and that the resulting measured errors at those points are close to or below the prescribed tolerance \(\eps\). The maximum absolute error over all target points is approximately \(1.2\times 10^{-4}\), \(5.0\times 10^{-8}\), and \(5.3\times 10^{-9}\) for the tolerances \(\eps=10^{-4},~10^{-6},~10^{-8}\), respectively.

As \(\eps\) decreases, the number of Gauss–Legendre panels \(\npan\) selected by the adaptive algorithm in the polar direction increases, as shown in Figure [fig:laplace95slp95s395histogram95nGL32]. Despite the close proximity of target points to the surface, the maximum \(\npan\) required to achieve the desired error tolerances \(\eps=10^{-4},~10^{-6},~10^{-8}\) was \(12\), \(13\) and \(14\), respectively. For \(\eps=10^{-8}\), this maximum was reached at only two target points situated approximately \(5.2 \cdot 10^{-5}\) from the surface.

8.3 Harmonic double layer potential: Type S10↩︎

We next demonstrate the performance of S3Q for a stronger singularity and illustrate the adaptive refinement behavior for a highly elongated geometry. Specifically, we consider the harmonic double layer potential, corresponding to \(p=3/2\), \[u(\xx) = \int_S \frac{\nn_\yy\cdot(\yy-\xx)~\sigma(\yy)}{\|\yy-\xx\|^3}\dS(\yy) \label{eq:laplace95double95layer}\tag{56}\] evaluated on a type S10 spheroid. Here, \(\nn_\yy\) denotes the outward-pointing normal vector at \(\yy\in S\). The prescribed tolerances are \(\eps=10^{-4},~10^{-6},~10^{-8}\), and the layer density is chosen as \[\sigma(\ggamma(\theta,\varphi))=1+\sin(6\varphi+\theta)\sin(\theta)^2. \label{eq:sigma95spheroid95s10}\tag{57}\] The base discretization uses \(n_t=160\), \(n_\varphi=100\), and we set \(\nGL=16\).

Figure 8 shows the result for error tolerance \(\eps=10^{-8}\). The measured absolute error is generally kept below the tolerance, with only 44 out of 30354 target points slightly exceeding it, and with a maximum error \(2.4\times 10^{-8}\) over all target points. For \(\eps=10^{-4}\) and \(\eps=10^{-6}\), the corresponding maximum absolute errors are \(8.3\times 10^{-5}\) and \(1.6\times 10^{-6}\), respectively. The figure also shows how \(\npan\) varies with the location of the target point relative to the surface.

This behavior is examined further in Figure 9, where the total number of quadrature points in the polar direction, \(\npan\times\nGL\), is plotted against the minimum distance \(d\) from the target point to the surface for the three prescribed tolerances. The required number of points grows approximately log-linearly with \(d\), and increases with stricter tolerances.

Note that the surface mesh shown in Figure 8 represents the base discretization with \(\nt=160\) and \(\nphi=100\), whereas each adaptively introduced Gauss–Legendre panel in the polar direction contains only \(\nGL=16\) points. Thus, \(\npan=10\) corresponds to the same number of quadrature points in the polar direction as the base discretization.

8.4 Helmholtz double layer potential: Multiple type S3↩︎

We next consider a Helmholtz double layer potential and illustrate how S3Q can be applied patchwise and coupled to fast summation by using the near-quadrature correction together with an FMM-accelerated standard quadrature. We consider the Helmholtz double layer potential on a collection of type S3 spheroids, denoted by \(\mathcal{S}\), given by \[u(\xx) = \frac{1}{4\pi}\int_{\mathcal{S}} \frac{\left(1-i\omega\|\yy-\xx\|\right)e^{i\omega\|\yy-\xx\|}~\nn_{\yy}\cdot(\xx-\yy)~\sigma(\yy)}{\|\yy-\xx\|^3}\dS(\yy),\] where \(\omega\) denotes the wavenumber and \(\nn_{\yy}\) the outward-pointing unit normal at \(\yy\in\mathcal{S}\).

To apply S3Q, we split the Helmholtz kernel into a non-smooth part and a smooth part, following [13]. More precisely, we write the Helmholtz double layer kernel as \[\mathcal{D}^H(\xx,\yy) = D(\xx,\yy)\left(H_1(\xx,\yy)+iH_2(\xx,\yy)\right),\] where \[D(\xx,\yy) = \frac{\nn_\yy\cdot(\xx-\yy)}{4\pi\|\yy-\xx\|^3}\] and \[H_1(\xx,\yy) = \cos(\omega\|\yy-\xx\|) + \omega\|\yy-\xx\|\sin(\omega\|\yy-\xx\|),\] \[H_2(\xx,\yy) = \sin(\omega\|\yy-\xx\|) - \omega\|\yy-\xx\|\cos(\omega\|\yy-\xx\|).\] The factor \(DH_1\) contains the non-smooth part of the kernel, while \(DH_2\) is smooth. Indeed, expanding \(H_1\) about \(\|\yy-\xx\|=0\) shows that \(DH_1\) consists of one contribution of type \(p=1/2\) and one contribution of type \(p=3/2\), and therefore falls within the class of nearly singular kernels treated by S3Q. Thus, S3Q is applied only to these two non-smooth contributions of \(DH_1\), whereas the smooth contribution \(DH_2\) is evaluated by the standard quadrature, accelerated by the FMM.⁹

The geometry consists of the cluster of 40 type S3 spheroids shown in Figure [fig:helmholtz95s3q95err]. The base discretization for each spheroid has three panels in the polar direction, each with \(\nt=10\), and \(\nphi=26\) points in the azimuthal direction. Figure [fig:helmholtz95spheroid] shows one such spheroid. In this experiment, the wavenumber is set to \(\omega=4.2\), the layer density is chosen as \[\sigma(\yy) = 2+\sin(y_3),\quad y=\left(y_1,y_2,y_3\right),\] and the prescribed tolerance is \(\eps=10^{-6}\). The full Helmholtz double layer potential is first evaluated using the FMM, after which S3Q is used as a local near-quadrature correction for the two non-smooth terms in \(DH_1\). Both the S3Q corrections and the quadrature error predictors of [32] are applied patchwise, so that each patch is treated independently when determining whether correction is required.

Figure [fig:helmholtz95s3q95err] shows the result. Of the 28800 target points in the grid, 27193 lie outside the particles, and 5257 of these exterior targets, or 19.3%, are classified as requiring S3Q. The maximum absolute error over all exterior target points is \(3.0\times 10^{-6}\). Excluding the cost of evaluating the quadrature error predictors, approximately 28% of the runtime is spent in the FMM and 72% in constructing the S3Q near-quadrature weights.

8.5 Stokes double layer potential: Axisymmetric geometry↩︎

We now evaluate the Stokes double layer potential \(\uu\), corresponding to \(p=5/2\), on an axisymmetric “peanut” geometry using S3Q. The three components of \(\uu\) are given by¹⁰, for \(i=1,2,3\), \[u_i(\xx) = \int_S \mathcal{T}_{ijk}(\xx-\yy)\sigma_j(\yy)n_k(\yy)\dS,\quad \mathcal{T}_{ijk}(\rr) = -6\frac{r_ir_jr_k}{\left\|\rr\right\|^5}, \label{eq:stokes95dlp}\tag{58}\] where \(\mathcal{T}_{ijk}\) is the stresslet tensor. We take the constant density \(\sigma_j(\yy)=1\), \(j=1,2,3\), for which the stresslet identity gives an exact reference solution. ¹¹ The prescribed tolerances are \(\eps=10^{-3}\) and \(\eps=10^{-6}\), and the base discretization uses \(\nt=320\), \(\nphi=180\), and we set \(\nGL=16\).

Figure 11 (a) and 11 (b) show the measured absolute error for \(\eps=10^{-3}\) and \(\eps=10^{-6}\), respectively. For \(\eps=10^{-3}\), the measured error remains close to the requested accuracy over nearly the entire target set, with a maximum absolute error of \(4.4\times 10^{-2}\), attained at the target point closest to the surface, for which \(d=7.9\times 10^{-5}\). For \(\eps=10^{-6}\), the same overall behavior is observed, except in a thin region nearest the surface where the prescribed accuracy is not achieved.

This loss of accuracy for very close targets is caused by catastrophic cancellation in the azimuthal SSQ step. The issue arises when the kernel has both a strong singularity and a near-vanishing numerator, and becomes more severe as the singularity strengthens. In this regime, the true value of the integral is small, while the SSQ quadrature sum is formed from much larger terms that nearly cancel, leading to substantial amplification of roundoff errors. This behavior was already noted in the original SSQ paper [21], and has recently been analyzed in detail in [24], where it is shown to be tied to the use of a fixed interpolation basis that does not reflect the local vanishing structure of the numerator.

To examine this more directly, let \(E_{\textrm{cancel}}\) denote the approximate cancellation error introduced in [24]. Figure 11 (c) shows that the targets for which \(E_{\textrm{cancel}}>\eps\) generally coincide with those whose errors exceed \(\eps\). Thus, the observed loss of accuracy for the closest targets is well explained by floating-point cancellation in the azimuthal SSQ step, rather than by a failure of the adaptive refinement or error-control mechanisms of S3Q. This is precisely the failure mode addressed by the translated-basis stabilization developed in [24].

9 Conclusions↩︎

We have presented an adaptive approach for evaluating nearly singular layer potentials on smooth axisymmetric surfaces, in which standard quadrature is used whenever possible and the proposed singularity swap surface quadrature (S3Q) is applied only where necessary. A key feature of the method is that all algorithmic choices and parameter values associated with the near-quadrature correction are determined automatically from the prescribed error tolerance.

The S3Q method combines singularity swap quadrature (SSQ) with a local adaptive discretization procedure. It requires only the prescribed error tolerance as essential user input, together with, optionally, the order of the Gauss–Legendre quadrature used for adaptive refinement in the polar direction. The near-quadrature correction employs SSQ in the azimuthal direction together with adaptive refinement in the polar direction, where either standard Gauss–Legendre quadrature or SSQ is used on the refined grid depending on the predicted quadrature error. After the azimuthal correction, the remaining integrand in the polar direction typically has a nearly singular behavior of one degree lower than that of the original layer potential kernel. This integral is then evaluated using SSQ only when the remaining near singularity demands it, and otherwise by standard Gauss–Legendre quadrature. A central component of the method is the use of quadrature and interpolation error predictors to control the local refinement. In particular, the interpolation error prediction extends the complex-variable framework of [31], previously used for quadrature error prediction in [32], and provides a practical mechanism for automatic parameter selection.

The numerical experiments show that this strategy provides reliable error control for harmonic single and double layer potentials, while requiring only moderate refinement in the polar direction even for targets close to the surface. They further show that the framework extends naturally to Helmholtz problems through kernel splitting, and that in the corresponding multi-body example the S3Q correction can be applied patchwise while the standard quadrature is accelerated through coupling to the fast multipole method.

For the Stokes double layer potential, the experiments show that the combination of standard quadrature and S3Q remains accurate over most of the close-evaluation region, but that a very thin innermost layer becomes dominated by floating-point cancellation in the azimuthal SSQ step. In the present experiments, this limitation is primarily visible for the \(1/r^5\) singularity and only for very close targets when high accuracy is requested.

A robust treatment of this regime is outside the scope of the present paper. However, the recent stabilization developed in [24] shows that this cancellation can be greatly reduced by replacing the standard SSQ basis with target-adapted translated bases. In the present surface setting, this suggests a natural next step: to combine the near-quadrature correction developed here with a stabilized azimuthal SSQ evaluator. This, together with a more optimized implementation of the method, is left for future work.

Acknowledgments↩︎

The authors would like to thank Ludvig af Klinteberg at the Department of Business and Mathematics, Mälardalen University (MDU), for helpful discussions. The authors acknowledge the support from the Swedish Research Council under grant 2023-04269.

10 Two proofs↩︎

11 Roots of planar squared-distance functions↩︎

This appendix derives the analytical formulas used in Section 6 for computing the complex root of the squared-distance function \(R_\lambda^2\) in spheroidal geometries. The key observation is that the spheroidal problem can be reduced to finding the roots of squared-distance functions for circles and ellipses in the plane.

We first recall the known result for a circle, which was previously derived in [32]. We then extend it to ellipses using a Joukowsky transform that maps ellipses to circles.

12 Stabilization of recurrence relations↩︎

As described in Remark [rem:mu95rec95stability], the forward recurrences for \(\mu_k^p(\chi)\) in ?? can become numerically unstable when \(\chi\) is small. The desired sequence is the minimal solution and decays approximately like \(\chi^k\). The recurrence also admits a complementary solution that grows with \(k\). In finite precision, even a small roundoff component in the growing solution can eventually dominate the computed values and obscure the desired solution.

When this instability is detected, we use the stabilization produce of [38], in which the three-term recurrence is recast as a homogeneous tridiagonal boundary-value problem of size \(k_0-1\). The endpoint \(k_0\) must be chosen large enough to include all Fourier modes required by the truncated expansion, and sufficiently far out that the minimal solution has decayed to a negligible size.

We use the following asymptotic model for the magnitude of the minimal solution: \[\mu_k^p(\chi)\approx \tilde{\mu}_k^p(\chi) \mathrel{\vcenter{:}}= C(p)\mu_0^p(\chi)\chi^k, \qquad 0<\chi<1. \label{eq:mu95est}\tag{77}\] The constants \(C(p)\) are chosen empirically for \(p=1/2,3/2,5/2\). This model is not used as a rigorous bound, but as a practical rule for choosing the size of the tridiagonal system and for deciding when the forward recurrence is expected to be reliable. The corresponding decay models and roundoff-growth indicators \(\Erecest(\chi)\) are listed in Table 1. In the implementation, the forward recurrence is used only if \(\Erecest(\chi)\), evaluated at the largest required mode, is below the requested recurrence tolerance. Otherwise, the tridiagonal stabilization is used.

Figure 12 compares the model 77 and the forward recurrence error predictor in Table 1 with accurate values of \(\mu_k^{3/2}(\chi)\) and with measured forward-recurrence errors. The agreement is sufficient for their intended roles: choosing the endpoint \(k_0\) and deciding when to switch from the forward recurrence to the stabilized tridiagonal solve. The same behavior is observed for \(p=1/2\) and \(p=5/2\).

It remains to choose the endpoint \(k_0\) of the tridiagonal system. Let \(\epsilon_{\textrm{rec}}\) denote the recurrence tolerance, which in the experiments is taken to be the requested accuracy tolerance. From the decay model 77 , define \[\kbar = \left\lceil\dfrac{\log\left(\dfrac{\epsilon_{\textrm{rec}}}{C(p)\mu_0^p(\chi)}\right)}{\log(\chi)}\right\rceil.\] Thus, according to the model, the minimal solution has then decayed to approximately \(\epsilon_{\textrm{rec}}\) by mode \(\kbar\). We then choose \[k_0 = \max\left\{\left\lceil S\kbar\right\rceil,\,\kmax+1,\,2\right\}, \label{eq:k0}\tag{78}\] where \(\kmax\) is the largest Fourier mode required by the truncated expansion 12 , and \(S\geq1\) is a safety factor. In the experiments we use \(S=1.5\).

References↩︎

---

1. Department of Mathematics, KTH Royal Institute of Technology, Stockholm, Sweden  (davkra@kth.se).↩︎

2. Department of Mathematics, KTH Royal Institute of Technology, Stockholm, Sweden  (akto@kth.se).↩︎

3. \(R^2(\theta,\varphi,\xx)\) is no longer a norm when evaluated with complex arguments, in which case we use the right-most expression in 6 .↩︎

4. For some target locations, such as points on the symmetry axis of an axisymmetric surface, no such complex roots exist because the distance function remains strictly positive for all complexified angles.↩︎

5. When \(\Gamma\) arises by fixing \(\theta\) on an axisymmetric surface, one simply has \(\boldsymbol{\xi}(\phi)=\ggamma(\theta,\phi)\).↩︎

6. The function \(\Fcal\) is defined in ?? , with \(\mu_k^p\) computed via the recurrence relations ?? . The initial values ?? ?? involve the complete elliptic integral of the first kind \(K(r^2)\). As the evaluation point \(\xx\) approaches the surface, the quantity \(\chi(\theta)=e^{-|\Im(\phiroot(\theta,\xx))|} \rightarrow 1^-\), and \(K(\chi^2)\) diverges logarithmically.↩︎

7. For Gauss–Legendre quadrature, singularities farther from \(E=[-1,1]\) are exponentially suppressed: if \(g\) extends analytically to a Bernstein ellipse of radius \(\varrho\), then \(\E_n[g]=\mathcal{O}(\varrho^{-2n})\) [36]. For example, with \(n=16\) and a singularity at \(z=0.5i\), a second singularity at \(2z\) contributes about \(10^{-6}\) times as much to the error. If two singularities are located at comparable distances from \(E\) and only one is retained in the predictor, the neglected contribution is typically of the same order of magnitude, resulting in a difference by at most a modest constant factor.↩︎

8. The FMM3D package can be found at https://github.com/flatironinstitute/FMM3D.↩︎

9. The Helmholtz single layer potential admits an analogous splitting; see [13]. The same idea also applies to Yukawa single and double layer kernels, which similarly enables the use of S3Q.↩︎

10. The Einstein summation convention is applied here, where summation over repeated indices in a term is implied.↩︎

11. For any constant density vector \(\tilde{\boldsymbol{\sigma}}\) the Stokes double layer potential satisfies \(\uu(\xx)=\mathbf{0}\) outside \(S\), \(\uu(\xx)=4\pi\tilde{\boldsymbol{\sigma}}\) on \(S\), and \(\uu(\xx)=8\pi\tilde{\boldsymbol{\sigma}}\) inside \(S\); see [47].↩︎