Estimation of quadrature errors for layer potentials evaluated near surfaces with spherical topology

Numerical simulations with rigid particles, drops, or vesicles constitute some examples that involve 3D objects with spherical topology. When the numerical method is based on boundary integral equations, the error in using a regular quadrature rule to approximate the layer potentials that appear in the formulation will increase rapidly as the evaluation point approaches the surface and the integrand becomes sharply peaked. To determine when the accuracy becomes insufficient, and a more costly special quadrature method should be used, error estimates are needed. In this paper, we present quadrature error estimates for layer potentials evaluated near surfaces of genus 0, parametrized using a polar and an azimuthal angle, discretized by a combination of the Gauss-Legendre and the trapezoidal quadrature rules. The error estimates involve no unknown coefficients, but complex-valued roots of a specified distance function. The evaluation of the error estimates in general requires a one-dimensional local root-finding procedure, but for specific geometries, we obtain analytical results. Based on these explicit solutions, we derive simplified error estimates for layer potentials evaluated near spheres; these simple formulas depend only on the distance from the surface, the radius of the sphere, and the number of discretization points. The usefulness of these error estimates is illustrated with numerical examples.


Introduction
We consider a generic layer potential over a regular surface S ⊂ R 3 , where 2p ∈ Z + and the evaluation (or target) point x ∈ R 3 is allowed to be close to, but not on, S. The functions k (x, y) and σ(y) as well as S are assumed to be smooth.When x is close to S, the integrand will be peaked around the point on S closest to x, implying that, while the integral is well defined analytically, it is difficult to resolve well numerically.
In paper [1], we derived estimates for the numerical errors that result when applying quadrature rules to such layer potentials.Specifically, we considered the panel based Gauss-Legendre quadrature rule and the global trapezoidal rule.The estimates that was derived have no unknown coefficients and can be efficiently evaluated given the discretization of the surface.The evaluation involves a local one-dimensional root finding procedure.In numerical experiments, we have found the estimates to be both sufficiently precise and computationally cheap to be practically useful.This means that they can be used to determine when the regular quadrature is insufficient for a required accuracy, and hence when a more costly special quadrature method must be invoked.In deriving these estimates, we assumed that the local (for Gauss-Legendre) or global (for trapezoidal rule) surface parametrization is such that the map between the parameter space and the surface coordinates is one-to-one.
For surfaces of genus 0, topologically equivalent to a sphere, it is quite common to use a global parametrization in two angles, i.e. spherical coordinates.At the poles of such a parametrization, the parameter to surface coordinate map is not one-to-one.Here, the derivative of the surface coordinate with respect to the azimuth angle vanish, and with it, the surface area element.For such parametrizations, the previously derived error estimates cannot be directly used for all evaluation points x.
It is quite common in applications of boundary integral equations to discretize surfaces with parametrizations based on spherical coordinates, and to that attach a Gaussian grid, i.e. a discretization with a Gauss-Legendre quadrature rule in the polar angle (or a mapping of the polar angle) and the trapezoidal rule in the periodic azimuthal angle.This has been used for simulations of Stokes flow with solid particles in e.g.[2,3], and with drops and vesicles [4][5][6][7].In the latter case, the deformable shapes are represented by spherical harmonics expansions.
In this paper, we consider the case of a general smooth surface of genus 0, parametrized using a polar and an azimuthal angle, discretized by a combination of the Gauss-Legendre quadrature rule and the trapezoidal rule as described above.Before we introduce the contributions of this paper, we will briefly describe the derivations of the estimates in paper [1], and the previous results that enabled that work.
Consider the following simple integrals with z 0 = a + ib, a, b ∈ R, b > 0, p ∈ Z + , and over a basic interval E, e.g.[−1, 1] for Gauss-Legendre and [0, 2π) for the trapezoidal rule.If b is small the integrands have two poles/one pole close to the integration interval along the real axis.The theory of Donaldson and Elliott [8] defines the quadrature error as a contour integral in the complex plane over the integrand multiplied with a so-called remainder function, that depends on the quadrature rule.Elliott et al. [9] derived error estimates for the error in the approximation of ( 2) with an n-point Gauss-Legendre quadrature rule.To estimate the contour integral, they used residue calculus for p = 1 and branch cuts for 0 < p < 1.In [10], af Klinteberg and Tornberg derived error estimates for both (2) and (3) for the Gauss-Legendre quadrature rule, for any p ∈ Z + .Corresponding results were derived also for the trapezoidal rule, but for integration over the unit circle.Previous studies on the trapezoidal rule include the survey by Trefethen and Weideman [11], and the error bound provided by Barnett in [12] for the quadrature error in evaluating the harmonic double layer potential.
In [13], a key step was taken to accurately estimate the quadrature errors for the Gauss-Legendre quadrature rule in the approximation of layer potentials in 2D, written in complex form.Introducing a parametrization of a smooth curve segment, γ(t) ∈ C, a typical form of an integral to evaluate is As compared to the estimates for the simple complex integral (3) above, the estimates derived for this integral require the knowledge of t 0 ∈ C such that γ(t 0 ) = z 0 .In practice, given the Gauss-Legendre points used to discretize the panel, a numerical procedure is used to compute t 0 .Note that not all layer potentials in 2D can be written in this form.Using the same techniques, remarkably accurate error estimates were derived also for layer potentials for the Helmholtz and Stokes equations, in [13] and [14], respectively.Now, let S in (1) be a curve Γ ∈ R 2 , Γ = γ(E), E ⊂ R. We can then write the layer potential in (1) in the equivalent form In the last step all the components that are assumed to be smooth have been collected in the function f (t), which has an implicit dependence on x.The location of the poles are in this case given by each t 0 ∈ C such that the denominator is zero.In [1], error estimates were derived for such integrals, both for the Gauss-Legendre and the trapezoidal quadrature rules.To evaluate these estimates, the pair of the complex conjugate roots {t 0 , t0 } closest to integration interval is needed, and is in practice found through a numerical root-finding procedure.The estimates will be stated in Section 3.
A key observation in [1] was that the error estimates for the numerical approximation of ( 5) can be derived the same way for Γ ∈ R d , d = 2, 3.The only difference is that R 2 (t, x) = γ(t) − x 2 will have three additive terms instead of two.Starting with the curve estimates in R 3 , error estimates for the prototype layer potential (1) were derived.
With S a two-dimensional surface in R 3 , parametrized by γ : All the components that are assumed to be smooth have here been collected in f (t, ϕ) which depends implicitly on x.In [1] error estimates were derived for the numerical approximation of ( 6) by composite Gauss-Legendre quadrature or global trapezoidal quadrature.Numerical examples were shown with tensor product quadrature rules based on Gauss-Legendre quadrature for surface discretizations of quadrilateral patches, and the tensor product trapezoidal rule for global discretizations.

