# Error Analysis for the Implicit Boundary Integral Method

Yimin Zhong Department of Mathematics and Statistics, Auburn University,
Auburn AL, USA<EMAIL_ADDRESS>Kui Ren Department of Applied Physics
and Applied Mathematics, Columbia University, New York NY, USA
<EMAIL_ADDRESS>Olof Runborg Department of Mathematics, KTH Royal
Institute of Technology, Stockholm, Sweden<EMAIL_ADDRESS>Richard Tsai
Department of Mathematics and Oden Institute for Computational Engineering and
Sciences, The University of Texas at Austin, Austin TX, USA
<EMAIL_ADDRESS>
###### Abstract

The implicit boundary integral method (IBIM) provides a framework to construct
quadrature rules on regular lattices for integrals over irregular domain
boundaries. This work provides a systematic error analysis for IBIMs on
uniform Cartesian grids for boundaries with different degree of regularities.
We first show that the quadrature error gains an addition order of
$\frac{d-1}{2}$ from the curvature for a strongly convex smooth boundary due
to the “randomness” in the signed distances. This gain is discounted for
degenerated convex surfaces. We then extend the error estimate to general
boundaries under some special circumstances, including how quadrature error
depends on the boundary’s local geometry relative to the underlying grid.
Bounds on the variance of the quadrature error under random shifts and
rotations of the lattices are also derived.

Key words. implicit boundary integral method, error analysis, level set,
solvent-excluded surface

##  1 Introduction

Let $\Omega\subset\mathbb{R}^{d}$ be a general simply connected domain with
$C^{\infty}$ boundary $\Gamma$. Without loss of generality, we assume that the
origin ${\mathbf{0}}\in\Omega$ and define the $\varepsilon$-tube
$T_{\varepsilon}$:

$T_{\varepsilon}:=\\{x\in\Omega\mid
0\leq\mathrm{dist}(\mathbf{x},\Gamma)\leq\varepsilon\\},$

where
$\mathrm{dist}(\mathbf{x},\Gamma)=\min_{\mathbf{y}\in\Gamma}|\mathbf{y}-\mathbf{x}|$
is the _unsigned_ distance function to the boundary. We make further
assumption that $\forall\mathbf{x}\in T_{\varepsilon}$, the projection
$P_{\Gamma}(\mathbf{x}):T_{\varepsilon}\mapsto\Gamma$:

$P_{\Gamma}(\mathbf{x})=\operatorname*{\arg\min}_{\mathbf{y}\in\Gamma}|\mathbf{y}-\mathbf{x}|$

is well-defined and the _signed_ distance function $d_{\Gamma}$ is defined by

$d_{\Gamma}(\mathbf{x}):=\mathbf{n}(P_{\Gamma}(\mathbf{x}))\cdot(\mathbf{x}-P_{\Gamma}(\mathbf{x})),$

where $\mathbf{n}$ denotes the outward unit normal vector on $\Gamma$. Let
$\sigma$ be the boundary Lebesgue measure defined on $\Gamma$. We can then use
the co-area formula to rewrite the following boundary integral for
$f(\mathbf{x})\in C^{\infty}(\Gamma)$

$\mathcal{I}(f):=\int_{\Gamma}f(\mathbf{x})d\sigma(\mathbf{x})$

into a volumetric integral in $T_{\varepsilon}$:

$\mathcal{I}(f)=\int_{T_{\varepsilon}}f(P_{\Gamma}(\mathbf{x}))\theta_{\varepsilon}(d_{\Gamma}(\mathbf{x}))J_{\varepsilon}(\mathbf{x},d_{\Gamma}(\mathbf{x}))d\mathbf{x}.$
(1)

The weight function
$\theta_{\varepsilon}(s):=\varepsilon^{-1}\theta(\varepsilon^{-1}s)$ is a
regularized 1D Dirac function such that

$\operatorname{supp}\theta=[-1,1]\;\text{ and }\;\int_{-1}^{1}\theta(s)ds=1.$

The function $J_{\varepsilon}(\cdot,\eta)$ is the Jacobian for the projection
$P_{\Gamma}$ on the level set surface $\\{d_{\Gamma}=\eta\\}$. For $d=2$,
$J_{\varepsilon}(\mathbf{x},\eta)=1-\eta\kappa(\mathbf{x})$ with $\kappa$ the
signed curvature of the level curve $\Gamma_{\eta}:=\\{\mathbf{x}\mid
d_{\Gamma}(\mathbf{x})=\eta\\}$. For $d=3$,
$J_{\varepsilon}(\mathbf{x},\eta)=1-2\eta H(\mathbf{x})+\eta^{2}G(\mathbf{x})$
with $H$ and $G$ denoting the mean curvature and Gaussian curvature of the
level surface $\Gamma_{\eta}$, respectively.

For the convenience of analysis, we introduce the function class
$\mathcal{W}_{q}$, $q\in\mathbb{N}$ that

$\mathcal{W}_{q}:=\Set{f\in C^{q-1}(\mathbb{R})\
}{\begin{array}[]{l}\operatorname{supp}f=[-\varepsilon,\varepsilon]\;\text{
and }\;\displaystyle\int_{-\varepsilon}^{\varepsilon}f(s)ds=1,\\\ \\\
f^{(q)}\text{ is a piecewise $C^{1}$ function on
$\mathbb{R}$.}\end{array}},\;\text{ if }q\geq 1$

and

$\mathcal{W}_{0}:=\Set{f\in L^{\infty}(\mathbb{R})\
}{\begin{array}[]{l}\operatorname{supp}f=[-\varepsilon,\varepsilon]\;\text{
and }\;\displaystyle\int_{-\varepsilon}^{\varepsilon}f(s)ds=1,\\\ \\\ f\text{
is a piecewise $C^{1}$ function on $\mathbb{R}$.}\end{array}}.$

Here $f^{(q)}$ denotes the $q$-th derivative of $f$. The index $q$ denotes the
order of regularity of the weight function. Usual choices such as

$\theta_{\varepsilon}^{\cos}(s)=\frac{1}{2\varepsilon}\left(1+\cos(\varepsilon^{-1}\pi
s)\right)\in\mathcal{W}_{2},\quad\theta_{\varepsilon}^{\Delta}(s)=\frac{1}{2\varepsilon}(1-\varepsilon^{-1}|s|)\in\mathcal{W}_{1}.$
(2)

Figure 1: Shape of the the weight functions $\theta_{0.5}^{\cos}$ and
$\theta_{0.5}^{\Delta}$ defined in (2).

have been adopted in the literature; see, for instance, references
[zhong2018implicit, kublik2013implicit, chen2017implicit,
engquist2005discretization, izzo2022corrected, izzo2023high]. These functions
are visualized in Figure 1. More choices of weight functions with higher
regularities can be found in [tornberg2004numerical].

The implicit boundary integral method (IBIM) [chen2017implicit,
kublik2013implicit, kublik2016integration, kublik2018extrapolative] uses the
summation over the regular lattice $(h\mathbb{Z})^{d}$ to approximate the
integral (1):

$\displaystyle\mathcal{I}_{h}(f)$
$\displaystyle=h^{d}\sum_{\mathbf{n}\in(h\mathbb{Z})^{d}\cap
T_{\varepsilon}}f(P_{\Gamma}(\mathbf{n}))\theta_{\varepsilon}(d_{\Gamma}(\mathbf{n}))J_{\varepsilon}(\mathbf{n},d_{\Gamma}(\mathbf{n})).$
(3)

We use $\Theta$ and $\mathcal{O}$ to denote the standard Landau notations for
big-Theta and big-O, respectively. When the width parameter
$\varepsilon=\Theta(h^{\alpha})$, $\alpha\in[0,1)$, the theoretical quadrature
error previously derived is at the order of $\mathcal{O}(h^{(q+1)(1-\alpha)})$
for a general domain with smooth boundary [engquist2005discretization].
However, the total number of lattice points is at the order
$\mathcal{O}(\varepsilon h^{-d})$, and the width parameter is often chosen as
$\varepsilon=\Theta(h)$. This turns the classical quadrature error into
$\mathcal{O}(1)$ in such a circumstance, a clearly undesired abrupt loss of
accuracy.

The objective of this paper is to perform a more careful investigation of the
above issue of vanishing accuracy. We perform a systematic error analysis
using the Poisson Summation Formula for boundaries of different orders of
regularity. For strongly convex smooth boundaries, we show that the quadrature
error gains an addition order of $\frac{d-1}{2}$ from the curvature of the
boundaries due to the “randomness” in the signed distance functions for such
boundaries. This additional gain avoids the accuracy catastrophe we see above.

The remainder of this paper is organized as follows. We first discuss in
Section 2 the error estimate of the implicit boundary integral method for
strongly convex boundaries and the error statistics. We then investigate the
issues and possible solutions for a general smooth convex boundary for 2D in
Section 3. Generalization of the estimate of implicit boundary integral to
general open curves in 2D is then presented in Section 4. Concluding remarks
are offered in Section LABEL:SEC:_5.

