Fast Galerkin Method for Solving Helmholtz Boundary Integral Equations on
Screens ††thanks: This work was partially funded by Fondecyt Regular 1171491
and grant Conicyt 21171479 Doctorado Nacional 2017.Fast Galerkin Method for
Solving Helmholtz BIEs on Screens

# Fast Galerkin Method for Solving Helmholtz Boundary Integral Equations on
Screens ††thanks: This work was partially funded by Fondecyt Regular 1171491
and grant Conicyt 21171479 Doctorado Nacional 2017.

Carlos Jerez-Hanckes, José Pinto  Facultad de Ingeniería y Ciencias,
Universidad Adolfo Ibáñez, Santiago, Chile. Email<EMAIL_ADDRESS>of Engineering, Pontificia Universidad Católica de Chile, Santiago, Chile,
Email<EMAIL_ADDRESS>
###### Abstract

We solve first-kind Fredholm boundary integral equations arising from
Helmholtz and Laplace problems on bounded, smooth screens in three-dimensions
with either Dirichlet or Neumann conditions. The proposed Galerkin-Bubnov
method takes as discretization elements pushed-forward weighted azimuthal
projections of standard spherical harmonics onto the canonical disk. We show
that these bases allow for spectral convergence and provide fully discrete
error analysis. Numerical experiments support our claims, with results
comparable to Nyström-type and $hp$-methods. Boundary integral equations,
spectral methods, wave scattering problems, screens problems, non-Lipschitz
domains.

##  1 Introduction

We study the solution of the Helmoltz and Laplace problems with Dirichlet or
Neumann conditions posed on an open orientable bounded surface
$\Gamma\subset\mathbb{R}^{3}$. These can be summarized as follows:

######  Problem 1.

Let $k\geq 0$, find $u$ defined in $\mathbb{R}^{3}$, such that

$\displaystyle-\Delta u-k^{2}u=0,\quad\text{on
}\mathbb{R}^{3}\setminus\overline{\Gamma},$ $\displaystyle
u=g_{d}\quad\text{or}\quad\partial_{n}u=g_{n}\quad\text{in }\Gamma,$
$\displaystyle\text{condition at infinity}(k),$

where $\partial_{n}u$ is the normal derivative of $u$ with the normal defined
via the parametrization used111A fixed orientation for the normal vector is
required, and one can choose either of the two orientations that yield
continuous normal vectors without loss of generality., $g_{d}$ and $g_{n}$ are
suitable Dirichlet and Neumann data, respectively, and the condition at
infinity reads

$\displaystyle\begin{cases}\lim_{\|\mathbf{x}\|\rightarrow\infty}\|\mathbf{x}\|\left(\partial_{\|\mathbf{x}\|}u-iku\right)=0,&k>0,\\\
u(\mathbf{x})=O(\|\mathbf{x}\|^{-1}),\text{ as
}\|\mathbf{x}\|\rightarrow\infty,&k=0,\end{cases}$

For a detailed discussion of these conditions see [25, Chap. 8 and 9]. The
Laplace case occurs when $k=0$.

Problem 1 and associated boundary integral equations (BIEs) have been
extensively studied [34, 33, 9, 13, 16, 22, 27, 17]. Indeed, BIEs rigorously
recast the volume problem onto the screen while taking into account the
corresponding condition at infinity. By meshing the open surface $\Gamma$,
standard local low-order approximations can be built and solved by today’s
compression algorithms [4, 40] to accelerate computations. However, solving
Problem 1 remains far from trivial as standard boundary element
implementations converge at the worst rate possible due to the solution’s
singular behavior near the screen surface. Thus, one has to resort to special
meshing techniques [14] or to increase the polynomial degree of the underlying
discretization in a suitable manner [15, 5] so as to recover better error
convergence rates. This however is a tedious procedure when solving for
multiple configurations as in uncertainty quantification or inverse problems.

In two-dimensional space, i.e. when $\Gamma$ is an open arc, an alternative
discretization of BIEs can be performed by means of weighted Chebyshev
polynomials –spectral discretization–, which explicitly capture the edge
singularity and allow for super-algebraic convergence whenever $g_{d}$ and
$g_{n}$ are smooth functions [29, 18, 19, 20]. Opposingly, for screens in
$\mathbb{R}^{3}$, to the best of our knowledge, no such equivalent
discretization exists. Moreover, though for closed surfaces an spectral BIE
method was presented by Sloan _et al._ [12], the presence of extra
singularities around the surface edges and the structure of the underlying
fundamental solution renders impossible the extension of this method to the
context of open surfaces directly.

Our main contribution is the derivation and analysis of spectral Galerkin-
Bubnov methods to tackle Laplace and Helmholtz BIEs on open surfaces in three-
dimensional space. In contrast to the Nystr̈om approach [6], our spectral
Galerkin method preserves all the theoretical properties of a Galerkin method
and its implementation relies in suitable quadrature rules and change of
variables without the need of special window functions, while maintaining the
convergence rate for smooth geometries. On the other hand, while special $p$\-
and $hp$-discretizations (see [14, 15, 5]) could also achieve super-algebraic
convergence rates, and being more flexible for non-smooths geometries, they
need a mesh which in practice could slowdown the implementation when solutions
for different geometries are needed as in uncertainty quantification or
inverse problem [21].

Our work is structured as follows. In Section 2.3 we define a family of
finite-dimensional functional spaces by projecting the spherical harmonics
into the unitary disk and consider their linear span. We then define the trial
and test spaces for the Galerkin method by mapping the functions to the screen
in consideration. From the definition of our trial space, we construct
adequate auxiliary functional spaces which allow for the associated error
analysis, and we describe how theses spaces are related to the standard
Sobolev one. This fact, in conjunction with standard Galerkin properties for
the variational versions of the BIEs (Section 3.1), lead to a semi-discrete a
priori error bound in Theorem 3. Then, in Section 5.1, we detail the
quadrature procedure used to approximate the action of the weakly- and hyper-
singular boundary integral operators (BIOs) on the trial space, and provide a
full error analysis which incorporate the quadrature error. Finally, in
Section 6 we show some numerical examples that exhibit the super-algebraic
convergence method.

##  2 Mathematical Tools

This section introduces general definitions needed for the analysis of Problem
1 and of the proposed Galerkin spectral method.

Given $\mathbf{x},\mathbf{x}^{\prime}\in\mathbb{C}^{d}$, $d\in\mathbb{N}$, the
standard dot product is denoted by
$\mathbf{x}\cdot\mathbf{x}^{\prime}=\sum_{j=1}^{d}x_{j}\overline{x^{\prime}}_{j}$,
and the Euclidean norm satisfies
$\|\mathbf{x}\|^{2}=\mathbf{x}\cdot\mathbf{x}$. For $d=3$ and real vectors, we
denote $\mathbf{x}\times\mathbf{x}^{\prime}$ the cross-product –vectorial
product. For $a,b\in\mathbb{R}$ we will make use of the notation $a\lesssim b$
whenever there is a positive number $C$ such that $a\leq Cb$, typically $C$ is
a constant not relevant for the underlying analysis. Also, we will denote by
$L^{2}(A)$ the classical Lebesgue space of square-integrable functions over a
measurable set $A\subset\mathbb{R}^{d}$, $d=2,3$.

###  2.1 Geometry

Throughout, we denote $\Gamma\subset\mathbb{R}^{3}$ a smooth orientable
connected surface with boundary, also called smooth screen or simply screen.
We assume that the screen is contained on a smooth orientable closed surface
which will be denoted $\widetilde{\Gamma}$. The canonical disk and spheres are
given, respectively, by

$\displaystyle\mathbb{D}:=\\{\mathbf{x}\in\mathbb{R}^{2}:\|\mathbf{x}\|\leq
1\\},\quad\mathbb{S}:=\\{\mathbf{x}\in\mathbb{R}^{3}:\|\mathbf{x}\|=1\\}.$

For any given screen, we assume that there is a parametrization
$\mathbf{r}:\mathbb{D}\rightarrow\Gamma$, being a smooth function such that in
every argument and coordinate can be extended to an analytic function on a
Bernstein ellipse in the complex plane of parameter222Given and interval
$[a,b]\subset\mathbb{R}$, the Bernstein ellipse of parameter $\rho$,
corresponds to an ellipse on the complex plane with foci at $a,b$, semi-major
axis $\frac{(b-a)(\rho+\rho^{-1})}{4}$, and semi-minor axis
$\frac{(b-a)(\rho-\rho^{-1})}{4}$.  $\rho>1$. The class of functions
portraying the regularity previously described are referred to as
$\rho-$analytic. Furthermore, we impose that the gradients of $\mathbf{r}$ (as
a matrix) have full rank for any point, and the Jacobian,
$J_{\mathbf{r}}(\mathbf{x})=\|\partial_{x_{1}}\mathbf{r}\times\partial_{x_{2}}\mathbf{r}\|$
is such that $J_{\mathbf{r}}(\mathbf{x})/\|\mathbf{x}\|$ is bounded and
nowhere null.

The direction of the unitary normal vector of $\Gamma$ is selected to be equal
to $\partial_{x_{1}}\mathbf{r}\times\partial_{x_{2}}\mathbf{r}$. This imply
that, rigorously speaking, a screen is characterized by the particular
parametrization $\mathbf{r}$ and no by the set of points that represent the
physical domain $\Gamma$.

On $\mathbb{S}$, we will often consider spherical coordinates characterized by
two angles $(\theta,\varphi)\in[0,2\pi]\times[0,\pi]$ . We will always assume
that $\theta$ corresponds to the polar angle of a given point when projected
to the $x_{3}=0$ plane, and $\varphi$ denotes the angle measured from the
$x_{3}-$axis. Similarly, in $\mathbb{D}$ we use polar coordinates
$(r,\theta)\in[0,1]\times[-\frac{\pi}{2},\frac{3\pi}{2}]$.

###  2.2 Classical Functional Spaces

Given a set $\mathcal{O}\subset\mathbb{R}^{d}$, $d\in\\{2,3\\}$, we denote by
$\mathcal{D}(\mathcal{O})$ the space of smooth functions with compact support
in $\mathcal{O}$ with topological dual $\mathcal{D}(\mathcal{O})^{\prime}$. If
$\mathcal{O}$ is open then $\mathcal{D}(\mathcal{O})$ are smooth functions
whose extension by zero is also smooth. If $\mathcal{O}$ is compact then
$\mathcal{D}(\mathcal{O})$ are the set of infinity differentiable functions.
We will require Sobolev spaces defined on open and closed surfaces. For
$\Gamma$ an open surface, the Sobolev spaces $H^{s}(\widetilde{\Gamma})$,
$s\in\mathbb{R}$, are defined by using local parametrizations and a partition
of unity [25, Chap. 2]. Spaces $H^{s}(\Gamma)$, and
$\widetilde{H}^{s}(\Gamma)$ , with $s\in\mathbb{R}$, are also standard and can
be defined as

$\displaystyle H^{s}(\Gamma)$
$\displaystyle:=\\{u\in\mathcal{D}^{\prime}(\Gamma):\exists\ U\in
H^{s}(\widetilde{\Gamma}),\ u=U|_{\Gamma}\\},$
$\displaystyle\widetilde{H}^{s}(\Gamma)$ $\displaystyle:=\\{u\in
H^{s}(\widetilde{\Gamma}):\ \mathrm{supp}(u)\subset\overline{\Gamma}\\},$

where the support and restriction have to be understood in the context of
distributions. One can identify the dual space of $H^{s}(\Gamma)$ with
$\widetilde{H}^{-s}(\Gamma)$, the corresponding duality product is denoted
$\langle\cdot,\cdot\rangle_{\Gamma}$. The duality product is an extension of
the $L^{2}(\Gamma)$ inner product as we work under the identification,
$H^{0}(\Gamma)=\widetilde{H}^{0}(\Gamma)=L^{2}(\Gamma)$. The following Lemma
is useful to retain control over the norms of functions defined on an
arbitrary screen.

######  Lemma 1 (Theorem 3.23 in [25]).

Let $s\in\mathbb{R}$, if $u\in H^{s}(\Gamma)$, we can define an equivalent
norm as

$\displaystyle\|u\circ\mathbf{r}\|_{H^{s}(\mathbb{D})}\cong\|u\|_{H^{s}(\Gamma)},$

where the unspecified constant depends on $\Gamma$ only. The same result holds
true when we change $H^{s}(\Gamma)$ for $\widetilde{H}^{s}({\Gamma})$ and
$H^{s}(\mathbb{D})$ for $\widetilde{H}^{s}(\mathbb{D})$, accordingly.

######  Remark 1.

Lemma 1 can be generalized to screens of restricted regularity by limiting the
range of $s$ depending on the regularity (cf. [25, Theorem 3.2.3] for
details).

Customarily, we extend the definitions of restrictions and normal derivatives
over $\Gamma$ to linear bounded maps in appropriate Sobolev spaces. In
particular, following [25, Chapter 2], we define the Dirichlet traces as the
maps
$\gamma_{d}^{\pm}:H^{s}(\mathbb{R}^{3}\setminus\overline{\Gamma})\rightarrow
H^{s-\frac{1}{2}}(\Gamma)$ which extensions of the following operator:

$\gamma_{d}^{\pm}u(\mathbf{x})=\lim_{\epsilon\downarrow
0}u(\mathbf{x}\pm\epsilon\widehat{\mathbf{n}}(\mathbf{x})),\quad\mathbf{x}\in\Gamma,$

for $s>\frac{1}{2}$, and where $\widehat{\mathbf{n}}$ denotes the unitary
normal vector whose direction depends on the parametrization $\mathbf{r}$ of
$\Gamma$. The Neumann trace $\gamma_{n}$ is defined as the extension of the
normal derivative:

$\gamma_{n}^{\pm}u(\mathbf{x})=\lim_{\epsilon\downarrow 0}\nabla
u(\mathbf{x}\pm\epsilon\widehat{\mathbf{n}}(\mathbf{x}))\cdot\widehat{\mathbf{n}}(\mathbf{x}),\quad\mathbf{x}\in\Gamma,$

where the extension can be done from
$H^{s}(\mathbb{R}^{3}\setminus\overline{\Gamma})\rightarrow
H^{s-\frac{1}{2}}(\Gamma)$ for $s>\frac{3}{2}$, or from
$H^{1}(\mathbb{R}^{3}\setminus\overline{\Gamma})\rightarrow
H^{-\frac{1}{2}}(\Gamma)$ for functions such that $\Delta u$ is locally in
$L^{2}(\mathbb{R}^{3}\setminus\Gamma)$ [25, Chapter 4]. Whenever
$\gamma^{+}_{d}=\gamma^{-}_{d}$ (resp. $\gamma^{+}_{n}=\gamma^{-}_{n}$) we
denote $\gamma_{d}=\gamma_{d}^{\pm}$ (resp. $\gamma_{n}=\gamma_{n}^{\pm}$).

###  2.3 Spherical Harmonics and Projected Basis

The standard spherical harmonics functions consist of a family of smooth
functions defined in $\mathbb{S}$ given by

$\displaystyle
Y^{l}_{m}(\theta,\varphi):=c_{m}\sqrt{\frac{(2l+1)(l-|m|)!}{4\pi(l+|m|)!}}\mathrm{P}^{l}_{|m|}(\cos\varphi)e^{im\theta},$

where $c_{m}=(-1)^{m}$ if $m<0$ or $c_{m}=1$ otherwise, with
$(\theta,\varphi)\in[0,2\pi]\times[0,\pi]$ are the spherical coordinates on
$\mathbb{S}$, $\mathrm{P}^{l}_{m}$ denotes the associated Legendre function,
and the indices $l\in\mathbb{N}$, $m\in\mathbb{Z}:|m|\leq l$. We will write
$Y^{l}_{m}(\mathbf{x})$ to denote the corresponding spherical harmonic
evaluated at a point $\mathbf{x}\in\mathbb{S}$. The spherical harmonics form
an orthonormal basis of $L^{2}(\mathbb{S})$ (_cf._[24] for more details).

Functions $Y^{l}_{m}$ can be projected onto $\mathbb{D}$, this fact being key
ingredient of our approximation basis. Let us start by defining two lifting
operators on $\mathcal{D}(\mathbb{D})$, the first is the _even lifting_ :

$\displaystyle L_{e}$
$\displaystyle:\mathcal{D}(\mathbb{D})\rightarrow\mathcal{D}(\mathbb{S}),$
$\displaystyle L_{e}(f)(\mathbf{x})$
$\displaystyle=f(x_{1},x_{2}),\quad\mathbf{x}\in\mathbb{S}\subset\mathbb{R}^{3},$

using polar-spherical coordinates we see that
$L_{e}(g)(\theta,\varphi)=g(\sin\varphi,\theta)$, where the arguments of $g$
are the distance to the origin and the polar angle on the plane $x_{3}=0$. The
odd lifting has to be defined on $\mathcal{D}_{0}(\mathbb{D})$ the space of
smooth functions on $\mathbb{D}$ which vanish on the unitary circle in the
plane $x_{3}=0$.

$\displaystyle L_{o}$
$\displaystyle:\mathcal{D}_{0}(\mathbb{D})\rightarrow\mathcal{D}(\mathbb{S}),$
$\displaystyle L_{o}(f)(\mathbf{x})$
$\displaystyle=\begin{cases}f(x_{1},x_{2})\quad x_{3}\geq 0\\\
-f(x_{1},x_{2})\quad x_{3}\leq
0\end{cases},\quad\mathbf{x}\in\mathbb{S}\subset\mathbb{R}^{3},$

Notice that every function in the image of $L_{e}$ (resp. $L_{o}$) is a even
(resp. odd) function on the $x_{3}$ variable. In particular, the spherical
harmonic $Y^{l}_{m}$ is even (resp. odd) when $m+l$ is even (resp. odd). The
definition of the lifting operators can be extended to distributions by
duality. In parallel, one can define a projection operator as

$\displaystyle\Pi_{\mathbb{S}}$
$\displaystyle:\mathcal{D}(\mathbb{S})\rightarrow\mathcal{D}(\mathbb{D}),$
$\displaystyle\Pi_{\mathbb{S}}(f)(\mathbf{x})$
$\displaystyle:=f\left(\mathbf{x},\sqrt{1-x_{1}^{2}-x_{2}^{2}}\right)\quad\mathbf{x}\in\mathbb{D}\subset\mathbb{R}^{2},$