Contributions and outline
In this paper, we will discuss the generalization of the results of paper [1] to the case of smooth surfaces topologically equivalent to a sphere, parametrized by γ A generic such surface can be represented by where C m ∈ C 3 and Y m (θ, ϕ) is the spherical harmonic function of degree and order m.
In section 3, we introduce the error estimates derived in [1] for the integral over a curve in R 2 or R 3 (5).In section 4 we then introduce the extension to a surface of genus 0 in R 3 for our specific discretizations, based on what was done in paper [1].In section 5 we derive analytical results for axisymmetric surfaces.We also show how, based on the explicit knowledge of the roots, it is possible to derive simplified error estimates for layer potentials evaluated near spheres.In section 6, we discuss how to numerically evaluate the estimates for a general surface with spherical topology, and in section 7 we show how the error estimates perform on different numerical examples.
3 Quadrature error estimates for curves in R 2  and R 3 Let us introduce the base interval E, which for the Gauss-Legendre quadrature will be [−1, 1] and for the trapezoidal rule [0, 2π).Consider an integral over such a base interval and an n-point quadrature rule with quadrature nodes {t } n =1 and corresponding quadrature weights {w } n =1 to approximate it, The error as a function of n will depend on the function g and the specific quadrature rule.We will consider closed curves for the trapezoidal rule and open curves (segments) with the Gauss-Legendre quadrature rule.We now introduce the squared distance function for a curve in R d (d = 2 or 3), to an evaluation point We will later evaluate this function also for t ∈ C, in which case we will use the right most expression.This expression can then evaluate as a complex number, and will no longer be a norm.Our integral of interest ( 5) can be written in the form and we want to estimate E n [Θ p ](x).As x is not on γ(t), we have R 2 (t, x) > 0 for t ∈ E. There will however be complex conjugate pairs of roots to R 2 (t, x), since R 2 (t, x) is real for real t.Let t 0 , t 0 be the pair closest to E, s.t.
Under the assumption that f is smooth, the region of analyticity of Θ p (t, x) is bounded by these roots.We will henceforth refer to them both as roots (of R 2 ) and singularities (of the integrand).They are in most applications not known a priori, but can be found numerically for a given target point x (see section 6.2).The quadrature error E n [Θ p ](x) can, following Donaldson and Elliott [8], be written as a contour integral in the complex plane over the integrand Θ p (t, x) multiplied with a so-called remainder function, that depends on the quadrature rule.If p is an integer, t 0 and t0 are pth order poles of the integrand.If p is a half-integer, the integrand has branch points at these singularities.
An important step in the derivation leading up to estimates of the quadrature error in [1], is to divide and multiply the integrand with the factor t − w, where w ∈ C is the singularity at the branch begin considered.This yields the singularity (t − w) −1 to consider in the complex plane, and introduces what is denoted the geometry factor.
The geometry factor G is, for an evaluation point x ∈ R d and w ∈ C a root of R 2 , defined as With these definitions, we are ready to state the error estimates from [1], for both the trapezoidal rule and the Gauss-Legendre quadrature rule.
Error estimate 1 (Trapezoidal rule) Consider the integral in (15) for an evaluation point x ∈ R d , with 2p ∈ Z + , where γ(E) is the parametrization of a smooth closed curve in R d where d = 2 or 3.The integrand is assumed to be periodic in t over the integration interval E = [0, 2π).The error in approximating the integral with the n-point trapezoidal rule can in the limit n → ∞ be estimated as Here, Γ(p) the gamma function, and the geometry factor G is defined in (17).The squared distance function is defined in (14), and t 0 , t 0 is the pair of complex conjugate roots of this R 2 (t, x) closest to the integration interval E.
Error estimate 2 (Gauss-Legendre rule) Consider the integral in (15) for an evaluation point x ∈ R d , with 2p ∈ Z + , where γ(E) is the parametrization of a smooth closed curve in The error in approximating the integral with the n-point Gauss-Legendre rule can in the limit n → ∞ be estimated as where √ z 2 − 1 is defined as Here, Γ(p) the gamma function, the geometry factor G is defined in (17), the squared distance function is defined in (14), and t 0 , t 0 is the pair of complex conjugate roots of this R 2 (t, x) closest to the integration interval E.
In paper [1], each estimate was written as two different estimates, one for positive integers p, and one for positive half-integers p.The derivation of the two estimates follows a different path, using residue calculus and branch cuts, respectively.However, using (p − 1)! = Γ(p) for integer p the two estimates can both be written in the form given above.
4 Quadrature errors near two-dimensional surfaces in R 3 Let us now consider the three-dimensional case, and the prototype layer potential (6).Here, S ⊂ R 3 is a two-dimensional surface parametrized by γ : As introduced in (8), we will specifically consider γ(t, ϕ) where In analogy to the squared distance function to a curve, as introduced in (14), we now introduce the squared distance function between the surface γ(t, ϕ) and the evaluation point x = (x, y, z), (20) Note that we will later evaluate R 2 also for complex arguments using the right most expression, in which case it is no longer a norm.With this, the integrand of ( 6) can be written The operators I[g], Q n [g] and E n [g] were introduced in the beginning of Section 3, with Here we use them with a subindex indicating if they are applied in the t or the ϕ direction, where it should be understood that Gauss-Legendre quadrature rule is applied in the t-direction, and the trapezoidal rule in the ϕ direction.For ease of notation, we will skip the brackets above, such that I t I ϕ Θ p means an integration of Θ p first in the ϕ and then in the t direction.We can then write the tensor product quadrature as and from here In this last step, we have neglected the quadratic error term, and used that E t,nt I ϕ = I ϕ E t,nt .This last fact follows from E t,nt I ϕ = (I t − Q t,nt )I ϕ combined with For a more detailed discussion, see [15].For some basic integrals, Elliott et al. [16] have shown that the quadratic error term that we here neglect can have an important contribution.As it is a higher order contribution, this is only true when the quadrature error is large, and we will derive our estimates without it, as was also done in [1].Explicitly writing out the first term in the right hand side of (23), we have where ϕ l = 2π(l − 1)/n ϕ , l = 1, . . ., n ϕ and w T Z l = 2π/n ϕ , ∀l.The term in the brackets (i.e.E ϕ,nϕ ) represents the quadrature error of the trapezoidal rule on the line L t that for a given t is defined as For short, we will denote this curve γ(t, •).For a fixed t, this is the quadrature error for the trapezoidal rule, for which an estimate is given in Error estimate 1.
The term I ϕ E t,nt Θ p can be written analogously to (24), simply swapping t and ϕ, introducing the Gauss-Legendre quadrature nodes and weights.The error estimate that needs to be integrated in this term is given in Error estimate 2.
To be able to distinguish if the geometry factor in the error estimate corresponds to γ(t, •) or γ(•, ϕ), we extend the definition of the geometry factor in (17) and denote We need to work with the absolute value of the error, and we will use the estimate Expanding back from this shorthand notation, this can be formulated as follows.
Error estimate 3 (Surface in R 3 ) Given an evaluation point x ∈ R 3 , consider the integral in ( 6) with 2p ∈ Z + , where S ⊂ R 3 is a two-dimensional smooth closed surface parametrized by γ : The integrand is assumed to be periodic in ϕ over the integration interval E 2 = [0, 2π).The error in approximating the integral with the n t point Gauss-Legendre rule in the t-direction and the nϕ-point trapezoidal rule in the ϕ direction is defined as where {t k , w GL k } nt k=1 and {ϕ l , w T Z l } nϕ l=1 are the Gauss-Legendre and trapezoidal rule quadrature nodes and weights.
Remark 1 For this error estimate to be useful, it should give a good approximation of the error already for moderate values of n t and nϕ.In practice we find this to be true as long as n t and nϕ are large enough for the surface to be well resolved.This will be discussed in the numerical results section.
Remark 2 Given x ∈ R 3 , it is not guaranteed that a root ϕ 0 (t, x) exists for all t ∈ [−1, 1], nor that t 0 (ϕ, x) exists for all ϕ ∈ [0, 2π).For example, for surfaces of spherical topology with a global parametrization as defined in Eq. ( 8), the squared distance function R 2 (t, ϕ, x) is independent of ϕ for t = −1, 1, and hence no root ϕ 0 exists.For axisymmetric surfaces, no root ϕ 0 exist for evaluation points along the axis of symmetry, interior or exterior to the surface.
The exposition in this section has followed what was done in [1], however combining integration by the Gauss-Legendre rule in one direction, and the trapezoidal rule in the other.In [1] discretizations of surfaces of genus 1 with either a global trapezoidal rule in both directions, or a panel based Gauss-Legendre rule were considered.When a panel based discretization is used, the error estimates for the panels closest to the evaluation point are added together.We cannot theoretically guarantee that the roots that we need for evaluation of the error estimate always exist, and there are in general no analytical formulas for the roots.In [1], the root finding is done numerically, and approximations to the integrals in (31) and (32) are made, using the fact that the error contribution is strongly localized to the region on the surface closest to the evaluation point.For each evaluation point x, only one root t 0 (ϕ * , x) and one root ϕ 0 (t * , x) are needed, where γ(t * , ϕ * ) is the grid point (quadrature node) on the surface closest to the evaluation point x.This approach works well apart from the rare occasions where the root finding algorithm fails for evaluation points quite far from the surface.
To understand how we can evaluate error estimates for surfaces of genus 0 with a global parametrization, we will first analytically consider the simpler case of an axisymmetric surface, at times further simplified to a sphere.We will then in Section 6 discuss the practical evaluation of the estimate, including how to approximate the remaining integrals and determine the roots as needed.