##  2 Error analysis for smooth convex boundaries in $\mathbb{R}^{d}$

We start with the simplest case when the boundary $\Gamma\in C^{\infty}$ is
strongly convex; that is, the Gaussian curvature is bounded away from zero
uniformly. To setup the notation, we denote by

$\mathcal{Q}(\mathbf{x}):=f(P_{\Gamma}(\mathbf{x}))\theta_{\varepsilon}(d_{\Gamma}(\mathbf{x}))J_{\varepsilon}(\mathbf{x},d_{\Gamma}(\mathbf{x}))$

the integrand in (1). We also denote the mollifier
$\psi_{\delta}(\mathbf{x})=\delta^{-d}\psi(\delta^{-1}\mathbf{x})$, where
$\psi(\mathbf{u})=\psi(|\mathbf{u}|)\in C_{0}^{\infty}([-1,1])$ is the bump
function that

$\int_{-1}^{1}\psi(s)ds=1\,.$

###  2.1 Poisson Summation Formula

We use the standard notation $\mathcal{S}(\mathbb{R}^{d})$ for the Schwartz
function space on $\mathbb{R}^{d}$. For $f\in\mathcal{S}(\mathbb{R}^{d})$, the
Fourier transform of $f$ is defined by

$\widehat{f}(\boldsymbol{\zeta})=\int_{\mathbb{R}^{d}}f(\mathbf{x})e^{-2\pi
i\mathbf{x}\cdot\boldsymbol{\zeta}}d\mathbf{x}.$

The following Poisson Summation Formula is a standard tool whose proof can be
found in [stein1993harmonic].

######  Lemma 2.1 (Poisson Summation Formula [stein1993harmonic]).

If $f$ satisfies that
$|f(\mathbf{u})|+|\widehat{f}(\mathbf{u})|=\mathcal{O}((1+|\mathbf{u}|)^{-d-\mu})$
for all $\mathbf{u}\in\mathbb{R}^{d}$ and certain $\mu>0$, then

$\sum_{\mathbf{n}\in\mathbb{Z}^{d}}f(\mathbf{n})=\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\widehat{f}(\mathbf{n}).$
(4)

In particular, the above formula holds for $f\in\mathcal{S}(\mathbb{R}^{d})$.

Let $h\in\mathbb{R}_{+}$ be the lattice resolution in the summation (3). If we
apply the Poisson Summation Formula to
$\mathcal{Q}_{h,\delta}(\mathbf{x}):=\mathcal{Q}(h{\mathbf{x}})\ast\psi_{\delta}\in\mathcal{S}(\mathbb{R}^{d})$,
we have

$\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\mathcal{Q}_{h,\delta}(\mathbf{n})=\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\widehat{\mathcal{Q}_{h,\delta}}(\mathbf{n})=h^{-d}\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\widehat{\mathcal{Q}}(h^{-1}{\mathbf{n}})\widehat{\psi}(\delta\mathbf{n})\,.$
(5)

For the term $\mathbf{n}={\mathbf{0}}$ in the last summation, it equals to

$h^{-d}\widehat{\mathcal{Q}}({\mathbf{0}})\widehat{\psi}({\mathbf{0}})=h^{-d}\int_{\mathbb{R}^{d}}\mathcal{Q}(\mathbf{x})d\mathbf{x}\int_{\mathbb{R}^{d}}\psi(\mathbf{x})d\mathbf{x}=h^{-d}\int_{\mathbb{R}^{d}}\mathcal{Q}(\mathbf{x})d\mathbf{x}.$
(6)

which produces the desired integral in (1). The objective of the rest of this
section is to provide an estimate for the remaining terms that
$\mathbf{n}\neq{\mathbf{0}}$.

###  2.2 Main results on strongly convex boundaries

We will need the following two lemmas whose proofs are given in Appendices
LABEL:PRF:_STATION and LABEL:PRF:_ESTIMATE, respectively.

######  Lemma 2.2.

Let $\Gamma\in C^{\infty}$ be closed and strongly convex. Then there exists a
constant $C>0$ such that when $\varepsilon\leq C$,
$|\widehat{\mathcal{Q}}(\boldsymbol{\zeta})|=\mathcal{O}(\varepsilon^{-(q+1)}|\boldsymbol{\zeta}|^{-(d+1)/2-q})$.

######  Lemma 2.3.

Let $\Gamma\in C^{\infty}$ be closed and strongly convex. Then the remaining
terms in (5) are bounded by

$\sum_{\mathbf{n}\in\mathbb{Z}^{d},\mathbf{n}\neq{\mathbf{0}}}\widehat{\mathcal{Q}}(h^{-1}{\mathbf{n}})\widehat{\psi}(\delta\mathbf{n})=\begin{cases}\mathcal{O}(\varepsilon^{-(q+1)}h^{\frac{d+1}{2}+q}\delta^{-\frac{d-1}{2}+q})\quad&\frac{d-1}{2}>q,\\\
\mathcal{O}(\varepsilon^{-(q+1)}h^{\frac{d+1}{2}+q}|\log\delta|)\quad&\frac{d-1}{2}=q,\\\
\mathcal{O}(\varepsilon^{-(q+1)}h^{\frac{d+1}{2}+q})\quad&\frac{d-1}{2}<q.\end{cases}$
(7)

Lemma 2.3 divides the regularity order $q$ into three classes: sub-critical
regularity ($q<\frac{d-1}{2}$), critical regularity ($q=\frac{d-1}{2}$), and
super-critical regularity ($q>\frac{d-1}{2}$). Base d on this result, we
prove, in the following subsections, the estimates of the error bound for
implicit boundary integral based on different regularity regimes of the weight
function $\theta_{\varepsilon}$. We summarize the resulting estimates in Tab
2.2.

Regularity regime | Error bound of $|\mathcal{I}f-\mathcal{I}_{h}f|$
---|---
$\min(2,\frac{d-1}{2})>q\geq 1$ | $\mathcal{O}(h^{\frac{2d-\alpha(d+1)}{d+1-2q}})$
$\frac{d-1}{2}>q\geq 2$ | $\mathcal{O}(h^{\frac{4d-2\alpha(d+1)}{d+3-2q}})$
$\frac{d-1}{2}=q$ | $\mathcal{O}(h^{\frac{d-1}{2}+(q+1)(1-\alpha)}|\log h|)$
$\frac{d-1}{2}<q$ | $\mathcal{O}(h^{\frac{d-1}{2}+(q+1)(1-\alpha)})$
|

####  2.2.1 Sub-critical regularity

We consider two separate cases.

######  Theorem 2.4.

Under the assumption that $\Gamma\in C^{\infty}$ is closed and strongly
convex, if $\varepsilon=\Theta(h^{\alpha})$ and $\min(2,\frac{d-1}{2})>q\geq
1$, then
$|\mathcal{I}(f)-\mathcal{I}_{h}(f)|=\mathcal{O}(h^{\frac{2d-\alpha(d+1)}{d+1-2q}})$.

###### Proof.

First, by Lemma 2.3, we have,
$\left|\mathcal{I}(f)-h^{d}\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\mathcal{Q}_{h,\delta}(\mathbf{n})\right|=\mathcal{O}(\varepsilon^{-q-1}h^{\frac{d+1}{2}+q}\delta^{-\frac{d-1}{2}+q})\,.$
(8) We also estimate the bound for
$\displaystyle\left|\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\mathcal{Q}(h\mathbf{n})-\mathcal{Q}_{h,\delta}(\mathbf{n})\right|$
$\displaystyle=\left|\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\int_{\mathbb{R}^{d}}\left(\mathcal{Q}(h\mathbf{n})-\mathcal{Q}(h(\mathbf{n}-\mathbf{y}))\right)\psi_{\delta}(\mathbf{y})d\mathbf{y}\right|\,.$
(9) The summation has at most $\mathcal{O}(\varepsilon h^{-d})$ terms, and we
also have the trivial estimate
$\mathcal{Q}(h(\mathbf{n}-\mathbf{y}))-\mathcal{Q}(h\mathbf{n})=\mathcal{O}(Lh\delta)\,,$
(10) where $L=\mathcal{O}(\varepsilon^{-2})$ is the Lipschitz constant of
$\mathcal{Q}$. This leads to
$\left|\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\mathcal{Q}(h\mathbf{n})-\mathcal{Q}_{h,\delta}(\mathbf{n})\right|=\mathcal{O}(\varepsilon^{-1}h^{-d+1}\delta)\,.$
(11) Combining the estimates (8) and (11) gives us
$\left|\mathcal{I}(f)-\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\mathcal{Q}(h\mathbf{n})\right|=\mathcal{O}\left(\varepsilon^{-q-1}h^{\frac{d+1}{2}+q}\delta^{-\frac{d-1}{2}+q}+\varepsilon^{-1}h\delta\right)\,.$
(12) Let $\varepsilon=\Theta(h^{\alpha})$ and balance the orders between $h$
and $\delta$ in (12), we conclude that the balance is attained at
$\delta=\mathcal{O}(h^{\frac{d-1-2q(\alpha-1)}{d+1-2q}})$, and
$\left|\mathcal{I}(g)-h^{d}\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\mathcal{Q}(h\mathbf{n})\right|=\mathcal{O}(h^{\frac{2d-\alpha(d+1)}{d+1-2q}})\,.$
(13) This completes the proof. ∎ In the second case, we have the following
result.