which is the inverse $L_{e}$ (resp. $L_{o}$) when restricted to the functions
with the corresponding symmetries on $\mathbb{S}$. Finally, following [39,
28], the projected basis are defined as

$\displaystyle
p^{l}_{m}(\mathbf{x}):=\sqrt{2}\Pi_{\mathbb{S}}(Y^{l}_{m})(\mathbf{x}),\quad
q^{l}_{m}(\mathbf{x}):=\frac{p^{l}_{m}(\mathbf{x})}{\sqrt{1-\|\mathbf{x}\|^{2}}},$

where $\mathbf{x}\in\mathbb{D}$. From the orthogonality property of spherical
harmonics, it holds that

$\displaystyle\int_{\mathbb{D}}p^{l}_{m}(\mathbf{x})\overline{q^{l^{\prime}}_{m^{\prime}}(\mathbf{x})}d\mathbf{x}=\delta_{m,m^{\prime}}\delta_{l,l^{\prime}}.$
(1)

Also, by symmetry of spherical harmonics, one can see that if $m+l$ is even,
then $p_{m}^{l}$ is a smooth function, while, if $m+l$ is odd, then
$q_{m}^{l}$ is smooth.

An explicit formula for the projected basis is given by:

$\displaystyle
p^{l}_{m}(r,\theta)=C^{l}_{m}\mathrm{P}^{l}_{|m|}(\sqrt{1-r^{2}})e^{im\theta},$
(2)

where $C^{l}_{m}=\sqrt{\frac{(2l+1)(l-|m|)!}{2\pi(l+|m|)!}}$ if $m\geq 0$ and
$C^{l}_{m}=(-1)^{m}C^{l}_{-m}$ if $m<0$.

To simplify notations, we sporadically use only one sub-index for the
projected basis. For $m+l$ even, we can reorder the basis with a one-
dimensional index defined as

$I_{e}(l,m):=\dfrac{l(l+1)+(l+m)}{2}.$

The even function indexed in this way will be denoted $p^{e}_{I_{e}(l,m)}$
(resp. $q^{e}_{I_{e}(l,m)}$) whereas for $m+l$ odd, we set

$I_{o}(l,m)=\frac{(l-1)l+m+l-1}{2},$

with functions denoted $p^{o}_{I_{o}(l,m)}$ (resp. $q^{o}_{I_{o}(l,m)}$). For
example, this leads to

$\displaystyle p^{e}_{0}$ $\displaystyle=p^{0}_{0},$ $\displaystyle\
p^{e}_{1}$ $\displaystyle=p^{1}_{-1},$ $\displaystyle\ p^{e}_{2}$
$\displaystyle=p^{1}_{1},$ $\displaystyle p^{e}_{3}$
$\displaystyle=p^{2}_{-2},$ $\displaystyle q^{o}_{0}$
$\displaystyle=q^{1}_{0}$ $\displaystyle\ q^{o}_{1}$
$\displaystyle=q^{2}_{-1}$ $\displaystyle\ q^{o}_{2}$
$\displaystyle=q^{2}_{1}$ $\displaystyle\ q^{o}_{3}$
$\displaystyle=q^{3}_{-2}.$

######  Lemma 2 (Proposition 2.1.20 in [37]).

The following inclusions are dense in the corresponding Sobolev spaces:

$\displaystyle\operatorname{span}\\{p_{l}^{e}\\}_{l\in\mathbb{N}}\subset{H}^{\frac{1}{2}}(\mathbb{D}),\quad$
$\displaystyle\operatorname{span}\\{p_{l}^{o}\\}_{l\in\mathbb{N}}\subset\widetilde{H}^{\frac{1}{2}}(\mathbb{D}),$
$\displaystyle\operatorname{span}\\{q_{l}^{e}\\}_{l\in\mathbb{N}}\subset\widetilde{H}^{-\frac{1}{2}}(\mathbb{D}),\quad$
$\displaystyle\operatorname{span}\\{q_{l}^{o}\\}_{l\in\mathbb{N}}\subset{H}^{-\frac{1}{2}}(\mathbb{D}).$

###  2.4 Auxiliary Functional Spaces

Classically, Sobolev spaces on $\mathbb{S}$ can be defined in terms of
functions expressed as an expansion of spherical harmonics in the following
way (_cf._ [26] or [3, Chapter 7]):

$\displaystyle
H^{s}(\mathbb{S}):=\left\\{u\in\mathcal{D}^{\prime}(\mathbb{S}):\|u\|_{H^{s}(\mathbb{S})}^{2}:=\sum_{l=0}^{\infty}(l+1)^{2s}\sum_{m=-l}^{l}|\langle
u,Y^{l}_{m}\rangle|^{2}<\infty\right\\},\quad s\in\mathbb{R},$

where $\langle u,Y^{l}_{m}\rangle$ is the extension by duality of the
$L^{2}(\mathbb{S})$-inner product. From this definition, given $u\in
H^{s}(\mathbb{S})$ and $N\in\mathbb{N}$, we can define a finite-dimensional
projection of $u$ as

$\displaystyle\Pi_{N}u=\sum_{l=0}^{N}\sum_{m=-l}^{l}\langle u,Y^{l}_{m}\rangle
Y^{l}_{m}.$

and using elementary properties we can bound the error of this projection as:

$\displaystyle\|u-u_{N}\|_{H^{s}(\mathbb{S})}\leq(N+1)^{s-r}\|u\|_{H^{r}(\mathbb{S})},$

for any reals $s,r$ such that $s<r$. The space $H^{s}_{e}(\mathbb{S})$ (resp.
$H^{s}_{o}(\mathbb{S})$) is defined as the completion of the even (resp. odd)
functions in the $x_{3}$ variable in $\mathcal{D}(\mathbb{S})$, with the same
norm. For these cases, the norms can be written as

$\displaystyle\|u\|_{H^{s}_{e}(\mathbb{S})}^{2}=\sum_{l=0}^{\infty}(n+1)^{2s}\sum_{\begin{subarray}{c}m=-l\\\
m+l\text{ even}\end{subarray}}^{l}|\langle u,Y^{l}_{m}\rangle|^{2},$
$\displaystyle\|u\|_{H^{s}_{o}(\mathbb{S})}^{2}=\sum_{l=0}^{\infty}(n+1)^{2s}\sum_{\begin{subarray}{c}m=-l\\\
m+l\text{ odd}\end{subarray}}^{l}|\langle u,Y^{l}_{m}\rangle|^{2}.$

We will also consider the special function
$w_{\mathbb{S}}:=\sqrt{1-x_{1}^{2}-x_{2}^{2}}$ defined in $\mathbb{S}$. In
spherical coordinates can be written as $w_{\mathbb{S}}=|\cos\varphi|$.

######  Remark 2.

As the function $w_{\mathbb{S}}$ do not depend on the $x_{3}$ coordinate, it
holds that

$\displaystyle w_{\mathbb{S}}L_{e}(u)=L_{e}(w_{\mathbb{D}}u),\quad
w_{\mathbb{S}}L_{o}(u)=L_{o}(w_{\mathbb{D}}u),$

where $w_{\mathbb{D}}(\mathbf{x})=\sqrt{1-\|\mathbf{x}\|^{2}}$, for
$\mathbf{x}\in\mathbb{D}$.

Let us introduce two families of auxiliary spaces on $\mathbb{D}$ associated
with even functions:

$\displaystyle P^{s}_{e}(\mathbb{D})$
$\displaystyle:=\\{u\in\mathcal{D}^{\prime}(\mathbb{D}):L_{e}(u)\in
H^{s}_{e}(\mathbb{S})\\}$ $\displaystyle Q^{s}_{e}(\mathbb{D})$
$\displaystyle:=\\{u\in\mathcal{D}^{\prime}(\mathbb{D}):w_{\mathbb{S}}L_{e}(u)\in
H^{s}_{e}(\mathbb{S})\\},$

with norms given by:

$\displaystyle\|u\|_{P^{s}_{e}(\mathbb{D})}$
$\displaystyle:=\|L_{e}(u)\|_{H^{s}_{e}(\mathbb{S})}=\left(\sum_{l=0}^{\infty}(n+1)^{2s}\sum_{\begin{subarray}{c}m=-l\\\
m+l\text{ even}\end{subarray}}^{l}|\langle
u,q^{l}_{m}\rangle_{\mathbb{D}}|^{2}\right)^{\frac{1}{2}},$
$\displaystyle\|u\|_{Q^{s}_{e}(\mathbb{D})}$
$\displaystyle:=\|w_{\mathbb{S}}L_{e}(u)\|_{H^{s}_{e}(\mathbb{S})}=\left(\sum_{l=0}^{\infty}(n+1)^{2s}\sum_{\begin{subarray}{c}m=-l\\\
m+l\text{ even}\end{subarray}}^{l}|\langle
u,p^{l}_{m}\rangle_{\mathbb{D}}|^{2}\right)^{\frac{1}{2}}.$

Odd function spaces $P^{s}_{o}(\mathbb{D})$, $Q^{s}_{o}(\mathbb{D})$ are
defined in a similar fashion. While the connection with standard spaces is not
as direct as the definition suggests, the next Lemma will be useful.333The
proof is given in Appendix A.

######  Lemma 3.

The following relations between auxiliary and classical spaces on $\mathbb{D}$
holds:

$\displaystyle P^{0}_{e}(\mathbb{D})$
$\displaystyle=L^{2}_{1/w_{\mathbb{D}}}(\mathbb{D}):=\left\\{u:\int_{\mathbb{D}}\frac{u\overline{u}}{w_{\mathbb{D}}}d\mathbf{x}<\infty\right\\}$
(3) $\displaystyle Q^{-\frac{1}{4}}_{e}(\mathbb{D})$
$\displaystyle\subset\widetilde{H}^{-\frac{1}{2}}(\mathbb{D})\subset
Q^{-\frac{1}{2}}_{e}(\mathbb{D})$ (4) $\displaystyle
P^{\frac{1}{2}}_{o}(\mathbb{D})$
$\displaystyle\subset\widetilde{H}^{\frac{1}{2}}(\mathbb{D})\subset
P^{\frac{1}{4}}_{o}(\mathbb{D}).$ (5)

with continuous inclusions.

######  Corollary 1.

The following inclusions are continuous:

$\displaystyle Q_{o}^{-\frac{1}{4}}(\mathbb{D})$ $\displaystyle\subset
H^{-\frac{1}{2}}(\mathbb{D})\subset Q_{o}^{-\frac{1}{2}}(\mathbb{D})$ (6)
$\displaystyle P_{e}^{\frac{1}{2}}(\mathbb{D})$ $\displaystyle\subset
H^{\frac{1}{2}}(\mathbb{D})\subset P_{e}^{\frac{1}{4}}(\mathbb{D}).$ (7)

###### Proof.

Both results are immediate consequences of the duality relation between
classical Sobolev spaces and Lemma 3. ∎

By definition of the auxiliary spaces as well the projector $\Pi_{N}$, we can
introduce four projectors onto the disk. Namely,

$\displaystyle\Pi^{e,p}_{N}:\sum_{l=0}^{\infty}u_{l}p^{e}_{l}\mapsto\sum_{l=0}^{I(N)}u_{l}p^{e}_{l},\quad$
$\displaystyle\Pi^{o,p}_{N}:\sum_{l=0}^{\infty}u_{l}p^{o}_{l}\mapsto\sum_{l=0}^{I(N)}u_{l}p^{o}_{l},$
$\displaystyle\Pi^{e,q}_{N}:\sum_{l=0}^{\infty}u_{l}q^{e}_{l}\mapsto\sum_{l=0}^{I(N)}u_{l}q^{e}_{l},\quad$
$\displaystyle\Pi^{o,q}_{N}:\sum_{l=0}^{\infty}u_{l}q^{o}_{l}\mapsto\sum_{l=0}^{I(N)}u_{l}q^{o}_{l},$

where $I(N)=I_{e}(N,N)=\frac{N(N+1)}{2}+N$. Using the norm definitions on the
corresponding spaces, it holds that

$\displaystyle\|\Pi^{e,p}_{N}u-u\|_{P^{s}_{e}(\mathbb{D})}\leq(N+1)^{s-r}\|u\|_{P^{r}_{e}(\mathbb{D})},$
(8)
$\displaystyle\|\Pi^{o,p}_{N}u-u\|_{P^{s}_{o}(\mathbb{D})}\leq(N+1)^{s-r}\|u\|_{P^{r}_{o}(\mathbb{D})},$
(9)
$\displaystyle\|\Pi^{e,q}_{N}u-u\|_{Q^{s}_{e}(\mathbb{D})}\leq(N+1)^{s-r}\|u\|_{Q^{r}_{e}(\mathbb{D})},$
(10)
$\displaystyle\|\Pi^{o,q}_{N}u-u\|_{Q^{s}_{e}(\mathbb{D})}\leq(N+1)^{s-r}\|u\|_{Q^{r}_{o}(\mathbb{D})},$
(11)

for any reals $s,r$ such that $s<r$.

##  3 Boundary Integral Formulation

We now reduce the original problem to the screen.

###  3.1 Boundary Integral Operators

Let us recall the definitions of the single and double layer potentials:

$\displaystyle(\mathcal{S}_{\Gamma}[k]\lambda)(\mathbf{x}):=\int_{\Gamma}G_{k}(\mathbf{x},\mathbf{x}^{\prime})\lambda(\mathbf{x}^{\prime})d\mathbf{x}^{\prime},\quad\mathbf{x}\in\mathbb{R}^{3}\setminus\overline{\Gamma},$
$\displaystyle(\mathcal{D}_{\Gamma}[k]\nu)(\mathbf{x}):=\int_{\Gamma}\gamma_{n,\mathbf{x}^{\prime}}G_{k}(\mathbf{x},\mathbf{x}^{\prime})\nu(\mathbf{x}^{\prime})d\mathbf{x}^{\prime},\quad\mathbf{x}\in\mathbb{R}^{3}\setminus\overline{\Gamma},$

respectively, where

$G_{k}(\mathbf{x},\mathbf{x}^{\prime})=\frac{e^{ik\|\mathbf{x}-\mathbf{x}^{\prime}\|}}{4\pi\|\mathbf{x}-\mathbf{x}^{\prime}\|}$

denotes the Helmholtz fundamental solution and
$\gamma_{n,\mathbf{x}^{\prime}}$ corresponds to the Neumann trace in the
specified variable.

Following the clasical formulation of Dirichlet (resp. Neumann) problems, we
can seek for solutions of the form
$u(\mathbf{x})=(\mathcal{S}_{\Gamma}[k]\lambda)(\mathbf{x})$ (resp.
$u(\mathbf{x})=(\mathcal{D}_{\Gamma}[k]\nu)(\mathbf{x})$), where $\lambda$
(resp. $\nu$) is an unknown density defined on $\Gamma$. By taking traces, one
naturally defines the weakly- and hyper-singular boundary integral operators
(BIOs),

$\displaystyle(V_{\Gamma}[k]\lambda)(\mathbf{x})$
$\displaystyle:=\gamma_{d}(\mathcal{S}_{\Gamma}[k]\lambda)(\mathbf{x})=\int_{\Gamma}G_{k}(\mathbf{x},\mathbf{x}^{\prime})\lambda(\mathbf{x}^{\prime})d\mathbf{x}^{\prime},\quad
x\in\Gamma,$ $\displaystyle(W_{\Gamma}[k]\nu)(\mathbf{x})$
$\displaystyle:=-\gamma_{n}(\mathcal{D}_{\Gamma}[k]\nu)(\mathbf{x})=-\gamma_{n}\int_{\Gamma}\gamma_{n,\mathbf{x}^{\prime}}G_{k}(\mathbf{x},\mathbf{x}^{\prime})\nu(\mathbf{x}^{\prime})d\mathbf{x}^{\prime},\quad
x\in\Gamma,$

where first BIO is defined as a Lebesgue integral, while the second one is
understood as a principal value [25, Chapter 5].

###  3.2 Boundary Integral Equations

With the above definitions, Problem 1 can be reduced to

######  Problem 2 (BIEs).

For $k\geq 0$, find $\lambda,\nu$
$\in\widetilde{H}^{-1/2}(\Gamma)\times\widetilde{H}^{1/2}(\Gamma)$ such that

$\displaystyle V_{\Gamma}[k]\lambda=g_{d},\quad\text{(Dirichlet BIE)},$ (12)
$\displaystyle W_{\Gamma}[k]\nu=g_{n},\quad\text{(Neumann BIE)}.$ (13)

The equivalence between these BIEs and their corresponding original problems
is established in [33].

######  Theorem 1 (Theorem 2.7, Lemma 2.8, in [32], Theorem 2.1.60 in [30]).

For any $k\geq 0$, $g_{d}\in H^{\frac{1}{2}}(\Gamma)$ (resp. $g_{n}\in
H^{-\frac{1}{2}}(\Gamma)$), there exists one solution
$\lambda\in\widetilde{H}^{-\frac{1}{2}}(\Gamma)$ (resp.
$\nu\in\widetilde{H}^{\frac{1}{2}}(\Gamma)$) for Problem 2. Moreover, the
solution operators are bounded:

$\displaystyle\|\lambda\|_{\widetilde{H}^{-\frac{1}{2}}(\Gamma)}\lesssim\|g_{d}\|_{H^{\frac{1}{2}}(\Gamma)},\quad\|\nu\|_{\widetilde{H}^{\frac{1}{2}}(\Gamma)}\lesssim\|g_{n}\|_{H^{-\frac{1}{2}}(\Gamma)},$

with unspecified constants depending on $\Gamma$ and $k$.

######  Theorem 2 (Corollary A.4 in [10]).

Given $g_{d}$, $g_{n}$ smooth functions on $\Gamma$, the pulled-back solutions
multiplied by the corresponding weight functions,
$(\lambda\circ\mathbf{r})w_{\mathbb{D}}$ and
$(\nu\circ\mathbf{r})(w_{\mathbb{D}})^{-1}$, are as smooth as
$g_{d}\circ\mathbf{r}$, and $g_{n}\circ\mathbf{r}$, respectively.

######  Remark 3.

The last theorem is far from trivial. In the two-dimensional case –open arcs–
solutions are also singular at the arc endpoints, which is proven
straightforwardly for any $k$ and arc as a perturbation of the simpler case
$k=0$, $\Gamma=(-1,1)\times\\{0\\}$. On the other hand, for screens in three-
dimensional space a careful analysis of the BIO symbols is needed (_cf._ [10]
and references therein).