Analytical derivations for axisymmetric and spherical surfaces
In this section we will consider an axisymmetric surface, as parametrized by with For some results, we will simplify further and set a(θ) = b(θ) = a and consider a sphere of radius a.
The parametrization γ(t, ϕ) relates to this parametrization through a mapping θ = θ(t) as given in Eq. ( 8), with γ • = γ •,A .The map θ(t) will change the parametrization of the surface in t and hence yield different locations of the Gauss-Legendre quadrature nodes on the surface.In this section, we will keep this choice open to the extent possible, and state most results in θ and ϕ.
Note that the axis of symmetry for γ •,A (θ, ϕ) is the z-axis.For a surface of a different shift and orientation, the evaluation point x can be translated and rotated into a local coordinate system of the particle, and the results stated below will apply.
Generally, the roots of R 2 cannot be found analytically, and we need to compute them using a root finding procedure.For the axisymmetric case, we can however analytically find the roots ϕ 0 given θ, and for a sphere, we can furthermore find the roots θ 0 given ϕ.
We will see that the estimate for the error incurred by the trapezoidal rule cannot be evaluated for an evaluation point at the symmetry axis, as the integrand in Eq. (31) becomes undefined.Using the analytical expressions for the roots, we can study appropriate limits as the evaluation point approaches the symmetry axis.We will also use these analytical results to derive a simplified error estimate for the sphere.