######  Theorem 2.5.

Under the assumption that $\Gamma\in C^{\infty}$ is closed and strongly
convex, if $\varepsilon=\Theta(h^{\alpha})$ and $\frac{d-1}{2}>q\geq 2$, then
$|\mathcal{I}(f)-\mathcal{I}_{h}(f)|=\mathcal{O}(h^{\frac{4d-2\alpha(d+1)}{d+3-2q}})$.

###### Proof.

Similar to the previous case, we use Lemma 2.3 to conclude that
$\left|\mathcal{I}(f)-h^{d}\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\mathcal{Q}_{h,\delta}(\mathbf{n})\right|=\mathcal{O}(\varepsilon^{-q-1}h^{\frac{d+1}{2}+q}\delta^{-\frac{d-1}{2}+q})\,.$
(14) Since $q\geq 2$, $\mathcal{Q}\in C^{1,1}(\mathbb{R}^{d})$, we can take
Taylor expansion
$\mathcal{Q}(h(\mathbf{n}-\mathbf{y}))-\mathcal{Q}(h\mathbf{n})=-h\mathbf{y}\cdot\nabla\mathcal{Q}(h\mathbf{n})+\mathcal{O}(Kh^{2}\delta^{2})$
(15) where $K=\mathcal{O}(\varepsilon^{-3})$ is the Lipschitz constant of
$\nabla\mathcal{Q}$. Using the spherical symmetry of $\psi_{\delta}$, the
first term on the right-hand side of (15) will be canceled. We thus have
$\left|\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\mathcal{Q}(h\mathbf{n})-\mathcal{Q}_{h,\delta}(\mathbf{n})\right|=\mathcal{O}(\varepsilon^{-2}h^{-d+2}\delta^{2})\,,$
(16) which leads to
$\left|\mathcal{I}(f)-\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\mathcal{Q}(h\mathbf{n})\right|=\mathcal{O}\left(\varepsilon^{-q-1}h^{\frac{d+1}{2}+q}\delta^{-\frac{d-1}{2}+q}+\varepsilon^{-2}h^{2}\delta^{2}\right)\,.$
(17) The bounds attain balance when
$\delta=\mathcal{O}(h^{\frac{d-1+(\alpha-1)(2-2q)}{d+3-2q}})$ at which point
we have —I(f) - ∑_n∈Z^d Q(hn)— = O(h^4d- 2α(d+1)d+3-2q) , which is the desired
result. ∎

####  2.2.2 Critical regularity

The same calculations can be performed in the case of critical regularity. We
have the following result.

######  Corollary 2.6.

Under the assumption that $\Gamma\in C^{\infty}$ is closed and strongly
convex, if $\varepsilon=\Theta(h^{\alpha})$ and $\frac{d-1}{2}=q$, then
$|\mathcal{I}(f)-\mathcal{I}_{h}(f)|=\mathcal{O}(h^{\frac{d-1}{2}+(q+1)(1-\alpha)}|\log
h|)$.

###### Proof.

Using the Lemma 2.3 for the critical value $q=\frac{d-1}{2}$, we get
$\left|\mathcal{I}(f)-h^{d}\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\mathcal{Q}_{h,\delta}(\mathbf{n})\right|=\mathcal{O}(\varepsilon^{-q-1}h^{\frac{d+1}{2}+q}|\log\delta|)\,.$
If $q<2$, we reuse (11), and otherwise we reuse (16), to get
$\left|\mathcal{I}(f)-h^{d}\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\mathcal{Q}(h\mathbf{n})\right|=\begin{cases}\mathcal{O}\left(\varepsilon^{-q-1}h^{\frac{d+1}{2}+q}|\log\delta|+\varepsilon^{-1}h\delta\right)\quad&q<2\,,\\\
\mathcal{O}\left(\varepsilon^{-q-1}h^{\frac{d+1}{2}+q}|\log\delta|+\varepsilon^{-2}h^{2}\delta^{2}\right)\quad&q\geq
2\,.\end{cases}$ (18) Let $\varepsilon=\Theta(h^{\alpha})$ and balance the
orders between $h$ and $\delta$ in (18), we find:

  1. (1) If $0\leq q<2$, the balance is achieved when $h^{\frac{d-1}{2}+q-q\alpha}|\log\delta|\approx\delta$, which means the error is $\mathcal{O}(h^{1-\alpha+\frac{d-1}{2}+q-q\alpha}|\log\delta|)$. Since $\frac{d-1}{2}+q-q\alpha>0$, $\delta=\mathcal{O}(h^{\frac{d-1}{2}+q-q\alpha-\tau})$ for any $\tau>0$. Hence $|\log\delta|=\mathcal{O}(|\log h|)$.
  2. (2) If $q\geq 2$, the balance is achieved when $h^{(-q+1)\alpha}h^{\frac{d-3}{2}+q}|\log\delta|\approx\delta^{2}$ and the quadrature error is $\mathcal{O}(h^{-(q+1)\alpha}h^{\frac{d+1}{2}+q}|\log\delta|)$. Because $\frac{d-3}{2}+q+(1-q)\alpha>0$, we find that $\delta^{2}=\mathcal{O}(h^{\frac{d-3}{2}+q+(1-q)\alpha-\tau})$ for certain $\tau>0$ and it implies $|\log\delta|=\mathcal{O}(|\log h|)$.

The error bounds in both cases are the same; that is,
$\left|\mathcal{I}(g)-h^{d}\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\mathcal{Q}(h\mathbf{n})\right|=\mathcal{O}(h^{\frac{d-1}{2}+(q+1)(1-\alpha)}|\log
h|)\,.$ (19) The finishes the proof. ∎

####  2.2.3 Super-critical regularity

The bound for the super-critical regularity case follows in the same way.

######  Corollary 2.7.

Under the assumption that $\Gamma\in C^{\infty}$ is closed and strongly
convex, if $\varepsilon=\Theta(h^{\alpha})$ and $\frac{d-1}{2}<q$, then
$|\mathcal{I}(f)-\mathcal{I}_{h}(f)|=\mathcal{O}(h^{\frac{d-1}{2}+(q+1)(1-\alpha)})$.

###### Proof.

Using Lemma 2.3, the estimate (18) now takes the form
$\left|\mathcal{I}(f)-h^{d}\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\mathcal{Q}(h\mathbf{n})\right|=\begin{cases}\mathcal{O}\left(\varepsilon^{-q-1}h^{\frac{d+1}{2}+q}+\varepsilon^{-1}h\delta\right)\quad&q<2\,,\\\
\mathcal{O}\left(\varepsilon^{-q-1}h^{\frac{d+1}{2}+q}+\varepsilon^{-2}h^{2}\delta^{2}\right)\quad&q\geq
2\,.\end{cases}$ (20) Therefore, the error estimate becomes
$\left|\mathcal{I}(g)-h^{d}\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\mathcal{Q}(h\mathbf{n})\right|=\mathcal{O}(h^{\frac{d-1}{2}+(q+1)(1-\alpha)})\,,$
which is what we need to prove. ∎ The tube width usually takes $\alpha=1$ for
practical applications of the implicit-boundary-integral method. If
$\theta_{\varepsilon}$ has sub-critical regularity $q<\frac{d-1}{2}$, we have
the error exponents as $\frac{d-1}{d+1-2q}$ and $\frac{2(d-1)}{d+3-2q}$ in
Table 2.2 for different regularity classes. These errors may still be far from
optimal, and the sharper exponents need more sophisticated estimates for the
remaining terms of the Poisson Summation Formula (4). The problem is closely
related to the generalized cases of the Gauss Circle Problem (GCP)
[berndt2018circle, heath1999lattice] that counts the number of lattice points
inside the domain. This can be reviewed as a special case of implicit boundary
integral by taking
$\theta_{\varepsilon}=\theta_{\varepsilon}^{C}=\frac{1}{2\varepsilon}\chi_{T_{\varepsilon}}$
as the characteristic function on the tube and $f\equiv 1$ with an
approximation of $J\approx 1$. In such case, a discontinuity of $\mathcal{Q}$
exists across the boundary $\partial T_{\varepsilon}$. It will involve
additional boundary contributions in (12) and (17). We will provide an
interesting example in Remark 2.8 below. The sharp discrepancy estimate of GCP
is still an open question. Recent developments in GCP and its variants can be
found in [heath1999lattice, berndt2018circle, bourgain2017mean,
huxley2003exponential, huxley2005exponential] and references therein.

######  Remark 2.8.