##  4 Spectral Discretizations

From the definition of auxiliary spaces Lemma 2, it is natural to consider a
collection of finite-dimensional spaces spanned by elements $q_{l}^{e}$ (resp.
$p_{l}^{o}$) for the discretization of the Dirichlet (resp. Neumann) BIEs.

For $N\in\mathbb{N}$, we set the following spaces defined over the disk:

$\displaystyle\mathbb{Q}^{N}_{e}(\mathbb{D})$
$\displaystyle:=\operatorname{span}\\{q^{e}_{l}\\}_{l=0}^{I(N)}=\left\\{u:\
u=\sum_{l=0}^{N}\sum_{\begin{subarray}{c}m=-l\\\ m+l\text{
even}\end{subarray}}^{l}u^{l}_{m}q^{l}_{m},\
u^{l}_{m}\in\mathbb{C}\right\\}\subset\widetilde{H}^{-\frac{1}{2}}(\mathbb{D}),$
$\displaystyle\mathbb{P}^{N}_{o}(\mathbb{D})$
$\displaystyle:=\operatorname{span}\\{p^{o}_{l}\\}_{l=0}^{I(N)}=\left\\{u:\
u=\sum_{l=0}^{N+1}\sum_{\begin{subarray}{c}m=-l\\\ m+l\text{
odd}\end{subarray}}^{l}u^{l}_{m}p^{l}_{m},\
u^{l}_{m}\in\mathbb{C}\right\\}\subset\widetilde{H}^{\frac{1}{2}}(\mathbb{D}),$

where again $I(N)=\frac{N(N+1)}{2}+N$. For an arbitrary smooth screen, we
define corresponding spaces through pullbacks as follows

$\displaystyle\mathbb{Q}^{N}_{e}(\Gamma)$ $\displaystyle:=\left\\{u:\
=\frac{(v\circ\mathbf{r}^{-1})\|\mathbf{r}^{-1}\|}{J_{\mathbf{r}}\circ\mathbf{r}^{-1}},\quad
v\in\mathbb{Q}^{N}_{e}(\mathbb{D})\right\\}\subset\widetilde{H}^{-\frac{1}{2}}(\Gamma),$
$\displaystyle\mathbb{P}^{N}_{o}(\Gamma)$ $\displaystyle:=\left\\{u:\
=v\circ\mathbf{r}^{-1},\quad
v\in\mathbb{P}^{N}_{o}(\mathbb{D})\right\\}\subset\widetilde{H}^{\frac{1}{2}}(\Gamma).$

where the inclusions can be easily shown. These spaces are spanned by the
basis functions:

$\displaystyle q^{e,r}_{n}$
$\displaystyle:=\frac{(q^{e}_{n}\circ\mathbf{r}^{-1})\|\mathbf{r}^{-1}\|}{J_{\mathbf{r}}\circ\mathbf{r}^{-1}},\quad
p^{o,r}_{n}:=p^{o}_{n}\circ\mathbf{r}^{-1},$

respectively, and where $n=0,\ldots,I(N)$.

###  4.1 Discrete Problem

Before we introduce discrete versions of Problem 2, we set forth the following
notation. Given $N\in\mathbb{N}$, vectors on $\mathbb{C}^{I(N)+1}$ are written
in bold symbols and superindex $N$. Associated to every vector, there are two
functions which are denoted with the same symbol (not in bold) and superindex
$N$, and an extra the sub-index $e$ if the function is to be understood in
$\mathbb{Q}^{N}_{e}(\Gamma)$, and $o$ if in $\mathbb{P}_{o}^{N}(\Gamma)$. For
example, given $\boldsymbol{\lambda}^{N}\in\mathbb{C}^{I(N)+1}$, one can write

$\displaystyle\lambda^{N}_{e}=\sum_{m=0}^{I(N)}\lambda^{N}_{m}q_{m}^{e,r},\quad\lambda^{N}_{o}=\sum_{m=0}^{I(N)}\lambda^{N}_{m}p_{m}^{o,r}.$

Given $N\in\mathbb{N}$, the Galerkin discretization of the BIE formulation
(12) reads as

######  Problem 3.

Seek $\mathbf{\boldsymbol{\lambda}}^{N}\in\mathbb{C}^{I(N)+1}$ (resp.
$\boldsymbol{\nu}^{N}\in\mathbb{C}^{I(N)+1}$) such that

$\displaystyle
V^{N}_{\Gamma}[k]\boldsymbol{\lambda}^{N}=\mathbf{g}_{d}^{N},\quad\text{(discrete
Dirichlet BIE)}.$ $\displaystyle
W^{N}_{\Gamma}[k]\boldsymbol{\nu}^{N}=\mathbf{g}_{n}^{N},\quad\text{(discrete
Neumann BIE)},$

where the respective discretization matrices elements are defined as

$\displaystyle(V^{N}_{\Gamma}[k])_{l,m}=\langle
V_{\Gamma}[k]q^{e,r}_{m},q_{l}^{e,r}\rangle_{\Gamma},\quad(W^{N}_{\Gamma}[k])_{l,n}=\langle
W_{\Gamma}[k]p^{o,r}_{m},p_{l}^{o,r}\rangle_{\Gamma},$

and the corresponding discrete right-hand sides are, for $l,m=0,\ldots,I(N)$,

$\displaystyle(g_{d}^{N})_{l}=\langle
g_{d},q^{e,r}_{l}\rangle_{\Gamma},\quad(g_{n}^{N})_{l}=\langle
g_{n},p^{o,r}_{l}\rangle_{\Gamma}.$

The discrete approximations for $\lambda,\nu$ (solutions of Problem 2) are
$\lambda^{N}_{e},\nu^{N}_{o}$ respectively. From standard Galerkin properties
we have the following result:

######  Lemma 4 (Theorem 4.29 in [30]).

There exists $N_{0}\in\mathbb{N}_{0}$ –potentially different for the weakly-
and hyper-singular BIEs– such that, for any $N\in\mathbb{N}$ with $N>N_{0}$,
the solutions $\boldsymbol{\lambda}^{N}$, and $\boldsymbol{\nu}^{N}$ of
Problem 3 exist, are unique and also we obtain quasi-optimality:

$\displaystyle\|\lambda-\lambda^{N}_{e}\|_{\widetilde{H}^{-\frac{1}{2}}(\Gamma)}$
$\displaystyle\lesssim\inf_{v\in\mathbb{Q}^{e}_{N}(\Gamma)}\|\lambda-v\|_{\widetilde{H}^{-\frac{1}{2}}(\Gamma)},\quad\|\nu-\nu^{N}_{o}\|_{\widetilde{H}^{\frac{1}{2}}(\Gamma)}\lesssim\inf_{v\in\mathbb{P}^{o}_{N}(\Gamma)}\|\nu-v\|_{\widetilde{H}^{\frac{1}{2}}(\Gamma)}.$

From these last estimates we obtain the rate of error convergence.

######  Theorem 3.

Given $g_{d},g_{n}\in\mathcal{D}(\Gamma)$, then for any regularity indices
$s_{d}>-\frac{1}{4}$, $s_{n}>\frac{1}{2}$, it holds that

$\displaystyle\|\lambda-\lambda^{N}_{e}\|_{\widetilde{H}^{-\frac{1}{2}}(\Gamma)}$
$\displaystyle\lesssim
N^{-\frac{1}{4}-s_{d}}\|\lambda\circ\mathbf{r}\|_{Q^{s_{d}}_{e}(\mathbb{D})},$
(14) $\displaystyle\|\nu-\nu^{N}_{o}\|_{\widetilde{H}^{\frac{1}{2}}(\Gamma)}$
$\displaystyle\lesssim
N^{\frac{1}{2}-s_{n}}\|\nu\circ\mathbf{r}\|_{P^{s_{n}}_{o}(\mathbb{D})}.$ (15)

###### Proof.

We prove only the Dirichlet case as the Neumann follows similar arguments.
Denote by $R$ a smooth non-zero function
$R(\mathbf{x}):=\dfrac{\|\mathbf{x}\|}{J_{\mathbf{r}}(\mathbf{x})}$. By using
Lemmas 1 and 4 one finds that

$\displaystyle\|\lambda-\lambda^{N}_{e}\|_{\widetilde{H}^{-\frac{1}{2}}(\Gamma)}$
$\displaystyle\lesssim\inf_{v\in\mathbb{Q}^{e}_{N}(\Gamma)}\|\lambda-v\|_{\widetilde{H}^{-\frac{1}{2}}(\Gamma)}=\inf_{v\in\mathbb{Q}^{e}_{N}(\mathbb{D})}\|\lambda\circ\mathbf{r}-Rv\|_{\widetilde{H}^{-\frac{1}{2}}(\mathbb{D})}.$

We can use the relation between Sobolev and auxiliary spaces (Lemma 3) to
obtain

$\displaystyle\|\lambda-\lambda^{N}_{e}\|_{\widetilde{H}^{-\frac{1}{2}}(\Gamma)}$
$\displaystyle\lesssim\inf_{v\in\mathbb{Q}^{e}_{N}(\mathbb{D})}\left\lVert\lambda\circ\mathbf{r}-Rv\right\rVert_{Q^{-\frac{1}{4}}_{e}(\mathbb{D})}\lesssim\inf_{v\in\mathbb{Q}^{e}_{N}(\mathbb{D})}\left\lVert\frac{\lambda}{R}-v\right\rVert_{Q^{-\frac{1}{4}}_{e}(\mathbb{D})}.$

Since $R^{-1}$ is smooth, the results follow by selecting
$v=\Pi^{e,q}_{N}\frac{\lambda}{R}$ and estimate (8) with $s=-1/4,r=s_{d}$. ∎

######  Remark 4.

Notice that the right-hand side in (14) is indeed finite for any
$s_{d}>-\frac{1}{4}$, $s_{n}>\frac{1}{2}$. This follows from the norm
definitions as

$\|\lambda\circ\mathbf{r}\|_{Q^{s_{d}}_{e}(\mathbb{D})}=\|L_{e}(w_{\mathbb{D}}\lambda\circ\mathbf{r})\|_{H^{s_{d}}(\mathbb{S})},$

and by Theorem 2, which ensures that $w_{\mathbb{D}}\lambda\circ\mathbf{r}$ is
smooth. Thus, implying a smooth lifting to the sphere. Using definitions of
the Sobolev spaces [25, Chapter 2], any Sobolev norm is finite. The same is
true for the Neumann case.

One concludes from Theorem 3 that the spectral method converges _super
algebraically_ , i.e. faster than any fixed negative power of $N$.

###  4.2 Matrix Computation

We now describe how matrix entries are computed. We start by detailing the
approximation of weakly-singular integrals that appear in the corresponding
matrix, namely

$\displaystyle V_{l,m}[k]=\langle
V_{\Gamma}[k]q^{e,r}_{m},q_{l}^{e,r}\rangle_{\Gamma}.$

Then, we briefly discuss how integrals for the hyper-singular BIO are obtained
using the same techniques, and also how regular entries are computed. Before
we proceed, we recall some notions of numerical quadrature and convergence.

####  4.2.1 Numerical Quadrature

Let $a<b$ two real numbers, and $f:(a,b)\rightarrow\mathbb{C}$. The Gauss-
Legendre quadrature rule approximate the integral of $f$ as a weighted sum of
point evaluations of $f$, the approximation is constructed as

$\displaystyle\int_{a}^{b}f(x)dx\approx\sum_{i=1}^{N_{q}}w^{L}_{i}f(x^{L}_{i}),$

where $N_{q}$ is the order of the quadrature, and
$(\\{w_{i}^{L}\\}_{i=1}^{N_{q}}$ ,$\\{x_{i}^{L}\\}_{i=1}^{N_{q}})$ are the
weights and points 444Notice that the points and weights depend on $a,b$, but
once given for a fixed interval they can be translated using a linear change
of variable. Consequently, we omit this dependence in notation.  of the
quadrature respectively. When $f$ is smooth, the quadrature converges with a
rate bounded as a function of $N_{q}$. In particular, by using the fact that
the quadrature rule is constructed to be exact for every polynomial up to some
degree and classical bounds for polynomial interpolation [36, Chap. 7 and 8],
one can establish that for $f$ with $(m+1)$ continuous derivatives555Better
convergence rates can be achieved for Gaussian quadrature rules [35, Theorem
3.6.24].:

$\displaystyle\left|\int_{a}^{b}f(x)dx-\sum_{i=1}^{N_{q}}w^{L}_{i}f(x^{L}_{i})\right|\lesssim
N_{q}^{-m}.$

Moreover, if $f$ admits and analytic extension to a Bernstein ellipse of
parameter $\rho$ in the complex plane, one can show that

$\displaystyle\left|\int_{a}^{b}f(x)dx-\sum_{i=1}^{N_{q}}w^{L}_{i}f(x^{L}_{i})\right|\lesssim\rho^{-N_{q}}.$

Whenever an integral of a function can be approximated with the same rates as
the last two, we say that the approximation is optimal. In particular, Gauss-
Legendre quadrature rules are optimal for any smooth function integral. Jacobi
quadrature rules are built as an approximation of the following family of
integrals:

$\displaystyle\int_{a}^{b}f(x)(x-a)^{\alpha}(b-x)^{\beta}dx\approx\sum_{i=1}^{N_{q}}w^{\alpha,\beta}_{i}f(x^{\alpha,\beta}_{i}),$

where, again, $N_{q}$ is the quadrature order,
$(\\{w_{i}^{\alpha,\beta}\\}_{i=1}^{N_{q}}$
,$\\{x_{i}^{\alpha,\beta}\\}_{i=1}^{N_{q}})$ are the pair of weights and
points, and $\alpha,\beta>-1$. This rule is also optimal, i.e. that the rate
of convergence is again $N_{q}^{-m}$ when $f$ has $(m+1)$ continuous
derivative and $\rho^{-N_{q}}$ when $f$ is $\rho-$analytic.

####  4.2.2 Approximation of weakly-singular integrals

Given a screen $\Gamma$ para-metrized by a function $\mathbf{r}$, we consider
the computation of integrals of the form:

$I_{m,l}[k]=\int_{\Gamma}\int_{\Gamma}\frac{\cos(k\|\mathbf{x}-\mathbf{x}^{\prime}\|)}{4\pi\|\mathbf{x}-\mathbf{x}^{\prime}\|}q_{m}^{e,r}(\mathbf{x}^{\prime})\overline{q_{l}^{e,r}(\mathbf{x})}d\mathbf{x}d\mathbf{x}^{\prime}.$

These integrals are associated with the real part of the weakly-singular BIO.
Its imaginary part is regular since the function
$\sin(k\|\mathbf{x}-\mathbf{x}^{\prime}\|)/\|\mathbf{x}-\mathbf{x}^{\prime}\|$
is smooth. Moreover, the cosine factor is smooth so from here onwards, we
assume that $k=0$ and denote $I_{m,l}[0]$, as $I_{m,l}$. By performing a
change of variable, this integral becomes

$I_{m,l}=\int_{\mathbb{D}}\int_{\mathbb{D}}\frac{1}{4\pi\|\mathbf{r}(\mathbf{x})-\mathbf{r}(\mathbf{x}^{\prime})\|}q^{e}_{m}(\mathbf{x}^{\prime})\overline{q^{e}_{l}(\mathbf{x})}d\mathbf{x}d\mathbf{x}^{\prime}.$

The first step is to take care of the kernel singularity, i.e. when
$\mathbf{x}=\mathbf{x}^{\prime}$. To this end, we make the following
parametrization:

$\displaystyle\mathbf{x}=r\mathbf{e}_{\theta},\quad\mathbf{x}^{\prime}=\mathbf{x}+\lambda
A(r,\beta)\mathbf{e}_{\theta+\beta},$ (16)

where

$\displaystyle A(r,\beta):=\sqrt{1-r^{2}\sin^{2}\beta}-r\cos\beta,$ (17)

represents the length of the segment whose direction is $\theta+\beta$, and
goes from point $\mathbf{x}$ to the boundary of the disk (see Figure 1).

Figure 1: Polar change of variables performed in (16).

The integral can be expressed as

$\displaystyle I_{m,l}$
$\displaystyle=\int_{-\frac{\pi}{2}}^{3\frac{\pi}{2}}\int_{0}^{1}\int_{-\frac{\pi}{2}}^{3\frac{\pi}{2}}\int_{0}^{1}G_{m,l}(\theta,r,\lambda,\beta)d\lambda
d\beta drd\theta,$ $\displaystyle G_{m,l}(\theta,r,\lambda,\beta)$
$\displaystyle:=\frac{\lambda
A^{2}(r,\beta)}{4\pi\|\mathbf{r}(\mathbf{x})-\mathbf{r}(\mathbf{x}+\lambda
A(r,\beta)\mathbf{e}_{\theta+\beta})\|}\frac{rP^{e}_{m}(\mathbf{x})\overline{P^{e}_{l}(\mathbf{x}^{\prime})}}{\sqrt{1-r^{2}}\sqrt{1-\|\mathbf{x}+\lambda
A(r,\beta)\mathbf{e}_{\theta+\beta}\|^{2}}}.$

Since $\mathbf{r}$ is smooth and its gradient has full rank, the factor

$\frac{\lambda
A^{2}(r,\beta)}{4\pi\|\mathbf{r}(\mathbf{x})-\mathbf{r}(\mathbf{x}+\lambda
A(r,\beta)\mathbf{e}_{\theta+\beta})\|},$

is at least bounded. The term $(1-r^{2})^{-\frac{1}{2}}$ while singular can be
tackled with the Jacobi rule. However, this is not enough since $A(r,\beta)$
and$\sqrt{1-\|\mathbf{x}+\lambda A(r,\beta)\mathbf{e}_{\theta+\beta}\|^{2}}$
also have singularities that prevent an optimal rate of convergence. The
following results characterize the behavior of these functions.

######  Lemma 5.

The function $A:[0,1]\times[\frac{-\pi}{2},\frac{3\pi}{2}]\rightarrow[0,2]$
has the following properties:

  1. (i)

$A(r,\beta)\geq 0$.

  2. (ii)

Partial derivatives discontinuities of $A(r,\beta)$ are located at
$\\{r=1,\beta=\frac{-\pi}{2},\frac{\pi}{2},\frac{3\pi}{2}\\}$.

  3. (iii)