Analytical roots to the squared distance function
We will start by finding the roots to the squared distance function defined with respect to a circle in the xy plane, and then extend this result.
Remark 3 Note that if α 0 is a root to R 2 (α, x), so is α 0 + 2πp for any p ∈ Z. Further, notice that we have where x = ρ cos( θ), y = ρ sin( θ) with θ = atan2(y, x) have been introduced in the last step.
Here, (a/ρ + ρ/a)/2 ≥ 1 with equality only when ρ = a.Since x is not on the curve, we cannot have z = 0 in this case, and hence λ > 1.
Remark 4 Note that the root (in Lemma 2) is not defined if , which is independent of ϕ and always positive, so no root can be found.Similarly, if θ = 0 or π, R 2 ( θ, ϕ, x) is again independent of ϕ, and no root can be found.
, and rewrite Eq. ( 35) as Eq. ( 41) describes a circle of radius ã in the xy plane at z = b, and so we can proceed similarly as we did for Lemma 1.In this case, we get , where With this λ, the expression of the analytical root for ϕ 0 ( θ) is of the same form as in Eq. ( 36), as given in Eq. (40).
Remark 5 Note that the root in Lemma 3 is not defined for the case when z = 0 and φ − atan2(y, x) = π/2 + pπ, p ∈ Z.In this case, R 2 (θ, φ, x) = a 2 + x 2 + y 2 , which is independent of θ and always positive, so no root can be found.
Corollary 1 Under the assumptions of Lemma 3, with the evaluation point at the z-axis inside or outside of the sphere, x = (0, 0, z) ∈ R

Error estimates for evaluation points close to and on the symmetry axis
With γ(t, ϕ) = γ •,A (θ(t), ϕ), E T Z (γ, f, p, n ϕ , x) from Error estimate 3 can be written as where Here ϕ 0 (θ, x) is the root associated with γ •,A (θ, ϕ).Similarly, G • γ,2 (θ, ϕ, x) is the second geometry factor associated with γ •,A (θ, ϕ), as will be explicitly defined below.Note that the differentiation in Eq. ( 27) is with respect to ϕ and hence the mapping of the t-coordinate yields no extra factor.
As was commented on in Remark 4, the root ϕ 0 (θ, x) is not defined for x = (0, 0, z), z = a.Furthermore, the geometry factor in E T Z f ac ((0, 0, z), θ), z = a is infinite.To proceed, we will derive the expression for the general E T Z f ac (x, θ), where x is not on the z-axis, and then consider the appropriate limit to understand the behavior for evaluation points along the z-axis.
We have λ > 1, and where the last step holds true for λ > λ 0 as this quantity approaches 1 from above as λ increases.Similarly Using these inequalities in Eq. ( 50), the result follows.
Remark 6 E T Z f ac (x, θ) uses the root ϕ 0 (θ, x), and is hence not defined for ρ 2 = x 2 + y 2 = 0 or θ = {0, π} (see the remark following Lemma 2).In the next theorem we bound the maximum value of E T Z f ac for a spherical surface as the evaluation point approaches the z-axis.
(56) It then follows where the constant C is the same as in Eq. (51).Note that the value of β does not affect the result due to the axisymmetry.
Proof (Theorem 2) The point on the surface of the sphere closest to x(α) is the point with θ = α at the same azimuthal angle (ϕ = β).Differentiating E T Z f ac (x(α), θ) as given in Eq. ( 50) with respect to θ, one can show that θ = α yields the maximum, as one would expect (the formulas do however get very long).
Starting with the formula in Eq. (51), we can evaluate the expression in the first denominator using Eq. ( 56) and θ = α, Corollary 2 Let all be as in Theorem 2, and specifically let the evaluation point x(α) be as in Eq. ( 56) with m = 1.As the evaluation point x approaches the z − axis the trapezoidal rule error E T Z γ (f, p, nϕ, x) in Eq. (48) vanishes, i.e. lim α→{0,π} Proof (Corollary 2) Using the definition in Eq. ( 48), Hence with C = 4π Γ(p) n p−1 E1 f t, ϕ 0 (θ(t), x) dt, it follows From Theorem 2 it follows that lim α→{0,π} E T Z f ac (x(α), α) = 0, and the result follows.To the left in Fig. 1, we plot E T Z f ac (x(α), θ) versus θ/π for several values of α for a sphere of radius 1, with p = 1/2.E T Z f ac (x(α), θ) peaks at θ = α, and decays rapidly away from θ = α.The closer the evaluation point is to the surface (i.e. the closer m is to 1), the larger the maximum magnitude for a fixed number of discretization points n ϕ .For a fixed m, the maximum magnitude decreases rapidly with n ϕ .In the middle figure, we plot E T Z f ac (x(α), α) versus α/π for the combinations n ϕ = 20 and 40 and m = 1.01 and 1, 05.In the rightmost plot, we zoom in to see the behavior for small values of α.As sin α ≈ α for small α, we expect to see a decay proportional to α 2nϕ , and we indicate these slopes in the plot for the two values of n ϕ .From Corollary 2, we have the result that the trapezoidal rule error E T Z γ is bounded by a constant times E T Z f ac (x(α), α), and in the plots we can see the fast decay of E T Z f ac (x(α), α) with decreasing α.Let us now consider the Gauss-Legendre rule error.Write E GL γ Eq. (32) from Error estimate 3 as where This term can be directly evaluated for points on the symmetry axis.For an evaluation point x = (0, 0, z), we find that for a sphere of radius a, where δ = |z| /a if |z| > a and δ = a/|z| if |z| < a.We use the result in Corollary 1 for this derivation.For details, see Section A. This result is independent of ϕ, and hence we have This means that the Gauss-Legendre error E GL will strongly dominate over the trapezoidal rule error E T Z for evaluation points close to the symmetry axis, and the trapezoidal rule error can hence safely be ignored.
For a general axisymmetric surface, we cannot follow the same approach as for the sphere, were we could identify the two parametrization angles for the closest point on the surface, and furthermore show that E T Z f ac (x, θ) attains its maximum value for that value of the polar angle θ.As an evaluation point sufficiently close to a general axisymmetric surface approaches the z-axis, the closest point on the surface will however be the north or the south pole.Hence, it is interesting to investigate this limit.
For general axisymmetric surfaces, we have not been able to prove that the trapezoidal error vanishes as the evaluation point approaches the z-axis, as we have done for the sphere in Corollary 2. From the theorem above, we have the result that E T Z f ac (x, θ) vanishes as x approaches the z-axis for some range of θ, under some assumption on the z-coordinate.We have however not proven for which value of θ the maximum of E T Z f ac (x, θ) is attained.We do conjecture that this θ → {0, π} as the evaluation point approaches the z-axis for z < 0 and z > 0, respectively, as this θ is the parameter for the point on the surface closest to the evaluation point.
Numerically, we note that the trapezoidal error contribution to the total error estimate decays rapidly as the evaluation point approaches the z-axis also for general axisymmetric surfaces.