One interesting case is when the weight function has a jump across the
boundary. For instance, the case where
$\theta^{C}_{\varepsilon}=\frac{1}{2\varepsilon}\chi_{T_{\varepsilon}}\in\mathcal{W}_{0}$
is the characteristic function on the tube. Due to the jump, we cannot
directly reuse (10) for (9). We can estimate the summation in three categories
$\left|\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\mathcal{Q}(h\mathbf{n})-\mathcal{Q}_{h,\delta}(\mathbf{n})\right|\leq\mathcal{V}_{in}+\mathcal{V}_{bd}+\mathcal{V}_{out},$
where the summation $\mathcal{V}_{in}$, $\mathcal{V}_{bd}$, and
$\mathcal{V}_{out}$ are respectively $\displaystyle\mathcal{V}_{in}$
$\displaystyle=\left|\sum_{\mathcal{S}_{in}}\mathcal{Q}(h\mathbf{n})-\mathcal{Q}_{h,\delta}(\mathbf{n})\right|,\quad\text{
where }\mathcal{S}_{in}:=\\{\mathbf{n}\mid h\mathbf{n}\in
T_{\varepsilon},\quad B_{h\mathbf{n}}(h\delta)\subset T_{\varepsilon}\\}$
$\displaystyle\mathcal{V}_{bd}$
$\displaystyle=\left|\sum_{\mathcal{S}_{bd}}\mathcal{Q}(h\mathbf{n})-\mathcal{Q}_{h,\delta}(\mathbf{n})\right|,\quad\text{
where }\mathcal{S}_{bd}:=\\{\mathbf{n}\mid h\mathbf{n}\in
T_{\varepsilon},\quad B_{h\mathbf{n}}(h\delta)\not\subset T_{\varepsilon}\\}$
$\displaystyle\mathcal{V}_{out}$
$\displaystyle=\left|\sum_{\mathcal{S}_{out}}\mathcal{Q}(h\mathbf{n})-\mathcal{Q}_{h,\delta}(\mathbf{n})\right|,\quad\text{
where }\mathcal{S}_{out}:=\\{\mathbf{n}\mid h\mathbf{n}\not\in
T_{\varepsilon},\quad B_{h\mathbf{n}}(h\delta)\cap
T_{\varepsilon}\neq\emptyset\\}$ We can estimate $\mathcal{V}_{in}$ by
$\mathcal{O}(\varepsilon h^{-d+2}\delta^{2})$ by the second order Taylor
expansion at $h\mathbf{n}$. For both $\mathcal{V}_{bd}$ and
$\mathcal{V}_{out}$, there are $\mathcal{O}(\delta h^{-d+1})$ lattice points,
and each summand is at order $\mathcal{O}(\varepsilon^{-1})$. Therefore, we
obtain the estimate
$\left|\mathcal{I}(f)-h^{d}\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\mathcal{Q}(h\mathbf{n})\right|=\mathcal{O}(\varepsilon^{-1}h^{\frac{d+1}{2}}\delta^{-\frac{d-1}{2}}+h^{2}\delta^{2}+\varepsilon^{-1}\delta
h)\,.$ (21) The balance is attained at
$\delta\sim\mathcal{O}(h^{\frac{d-1}{d+1}})$ and error is
$\mathcal{O}(h^{\frac{2d-\alpha(d+1)}{d+1}})$ for
$\varepsilon=\Theta(h^{\alpha})$. In particular, when $\alpha=1$, we find that
the error is bounded by $\mathcal{O}(h^{\frac{d-1}{d+1}})$. As a comparison,
in this case, a standard Monte Carlo integration for $\alpha=1$ gives a
smaller error, bounded by $\mathcal{O}(h^{\frac{d-1}{2}})$. For fixed tube
width of $\mathcal{O}(1)$ ($\alpha=0$), the error estimates are
$\mathcal{O}(h^{\frac{2d}{d-1}})$ for $\theta_{\varepsilon}^{\Delta}$ ($q=1$
and $d>3$), which is almost second order convergence in high dimensions. It
exceeds the theoretical estimate $\mathcal{O}(h)$. For
$\theta_{\varepsilon}^{\cos}$ ($q=2$ and $d>5$), the error becomes
$\mathcal{O}(h^{\frac{4d}{d-1}})$, which is almost 4th order when the
dimension is high, exceeding the theoretical estimate $\mathcal{O}(h^{2})$. We
gain these additional convergence rates from the strong convexity in high
dimensions. However, for low dimensional surfaces of $d=2$ or $d=3$, once the
regularity $q\geq 1\geq\frac{d-1}{2}$, we have a generic error estimate
$\mathcal{O}(h^{\frac{d+1}{2}+q}|\log h|)$ or
$\mathcal{O}(h^{\frac{d+1}{2}+q})$.

###  2.3 Variance of error under random translations

We now characterize the variance of the quadrature error caused by random
shifts of the lattice, an important issue for practical applications
[izzo2022corrected, izzo2023high, engquist2005discretization]. We denote by
$\boldsymbol{\zeta}\in\mathbb{R}^{d}$ the random shift, and denote by
$\mathcal{I}_{h}(f;\boldsymbol{\zeta})$ the implicit-boundary-integral
quadrature with the random shift; that is, I_h(f;ξ) := h^d ∑_n∈Z^d Q(h(n\+ ξ))
. It is clear that the quadrature error is $h$-periodic, hence we only have to
consider a random offset $\boldsymbol{\xi}\in[0,1]^{d}$. First, we note that
$\mathbb{E}_{\boldsymbol{\xi}}\mathcal{I}_{h}(f;\boldsymbol{\xi})=\mathcal{I}(f)$,
since
$\mathbb{E}_{\boldsymbol{\xi}}\mathcal{I}_{h}(f;\boldsymbol{\xi})=h^{d}\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\int_{[0,1]^{d}}\mathcal{Q}(h(\mathbf{n}+\boldsymbol{\xi}))d\boldsymbol{\xi}=\int_{\mathbb{R}^{d}}\mathcal{Q}(\mathbf{x})d\mathbf{x}=\mathcal{I}(f).$
We therefore define the variance of $\mathcal{I}_{h}(f;\boldsymbol{\zeta})$ as
$\textrm{Var}_{\boldsymbol{\xi}}[\mathcal{I}_{h}(f;\boldsymbol{\xi})]:=\int_{\boldsymbol{\xi}\in[0,1]^{d}}\left|\mathcal{I}(f)-h^{d}\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\mathcal{Q}(h(\mathbf{n}+\boldsymbol{\xi}))\right|^{2}d\boldsymbol{\xi}\,.$
(22) We are interested in providing an estimate for this variance. The
following theorem follows the idea of the classical result by Kendall
[kendall1948number]. General $L^{p}$ estimates can be derived similarly; see,
for instance, [brandolini2015lp, huxley2014fourth].

######  Theorem 2.9.

Let $\Gamma\in C^{\infty}$ be closed and strongly convex. If
$\varepsilon=\Theta(h^{\alpha})$, then
$\mathrm{Var}_{\boldsymbol{\xi}}\left[\mathcal{I}_{h}(f;\boldsymbol{\xi})\right]=\mathcal{O}(h^{d-1+2(q+1)(1-\alpha)})$.

###### Proof.

Let
$\mathcal{A}(h,\boldsymbol{\zeta}):=h^{d}\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\mathcal{Q}(h(\mathbf{n}+\boldsymbol{\zeta}))$.
It is then clear that $\mathcal{A}(h,\boldsymbol{\zeta})$ has a Fourier series
expansion that converges in $L^{2}([0,1]^{d})$:
$\mathcal{A}(h,\boldsymbol{\zeta})\sim\sum_{\boldsymbol{\xi}\in\mathbb{Z}^{d}}a_{\boldsymbol{\xi}}(h)e^{2\pi
i(\boldsymbol{\zeta}\cdot\boldsymbol{\xi})}$ where the Fourier coefficients
are given as $\displaystyle a_{\boldsymbol{\xi}}(h)$
$\displaystyle=\int_{[0,1]^{d}}\mathcal{A}(h,\boldsymbol{\zeta})e^{-2\pi
i(\boldsymbol{\zeta}\cdot\boldsymbol{\xi})}d\boldsymbol{\zeta}$
$\displaystyle=h^{d}\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\int_{[0,1]^{d}}\mathcal{Q}(h(\mathbf{n}+\boldsymbol{\zeta}))e^{-2\pi
i(\boldsymbol{\xi}\cdot\boldsymbol{\zeta})}d\boldsymbol{\zeta}$
$\displaystyle=h^{d}\sum_{\mathbf{n}\in\mathbb{Z}^{d}}\int_{[0,1]^{d}}\mathcal{Q}(h(\mathbf{n}+\boldsymbol{\zeta}))e^{-2\pi
i(\boldsymbol{\xi}\cdot\boldsymbol{\zeta})}d\boldsymbol{\zeta}$
$\displaystyle=\int_{\mathbb{R}^{d}}\mathcal{Q}(\mathbf{x})e^{-2\pi
i(h^{-1}\boldsymbol{\xi}\cdot\mathbf{x})}d\mathbf{x}$
$\displaystyle=\widehat{\mathcal{Q}}(h^{-1}\boldsymbol{\xi})\,.$ Notice that
$a_{{\mathbf{0}}}=\mathcal{I}(f)$. Therefore, we have
$\mathrm{Var}_{\boldsymbol{\xi}}\left[\mathcal{I}_{h}(f;\boldsymbol{\xi})\right]=\sum_{\boldsymbol{\xi}\in\mathbb{Z}^{d},|\boldsymbol{\xi}|\neq
0}|a_{\boldsymbol{\xi}}(h)|^{2}=\sum_{\boldsymbol{\xi}\in\mathbb{Z}^{d},|\boldsymbol{\xi}|\neq
0}|\widehat{\mathcal{Q}}(h^{-1}\boldsymbol{\xi})|^{2}\,.$ Using Lemma 2.2,
with $q\geq 0$, we get
$\displaystyle\mathrm{Var}_{\boldsymbol{\xi}}\left[\mathcal{I}_{h}(f;\boldsymbol{\xi})\right]$
$\displaystyle=\mathcal{O}(\varepsilon^{-2(q+1)}h^{d+1+2q})\sum_{\boldsymbol{\xi}\in\mathbb{Z}^{d},|\boldsymbol{\xi}|\neq
0}|\boldsymbol{\xi}|^{-(d+1)-2q}$
$\displaystyle=\mathcal{O}(h^{d-1+2(q+1)(1-\alpha)})\,.$ The proof is
complete. ∎ A probabilistic consequence of the above variance characterization
is the following. If one randomly chooses a shift
$\boldsymbol{\xi}\in\mathbb{R}^{d}$ in the implicit boundary integration, then
by the Chebyshev inequality, there exists a constant $C$, depending on both
$f$ and the surface $\Gamma$, that
$\mathbb{P}\Big{(}|\mathcal{I}_{h}(f;\boldsymbol{\xi})-\mathcal{I}(f)|\geq
kh^{\frac{d-1}{2}+(q+1)(1-\alpha)}\Big{)}<Ck^{-2},\ \ k>0\,.$ Hence, for
weight functions with sub-critical regularity $q<\frac{d-1}{2}$, this implies
that it is rare to achieve the possible worst bounds shown in Theorem 2.4 and
Theorem 2.5.