For $|\beta|\leq\frac{\pi}{2}$, $A(r,\beta)\cos\beta\leq 1-r$.

###### Proof.

The first item is immediate by definition. For the second one, notice that the
discontinuities occur when the square-root parts vanish, i.e.

$(1-r^{2})+r^{2}\cos^{2}\beta=0.$

As a sum of two non negative terms, the singularities occurs when both terms
vanish, thus leading to our result directly. The last item is obtained by
direct evaluation at extreme points $\beta=\pm\frac{\pi}{2}$ and critical
values of $\beta$ such that $\dfrac{\partial
A(r,\beta)\cos\beta}{\partial\beta}=0$, or by directly using elementary
geometrical proprieties. ∎

######  Lemma 6.

Discontinuities of any partial derivative of the term

$\sqrt{1-\|\mathbf{x}+\lambda A(r,\beta)\mathbf{e}_{\theta+\beta}\|^{2}},$

occur at
$\\{\lambda=1\\}\cup\\{r=1,\lambda=0\\}\cup\\{r=1,\beta=\pm\frac{\pi}{2},\frac{3\pi}{2}\\}$.

###### Proof.

Again, critical points occur only when the term inside the square-root
vanishes. This term can be expressed as

$\displaystyle 1-\|\mathbf{x}+\lambda
A(r,\beta)\mathbf{e}_{\theta+\beta}\|^{2}=(1-\lambda)\left[(1+\lambda)(1-r^{2})-2\lambda
r\cos(\beta)A(r,\beta)\right],$ (18)

for which all previously listed points are critical. Hence, we are left to
check that no other critical points exist.

First, let us consider the case $\frac{\pi}{2}\leq\beta\leq\frac{3\pi}{2}$ so
that $-2\lambda r\cos\beta A(r,\beta)\geq 0$. Thus, the singularities are
characterized by $\lambda=1$ –because of the first factor–, and also the
points where

$\displaystyle(1+\lambda)(1-r^{2})=0,$ $\displaystyle 2\lambda r\cos\beta
A(r,\beta)=0.$

From this condition and Lemma 5 it is easy to see that no further critical
points take place. Now, if $|\beta|\leq\frac{\pi}{2}$, by the third item of
Lemma 5 we have that

$\displaystyle((1+\lambda)(1-r^{2})-2\lambda r\cos\beta
A(r,\beta))\geq(1+\lambda)(1-r^{2})-2\lambda r(1-r)),$

which implies that the singularities can only occur if $\lambda=1$ –because of
the first term–, or $r=1$, and the result follows. ∎

Based on the above considerations, we split the integrals $I_{m,l}$ according
to the singularities and critical points as follows. The detailed analysis for
the last three types is given in the Appendix B.

  * (a)

The first type is

$\displaystyle
I^{a}_{m,l}:=\int_{-\frac{\pi}{2}}^{\frac{3\pi}{2}}\int_{0}^{\frac{\sqrt{3}}{2}}\int_{-\frac{\pi}{2}}^{\frac{3\pi}{2}}\int_{0}^{1}G_{m,l}(\theta,r,\lambda,\beta)d\lambda
d\beta drd\theta,$

for which only a singularity of the form $(1-\lambda)^{-\frac{1}{2}}$ occurs.
This can numerically be treated using the corresponding Jacobi rule when
integrating in $\lambda$ and Gauss-Legendre for the integrals in
$\theta,r,\beta$, resulting in an optimal approximation.

  * (b)

The second type is

$\displaystyle
I^{b}_{m,l}:=\int_{-\frac{\pi}{2}}^{\frac{3\pi}{2}}\int_{\frac{\sqrt{3}}{2}}^{1}\int_{(\frac{-\pi}{3},\frac{\pi}{3})\cup(\frac{2\pi}{3},\frac{4\pi}{3})}\int_{\frac{1}{2}}^{1}G_{m,l}(\theta,r,\lambda,\beta)d\lambda
d\beta drd\theta,$

which has singularities of the form $(1-r)^{-\frac{1}{2}}$, and
$(1-\lambda)^{-\frac{1}{2}}$, so we use a Jacobi rule in these two variables
and Gauss-Legendre for $\theta$ and $\beta$.

  * (c)

The third one has critical points occurring as a combination of variables:

$\displaystyle
I^{c}_{m,l}:=\int_{-\frac{\pi}{2}}^{\frac{3\pi}{2}}\int_{\frac{\sqrt{3}}{2}}^{1}\int_{(\frac{-\pi}{2},\frac{-\pi}{3})\cup(\frac{\pi}{3},\frac{\pi}{2})\cup(\frac{\pi}{2},\frac{2\pi}{3})\cup(\frac{4\pi}{3},\frac{3\pi}{2})}\int_{0}^{\frac{1}{2}}G_{m,l}(\theta,r,\lambda,\beta)d\lambda
d\beta drd\theta,$ (19)

Here, critical points lie in $(\lambda,r,\beta)=(0,1,\beta)$, and
$(\lambda,r,\beta)=(\lambda,1,\beta_{0})$,
$\beta_{0}\in\\{\pm\frac{\pi}{2},\frac{3\pi}{2}\\}$. To tackle this, we will
use two polar change of variables (see B.2).

  * (d)

The fourth one is of the form

$\displaystyle
I^{d}_{m,l}:=\int_{-\frac{\pi}{2}}^{\frac{3\pi}{2}}\int_{\frac{\sqrt{3}}{2}}^{1}\int_{(\frac{-\pi}{3},\frac{\pi}{3})\cup(\frac{2\pi}{3},\frac{4\pi}{3})}\int_{0}^{\frac{1}{2}}G_{m,l}(\theta,r,\lambda,\beta)d\lambda
d\beta drd\theta,$ (20)

being sightly simpler that the previous case and only requiring one polar
change of variable in $\lambda$ and $r$ variables.

  * (e)

Finally, the last integral type is given by

$\displaystyle
I^{e}_{m,l}:=\int_{-\frac{\pi}{2}}^{\frac{3\pi}{2}}\int_{\frac{\sqrt{3}}{2}}^{1}\int_{(\frac{-\pi}{2},\frac{-\pi}{3})\cup(\frac{\pi}{3},\frac{\pi}{2})\cup(\frac{\pi}{2},\frac{2\pi}{3})\cup(\frac{4\pi}{3},\frac{3\pi}{2})}\int_{\frac{1}{2}}^{1}G_{m,l}(\theta,r,\lambda,\beta)d\lambda
d\beta drd\theta,$ (21)

which needs one polar change of variable in $(r,\cos\beta)$ and application of
a Jacobi rule for the integral in the $\lambda$ variable.

####  4.2.3 Approximation of smooth integrals

Smooth integrals can be of two forms:

$\displaystyle I^{f}_{m,l}$
$\displaystyle:=\int_{\Gamma}\int_{\Gamma}{G_{reg}}(\mathbf{x},\mathbf{x}^{\prime})q^{e,r}_{m}(\mathbf{x}^{\prime})\overline{q_{l}^{e,r}(\mathbf{x})}d\mathbf{x}d\mathbf{x}^{\prime},\quad
I^{g}_{m,l}:=\int_{\Gamma}g(\mathbf{x})\overline{q^{e,r}_{l}(\mathbf{x})}d\mathbf{x},$

where ${G_{reg}}$, and $g$ are smooth functions. The former comes from the
imaginary part of the fundamental solution, whereas the latter from testing
the right-hand side. We focus only in the second case, as the first one is
just a tensorisation. By using the screen parametrization, we obtain

$\displaystyle
I^{g}_{m,l}=\int_{-\frac{\pi}{2}}^{\frac{3\pi}{2}}\int_{0}^{1}g(\mathbf{r}(\mathbf{x}))\frac{\overline{p^{e}_{l}(\mathbf{x})}}{\sqrt{1-r^{2}}}rdrd\theta,$

where $\mathbf{x}=r\mathbf{e}_{\theta}$. Using the Gauss-Legendre rule when
integrating in $\theta$ and Jacobi for $r$ we obtain optimal rates of
convergence.

####  4.2.4 Approximation of hyper-singular Integrals

As it was pointed out, the discretization basis for the hyper-singular
operator BIO ($p^{o,r}_{m}$) vanishes on the boundary $\partial\Gamma$, and
consequently, the entries of the corresponding matrix can be computed using
the integration-by-parts formula [30, Corollary 3.3.24]:

$\displaystyle\begin{split}\langle
W_{\Gamma}[k]p^{o,r}_{m},p^{o,r}_{l}\rangle=\int_{\Gamma}\int_{\Gamma}\frac{e^{ik\|\mathbf{x}-\mathbf{x}^{\prime}\|}}{4\pi\|\mathbf{x}-\mathbf{x}^{\prime}\|}curl_{\Gamma}p^{o,r}_{m}(\mathbf{x}^{\prime})\cdot\overline{curl_{\Gamma}p^{o,r}_{l}(\mathbf{x})}d\mathbf{x}^{\prime}d\mathbf{x}\\\
-k^{2}\int_{\Gamma}\int_{\Gamma}\frac{e^{ik\|\mathbf{x}-\mathbf{x}^{\prime}\|}}{4\pi\|\mathbf{x}-\mathbf{x}^{\prime}\|}\widehat{\mathbf{n}}(\mathbf{x})\cdot\widehat{\mathbf{n}}(\mathbf{x}^{\prime})p^{o,r}_{m}(\mathbf{x}^{\prime})\overline{p^{o,r}_{l}(\mathbf{x})}d\mathbf{x}^{\prime}d\mathbf{x},\end{split}$
(22)

where $\widehat{\mathbf{n}}$ denote the unitary normal vector to $\Gamma$
(with a fixed orientation) and
$curl_{\Gamma}f=\widehat{\mathbf{n}}\times\nabla f$, whenever $f$ can be
extended to a neighborhood of $\Gamma$. We start by considering the second
integral on the right-hand side, we reduce the computations to $\mathbb{D}$
and obtain

$\displaystyle
k^{2}\int_{\mathbb{D}}\int_{\mathbb{D}}\frac{e^{ik\|\mathbf{r}(\mathbf{x})-\mathbf{r}(\mathbf{x}^{\prime})\|}}{4\pi\|\mathbf{r}(\mathbf{x})-\mathbf{r}(\mathbf{x}^{\prime})\|}\widehat{\mathbf{n}}(\mathbf{r}(\mathbf{x}))\cdot\widehat{\mathbf{n}}(\mathbf{r}(\mathbf{x}^{\prime}))p^{o}_{m}(\mathbf{x}^{\prime})\overline{p^{o}_{l}(\mathbf{x})}\frac{J_{\mathbf{r}}(\mathbf{x})}{\|\mathbf{x}\|}\frac{J_{\mathbf{r}}(\mathbf{x}^{\prime})}{\|\mathbf{x}^{\prime}\|}d\mathbf{x}^{\prime}d\mathbf{x},$

Since functions $p^{o}_{l}$ can be characterized as the product between a
smooth function and the weight function $w_{\mathbb{D}}$ the same change of
variables used for the weakly-singular case works for this integral666It is
necessary to change the parameters of the Jacobi quadrature rule, as now the
singularity is of the form $\sqrt{1-x^{2}}$, instead of $1/\sqrt{1-x^{2}}$..
For the first integral in the right-hand side in (22), we compute the surface
$curl$ operators. Using the parametrization of $\Gamma$, the explicit
expression

$(curl_{\Gamma}f)(\mathbf{r})=\frac{1}{J_{\mathbf{r}}}\left(\partial_{x_{2}}\mathbf{r}\partial_{x_{1}}(f\circ
r)-\partial_{x_{1}}\mathbf{r}\partial_{x_{2}}(f\circ\mathbf{r})\right)$

arises, where $\partial_{x_{1}},\partial_{x_{2}}$ denote the partial
derivatives with respect to the arguments of the parametrization $\mathbf{r}$.
In our implementation, $\mathbf{r}$ is given in polar coordinates and so
$x_{1}$ and $x_{2}$ are the radial and angular variables, respectively. The
function $f\circ\mathbf{r}$ corresponds to $p^{o}_{m}$ for some $m$,and thus,
we require their partial derivatives. Moreover, using the two-indices
representation we can write $p^{o}_{m}=p^{\mathcal{l}}_{\mathcal{m}}$, for a
pair of integers $\mathcal{l},\mathcal{m}$ such that $\mathcal{m}+\mathcal{l}$
is odd. The angular derivative is given by

$\partial_{\theta}p^{\mathcal{l}}_{\mathcal{m}}=\partial_{\theta}\left(C^{\mathcal{l}}_{\mathcal{m}}\mathrm{P}^{\mathcal{l}}_{|\mathcal{m}|}(\sqrt{1-r^{2}})e^{icm\theta}\right)=i\mathcal{m}C^{\mathcal{l}}_{\mathcal{m}}\mathrm{P}^{\mathcal{l}}_{|\mathcal{m}|}(\sqrt{1-r^{2}})e^{i\mathcal{m}\theta}=i\mathcal{m}p^{\mathcal{l}}_{\mathcal{m}},$

where the first equality is the explicit definition of the projected basis
(2). We can express this in terms of the basis used for the discretization of
the weakly-singular BIO as we have the following recursive relation (see [2,
15.87]):

$\displaystyle\mathrm{P}^{\mathcal{l}}_{\mathcal{m}}(x)=\frac{-\sqrt{1-x^{2}}}{2\mathcal{m}x}\left(\mathrm{P}^{\mathcal{l}}_{\mathcal{m}+1}(x)+(\mathcal{l}+\mathcal{m})(\mathcal{l}-\mathcal{m}+1)\mathrm{P}^{\mathcal{l}}_{\mathcal{m}-1}(x)\right),\quad\mathcal{m}\neq
0$

so we conclude that

$\displaystyle\partial_{\theta}p^{\mathcal{l}}_{\mathcal{m}}=\frac{-ir}{2}\left(a^{\mathcal{l}}_{\mathcal{m}}q^{\mathcal{l}}_{\mathcal{m}+1}+b^{\mathcal{l}}_{\mathcal{m}}q^{\mathcal{l}}_{\mathcal{m}-1}\right),$
(23)

where

$\displaystyle
a^{\mathcal{l}}_{\mathcal{m}}:=e^{-i\theta}\begin{cases}\sqrt{(\mathcal{l}-|\mathcal{m}|)(\mathcal{l}+|\mathcal{m}|+1)},\quad\mathcal{m}\geq
0\\\
-\sqrt{(\mathcal{l}+|\mathcal{m}|)(\mathcal{l}-|\mathcal{m}|+1)},\quad\mathcal{m}<0\end{cases}$
(24) $\displaystyle
b^{\mathcal{l}}_{\mathcal{m}}:=e^{i\theta}\begin{cases}\sqrt{(\mathcal{l}+|\mathcal{m}|)(\mathcal{l}-|\mathcal{m}|+1)},\quad\mathcal{m}\geq
0\\\
-\sqrt{(\mathcal{l}-|\mathcal{m}|)(\mathcal{l}+|\mathcal{m}|+1)},\quad\mathcal{m}<0\end{cases}$
(25)

Notice that, for $\mathcal{m}=0$, the derivative is zero and this is also true
for expression (23). Since $\mathcal{m}+\mathcal{l}$ is odd, we have that
$q^{\mathcal{l}}_{\mathcal{m}+1}$ and $q^{\mathcal{l}}_{\mathcal{m}-1}$ are
even functions. For the derivative with respect to $r$, we need the derivative
of the associated Legendre functions, given by the following recursion (see
[2, 15.91]):

$\displaystyle\frac{d\mathrm{P}^{\mathcal{l}}_{\mathcal{m}}(x)}{dx}=\frac{1}{2\sqrt{1-x^{2}}}\left(-\mathrm{P}^{\mathcal{l}}_{\mathcal{m}+1}(x)+(\mathcal{l}+\mathcal{m})(\mathcal{l}-\mathcal{m}+1)\mathrm{P}^{\mathcal{l}}_{\mathcal{m}-1}(x)\right),$

Hence, we obtain

$\partial_{r}p^{\mathcal{l}}_{\mathcal{m}}=\frac{1}{2}\left(a^{\mathcal{l}}_{\mathcal{m}}q^{\mathcal{l}}_{\mathcal{m}+1}-b^{\mathcal{l}}_{\mathcal{m}}q^{\mathcal{l}}_{\mathcal{m}-1}\right),$

where $a^{\mathcal{l}}_{\mathcal{m}},b^{\mathcal{l}}_{\mathcal{m}}$ are
defined as in (24). Again, we have expressed the derivative in terms of the
basis of the weakly-singular case. We conclude that for the computation of the
hyper-singular BIO only minor modifications respect to the weakly-singular one
are needed. These modifications are the change of the parameter of the
respective Jacobi rule, the inclusion of the product of the normal vectors,
and the extra smooth factors
$e^{\pm{\imath}\theta},\partial_{x_{j}}\mathbf{r}\cdot\partial_{x_{j}^{\prime}}\mathbf{r}$
($j,j^{\prime}=1,2$) that have to be included in the kernel function. We have
omitted the details of smooth integrals associated with the hyper-singular BIO
as they do no present any extra challenge.

###  4.3 Numerical Implementation

Throughout this section we denote by
$\widehat{N}=I(N)+1=N+\frac{N(N+1)}{2}+1$, the number of degrees of freedom
when using the spaces $\mathbb{Q}^{N}_{e}(\Gamma)$,
$\mathbb{P}^{N}_{o}(\Gamma)$ for the discretization of the underlying integral
equations. Every integral needed to assembly the matrix discretization of the
weakly-singular BIO ($I^{a}$, $I^{b}$, $I^{c}$, $I^{d}$, $I^{e}$, $I^{f}$ ) is
a four-dimensional integral. The total number of these integral computations
can be reduced using the matrix symmetries777The matrix associated with the
real part of the fundamental solution ($I^{\alpha},\alpha\in\\{a,b,c,d,e\\}$)
is symmetric, while the imaginary one ($I^{f}$) is anti-symmetric..
Consequently, only $\frac{\widehat{N}(\widehat{N}+1)}{2}$ combinations are
needed instead of $\widehat{N}^{2}$. Furthermore, since
$p_{-m}^{l}=(-1)^{m}\overline{p_{m}^{l}}$ the actual number of interactions
needed to compute the weakly-singular BIO is
$\frac{1}{4}\left(\widehat{N}+\frac{N+2}{2}\right)\left(\frac{3}{2}\widehat{N}-\frac{N-2}{4}\right)$,
assuming $N$ even888If $N$ is odd, the computational cost is
$\frac{1}{4}\left(\widehat{N}+\frac{N+1}{2}\right)\left(\frac{3}{2}\widehat{N}-\frac{N-3}{4}\right)$..