Simplified error estimate for a spherical surface
In this section, we will present a simplified error estimate for a spherical surface, and then discuss the derivation of it, starting from Error estimate 3.
Introduce ζ = x = x 2 + y 2 + z 2 and an even integer n.The error in approximating the integral with the n t = n/2 point Gauss-Legendre rule in the t-direction and the nϕ = n-point trapezoidal rule in the ϕ direction can be estimated as where Remark 7 Note this error estimate only depends on the evaluation point through ζ = x .This means that the error is estimated to decay equally in all directions as we move away from the sphere.This is only a good approximation under the map t = − cos(θ), and not under the linear map, as will be discussed in section 7.1.
To derive the simplified estimate, we will use the assumption that the error only depends on the distance to the sphere, and will pick an evaluation point that yields the simplest expressions to work with.We will start by considering the trapezoidal rule error at a point x = (x, y, 0); in this case ζ corresponds to the distance from the z-axis that was previously denoted by ρ, so we will continue with this notation, assuming ρ 2 = x 2 + y 2 = a.The expression for E T Z f ac (x, θ) is given in Eq. ( 50), but we will start with the equivalent expression in Eq. ( 53).For this case we obtain λ = (ρ/a + a/ρ)/(2 sin θ) = (δ + 1/δ)/(2 sin θ) where we let δ = ρ/a if ρ > a and δ = a/ρ if ρ < a, such that δ > 1.With this, we have From Theorem 2, we know that E T Z f ac (x, θ) will attain its maximum at θ = π/2 for the chosen evaluation point.At this value of θ, the last term vanishes.Ignoring that term we get λ + √ λ 2 − 1 ≈ δ/ sin θ with equality at θ = π/2.If we use this approximation, and evaluate the part of E T Z f ac that is taken to the power of p at θ = π/2, we obtain Now, under the map t = − cos(θ) an approximation to the trapezoidal rule error will be Under the coordinate transformation t = sin β, the integral can be written as Using that for q > 1, and in total we get We will now continue by estimating the size of the Gauss-Legendre rule error.We write E GL γ from Error estimate 3 as in Eq. ( 57), with E GL f ac (x, ϕ) defined in Eq. ( 58).
Again, picking an evaluation point at the equator, we can derive the approximation where ρ = x .The details of this derivation is given in Section A.
Comparing Eq. (68) (n = n ϕ ) and Eq.(69) (n = 2n t ), the expressions are very similar.The ratio is ).This means that the contributions of the two errors are of about equal size for evaluation points in the xy-plane at z = 0, as opposed to the case of evaluation points at the z-axis, where the contribution for the trapezoidal rule E T Z γ vanishes.The total error is however approximately equal for the same ζ = x .Since there is some overestimation of the errors, we have chosen to estimate the full error as a function of ζ by the derived expression for the trapezoidal rule error at the equator, as is given in (63) in Error estimate 4. The accuracy of this estimate will be numerically evaluated in Section 7.1.