###  2.4 Numerical experiments

In this section, we numerically verify (i) the quadrature error bounds shown
in Tab 2.2 on circles in 2D and spheres in 3D; and (ii) the variance of
quadrature error under random translations. The reference values of the
boundary integrals are computed by Mathematica to machine precision. The
MATLAB source code of all numerical experiments is hosted in the GitHub
repository111https://github.com/lowrank/ibim-error-analysis-experiments.

####  2.4.1 2D circle

The 2D experiments are performed on the circle $\Gamma$:
$\frac{(x-x_{0})^{2}}{r^{2}}+\frac{(y-y_{0})^{2}}{r^{2}}=1$ with
$r=\frac{3}{4}$ and $(x_{0},y_{0})$ a randomly sampled point. The test
integrand function $f:\mathbb{R}^{2}\mapsto\mathbb{R}$ is given as
$f(x,y)=\cos(x^{2}-y)\sin(y^{2}-x^{3})\,.$ With the weight function
$\theta_{\varepsilon}^{\Delta}\in\mathcal{W}_{1}$, the quadrature error has a
upper bound estimate $\mathcal{O}(h^{\frac{1}{2}+2(1-\alpha)})$. As we can see
in Figure 2, for $\alpha=1$, $\alpha=\frac{1}{2}$ and $\alpha=0$, the
estimated error bounds $\mathcal{O}(h^{\frac{1}{2}})$,
$\mathcal{O}(h^{\frac{3}{2}})$ and $\mathcal{O}(h^{\frac{5}{2}})$ appear
consistent with the numerical results. Figure 2: Quadrature error with weight
function $\theta_{\varepsilon}^{\Delta}$ with respect to grid size $h$. Left:
tube width $\varepsilon=2h$. Middle: tube width
$\varepsilon=2h^{\frac{1}{2}}$. Right: tube width $\varepsilon=0.1$. With the
weight function $\theta_{\varepsilon}^{\cos}\in\mathcal{W}_{2}$, the
quadrature error has an upper bound estimate
$\mathcal{O}(h^{\frac{1}{2}+3(1-\alpha)})$. In Figure 3, we find that for
$\alpha=1$, $\alpha=\frac{1}{2}$ and $\alpha=0$, the estimated error bounds
$\mathcal{O}(h^{\frac{1}{2}})$, $\mathcal{O}(h^{2})$ and
$\mathcal{O}(h^{\frac{7}{2}})$ match the numerical results approximately.
Figure 3: Quadrature error with weight function $\theta_{\varepsilon}^{\cos}$
with respect to grid size $h$. Left: tube width $\varepsilon=2h$. Middle: tube
width $\varepsilon=2h^{\frac{1}{2}}$. Right: tube width $\varepsilon=0.1$. The
variances of error are plotted in Figure 4 and Figure 5 with 32 random shifts,
consistent with the result of Theorem 2.9. Figure 4: Variance of quadrature
error with weight function $\theta_{\varepsilon}^{\Delta}$ with respect to
grid size $h$. Left: tube width $\varepsilon=2h$. Middle: tube width
$\varepsilon=2h^{\frac{1}{2}}$. Right: tube width $\varepsilon=0.1$.  Figure
5: Variance of quadrature error with weight function
$\theta_{\varepsilon}^{\cos}$ with respect to grid size $h$. Left: tube width
$\varepsilon=2h$. Middle: tube width $\varepsilon=2h^{\frac{1}{2}}$. Right:
tube width $\varepsilon=0.1$.

####  2.4.2 3D sphere

The 3D experiments are performed on the sphere $\Gamma$:
$\frac{(x-x_{0})^{2}}{r^{2}}+\frac{(y-y_{0})^{2}}{r^{2}}+\frac{(z-z_{0})^{2}}{r^{2}}=1$
with $r=\frac{3}{4}$ and $(x_{0},y_{0},z_{0})$ a randomly sampled point. The
test integrand function $f:\mathbb{R}^{2}\mapsto\mathbb{R}$ in this case is
$f(x,y,z)=\cos(x^{2}-y-z^{3})\sin(y^{2}-x^{3}-z).$ With the weight function
$\theta_{\varepsilon}^{\Delta}\in\mathcal{W}_{1}$ (resp.
$\theta_{\varepsilon}^{\cos}\in\mathcal{W}_{2}$), the quadrature errors are
shown in Figure 6 (resp. Figure 7) for $\alpha=1$, $\alpha=\frac{1}{2}$ and
$\alpha=0$. The estimated error bounds are still consistent with the rate
$\mathcal{O}(h^{1+2(1-\alpha)})$ (resp. $\mathcal{O}(h^{1+3(1-\alpha)})$). The
corresponding variance of quadrature error is plotted in Figure 8 (resp.
Figure 9). Figure 6: Quadrature error with weight function
$\theta_{\varepsilon}^{\Delta}$ with respect to grid size $h$. Left: tube
width $\varepsilon=2h$. Middle: tube width $\varepsilon=2h^{\frac{1}{2}}$.
Right: tube width $\varepsilon=0.1$.  Figure 7: Quadrature error with weight
function $\theta_{\varepsilon}^{\cos}$ with respect to grid size $h$. Left:
tube width $\varepsilon=2h$. Middle: tube width
$\varepsilon=2h^{\frac{1}{2}}$. Right: tube width $\varepsilon=0.1$.  Figure
8: Variance of quadrature error with weight function
$\theta_{\varepsilon}^{\Delta}$ with respect to grid size $h$. Left: tube
width $\varepsilon=2h$. Middle: tube width $\varepsilon=2h^{\frac{1}{2}}$.
Right: tube width $\varepsilon=0.1$.  Figure 9: Variance of quadrature error
with weight function $\theta_{\varepsilon}^{\cos}$ with respect to grid size
$h$. Left: tube width $\varepsilon=2h$. Middle: tube width
$\varepsilon=2h^{\frac{1}{2}}$. Right: tube width $\varepsilon=0.1$.

######  Remark 2.10.

For the efficiency of computation, when the tube width
$\varepsilon=\Theta(h)$, one can apply the approximated Jacobian factor
$J_{\varepsilon}\approx 1$ instead of the accurate one if the dimension
$d=2,3$, without altering the error estimate $\mathcal{O}(h^{\frac{d-1}{2}})$.
If the dimension $d=4,5$, one can alternatively use the approximation
$J_{\varepsilon}(\mathbf{x},\eta)\approx 1-d_{\Gamma}(\mathbf{x})\Delta
d_{\Gamma}(\mathbf{x})$, which can be effectively computed by a local central
difference scheme.

######  Remark 2.11.

It is worth noting that when the boundary is strongly convex, the quadrature
error $\mathcal{O}(h^{\frac{d-1}{2}})$ for a thin tube with width
$\varepsilon=\Theta(h)$ is somewhat equivalent to the error using a standard
Monte Carlo method. This seems to imply certain equidistributed “randomness”
of the lattice points inside the tube. However, this is still an open question
even for spheres unless the tube width $\varepsilon$ is at least
$\Theta(h^{\frac{11}{16}})$, which can be derived by combining the Theorem 1
of [duke1990representation] and the lattice count theorem in
[heath1999lattice].