A four-dimensional integral computed by tensorized 1D-quadrature rules, with
parameters $N_{q}^{1},N_{q}^{2}$, $N_{q}^{3},N_{q}^{4}$, has a computational
cost of $\mathcal{O}(\Pi_{j=1}^{4}N_{q}^{j})$ operations and evaluations. To
compute integrals arising from the weakly singular BIO
$I^{\alpha}\in\\{b,c,d,e\\}$, we denote by $N_{q}^{\theta}$ the number of
points for the $\theta$ variable, and $N_{q}^{\alpha}$ the number of points
for the other three variables depending on $\alpha$. For $\alpha=a$, we use
$N_{q}^{\theta}$ for variables $\theta$ and $\beta$, and $N^{a}_{q}$ for the
rest. For the hyper-singular case, same rules apply.

For the smooth integrals $\alpha\in\\{f,g\\}$ we could use $N_{q}^{\theta}$
for the $\theta$ and $\theta^{\prime}$ variables, and
$N^{\alpha}_{q},\alpha\in\\{f,g\\}$ for $\mathbf{r}$, and
$\mathbf{r}^{\prime}$. However, in practice it is better to reformulate the
integrals onto the sphere, where the basis correspond to spherical harmonics,
and approximate the integrals using the spherical harmonics transforms. In
particular we use the implementation detailed in [31].

##  5 Full Discretizaton Error Analysis

The rate of convergence of the spectral Galerkin discretization method was
already established in Theorem 3. Yet, this does not illustrate the
performance of the fully discrete method as extra error terms appears as
consequence of the quadrature approximation of integral terms. Thus, we first
measure the perturbation in error convergence rates due to quadrature error.

###  5.1 Quadrature Error

For the sake of brevity, we denote $Q^{\alpha}$ the quadrature approximation
of $I^{\alpha}$, $\alpha\in\\{a,b,c,d,e,f,g\\}$, defined in Section 4.2.2.
Since we assume that the screen is parametrized by a $\rho-$analytic function,
and the approximation is optimal, we have that

$\displaystyle|I^{\alpha}-Q^{\alpha}|\lesssim\rho^{-N_{q}^{\theta}}+\rho^{-N_{q}^{\alpha}}.$
(26)

While this bound is precise in terms of how the quadrature error decreases
with increasing number of quadrature points, the unspecified constant depends
on the trial and test basis indices. Hence, since the rate of convergence
depends on the number of trial functions, we need a more detailed quadrature
error analysis considering the exact index of the trial and test basis. For
this, let us consider the canonical integral

$\displaystyle
I^{l}_{m}=\int_{0}^{1}\int_{-\frac{\pi}{2}}^{\frac{3\pi}{2}}g(\mathbf{x})p^{l}_{m}(\mathbf{x})rd\theta
dr,$

where $\mathbf{x}=r\mathbf{e}_{\theta}$, $g$ is $\rho-$analytic in
$(r,\theta)$ and $l+m$ is even. It is enough to consider this case, as all
integrals discussed in Section 4.2.2 can be reduced to this form by means of
analytic change of variables to tensorization of integral. Denote by
$Q^{l}_{m}$ the quadrature approximation of $I^{l}_{m}$ obtained by a Gauss-
Legendre rule in both variables, with $N^{\theta}_{q}$ points in the $\theta$
variable, and $N_{q}$ points in the $r$ variable.

We denote by $\mathcal{E}_{\rho}[a,b]$ the region enclosed by the Bernstein
ellipse in the complex plane with foci in $a,b$ and parameter $\rho$. Now, we
recall the classical error bound for analytic integrands.

######  Theorem 4 (Theorem 5.3.13 in [30]).

If $m+l$ is even, for $g$ $\rho$-analytic in $[0,1]$ in the radial variable
and $\rho$-analytic in $[-\frac{\pi}{2},\frac{3\pi}{2}]$ for the angular one,
it holds that

$\displaystyle\begin{split}|I^{l}_{m}-Q^{l}_{m}|\lesssim(2\rho)^{-2N_{q}^{\theta}}\max_{r\in[0,1]}\max_{z\in\partial\mathcal{E}_{\rho}[\frac{-\pi}{2},\frac{3\pi}{2}]}|g(r,z)p^{l}_{m}(r,z)|\\\
+(2\rho)^{-2N_{q}}\max_{\theta\in[\frac{-\pi}{2},\frac{3\pi}{2}]}\max_{z\in\partial\mathcal{E}_{\rho}[0,1]}|g(z,\theta)p^{l}_{m}(z,\theta)|\end{split}$

where the unspecified constant does not depend on the integrand of
$I^{l}_{m}$.

Since $g$ is assumed to be known, we can further simplify the error bound as

$\displaystyle\begin{split}|I^{l}_{m}-Q^{l}_{m}|\lesssim(2\rho)^{-2N_{q}^{\theta}}\max_{r\in[0,1]}\max_{z\in\partial\mathcal{E}_{\rho}[\frac{-\pi}{2},\frac{3\pi}{2}]}|p^{l}_{m}(r,z)|\\\
+(2\rho)^{-2N_{q}}\max_{\theta\in[\frac{-\pi}{2},\frac{3\pi}{2}]}\max_{z\in\partial\mathcal{E}_{\rho}[0,1]}|p^{l}_{m}(z,\theta)|,\end{split}$

with a constant that now depends on $g$.

######  Corollary 2.

Under hypothesis of Theorem 4, for the integral

$\displaystyle\widetilde{I}^{l}_{m}:=\int_{0}^{1}\int_{-\frac{\pi}{2}}^{\frac{3\pi}{2}}g(\mathbf{x})q^{l}_{m}(\mathbf{x})rd\theta
dr,$

it holds that

$\displaystyle\begin{split}|\widetilde{I}^{l}_{m}-\widetilde{Q}^{l}_{m}|\lesssim(2\rho)^{-2N_{q}^{\theta}}\max_{r\in[0,1]}\max_{z\in\partial\mathcal{E}_{\rho}[\frac{-\pi}{2},\frac{3\pi}{2}]}|p^{l}_{m}(r,z)|\\\
+(2\rho)^{-2N_{q}}\max_{\theta\in[\frac{-\pi}{2},\frac{3\pi}{2}]}\max_{z\in\partial\mathcal{E}_{\rho}[0,1]}|p^{l}_{m}(z,\theta)|,\end{split}$

where $\widetilde{Q}^{l}_{m}$ denotes the quadrature approximation using a
Gauss-Legendre rule with $N_{q}^{\theta}$ points in $\theta$, and a Jacobi
rule with $N_{q}$ points in the $r$ variable.

We now proceed to estimate maxima of analytic extensions for the functions
$p^{l}_{m}$. Remember that the explicit definition (2) (see Section 2.3):

$\displaystyle
p^{l}_{m}(r,\theta)=C^{l}_{m}\mathrm{P}^{l}_{|m|}(\sqrt{1-r^{2}})e^{im\theta},$

where $|C^{l}_{m}|=\sqrt{\frac{(2l+1)(l-|m|)!}{2\pi(l+|m|)!}}$, with
$\mathrm{P}^{l}_{|m|}(\sqrt{1-r^{2}})$ smooth in the $r$ variable. By using
[23, Theorem 3], one deduces that

$\displaystyle\max_{r\in[0,1]}\max_{z\in\partial\mathcal{E}_{\rho}[\frac{-\pi}{2},\frac{3\pi}{2}]}|p^{l}_{m}(r,z)|\lesssim\sqrt{2l+1}e^{|m|\frac{\pi}{2}(\rho-\rho^{-1})}\lesssim\sqrt{2l+1}e^{l\pi\rho},$

where the last inequality follows for $|m|\leq l$ and $\rho>1$. On the other
hand, the second term can be bounded as

$\displaystyle\max_{\theta\in[\frac{-\pi}{2},\frac{3\pi}{2}]}\max_{z\in\partial\mathcal{E}_{\rho}[0,1]}|p^{l}_{m}(z,\theta)|\leq|C^{l}_{m}|\max_{z\in\partial\mathcal{E}_{\rho}[0,1]}|\mathrm{P}^{l}_{|m|}(\sqrt{1-z^{2}})|.$

We can use the Rodríguez formula to express the associated Legendre function
in terms of Legendre polynomials:

$\displaystyle|\mathrm{P}^{l}_{|m|}(\sqrt{1-z^{2}})|=|z^{m}|\left|\left(\frac{d^{m}}{dx^{m}}\mathrm{P}^{l}(x)\right)|_{x=\sqrt{1-z^{2}}}\right|,$

where $\mathrm{P}^{l}$ denotes the $l$th Legendre polynomial. Obviously,
$|z|<\frac{\rho+\rho^{-1}}{4}$, for every
$z\in\partial\mathcal{E}_{\rho}[0,1]$. Moreover,

$\displaystyle\max_{z\in\partial\mathcal{E}_{\rho}[0,1]}\left|\left(\frac{d^{m}}{dx^{m}}\mathrm{P}^{l}(x)\right)\Big{|}_{x=\sqrt{1-z^{2}}}\right|=\max_{z\in
A_{\rho}}\left|\left(\frac{d^{m}}{dx^{m}}\mathrm{P}^{l}(x)\right)\Big{|}_{x=z}\right|,$

where $A_{\rho}$ is the image of $\partial\mathcal{E}_{\rho}[0,1]$ under the
transformation $(1-x^{2})^{\frac{1}{2}}$, where the square root has to be
understood as the pre-image of the square function. Since $\mathrm{P}^{l}$ are
polynomials and using the maximum modulus principle, there exists
$\widehat{\rho}>1$ such that

$\displaystyle\max_{z\in
A_{\rho}}\left|\left(\frac{d^{m}}{dx^{m}}\mathrm{P}^{l}(x)\right)\Big{|}_{x=z}\right|\leq\max_{z\in\partial\mathcal{E}_{\widehat{\rho}}[-1,1]}\left|\left(\frac{d^{m}}{dx^{m}}\mathrm{P}^{l}(x)\right)\Big{|}_{x=z}\right|.$

Furthermore, using the Cauchy integral formula we have that

$\displaystyle\frac{d^{m}}{dz^{m}}\mathrm{P}^{l}(z)=\frac{m!}{2i\pi}\int_{\partial\mathcal{E}_{2\widehat{\rho}}[-1,1]}\frac{\mathrm{P}^{l}(x)}{(z-x)^{m+1}}dx,\quad\forall
z\in\partial\mathcal{E}_{\widehat{\rho}}[-1,1].$

Hence, by using [38][Theorem 4.1] we have

$\displaystyle\max_{z\in\partial\mathcal{E}_{\widehat{\rho}}[-1,1]}\left|\left(\frac{d^{m}}{dx^{m}}\mathrm{P}^{l}(x)\right)|_{x=z}\right|\leq\frac{m!}{2\pi}L(\partial\mathcal{E}_{2\widehat{\rho}}[-1,1])\mathrm{P}^{l}(\widehat{\rho}+(4\widehat{\rho})^{-1})\max_{x\in\partial\mathcal{E}_{2\widehat{\rho}}[-1,1]}\frac{1}{|z-x|^{m+1}},$

where $L(\partial\mathcal{E}_{2\widehat{\rho}}[-1,1])$ is the length of the
corresponding ellipse, and as such, it can be estimated as
$L(\partial\mathcal{E}_{2\widehat{\rho}}[-1,1])\lesssim\widehat{\rho}$. Also,
notice that the minimum distance between
$\partial\mathcal{E}_{2\widehat{\rho}}[-1,1]$ and
$\partial\mathcal{E}_{\widehat{\rho}}[-1,1]$ is larger999This can be shown
using elementary geometrical computations. than $\widehat{\rho}$. Thus, one
has

$\displaystyle\max_{z\in\partial\mathcal{E}_{\widehat{\rho}}[-1,1]}\left|\left(\frac{d^{m}}{dx^{m}}\mathrm{P}^{l}(x)\right)|_{x=z}\right|\lesssim\frac{m!}{\widehat{\rho}^{m}}\mathrm{P}^{l}(\widehat{\rho}+(4\widehat{\rho})^{-1})\lesssim\frac{m!}{2^{l}}\widehat{\rho}^{l-m},$

wherein the rightmost inequality follows for $l$ is large as the polynomial is
dominated by the monomial of greatest degree. Finally, combining all the
computations for $l$ large enough, we have that

$\displaystyle\max_{\theta\in[\frac{-\pi}{2},\frac{3\pi}{2}]}\max_{z\in\partial\mathcal{E}_{\rho}[0,1]}|p^{l}_{m}(z,\theta)|\lesssim|C^{l}_{|m|}|m!\widehat{\rho}^{l}\lesssim\sqrt{\frac{(2l+1)(l!)^{2}}{(2l)!}}\widehat{\rho}^{l}\lesssim\sqrt{2l+1}\widehat{\rho}^{l}$

and the quadrature error is then bounded as

$\displaystyle|{I}^{l}_{m}-{Q}^{l}_{m}|\lesssim\sqrt{2l+1}\left[(2\rho)^{-2N_{q}^{\theta}}e^{\pi
l\rho}+(2\rho)^{-2N_{q}}\widehat{\rho}^{l}\right].$

This bound be further simplified to

$\displaystyle|{I}^{l}_{m}-{Q}^{l}_{m}|\lesssim\sqrt{2l+1}\left[(2\rho)^{-2N_{q}^{\theta}}+(2\rho)^{-2N_{q}}\right]\widetilde{\rho}^{l}$

where $\widetilde{\rho}>1$. We can summarize quadrature error bound in the
following result:

######  Lemma 7.

There is $N_{0}\in\mathbb{N}$ such that given integers $m,l$, $|m|<l$ and
$l>N_{0}$, and
$g:[0,1]\times[\frac{-\pi}{2},\frac{3\pi}{2}]\rightarrow\mathbb{C}$
$\rho-$analytic in both variables. For the integral

$I^{l}_{m}:=\int_{0}^{1}\int_{\frac{-\pi}{2}}^{\frac{3\pi}{2}}g(\mathbf{x})p^{l}_{m}(\mathbf{x})rd\theta
dr,$

we have the error bound

$\displaystyle|{I}^{l}_{m}-{Q}^{l}_{m}|\lesssim\sqrt{2l+1}\left((2\rho)^{-2N_{q}^{\theta}}+(2\rho)^{-2N_{q}}\right)\widetilde{\rho}^{l}$

where $Q^{l}_{m}$ denotes the approximation by Gauss-Legendre of order
$N_{q},N_{q}^{\theta}$, for both variables accordingly, $\widetilde{\rho}>1$,
and the unspecified constant does not depend on $l,m$.

######  Corollary 3.

Under the hypothesis and notations of Lemma 7, given another pair of integers
$l^{\prime},m^{\prime}$ such that $|m^{\prime}|<l^{\prime}$ and
$l^{\prime}>N_{0}$, and
$G(\mathbf{x},\mathbf{x}^{\prime}):\left([0,1]\times[\frac{-\pi}{2},\frac{3\pi}{2}]\right)^{2}\rightarrow\mathbb{C}$,
for the integral

$\displaystyle
I^{l,l^{\prime}}_{m,m^{\prime}}:=\int_{0}^{1}\int_{-\frac{\pi}{2}}^{\frac{3\pi}{2}}\int_{0}^{1}\int_{-\frac{\pi}{2}}^{\frac{3\pi}{2}}G(\mathbf{x},\mathbf{x}^{\prime})p^{l}_{m}(\mathbf{x})p^{l^{\prime}}_{m^{\prime}}(\mathbf{x}^{\prime})rr^{\prime}d\theta
drd\theta^{\prime}dr^{\prime},$

it holds that

$\displaystyle|I^{l,l^{\prime}}_{m,m^{\prime}}-Q^{l,l^{\prime}}_{m,m^{\prime}}|\lesssim\sqrt{(2l+1)(2l^{\prime}+1)}\left((2\rho)^{-2N_{q}^{\theta}}+(2\rho)^{-2N_{q}}\right)\widetilde{\rho}^{l}\widetilde{\rho}^{\prime
l^{\prime}}$

where $Q^{l,l^{\prime}}_{m,m^{\prime}}$ denotes the Gauss-Legendre quadrature
of orders $N_{q}$ in $r,r^{\prime}$ and $N_{q}^{\theta}$ for
$\theta,\theta^{\prime}$, $\widetilde{\rho}>1$ and
$\widetilde{\rho}^{\prime}>1$, and the unspecified constant does not depend of
$l,m,l^{\prime},m^{\prime}$.

###  5.2 Fully Discrete Error Analysis:

Recall Problem 3, where the unknowns are vectors of dimension
$I(N)=\frac{N(N+1)}{2}+N$. Let us now consider the same problem where the
corresponding matrices and righ-hand sides are approximated with the
quadrature method detailed in Section 4.2.

######  Problem 4.

Find $\boldsymbol{\lambda}^{N,q},\boldsymbol{\nu}^{N,q}\in\mathbb{C}^{I(N)+1}$
such that

$\displaystyle
V^{N,q}[k]\boldsymbol{\lambda}^{N,q}=\mathbf{g}_{d}^{N,q},\quad\text{(Fully
discrete Dirichlet BIE)},$ (27) $\displaystyle
W^{N,q}[k]\boldsymbol{\nu}^{N,q}=\mathbf{g}_{n}^{N,q},\quad\text{(Fully
discrete Neumann BIE)},$ (28)

where $V^{N,q}[k]$ (resp.
$W^{N,q}[k],\mathbf{g}_{d}^{N,q},\mathbf{g}_{n}^{N,q}$) is the quadrature
approximation to $V^{N}[k]$ (resp.
$W^{N}[k],\mathbf{g}_{d}^{N},\mathbf{g}_{n}^{N}$) constructed as described in
Section 4.2.

The approximations obtained from the fully discrete problems are written