Numerical evaluation of the error estimate
An estimate E EST γ for the quadrature error E Q γ in the evaluation of a layer potential was given in Error estimate 3. The integrals in (31)-(32) adding up to E EST γ can however not be evaluated analytically, and we need to introduce an approximation that is sufficiently precise and computationally cheap to evaluate.For general surfaces, we do not have access to analytical expressions for the roots appearing in the estimate.Hence, they must be computed using a numerical root finding procedure, and it will be of interest to minimize the number of root evaluations.
The integrand in (31) is not well defined for evaluation points along the symmetry axis for an surface, as was discussed in section 5.2.We proved that the contribution from the trapezoidal error (31) vanishes in the limit of the target point approaching the symmetry axis of a sphere (Corollary 2).For a general axisymmetric surface, we were able to prove only a weaker result but we quantified the contribution of the trapezoidal error (31) compared to the contribution of the Gauss-Legendre error (32) by numerical experiments.Also here we find that the first quantity decays very rapidly as the evaluation point approaches the z-axis, depending also on the distance of the evaluation point from the surface.To be more precise, the problematic target points lie in the region where ρ → 0 or/and a(θ) sin(θ) → 0 (see Theorem 3) geometrically represented by the cones with apices at the poles and increasing width inside or outside the surface respectively for the interior or exterior problem.To see why, consider a number of target points placed in the normal direction starting from a grid point γ(t * , ϕ * ); it is clear that they will all refer to the same closest grid point, but ρ (the distance from the z-axis) will increase together with the distance from the surface.In practical applications we will safely ignore the trapezoidal rule contribution for the evaluation points x such that where K c is a fixed constant, ρ the distance to the z-axis and A a representative radius/length scale.This procedure is not very sensitive to the choice of A and K c , as there is a rather wide range where the trapezoidal rule error contribution is negligible but it is still numerically stable to compute the error estimate.Practically, to evaluate the distance between the target point and the surface it is sufficient to approximate the minimum in (70) by the minimum over the surface grid points only.When considering non-axisymmetric surfaces, the situation is much less predictable and we need a different strategy.This is based on using a local approximation of the surface centered away from the poles; in this way all the quantities needed for the estimate evaluation are locally well defined, and the singularities mentioned above are eluded.How to numerically compute such an approximation will be further discussed in the next subsections, where we describe how to practically evaluate the two error contributions (31) and (32), that add up to the total error.We start by the approximation of the integrals before we discuss the root finding.