###  2.5 Convex but not strongly convex boundaries

It turns out that the results in this section can be modified slightly for
closed, convex but not strongly convex, surface $\Gamma\in C^{\infty}$ with at
least one positive principal curvature. In this case, the stationary phase
estimate used in proving Lemma 2.2 gives a degenerate result (see Remark
LABEL:REMARK:Convex), which leads to a modified version of Lemma 2.2 in which
one replaces $d$ with $\Lambda+1$ ($\Lambda$ being the number of positive
principle curvature of $\Gamma$). Based on this revised version of Lemma 2.2,
we can reproduce all the theorems in this section to have the following
result.

######  Corollary 2.12.

Let $\Gamma\in C^{\infty}$ be closed and convex with at least $\Lambda>0$
strictly positive principal curvatures everywhere. Then every theorem of
Section 2.2 and Section 2.3 holds with $d$ replaced by $(\Lambda+1)$. As an
extreme example, in the following Figure 10, we show that if the boundary
$\Gamma$ has a segment component parallel to $y=x$, where $\Lambda=0$,
numerical quadrature with tube width $\varepsilon=\Theta(h)$ will introduce an
$\mathcal{O}(1)$ error even when integrating a constant function. In addition,
if one computes the quadrature with random shifts, the standard deviation is
still $\mathcal{O}(1)$. Such an argument can be easily adapted to segments
with rational slopes. Figure 10: $\Gamma$ (red solid line) has a segment
component as $y=x$. The tube, denoted by the red dash line, has a width
$\varepsilon=\frac{3}{2}h$. The green dots are lattice points, and the blue
dash lines are projections to $\Gamma$.

######  Remark 2.13.

Such a difference between smooth convex sets and polytopes has already been
spotted in the famous lattice point problem. The early works trace back to
Hardy and Littlewood [hardy1922some, hardy1922some2] over a century ago. The
lattice point problem studies the discrepancy —t B ∩Z^d— - meas(B) t^d. The
discrepancy could be regarded as random fluctuations if
$B\subset\mathbb{R}^{d}$ has a smooth convex boundary. In contrast, if $B$ is
a polytope, the leading term in the discrepancy may still behave polynomially
[macdonald1971polynomials] in the variable $t$, which conceals the randomness
behind the successive order. In addition, the fluctuations of randomness are
closely related to the Diophantine approximation problem; see
[beck1994probabilistic, matousek1999geometric] and references therein.

##  3 Theories for smooth closed curves in $\mathbb{R}^{2}$

In this section, we will establish a general theory for the implicit boundary
integral on smooth closed curves in $\mathbb{R}^{2}$. The case when the curve
is closed and strongly convex is only a special case of the results in Section
2. We therefore focus only on the cases where the curvature vanishes at
certain points of the curves.

###  3.1 Vanishing curvature at isolated points

We start with a generic case where the curvature vanishes at a finite number
of points. It suffices to consider only one point with vanishing curvature
using a partition of unity argument. We set this point $\mathbf{z}$ as the
origin. We assume that the curvature at $\mathbf{z}$ vanishes to order
$\kappa-2$ for certain $\kappa>2$, which means _locally_ we can arrange the
frame and represent $\Gamma$ as $(x,g(x))\subset\mathbb{R}^{2}$ such that
$g(x)=|x|^{\kappa}h(x),\quad h\in C^{\infty}[-r,r],\quad h(x)\neq
0\quad\forall x\in[-r,r]$ (23) for certain $r>0$. Then in the proof of Lemma
2.2 (see Appendix LABEL:PRF:_STATION), instead of gaining an additional decay
factor $\mathcal{O}(|\boldsymbol{\zeta}|^{1/2})$ for (LABEL:EQ:_STA_PHASE),
the factor becomes $\mathcal{O}(|\boldsymbol{\zeta}|^{1/\kappa})$ in the
spirit of Lemma LABEL:LEM:_van_der_Corput_Revised (a generalized version of
the van der Corput Lemma [stein1993harmonic]). Repeating the proofs in Section
2.2, one can obtain the following result.

######  Corollary 3.1.

Let $d=2$, and $\Gamma$ be a smooth closed curve with vanishing curvatures to
order $\kappa-2$ ($\kappa>2$) at finitely many points. Then, when the tube
width is $\varepsilon=\mathcal{O}(h^{\alpha})$, the quadrature error scales as
—If - I_h f—=O(h^1κ+(q+1)(1-α)) . When $\alpha=1$, the quadrature error is
bounded by $\mathcal{O}(h^{\frac{1}{\kappa}})$ which degenerates to
$\mathcal{O}(1)$ as $\kappa\to\infty$ (that is when the curve locally becomes
straight, as in Figure 10). While this worst quadrature error may be
unsatisfactory for practical uses, under random rigid transformations, the
“average” quadrature error could be greatly improved over this bound.

###  3.2 Variance of error under random rigid transformations

The discussions in Section 2.3 suggest that the distribution of signed
distances on lattice points to a strongly convex smooth boundary behaves
“almost” random so that most rigid transformations (rotations and
translations) of the lattice points will not change the quadrature error too
much. In this section, we show that under random rigid transformations, even
with a finite number of points on $\Gamma$ with vanishing curvatures will
_not_ alter the “average” quadrature error
$\mathcal{O}(h^{\frac{1}{2}+(q+1)(1-\alpha)})$.

######  Theorem 3.2.

Let $\kappa\in(2,\infty)$ and $\Gamma$ be a convex closed $C^{\infty}$ curve
with finitely many points of vanishing curvature of maximal order
$(\kappa-2)$. For any rotation $\eta\in\mathrm{SO}(2)$ and translation
$\boldsymbol{\xi}\in[0,1]^{2}$, let $\mathcal{I}_{h}(f;\eta,\boldsymbol{\xi})$
denote the implicit boundary integral with tube width
$\varepsilon=\Theta(h^{\alpha})$, $\alpha\in[0,1]$ on the transformed boundary
$\eta\Gamma+h\boldsymbol{\xi}$. Then, we have that
$\int_{[0,1]^{2}}\int_{\eta\in\mathrm{SO}(2)}|\mathcal{I}_{h}(f;\eta,\boldsymbol{\xi})-\mathcal{I}(f)|^{2}d\eta
d\boldsymbol{\xi}=\mathcal{O}(h^{1+2(q+1)(1-\alpha)})\,,$ (24) where $d\eta$
is the normalized Haar measure on $\mathrm{SO}(2)$ and $q\geq 0$ is the
regularity order of the weight function.

###### Proof.