$\displaystyle\lambda^{N,q}_{e}=\sum_{m=0}^{I(N)}\lambda^{N,q}_{m}q_{m}^{e,r},\quad\nu^{N,q}_{o}=\sum_{m=0}^{I(N)}\nu^{N,q}_{m}p_{m}^{o,r}.$

We refer to Problem 4 as fully discrete problems. For their analysis, we
detail the Dirichlet case as the Neumann case follows similar ideas with
results for both is reported at the end of the section.

We first estimate the quadrature error between matrices $V^{N}[k]-V^{N,q}[k]$.
Recall the notation introduced at beginning of Section 4.1. For simplicity,
the number quadrature points is determined by only two variables
$N_{q}^{\theta}$ and $N_{q}$. Following the notation in Section 4.3, this
corresponds to the simpler case $N_{q}^{\alpha}=N_{q}$ for every
$\alpha\in\\{a,b,c,d,e,f\\}$. For different values of $N_{q}^{\alpha}$, the
results are essentially the similar but we forgo this analysis for the sake of
brevity.

######  Lemma 8.

There is $N_{0}\in\mathbb{N}$, such that for $N>N_{0}$, given vectors
$\boldsymbol{\mu}^{N}$, $\boldsymbol{\psi}^{N}\in\mathbb{C}^{I(N)+1}$, we have
that

$\displaystyle\begin{split}\left|\boldsymbol{\psi}^{N}\cdot\left(V^{N}[k]-V^{N,q}[k]\right)\boldsymbol{\mu}^{N}\right|\lesssim\widetilde{\rho}^{2N}N^{4}\left[(2\rho)^{-2N_{q}}+(2\rho)^{-2N^{\theta}_{q}}\right]\\\
\times\|\mu_{e}^{N}\|_{\widetilde{H}^{-\frac{1}{2}}(\Gamma)}\|\psi_{e}^{N}\|_{\widetilde{H}^{-\frac{1}{2}}(\Gamma)},\end{split}$

where, as before, $\rho>1,\widetilde{\rho}>1$.

###### Proof.

Expanding matrix products yields

$\displaystyle\begin{split}&\left|\boldsymbol{\psi}^{N}\cdot\left(V^{N}[k]-V^{N,q}[k]\right)\boldsymbol{\mu}^{N}\right|=\\\
&\left|\sum_{l^{\prime}=0}^{N}\sum_{\begin{subarray}{c}m^{\prime}=-l^{\prime}\\\
m^{\prime}+l^{\prime}\text{
even}\end{subarray}}^{l^{\prime}}\overline{\psi^{N}_{I_{e}(l^{\prime},m^{\prime})}}\sum_{l=0}^{N}\sum_{\begin{subarray}{c}m=-l\\\
m+l\text{
even}\end{subarray}}^{l}\left(V^{N}[k]-V^{N,q}[k]\right)_{I_{e}(l^{\prime},m^{\prime}),I_{e}(l,m)}\mu^{N}_{I_{e}(l,m)}\right|,\end{split}$

where the index $I_{e}(l,m)=\frac{l(l+1)+(l+m)}{2}$ was defined in Section
2.3. Now, $V^{N}[k]$ (see Section 4.1) is a matrix whose entries can be
written as the sum of integrals $I^{a},I^{b},I^{c},I^{d},I^{e},I^{f}$, and
$V^{N,q}[k]$ is the sum corresponding to quadratures approximations
$Q^{a},Q^{b},Q^{c},Q^{d},Q^{e},Q^{f}$ as described in Section 4.2.2. By
construction, we can use101010Each of the integrals
$I^{\alpha},\alpha\in\\{a,b,c,d,e,f\\}$ have different integration domain but
they could be fixed using smooth window functions that extend the integrands
to the required domain so the application of the Corollary directly.
Corollary 3 and obtain

$\displaystyle\begin{split}&\left|\boldsymbol{\psi}^{N}\cdot\left(V^{N}[k]-V^{N,q}[k]\right)\boldsymbol{\mu}^{N}\right|\lesssim\left[(2\rho)^{-2N_{q}^{\theta}}+(2\rho)^{-2N_{q}}\right]\\\
&{\times}\left|\sum_{l^{\prime}=0}^{N}\sum_{\begin{subarray}{c}m^{\prime}=-l^{\prime}\\\
m^{\prime}+l^{\prime}\text{
even}\end{subarray}}^{l^{\prime}}\sum_{l=0}^{N}\sum_{\begin{subarray}{c}m=-l\\\
m+l\text{
even}\end{subarray}}^{l}\sqrt{(2l+1)(2l^{\prime}+1)}\left(\widetilde{\rho}^{{}^{\prime}l^{\prime}}\widetilde{\rho}^{l}\right)\overline{\psi^{N}_{I_{e}(l^{\prime},m^{\prime})}}\mu^{N}_{I_{e}(l,m)}\right|\end{split}.$

Redefining
$\widetilde{\rho}:=\max\\{\widetilde{\rho}^{\prime},\widetilde{\rho}\\}$ leads
to

$\displaystyle\begin{split}&|\boldsymbol{\psi}^{N}\cdot\left(V^{N}[k]-V^{N,q}[k]\right)\boldsymbol{\mu}^{N}|\lesssim\widetilde{\rho}^{2N}\left[(2\rho)^{-2N_{q}^{\theta}}+(2\rho)^{-2N_{q}})\right]\\\
&\times\left|\sum_{l^{\prime}=0}^{N}\sum_{\begin{subarray}{c}m^{\prime}=-l^{\prime}\\\
m^{\prime}+l^{\prime}\text{
even}\end{subarray}}^{l^{\prime}}\sum_{l=0}^{N}\sum_{\begin{subarray}{c}m=-l\\\
m+l\text{
even}\end{subarray}}^{l}\sqrt{(2l+1)(2l^{\prime}+1)}\overline{\psi^{N}_{I_{e}(l^{\prime},m^{\prime})}}\mu^{N}_{I_{e}(l,m)}\right|.\end{split}$

Applying the Cauchy-Schwartz inequality one obtains

$\displaystyle\begin{split}&|\boldsymbol{\psi}^{N}\cdot\left(V^{N}[k]-V^{N,q}[k]\right)\boldsymbol{\mu}^{N}|\lesssim\widetilde{\rho}^{2N}\left[(2\rho)^{-2N_{q}^{\theta}}+(2\rho)^{-2N_{q}})\right]\\\
&\times
N^{4}\left(\sum_{l^{\prime}=0}^{N}\sum_{\begin{subarray}{c}m^{\prime}=-l^{\prime}\\\
m^{\prime}+l^{\prime}\text{
even}\end{subarray}}^{l^{\prime}}\frac{(\psi^{N}_{I_{e}(l^{\prime},m^{\prime})})^{2}}{l^{\prime}+1}\right)^{\frac{1}{2}}\left(\sum_{l=0}^{N}\sum_{\begin{subarray}{c}m=-l\\\
m+l\text{
even}\end{subarray}}^{l}\frac{(\mu^{N}_{I_{e}(l,m)})^{2}}{l+1}\right)^{\frac{1}{2}}\end{split},$

and by the auxiliary spaces norms definitions, it holds that

$\displaystyle\begin{split}|\boldsymbol{\psi}^{N}\cdot\left(V^{N}[k]-V^{N,q}[k]\right)\boldsymbol{\mu}^{N}|&\lesssim\widetilde{\rho}^{2N}\left[(2\rho)^{-2N_{q}^{\theta}}+(2\rho)^{-2N_{q}})\right]\\\
&\times
N^{4}\|\psi^{N}\circ\mathbf{r}\|_{Q_{e}^{-\frac{1}{2}}(\mathbb{D})}\|\mu^{N}\circ\mathbf{r}\|_{Q_{e}^{-\frac{1}{2}}(\mathbb{D})}.\end{split}$

The results follows directly from Lemmas 3 and 1. ∎

We notice that for any fixed value of $N$ the error term

$\displaystyle\widetilde{\rho}^{2N}\left[(2\rho)^{-2N_{q}^{\theta}}+(2\rho)^{-2N_{q}})\right]N^{4}$

goes to zero as the quadrature order increases. Thus, we can use the standard
Strang’s Lemma to bound the fully discrete error.

######  Theorem 5 (Theorem 4.2.11 in [30]).

There exists $N_{0}>0$ such that for every $N>N_{0}$, there is a $N_{q,0}$
that depends on $N$ such that, if $N_{q}^{\theta}$, $N_{q}$ are both greater
than $N_{q,0}$, the fully discrete Dirichlet problem 4 has a unique solution,
and the following error estimate follows for $s>-\frac{1}{4}$:

$\displaystyle\begin{split}\|\lambda-\lambda_{e}^{N,q}\|_{\widetilde{H}^{-\frac{1}{2}}(\Gamma)}&\lesssim
N^{-\frac{1}{4}-s}\|\lambda\circ\mathbf{r}\|_{Q^{s}_{e}(\mathbb{D})}\\\
&+\widetilde{\rho}^{2N}\left[(2\rho)^{-2N_{q}^{\theta}}+(2\rho)^{-2N_{q}})\right]N^{4}\|\lambda\|_{\widetilde{H}^{-\frac{1}{2}}(\Gamma)}\\\
&+\widetilde{\rho}^{N}\left[(2\rho)^{-2N_{q}^{\theta}}+(2\rho)^{-2N_{q}})\right]N^{2},\end{split}$
(29)

where $\lambda$ denotes the continuous solution of the Dirichlet BIE (12).

###### Proof.

Existence and uniqueness are obtained following [30, Theorem 4.2.11]. Moreover
from the same reference, it holds that

$\displaystyle\begin{split}\|\lambda-\lambda^{N,q}_{e}\|_{\widetilde{H}^{-\frac{1}{2}}(\Gamma)}\lesssim\inf_{v^{N}\in\mathbb{Q}^{e}_{N}(\Gamma)}\left(\|\lambda-v^{N}\|_{\widetilde{H}^{-\frac{1}{2}}(\Gamma)}+c(N_{q},N)\|v^{N}\|_{\widetilde{H}^{-\frac{1}{2}}(\Gamma)}\right)\\\
+\sup_{v^{N}\in\mathbb{Q}^{e}_{N}(\Gamma)}\frac{\boldsymbol{v}^{N}\cdot(\mathbf{g}_{d}^{N}-\mathbf{g}_{d}^{N,q})}{\|v^{N}\|_{\widetilde{H}^{-\frac{1}{2}}(\Gamma)}},\end{split}$

where
$c(N_{q},N){:=}\widetilde{\rho}^{2N}\left[(2\rho)^{-2N_{q}^{\theta}}+(2\rho)^{-2N_{q}})\right]N^{4}$.
The second term on the right-hand side can be estimated following the same
ideas of Lemma 8, and we obtain:

$\displaystyle\sup_{v^{N}\in\mathbb{Q}^{e}_{N}(\Gamma)}\frac{\boldsymbol{v}^{N}\cdot(\mathbf{g}_{d}^{N}-\mathbf{g}_{d}^{N,q})}{\|v^{N}\|_{\widetilde{H}^{-\frac{1}{2}}(\Gamma)}}\lesssim\widetilde{\rho}^{N}\left[(2\rho)^{-2N_{q}^{\theta}}+(2\rho)^{-2N_{q}})\right]N^{2}.$

With this, the result is obtained following Theorem 3. ∎

######  Remark 5.

In the right-hand side of (29), the last term grows as
$\widetilde{\rho}^{N}N^{2}$ due a to single numeric quadrature whereas the
second term does at rate of $\widetilde{\rho}^{2N}N^{4}$ as it arises from the
tensorization of two quadrature rules.

In a similar fashion, we obtain the equivalent result for the Neumann problem.

######  Theorem 6.

There exists $N_{0}>0$ such that for every $N>N_{0}$, there is a $N_{q,0}$
that depends on $N$ such that, if $N_{q}^{\theta}$, $N_{q}$ is both greater
than $N_{q,0}$, the fully discrete Neumann Problem 4 has a unique solution,
and the following error estimate follows for $s>\frac{1}{2}$:

$\displaystyle\begin{split}\|\nu-\nu^{N,q}_{o}\|_{\widetilde{H}^{\frac{1}{2}}(\Gamma)}\lesssim&\
N^{\frac{1}{2}-s}\|\nu\circ\mathbf{r}\|_{P^{s}_{o}(\mathbb{D})}+\widetilde{\rho}^{2N}\left[(2\rho)^{-2N_{q}^{\theta}}+(2\rho)^{-2N_{q}})\right]N^{2}\|\nu\|_{\widetilde{H}^{\frac{1}{2}}(\Gamma)}\\\
&+\widetilde{\rho}^{N}\left[(2\rho)^{-2N_{q}^{\theta}}+(2\rho)^{-2N_{q}})\right]N.\end{split}$
(30)

######  Remark 6.

The differences between the second (resp. third) term in the right hand-sides
of (29) and (30) is caused for the norm used to measure the error in each
case.

######  Remark 7.

From the fully discrete error analysis, one concludes that the number of
quadrature points should be a linear function of the parameter $N$ so as to
obtain good approximations of the BIOs.

##  6 Numerical Results

In what follows, we conduct a series of numerical experiments to verify our
claims, showcase insights and show limitations of the provided results. These
computational results were carried out on a desktop PC I7-4790k with 8Gb of
RAM.

###  6.1 Quadrature Results

(a)  $I^{a},I^{b},I^{c}$

(b)  $I^{d},I^{e}$

Figure 2:  Quadrature errors for $k=2.8$, computed against an overkill with
$N_{q}=76$. The error is computed by taking the 2-norm
($\|A\|_{2}=\sup_{\mathbf{v}\neq 0}\frac{\mathbf{v}\cdot
A\mathbf{v}}{\|\mathbf{v}\|}$) we fix $N=8$, (45 degrees of freedom).

Before studying the accuracy of the full spectral Galerkin method, we consider
the performance of the quadrature procedure detailed in Section 4.2.2. To this
end, we consider the screen $\Gamma$ given by the following parametrization:

$\mathbf{r}(r,\theta)=r(\cos\theta,2.8\sin\theta,-0.56r).$

We compute the integrals $I^{\alpha}$, $\alpha\in\\{a,b,c,d,e\\}$ for an
increasing number of quadrature points. In particular, we only select the
variable $N_{q}^{\alpha}$ and fix $N_{q}^{\theta}=N_{q}^{\alpha}+12$. Results
reported in Figure 2 show that quadrature errors decay linearly in the log-
linear scale –exponential decay– as described in (26).

As a second test, we consider the same screen and wavenumber for an increasing
number of functions and quadrature points. As before, we only modify the
variable $N_{q}^{\alpha}$, $\alpha\in{a,b,c,d,e,f}$ and fix
$N_{q}^{\theta}=N_{q}^{\alpha}+12$. The results are presented in Table 1. In
contrast to the previous experiment, we show the error for the consolidated
variable $V^{N}[k]-V^{N,q}[k]$ and observe that the quadrature error is almost
constant when the quadrature points increase linearly with $N$ as stated in
Remark 7.

Table 1: Quadrature error for the weakly-singular operator with $k=2.8$ computed against an overkill with $N_{q}=92,N=24$. The rule for increasing the quadrature points with $N$ is $N_{q}(N)=1.75N+50$. The error again is computed as the 2-norm of the approximation of the weakly-singular operator matrix. $N$ | 2 | 6 | 14 | 16 | 18 | 20
---|---|---|---|---|---|---
$N_{q}$ | 18 | 25 | 39 | 42 | 46 | 49
Error | 6.20e-15 | 4.89e-15 | 8.88e-15 | 9.60e-15 | 7.14e-15 | 8.39e-15

###  6.2 Code Validation

We now show that our method is correctly implemented for the case of the
Laplace Dirichlet and Neumann problem for the disk $\mathbb{D}$. Recall the
closed form of the matrix entries in (32). We consider $g_{d}$ (resp. $g_{n}$)
as the Dirichlet (resp. Neumann) traces of a plane wave:

$\displaystyle\exp(\imath k_{0}\mathbf{x}\cdot\mathbf{d})$

where $\mathbf{d}$ is an unitary vector which can be characterized in terms of
two angles:

$\mathbf{d}=(\cos\theta_{0}\cos\varphi_{0},\sin\theta_{0}\cos\varphi_{0},\sin\varphi_{0}),$

and $k_{0}$ is the wavenumber of the plane wave. In Figure 3(b) we show the
error of the approximated solution with respect to the semi-analytic solution
–the one obtained using (32) for the matrix computations and a fixed value of
$N$. Notice that the super-algebraic convergence rate is achieved, i.e. linear
behavior in log-linear scale.

(a)  $\Gamma_{1}$

(b)  right hand-side

Figure 3: (a) Convergence curves for the Laplace problems on the disk.
Parameters for the plane wave are $k_{0}=2.8$, $\theta_{0}=\pi/3$,
$\varphi_{0}=\frac{\pi}{4}$. Running times are for $N=14$ are 181s, 730s for
Dirichlet and Neumann cases, respectively. (b) Maximum value of absolute
values in the right-hand side for the corresponding $N$.

As in the two-dimensional case [19], Figure 3(b) suggests that one could
deduce the decaying behavior of solution coefficients from the right-hand side
coefficients. For the Laplace case on a disk, this comes directly from (32)
but for more general cases, to the best of our knowledge, this has not been
done. One could prove results in this context by establishing a complete
theory of pseudo-differential operators on screens acting in the auxiliary
spaces, in a similar fashion [1] presented for the two-dimensional case.

###  6.3 More complex screens

We fix the wavenumbers $k=2.8$, and for right-hand side we use a plane wave
with $k_{0}=k$, $\theta_{0}=\frac{\pi}{3}$, $\varphi_{0}=\frac{\pi}{4}$. Let
us first consider two distinct geometries: $\Gamma_{1}$ truncated elliptic
screen given by

$\mathbf{r}_{1}(r,\theta)=r(\cos\theta,2.8\sin\theta,0)$

and $\Gamma_{2}$ a truncated paraboloid,

$\mathbf{r}_{2}(r,\theta)=r(\cos\theta,2.8\sin\theta,-0.56r).$

Results for $\Gamma_{1},\Gamma_{2}$ are presented in Figure 4, wherein we
obtain super-algebraic convergence –linear convergence in log-linear scale– as
stated in Theorems 5, 6.

(a)  $\Gamma_{1}$

(b)  $\Gamma_{2}$