Approximation of integrals in the error estimate
In this section, we discuss how to approximate the integrals in (31)-(32), to be able to efficiently compute a sufficiently precise estimate E EST γ for the quadrature error E Q γ as defined in Error estimate 3. Given an evaluation point x ∈ R 3 , we start by identifying x * ∈ R 3 , the closest discrete point on the surface γ, and the parameters t * , ϕ * such that γ(t * , ϕ * ) = x * .This means that t * will be of the n t Gauss-Legendre quadrature nodes, and ϕ * will be one of the n ϕ (equidistant) trapezoidal rule quadrature nodes.Loosely speaking, the contribution to the quadrature error will have a peak around x * .What this means is that est T Z (ϕ 0 (t, x), n ϕ , p) in the integral over t in (31) will have a peak close to t = t * , and decay rapidly away from t * due to the variation in ϕ 0 (t, x).Similarly est GL (t 0 (ϕ, x), n t , p) will have a peak close to ϕ = ϕ * , decaying rapidly away from ϕ * due to the variation in t 0 (ϕ, x).Let us here remind that the roots t 0 (ϕ, x) and ϕ 0 (t, x) are roots to R 2 (t, ϕ, x) with one variable kept fixed, as defined in Error estimate 3, and further denote We now assume that est T Z is the most rapidly varying factor in the integrand in (31), and similarly for est GL in (32), and approximate This is typically a good approximation, unless the surface grid is very stretched such that the grid resolutions on the surface in the two directions are very different, in which case the geometry factors can vary rapidly as well.Now, we want to find a simple expression for how t 0 (ϕ, x) varies with ϕ around ϕ = ϕ * .If we replace γ with its bivariate linear approximation around (t * , ϕ * ) in the definition of the squared distance function (20), we get a quadratic equation that we can solve to find the root.From here, we find the approximation where, with ∆ϕ = ϕ − ϕ * , and The root t * 0 is by definition the root at ϕ = ϕ * , and in practice, as will be discussed in the next subsection, at least an accurate approximation to it, while t L 0 (ϕ * ) is a simpler approximation.In order to better capture the magnitude of the peak of est GL (t 0 (ϕ, x), n t , p), we want to this more accurate value, but we also want to use the simple dependence on ϕ.This leads us to define See [1] for more details and a discussion regarding the effect of making these approximations.
The same approximations can naturally be made to define φ0 (t), an approximation to ϕ 0 (t, x).Away from t = t * , we have Im φ0 Given this decay, it is a reasonable approximation to expand the interval of integration in (72 as the tails will be negligible, and we approximate the integral in (72) by With the variable transformation x = n ϕ ks, we can write where h T Z ± (x) = est T Z φ0 (t * ± x/(n ϕ k)) e x .Gauss-Laguerre quadrature is a Gaussian quadrature for integrals of this type [17, S 3 .5(v)],and we find that is is sufficiently accurate to evaluate each of the integrals of h T Z + (x) and h T Z − (x) with 8 quadrature nodes.
For the Gauss-Legendre estimate, we have the bound [10] est Hence, with t 0 (ϕ, x) approximated with t0 (ϕ), we have an estimated decay e −2ntk −1 |ϕ−ϕ * | .Based on this decay we use the variable transformation x = 2n t k −1 s, and write the approximation of the integral in (73) as where h GL ± (x) = est GL t0 (ϕ * ± xk/(2n t )), n t p e x .Again, each of these integrals is approximated with an 8-point Gauss-Laguerre quadrature rule.
In [1], we used this strategy for the global trapezoidal rule discretization, and a different strategy for the panel based Gauss-Legendre quadrature.Here, we have a mix of the two quadrature rules, but both are used globally on the surface, and we have extended this approach to be used for both contributions to the error estimate.It remains now to discuss the root finding, and the evaluation of the factors in front of the integrals in (72)-(73).

Root finding
To evaluate the error estimate as described in the previous subsection, we need to determine t * 0 = t 0 (ϕ * , x) and ϕ * 0 = ϕ 0 (t * , x) as defined in (71).For spherical topologies, these roots are not always well defined, as introduced in 5.2 for axisymmetric geometries (where ϕ 0 (θ, x) is not defined for evaluation points on the symmetry axis) and further discussed at the beginning of section 6 for general surfaces.For axisymmetric geometries we identified the problematic region in the cone defined by eq. ( 70): for these evaluation points we can ignore the trapezoidal error contribution and then we do not need to compute the roots ϕ 0 (θ, x).We will then not consider these points in the following discussion.For a general surface, it is not so easy the determine a similar set, and we will proceed with a discrete approach as later discussed.
In Section 5, we derived analytical expressions for the roots of the squared distance function for special geometries: we have analytical expressions for θ 0 only for a spherical surface, and for ϕ 0 for any axisymmetric surface.Hence, in general we need a numerical procedure to determine the roots, and we will define this procedure in the (θ, ϕ) coordinate system.Given an evaluation point x = (x, y, z), we define Given a parametrization γ • (θ, ϕ), it is easy to solve R 2 (θ, ϕ, x) = 0 using a one dimensional Newton's method to find a root ϕ 0 given θ, or similarly, to find a root θ 0 given ϕ.Specifically, in our setting, we need to determine θ * 0 = θ 0 (ϕ * , x) and ϕ * 0 = ϕ 0 (θ * , x), where θ * = θ(t * ).We typically use an initial guess of θ * +i/10 for θ * 0 , and correspondingly for ϕ * 0 .The iterations then usually converge with a strict tolerance in less than five iterations.However in rare cases, usually for evaluation points far away from the surface, it may happen that the iterations fail to converge, but most of the times it is sufficient to increase the magnitude of the imaginary part of the initial guess for the iterations to converge well.
If no parametrization is available and we know only the quadrature node values of γ • at n t × n ϕ nodes, we need to define an approximation γ• (θ, ϕ) that can be evaluated at different arguments of θ and ϕ, and that allows for differentiation to formulate Newton's method.One way is to compute the spherical harmonics coefficients using a discrete transform [18].Evaluating the spherical harmonics expansion is however a global procedure with O(n t n ϕ ) cost for arbitrary arguments, and this would need to be done at each step in the Newton iteration.
We can however use a local approach.In the Newton iteration, one of the variables, θ or ϕ will be fixed.Assume that ϕ = ϕ * and introduce a qth order Taylor expansion in θ, This approximation can be used in Newton's method to determine an approximation to the root θ * 0 = θ 0 (ϕ * , x).Similarly, a Taylor expansion in ϕ can be introduced to obtain an approximation to ϕ * 0 = ϕ 0 (θ * , x).This is a solid strategy that eludes the pole singularities for any kind of geometry.Indeed, since the discretization is based on Gauss-Legendre nodes in the polar angle, the closest grid point where the expansion (82) is centered, will never be a pole, avoiding the above mentioned problems.
For this reason, we will use this approach for non-axisymmetric surfaces, even if we have a parametrization of γ • available.In this case, the analytical expression for γ • can however be used to determine the q first partial derivatives of γ • .When this is not the case, the derivatives both with respect to θ and φ, can be evaluated e.g from a spherical harmonics expansion.The advantage compared to using a global expansion is that we can evaluate these derivatives at all grid points in one sweep [4,19], and then use different local expansions when estimating the quadrature error for different evaluation points.
Once the roots have been determined, these Taylor expansions can also be used to evaluate the geometry factors in (72)-(73) as defined in (26)-( 27).The roots t * 0 and ϕ * 0 are also needed to evaluate f in (72)-(73).We recall that f depends on the density σ, which may not be known analytically but be available at the grid points only (e.g. if σ is a solution to a discretized integral equation).In this case, again, f can either be approximated locally by a Taylor expansion, or globally by a spherical harmonics expansion.

Numerical experiments
In this section, we will compare the quadrature error estimate E EST γ defined (30), and evaluated as discussed in the previous section, to the measured error E Q γ as defined in (29) for some different examples.The measured error E Q γ will be computed by using a reference solution on an upsampled grid with upsampling rate set to five.
As outlined in Remark 1, we will see that the estimate provides a good approximation of the error also for moderate values of n t and n ϕ as long as the geometry and the layer density are well resolved.We choose n t and n ϕ so that this is true for all the following numerical examples.For the rootfinding procedure, we follow the strategy presented in the previous section: the analytical expression for γ• is used when dealing with axisymmetric surfaces (excluding the trapezoidal error contribution for target points defined in (70)), and γ • is locally approximated with a Taylor expansion for other geometries.In the first case, the parameters defining the cone in (70) will be kept fixed to A = 1 and K C = 10.In the latter case, the order of the Taylor expansion will be fixed as q = 4 (see eq. ( 82)), which is accurate enough for our purposes; an analysis of how the choice of q can affect the accuracy of the roots can be found in [1].In all the presented numerical tests we will use the analytical expression for the density.

A sphere
In the first example we consider the harmonic single layer potential evaluated near a sphere of radius a = 1 with unit density, σ(x) = 1.We want to compare the estimated error E EST γ to the actual measured error E Q γ , using the full error estimate in (30), approximated as described in the previous section, and, for the cosine map, also the simplified error estimate (4), derived in Section 5.3.For this case, referring to eq. ( 30), p = 1/2 and f simplifies to 2 ), if using the linear mapping (9) 1, if using the cosine mapping (10).Fig. 2 shows the resulting surface grids and the different behavior of the quadrature error exterior to the sphere with discretizations using the linear and non-linear mapping θ(t) as given in ( 9)- (10), with n t = 30 and n ϕ = 60 points.The linear map clusters the grid points more towards the poles, and at a fixed distance from the sphere, the error is smaller in these regions, while the cosine mapping yields a more even error.In both cases we see oscillations in the error on a length scale of the grid size.Fig. 2: Error E Q γ (log 10 scale) in computing the harmonic single layer potential with unit density on a at y = 0 cutting a sphere of radius 1 discretized by using (a) the linear mapping, (b) the cosine mapping.In both cases, n t = 30 and n ϕ = 60.
In Fig. 3, the full error estimate E EST γ from Error estimate 3, approximated as described in the previous section, is compared to the measured errors E Q γ .The measured errors (log 10 scale) are shown in one selected plane as colored fields, with the contours of the estimates drawn in black.This is done for eval- uation points both interior and exterior to the sphere, and the error estimates can be seen to work The contours of the estimate are smoothly enclosing the oscillatory error, with an over-estimation that is larger further out exterior to the sphere.Note however that the last contour is at an error level of 10 −14 , which is a very low level.
In Error estimate 4, we derived a simplified error estimate applicable to this case (when using the cosine map) that depends only on the radius a, the distance from the surface, the number of discretization points and the halfinteger p, where p = 1/2 for the integral in (83).In Fig. 4, we compare the error measured in a set of evaluation points to the simplified error estimate, both versus the distance to the surface (here a negative distance d is used for interior point) for a fixed grid resolution, and versus the grid resolution for points at a fixed distance to the sphere.The discrete points are set using the parametrization for a sphere with radius (1 + d), over the full range of the polar angle 0 ≤ θ ≤ π, but only over an angle sector in the azimuthal angle, 0 ≤ ϕ ≤ π/n ϕ .The actual errors do not depend on ϕ, more than that there is an oscillation determined by the grid size, and this range is sufficient to cover the range of errors.Under the cosine-map, the error is much less dependent on θ compared to the linear mapping, but there is a variation, and we include the full range here.From the discrete dots, each representing a different evaluation point, we can see the range of errors for evaluation points at the same distance to the sphere.The simplified error estimate works better than we could expect, and gives a rather tight upper bound of the error.