Let
$\mathcal{Q}_{\Gamma}(\mathbf{x}):=f(P_{\Gamma}\mathbf{x})\theta_{\varepsilon}(d_{\Gamma}(\mathbf{x}))$.
The quadrature $\mathcal{I}_{h}(f;\eta,\boldsymbol{\xi})$ can be written as
$\mathcal{I}_{h}(f;\eta)=h^{2}\sum_{\mathbf{w}\in\mathbb{Z}^{2}}\mathcal{Q}_{\eta\Gamma+h\boldsymbol{\xi}}(h\mathbf{w})=h^{2}\sum_{\mathbf{w}\in\mathbb{Z}^{2}}\mathcal{Q}_{\Gamma}(h\eta^{-1}(\mathbf{w}-\boldsymbol{\xi})).$
Similar to that in the proof of Theorem 2.9, we have
$\displaystyle\int_{\mathrm{SO}(2)}\int_{[0,1]^{2}}|\mathcal{I}(f;\eta,\boldsymbol{\xi})-\mathcal{I}(f)|^{2}d\boldsymbol{\xi}d\eta$
$\displaystyle=\int_{\mathrm{SO}(2)}\sum_{\mathbf{w}\in\mathbb{Z}^{2},\mathbf{w}\neq{\mathbf{0}}}|\widehat{\mathcal{Q}}_{\Gamma}(h^{-1}\eta^{-1}\mathbf{w})|^{2}d\eta.$
Next, we provide a slightly different estimate of
$\widehat{\mathcal{Q}}_{\Gamma}$ from the one in the proof of Lemma 2.2.
Without loss of generality, we choose a partition of unity
$\\{\phi_{j}\\}_{j=1}^{N}$ for $\Gamma$ such that each support
$\operatorname{supp}\phi_{j}\subset\Gamma$ contains exactly one point with
vanishing curvature to order $\kappa-2$. The case where
$\operatorname{supp}\phi_{j}$ does not contain any point with vanishing
curvature is already handled in Theorem 2.9. Locally, the support of
$\phi_{j}$ is represented by suppϕ_j = {(x, g_j(x)) ∣x∈(-r, r)} , with
$g_{j}(x)$ taking the form of (23). Let
$\mathbf{w}=|\mathbf{w}|(\sin\delta,\cos\delta)$. Then, according to
(LABEL:EQ:_STA_PHASE), if we set
$\mathbf{x}=(x(\delta),g_{j}(\delta))\in\Gamma$, then
$\mathbf{x}+s\mathbf{n}(\mathbf{x})$ becomes a stationary point if
$\sin(\delta)+g^{\prime}_{j}(x(\delta))\cos(\delta)-\frac{sg^{\prime\prime}_{j}(x(\delta))}{(1+|g^{\prime}_{j}(x(\delta))|^{2}))^{3/2}}\left(\sin(\delta)+g^{\prime}_{j}(x(\delta))\cos(\delta)\right)=0\,.$
(25) When $\varepsilon$ is small enough, (25) is equivalent to
$\sin\delta+g^{\prime}_{j}(x(\delta))\cos\delta=0\,.$ This leads to
$|x(\delta)|=\Theta(|\delta|^{1/(\kappa-1)})$ since locally
$g^{\prime}_{j}(x)=\Theta(|x|^{\kappa-1})$. We estimate
$\widehat{\mathcal{Q}}_{\Gamma}$ in three different cases:

  1. 1. If $|x(\delta)|>r$, then there are no stationary points near the point with vanishing curvature. The estimate will be the same as that in Lemma 2.2, that is, $|\widehat{\mathcal{Q}_{\Gamma}\phi_{j}}(\mathbf{w})|=\mathcal{O}(\varepsilon^{-q-1}|\mathbf{w}|^{-1/2-(q+1)})\,.$
  2. 2. If $|x(\delta)|<|\mathbf{w}|^{-\frac{\kappa-2}{\kappa(\kappa-1)}}$, the point is close to the origin. We can then use Lemma LABEL:LEM:_van_der_Corput_Revised) (the revised van der Corput Lemma) and Lemma 2.2 to conclude that $|\widehat{\mathcal{Q}_{\Gamma}\phi_{j}}(\mathbf{w})|=\mathcal{O}(\varepsilon^{-q-1}|\mathbf{w}|^{-1/\kappa-(q+1)})\,.$
  3. 3. If $r\geq|x(\delta)|\geq|\mathbf{w}|^{-\frac{\kappa-2}{\kappa(\kappa-1)}}$, the point’s Hessian is at order $\mathcal{O}(|x(\delta)|^{\kappa-2})$. We apply the standard stationary phase approximation with the additional factor $|x(\delta)|^{-(\kappa-2)/2}=\Theta(|\delta|^{-\frac{\kappa-2}{2\kappa-2}})$ to get $|\widehat{\mathcal{Q}_{\Gamma}\phi_{j}}(\mathbf{w})|=\mathcal{O}(\varepsilon^{-q-1}|\mathbf{w}|^{-1/2-(q+1)}|\delta|^{-\frac{\kappa-2}{2\kappa-2}})\,.$

Therefore, we have
$\displaystyle\int_{\mathrm{SO(2)}}|\widehat{\mathcal{Q}_{\Gamma}\phi_{j}}(h^{-1}\eta^{-1}\mathbf{w})|^{2}d\eta$
$\displaystyle=\mathcal{O}\left(\left|\varepsilon^{-q-1}|h^{-1}\mathbf{w}|^{-1/2-(q+1)}\right|^{2}\right)$
$\displaystyle\quad+\mathcal{O}\left(|h^{-1}\mathbf{w}|^{-\frac{\kappa-2}{\kappa}}\left|\varepsilon^{-q-1}|h^{-1}\mathbf{w}|^{-1/\kappa-(q+1)}\right|^{2}\right)$
$\displaystyle\quad+\int_{|h^{-1}\mathbf{w}|^{-\frac{\kappa-2}{\kappa}}}^{\mathcal{O}(1)}\mathcal{O}\left(\left|\varepsilon^{-q-1}|h^{-1}\mathbf{w}|^{-1/2-(q+1)}|\delta|^{-\frac{\kappa-2}{2\kappa-2}}\right|^{2}\right)d\delta\,.$
The first two terms are both bounded by
$\mathcal{O}\left(\varepsilon^{-2q-2}h^{1+2(q+1)}|\mathbf{w}|^{-1-2(q+1)}\right)$,
and the last integral is bounded by
$\mathcal{O}(\varepsilon^{-2q-2}|\mathbf{w}|^{-1-2(q+1)}h^{1+2(q+1)})\int_{|h^{-1}\mathbf{w}|^{-\frac{\kappa-2}{\kappa}}}^{r}|\delta|^{-\frac{\kappa-2}{\kappa-1}}d\delta=\mathcal{O}(\varepsilon^{-2q-2}|\mathbf{w}|^{-1-2(q+1)}h^{1+2(q+1)})\,.$
Taking $\varepsilon=\mathcal{O}(h^{\alpha})$ then gives us
$\displaystyle\int_{\mathrm{SO}(2)}\sum_{\mathbf{w}\in\mathbb{Z}^{2},\mathbf{w}\neq{\mathbf{0}}}|\widehat{\mathcal{Q}_{\Gamma}\phi_{j}}(h^{-1}\eta^{-1}\mathbf{w})|^{2}d\eta$
$\displaystyle=\sum_{\mathbf{w}\in\mathbb{Z}^{2},\mathbf{w}\neq{\mathbf{0}}}\mathcal{O}(\varepsilon^{-2q-2}|\mathbf{w}|^{-1-2(q+1)}h^{1+2(q+1)})$
$\displaystyle=\mathcal{O}(h^{1+2(q+1)(1-\alpha)})\,.$ The last step is to use
the Cauchy-Schwartz inequality and the finiteness of $N$ to conclude that
$\displaystyle\int_{\mathrm{SO}(2)}\sum_{\mathbf{w}\in\mathbb{Z}^{2},\mathbf{w}\neq{\mathbf{0}}}|\widehat{\mathcal{Q}_{\Gamma}}(h^{-1}\eta^{-1}\mathbf{w})|^{2}d\eta$
$\displaystyle=\int_{\mathrm{SO}(2)}\sum_{\mathbf{w}\in\mathbb{Z}^{2},\mathbf{w}\neq{\mathbf{0}}}|\sum_{j=1}^{N}\widehat{\mathcal{Q}_{\Gamma}\phi_{j}}(h^{-1}\eta^{-1}\mathbf{w})|^{2}d\eta$
$\displaystyle\leq
N\int_{\mathrm{SO}(2)}\sum_{\mathbf{w}\in\mathbb{Z}^{2},\mathbf{w}\neq{\mathbf{0}}}|\widehat{\mathcal{Q}_{\Gamma}\phi_{j}}(h^{-1}\eta^{-1}\mathbf{w})|^{2}d\eta$
$\displaystyle=\mathcal{O}(h^{1+2(q+1)(1-\alpha)})\,.$ The proof is complete.
∎ In fact, one can work a little harder to remove the convexity requirement in
the above theorem. The main idea is that, instead of using the representation
(23) for the curve, which sustains the convexity, we can represent the curve
locally as $(x,g(x))\subset\mathbb{R}^{2}$ for $x\in(-r,r)$ that
$g(x)=\operatorname{sgn}(x)|x|^{\kappa}h(x),\quad h(x)\neq 0\quad\forall
x\in(-r,r),\quad h\in C^{2}(-r,r)\,.$ (26) The proof of Lemma
LABEL:LEM:_van_der_Corput_Revised (the revised van der Corput Lemma) shows
both (23) and (26) can be handled by the same procedure; see Remark
LABEL:RMK:Nonconvex. We can thus have the following corollary.

######  Corollary 3.3.

The same estimate in Theorem 3.2 holds for any closed boundary $\Gamma\in
C^{\infty}$ if it only has finitely many points with vanishing curvature,
whose maximal order is $(\kappa-2)$ for $\kappa\in(2,\infty)$.

###  3.3 Numerical experiments

Here, we present some numerical simulations to verify the variance estimates
in Theorem 3.2 and Corollary 3.3 for the quadrature error under random rigid
transformations.

####  3.3.1 Convex curve with vanishing curvatures

To numerically verify the variance estimate in Theorem 3.2, we consider the
convex curve $\frac{(x-x_{0})^{4}}{r^{4}}+\frac{(y-y_{0})^{2}}{r^{2}}=1$ (27)
with $r=\frac{3}{4}$ and a random center $(x_{0},y_{0})$. The curvature
vanishes to order two at $(x_{0},y_{0}\pm r)$; See Figure 11. The integrand
function is selected as $f(x,y)=\cos(x^{2}-y)\sin(y^{2}-x^{3}).$ Figure 11:
The smooth boundary with vanishing curvatures to order two in (27). The weight
function is $\theta_{\varepsilon}^{\Delta}\in\mathcal{W}_{1}$. Each experiment
is performed with 32 random rigid transformations independently. The decay
rates of $\mathcal{O}(h^{5-4\alpha})$ and $\mathcal{O}(h^{7-6\alpha})$ have
been observed in Figure 12 and Figure 13 for different tube widths
$\varepsilon=\Theta(h^{\alpha})$, respectively. Figure 12: Variance of
quadrature error for weight function $\theta_{\varepsilon}^{\Delta}$. Left:
$\varepsilon=2h$. Middle: $\varepsilon=2h^{\frac{1}{2}}$. Right:
$\varepsilon=0.1$. Figure 13: Variance of quadrature error for weight function
$\theta_{\varepsilon}^{\cos}$. Left: $\varepsilon=2h$. Middle:
$\varepsilon=2h^{\frac{1}{2}}$. Right: $\varepsilon=0.1$.