Figure 4:  Error in the $Q^{-\frac{1}{2}}_{e}$-Norm for the Dirichlet
problems, and $P_{o}^{\frac{1}{2}}$ for the Neumann. Overkill solutions are
computed using $N=20$. Run times for Dirichlet and Neumann problem (with
$N=20$) are 749s, 2087s respectively.

Next, we show the impact of critical points through the following screens.
$\Gamma_{3}$, described by

$\mathbf{r}_{3}(r,\theta)=r(1-0.2r^{3}\cos(3\theta))(\cos\theta,\sin\theta,0),$

and $\Gamma_{4}$, with parametrization

$\mathbf{r}_{4}(r,\theta)=r(1-0.3r^{3}\cos(3\theta))(\cos\theta,\sin\theta,0).$

The case $k=2.8$ is depicted in Figure 5. While these last two screens are
similar, one can see that the error convergence for $\Gamma_{4}$ is worse than
for the other cases. This is explained by the Jacobians’ behavior: for
$\Gamma_{3}$, one has

$|r(1-0.2r^{3}\cos(3\theta))(1-0.8r^{3}\cos\theta)|,$

which is no-where null, whereas for $\Gamma_{4}$

$|r(1-0.3r^{3}\cos(3\theta))(1-1.2r^{3}\cos\theta)|,$

is zero on the curve $1-1.2r^{3}\cos\theta=0$.

(a)  $\Gamma_{3}$

(b)  $\Gamma_{4}$

Figure 5:  Error in the $Q^{-\frac{1}{2}}_{e}$-Norm for the Dirichlet
problems, and $P_{o}^{\frac{1}{2}}$ for the Neumann. Overkill solutions are
computed using $N=28$. Run times for Dirichlet and Neumann problem (with
$N=28$) are 2800s, 6198s respectively.

Finally, and for illustration purposes, we consider a highly complex screen,
$\Gamma_{5}$, along with error convergence and volume solution plots (see
Figures 6(a)–(d)). The screen is given by the parametrization:

$\displaystyle\mathbf{r}_{5}(r,\theta)=(x_{5}(r,\theta),y_{5}(r,\theta),z_{5}(r,\theta)),$

where
$(x_{5}(r,\theta),y_{5}(r,\theta)):=r(1-0.2r^{3}\cos(3\theta))(\cos\theta,\sin\theta)$
and $z_{5}(r,\theta):=0.2r\cos(7y_{5}(r,\theta))$.

(a)  $\Gamma_{5}$

(b)  $\Gamma_{5}$

(c)  Convergence

(d)  Volume solution

Figure 6:  (a)-(b) Screen $\Gamma_{5}$ profiles. (c) Convergence plots for
Dirichlet and Neumann problems computed against an overkill solution computed
with $N=28$, for $k=2.8$, $k_{0}=k$, $\theta_{0}=\frac{\pi}{2}$, and
$\phi_{0}=0$. (d) Plot of the volume solution for the Dirichlet problem.

##  7 Concluding Remarks

The present work presents a high-order discretization method for open screens
based on weighted projections of the spherical harmonics functions. We have
proved that the method converges super-algebraically and also described
implementation details. As an efficient solver for the forward problem, future
efforts are directed towards using our method to solve shape optimization or
inverse problems.

While the analysis presented here describes how approximation errors behave as
a function of the number of degrees of freedom, it does not provide any
insight on how to choose this number in terms of the problem parameters, in
particular, the wavenumber $k$. We refer to [7] for a general review on this
topic and [8] for results related to screens.

## References

  * [1] François Alouges and Martin Averseng.  New preconditioners for laplace and helmholtz integral equations on open curves: Analytical framework and numerical results.  2019\.
  * [2] George B. Arfken, Hans J. Weber, and Frank E. Harris.  Chapter 15 - Legendre functions.  In George B. Arfken, Hans J. Weber, and Frank E. Harris, editors, Mathematical Methods for Physicists (Seventh Edition), pages 715 – 772\. Academic Press, Boston, seventh edition edition, 2013.
  * [3] K. Atkinson and W. Han.  Theoretical Numerical Analysis: A Functional Analysis Framework.  Texts in Applied Mathematics. Springer New York, 2001.
  * [4] M. Bebendorf.  Hierarchical Matrices: A Means to Efficiently Solve Elliptic Boundary Value Problems.  Springer Science & Business Media, 2008.
  * [5] Alexei Bespalov and Norbert Heuer.  The p-version of the boundary element method for weakly singular operators on piecewise plane open surfaces.  Numerische Mathematik, 106:69–97, 03 2007.
  * [6] Oscar P. Bruno and Stéphane K. Lintner.  A high-order integral solver for scalar problems of diffraction by screens and apertures in three-dimensional space.  Journal of Computational Physics, 252:250 – 274, 2013.
  * [7] Simon Chandler-Wilde, Ivan Graham, Stephen Langdon, and Euan Spence.  Numerical-asymptotic boundary integral methods in high-frequency acoustic scattering.  Acta Numerica, 21:89 – 305, 04 2012.
  * [8] Simon Chandler-Wilde and David Hewett.  Wavenumber-explicit continuity and coercivity estimates in acoustic scattering by planar screens.  Integral Equations and Operator Theory, 82:423–449, 07 2015.
  * [9] Martin Costabel and Monique Dauge.  Crack singularities for general elliptic systems.  Mathematische Nachrichten, 235(1):29–49, 2002.
  * [10] Martin Costabel, Monique Dauge, and Roland Duduchava.  Asymptotics without logarithmic terms for crack problems.  Communications in Partial Differential Equations, 28(5-6):869–926, 2003.
  * [11] Walter Gautschi.  Some elementary inequalities relating to the gamma and incomplete gamma function.  Journal of Mathematics and Physics, 38(1-4):77–81, 1959.
  * [12] I.G. Graham and Ian Sloan.  Fully discrete spectral boundary integral methods for Helmholtz problems on smooth closed surfaces in $\mathbb{R}^{3}$.  Numerische Mathematik, 92:289–323, 01 2002.
  * [13] Tuong Ha-Duong.  On the transient acoustic scattering by a flat object.  Japan Journal of Applied Mathematics, 7, 10 1990.
  * [14] Norbert Heuer, Matthias Maischak, and Ernst Stephan.  Exponential convergence of the hp-version for the boundary element method on open surfaces.  Numerische Mathematik, 83, 02 1970.
  * [15] Norbert Heuer, Mario E. Mellado, and Ernst P. Stephan.  A p-adaptive algorithm for the bem with the hypersingular operator on the plane screen.  International Journal for Numerical Methods in Engineering, 53(1):85–104, 2002.
  * [16] R. Hiptmair, C. Jerez-Hanckes, and C. Urzúa-Torres.  Closed-form exact inverses of the weakly singular and hypersingular operators on disks.  Integral Equations and Operator Theory, 90(4):1–14, 2018.
  * [17] R. Hiptmair, C. Jerez-Hanckes, and C. Urzúa-Torres.  Optimal operator preconditioning for Galerkin boundary element methods on 3-dimensional screens.  SIAM Journal on Numerical Analysis, 58(1):834–857, 2020.
  * [18] Carlos Jerez-Hanckes and Jean-Claude Nédélec.  Explicit variational forms for the inverses of integral logarithmic operators over an interval.  SIAM Journal on Mathematical Analysis, 44(4):2666–2694, jan 2012\.
  * [19] Carlos Jerez-Hanckes and José Pinto.  High-order Galerkin method for Helmholtz and Laplace problems on multiple open arcs.  ESAIM: M2AN, 54(6):1975–2009, 2020.
  * [20] Carlos Jerez-Hanckes and José Pinto.  Spectral Galerkin method for solving Helmholtz and Laplace Dirichlet problems on multiple open arcs.  In Spencer J. Sherwin, David Moxey, Joaquim Peiró, Peter E. Vincent, and Christoph Schwab, editors, Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2018, pages 383–393, Cham, 2020\. Springer International Publishing.
  * [21] Rainer Kress.  Inverse scattering from an open arc.  Mathematical Methods in the Applied Sciences, 18(4):267–293, 1995\.
  * [22] Stéphane K. Lintner and Oscar P. Bruno.  A generalized Calderón formula for open-arc diffraction problems: theoretical considerations.  Proceedings of the Royal Society of Edinburgh: Section A Mathematics, 145(2):331–364, 2015.
  * [23] G Lohöfer.  Inequalities for Legendre functions and Gegenbauer functions.  Journal of Approximation Theory, 64(2):226 – 234, 1991.
  * [24] T.M. MacRobert.  Spherical Harmonics: An Elementary Treatise on Harmonic Functions, with Applications.  Chronic illness in the United States. Dover Publications, 1948.
  * [25] W. McLean.  Strongly Elliptic Systems and Boundary Integral Equations.  Cambridge University Press, 2000.
  * [26] Duong Pham, T. Tran, and Alexey Chernov.  Pseudodifferential equations on the sphere with spherical splines.  Mathematical Models and Methods in Applied Sciences, 21:1933–1959, 09 2011.
  * [27] Pedro Ramaciotti and Jean-Claude Nédélec.  About some boundary integral operators on the unit disk related to the Laplace equation.  SIAM J. Numer. Anal., 55(4):1892–1914, 2017.
  * [28] Pedro Ramaciotti Morales.  Theoretical and numerical aspects of wave propagation phenomena in complex domains and applications to remote sensing.  PhD thesis, 10 2016.
  * [29] J. Saranen and G. Vainikko.  Periodic Integral and Pseudodifferential Equations with Numerical Approximation.  Springer Monographs in Mathematics. Springer Berlin Heidelberg, 2013.
  * [30] S. A. Sauter and Ch. Schwab.  Boundary Element Methods, volume 39.  Springer Series in Computational Mathematics, 2011.
  * [31] Nathanaël Schaeffer.  Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations.  Geochemistry, Geophysics, Geosystems, 14(3):751–758, 2013.
  * [32] E. P. Stephan.  A boundary integral equation method for three-dimensional crack problems in elasticity.  Math. Methods Appl. Sci., 8(4):609–623, 1986.
  * [33] Ernst P. Stephan.  Boundary integral equations for magnetic screens in $\mathbb{R}^{3}$.  Proceedings of the Royal Society of Edinburgh: Section A Mathematics, 102:189 – 210, 01 1986.
  * [34] Ernst P. Stephan.  Boundary integral equations for screen problems in ${\mathbb{R}}^{3}$.  Integral Equations and Operator Theory, 10(2):236–257, Mar 1987\.
  * [35] J. Stoer and R. Bulirsch.  Introduction to numerical analysis.  Springer-Verlag, 1980.
  * [36] L.N. Trefethen.  Approximation Theory and Approximation Practice.  Other Titles in Applied Mathematics. SIAM, 2013.
  * [37] Carolina Urzúa-Torres.  Operator Preconditioning for Galerkin Boundary Element Methods on Screens.  PhD thesis, 2018.
  * [38] Haiyong Wang and Lun Zhang.  Jacobi polynomials on the Bernstein ellipse.  J. Sci. Comput., 75(1):457–477, April 2018.
  * [39] P. Wolfe.  Eigenfunctions of the Integral Equation for the Potential of the Charged Disk.  Journal of Mathematical Physics, 12:1215–1218, July 1971.
  * [40] L. Yijun.  Fast Multipole Boundary Element Method: Theory and Applications in Engineering.  Cambridge University, 2009.

##  Appendix A Proof of Lemma 3

Recall the weakly- and hyper-singular BIOs defined in Section 3.1. For $k=0$
and $\Gamma=\mathbb{D}$, we write

$\displaystyle V_{\mathbb{D}}u(\mathbf{x})$
$\displaystyle:=\int_{\mathbb{D}}\frac{1}{4\pi\|\mathbf{x}-\mathbf{x}^{\prime}\|_{2}}u(\mathbf{x}^{\prime})d\mathbf{x}^{\prime},$
$\displaystyle W_{\mathbb{D}}u(\mathbf{x})$
$\displaystyle:=-\gamma_{n,\mathbf{x}}\int_{\mathbb{D}}\gamma_{n,\mathbf{x}^{\prime}}\frac{1}{4\pi\|\mathbf{x}-\mathbf{x}^{\prime}\|_{2}}u(\mathbf{x}^{\prime})d\mathbf{x}^{\prime},$

These operators have two key properties. First, they are continuous and
elliptic so they can be used to define equivalent norms. In fact, by [30,
Theorem 3.5.9] we have that

$\displaystyle\|u\|^{2}_{\widetilde{H}^{-\frac{1}{2}}(\mathbb{D})}\cong\langle
V_{\mathbb{D}}u,u\rangle_{\mathbb{D}},\quad\|u\|^{2}_{\widetilde{H}^{\frac{1}{2}}(\mathbb{D})}\cong\langle
W_{\mathbb{D}}u,u\rangle_{\mathbb{D}}.$ (31)

Secondly, we have a characterization of the eigenvalues of these two
operators111111See [39] for the original proof of the weakly-singular BIO and
[27, Theorem 2.7.1] for the hyper-singular case:

$\displaystyle V_{\mathbb{D}}q^{l}_{m}=\frac{1}{4}\Lambda_{l,m}p^{l}_{m},\quad
l+m\text{ even},\quad
W_{\mathbb{D}}p^{l}_{m}=\frac{1}{\Lambda_{l,m}}q^{l}_{m},\quad l+m\text{
odd},$ (32)

where

$\displaystyle\Lambda_{l,m}=\frac{\Gamma\left(\frac{l+|m|+1}{2}\right)\Gamma\left(\frac{l-|m|+1}{2}\right)}{\Gamma\left(\frac{l+|m|+2}{2}\right)\Gamma\left(\frac{l-|m|+2}{2}\right)},$

here $\Gamma$ denotes the Gamma function and not an screen. Using Gautschi’s
inequality [11], it holds that

$\displaystyle\frac{1}{l+1}\lesssim\Lambda_{l,m}\lesssim\frac{1}{\sqrt{l+1}}.$
(33)

With these elements, we proceed with the proof of Lemma 3. Pick any smooth
function $u$ on $\mathbb{D}$ and consider its even lifting to $\mathbb{S}$.
Since spherical harmonics are dense on smooth functions defined on the sphere,
we can approximate the lifting $u$ by even spherical harmonics. Thus, $u$ can
be expanded as

$u=\sum_{l=0}^{\infty}\sum_{\begin{subarray}{c}m=-l\\\ m+l\text{
even}\end{subarray}}^{l}\langle u,q_{m}^{l}\rangle p_{m}^{l},$

now (3) follows from the orthogonality relation (1). To prove (4), we use the
density of functions $q_{l}^{e}$ in
$\widetilde{H}^{-\frac{1}{2}}(\mathbb{D})$, i.e.

$u=\sum_{l=0}^{\infty}\sum_{\begin{subarray}{c}m=-l\\\ m+l\text{
even}\end{subarray}}^{l}u^{l}_{m}q^{l}_{m},$

wherein, by orthogonality it holds that $u^{l}_{m}=\langle
u,p^{l}_{m}\rangle$. Hence, computing the norm of $u$ using the equivalence
(31), the relation (32) and the estimate (33) yields

$\displaystyle\|u\|_{Q_{e}^{-\frac{1}{2}}(\mathbb{D})}\lesssim\|u\|_{\widetilde{H}^{-\frac{1}{2}}(\mathbb{D})}\lesssim\|u\|_{Q_{e}^{-\frac{1}{4}}(\mathbb{D})},$

which implies the result. Similar ideas are used to show (5). ∎

##  Appendix B Singular Integrals Analysis

We consider the integrals $I^{c},I^{d},I^{e}$ defined in Section 4.2.2. These
integrals have in common that the singularities occur in specifics points in
three- or two-dimensional spaces, i.e. when three or two variables take a
specific value. In contrast, $I^{a},I^{b}$, singularities occur when one
variable takes a specific value regardless of the other.

###  B.1 General idea

Let us start with the simpler case of an integral which has a singularity in
2D:

$\displaystyle I:=\int_{0}^{1}\int_{0}^{1}\frac{1}{\sqrt{x+y}}dydx,$

the integrand has a singularity at $(x,y)=(0,0)$. Performing the polar change
of variables $x=\rho\cos\alpha$, $y=\rho\sin\alpha$, we have that

$\displaystyle\begin{split}I=\int_{0}^{\frac{\pi}{4}}\int_{0}^{\frac{1}{\cos\alpha}}\frac{\sqrt{\rho}}{\sqrt{\cos\alpha+\sin\alpha}}d\rho
d\alpha+\int_{\frac{\pi}{4}}^{\frac{\pi}{2}}\int_{0}^{\frac{1}{\sin\alpha}}\frac{\sqrt{\rho}}{\sqrt{\cos\alpha+\sin\alpha}}d\rho
d\alpha.\end{split}$

Moreover, one can fix the integration domain for the $\rho$ variable by doing
a linear change of variable

$\displaystyle\begin{split}I=\int_{0}^{\frac{\pi}{4}}\int_{0}^{1}\frac{\sqrt{t}}{(\cos\alpha)^{3/2}\sqrt{\cos\alpha+\sin\alpha}}dtd\alpha+\\\
\int_{\frac{\pi}{4}}^{\frac{\pi}{2}}\int_{0}^{1}\frac{\sqrt{t}}{(\sin\alpha)^{3/2}\sqrt{\cos\alpha+\sin\alpha}}dtd\alpha.\end{split}$

These last two integrals can be straightforwardly computed by using a Jacobi
rule in $t$ and Gauss-Legendre in $\alpha$, resulting in an optimal
convergence rate –exponential in this particular case. The idea is in fact
very simple: use polar coordinates with the origin in the point where the
singularity occurs, transferring the multidimensional singularity to the
radial coordinate only. The rest of this section gives the detail on each
change of variable that is needed and also proving that the resulting
integrals are computed with optimal rates.

###  B.2 Integral $I^{c}$

Consider the integral of the form

$\displaystyle\begin{split}J^{c}:=\int_{\frac{\sqrt{3}}{2}}^{1}\int_{\frac{-\pi}{2}}^{\frac{-\pi}{3}}\int_{0}^{\frac{1}{2}}\frac{\lambda
A^{2}(r,\beta)}{4\pi\|\mathbf{r}(\mathbf{x})-\mathbf{r}(\mathbf{x}+\lambda
A(r,\beta)\mathbf{e}_{\theta+\beta})\|}\\\ \frac{rd\lambda d\beta
dr}{\sqrt{1-r^{2}}\sqrt{1-\|\mathbf{x}+\lambda
A(r,\beta)\mathbf{e}_{\theta+\beta}\|^{2}}},\end{split}$ (34)