A prolate spheroid
In the second example we consider an axisymmetric ellipsoid, a prolate spheroid, with ratio 3-1 between the long and short axes.Here the density function is given by σ(θ, ϕ) = 1 + sin(6ϕ + θ) sin 2 (θ). (84) and is in Fig. 5a visualized on the surface by the black and white colormap.We can see how the varying density breaks the geometric symmetry of the problem.We first consider the quadrature error for evaluation points on a vertical wall placed at y = 1.02,Fig. 5a-6a, and then we place random evaluation points around the spheroid (Fig. 5b), and plot the error vs the estimate in Fig.In the first case n t = 40, in the second case n t = 60, and for both we set n ϕ = 2n t .In both cases the estimates can predict very well the actual error.Moreover, it is clear that the density has an effect on the error, but still the simplification made in (72)-( 73) is good enough to capture the behavior of the overall error.

Non-axisymmetric geometry
In the last example we consider a non-axisymmetric geometry given by We use the cosine mapping and define γ(t, ϕ) = γ • (cos −1 (−t), ϕ).
For this case, referring to eq. ( 30), p = 1/2 and f (t, ϕ) = e −ω γ(t,ϕ)−x σ(t, ϕ) ∂ γ ∂t × ∂ γ ∂ϕ .We consider the case ω = 3 and the density function σ given by eq. ( 84).In Fig. 8a we show the error and the estimates computed on the whole spherical shell, plotted here with the horizontal axis being the the azimuthal angle and the vertical axis being the polar angle.In Fig. 8b we zoom in on the white rectangle highlighted in Fig. 8a, to better show the agreement between estimate and error.

Conclusions
In this paper, we studied the error incurred by numerically approximating layer potentials over surfaces of spherical topology.We have derived error estimates for discretizations with the trapezoidal rule in the azimuthal angle, and a Gauss-Legendre rule in a variable that maps to the polar angle.The framework for the derivation of the error estimates, and the practical evaluation there of, were introduced in [1] for surfaces of genus 1, discretized by either a global trapezoidal rule in both directions, or a panel based Gauss-Legendre rule.Here we extended this approach with special attention to the global parametrization.There is one component of the error estimate that cannot be directly evaluated for evaluation points on the symmetry axis of an axisymmetric surface.Starting by deriving analytical expressions for the roots of a squared distance function for a sphere, we were able to prove that this contribution vanishes at these points.We could also derive a simplified error estimate for the sphere, that shows the decay in error with the distance of the evaluation point to the sphere with a simple formula.Some analytical results were also extended to the more general case of an axisymmetric surface, and we devised a strategy for evaluating the error estimate for general surfaces, avoiding the difficulties associated with the discretization around the poles.

Fig. 3 :
Fig. 3: Error E Q γ (colors) and estimates E EST γ (black lines) plotted in log 10 scale when computing the harmonic single layer potential with unit density over a sphere of radius 1 discretized by using (a) the linear mapping, (b) the cosine mapping.The results are shown for evaluation points in the xz-plane for y = 0, both inside and outside of the sphere that is indicated with a dashed white line.

Fig. 4 :
Fig. 4: Measured quadrature errors E Q γ in evaluating (83) (using the cosine map) at a set of discrete evaluation points x ∈ R 3 (blue dots) and the simplified estimate (red lines); (a) varying distance to the sphere (negative values for interior) and fixed n t = 30 and (b) varying n t at the fixed distance d = 0.1.Note that n = n ϕ = 2n t in the simplified estimate.

Fig. 5 :
Fig.5: Target points considered when computing the harmonic single (a) and the double (b) layer potentials evaluated near a prolate spheroid.The and white colormap represents the density given by eq.(84).The red and blue colormap represents the actual error evaluated on the target wall.

Fig. 6 :
Fig. 6: (a) Error E Q γ (colors) and estimates E EST γ (black lines) plotted in log 10 scale when computing the harmonic single layer potential near the ellipsoid on a plane at y = 1.02, shown in Fig. 5a.(b) Estimates vs error in computing the harmonic double layer potential at random evaluation points showed in Fig. 5b.The three lines from top to bottom indicate where the estimate of E EST γ is a factor of 10, 1 and 1/10 times the measured value of E Q γ .If the estimate was never underestimating the error, no dot would fall below the red line.

Fig. 7 :
Fig.7: Half of the spherical shell enclosing the non axisymmetric shape defined in (86).The black and white colormap represents the density given by eq.(84).