####  3.3.2 Star shape

For the case of non-convex curves, we take the star-shaped curve in polar
coordinates: $\rho(\theta)=R+r\cos(m\theta)$ (28) with parameters $R=0.75$,
$r=0.2$, and $m=3$; see Figure 14. Figure 14: Star-shape boundary with
parameter $R=0.75$, $r=0.2$, and $m=3$ in (28). The curve is non-convex but
consists of finitely many points with vanishing curvatures. We use it to
verify the result in Corollary 3.3. The integrand function is selected as
$f(x,y)=\cos(x^{2}-y)\sin(y^{2}-x^{3}).$ The weight function is
$\theta_{\varepsilon}^{\Delta}\in\mathcal{W}_{1}$. Each experiment is
performed with 32 random rigid transformations independently. The decay rate
of $\mathcal{O}(h^{5-4\alpha})$ is observed in Figure 15. Figure 15: Variance
of quadrature error for weight function $\theta_{\varepsilon}^{\Delta}$. Left:
$\varepsilon=2h$. Middle: $\varepsilon=2h^{\frac{1}{2}}$. Right:
$\varepsilon=0.1$.

##  4 Further discussions on open curves

In this section, we consider the quadrature error on open curves instead of
closed ones. For technical reasons in the derivation of the theory, mainly the
fact that in the case of open curves, we cannot construct a partition of unity
with $C_{0}^{\infty}$ functions, we have only weaker results compared to the
ones in Section 2 and Section 3. We consider two groups of curves: (i) curves
with finitely many inflection points (curvature zero) and (ii) line segments.

###  4.1 General curves

Let $\Gamma$ be a collection of curves that
$\Gamma=\bigcup_{k=1}^{m}\Gamma_{k}$, where each curve $\Gamma_{k}\in
C^{\infty}$ has only a finite number of points with vanishing curvature, whose
maximal order is $(\kappa-2)$, $\kappa\in(2,\infty)$. Let
$\mathcal{Q}_{k}(\mathbf{x}):=f(P_{\Gamma_{k}}\mathbf{x})\theta_{\varepsilon}(d_{\Gamma_{k}}(\mathbf{x}))J_{\varepsilon}(\mathbf{z},d_{\Gamma_{k}}(\mathbf{x}))$,
we define generalized implicit boundary integral on the collection $\Gamma$ as
$\displaystyle\mathcal{I}_{h}(f):=h^{2}\sum_{k=1}^{m}\sum_{\mathbf{z}\in(h\mathbb{Z})^{2}\cap
T_{k,\varepsilon}}\mathcal{Q}_{k}(\mathbf{z}),$ (29) where
$T_{k,\varepsilon}=\\{\mathbf{x}+t\mathbf{n}_{k}(\mathbf{x})\mid(\mathbf{x},t)\in\Gamma_{k}\times[-\varepsilon,\varepsilon]\\}$
is the segmented tube associated with $\Gamma_{k}$ and
$\mathbf{n}_{k}(\mathbf{x})$ is the unit normal vector at
$\mathbf{x}\in\Gamma_{k}$. Clearly, the smooth closed boundary in Corollary
3.3 is a special case for $m=1$.

######  Theorem 4.1.

For any rotation $\eta\in\mathrm{SO}(2)$ and translation
$\boldsymbol{\xi}\in[0,1]^{2}$, we denote the generalized implicit boundary
integral with tube width $\varepsilon=\Theta(h^{\alpha})$ on the transformed
curve collection $\eta\Gamma+h\boldsymbol{\xi}$ as
$\mathcal{I}_{h}(f;\eta,\boldsymbol{\xi})$, then if there are no points with
vanishing curvature or there exist finitely many points with vanishing
curvature to the maximal order $\kappa-2$ for
$\kappa\in[2,\frac{3+\sqrt{5}}{2})$, then ∫_SO(2)∫_[0, 1]^2 —I(f) - I_h(f; η,
ξ)—^2 dξdη= O(h^2κ + 2κ- 1(1-α) ).

###### Proof.

Without loss of generality, we may assume $\Gamma_{k}$ has at most one point with vanishing curvature at its endpoint, otherwise one can split the curve into shorter pieces. Let Γ_k ≔{(x, g_k(x)) ∣x∈[0, r]} such that $g_{k}(x)=|x|^{\kappa}h(x)$, $h\in C^{\infty}[0,r]$ and $h(x)>0$ for all $x\in[0,r]$. When $\kappa=2$, then there are no points with vanishing curvature. Using the derivation in Appendix LABEL:PRF:_STATION, the Fourier transform of $\mathcal{Q}_{k}\chi_{T_{k,\varepsilon}}$ becomes $\displaystyle\widehat{\mathcal{Q}_{k}\chi_{T_{k,\varepsilon}}}(\boldsymbol{\zeta})$ $\displaystyle=\int_{T_{k,\varepsilon}}\mathcal{Q}_{k}(\mathbf{x})e^{-2\pi i\mathbf{x}\cdot\boldsymbol{\zeta}}d\mathbf{x}$ (30) $\displaystyle=\int_{-\varepsilon}^{\varepsilon}\theta_{\varepsilon}(s)\int_{\Gamma_{k}}f(\mathbf{x}^{\prime})e^{-2\pi i(\mathbf{x}^{\prime}+s\mathbf{n}(\mathbf{x}^{\prime}))\cdot\boldsymbol{\zeta}}d\sigma ds$ (30) $\displaystyle=\int_{-\varepsilon}^{\varepsilon}\theta_{\varepsilon}(s)\int_{0}^{r}f(x)e^{-2\pi i\phi(x,\zeta_{1},\zeta_{2},s)}\sqrt{1+|g^{\prime}_{k}(x)|^{2}}dxds$ (30) where $d\sigma$ is the Lebesgue measure on $\Gamma_{k}$, $\chi_{T_{k,\varepsilon}}$ denotes the characteristic function on $T_{k,\varepsilon}$, and $\phi(x,\zeta_{1},\zeta_{2},s)=x\zeta_{1}+g_{k}(x)\zeta_{2}+s\frac{-g^{\prime}_{k}(x)\zeta_{1}+\zeta_{2}}{\sqrt{1+|g^{\prime}_{k}(x)|^{2}}}\,.$ Let $\\{\chi_{l}\\}_{l\geq 1}$ be a partition of unity over the interval $[0,r]$. It is clear that we only have to estimate $\widehat{\mathcal{Q}_{k}\chi_{T_{k,\varepsilon}}}$ for each integrand $f(x)\chi_{l}$ instead of $f(x)$. Therefore, we may further assume $f(r)=0$. Let us define the auxiliary function $G(x,s)$ that $G(x,s)=\int_{0}^{x}e^{-2\pi i\phi(t,\zeta_{1},\zeta_{2},s)}dt\,.$ Then, using integration by parts on (30), as well as the the assumption that $f(r)=0$, we obtain $\displaystyle\widehat{\mathcal{Q}_{k}\chi_{T_{k,\varepsilon}}}(\boldsymbol{\zeta})$ $\displaystyle=\int_{-\varepsilon}^{\varepsilon}\theta_{\varepsilon}(s)\int_{0}^{r}f(x)e^{-2\pi i\phi(x,\zeta_{1},\zeta_{2},s)}\sqrt{1+|g^{\prime}_{k}(x)|^{2}}dxds$ (31) $\displaystyle=\int_{-\varepsilon}^{\varepsilon}\theta_{\varepsilon}(s)\left[-\int_{0}^{r}G(x,s)\frac{d}{dx}\left(f(x)\sqrt{1+|g^{\prime}_{k}(x)|^{2}}\right)dx\right]ds\,.$ (31) Let $\mathbf{x}^{\prime}=(x,g_{k}(x))$ and $\mathbf{n}(\mathbf{x}^{\prime})=(\sin\psi(x),\cos\psi(x))$, then $\psi(x)=-\arctan(g^{\prime}(x))=\Theta(|x|^{\kappa-1})$ and $\psi\in[-\psi_{0},0]$ for certain $\psi_{0}\in(0,\frac{\pi}{2})$. In the following, we use the notation $\boldsymbol{\zeta}=|\boldsymbol{\zeta}|(\sin\delta,\cos\delta)$ and estimate the integral in the following form $\mathcal{E}:=\int_{-\varepsilon}^{\varepsilon}\theta_{\varepsilon}(s)\int_{0}^{r}G(x,s)\frac{d}{dx}\left(f(x)\sqrt{1+|g^{\prime}_{k}(x)|^{2}}\right)dxds.$ (32) We split the integral in $x$ into two different regimes. |
Table 1: Error estimates for implicit boundary integral method with
$\varepsilon=\Theta(h^{\alpha})$, $\alpha\in[0,1]$ on strongly convex
surfaces.