where in comparison to (19), we restrict to the first interval for the $\beta$
variable as the other cases follows similarly, and we have omitted the
integral in the $\theta$ variable as it is not relevant to the singularity
analysis and can be treated with Gauss-Legendre quadrature having not effect
on the rate of convergence. Since trial and test polynomials are smooth
functions they have also been neglected in the ensuing singularity analysis.
Consider the following change of variables: $u^{2}=1-r^{2}$, $\cos\beta=v$.
Define $d:=4\pi\|\mathbf{r}(\mathbf{x})-\mathbf{r}(\mathbf{x}+\lambda
A(r,\beta)\mathbf{e}_{\theta+\beta})\|$ and use expansion (18) so that (34)
becomes

$\displaystyle
J^{c}=\int_{0}^{\frac{1}{2}}\int_{0}^{\frac{1}{2}}\int_{0}^{\frac{1}{2}}\frac{\lambda
d^{-1}A^{2}dudvd\lambda}{\sqrt{1-\lambda}\sqrt{1-v^{2}}\sqrt{u^{2}(1+\lambda)-2\lambda\sqrt{1-u^{2}}vA}},$

wherein by definition of $A(r,\beta)$ (see (17)), we have that

$\displaystyle A=\left(\sqrt{u^{2}+(1-u^{2})v^{2}}-\sqrt{1-u^{2}}v\right).$

Apply the first polar change of variables $u=\rho\cos\alpha$,
$v=\rho\sin\alpha$ so that

$\displaystyle
A=\rho\left(\sqrt{1-\rho^{2}\cos^{2}\alpha\sin^{2}\alpha}-\sin\alpha\sqrt{1-\rho^{2}\cos^{2}\alpha}\right)=:\rho\widetilde{A},$

and we obtain two integrals

$\displaystyle
J^{c}_{1}:=\int_{0}^{\frac{1}{2}}\frac{\lambda}{\sqrt{1-\lambda}}\int_{0}^{\frac{\pi}{4}}\int_{0}^{\frac{1}{2\cos\alpha}}\frac{(1-\rho^{2}\sin^{2}\alpha)^{-1}\rho^{2}d^{-1}\widetilde{A}^{2}d\rho
d\alpha
d\lambda}{\sqrt{\cos^{2}\alpha(1+\lambda)-2\lambda\sqrt{1-\rho^{2}\cos^{2}\alpha}\sin\alpha\widetilde{A}}},$
(35) $\displaystyle
J^{c}_{2}:=\int_{0}^{\frac{1}{2}}\frac{\lambda}{\sqrt{1-\lambda}}\int_{\frac{\pi}{4}}^{\frac{\pi}{2}}\int_{0}^{\frac{1}{2\sin\alpha}}\frac{(1-\rho^{2}\sin^{2}\alpha)^{-1}\rho^{2}d^{-1}\widetilde{A}^{2}d\rho
d\alpha
d\lambda}{\sqrt{\cos^{2}\alpha(1+\lambda)-2\lambda\sqrt{1-\rho^{2}\cos^{2}\alpha}\sin\alpha\widetilde{A}}}.$
(36)

We apply a linear change of variables to fix the integration domain of the
$\rho$ variable. For the first integral, it holds that

$\displaystyle
J^{c}_{1}=\int_{0}^{\frac{1}{2}}\frac{\lambda}{\sqrt{1-\lambda}}\int_{0}^{\frac{\pi}{4}}\int_{0}^{\frac{1}{2}}\frac{t^{2}d^{-1}\widetilde{A}^{2}\cos\alpha^{-3}dtd\alpha
d\lambda}{\sqrt{1-t^{2}\tan^{2}\alpha}\sqrt{\cos^{2}\alpha(1+\lambda)-2\lambda\sqrt{1-t^{2}}\sin\alpha\widetilde{A}}},$

wherein
$\widetilde{A}=\left(\sqrt{1-t^{2}\sin^{2}\alpha}-\sin\alpha\sqrt{1-t^{2}}\right)$.
One can see that this integral converges at the optimal rate when we use the
Gauss-Legendre rule for all the variables as neither of the terms inside the
square roots vanish in the integration domain. For the second integral we have

$\displaystyle
J^{c}_{2}=\int_{0}^{\frac{1}{2}}\frac{\lambda}{\sqrt{1-\lambda}}\int_{\frac{\pi}{4}}^{\frac{\pi}{2}}\int_{0}^{\frac{1}{2}}\frac{t^{2}d^{-1}\widetilde{A}^{2}\sin\alpha^{-3}dtd\alpha
d\lambda}{\sqrt{1-t^{2}}\sqrt{\cos^{2}\alpha(1+\lambda)-2\lambda\sqrt{1-t^{2}\cot^{2}\alpha}\sin\alpha\widetilde{A}}},$

and where
$\widetilde{A}=\left(\sqrt{1-t^{2}\cos^{2}\alpha}-\sin\alpha\sqrt{1-t^{2}\cot^{2}\alpha}\right)$.
In contrast to $J^{c}_{1}$, one can easily verify that the integrated has one
singularity in $(\lambda,\alpha)=(0,\frac{\pi}{2})$, so further
transformations are needed. In particular, we use $z=\cos\alpha$,
$x^{2}=\lambda$. Thus, we have that

$\displaystyle
J^{c}_{2}=\int_{0}^{1/\sqrt{2}}\frac{2x^{3}}{\sqrt{1-x^{2}}}\int_{0}^{1/\sqrt{2}}\int_{0}^{\frac{1}{2}}\frac{t^{2}d^{-1}\widetilde{A}^{2}(1-z^{2})^{-2}dtdzdx}{\sqrt{1-t^{2}}\sqrt{z^{2}(1+x^{2})-2x^{2}\sqrt{1-z^{2}-t^{2}z^{2}}\widetilde{A}}},$

with
$\widetilde{A}=\left(\sqrt{1-t^{2}z^{2}}-\sqrt{1-z^{2}-t^{2}z^{2}}\right).$
Once again we make a polar change of variables $x=\sigma\cos(\phi),$
$z=\sigma\sin(\phi)$, and we obtain two integrals

$\displaystyle\begin{split}J^{c}_{2,1}&:=\int_{0}^{\frac{1}{2}}\int_{0}^{\frac{\pi}{4}}\int_{0}^{1/(\sqrt{2}\cos\phi)}\frac{2(\sigma\cos\phi)^{3}t^{2}d^{-1}\widetilde{A}^{2}(1-(\sigma\sin\phi)^{2})^{-2}}{\sqrt{1-(\sigma\cos\phi)^{2}}\sqrt{1-t^{2}}}\\\
&\times\frac{d\sigma d\phi
dt}{\sqrt{\sin^{2}\phi(1+(\sigma\cos\phi)^{2})-2\cos^{2}\phi\sqrt{1-(\sigma\sin\phi)^{2}(1+t^{2})}\widetilde{A}}},\\\
\end{split}$
$\displaystyle\begin{split}J^{c}_{2,2}&:=\int_{0}^{\frac{1}{2}}\int_{\frac{\pi}{4}}^{\frac{\pi}{2}}\int_{0}^{1/(\sqrt{2}\sin\phi)}\frac{2(\sigma\cos\phi)^{3}t^{2}d^{-1}\widetilde{A}^{2}(1-(\sigma\sin\phi)^{2})^{-2}}{\sqrt{1-(\sigma\cos\phi)^{2}}\sqrt{1-t^{2}}}\\\
&\times\frac{d\sigma d\phi
dt}{\sqrt{\sin^{2}\phi(1+(\sigma\cos\phi)^{2})-2\cos^{2}\phi\sqrt{1-(\sigma\sin\phi)^{2}(1+t^{2})}\widetilde{A}}},\\\
\end{split}$

where
$\widetilde{A}=\left(\sqrt{1-t^{2}(\sigma\sin\phi)^{2}}-\sqrt{1-(\sigma\sin\phi)^{2}(1+t^{2})}\right)$.
Finally, we take the corresponding change of variables needed to fix the
integration domain for the $\sigma$ variable, and we obtain

$\displaystyle\begin{split}J^{c}_{2,1}&=\int_{0}^{\frac{1}{2}}\int_{0}^{\frac{\pi}{4}}\int_{0}^{1/\sqrt{2}}\frac{2s^{3}t^{2}d^{-1}\widetilde{A}^{2}(1-(s\tan\phi)^{2})^{-2}}{\cos\phi\sqrt{1-s^{2}}\sqrt{1-t^{2}}}\\\
&\times\frac{dsd\phi
dt}{\sqrt{\sin^{2}\phi(1+s^{2})-2\cos^{2}\phi\sqrt{1-(s\tan^{2}\phi)(1+t^{2})}\widetilde{A}}},\\\
\end{split}$

with
$\widetilde{A}=\left(\sqrt{1-t^{2}(s\tan\phi))^{2}}-\sqrt{1-(s\tan\phi)^{2}(1+t^{2})}\right).$
It is easy to see that the integrand has not singularities and as so it can be
integrated with optimum rate using a tenzorisation of the Gauss-Legendre rule.
For the second integral we have

$\displaystyle\begin{split}J^{c}_{2,2}&=\int_{0}^{\frac{1}{2}}\int_{\frac{\pi}{4}}^{\frac{\pi}{2}}\int_{0}^{1/\sqrt{2}}\frac{2(s\cot\phi)^{3}t^{2}d^{-1}\widetilde{A}^{2}(1-s^{2})^{-2}}{\sin\phi\sqrt{1-(s\cot\phi)^{2}}\sqrt{1-t^{2}}}\\\
&\times\frac{dsd\phi
dt}{\sqrt{\sin^{2}\phi(1+(s\cot\phi)^{2})-2\cos^{2}\phi\sqrt{1-s^{2}(1+t^{2})}\widetilde{A}}},\\\
\end{split}$

where
$\widetilde{A}=\left(\sqrt{1-t^{2}s^{2}}-\sqrt{1-s^{2}(1+t^{2})}\right)$. Once
again the integrand does not has singularities and the optimum rate of
convergence is retrieved.

###  B.3 Integral $I^{d}$

As in the previous case we simplify the Integral $I^{d}$ in (20) to the
following integral

$\displaystyle
J^{d}:=\int_{\frac{\sqrt{3}}{2}}^{1}\int_{\frac{-\pi}{3}}^{\frac{\pi}{3}}\int_{0}^{\frac{1}{2}}\frac{\lambda
A^{2}(r,\beta)}{4\pi\|\mathbf{r}(\mathbf{x})-\mathbf{r}(\mathbf{x}+\lambda
A(r,\beta)\mathbf{e}_{\theta+\beta})\|}\frac{(1-r^{2})^{-\frac{1}{2}}rd\lambda
d\beta dr}{\sqrt{1-\|\mathbf{x}+\lambda
A(r,\beta)\mathbf{e}_{\theta+\beta}\|^{2}}},$ (37)

where we have only take the first interval for the $\beta$ variable as the
other case is similar. We have to make a transformation in $\lambda,r$: start
by fixing the singularity $(\lambda,r)=(0,1)$ to the origin of the new
variables $r^{2}=1-u^{2}$, $v^{2}=\lambda$, by doing so we obtain

$\displaystyle
J^{d}:=\int_{0}^{\frac{1}{2}}\int_{\frac{-\pi}{3}}^{\frac{\pi}{3}}\int_{0}^{1/\sqrt{2}}\frac{2v^{3}d^{-1}A^{2}dvd\beta
dr}{\sqrt{1-v^{2}}\sqrt{(u^{2}(1+v^{2})-2v^{2}\sqrt{1-u^{2}}\cos\beta A}},$

where $A=\sqrt{u^{2}+(1-u^{2})\cos^{2}\beta}-\sqrt{1-u^{2}}\cos\beta$. Now, we
make the polar change of variables $u=\rho\cos\alpha$, $v=\rho\sin\alpha$,
which leads to

$\displaystyle\begin{split}J^{d}_{1}&:=\int_{\frac{-\pi}{3}}^{\frac{\pi}{3}}\int_{0}^{\arctan(\sqrt{2})}\int_{0}^{1/(2\cos\alpha)}\frac{d^{-1}2(\rho\sin\alpha)^{3}A^{2}}{\sqrt{1-(\rho\sin\alpha)^{2}}}\\\
&\times\frac{d\rho d\alpha
d\beta}{\sqrt{\cos^{2}\alpha(1+(\rho\sin\alpha)^{2})-2\sin^{2}\alpha\sqrt{1-(\rho\cos\alpha)^{2}}\cos\beta
A}},\end{split}$
$\displaystyle\begin{split}J^{d}_{2}&:=\int_{\frac{-\pi}{3}}^{\frac{\pi}{3}}\int_{\arctan(\sqrt{2})}^{\frac{\pi}{2}}\int_{0}^{1/(\sqrt{2}\sin\alpha)}\frac{d^{-1}2(\rho\sin\alpha)^{3}A^{2}}{\sqrt{1-(\rho\sin\alpha)^{2}}}\\\
&\times\frac{d\rho d\alpha
d\beta}{\sqrt{\cos^{2}\alpha(1+(\rho\sin\alpha)^{2})-2\sin^{2}\alpha\sqrt{1-(\rho\cos\alpha)^{2}}\cos\beta
A}},\end{split}$

with
$A=\sqrt{(\rho\cos\alpha)^{2}+(1-(\rho\cos\alpha)^{2})\cos^{2}\beta}-\sqrt{1-(\rho\cos\alpha)^{2}}\cos\beta.$

Finally we apply the corresponding linear transformation to fix the $\rho$
integration domain,

$\displaystyle\begin{split}J^{d}_{1}&:=\int_{\frac{-\pi}{3}}^{\frac{\pi}{3}}\int_{0}^{\arctan(\sqrt{2})}\int_{0}^{\frac{1}{2}}\frac{d^{-1}2(t\tan\alpha)^{3}A^{2}}{\cos\alpha\sqrt{1-(t\tan\alpha)^{2}}}\\\
&\times\frac{dtd\alpha
d\beta}{\sqrt{\cos^{2}\alpha(1+(t\tan\alpha)^{2})-2\sin^{2}\alpha\sqrt{1-t^{2}}\cos\beta
A}},\end{split}$

where $A=\sqrt{t^{2}+(1-t^{2})\cos^{2}\beta}-\sqrt{1-t^{2}}\cos\beta$. This
integrand is smooth: only possible singularities can occur when $\rho=0$ but
are eliminated by the numerator, and so the rate of convergence is optimal.
Similarly for the second part as

$\displaystyle\begin{split}J^{d}_{2}&:=\int_{\frac{-\pi}{3}}^{\frac{\pi}{3}}\int_{\arctan(\sqrt{2})}^{\frac{\pi}{2}}\int_{0}^{1/\sqrt{2}}\frac{d^{-1}2t^{3}A^{2}}{\sin\alpha\sqrt{1-t^{2}}}\\\
&\times\frac{d\rho d\alpha
d\beta}{\sqrt{\cos^{2}\alpha(1+t^{2})-2\sin^{2}\alpha\sqrt{1-(t\cot\alpha)^{2}}\cos\beta
A}},\end{split}$

where the new change of variables leads to

$\displaystyle
A=\sqrt{(t\cot\alpha)^{2}+(1-(t\cot\alpha)^{2})\cos^{2}\beta}-\sqrt{1-(t\cot\alpha)^{2}}\cos\beta.$

Once again the integrand is smooth and the optimal convergence rate is
retrieved.

###  B.4 Integral $I^{e}$

Finally, we consider the simplification of $I^{e}$ defined as in (21):

$\displaystyle
J^{e}:=\int_{\frac{\sqrt{3}}{2}}^{1}\int_{\frac{-\pi}{2}}^{\frac{-\pi}{3}}\int_{\frac{1}{2}}^{1}\frac{\lambda
A^{2}(r,\beta)}{4\pi\|\mathbf{r}(\mathbf{x})-\mathbf{r}(\mathbf{x}+\lambda
A(r,\beta)\mathbf{e}_{\theta+\beta})\|}\frac{rd\lambda d\beta
dr}{\sqrt{1-r^{2}}\sqrt{1-\|\mathbf{x}+\lambda
A(r,\beta)\mathbf{e}_{\theta+\beta}\|^{2}}},$

we have again fixed the value of $\beta$ to the first interval as other cases
are similar. The analysis here follows that of the first part of $J_{1}$ (up
until the first polar change of variables). We obtain two integrals, the first
one being

$\displaystyle
J^{e}_{1}=\int_{\frac{1}{2}}^{1}\frac{\lambda}{\sqrt{1-\lambda}}\int_{0}^{\frac{\pi}{4}}\int_{0}^{\frac{1}{2}}\frac{t^{2}d^{-1}\widetilde{A}^{2}\cos^{-3}\alpha
dtd\alpha
d\lambda}{\sqrt{1-t^{2}\tan^{2}\alpha}\sqrt{\cos^{2}\alpha(1+\lambda)-2\lambda\sqrt{1-t^{2}}\sin\alpha\widetilde{A}}},$

with
$\widetilde{A}=\left(\sqrt{1-t^{2}\sin^{2}\alpha}-\sin\alpha\sqrt{1-t^{2}}\right)$,
and the second one

$\displaystyle
J^{e}_{2}=\int_{\frac{1}{2}}^{1}\frac{\lambda}{\sqrt{1-\lambda}}\int_{\frac{\pi}{4}}^{\frac{\pi}{2}}\int_{0}^{\frac{1}{2}}\frac{t^{2}d^{-1}\widetilde{A}^{2}\sin\alpha^{-3}dtd\alpha
d\lambda}{\sqrt{1-t^{2}}\sqrt{\cos^{2}\alpha(1+\lambda)-2\lambda\sqrt{1-t^{2}\cot^{2}\alpha}\sin\alpha\widetilde{A}}},$

where
$\widetilde{A}=\left(\sqrt{1-t^{2}\cos^{2}\alpha}-\sin\alpha\sqrt{1-t^{2}\cot^{2}\alpha}\right)$.
Contrary to the first case, these integrands are not smooth due to the term
$(1-\lambda)^{-1}$. However, by using a Jacobi rule in $\lambda$ one recovers
the optimal convergence rate.