# A Second-Order Nonlocal Approximation to Manifold Poisson Models with
Neumann Boundary ††thanks: This work was supported by NSFC grant 12071244,
11671005.

Yajie Zhang Department of Statistics and Mathematics, Zhongnan University of
Economics and Laws, Wuhan, China, 430000 Email<EMAIL_ADDRESS>Zuoqiang
Shi Corresponding Author, Department of Mathematical Sciences, Tsinghua
University, Beijing, China, 100084\. Email<EMAIL_ADDRESS>
###### Abstract

In this paper, we propose a class of nonlocal models to approximate the
Poisson model on manifolds with homogeneous Neumann boundary condition, where
the manifolds are assumed to be embedded in high dimensional Euclid spaces. In
comparison to the existing nonlocal approximation of Poisson models with
Neumann boundary, we optimize the truncation error of model by adding an
augmented term along the $2\delta$ layer of boundary, with $2\delta$ be the
nonlocal interaction horizon. Such term is formulated by the integration of
the second order normal derivative of solution through the boundary, while the
second order normal derivative is expressed as the difference between the
interior Laplacian and the boundary Laplacian. The concentration of our paper
is on the construction of nonlocal model, the well-posedness of model, and its
second-order convergence rate to its local counterpart. The localization rate
of our nonlocal model is currently optimal among all related works even for
the case of high dimensional Euclid spaces.

###### keywords:

Manifold Poisson equation, Neumann boundary, nonlocal approximation, well-
posedness, second order convergence.

45P05; 45A05; 35A15; 46E35

##  1 Introduction

Partial differential equations on manifolds have never failed to attract
researchers in the past decades. Its application includes material science [9]
[19] , fluid flow [21] [23], and biology physics [6] [20] [31]. In addition,
many recent work on machine learning [7] [11] [26] [29] [35] and image
processing [10] [22] [24] [25] [32] [34] [44] relates their model with
manifold PDEs. One frequently used manifold model is the manifold Poisson
equation. Apart from its broad application, such equation is mathematically
interesting by itself since it reveals much information of the manifold. Among
all the recent study on the numerical analysis of the Poisson models, one
approach is its nonlocal approximation. The main advantage for nonlocal model
is that the use of explicit spatial differential operators is avoided, hence
new numerical schemes can be explored, for instance, the point integral
method(PIM). In spite of the large quantity of articles on the nonlocal
approximation of Poisson model on Euclid spaces, less work have been done on
its extension from Euclid spaces into manifolds. Due to the demand of
numerical methods for manifold PDEs, it is necessary to construct one nonlocal
model that approximates the manifold Poisson equation with high accuracy,
hence being able to solve it by proper numerical scheme.

In the past decade, researchers devoted on the nonlocal approximation of local
models on Euclid domains, among which the most commonly studied equation is

$\frac{1}{\delta^{2}}\int_{\Omega}(u_{\delta}({x})-u_{\delta}({y}))R_{\delta}({x},{y})d{y}=f({x}),\qquad{x}\in\Omega.$
(1.1)

Here $\Omega\subset\mathbb{R}^{k}$ is a bounded Euclid domain with smooth
boundary, $f\in H^{2}(\Omega)$, $\delta$ is the nonlocal horizon parameter
that describes the range of nonlocal interaction,
$R_{\delta}({x},{y})=C_{\delta}R\big{(}\frac{|{x}-{y}|^{2}}{4\delta^{2}}\big{)}$
is the nonlocal kernel function, where $R\in C^{2}(\mathbb{R}^{+})\cap
L^{1}[0,\infty)$ is a properly-chosen non-negative function with compact
support, and $C_{\delta}=\frac{1}{(4\pi\delta^{2})^{k/2}}$ is the
normalization factor. Such equation usually appeared in the discussion of
peridynamics models [3] [8] [13] [33] [40] [42]. The truncation error of (1.1)
to its local counterpart $\Delta u=f$ has been proved to be
$\mathcal{O}(\delta^{2})$ in the interior region away from boundary, and
$\mathcal{O}(\delta^{-1})$ in the $2\delta-$layer along the boundary. Efforts
have been made to reduce the order of truncation error by making modification
of (1.1) along the boundary layer. See [4] [5] [12] [14] for Neumann boundary
condition and [1] [2] [15] [30] [49] for other types of boundary conditions.
Other related modification of (1.1) involves fractional nonlocal kernel [17]
or nonlocal gradient terms [16]. Those studies yield to at most
$\mathcal{O}(\delta)$ convergence rate from $u_{\delta}$ to $u$.

As a breakthrough, in one dimensional [41] and two dimensional [45] cases, the
nonlocal models with $\mathcal{O}(\delta^{2})$ convergence rate to its local
counterpart were successfully constructed under Neumann boundary condition.
One year later, Lee H. and Du Q. in [27] introduced a nonlocal model under
Dirichlet boundary condition by imposing a special volumetric constraint along
the boundary layer, which assures $\mathcal{O}(\delta^{2})$ convergence rate
in 1d segment and 2d plain disk. In 2022, Yin X. in [18] introduced a nonlocal
Poisson model with polygonally-supported interaction kernel instead of the
classical kernels with fixed horizon, which brings much convenience on the
numerical simulation of nonlocal models.

Nevertheless, in comparison to the case of Euclid domains, less work have been
done on the nonlocal approximation of Poisson model on manifolds. Such
deficiency restricted the exploration of new numerical methods for manifold
PDEs. In 2018, nonlocal Poisson model on manifolds were first studied in [38]
in the case of Neumann boundary. Its subsequent works describe the manifold
nonlocal models with Dirichlet boundary [36] [43] , with interface [46] , and
with anisotropic coefficients [37]. As a representative of these works, in
[38], a nonlocal Poisson model was constructed by utilizing the following
approximation:

$\\!\\!\\!\\!\\!\\!\\!\\!-\int_{\mathcal{M}}\Delta_{\mathcal{M}}u({y})\bar{R}_{\delta}({x},{y})d\mu_{{y}}\approx\int_{\mathcal{M}}\frac{1}{\delta^{2}}R_{\delta}({x},{y})(u({x})-u({y}))d\mu_{{y}}-2\int_{\partial\mathcal{M}}\bar{R}_{\delta}({x},{y})\frac{\partial
u}{\partial\bm{n}}({y})d\tau_{y},$ (1.2)

Here $\mathcal{M}$ is a compact, smooth $m$ dimensional manifold embedded in
$\mathbb{R}^{d}$, with $\partial\mathcal{M}$ a smooth $(m-1)$ dimensional
submanifold. $\Delta_{\mathcal{M}}$ is the Laplace-Beltrami operator on
$\mathcal{M}$. We omit the definition of $\Delta_{\mathcal{M}}$ and
$\nabla_{\mathcal{M}}$, see page $2-3$ of [38] for reference. $u\in
H^{4}(\mathcal{M})$, $\frac{\partial
u}{\partial\bm{n}}=\nabla_{\mathcal{M}}u\cdot\bm{n}$, $d\mu_{y}$ and
$d\tau_{y}$ are the volume form of $\mathcal{M}$ and $\partial\mathcal{M}$.
The kernel function $R_{\delta}$ has the same definition as $R_{\delta}$ in
(1.1), while
$\bar{R}_{\delta}({x},{y})=C_{\delta}\bar{R}\big{(}\frac{|{x}-{y}|^{2}}{4\delta^{2}}\big{)}$
and $\bar{R}(r)=\int_{r}^{+\infty}R(s)ds$. It can be derived that
$2\delta^{2}\nabla^{{x}}_{\mathcal{M}}\bar{R}_{\delta}({x},{y})\approx({y}-{x})R_{\delta}({x},{y})$
for ${x},{y}\in\mathcal{M}$.

The error of (1.2) has been rawly analyzed in [38], and the nonlocal model
with Neumann boundary was set up by eliminating the last term of (1.2) and
equalizing both sides. Such model converges to its local counterpart in a rate
of $\mathcal{O}(\delta)$. The same rate was reached in [28] in the case of
Dirichlet boundary. In such model, the Dirichlet boundary was approximated by
Robin condition so that the term $\frac{\partial u}{\partial\bm{n}}$ in (1.2)
can be replaced by $\lambda u$ for some constant $\lambda$.

Due to the demand of high accuracy for nonlocal approximation, one natural
question to consider is the existence of nonlocal Poisson model that
approaches its local counterpart in a rate of $\mathcal{O}(\delta^{2})$. In
our previous work [47] [48], such rate is attained in the case of Dirichlet
boundary. The crucial part in the construction of model is that we improve the
accuracy of approximation (1.2) by adding a term involving the second-order
normal derivative of $u$ along the boundary:

$\begin{split}\\!\\!\\!\\!\\!\\!\\!\\!\\!-\int_{\mathcal{M}}\Delta_{\mathcal{M}}u({y})\bar{R}_{\delta}({x},{y})d\mu_{y}=\
&\frac{1}{\delta^{2}}\int_{\mathcal{M}}(u({x})-u({y}))R_{\delta}({x},{y})d\mu_{y}-2\int_{\partial\mathcal{M}}\frac{\partial
u}{\partial\bm{n}}({y})\bar{R}_{\delta}({x},{y})d\tau_{y}\\\
&-\int_{\partial\mathcal{M}}({x}-{y})\cdot\bm{n}({y})\frac{\partial^{2}u}{\partial^{2}\bm{n}}({y})\bar{R}_{\delta}({x},{y})d\tau_{y}+r_{in}({x}),\
{x}\in\mathcal{M},\end{split}$ (1.3)

where $r_{in}$ is an $\mathcal{O}(\delta^{2})$ function and has been given a
complete control in Theorem 3.1 of [48]. Based on (1.3), we constructed a
nonlocal Poisson model in [47] with $\mathcal{O}(\delta^{2})$ convergence rate
to its local counterpart in the case of Dirichlet boundary, where the terms
involving the normal derivatives of $u$ were approximated by some integration
of $u$ along the boundary layer.

In the case of Neumann boundary, since no information is given regarding the
second order normal derivative, it is difficult to convert (1.3) into a
nonlocal model. In [45] under the case of 2d Euclid domain, such derivative is
expressed by the integration of $u$ through level curves that are parallel to
the boundary curve. The corresponding model is well-posed and has
$\mathcal{O}(\delta^{2})$ localization rate. Nevertheless, such model cannot
be extended into the case when the domain is a high dimensional manifold.

In this paper, to overcome the issue regarding the second order normal
derivative, we introduce an augmented function that approximates the value of
$u$ along the boundary to guarantee the coercivity of model. Such derivative
can be approximated by the derivative of the augmented function, which equals
to the difference between the interior Laplacian and boundary Laplacian
according to our previous result [48]. Since the first order normal derivative
vanishes, we are able to rewrite the approximation (1.3) and set up a nonlocal
model accordingly. The main contribution of this work is the construction of
nonlocal model, the well-posedness of model and its second order convergence
rate to its local counterpart. Our results can be easily generalized into the
case with non-homogeneous Neumann boundary condition. To the author’s best
knowledge, even in Euclid spaces with Neumann boundary, no work has ever
appeared on the construction of nonlocal Poisson model with
$\mathcal{O}(\delta^{2})$ convergence under dimension $d\geq 3$, having such
model will result in higher efficiency in the numerical implementation. In
addition, our model can be easily discretized and be solved by meshless scheme
like point integral method(PIM, see [28]) , which brings more convenience on
solving manifold models numerically. On the contrary, we can not apply the
finite-element method(FEM) on high dimensional manifold due to the
complication of mesh construction. Besides, this work gives an orientation on
how to handle nonlocal problems on manifolds.

The paper is organized as follows: we introduce the idea on the formulation of
nonlocal Poisson model in section 2. In section 3, we present the well-
posedness of model and prove the solution to be in a more regular Hilbert
space. The convergence rate of model to its local counterpart is given in
section 4. In section 5, we make a conclusion on our paper and discuss our
future work. The control of an error function have been moved into appendix
due to its complication.

##  2 Construction of Nonlocal Model

In this work, our purpose is to find a second-order nonlocal approximation to
the following Poisson problem:

$\begin{cases}-\Delta_{\mathcal{M}}u({x})=f({x})&{x}\in\mathcal{M},\\\
\frac{\partial u}{\partial\bm{n}}({x})=0&{x}\in\partial\mathcal{M},\\\
\int_{\mathcal{M}}u({x})d\mu_{{x}}=0;\end{cases}$ (2.1)

here $\mathcal{M}$ is a compact, smooth $m$ dimensional manifold embedded in
$\mathbb{R}^{d}$, with $\partial\mathcal{M}$ a smooth $(m-1)$ dimensional
submanifold. $\Delta_{\mathcal{M}}$ is the Laplace-Beltrami operator on
$\mathcal{M}$. $f$ is an $H^{2}$ function on $\mathcal{M}$ with
$\int_{\mathcal{M}}f({x})d\mu_{{x}}=0$.

In order to construct a nonlocal model, we will first introduce the
interaction kernel function
$R_{\delta}({x},{y})=C_{\delta}R\big{(}\frac{|{x}-{y}|^{2}}{4\delta^{2}}\big{)}$
that appeared in (1.2). A better choice of $R$ will result in much less work
in the error analysis. To this end, we make the following assumptions on
$R(r)$:

  1. 1.

Smoothness: $R\in W^{2,+\infty}[0,+\infty)$, i.e., for any $r\geq 0$ we have
$\big{|}D^{2}R(r)\big{|}\leq C$;

  2. 2.

Nonnegativity: $R(r)>0$ for any $r\geq 0$;

  3. 3.

Compact support: $R(r)=0$ for any $r>1$;

  4. 4.

Nondegenearcy: $\exists\ \delta_{0}>0$ so that $R(r)\geq\delta_{0}>0$ for
$0\leq r\leq 1/2$.

These assumptions are very mild, almost all the smooth functions in the
literature satisfy these conditions. The first condition makes $R$ not blow up
near the origin, while the fourth condition assures a non-degenerate interval
for $R$. Despite that the removal of third condition will make the problem
more interesting, we keep this condition to avoid complication in the error
analysis. A frequently used kernel function $R$ is the following:

$R(r)=\begin{cases}\frac{1}{2}(1+\cos\pi r),&0\leq r\leq 1,\\\
0,&r>1.\end{cases}$ (2.2)

Based on $R$, we write down the other three related functions
${R}(r)=\int_{r}^{+\infty}\underline{R}(s)ds$,
$\bar{R}(r)=\int_{r}^{+\infty}R(s)ds$,
$\overset{=}{R}(r)=\int_{r}^{+\infty}\bar{R}(s)ds$ and their corresponding
kernel functions

$\underline{R}_{\delta}({x},{y})=C_{\delta}\underline{R}\big{(}\frac{|{x}-{y}|^{2}}{4\delta^{2}}\big{)},\
\bar{R}_{\delta}({x},{y})=C_{\delta}\bar{R}\big{(}\frac{|{x}-{y}|^{2}}{4\delta^{2}}\big{)},\
\overset{=}{R}_{\delta}({x},{y})=C_{\delta}\overset{=}{R}\big{(}\frac{|{x}-{y}|^{2}}{4\delta^{2}}\big{)}.$
(2.3)

It is clear that $\underline{R}$, $\bar{R}$ and $\overset{=}{R}$ satisfy the
assumptions 1-4 as well by replacing $\delta_{0}$ with some other positive
constants in the fourth condition.

Now we start the derivation of nonlocal model. Let $u$ be the unique solution
of (2.1), then according to the boundary condition $\frac{\partial
u}{\partial\bm{n}}\equiv 0$, we can rewrite (1.3) into

$\begin{split}\frac{1}{\delta^{2}}\int_{\mathcal{M}}(u({x})-u({y}))R_{\delta}({x},{y})d\mu_{y}-\int_{\partial\mathcal{M}}({x}-{y})\cdot\bm{n}({y})\frac{\partial^{2}u}{\partial^{2}\bm{n}}({y})\bar{R}_{\delta}({x},{y})d\tau_{y}\\\
=\int_{\mathcal{M}}f({y})\bar{R}_{\delta}({x},{y})d\mu_{y}+r_{1}({x}),\qquad{x}\in\mathcal{M};\end{split}$
(2.4)

where the error function $r_{1}$ is exactly $r_{in}$ in the equation (1.3).
For the term $\frac{\partial^{2}u}{\partial^{2}\bm{n}}$, the equation $(4.15)$
of [48], gives:

$\begin{split}&\Delta_{\mathcal{M}}u({x})=\sum\limits_{i=1}^{m-1}\bm{H}_{u}({x}):\bm{\tau_{i}}\otimes\bm{\tau_{i}}+\frac{\partial^{2}u}{\partial\bm{n}^{2}}({x})=\Delta_{\partial\mathcal{M}}u({x})+\frac{\partial^{2}u}{\partial\bm{n}^{2}}({x}),\qquad{x}\in\partial\mathcal{M},\end{split}$
(2.5)

where $\\{\bm{\tau}_{i}\\}_{i=1}^{m-1}$ is an orthonormal basis of the tangent
space of $\partial\mathcal{M}$ at ${x}$. We then apply (2.5) to (2.4) to
obtain

$\begin{split}\frac{1}{\delta^{2}}\int_{\mathcal{M}}(u({x})-u({y}))R_{\delta}({x},{y})d\mu_{y}+\int_{\partial\mathcal{M}}\Delta_{\partial\mathcal{M}}u({y})({x}-{y})\cdot\bm{n}({y})\bar{R}_{\delta}({x},{y})d\tau_{y}\\\
=\int_{\mathcal{M}}f({y})\bar{R}_{\delta}({x},{y})d\mu_{y}+\int_{\partial\mathcal{M}}({x}-{y})\cdot\bm{n}({y})f({y})\bar{R}_{\delta}({x},{y})d\tau_{y}+r_{1}({x}),\
{x}\in\mathcal{M}.\end{split}$ (2.6)

Now we aim to seek for a term to approximate the second term of (2.6) with
error $r_{2}$. Let the operator $P:L^{2}(\mathcal{M})\to
L^{2}(\partial\mathcal{M})$ be defined as

$\begin{split}&Pu({y})=\int_{\mathcal{M}}u({x})({x}-{y})\cdot\bm{n}({y})\bar{R}_{\delta}({x},{y})d\mu_{x},\qquad{y}\in\partial\mathcal{M},\end{split}$
(2.7)

and

$\begin{split}&\omega({y})=\int_{\mathcal{M}}({x}-{y})\cdot\bm{n}({y})\bar{R}_{\delta}({x},{y})d\mu_{x},\qquad{y}\in\partial\mathcal{M},\end{split}$
(2.8)

and we define the error function

$r_{2}({x})=\int_{\partial\mathcal{M}}(\Delta_{\partial\mathcal{M}}u({y})-\nabla_{\partial\mathcal{M}}\cdot\frac{\nabla_{\partial\mathcal{M}}Pu({y})}{\omega({y})})\
\bm{n}({y})\cdot({x}-{y})\bar{R}_{\delta}({x},{y})d\tau_{y},\qquad{x}\in\mathcal{M},$
(2.9)

then (2.6) can be written as

$\begin{split}&\frac{1}{\delta^{2}}\int_{\mathcal{M}}(u({x})-u({y}))R_{\delta}({x},{y})d\mu_{y}+\int_{\partial\mathcal{M}}(\nabla_{\partial\mathcal{M}}\cdot\frac{\nabla_{\partial\mathcal{M}}Pu({y})}{\omega({y})})\bm{n}({y})\cdot({x}-{y})\bar{R}_{\delta}({x},{y})d\tau_{y}\\\
&=\int_{\mathcal{M}}f({y})\bar{R}_{\delta}({x},{y})d\mu_{y}+\int_{\partial\mathcal{M}}({x}-{y})\cdot\bm{n}({y})f({y})\bar{R}_{\delta}({x},{y})d\tau_{y}+r_{1}({x})-r_{2}({x}),\qquad{x}\in\mathcal{M}.\end{split}$
(2.10)

Our next goal is to make the average of either side of (2.10) on $\mathcal{M}$
to be $0$. Let the functional $\bar{P}:L^{2}(\mathcal{M})\to\mathbb{R}$ be
defined as

$\bar{P}(u)=\frac{1}{|\mathcal{M}|}\int_{\mathcal{M}}\int_{\partial\mathcal{M}}(\nabla_{\partial\mathcal{M}}\cdot\frac{\nabla_{\partial\mathcal{M}}Pu({y})}{\omega({y})})\bm{n}({y})\cdot({x}-{y})\bar{R}_{\delta}({x},{y})d\tau_{y}d\mu_{{x}},$
(2.11)

and the functional $\bar{Q}:L^{2}(\mathcal{M})\to\mathbb{R}$ be defined as

$\bar{Q}(f)=\frac{1}{|\mathcal{M}|}\int_{\mathcal{M}}(\int_{\mathcal{M}}f({y})\bar{R}_{\delta}({x},{y})d\mu_{y}+\int_{\partial\mathcal{M}}({x}-{y})\cdot\bm{n}({y})f({y})\bar{R}_{\delta}({x},{y})d\tau_{y})d\mu_{x},$
(2.12)

then we rewrite (2.10) into

$\begin{split}&\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\frac{1}{\delta^{2}}\int_{\mathcal{M}}(u({x})-u({y}))R_{\delta}({x},{y})d\mu_{y}+\int_{\partial\mathcal{M}}(\nabla_{\partial\mathcal{M}}\cdot\frac{\nabla_{\partial\mathcal{M}}Pu({y})}{\omega({y})})\bm{n}({y})\cdot({x}-{y})\bar{R}_{\delta}({x},{y})d\tau_{y}-\bar{P}(u)\\\
=&\int_{\mathcal{M}}f({y})\bar{R}_{\delta}({x},{y})d\mu_{y}+\int_{\partial\mathcal{M}}({x}-{y})\cdot\bm{n}({y})f({y})\bar{R}_{\delta}({x},{y})d\tau_{y}-\bar{Q}(f)\\\
&+(\bar{Q}(f)-\bar{P}(u)+r_{1}({x})-r_{2}({x}))\qquad\forall\
{x}\in\mathcal{M}.\end{split}$ (2.13)

Finally, we remove the error terms
$(\bar{Q}(f)-\bar{P}(u)+r_{1}({x})-r_{2}({x}))$ in (2.13), and replace the
functions $u$ into $u_{\delta}$, to eventually establish the following
nonlocal model:

$\begin{split}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\frac{1}{\delta^{2}}\int_{\mathcal{M}}(u_{\delta}({x})-u_{\delta}({y}))R_{\delta}({x},{y})d\mu_{y}+\int_{\partial\mathcal{M}}(\nabla_{\partial\mathcal{M}}\cdot\frac{\nabla_{\partial\mathcal{M}}Pu_{\delta}({y})}{\omega({y})})\bm{n}({y})\cdot({x}-{y})\bar{R}_{\delta}({x},{y})d\tau_{y}-\bar{P}(u_{\delta})\\\
=\int_{\mathcal{M}}f({y})\bar{R}_{\delta}({x},{y})d\mu_{y}+\int_{\partial\mathcal{M}}({x}-{y})\cdot\bm{n}({y})f({y})\bar{R}_{\delta}({x},{y})d\tau_{y}-\bar{Q}(f),\qquad\forall\
{x}\in\mathcal{M};\end{split}$ (2.14)

where $P,\omega,\bar{P},\bar{Q}$ are defined in (2.7), (2.8), (2.11) and
(2.12).

Since the LHS of (2.14) is an operator on $u_{\delta}$, while the RHS of
(2.14) is an operator on $f$, we will frequently denote the nonlocal model
(2.14) as

$\mathcal{L}_{\delta}u_{\delta}({x})=f_{\delta}({x})\qquad\forall\
{x}\in\mathcal{M},$ (2.15)

for brevity’s sake in the following content. Next, we are ready to state our
main theorem in this paper.

######  Theorem 2.1.

  1. 1.

(Well-Posedness) Let $\mathcal{M}$ be a $C^{2}$ smooth, $m$ dimensional
manifold embedded in the Euclid space $\mathbb{R}^{n}$, while its boundary
$\partial\mathcal{M}$ be $C^{2}$ smooth, $(m-1)$ dimensional submanifold. Then
for any $f\in H^{1}(\mathcal{M})$ with $\int_{\mathcal{M}}f({x})d\mu_{{x}}=0$,
there is a unique function $u_{\delta}\in L^{2}(\mathcal{M})$ such that
$\int_{\mathcal{M}}u_{\delta}({x})d\mu_{x}=0$ and $u_{\delta}$ solves (2.14)
for $a.e.\ {x}\in\mathcal{M}$, with the estimate

$\left\|u_{\delta}\right\|_{H^{1}(\mathcal{M})}\leq
C\left\|f\right\|_{H^{1}(\mathcal{M})},$ (2.16)

where the constant $C$ depend only on the shape of $\mathcal{M}$ and
$\partial\mathcal{M}$.

  2. 2.

(Vanishing Nonlocality) If in addition to part 1, $f\in H^{2}(\mathcal{M})$,
and $u\in H^{4}(\mathcal{M})$ be the unique solution to the model (2.1), then
we have

$\left\|u-u_{\delta}\right\|_{H^{1}(\mathcal{M})}\leq
C\delta^{2}\left\|f\right\|_{H^{2}(\mathcal{M})},$ (2.17)

where the constant $C$ depend only on the shape of $\mathcal{M}$ and
$\partial\mathcal{M}$.

##  3 Well-Posedness

Our purpose here in this section is to prove the first part of Theorem 2.1. To
begin with, we denote the following function space

$L_{0}^{2}(\mathcal{M})=\\{q:\mathcal{M}\to\mathbb{R}|\ q\in
L^{2}(\mathcal{M})\mbox{ and }\int_{\mathcal{M}}q({x})d\mu_{{x}}=0\\}.$ (3.1)

Apparently, we can check that $f_{\delta}\in L_{0}^{2}(\mathcal{M})$ and
$\mathcal{L}_{\delta}q\in L_{0}^{2}(\mathcal{M})$ for all $q\in
L_{0}^{2}(\mathcal{M})$.

To prove the existence of solution $u_{\delta}\in L_{0}^{2}(\mathcal{M})$ to
(2.15), we need to employ Lax-Milgram theorem in the function space
$L_{0}^{2}(\mathcal{M})$. For all $p_{\delta},q_{\delta}\in
L_{0}^{2}(\mathcal{M})$, we define the following inner product

$(p_{\delta},q_{\delta})=\int_{\mathcal{M}}p_{\delta}({x})q_{\delta}({x})d\mu_{{x}};$
(3.2)

and the bilinear form

$\begin{split}&B_{\delta}[p_{\delta},q_{\delta}]=(\mathcal{L}_{\delta}p_{\delta},q_{\delta})=\frac{1}{\delta^{2}}\int_{\mathcal{M}}\int_{\mathcal{M}}q_{\delta}({x})(p_{\delta}({x})-p_{\delta}({y}))R_{\delta}({x},{y})d\mu_{y}d\mu_{x}\\\
&\\!\\!+\int_{\mathcal{M}}\int_{\partial\mathcal{M}}(\nabla_{\partial\mathcal{M}}\cdot\frac{\nabla_{\partial\mathcal{M}}Pp_{\delta}({y})}{\omega({y})})q_{\delta}({x})\bm{n}({y})\cdot({x}-{y})\bar{R}_{\delta}({x},{y})d\tau_{y}d\mu_{x}-\int_{\mathcal{M}}\bar{P}(p_{\delta})q_{\delta}({x})d\mu_{x}.\end{split}$
(3.3)

Since $f_{\delta}\in L_{0}^{2}(\mathcal{M})$, the nonlocal model (2.15) is
equivalent to the following weak formulation: find $u_{\delta}\in
L_{0}^{2}(\mathcal{M})$ such that

$B_{\delta}[u_{\delta},q_{\delta}]=(f_{\delta},q_{\delta})\qquad\forall\
q_{\delta}\in L_{0}^{2}(\mathcal{M}).$ (3.4)

Next, we present a lemma to illustrate that $B_{\delta}[\cdot,\cdot]$ satisfiy
the Lax-Milgram conditions.

######  Lemma 3.1.

For the bilinear form $B_{\delta}[\cdot,\cdot]$, we have

$B_{\delta}[p_{\delta},q_{\delta}]\leq
C_{\delta}\left\|p_{\delta}\right\|_{L^{2}(\mathcal{M})}\left\|q_{\delta}\right\|_{L^{2}(\mathcal{M})},\qquad\forall\
p_{\delta},q_{\delta}\in L_{0}^{2}(\mathcal{M});$ (3.5)

and

$B_{\delta}[p_{\delta},p_{\delta}]\geq
C\left\|p_{\delta}\right\|^{2}_{L^{2}(\mathcal{M})},\qquad\forall\
p_{\delta}\in L_{0}^{2}(\mathcal{M});$ (3.6)

for some constant $C$ depend only on the shape of $\mathcal{M}$, and some
constant $C_{\delta}$ depend on $\mathcal{M}$ and $\delta$.

######  Proof 3.2.

First of all, the formulation (3.3) gives

$\begin{split}&\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!B_{\delta}[p_{\delta},q_{\delta}]=\frac{1}{\delta^{2}}\int_{\mathcal{M}}\int_{\mathcal{M}}q_{\delta}({x})(p_{\delta}({x})-p_{\delta}({y}))R_{\delta}({x},{y})d\mu_{y}d\mu_{x}\\\
&+\int_{\partial\mathcal{M}}(\nabla_{\partial\mathcal{M}}\cdot\frac{\nabla_{\partial\mathcal{M}}Pp_{\delta}({y})}{\omega({y})})Pq_{\delta}({y})d\tau_{y}+\bar{P}(p_{\delta})\int_{\mathcal{M}}q_{\delta}({x})d\mu_{x}\\\
=&\frac{1}{2\delta^{2}}\int_{\mathcal{M}}\int_{\mathcal{M}}(q_{\delta}({x})-q_{\delta}({y}))(p_{\delta}({x})-p_{\delta}({y}))R_{\delta}({x},{y})d\mu_{y}d\mu_{x}\\\
&+\int_{\partial\mathcal{M}}\frac{\nabla_{\partial\mathcal{M}}Pp_{\delta}({y})\nabla_{\partial\mathcal{M}}Pq_{\delta}({y})}{-\omega({y})}d\tau_{y}\\\
\leq&C_{\delta}\int_{\mathcal{M}}|q_{\delta}({x})|(|p_{\delta}({x})|+|p_{\delta}({y})|)d\mu_{y}d\mu_{x}+C_{\delta}\int_{\partial\mathcal{M}}(\int_{\mathcal{M}}|p_{\delta}({x})|d\mu_{{x}}\int_{\mathcal{M}}|q_{\delta}({x})|d\mu_{{x}})d\mu_{{y}}\\\
\leq&C_{\delta}\left\|p_{\delta}\right\|_{L^{2}(\mathcal{M})}\left\|q_{\delta}\right\|_{L^{2}(\mathcal{M})},\end{split}$
(3.7)

what left in this lemma is the estimate (3.6). In fact, we replace
$q_{\delta}$ by $p_{\delta}$ in (3.7) to find

$\begin{split}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!B_{\delta}[p_{\delta},p_{\delta}]=\frac{1}{2\delta^{2}}\int_{\mathcal{M}}\int_{\mathcal{M}}(p_{\delta}({x})-p_{\delta}({y}))^{2}R_{\delta}({x},{y})d\mu_{y}d\mu_{x}+\int_{\partial\mathcal{M}}\frac{\nabla_{\partial\mathcal{M}}Pp_{\delta}({y})\cdot\nabla_{\partial\mathcal{M}}Pp_{\delta}({y})}{-\omega({y})}d\tau_{y}.\end{split}$
(3.8)

Recall the definition of $\omega$ in (2.8), we have
$C_{1}\delta\leq-\omega({y})\leq C_{2}\delta$ for some constant
$C_{1},C_{2}>0$ and all ${y}\in\partial\mathcal{M}$. The second term of (3.8)
is apparently non-negative. For the first term of (3.8), we apply Lemma 2 and
Lemma 3 of [38] to obtain

$\frac{1}{2\delta^{2}}\int_{\mathcal{M}}\int_{\mathcal{M}}(p_{\delta}({x})-p_{\delta}({y}))^{2}R_{\delta}({x},{y})d\mu_{y}d\mu_{x}\geq
C\left\|\nabla\tilde{p}_{\delta}\right\|_{L^{2}(\mathcal{M})}^{2},$ (3.9)

where $\tilde{p}_{\delta}$ is the weighted average of $p_{\delta}$ with
respect to $R_{\delta}$ on $\mathcal{M}$:

$\tilde{p}_{\delta}({x})=\frac{1}{\omega_{1}({x})}\int_{\mathcal{M}}p_{\delta}({y})R_{\delta}({x},{y})d\mu_{y},\
\ \omega_{1}({x})=\int_{\mathcal{M}}R_{\delta}({x},{y})d\mu_{y},\qquad\forall\
{x}\in\mathcal{M}.$ (3.10)

Besides, we apply Poincare inequality on $\bar{q}_{\delta}$:

$\left\|\nabla\tilde{p}_{\delta}\right\|^{2}_{L^{2}(\mathcal{M})}+(\int_{\mathcal{M}}\tilde{p}_{\delta}({x})d\mu_{{x}})^{2}\geq
C\left\|\tilde{p}_{\delta}\right\|^{2}_{L^{2}(\mathcal{M})}.$ (3.11)

For terms in (3.11), we apply Cauchy’s inequality to obtain the following two
controls:

$\begin{split}&\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!(\int_{\mathcal{M}}\tilde{p}_{\delta}({x})d\mu_{{x}})^{2}=(\int_{\mathcal{M}}\tilde{p}_{\delta}({x})-p_{\delta}({x})d\mu_{{x}})^{2}=\big{(}\int_{\mathcal{M}}\int_{\mathcal{M}}\frac{1}{\omega_{1}({x})}(p_{\delta}({y})-p_{\delta}({x}))R_{\delta}({x},{y})d\mu_{y}d\mu_{{x}}\big{)}^{2}\\\
&\leq
C\int_{\mathcal{M}}\int_{\mathcal{M}}(p_{\delta}({y})-p_{\delta}({x}))^{2}R_{\delta}({x},{y})d\mu_{y}d{\mu}_{{x}}\int_{\mathcal{M}}\int_{\mathcal{M}}\frac{1}{\omega^{2}_{1}({x})}R_{\delta}({x},{y})d\mu_{y}d\mu_{{x}}\\\
&\leq
C\int_{\mathcal{M}}\int_{\mathcal{M}}(p_{\delta}({y})-p_{\delta}({x}))^{2}R_{\delta}({x},{y})d\mu_{y}d{\mu}_{{x}}\leq
C\delta^{2}B_{\delta}[p_{\delta},p_{\delta}],\end{split}$ (3.12)

$\begin{split}&\big{|}\left\|\tilde{p}_{\delta}\right\|^{2}_{L^{2}(\mathcal{M})}-\left\|{p}_{\delta}\right\|^{2}_{L^{2}(\mathcal{M})}\big{|}\leq\left\|{p}_{\delta}-\tilde{p}_{\delta}\right\|^{2}_{L^{2}(\mathcal{M})}=\int_{\mathcal{M}}(\tilde{p}_{\delta}({x})-p_{\delta}({x}))^{2}d\mu_{{x}}\\\
&=\big{(}\int_{\mathcal{M}}(\int_{\mathcal{M}}\frac{1}{\omega_{1}({x})}(p_{\delta}({y})-p_{\delta}({x}))R_{\delta}({x},{y})d\mu_{y}\big{)}^{2}d\mu_{{x}}\\\
&\leq
C\int_{\mathcal{M}}(\int_{\mathcal{M}}(p_{\delta}({y})-p_{\delta}({x}))^{2}R_{\delta}({x},{y})d\mu_{y})(\int_{\mathcal{M}}\frac{1}{\omega^{2}_{1}({x})}R_{\delta}({x},{y})d\mu_{y})d\mu_{{x}}\\\
&\leq
C\int_{\mathcal{M}}\int_{\mathcal{M}}(p_{\delta}({y})-p_{\delta}({x}))^{2}R_{\delta}({x},{y})d\mu_{y}d{\mu}_{{x}}\leq
C\delta^{2}B_{\delta}[p_{\delta},p_{\delta}].\end{split}$ (3.13)

We take the above estimates on
$\int_{\mathcal{M}}\tilde{p}_{\delta}({x})d\mu_{{x}}$ and
$\left\|\tilde{p}_{\delta}\right\|^{2}_{L^{2}(\mathcal{M})}$ into (3.11) to
find

$\delta^{2}B_{\delta}[p_{\delta},p_{\delta}]+\left\|\nabla\tilde{p}_{\delta}\right\|^{2}_{L^{2}(\mathcal{M})}\geq
C\left\|{p}_{\delta}\right\|^{2}_{L^{2}(\mathcal{M})}$ (3.14)

We then derive (3.6) by utilizing the estimate (3.9).

Lemma 3.1 implies that $B_{\delta}[\cdot,\cdot]$ satisfy all the Lax-Milgram
conditions. Under such condition, we are able to conclude that there is a
unique solution $u_{\delta}\in L_{0}^{2}(\mathcal{M})$ to the weak formulation
(3.4), with the estimate

$\left\|u_{\delta}\right\|_{L^{2}(\mathcal{M})}\leq
C\left\|f_{\delta}\right\|_{L^{2}(\mathcal{M})};$ (3.15)

and equivalently, $u_{\delta}\in L_{0}^{2}(\mathcal{M})$ uniquely solves
(2.15) as well.

Now let us control the $L^{2}$ norm of $u_{\delta}$. Since we have

$\begin{split}&\\!\\!\\!\\!\\!\\!\\!\\!\\!\left\|f_{\delta}\right\|_{L^{2}(\mathcal{M})}=\big{(}\int_{\mathcal{M}}(\int_{\mathcal{M}}f({y})\bar{R}_{\delta}({x},{y})d\mu_{y}+\int_{\partial\mathcal{M}}({x}-{y})\cdot\bm{n}({y})f({y})\bar{R}_{\delta}({x},{y})d\tau_{y}-\bar{f})^{2}d\mu_{{x}}\big{)}^{\frac{1}{2}}\\\
=&\big{(}\int_{\mathcal{M}}(\int_{\mathcal{M}}f({y})\bar{R}_{\delta}({x},{y})d\mu_{y}+\int_{\partial\mathcal{M}}({x}-{y})\cdot\bm{n}({y})f({y})\bar{R}_{\delta}({x},{y})d\tau_{y})^{2}d\mu_{{x}}\big{)}^{\frac{1}{2}}\\\
\leq&C\big{(}\int_{\mathcal{M}}(\int_{\mathcal{M}}f^{2}({y})\bar{R}_{\delta}({x},{y})d\mu_{y})(\int_{\mathcal{M}}\bar{R}_{\delta}({x},{y})d\mu_{y})d\mu_{{x}}\big{)}^{\frac{1}{2}}\\\
&+C\big{(}\int_{\mathcal{M}}(\int_{\partial\mathcal{M}}\delta^{2}f^{2}({y})\bar{R}_{\delta}({x},{y})d\tau_{y})(\int_{\partial\mathcal{M}}\bar{R}_{\delta}({x},{y})d\tau_{y})d\mu_{{x}}\big{)}^{\frac{1}{2}}\\\
\leq&C\big{(}(\int_{\mathcal{M}}\int_{\mathcal{M}}f^{2}({y})\bar{R}_{\delta}({x},{y})d\mu_{y}d\mu_{{x}})^{\frac{1}{2}}+\delta(\int_{\mathcal{M}}\int_{\partial\mathcal{M}}f^{2}({y})\bar{R}_{\delta}({x},{y})d\tau_{y}d\mu_{{x}})^{\frac{1}{2}}\big{)}\\\
\leq&C(\left\|f\right\|_{L^{2}(\mathcal{M})}+\delta^{\frac{1}{2}}\left\|f\right\|_{L^{2}(\partial\mathcal{M})})\leq
C\left\|f\right\|_{H^{1}(\mathcal{M})},\end{split}$ (3.16)

this implies $\left\|u_{\delta}\right\|_{L^{2}(\mathcal{M})}\leq
C\left\|f\right\|_{H^{1}(\mathcal{M})}$.

What left in the first part of Theorem 2.1 is to present that $u_{\delta}\in
H^{1}(\mathcal{M})$. To this end, we will need the following lemma.

######  Lemma 3.3.

Let $r_{\delta}\in H^{1}(\mathcal{M})\cap L_{0}^{2}(\mathcal{M})$, and
$p_{\delta}\in L_{0}^{2}(\mathcal{M})$ solves

$B_{\delta}[p_{\delta},q_{\delta}]=(r_{\delta},q_{\delta})\qquad\forall\
q_{\delta}\in L_{0}^{2}(\mathcal{M}),$ (3.17)

where the function $r_{\delta}$ satisfies

$\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!(r_{\delta},q_{\delta})\leq
C\left\|r_{1\delta}\right\|_{H^{m}(\mathcal{M})}(\left\|q_{\delta}\right\|_{H^{1}(\mathcal{M})}+\left\|\bar{q}_{\delta}\right\|_{H^{1}(\mathcal{M})}+\left\|\overset{=}{q}_{\delta}\right\|_{H^{1}(\mathcal{M})}),\
\forall\ q_{\delta}\in H^{1}(\mathcal{M})\cap L_{0}^{2}(\mathcal{M})$ (3.18)

for some $m\in\mathbb{N}^{+}$ and some function $r_{1\delta}\in
H^{m}(\mathcal{M})$. Here the notations

$\bar{q}_{\delta}({x})=\frac{1}{\omega_{2}({x})}\int_{\mathcal{M}}q_{\delta}({y})\bar{R}_{\delta}({x},{y})d\mu_{{x}},\qquad\omega_{2}({x})=\int_{\mathcal{M}}\bar{R}_{\delta}({x},{y})d\mu_{{x}},$
(3.19)

and

$\overset{=}{q}_{\delta}({x})=\frac{1}{\omega_{3}({x})}\int_{\mathcal{M}}q_{\delta}({y})\overset{=}{R}_{\delta}({x},{y})d\mu_{{x}},\qquad\omega_{3}({x})=\int_{\mathcal{M}}\overset{=}{R}_{\delta}({x},{y})d\mu_{{x}}.$
(3.20)

Then we have $p_{\delta}\in H^{1}(\mathcal{M})$, with the estimate

$\left\|p_{\delta}\right\|_{H^{1}(\mathcal{M})}\leq
C(\left\|r_{1\delta}\right\|_{H^{m}(\mathcal{M})}+\delta^{2}\left\|\nabla_{\mathcal{M}}r_{\delta}\right\|_{L^{2}(\mathcal{M})}).$
(3.21)

######  Proof 3.4.

As we mentioned, (3.17) is equivalent to the following equation:

$\mathcal{L}_{\delta}p_{\delta}({x})=r_{\delta}({x})\qquad\forall\
{x}\in\mathcal{M}.$ (3.22)

We take the derivative of (3.22) to obtain

$\nabla^{{x}}_{\mathcal{M}}\mathcal{L}_{\delta}p_{\delta}({x})=\nabla^{{x}}_{\mathcal{M}}r_{\delta}({x})\qquad\forall\
{x}\in\mathcal{M},$ (3.23)

recall the definition of $\mathcal{L}_{\delta}$ in (2.14), we write down the
complete form of (3.23) as

$\begin{split}&\frac{1}{\delta^{2}}\nabla^{{x}}_{\mathcal{M}}\int_{\mathcal{M}}p_{\delta}({x})R_{\delta}({x},{y})d\mu_{y}-\frac{1}{\delta^{2}}\nabla^{{x}}_{\mathcal{M}}\int_{\mathcal{M}}p_{\delta}({y})R_{\delta}({x},{y})d\mu_{y}\\\
&+\nabla^{{x}}_{\mathcal{M}}\int_{\partial\mathcal{M}}(\nabla^{{y}}_{\partial\mathcal{M}}\cdot\frac{\nabla^{{y}}_{\partial\mathcal{M}}Pp_{\delta}({y})}{\omega({y})})\bm{n}({y})\cdot({x}-{y})\bar{R}_{\delta}({x},{y})d\tau_{y}=\nabla^{{x}}_{\mathcal{M}}r_{\delta}({x})\qquad\forall\
{x}\in\mathcal{M}.\end{split}$ (3.24)

To further simplify (3.24), we introduce again the auxiliary functions
$\tilde{p}_{\delta}$ and $\omega_{1}$ defined in (3.10), and do integration by
parts to the boundary term that involves $Pp_{\delta}$, to rewrite (3.24) as

$\begin{split}&\nabla^{{x}}_{\mathcal{M}}(p_{\delta}({x})\omega_{1}({x}))-\nabla^{{x}}_{\mathcal{M}}(\tilde{p}_{\delta}({x})\omega_{1}({x}))\\\
&-\delta^{2}\nabla^{{x}}_{\mathcal{M}}\int_{\partial\mathcal{M}}(\frac{\nabla^{{y}}_{\partial\mathcal{M}}Pp_{\delta}({y})}{\omega({y})})\cdot\nabla^{{y}}_{\partial\mathcal{M}}(\bm{n}({y})\cdot({x}-{y})\bar{R}_{\delta}({x},{y}))d\tau_{y}=\delta^{2}\nabla^{{x}}_{\mathcal{M}}r_{\delta}({x}),\
\forall\ {x}\in\mathcal{M},\end{split}$ (3.25)

we then obtain the following equality for $\nabla_{\mathcal{M}}p_{\delta}$:

$\begin{split}\nabla^{{x}}_{\mathcal{M}}p_{\delta}({x})=&\nabla^{{x}}_{\mathcal{M}}\tilde{p}_{\delta}({x})-\frac{\nabla^{{x}}_{\mathcal{M}}\omega_{1}({x})}{\omega_{1}({x})}(p_{\delta}({x})-\tilde{p}_{\delta}({x}))+\delta^{2}\frac{1}{\omega_{1}({x})}\nabla^{{x}}_{\mathcal{M}}r_{\delta}({x})\\\
&+\delta^{2}\int_{\partial\mathcal{M}}\frac{\nabla^{{y}}_{\partial\mathcal{M}}Pp_{\delta}({y})}{\omega({y})}\cdot\frac{\nabla^{{y}}_{\partial\mathcal{M}}\nabla^{{x}}_{\mathcal{M}}(\bm{n}({y})\cdot({x}-{y})\bar{R}_{\delta}({x},{y}))}{\omega_{1}({x})}d\tau_{y},\
\forall\ {x}\in\mathcal{M}.\end{split}$ (3.26)

The $L^{2}(\mathcal{M})$ norm of the first $3$ terms of (3.26) are comparably
easy to control. By the definition of $\omega_{1}$ in (3.10), we have
$C_{3}\leq\omega_{1}({x})\leq C_{4}$ for some constant $C_{3},C_{4}>0$ and all
${x}\in\mathcal{M}$, with

$\begin{split}|\nabla^{{x}}_{\mathcal{M}}\omega_{1}({x})|=&|\int_{\mathcal{M}}\nabla^{{x}}_{\mathcal{M}}R_{\delta}({x},{y})d\mu_{{y}}|=|\frac{1}{2\delta^{2}}\int_{\mathcal{M}}\underline{R}_{\delta}({x},{y})\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})d\mu_{{y}}|\\\
&\leq
C\delta^{-1}\int_{\mathcal{M}}\underline{R}_{\delta}({x},{y})d\mu_{{y}}\leq
C\delta^{-1},\end{split}$ (3.27)

where $\underline{R}_{\delta}:\mathcal{M}\times\mathcal{M}\to\mathbb{R}$ is
the kernel function defined in (2.3). We then combine the control of
$\left\|\nabla_{\mathcal{M}}\tilde{p}_{\delta}\right\|_{L^{2}(\mathcal{M})}$
and $\left\|p_{\delta}-\tilde{p}_{\delta}\right\|_{L^{2}(\mathcal{M})}$ in
(3.9) (3.13) to conclude

$\begin{split}&\left\|\nabla^{{x}}_{\mathcal{M}}\tilde{p}_{\delta}({x})-\frac{\nabla^{{x}}_{\mathcal{M}}\omega_{1}({x})}{\omega_{1}({x})}(p_{\delta}({x})-\tilde{p}_{\delta}({x}))+\delta^{2}\frac{1}{\omega_{1}({x})}\nabla^{{x}}_{\mathcal{M}}r_{\delta}({x})\right\|^{2}_{L_{{x}}^{2}(\mathcal{M})}\\\
&\leq
C(B_{\delta}[p_{\delta},p_{\delta}]+\delta^{-2}\delta^{2}B_{\delta}[p_{\delta},p_{\delta}]+\delta^{4}\left\|\nabla_{\mathcal{M}}r_{\delta}\right\|^{2}_{L^{2}(\mathcal{M})})\\\
&\leq
CB_{\delta}[p_{\delta},p_{\delta}]+C\delta^{4}\left\|\nabla_{\mathcal{M}}r_{\delta}\right\|^{2}_{L^{2}(\mathcal{M})}.\end{split}$
(3.28)

The last term of (3.26) is more involved, we will use the notations in
Appendix A. First, we notice
$\nabla^{{y}}_{\partial\mathcal{M}}Pp_{\delta}({y})\in\mathcal{T}_{\partial\mathcal{M}}^{{y}}$.
On the other hand, for any unit vector
$\bm{\tau}\in\mathcal{T}_{\partial\mathcal{M}}^{{y}}$, we can calculate

$\begin{split}&\bm{\tau}\cdot\nabla^{{y}}_{\partial\mathcal{M}}\big{(}\nabla^{{x}}_{\mathcal{M}}(\bm{n}({y})\cdot({x}-{y})\bar{R}_{\delta}({x},{y}))\big{)}\\\
=&\tau^{i}\partial^{i}_{{y}}\big{(}n({y})\bar{R}_{\delta}({x},{y})-({x}^{l}-{y}^{l})n^{l}({y})\frac{1}{2\delta^{2}}(\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})\big{)}\\\
=&\tau^{i}\partial^{i}_{{y}}\bm{n}({y})\bar{R}_{\delta}({x},{y})+\bm{n}({y})\frac{1}{2\delta^{2}}\tau^{i}(\mathcal{T}_{\partial\mathcal{M}}^{{y}}({x}-{y}))^{i}{R}_{\delta}({x},{y})\\\
&-({x}^{l}-{y}^{l})\tau^{i}\partial^{i}_{{y}}n^{l}({y})\frac{1}{2\delta^{2}}(\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})+({x}^{l}-{y}^{l})n^{l}({y})\frac{1}{2\delta^{2}}(\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}){R}_{\delta}({x},{y})\\\
&-({x}^{l}-{y}^{l})n^{l}({y})\frac{1}{2\delta^{2}}(\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}))\frac{1}{2\delta^{2}}\tau^{i}(\mathcal{T}_{\partial\mathcal{M}}^{{y}}({x}-{y}))^{i}\underline{R}_{\delta}({x},{y})\\\
=&\textbf{p}_{1\delta}({x},{y},\bm{\tau})\bar{R}_{\delta}({x},{y})+\textbf{p}_{2\delta}({x},{y},\bm{\tau}){R}_{\delta}({x},{y})+\textbf{p}_{3\delta}({x},{y},\bm{\tau})\underline{R}_{\delta}({x},{y}),\end{split}$
(3.29)

where $|\textbf{p}_{k\delta}({x},{y},\bm{\tau})|\leq C\delta^{-1},\ k=1,2,3$
by a simple observation on each term that composes (3.29). For brevity’s sake,
we denote
$R_{1\delta}=\bar{R}_{\delta},R_{2\delta}=R_{\delta},R_{3\delta}=\underline{R}_{\delta}$.
Next, we take
$\bm{\tau}=\nabla^{{y}}_{\partial\mathcal{M}}Pp_{\delta}({y})/|\nabla^{{y}}_{\partial\mathcal{M}}Pp_{\delta}({y})|$,
and apply (3.29) to the following calculation:

$\begin{split}\int_{\mathcal{M}}&(\delta^{2}\int_{\partial\mathcal{M}}\frac{\nabla^{{y}}_{\partial\mathcal{M}}Pp_{\delta}({y})}{\omega({y})}\cdot\frac{\nabla^{{y}}_{\partial\mathcal{M}}\nabla^{{x}}_{\mathcal{M}}(\bm{n}({y})\cdot({x}-{y})\bar{R}_{\delta}({x},{y}))}{\omega_{1}({x})}d\tau_{y})^{2}d\mu_{{x}}\\\
\leq&\sum\limits_{k=1}^{3}\int_{\mathcal{M}}(\delta^{2}\int_{\partial\mathcal{M}}\frac{|\nabla^{{y}}_{\partial\mathcal{M}}Pp_{\delta}({y})|}{-\omega({y})}\frac{|\textbf{p}_{k\delta}({x},{y},\bm{\tau})|{R}_{k\delta}({x},{y})}{\omega_{1}({x})}d\tau_{y})^{2}d\mu_{{x}}\\\
\leq&C\sum\limits_{k=1}^{3}\int_{\mathcal{M}}(\int_{\partial\mathcal{M}}|\nabla^{{y}}_{\partial\mathcal{M}}Pp_{\delta}({y})|{R}_{k\delta}({x},{y})d\tau_{y})^{2}d\mu_{{x}}\\\
\leq&C\sum\limits_{k=1}^{3}\int_{\mathcal{M}}(\int_{\partial\mathcal{M}}|\nabla^{{y}}_{\partial\mathcal{M}}Pp_{\delta}({y})|^{2}{R}_{k\delta}({x},{y})d\tau_{y})(\int_{\partial\mathcal{M}}{R}_{k\delta}({x},{y})d\tau_{y})d\mu_{{x}}\\\
\leq&C\delta^{-1}\sum\limits_{k=1}^{3}\int_{\mathcal{M}}\int_{\partial\mathcal{M}}|\nabla^{{y}}_{\partial\mathcal{M}}Pp_{\delta}({y})|^{2}{R}_{k\delta}({x},{y})d\tau_{y}d\mu_{{x}}\\\
\leq&C\delta^{-1}\int_{\partial\mathcal{M}}|\nabla^{{y}}_{\partial\mathcal{M}}Pp_{\delta}({y})|^{2}d\tau_{y}\leq
C\int_{\partial\mathcal{M}}\frac{\nabla_{\partial\mathcal{M}}Pp_{\delta}({y})\cdot\nabla_{\partial\mathcal{M}}Pp_{\delta}({y})}{-\omega({y})}d\tau_{y}\leq
CB_{\delta}[p_{\delta},p_{\delta}].\end{split}$ (3.30)

We then employ (3.28) and (3.30) to (3.26) to discover

$\left\|\nabla_{\mathcal{M}}p_{\delta}\right\|^{2}_{L^{2}(\mathcal{M})}\leq
C(B_{\delta}[p_{\delta},p_{\delta}]+\delta^{4}\left\|\nabla_{\mathcal{M}}r_{\delta}\right\|^{2}_{L^{2}(\mathcal{M})});$
(3.31)

together with the estimate (3.6) on
$\left\|p_{\delta}\right\|^{2}_{L^{2}(\mathcal{M})}$, we conclude that

$\left\|p_{\delta}\right\|^{2}_{H^{1}(\mathcal{M})}\leq
C_{1}(B_{\delta}[p_{\delta},p_{\delta}]+\delta^{4}\left\|\nabla_{\mathcal{M}}r_{\delta}\right\|^{2}_{L^{2}(\mathcal{M})})$
(3.32)

for some constant $C_{1}$ depend only on the shape of $\mathcal{M}$.

Finally, let us return to the assumption of our lemma to study
$B_{\delta}[p_{\delta},p_{\delta}]$. In fact, we apply (3.17) and (3.18) to
find

$B_{\delta}[p_{\delta},p_{\delta}]=(r_{\delta},p_{\delta})\leq
C\left\|r_{1\delta}\right\|_{H^{m}(\mathcal{M})}(\left\|p_{\delta}\right\|_{H^{1}(\mathcal{M})}+\left\|\bar{p}_{\delta}\right\|_{H^{1}(\mathcal{M})}+\left\|\overset{=}{p}_{\delta}\right\|_{H^{1}(\mathcal{M})}),$
(3.33)

here $\bar{p}_{\delta}$ and $\overset{=}{p}_{\delta}$ are defined in (3.19)
and (3.20) respectively after replacing $q_{\delta}$ by $p_{\delta}$. By
substituting $R_{\delta}$ into $\bar{R}_{\delta}$ in (3.9) and (3.13), we have

$\begin{split}\left\|\bar{p}_{\delta}-p_{\delta}\right\|^{2}_{L^{2}(\mathcal{M})}\leq
C\int_{\mathcal{M}}\int_{\mathcal{M}}(p_{\delta}({y})-p_{\delta}({x}))^{2}\bar{R}_{\delta}({x},{y})d\mu_{y}d{\mu}_{{x}};\end{split}$
(3.34)
$\frac{1}{2\delta^{2}}\int_{\mathcal{M}}\int_{\mathcal{M}}(p_{\delta}({x})-p_{\delta}({y}))^{2}\bar{R}_{\delta}({x},{y})d\mu_{y}d\mu_{x}\geq
C\left\|\nabla\bar{p}_{\delta}\right\|_{L^{2}(\mathcal{M})}^{2}.$ (3.35)

Besides, the equation (3.12) of [47] gives that

$\int_{\mathcal{M}}\int_{\mathcal{M}}(p_{\delta}({x})-p_{\delta}({y}))^{2}\bar{R}_{\delta}({x},{y})d\mu_{y}d\mu_{x}\leq
C\int_{\mathcal{M}}\int_{\mathcal{M}}(p_{\delta}({x})-p_{\delta}({y}))^{2}{R}_{\delta}({x},{y})d\mu_{y}d\mu_{x}$
(3.36)

We combine the estimates (3.34) (3.35) (3.36) and (3.6) to obtain

$\left\|\bar{p}_{\delta}\right\|^{2}_{H^{1}(\mathcal{M})}\leq
C_{2}\frac{1}{\delta^{2}}\int_{\mathcal{M}}\int_{\mathcal{M}}(p_{\delta}({y})-p_{\delta}({x}))^{2}R_{\delta}({x},{y})d\mu_{y}d{\mu}_{{x}}\leq
C_{2}B_{\delta}[p_{\delta},p_{\delta}]$ (3.37)

for some constant $C_{2}>0$. By the same argument as $\bar{p}_{\delta}$, we
have

$\left\|\overset{=}{p}_{\delta}\right\|^{2}_{H^{1}(\mathcal{M})}\leq
C_{3}B_{\delta}[p_{\delta},p_{\delta}]$ (3.38)

for some constant $C_{3}>0$ as well. We apply such estimates of
$\bar{p}_{\delta}$ and $\overset{=}{p}_{\delta}$, together with the estimates
of $p_{\delta}$ in (3.32) into (3.33) to discover

$\begin{split}&\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!B_{\delta}[p_{\delta},p_{\delta}]\leq
C\left\|r_{1\delta}\right\|_{H^{m}(\mathcal{M})}(\left\|p_{\delta}\right\|_{H^{1}(\mathcal{M})}+\left\|\bar{p}_{\delta}\right\|_{H^{1}(\mathcal{M})}+\left\|\overset{=}{p}_{\delta}\right\|_{H^{1}(\mathcal{M})})\\\
\leq&(C^{2}C_{1}\left\|r_{1\delta}\right\|^{2}_{H^{m}(\mathcal{M})}+\frac{1}{4C_{1}}\left\|p_{\delta}\right\|^{2}_{H^{1}(\mathcal{M})})+(C^{2}C_{2}\left\|r_{1\delta}\right\|^{2}_{H^{m}(\mathcal{M})}+\frac{1}{4C_{2}}\left\|\bar{p}_{\delta}\right\|^{2}_{H^{1}(\mathcal{M})})\\\
&+(C^{2}C_{3}\left\|r_{1\delta}\right\|^{2}_{H^{m}(\mathcal{M})}+\frac{1}{4C_{3}}\left\|\overset{=}{p}_{\delta}\right\|^{2}_{H^{1}(\mathcal{M})})\\\
\leq&C^{2}(C_{1}+C_{2}+C_{3})\left\|r_{1\delta}\right\|^{2}_{H^{m}(\mathcal{M})}+\frac{3}{4}B_{\delta}[p_{\delta},p_{\delta}]+\frac{1}{4}\delta^{4}\left\|\nabla_{\mathcal{M}}r_{\delta}\right\|^{2}_{L^{2}(\mathcal{M})}.\end{split}$
(3.39)

We then subtract both sides of (3.39) by
$\frac{3}{4}B_{\delta}[p_{\delta},p_{\delta}]$ to obtain

$B_{\delta}[p_{\delta},p_{\delta}]\leq
C(\left\|r_{1\delta}\right\|^{2}_{H^{m}(\mathcal{M})}+\delta^{4}\left\|\nabla_{\mathcal{M}}r_{\delta}\right\|^{2}_{L^{2}(\mathcal{M})}).$
(3.40)

Finally, we apply (3.32) into (3.40) to conclude

$\left\|p_{\delta}\right\|_{H^{1}(\mathcal{M})}\leq
C(\left\|r_{1\delta}\right\|_{H^{m}(\mathcal{M})}+\delta^{2}\left\|\nabla_{\mathcal{M}}r_{\delta}\right\|_{L^{2}(\mathcal{M})}).$
(3.41)

We then complete the proof of lemma.

We directly apply Lemma 3.3 on the bilinear form (3.4) to obtain

$\left\|u_{\delta}\right\|_{H^{1}(\mathcal{M})}\leq
C(\left\|f_{\delta}\right\|_{L^{2}(\mathcal{M})}+\delta^{2}\left\|\nabla_{\mathcal{M}}f_{\delta}\right\|_{L^{2}(\mathcal{M})}).$
(3.42)

The $L^{2}$ norm of $f_{\delta}$ has been controlled in (3.16), for
$\nabla_{\mathcal{M}}f_{\delta}$, we calculate

$\begin{split}|\nabla_{\mathcal{M}}&f_{\delta}({x})|=\big{|}\nabla_{\mathcal{M}}^{{x}}\int_{\mathcal{M}}f({y})\bar{R}_{\delta}({x},{y})d\mu_{y}+\nabla_{\mathcal{M}}^{{x}}\int_{\partial\mathcal{M}}({x}-{y})\cdot\bm{n}({y})f({y})\bar{R}_{\delta}({x},{y})d\tau_{y}\big{|}\\\
=&\big{|}\frac{1}{2\delta^{2}}\int_{\mathcal{M}}f({y})\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}){R}_{\delta}({x},{y})d\mu_{y}+\int_{\partial\mathcal{M}}\bm{n}({y})f({y})\bar{R}_{\delta}({x},{y})d\tau_{y}\\\
&+\frac{1}{2\delta^{2}}\int_{\partial\mathcal{M}}({x}-{y})\cdot\bm{n}({y})f({y})\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}){R}_{\delta}({x},{y})d\tau_{y}\big{|}\\\
\leq&C(\delta^{-1}\int_{\mathcal{M}}|f({y})|{R}_{\delta}({x},{y})d\mu_{y}+\int_{\partial\mathcal{M}}|f({y})|\bar{R}_{\delta}({x},{y})d\tau_{y}+\int_{\partial\mathcal{M}}|f({y})|{R}_{\delta}({x},{y})d\tau_{y}),\end{split}$
(3.43)

By the same arguments as (3.16), we can calculate

$\begin{split}\left\|\nabla_{\mathcal{M}}f_{\delta}\right\|_{L^{2}(\mathcal{M})}\leq
C(\delta^{-1}\left\|f\right\|_{L^{2}(\mathcal{M})}+\delta^{-1/2}\left\|f\right\|_{L^{2}(\partial\mathcal{M})})\leq
C\delta^{-1}\left\|f\right\|_{H^{1}(\mathcal{M})}.\end{split}$ (3.44)

We then apply the estimates (3.16) (3.44) into (3.42) to conclude

$\left\|u_{\delta}\right\|_{H^{1}(\mathcal{M})}\leq
C\left\|f\right\|_{H^{1}(\mathcal{M})}.$ (3.45)

##  4 Vanishing Nonlocality

In this section, our purpose is to prove the second part of Theorem 2.1. In
fact, we let

$e_{\delta}=u-u_{\delta},$ (4.1)

then $e_{\delta}\in L_{0}^{2}(\mathcal{M})$, and we subtract the equation
(2.14) by (2.13) to discover

$\begin{split}\mathcal{L}_{\delta}e_{\delta}({x})=\bar{Q}(f)-\bar{P}(u)+r_{1}({x})-r_{2}({x}),\qquad\forall\
{x}\in\mathcal{M},\end{split}$ (4.2)

where $r_{1},r_{2}$ are the functions defined in (1.3) (2.9),
$\bar{P}(u),\bar{Q}(f)$ are the functionals defined in (2.11) (2.12)
respectively. This equation implies

$B[e_{\delta},q_{\delta}]=(\mathcal{L}_{\delta}e_{\delta},q_{\delta})=(r_{1}-r_{2},q_{\delta}),\qquad\forall\
q_{\delta}\in L_{0}^{2}(\mathcal{M}).$ (4.3)

According to Lemma 3.3, in condition of (4.3), if we can present

$\left\|\nabla_{\mathcal{M}}r_{1}\right\|\leq
C\delta^{\frac{1}{2}}\left\|u\right\|_{H^{4}(\mathcal{M})},\qquad\left\|\nabla_{\mathcal{M}}r_{2}\right\|\leq
C\delta^{\frac{1}{2}}\left\|u\right\|_{H^{4}(\mathcal{M})},$ (4.4)

and

$\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!(r_{k},q_{\delta})\leq
C\delta^{2}\left\|u\right\|_{H^{4}(\mathcal{M})}(\left\|q_{\delta}\right\|_{H^{1}(\mathcal{M})}+\left\|\bar{q}_{\delta}\right\|_{H^{1}(\mathcal{M})}+\left\|\overset{=}{q}_{\delta}\right\|_{H^{1}(\mathcal{M})}),\
\forall\ q_{\delta}\in H^{1}(\mathcal{M})\cap L_{0}^{2}(\mathcal{M})$ (4.5)

for $k=1,2$, then we will have

$\left\|e_{\delta}\right\|_{H^{1}(\mathcal{M})}\leq
C\delta^{2}\left\|u\right\|_{H^{4}(\mathcal{M})}+C\delta^{2}\left\|\nabla_{\mathcal{M}}(r_{1}-r_{2})\right\|_{L^{2}(\mathcal{M})}\leq
C\delta^{2}\left\|u\right\|_{H^{4}(\mathcal{M})},$ (4.6)

which is exactly the second part of Theorem 2.1. Our task is then reduced to
the proof of (4.4) and (4.5).

The term $r_{1}$ is exactly the term $r_{in}$ in [48]. In [48], we have
devoted the whole content on the control of $r_{in}$, and it satisfy the
conditions (4.4) and (4.5) according to Theorem 3.1 of [48].

For the term $r_{2}$, we recall its definition (2.9) to write

$\begin{split}r_{2}({x})=\int_{\partial\mathcal{M}}b({y})\bm{n}({y})\cdot({x}-{y})\bar{R}_{\delta}({x},{y})d\tau_{y},\qquad{x}\in\mathcal{M},\end{split}$
(4.7)

where

$b({y})=\Delta_{\partial\mathcal{M}}u({y})-\nabla_{\partial\mathcal{M}}\cdot\frac{\nabla_{\partial\mathcal{M}}Pu({y})}{\omega({y})},\qquad{y}\in\partial\mathcal{M}.$
(4.8)

By a simple observation, the estimates (4.4) and (4.5) on $r_{2}$ are direct
corollaries of the following two lemmas:

######  Lemma 4.1.

Let the function $b:\partial\mathcal{M}\to\mathbb{R}$ be defined in (4.8),
then

$\left\|b\right\|_{L^{2}(\partial\mathcal{M})}\leq
C\delta\left\|u\right\|_{H^{4}(\mathcal{M})}.$ (4.9)

######  Lemma 4.2.

Let $r_{2}$ be defined in (2.9), then we have the following two estimates:

$\left\|\nabla_{\mathcal{M}}r_{2}\right\|_{L^{2}(\mathcal{M})}\leq
C\delta^{-1/2}\left\|b\right\|_{L^{2}(\partial\mathcal{M})};$ (4.10)

and

$\begin{split}\int_{\mathcal{M}}r_{2}({x})q_{\delta}({x})d\mu_{x}\leq
C\delta\left\|b\right\|_{L^{2}(\partial\mathcal{M})}\left\|q_{\delta}\right\|_{H^{1}(\mathcal{M})},\qquad\forall\
q_{\delta}\in H^{1}(\mathcal{M})\cap L_{0}^{2}(\mathcal{M}).\end{split}$
(4.11)

We have moved the proof of Lemma 4.1 into Appendix due to its complication.

######  Proof 4.3 (Proof of Lemma 4.2).

We can calculate

$\begin{split}&\left\|\nabla_{\mathcal{M}}r_{2}\right\|^{2}_{L^{2}(\mathcal{M})}=\int_{\mathcal{M}}|\nabla^{{x}}_{\mathcal{M}}\int_{\partial\mathcal{M}}b({y})\bm{n}({y})\cdot({x}-{y})\bar{R}_{\delta}({x},{y})d\tau_{y}|^{2}d\mu_{{x}}\\\
&=\int_{\mathcal{M}}|\int_{\partial\mathcal{M}}b({y})\bm{n}({y})\bar{R}_{\delta}({x},{y})-b({y})\bm{n}({y})\cdot({x}-{y})\frac{1}{2\delta^{2}}\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}){R}_{\delta}({x},{y})d\tau_{y}|^{2}d\mu_{{x}}\\\
&\leq
C\int_{\mathcal{M}}(\int_{\partial\mathcal{M}}|b({y})|(\bar{R}_{\delta}({x},{y})+{R}_{\delta}({x},{y}))d\tau_{y})^{2}d\mu_{{x}}\\\
&\leq
C\int_{\mathcal{M}}(\int_{\partial\mathcal{M}}b^{2}({y})(\bar{R}_{\delta}({x},{y})+{R}_{\delta}({x},{y}))d\tau_{y})(\int_{\partial\mathcal{M}}(\bar{R}_{\delta}({x},{y})+{R}_{\delta}({x},{y}))d\tau_{y})d\mu_{{x}}\\\
&\leq\frac{C}{\delta}\int_{\mathcal{M}}(\int_{\partial\mathcal{M}}b^{2}({y})\bar{R}_{\delta}({x},{y})+{R}_{\delta}({x},{y}))d\tau_{y}d\mu_{{x}}\leq\frac{C}{\delta}\int_{\partial\mathcal{M}}b^{2}({y})d\tau_{y}\leq\frac{C}{\delta}\left\|b\right\|^{2}_{L^{2}(\partial\mathcal{M})},\end{split}$
(4.12)

and

$\begin{split}&\int_{\mathcal{M}}r_{2}({x})q_{\delta}({x})d\mu_{x}=\int_{\mathcal{M}}\int_{\partial\mathcal{M}}q_{\delta}({x})b({y})\bm{n}({y})\cdot({x}-{y})\bar{R}_{\delta}({x},{y})d\tau_{y}d\mu_{{x}}\\\
=&-\int_{\mathcal{M}}\int_{\partial\mathcal{M}}q_{\delta}({x})b({y})\bm{n}({y})2\delta^{2}\nabla^{{x}}_{\mathcal{M}}\overset{=}{R}_{\delta}({x},{y})d\tau_{y}d\mu_{{x}}\\\
&+\int_{\mathcal{M}}\int_{\partial\mathcal{M}}q_{\delta}({x})b({y})\bm{n}({y})\cdot({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}))\bar{R}_{\delta}({x},{y})d\tau_{y}d\mu_{{x}}\\\
=&-\int_{\partial\mathcal{M}}\int_{\partial\mathcal{M}}q_{\delta}({x})b({y})\bm{n}({y})\cdot\bm{n}({x})2\delta^{2}\overset{=}{R}_{\delta}({x},{y})d\tau_{y}d\tau_{{x}}\\\
&+\int_{\mathcal{M}}\int_{\partial\mathcal{M}}(\nabla_{\mathcal{M}}^{{x}}q_{\delta}({x}))b({y})\bm{n}({y})2\delta^{2}\overset{=}{R}_{\delta}({x},{y})d\tau_{y}d\mu_{{x}}\\\
&+\int_{\mathcal{M}}\int_{\partial\mathcal{M}}q_{\delta}({x})b({y})\bm{n}({y})\cdot({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}))\bar{R}_{\delta}({x},{y})d\tau_{y}d\mu_{{x}}\\\
\leq&C\delta^{2}(\int_{\partial\mathcal{M}}\int_{\partial\mathcal{M}}|q_{\delta}({x})b({y})|\overset{=}{R}_{\delta}({x},{y})d\tau_{y}d\tau_{{x}}+\int_{\mathcal{M}}\int_{\partial\mathcal{M}}|\nabla_{\mathcal{M}}^{{x}}q_{\delta}({x})b({y})|\overset{=}{R}_{\delta}({x},{y})d\tau_{y}d\mu_{{x}}\\\
&+\int_{\mathcal{M}}\int_{\partial\mathcal{M}}|q_{\delta}({x})b({y})|\bar{R}_{\delta}({x},{y})d\tau_{y}d\mu_{{x}}).\end{split}$
(4.13)

Using Cauchy’s inequality, we have

$\begin{split}\int_{\partial\mathcal{M}}&\int_{\partial\mathcal{M}}|q_{\delta}({x})b({y})|\overset{=}{R}_{\delta}({x},{y})d\tau_{y}d\tau_{{x}}\\\
\leq&C(\int_{\partial\mathcal{M}}\int_{\partial\mathcal{M}}q_{\delta}^{2}({x})\overset{=}{R}_{\delta}({x},{y})d\tau_{y}d\tau_{{x}}\int_{\partial\mathcal{M}}\int_{\partial\mathcal{M}}b^{2}({y})\overset{=}{R}_{\delta}({x},{y})d\tau_{y}d\tau_{{x}})^{\frac{1}{2}}\\\
\leq&C(\frac{1}{\delta^{2}}\int_{\partial\mathcal{M}}q_{\delta}^{2}({x})d\tau_{{x}}\int_{\partial\mathcal{M}}b^{2}({y})d\tau_{y})^{\frac{1}{2}}\leq\frac{C}{\delta}\left\|q_{\delta}\right\|_{L^{2}(\partial\mathcal{M})}\left\|b\right\|_{L^{2}(\partial\mathcal{M})}\\\
\leq&\frac{C}{\delta}\left\|q_{\delta}\right\|_{H^{1}(\mathcal{M})}\left\|b\right\|_{L^{2}(\partial\mathcal{M})}.\end{split}$
(4.14)

By a similar argument as (4.14), the last two terms of (4.13) have the same
bound as (4.14). This indicates

$\int_{\mathcal{M}}r_{2}({x})q_{\delta}({x})d\mu_{x}\leq
C\delta\left\|q_{\delta}\right\|_{H^{1}(\mathcal{M})}\left\|b\right\|_{L^{2}(\partial\mathcal{M})}.$
(4.15)

We then complete the proof of Lemma 4.2.

Consequently, the error functions $r_{1}$ and $r_{2}$ both satisfy the
estimates (4.4) (4.5). According to Lemma 3.3, the second part of Theorem 2.1
holds.

##  5 Conclusion and Future Work

In this paper, we propose a nonlocal approximation of manifold Poisson model
with Neumann boundary. The well-posedness and the second order localization
rate is analyzed. The application of our nonlocal model is that, such model
can be easily discretized and be solved numerically by meshless scheme like
point integral method(PIM, see [28] [39]), while this solution can be regarded
as the numerical solution of the local Poisson model as well. Compared to the
traditional manifold finite element method(MFEM), one advantage of PIM is that
only local mesh is required so that we do not need a global mesh like MFEM. In
our subsequent work, we will implement several numerical examples to
illustrate the rate of convergence on the other side. Besides, we will try to
seek for accurate nonlocal approximations of some other PDEs, for instance,
the elliptic transmission problems, the Stokes problems, and the Poisson
problems with polygonal domains.

##  Appendix A Proof of Lemma 4.1

Our goal in this section is to prove

$\left\|b\right\|_{L^{2}(\partial\mathcal{M})}\leq
C\delta\left\|u\right\|_{H^{4}(\mathcal{M})}.$ (A.1)

First of all, let us do the following decomposition:

$\begin{split}b({y})=&\Delta_{\partial\mathcal{M}}u({y})-\nabla_{\partial\mathcal{M}}\cdot\frac{\nabla_{\partial\mathcal{M}}Pu({y})}{\omega({y})}=\nabla_{\partial\mathcal{M}}\cdot(\nabla_{\partial\mathcal{M}}u({y})-\frac{\nabla_{\partial\mathcal{M}}Pu({y})}{\omega({y})})\\\
=&\nabla_{\partial\mathcal{M}}\cdot\bigg{(}\frac{1}{\omega({y})}\Big{(}\nabla_{\partial\mathcal{M}}(u({y})\omega({y})-Pu({y}))-\
u({y})\nabla_{\partial\mathcal{M}}\omega({y})\Big{)}\bigg{)}\\\
=&\nabla_{\partial\mathcal{M}}\cdot\bigg{(}\frac{1}{\omega({y})}\Big{(}\nabla_{\partial\mathcal{M}}s({y})-\
u({y})\nabla_{\partial\mathcal{M}}\omega({y})\Big{)}\bigg{)}\\\
=&\frac{1}{\omega({y})}\Delta_{\partial\mathcal{M}}s({y})-\frac{\nabla_{\partial\mathcal{M}}\omega({y})}{\omega^{2}({y})}\cdot\nabla_{\partial\mathcal{M}}s({y})-\frac{\Delta_{\partial\mathcal{M}}\omega({y})}{\omega({y})}u({y})\\\
&-\frac{\nabla_{\partial\mathcal{M}}\omega({y})\cdot\nabla_{\partial\mathcal{M}}\omega({y})}{\omega^{2}({y})}u({y})+\frac{\nabla_{\partial\mathcal{M}}\omega({y})}{\omega({y})}\cdot\nabla_{\partial\mathcal{M}}u({y});\end{split}$
(A.2)

where the function

$s({y})=u({y})\omega({y})-Pu({y})=-\int_{\mathcal{M}}(u({x})-u({y}))(({x}-{y})\cdot\bm{n}({y}))\bar{R}_{\delta}({x},{y})d\mu_{x},$
(A.3)

and $\omega$ is defined in (2.8). Apparently, there are constants
$C_{1},C_{2}>0$ such that $C_{1}\delta\leq-\omega({y})\leq C_{2}\delta$ for
all ${y}\in\partial\mathcal{M}$, therefore, if we can present the following 3
estimates:

  1. 1.

$|\nabla_{\partial\mathcal{M}}\omega({y})|\leq C\delta^{2},\qquad\forall\
{y}\in\partial\mathcal{M}$;

  2. 2.

$|\Delta_{\partial\mathcal{M}}\omega({y})|\leq C\delta^{2},\qquad\forall\
{y}\in\partial\mathcal{M}$;

  3. 3.

$\left\|\nabla_{\partial\mathcal{M}}s\right\|_{L^{2}(\partial\mathcal{M})}+\left\|\Delta_{\partial\mathcal{M}}s\right\|_{L^{2}(\partial\mathcal{M})}\leq
C\delta^{2}\left\|u\right\|_{H^{4}(\mathcal{M})};$

then (A.2) will indicate that

$\begin{split}&\left\|b\right\|_{L^{2}(\partial\mathcal{M})}=\left\|\Delta_{\partial\mathcal{M}}u({y})-\nabla_{\partial\mathcal{M}}\cdot\frac{\nabla_{\partial\mathcal{M}}Pu({y})}{\omega({y})}\right\|_{L_{{y}}^{2}(\partial\mathcal{M})}\\\
&\leq
C(\delta^{-1}\left\|\Delta_{\partial\mathcal{M}}s\right\|_{L^{2}(\partial\mathcal{M})}+\left\|\nabla_{\partial\mathcal{M}}s\right\|_{L^{2}(\partial\mathcal{M})}+\delta\left\|u\right\|_{L^{2}(\partial\mathcal{M})}+\delta\left\|\nabla_{\partial\mathcal{M}}u\right\|_{L^{2}(\partial\mathcal{M})})\\\
&\leq C\delta\left\|u\right\|_{H^{4}(\mathcal{M})};\end{split}$ (A.4)

consequently, the control of $r_{2}$ will be completed. In the following 3
subsections, we will give a proof for the above 3 estimates.

###  A.1 Control on $\nabla_{\partial\mathcal{M}}\omega$

We start from the control of the term
$|\nabla_{\partial\mathcal{M}}\omega({y})|$. For simplicity, for any point
${x}\in\mathcal{M},\ {y}\in\partial\mathcal{M}$, let us denote
$\mathcal{T}_{\partial\mathcal{M}}^{{y}}$ to be the tangent space of
$\partial\mathcal{M}$ at ${y}$, $\mathcal{T}_{\mathcal{M}}^{{x}}$ to be the
tangent space of $\mathcal{M}$ at ${x}$. Moreover, for any
$\textbf{v}\in\mathbb{R}^{d}$, we denote
$\mathcal{T}_{\partial\mathcal{M}}^{{y}}(\textbf{v})$ as the projection of v
onto $\mathcal{T}_{\partial\mathcal{M}}^{{y}}$,
$\mathcal{T}_{\mathcal{M}}^{{x}}(\textbf{v})$ as the projection of v onto
$\mathcal{T}_{\mathcal{M}}^{{x}}$. Finally, we denote $n^{j},x^{j},y^{j}$ to
be the $j^{th}$ component of $\bm{n},{x}$ and ${y}$; $\partial_{{y}}^{i}$ to
be the $i^{th}$ component of $\nabla_{\partial\mathcal{M}}^{{y}}$,
$\partial_{{x}}^{i}$ to be the $i^{th}$ component of
$\nabla_{\mathcal{M}}^{{x}}$.

Apparently,
$\nabla_{\partial\mathcal{M}}\omega({y})\in\mathcal{T}_{\partial\mathcal{M}}^{{y}}$.
In order to seek for a bound on $\nabla_{\partial\mathcal{M}}\omega({y})$, we
only need to control $\nabla_{\partial\mathcal{M}}\omega({y})\cdot\bm{\tau}$
for arbitrary unit vector
$\bm{\tau}\in\mathcal{T}_{\partial\mathcal{M}}^{{y}}$. To this end, we fix a
unit vector $\bm{\tau}\in\mathcal{T}_{\partial\mathcal{M}}^{{y}}$ and denote
$\tau^{j}$ to be the $j^{th}$ component of $\bm{\tau}$. We then calculate

$\begin{split}&\nabla_{\partial\mathcal{M}}\omega({y})\cdot\bm{\tau}=\bm{\tau}\cdot\nabla^{{y}}_{\partial\mathcal{M}}\int_{\mathcal{M}}({x}-{y})\cdot\bm{n}({y})\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
=&\int_{\mathcal{M}}\tau^{i}\partial_{{y}}^{i}(x^{j}-y^{j})n^{j}({y})\bar{R}_{\delta}({x},{y})d\mu_{x}+\int_{\mathcal{M}}\tau^{i}(x^{j}-y^{j})\partial_{{y}}^{i}n^{j}({y})\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
&+\int_{\mathcal{M}}(x^{j}-y^{j})n^{j}({y})\tau^{i}\partial_{{y}}^{i}\bar{R}_{\delta}({x},{y})d\mu_{x}.\end{split}$
(A.5)

Now we analyze each term that composes (A.5), first,

$\int_{\mathcal{M}}\tau^{i}\partial_{{y}}^{i}(x^{j}-y^{j})n^{j}({y})\bar{R}_{\delta}({x},{y})d\mu_{x}=\int_{\mathcal{M}}\tau^{i}n^{i}({y})\bar{R}_{\delta}({x},{y})d\mu_{x}=0.$
(A.6)

For the second term of (A.5), we do the following decomposition and
integration by parts

$\begin{split}&\int_{\mathcal{M}}\tau^{i}(x^{j}-y^{j})\partial_{{y}}^{i}n^{j}({y})\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
=&-\int_{\mathcal{M}}\partial_{{y}}^{i}n^{j}({y})\tau^{i}(({y}^{j}-{x}^{j})\bar{R}_{\delta}({x},{y})-2\delta^{2}\partial_{{x}}^{j}\overset{=}{R}_{\delta}({x},{y}))d\mu_{x}\\\
&-2\delta^{2}\int_{\mathcal{M}}\partial_{{y}}^{i}n^{j}({y})\tau^{i}\partial_{{x}}^{j}\overset{=}{R}_{\delta}({x},{y})d\mu_{x}\\\
=&-\int_{\mathcal{M}}\partial_{{y}}^{i}n^{j}({y})\tau^{i}(({y}^{j}-{x}^{j})\bar{R}_{\delta}({x},{y})-2\delta^{2}\partial_{{x}}^{j}\overset{=}{R}_{\delta}({x},{y}))d\mu_{x}\\\
&+2\delta^{2}\int_{\mathcal{M}}div_{{x}}(\partial_{{y}}^{i}\bm{n}({y})\tau^{i})\overset{=}{R}_{\delta}({x},{y})d\mu_{x}-2\delta^{2}\int_{\partial\mathcal{M}}\partial_{{y}}^{i}n^{j}({y})\tau^{i}n^{j}({x})\overset{=}{R}_{\delta}({x},{y})d\tau_{x}\\\
=&\int_{\mathcal{M}}\tau^{i}\partial_{{y}}^{i}\bm{n}({y})\cdot({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}))\bar{R}_{\delta}({x},{y})d\mu_{x}+2\delta^{2}\int_{\mathcal{M}}div_{{x}}(\partial_{{y}}^{i}\bm{n}({y})\tau^{i})\overset{=}{R}_{\delta}({x},{y})d\mu_{x}\\\
&+2\delta^{2}\int_{\partial\mathcal{M}}\tau^{i}(\partial_{{y}}^{i}\bm{n}({y})\cdot(\bm{n}({y})-\bm{n}({x}))\overset{=}{R}_{\delta}({x},{y})d\tau_{x}=r_{21}({y})+r_{22}({y})+r_{23}({y});\end{split}$
(A.7)

here in the second last equality of (A.7) we have used the fact that
$\partial_{{y}}^{i}\bm{n}({y})\cdot\bm{n}({y})=0$ and
$\nabla^{{x}}_{\mathcal{M}}\bar{R}_{\delta}({x},{y})=-\frac{1}{2\delta^{2}}\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})R_{\delta}({x},{y})$.

Finally, for the last term of (A.5), since we can do the following
decomposition

$\begin{split}\tau^{j}\partial^{j}_{{y}}\bar{R}_{\delta}({x},{y})=&\bm{\tau}\cdot\nabla_{\partial\mathcal{M}}^{{y}}\bar{R}_{\delta}({x},{y})=\frac{1}{2\delta^{2}}\bm{\tau}\cdot\mathcal{T}_{\partial\mathcal{M}}^{{y}}({x}-{y}){R}_{\delta}({x},{y})=\frac{1}{2\delta^{2}}\bm{\tau}\cdot({x}-{y}){R}_{\delta}({x},{y})\\\
=&\frac{1}{2\delta^{2}}\bm{\tau}\cdot({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})+\frac{1}{2\delta^{2}}\bm{\tau}\cdot\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}){R}_{\delta}({x},{y})\\\
=&\frac{1}{2\delta^{2}}\bm{\tau}\cdot({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})-\bm{\tau}\cdot\nabla_{\mathcal{M}}^{{x}}\bar{R}_{\delta}({x},{y})\\\
=&\frac{1}{2\delta^{2}}\bm{\tau}\cdot({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})-\tau^{j}\partial_{{x}}^{j}\bar{R}_{\delta}({x},{y});\end{split}$
(A.8)

we then apply (A.8) and do integration by parts to find

$\begin{split}&\int_{\mathcal{M}}(x^{j}-y^{j})n^{j}({y})\tau^{i}\partial_{{y}}^{i}\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
=&-\int_{\mathcal{M}}(x^{j}-y^{j})n^{j}({y})(\bm{\tau}\cdot\frac{1}{2\delta^{2}}({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}))-\partial_{{x}}^{i}\bar{R}_{\delta}({x},{y}))d\mu_{x}\\\
=&\int_{\mathcal{M}}(x^{j}-y^{j})n^{j}({y})\
\bm{\tau}\cdot\frac{1}{2\delta^{2}}({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})d\mu_{x}\\\
&-\int_{\partial\mathcal{M}}(x^{j}-y^{j})n^{j}({y})\tau^{i}n^{i}({x})\bar{R}_{\delta}({x},{y})d\tau_{x}+\int_{\mathcal{M}}\partial_{{x}}^{i}(x^{j}-y^{j})n^{j}({y})\tau^{i}\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
=&\int_{\mathcal{M}}(x^{j}-y^{j})n^{j}({y})\
(\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\cdot\frac{1}{2\delta^{2}}({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})d\mu_{x}\\\
&-\int_{\partial\mathcal{M}}[({x}-{y}-\mathcal{T}_{\partial\mathcal{M}}^{{y}}({x}-{y}))\cdot\bm{n}({y})]\
[\bm{\tau}\cdot(\bm{n}({x})-\bm{n}({y}))]\bar{R}_{\delta}({x},{y})d\tau_{x}\\\
&+\int_{\mathcal{M}}(\bm{n}({y})-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{n}({y}))\cdot(\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\bar{R}_{\delta}({x},{y})d\mu_{x}=r_{24}({y})+r_{25}({y})+r_{26}({y});\end{split}$
(A.9)

here in the second last equality of (A.9) we have used the fact that
$\bm{n}({y})\cdot\bm{\tau}=0$,
$\bm{n}({y})\cdot\mathcal{T}_{\partial\mathcal{M}}^{{y}}({x}-{y})=0$,
$\partial_{{x}}^{i}(x^{j}-y^{j})n^{j}({y})=(\bm{n}({y})-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{n}({y}))^{i}$;
and
$(\textbf{u}-\mathcal{T}_{\mathcal{M}}^{{x}}(\textbf{u}))\cdot\mathcal{T}_{\mathcal{M}}^{{x}}(\textbf{v})=0$
for all $\textbf{u},\textbf{v}\in\mathbb{R}^{d}$.

In conclusion, we have

$\nabla_{\partial\mathcal{M}}\omega({y})\cdot\bm{\tau}=r_{21}({y})+r_{22}({y})+r_{23}({y})+r_{24}({y})+r_{25}({y})+r_{26}({y}).$
(A.10)

Here $\\{r_{2t}\\}_{t=1}^{6}$ is defined in (A.7) and (A.9). By noticing that
$\int_{\mathcal{M}}\bar{R}_{\delta}({x},{y})d\mu_{x}$ is $\mathcal{O}(1)$,
$\int_{\partial\mathcal{M}}\bar{R}_{\delta}({x},{y})d\tau_{x}$ is
$\mathcal{O}(\delta^{-1})$, $\partial_{{y}}^{i}\bm{n}({y}),\
div_{{x}}(\partial_{{y}}^{i}\bm{n}({y})\tau^{i})$ is $\mathcal{O}(1)$,
$(\bm{n}({y})-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{n}({y})),\
(\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})$ is
$\mathcal{O}(\delta)$, $({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}))$ is
$\mathcal{O}(\delta^{2})$, we are able to conclude that

$|\nabla_{\partial\mathcal{M}}\omega({y})\cdot\bm{\tau}|=|\sum_{t=1}^{6}r_{2t}({y})|\leq
C\delta^{2}\qquad\forall\ {y}\in\partial\mathcal{M},\
\bm{\tau}\in\mathcal{T}_{\partial\mathcal{M}}^{{y}}.$ (A.11)

We then derive $|\nabla_{\partial\mathcal{M}}\omega({y})|\leq C\delta^{2}$ by
the randomness of unit vector $\bm{\tau}$.

###  A.2 Control on $\Delta_{\partial\mathcal{M}}\omega$

For the second step, we start to control $\Delta_{\partial\mathcal{M}}\omega$.
The following equality for Laplace-Beltrami operator is presented in $(4.15)$
of our previous paper [48] (which we also mentioned in (2.5)):

$\begin{split}&\Delta_{\partial\mathcal{M}}\omega({y})=\sum\limits_{l=1}^{m-1}\bm{H}_{\omega}({y}):\bm{\tau_{l}}\otimes\bm{\tau_{l}}=\sum\limits_{l=1}^{m-1}(\partial^{j}_{{y}}\partial^{k}_{{y}}\omega({y}))\tau_{l}^{j}\tau_{l}^{k},\end{split}$
(A.12)

where $\\{\bm{\tau}_{l}\\}_{l=1}^{m-1}$ is an orthonormal basis of the space
$\mathcal{T}_{\partial\mathcal{M}}^{{y}}$. Therefore, to seek for a bound of
$\Delta_{\partial\mathcal{M}}\omega({y})$, we only need to control
$(\partial^{j}_{{y}}\partial^{k}_{{y}}\omega({y}))\tau^{j}\tau^{k}$ for
arbitrary unit vector $\bm{\tau}\in\mathcal{T}_{\partial\mathcal{M}}^{{y}}$.
To this end, we fix the unit vector
$\bm{\tau}\in\mathcal{T}_{\partial\mathcal{M}}^{{y}}$ and concentrate on the
control of $(\partial^{j}_{{y}}\partial^{k}_{{y}}\omega({y}))\tau^{j}\tau^{k}$
in the following content. In fact, (A.10) gives

$\begin{split}(\partial^{j}_{{y}}\partial^{k}_{{y}}\omega({y}))\tau^{j}\tau^{k}=\tau^{j}\partial^{j}_{{y}}(\tau^{k}\partial^{k}_{{y}}\omega({y}))=\tau^{j}\partial^{j}_{{y}}(\nabla_{\partial\mathcal{M}}\omega({y})\cdot\bm{\tau})=\sum\limits_{t=1}^{6}\tau^{j}\partial^{j}_{{y}}r_{2t}({y}).\end{split}$
(A.13)

Now our purpose becomes the control of $\tau^{j}\partial^{j}_{{y}}r_{2t},\
t=1,2,...,6$ where $r_{2t}$ is defined in (A.7) and (A.9). For brevity’s sake,
we write $\mathcal{O}(\delta^{2})$ instead of the terms that are no larger
than $C\delta^{2}$ for some constant $C$ depend only on the shape of
$\mathcal{M}$.

To simplify the calculation, let us introduce the following lemma.

######  Lemma A.1.

  1. 1.

for any $C^{1}$ function
$g_{1}:\mathcal{M}\times\partial\mathcal{M}\to\mathbb{R}$ with
$|g_{1}({x},{y})|\leq C\delta^{2}$ for all $|{x}-{y}|\leq 2\delta$, we have

$\begin{split}\tau^{j}&\partial_{{y}}^{j}\int_{\mathcal{M}}g_{1}({x},{y})\bar{R}({x},{y})d\mu_{{x}}=\int_{\mathcal{M}}\tau^{j}((\partial_{{y}}^{j}+\partial_{{x}}^{j})g_{1}({x},{y}))\bar{R}({x},{y})d\mu_{{x}}+\mathcal{O}(\delta^{2});\end{split}$
(A.14)

  2. 2.

for any $C^{1}$ function
$g_{2}:\partial\mathcal{M}\times\partial\mathcal{M}\to\mathbb{R}$ with
$|g_{2}({x},{y})|\leq C\delta^{3}$ for all $|{x}-{y}|\leq 2\delta$, we have

$\begin{split}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\tau^{j}&\partial_{{y}}^{j}\int_{\partial\mathcal{M}}g_{2}({x},{y})\bar{R}({x},{y})d\mu_{{x}}=\int_{\partial\mathcal{M}}\tau^{j}((\partial_{{y}}^{j}+\partial_{{x},\partial\mathcal{M}}^{j})g_{2}({x},{y}))\bar{R}({x},{y})d\mu_{{x}}+\mathcal{O}(\delta^{2}),\end{split}$
(A.15)

where $\partial_{{x},\partial\mathcal{M}},{x}\in\partial\mathcal{M}$ is the
gradient operator on $\partial\mathcal{M}$.

The kernel function $\bar{R}_{\delta}$ in (A.14) and (A.15) can be replaced by
the other kernel functions $R_{\delta}$ and $\overset{=}{R}_{\delta}$.

######  Proof A.2.

we apply (A.8) and do integration by parts to obtain

$\begin{split}&\tau^{j}\partial_{{y}}^{j}\int_{\mathcal{M}}g_{1}({x},{y})\bar{R}({x},{y})d\mu_{{x}}\\\
=&\int_{\mathcal{M}}(\tau^{j}\partial_{{y}}^{j}g_{1}({x},{y}))\bar{R}({x},{y})d\mu_{{x}}+\int_{\mathcal{M}}g_{1}({x},{y})(\tau^{j}\partial_{{y}}^{j}\bar{R}({x},{y}))d\mu_{{x}}\\\
=&\int_{\mathcal{M}}(\tau^{j}\partial_{{y}}^{j}g_{1}({x},{y}))\bar{R}({x},{y})d\mu_{{x}}-\int_{\mathcal{M}}g_{1}({x},{y})(\tau^{j}\partial_{{x}}^{j}\bar{R}({x},{y}))d\mu_{{x}}\\\
&+\int_{\mathcal{M}}g_{1}({x},{y})\frac{1}{2\delta^{2}}\bm{\tau}\cdot({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})d\mu_{{x}}\\\
=&\int_{\mathcal{M}}\tau^{j}((\partial_{{y}}^{j}+\partial_{{x}}^{j})g_{1}({x},{y}))\bar{R}({x},{y})d\mu_{{x}}-\int_{\partial\mathcal{M}}g_{1}({x},{y})\bar{R}({x},{y})\tau^{j}n^{j}({x})d\tau_{{x}}\\\
&+\int_{\mathcal{M}}g_{1}({x},{y})\frac{1}{2\delta^{2}}\bm{\tau}\cdot({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})d\mu_{{x}}\end{split}$
(A.16)

since $g_{1}({x},{y})$ and
$({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}))$ is
$\mathcal{O}(\delta^{2})$,
$\tau^{j}n^{j}({x})=\tau^{j}(n^{j}({x})-n^{j}({y}))$ is $\mathcal{O}(\delta)$,
we can derive that the last two terms of (A.16) is $\mathcal{O}(\delta^{2})$.
This is exactly (A.14).

By a similar argument, we obtain

$\begin{split}\tau^{j}&\partial_{{y}}^{j}\int_{\partial\mathcal{M}}g_{2}({x},{y})\bar{R}({x},{y})d\mu_{{x}}=\int_{\partial\mathcal{M}}\tau^{j}((\partial_{{y}}^{j}+\partial_{{x},\partial\mathcal{M}}^{j})g_{2}({x},{y}))\bar{R}({x},{y})d\tau_{{x}}\\\
&+\int_{\partial\mathcal{M}}g_{2}({x},{y})\frac{1}{2\delta^{2}}\bm{\tau}\cdot({x}-{y}-\mathcal{T}_{\partial\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})d\tau_{{x}};\end{split}$
(A.17)

the last of (A.17) is $\mathcal{O}(\delta^{2})$ since $g_{2}({x},{y})$ is
$\mathcal{O}(\delta^{3})$ and
$({x}-{y}-\mathcal{T}_{\partial\mathcal{M}}^{{x}}({x}-{y}))$ is
$\mathcal{O}(\delta^{2})$. we then obtain (A.15).

Now we move back to control each term of (A.13) by utilizing Lemma A.1. Still,
we divide our calculation into 6 parts.

  1. 1.

To control the first term of (A.13), we need to analyze the following term in
front:

$\begin{split}&\tau^{j}(\partial^{j}_{{y}}+\partial^{j}_{{x}})({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}))\\\
=&\tau^{j}\partial^{j}_{{y}}({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}))+\tau^{j}\partial^{j}_{{x}}({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}))\\\
=&(-\bm{\tau}+\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})+(\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}-\tau^{j}\partial^{j}_{{x}}(\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})))\\\
=&-\bm{\tau}+\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}-({x}^{j}-{y}^{j})\partial^{j}_{{x}}(\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})=-\bm{\Phi}({y})+\bm{\Phi}({x})-({x}-{y})\cdot\nabla_{\mathcal{M}}^{{x}}\bm{\Phi}({x})\end{split}$
(A.18)

here
$\bm{\Phi}({x})=\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau},{x}\in\mathcal{M}$ is
a $C^{2}$ vector-valued function, and the second last equality of (A.18) is
from the following calculation: for any $\textbf{v}\in\mathbb{R}^{d}$,

$\begin{split}&\tau^{j}\partial^{j}_{{x}}(\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}))\cdot\textbf{v}=\tau^{j}\partial^{j}_{{x}}(\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})\cdot\textbf{v})=\tau^{j}\partial^{j}_{{x}}(({x}-{y})\cdot\mathcal{T}_{\mathcal{M}}^{{x}}\textbf{v})\\\
&=\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}\cdot\mathcal{T}_{\mathcal{M}}^{{x}}\textbf{v}+({x}-{y})\cdot\tau^{j}\partial^{j}_{{x}}\mathcal{T}_{\mathcal{M}}^{{x}}\textbf{v}=\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}\cdot\textbf{v}+({x}-{y})\otimes\bm{\tau}:\nabla^{{x}}_{\mathcal{M}}\nabla^{{x}}_{\mathcal{M}}(\textbf{v}\cdot{x})\\\
&=\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}\cdot\textbf{v}+\bm{\tau}\cdot({x}^{j}-{y}^{j})\partial^{j}_{{x}}\mathcal{T}_{\mathcal{M}}^{{x}}\textbf{v}=\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}\cdot\textbf{v}+({x}^{j}-{y}^{j})\partial^{j}_{{x}}(\mathcal{T}_{\mathcal{M}}^{{x}}\textbf{v}\cdot\bm{\tau})\\\
&=(\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}+({x}^{j}-{y}^{j})\partial^{j}_{{x}}(\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}))\cdot\textbf{v}.\end{split}$
(A.19)

We then conclude that
$\tau^{j}(\partial^{j}_{{y}}+\partial^{j}_{{x}})({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}))$
is $\mathcal{O}(\delta^{2})$. Thereafter, we apply (A.14) on
$\tau^{j}\partial^{j}_{{y}}r_{21}$ to find

$\begin{split}\tau^{j}\partial^{j}_{{y}}&r_{21}({y})=\tau^{j}\partial^{j}_{{y}}\int_{\mathcal{M}}\tau^{i}\partial_{{y}}^{i}\bm{n}({y})\cdot({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}))\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
=&\int_{\mathcal{M}}\tau^{i}\tau^{j}\partial_{{y}}^{i}\partial^{j}_{{y}}\bm{n}({y})\cdot({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}))\bar{R}_{\delta}({x},{y})d\mu_{x}+\mathcal{O}(\delta^{2})\\\
&+\int_{\mathcal{M}}\tau^{i}\partial_{{y}}^{i}\bm{n}({y})\cdot\tau^{j}(\partial^{j}_{{y}}+\partial^{j}_{{x}})({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}))\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
=&\mathcal{O}(\delta^{2})+\mathcal{O}(\delta^{2})+\mathcal{O}(\delta^{2})=\mathcal{O}(\delta^{2}),\end{split}$
(A.20)

We then complete the control of $\tau^{j}\partial^{j}_{{y}}r_{21}$.

  2. 2.

For the term $\tau^{j}\partial^{j}_{{y}}r_{22}$, we apply again (A.14) to
discover

$\begin{split}&\tau^{j}\partial^{j}_{{y}}r_{22}=2\delta^{2}\tau^{j}\partial^{j}_{{y}}\int_{\mathcal{M}}div_{{x}}(\partial_{{y}}^{i}\bm{n}({y})\tau^{i})\overset{=}{R}_{\delta}({x},{y})d\mu_{x}\\\
=&2\delta^{2}\int_{\mathcal{M}}div_{{x}}(\partial_{{y}}^{i}\partial^{j}_{{y}}\bm{n}({y})\tau^{i}\tau^{j})\overset{=}{R}_{\delta}({x},{y})d\mu_{{x}}+\mathcal{O}(\delta^{2})\\\
&+2\delta^{2}\int_{\mathcal{M}}\partial^{j}_{{x}}(div_{{x}}(\partial_{{y}}^{i}\bm{n}({y})\tau^{i}))\tau^{j}\overset{=}{R}_{\delta}({x},{y})d\mu_{x}=\mathcal{O}(\delta^{2})\end{split}$
(A.21)

here the $\mathcal{O}(1)$ control on the terms
$\partial^{j}_{{x}}(div_{{x}}\textbf{v})$ and $div_{{x}}\textbf{v}$ for
$\textbf{v}=\partial_{{y}}^{i}\bm{n}({y})\tau^{i}$ is by the smoothness of
$\mathcal{M}$. We thereby complete the control of
$\tau^{j}\partial^{j}_{{y}}r_{22}$.

  3. 3.

To estimate the next term $\tau^{j}\partial^{j}_{{y}}r_{23}$, we apply (A.15)
to discover

$\begin{split}\tau^{j}\partial^{j}_{{y}}r_{23}=&2\delta^{2}\tau^{j}\partial_{{y}}^{j}\int_{\partial\mathcal{M}}\tau^{i}(\partial_{{y}}^{i}\bm{n}({y})\cdot(\bm{n}({y})-\bm{n}({x}))\overset{=}{R}_{\delta}({x},{y})d\tau_{x}\\\
=&2\delta^{2}\int_{\partial\mathcal{M}}\tau^{i}\tau^{j}\partial_{{y}}^{i}\partial_{{y}}^{j}\bm{n}({y})\cdot(\bm{n}({y})-\bm{n}({x}))\overset{=}{R}_{\delta}({x},{y})d\tau_{x}+\mathcal{O}(\delta^{2})\\\
&+2\delta^{2}\int_{\partial\mathcal{M}}(\tau^{i}\partial_{{y}}^{i}\bm{n}({y})-\tau^{j}\partial_{{x},\partial\mathcal{M}}^{j}\bm{n}({x}))\cdot(\tau^{j}\partial_{{y}}^{j}\bm{n}({y}))\overset{=}{R}_{\delta}({x},{y})d\tau_{x}\end{split}$
(A.22)

by the ${C}^{2}$ smoothness of $\bm{n}$ through $\partial\mathcal{M}$, we
obtain $\tau^{j}\partial^{j}_{{y}}r_{23}$ is $\mathcal{O}(\delta^{2})$.

  4. 4.

Now we move on to the control on $\tau^{j}\partial^{j}_{{y}}r_{24}$. We apply
(A.14) to find

$\begin{split}&\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\tau^{j}\partial^{j}_{{y}}r_{24}({y})\\\
&\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=\tau^{j}\partial^{j}_{{y}}\int_{\mathcal{M}}(x^{i}-y^{i})n^{i}({y})\
(\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\cdot\frac{1}{2\delta^{2}}({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})d\mu_{x}\\\
&\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=\int_{\mathcal{M}}\tau^{j}((\partial^{j}_{{y}}+\partial^{j}_{{x}})(x^{i}-y^{i}))n^{i}({y})\
(\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\cdot\frac{1}{2\delta^{2}}({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})d\mu_{x}\\\
&\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+\int_{\mathcal{M}}(x^{i}-y^{i})\tau^{j}\partial^{j}_{{y}}n^{i}({y})\
(\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\cdot\frac{1}{2\delta^{2}}({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})d\mu_{x}\\\
&\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-\int_{\mathcal{M}}(x^{i}-y^{i})n^{i}({y})\tau^{j}(\partial^{j}_{{x}}\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\cdot\frac{1}{2\delta^{2}}({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})d\mu_{x}+\mathcal{O}(\delta^{2})\\\
&\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+\int_{\mathcal{M}}(x^{i}-y^{i})n^{i}({y})\
(\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\cdot\frac{1}{2\delta^{2}}\tau^{j}(\partial^{j}_{{y}}+\partial^{j}_{{x}})({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})d\mu_{x}\\\
&\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=\
r_{241}({y})+r_{242}({y})+r_{243}({y})+r_{244}({y})+\mathcal{O}(\delta^{2}),\end{split}$
(A.23)

here $r_{242}({y})$ is apparently $\mathcal{O}(\delta^{2})$, and

$\begin{split}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!r_{241}({y})=\int_{\mathcal{M}}\tau^{j}(n^{j}({y})-\mathcal{T}_{\mathcal{M}}^{{x}}n^{j}({y}))\
(\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\cdot\frac{1}{2\delta^{2}}({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})d\mu_{x},\end{split}$
(A.24)

which is $\mathcal{O}(\delta^{2})$ as well. In addition, since for any
$\textbf{v}\in\mathbb{R}^{d}$,

$\begin{split}(\tau^{j}\partial_{{x}}^{j}\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\cdot\textbf{v}=\bm{\tau}\otimes\textbf{v}:\nabla^{{x}}_{\mathcal{M}}\nabla^{{x}}_{\mathcal{M}}(\bm{\tau}\cdot{x})=(v^{j}\partial_{{x}}^{j}\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\cdot\bm{\tau}=v^{j}\partial_{{x}}^{j}(\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}\cdot\bm{\tau})\\\
=-v^{j}\partial_{{x}}^{j}((\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}-\bm{\tau})\cdot(\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}-\bm{\tau}))=-2v^{j}(\partial_{{x}}^{j}\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\cdot(\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}-\bm{\tau}),\end{split}$
(A.25)

while
$\partial_{{x}}^{j}\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}=\partial^{j}_{{x}}\bm{\Phi}({x})\in
C^{1}(\mathcal{M})$, we can derive
$\tau^{j}\partial_{{x}}^{j}\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}$ is
$\mathcal{O}(\delta)$, which indicates that $r_{243}({y})$ is
$\mathcal{O}(\delta^{2})$. Moreover, we apply (A.18) to find

$\begin{split}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!r_{244}({y})=\int_{\mathcal{M}}(x^{i}-y^{i})n^{i}({y})\
(\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\cdot\frac{1}{2\delta^{2}}(-\bm{\Phi}({y})+\bm{\Phi}({x})-({x}-{y})\cdot\nabla_{\mathcal{M}}^{{x}}\bm{\Phi}({x})){R}_{\delta}({x},{y})d\mu_{x},\end{split}$
(A.26)

hence $r_{244}({y})$ is $\mathcal{O}(\delta^{2})$ as well. We then conclude
that $\tau^{j}\partial^{j}_{{y}}r_{24}({y})$ is $\mathcal{O}(\delta^{2})$.

  5. 5.

Next, for $\tau^{j}\partial_{{y}}^{j}r_{25}$, we apply (A.15) to obtain

$\begin{split}&\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\tau^{j}\partial_{{y}}^{j}r_{25}({y})=\tau^{j}\partial_{{y}}^{j}\int_{\partial\mathcal{M}}[({x}-{y}-\mathcal{T}_{\partial\mathcal{M}}^{{y}}({x}-{y}))\cdot\bm{n}({y})]\
[\bm{\tau}\cdot(\bm{n}({x})-\bm{n}({y}))]\bar{R}_{\delta}({x},{y})d\tau_{x}\\\
&\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=\int_{\partial\mathcal{M}}[\tau^{j}(\partial_{{y}}^{j}+\partial_{{x},\partial\mathcal{M}}^{j})({x}-{y}-\mathcal{T}_{\partial\mathcal{M}}^{{y}}({x}-{y}))]\cdot\bm{n}({y})\
[\bm{\tau}\cdot(\bm{n}({x})-\bm{n}({y}))]\bar{R}_{\delta}({x},{y})d\tau_{x}\\\
&\\!\\!\\!\\!\\!+\int_{\partial\mathcal{M}}[({x}-{y}-\mathcal{T}_{\partial\mathcal{M}}^{{y}}({x}-{y}))\cdot\tau^{j}\partial_{{y}}^{j}\bm{n}({y})]\
[\bm{\tau}\cdot(\bm{n}({x})-\bm{n}({y}))]\bar{R}_{\delta}({x},{y})d\tau_{x}+\mathcal{O}(\delta^{2})\\\
&\\!\\!\\!\\!\\!+\int_{\partial\mathcal{M}}[({x}-{y}-\mathcal{T}_{\partial\mathcal{M}}^{{y}}({x}-{y}))\cdot\bm{n}({y})]\
[\bm{\tau}\cdot(\tau^{j}\partial_{{x},\partial\mathcal{M}}^{j}\bm{n}({x})-\tau^{j}\partial_{{y}}^{j}\bm{n}({y}))]\bar{R}_{\delta}({x},{y})d\tau_{x}\\\
&\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=r_{251}({y})+r_{252}({y})+r_{253}({y})+\mathcal{O}(\delta^{2});\end{split}$
(A.27)

here $r_{252}({y})$ is apparently $\mathcal{O}(\delta^{2})$, $r_{253}({y})$ is
$\mathcal{O}(\delta^{2})$ by the continuity of the vector function
$\tau^{j}\partial_{{y}}^{j}\bm{n}({y})$; for $r_{251}({y})$, we repeat the
calculation (A.18) on the $(m-1)$ dimensional manifold $\partial\mathcal{M}$
instead of $\mathcal{M}$ to obtain

$\begin{split}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!r_{251}({y})=\int_{\partial\mathcal{M}}[-\bm{\phi}({y})+\bm{\phi}({x})-({x}-{y})\cdot\nabla_{\partial\mathcal{M}}^{{y}}\bm{\phi}({y})]\cdot\bm{n}({y})\
[\bm{\tau}\cdot(\bm{n}({x})-\bm{n}({y}))]\bar{R}_{\delta}({x},{y})d\tau_{x},\end{split}$
(A.28)

where
$\bm{\phi}({x})=\mathcal{T}_{\partial\mathcal{M}}^{{x}}\bm{\tau},{x}\in\partial\mathcal{M}$;
this is obviously $\mathcal{O}(\delta^{2})$ by the smoothness of $\bm{\phi}$.
We thereby conclude that $\tau^{j}\partial_{{y}}^{j}r_{25}({y})$ is
$\mathcal{O}(\delta^{2})$.

  6. 6.

Finally, we move on to the control of $\tau^{j}\partial_{{y}}^{j}r_{26}$.
Still, we apply (A.14) to find

$\begin{split}\tau^{j}\partial_{{y}}^{j}&r_{26}({y})=\tau^{j}\partial_{{y}}^{j}\int_{\mathcal{M}}(\bm{n}({y})-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{n}({y}))\cdot(\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
=&\int_{\mathcal{M}}[\tau^{j}\partial_{{y}}^{j}\bm{n}({y})-\mathcal{T}_{\mathcal{M}}^{{x}}(\tau^{j}\partial_{{y}}^{j}\bm{n}({y}))]\cdot(\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
&+\mathcal{O}(\delta^{2})+\int_{\mathcal{M}}-\tau^{j}\partial_{{x}}^{j}\mathcal{T}_{\mathcal{M}}^{{x}}\bm{n}({y})\cdot(\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
&+\int_{\mathcal{M}}(\bm{n}({y})-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{n}({y}))\cdot(-\tau^{j}\partial_{{x}}^{j}\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
=&r_{261}({y})+r_{262}({y})+r_{263}({y})+\mathcal{O}(\delta^{2}),\end{split}$
(A.29)

for $r_{262}$, we apply the following equality: for any
$\textbf{v}\in\mathbb{R}^{d}$,

$\begin{split}(\tau^{j}\partial_{{x}}^{j}\mathcal{T}_{\mathcal{M}}^{{x}}\bm{n}({y}))\cdot\textbf{v}=v^{j}\partial_{{x}}^{j}(\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}\cdot\bm{n}({y}))=-v^{j}\partial_{{x}}^{j}((\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}-\bm{\tau})\cdot(\mathcal{T}_{\mathcal{M}}^{{x}}\bm{n}({y})-\bm{n}({y})))\\\
=-v^{j}(\partial_{{x}}^{j}\mathcal{T}_{\mathcal{M}}^{{x}}\bm{n}({y}))\cdot(\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}-\bm{\tau})-v^{j}(\partial_{{x}}^{j}\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\cdot(\mathcal{T}_{\mathcal{M}}^{{x}}\bm{n}({y})-\bm{n}({y})),\end{split}$
(A.30)

hence $\tau^{j}\partial_{{x}}^{j}\mathcal{T}_{\mathcal{M}}^{{x}}\bm{n}({y})$
is $\mathcal{O}(\delta)$ by the smoothness of
$\partial_{{x}}^{j}\mathcal{T}_{\mathcal{M}}^{{x}}\bm{n}({y})$ and
$\partial_{{x}}^{j}\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}$. We then conclude
that $r_{262}({y})$ is $\mathcal{O}(\delta^{2})$. Similarly, we apply (A.25)
to obtain that $r_{263}({y})$ is $\mathcal{O}(\delta^{2})$. The term
$r_{261}({y})$ is more involved, we can calculate

$\begin{split}&r_{261}({y})=\int_{\mathcal{M}}[\tau^{j}\partial_{{y}}^{j}\bm{n}({y})-\mathcal{T}_{\mathcal{M}}^{{x}}(\tau^{j}\partial_{{y}}^{j}\bm{n}({y}))]\cdot(\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
=&\int_{\mathcal{M}}\tau^{j}\partial_{{y}}^{j}\bm{n}({y})\cdot(\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
=&\int_{\mathcal{M}}\tau^{j}\partial_{{y}}^{j}\bm{n}({y})\cdot(\mathcal{T}_{\mathcal{M}}^{{y}}\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau}-({y}^{i}-{x}^{i})\partial_{{y}}^{i}\mathcal{T}_{\mathcal{M}}^{{y}}\bm{\tau})\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
&+\int_{\mathcal{M}}\tau^{j}\partial_{{y}}^{j}\bm{n}({y})\cdot(\partial_{{y}}^{i}\mathcal{T}_{\mathcal{M}}^{{y}}\bm{\tau})({y}^{i}-{x}^{i})\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
=&\int_{\mathcal{M}}\tau^{j}\partial_{{y}}^{j}\bm{n}({y})\cdot(\bm{\Phi}({y})-\bm{\Phi}({x})-({y}-{x})\cdot\nabla_{\mathcal{M}}^{{y}}\bm{\Phi}({y}))\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
&+\int_{\mathcal{M}}\tau^{j}\partial_{{y}}^{j}\bm{n}({y})\cdot(\partial_{{y}}^{i}\mathcal{T}_{\mathcal{M}}^{{y}}\bm{\tau})(({y}^{i}-{x}^{i})\bar{R}_{\delta}({x},{y})-\partial^{i}_{{x}}\overset{=}{R}_{\delta}({x},{y}))d\mu_{x}\\\
&+\int_{\mathcal{M}}\tau^{j}\partial_{{y}}^{j}\bm{n}({y})\cdot(\partial_{{y}}^{i}\mathcal{T}_{\mathcal{M}}^{{y}}\bm{\tau})\partial^{i}_{{x}}\overset{=}{R}_{\delta}({x},{y})d\mu_{x}\\\
=&\mathcal{O}(\delta^{2})+\int_{\mathcal{M}}\tau^{j}\partial_{{y}}^{j}\bm{n}({y})\cdot(\partial_{{y}}^{i}\mathcal{T}_{\mathcal{M}}^{{y}}\bm{\tau})(({y}-{x}-\mathcal{T}_{\mathcal{M}}^{{y}}({y}-{x}))^{i}\bar{R}_{\delta}({x},{y}))d\mu_{x}\\\
&+2\delta^{2}\int_{\partial\mathcal{M}}\tau^{j}\partial_{{y}}^{j}\bm{n}({y})\cdot(\partial_{{y}}^{i}\mathcal{T}_{\mathcal{M}}^{{y}}\bm{\tau})n^{i}({x})\overset{=}{R}_{\delta}({x},{y})d\tau_{x}\\\
&-2\delta^{2}\int_{\mathcal{M}}div(\tau^{j}\partial_{{y}}^{j}\bm{n}({y})\cdot(\partial_{{y}}\mathcal{T}_{\mathcal{M}}^{{y}}\bm{\tau}))\overset{=}{R}_{\delta}({x},{y})d\mu_{x}\\\
=&\mathcal{O}(\delta^{2})+\mathcal{O}(\delta^{2})+2\delta^{2}\int_{\partial\mathcal{M}}\tau^{j}\partial_{{y}}^{j}\bm{n}({y})\cdot(\partial_{{y}}^{i}\mathcal{T}_{\mathcal{M}}^{{y}}\bm{\tau})(n^{i}({x})-n^{i}({y}))\overset{=}{R}_{\delta}({x},{y})d\tau_{x}\\\
=&\mathcal{O}(\delta^{2}).\end{split}$ (A.31)

Hence we complete the control of $\tau^{j}\partial_{{y}}^{j}r_{26}({y})$.

We then combine all the above estimates on
$\tau^{j}\partial_{{y}}^{j}r_{2t}({y}),t=1,2,...,6$ to conclude that
$|\Delta_{\partial\mathcal{M}}\omega({y})|\leq C\delta^{2}$ according to
(A.12) and (A.13).

###  A.3 Control of
$\nabla_{\partial\mathcal{M}}s,\Delta_{\partial\mathcal{M}}s.$

Next, we move on to the control of $\nabla_{\partial\mathcal{M}}s({y})$ and
$\Delta_{\partial\mathcal{M}}s({y})$ where $s({y})$ is defined in (A.3).
First, for any unit vector
$\bm{\tau}\in\mathcal{T}_{\partial\mathcal{M}}^{{y}}$, we apply (A.16) in
Lemma A.1 to derive

$\begin{split}\bm{\tau}\cdot&\nabla_{\partial\mathcal{M}}s({y})=-\tau^{i}\partial_{{y}}^{i}\int_{\mathcal{M}}(u({x})-u({y}))(({x}-{y})\cdot\bm{n}({y}))\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
=&-\int_{\mathcal{M}}\tau^{i}(\partial_{{x}}^{i}u({x})-\partial_{{y}}^{i}u({y}))(({x}-{y})\cdot\bm{n}({y}))\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
&+\int_{\mathcal{M}}(u({x})-u({y}))(\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\cdot\bm{n}({y})\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
&-\int_{\mathcal{M}}(u({x})-u({y}))({x}-{y})\cdot(\tau^{i}\partial_{{y}}^{i}\bm{n}({y}))\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
&-\int_{\partial\mathcal{M}}(u({x})-u({y}))(({x}-{y})\cdot\bm{n}({y}))\bar{R}_{\delta}({x},{y})\tau^{i}n^{i}({x})d\tau_{x}\\\
&+\int_{\mathcal{M}}(u({x})-u({y}))(({x}-{y})\cdot\bm{n}({y}))\
\bm{\tau}\cdot\frac{1}{2\delta^{2}}({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})d\mu_{x}\\\
=&r_{31}({y})+r_{32}({y})+r_{33}({y})+r_{34}({y})+r_{35}({y}).\end{split}$
(A.32)

Secondly, we study the following term

$\begin{split}&\tau^{i}\tau^{j}\partial_{{y}}^{i}\partial_{{y}}^{j}s({y})=\tau^{j}\partial_{{y}}^{j}(\bm{\tau}\cdot\nabla_{\partial\mathcal{M}}s({y}))=\tau^{j}\partial_{{y}}^{j}(r_{31}({y})+r_{32}({y})+r_{33}({y})+r_{34}({y})+r_{35}({y}));\end{split}$
(A.33)

we apply (A.16) and (A.17) in Lemma A.1 again on each term that composes
(A.33) and utilize the results (A.18) (A.25) (A.28) (A.30) to obtain

$\begin{split}&\tau^{i}\tau^{j}\partial_{{y}}^{i}\partial_{{y}}^{j}s({y})=\tau^{j}\partial_{{y}}^{j}(r_{31}({y})+r_{32}({y})+r_{33}({y})+r_{34}({y})+r_{35}({y}))\\\
&=-\int_{\mathcal{M}}\tau^{i}\tau^{j}(\partial_{{x}}^{i}\partial_{{x}}^{j}u({x})-\partial_{{y}}^{i}\partial_{{y}}^{j}u({y}))(({x}-{y})\cdot\bm{n}({y}))\bar{R}_{\delta}({x},{y})d\mu_{x}\\\
&+2\int_{\mathcal{M}}\tau^{i}(\partial_{{x}}^{i}u({x})-\partial_{{y}}^{i}u({y}))[(\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\cdot\bm{n}({y})\bar{R}_{\delta}({x},{y})\\\
&-({x}-{y})\cdot(\tau^{i}\partial_{{y}}^{i}\bm{n}({y}))\bar{R}_{\delta}({x},{y})+(({x}-{y})\cdot\bm{n}({y}))\
\bm{\tau}\cdot\frac{1}{2\delta^{2}}({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})]d\mu_{x}\\\
&-\int_{\partial\mathcal{M}}\tau^{i}(\partial_{{x}}^{i}u({x})+\partial_{{x},\partial\mathcal{M}}^{i}u({x})-2\partial_{{y}}^{i}u({y}))(({x}-{y})\cdot\bm{n}({y}))\bar{R}_{\delta}({x},{y})\tau^{j}(n^{j}({x})-n^{j}({y}))d\tau_{x}\\\
&+2\int_{\mathcal{M}}(u({x})-u({y}))[(\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\cdot\tau^{j}\partial_{{y}}^{j}\bm{n}({y})\bar{R}_{\delta}({x},{y})\\\
&+(({x}-{y})\cdot(\tau^{i}\partial_{{y}}^{i}\bm{n}({y}))-(\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\cdot\bm{n}({y}))(\bm{\tau}\cdot\frac{1}{2\delta^{2}}({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})]d\mu_{x}\\\
&-\int_{\partial\mathcal{M}}(u({x})-u({y}))\tau^{j}(n^{j}({x})-n^{j}({y}))[2(({x}-{y})\cdot\tau^{i}\partial_{{y}}^{i}\bm{n}({y})-(\bm{\tau}-\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\cdot\bm{n}({y}))\bar{R}_{\delta}({x},{y})\\\
&+(({x}-{y})\cdot\bm{n}({y}))\bm{\tau}\cdot\frac{1}{2\delta^{2}}(2({x}-{y})-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y})-\mathcal{T}_{\partial\mathcal{M}}^{{x}}({x}-{y})){R}_{\delta}({x},{y})]d\tau_{x}\\\
&-\int_{\mathcal{M}}(u({x})-u({y}))[(\tau^{j}(\partial_{{x}}^{j}\mathcal{T}_{\mathcal{M}}^{{x}}\bm{\tau})\cdot\bm{n}({y})+({x}-{y})\cdot\tau^{i}\tau^{j}\partial_{{y}}^{i}\partial_{{y}}^{j}\bm{n}({y}))\bar{R}_{\delta}({x},{y})+({x}-{y})\cdot\bm{n}({y})\\\
&((\frac{\bm{\tau}}{2\delta^{2}}\cdot({x}-{y}-\mathcal{T}_{\mathcal{M}}^{{x}}({x}-{y}))^{2}\underline{R}_{\delta}({x},{y})+\frac{\bm{\tau}}{2\delta^{2}}\cdot(\bm{\Phi}({y})-\bm{\Phi}({x})-({y}-{x})\cdot\nabla_{\mathcal{M}}^{{x}}\bm{\Phi}({x})){R}_{\delta}({x},{y}))]\\\
&d\mu_{x}-\int_{\partial\mathcal{M}}(u({x})-u({y}))(({x}-{y})\cdot\bm{n}({y}))\bar{R}_{\delta}({x},{y})\tau^{i}\tau^{j}(\partial_{{x},\partial\mathcal{M}}^{j}n^{i}({x})-\partial_{{y}}^{j}n^{i}({y}))d\tau_{x}.\end{split}$
(A.34)

Now let us introduce $\nabla_{\partial\mathcal{M}}s$ and
$\Delta_{\partial\mathcal{M}}s$. By the equality (A.12) and the fact
$\nabla_{\partial\mathcal{M}}s({y})\in\mathcal{T}_{\partial\mathcal{M}}^{{y}}$,
let $\\{\bm{\tau}_{l}\\}_{l=1}^{m-1}$ be an orthonormal basis of the space
$\mathcal{T}_{\partial\mathcal{M}}^{{y}}$, then

$\nabla_{\partial\mathcal{M}}s({y})=\sum\limits_{l=1}^{m-1}(\bm{\tau}_{l}\cdot\nabla_{\partial\mathcal{M}}s({y}))\bm{\tau}_{l},\qquad\Delta_{\partial\mathcal{M}}s({y})=\sum\limits_{l=1}^{m-1}\tau^{i}_{l}\tau^{j}_{l}\partial_{{y}}^{i}\partial_{{y}}^{j}s({y}).$
(A.35)

Consequently, we can write down the complete formulation of
$\nabla_{\partial\mathcal{M}}s({y})$ and $\Delta_{\partial\mathcal{M}}s({y})$
by plugging (A.32) and (A.34) into the equation (A.35).

Our next goal is the control of (A.35). In fact, by observing (A.32) (A.34)
(A.35) and applying the equality

$\tau^{i}\partial_{{y}}^{i}u({y})=\tau^{i}\partial_{{y},\mathcal{M}}^{i}u({y}),\
\tau^{i}\tau^{j}\partial_{{y}}^{i}\partial_{{y}}^{j}u({y})=\tau^{i}\tau^{j}\partial_{{y},\mathcal{M}}^{i}\partial_{{y},\mathcal{M}}^{j}u({y}),\qquad\forall\
\bm{\tau}\in\mathcal{T}_{\partial\mathcal{M}}^{{y}},$

we can check that the terms that compose $\nabla_{\partial\mathcal{M}}s({y})$
and $\Delta_{\partial\mathcal{M}}s({y})$ are in one of the following two
forms:

  1. 1.

$\int_{\mathcal{M}}g_{1}({x},{y},\bm{\tau}_{l})(v_{1}({x})-v_{1}({y}))\tilde{R}_{\delta}({x},{y})d\mu_{{x}},$
(A.36)

  2. 2.

$\int_{\partial\mathcal{M}}g_{2}({x},{y},\bm{\tau}_{l})(v_{2}({x})-v_{2}({y}))\tilde{R}_{\delta}({x},{y})d\tau_{{x}};$
(A.37)

where the function $v_{1}$ is one of $u$, $\nabla_{\mathcal{M}}^{i}u$ and
$\nabla_{\mathcal{M}}^{i}\nabla_{\mathcal{M}}^{j}u$; $v_{2}$ is one of $u$,
$\nabla_{\mathcal{M}}^{i}u$ and $\nabla_{\partial\mathcal{M}}^{i}u$; the
kernel function $\tilde{R}_{\delta}$ is one of $\bar{R}_{\delta},\ R_{\delta}$
and $\underline{R}_{\delta}$; $l\in\\{1,2,...,m-1\\}$; in addition,
$g_{i}:\mathcal{M}\times\partial\mathcal{M}\times\mathbb{R}^{d}\to\mathbb{R}\
\mbox{or}\ \mathbb{R}^{d}$ is a function such that
$|g_{i}({x},{y},\bm{\tau})|\leq C\delta^{i};\ i=1,2$ when $|{x}-{y}|\leq
2\delta$ and $\bm{\bm{\tau}}\in\mathcal{T}_{\partial\mathcal{M}}^{{y}}$ .

Such observation implies that

$\begin{split}|\nabla_{\partial\mathcal{M}}s({y})|+|\Delta_{\partial\mathcal{M}}s({y})|\leq&\sum_{i=1}^{s_{1}}\int_{\mathcal{M}}\delta|v_{1}({x})-v_{1}({y})|\tilde{R}_{\delta}({x},{y})d\mu_{{x}}\\\
&+\sum_{i=1}^{s_{2}}\int_{\partial\mathcal{M}}\delta^{2}|v_{2}({x})-v_{2}({y})|\tilde{R}_{\delta}({x},{y})d\tau_{{x}}.\end{split}$
(A.38)

where $s_{1}$ and $s_{2}$ are the number of terms that composes
$\nabla_{\partial\mathcal{M}}s({y})$ and $\Delta_{\partial\mathcal{M}}s({y})$.

Now we aim to control the $L_{{y}}^{2}(\partial\mathcal{M})$ norm of the right
hand side of (A.38). We will need the following lemma.

######  Lemma A.3.

  1. 1.

let the function $v_{1}\in H^{2}(\mathcal{M})$,
$d_{1}:\mathcal{M}\times\partial\mathcal{M}\to\mathbb{R}$ be defined as
$d_{1}({x},{y})=v_{1}({x})-v_{1}({y})-({x}-{y})\cdot\nabla_{\mathcal{M}}^{{y}}v_{1}({y}),$
then we have

$\int_{\partial\mathcal{M}}\int_{\mathcal{M}}|d_{1}({x},{y})|^{2}\bar{R}_{\delta}({x},{y})d\mu_{{x}}d\tau_{{y}}\leq
C\delta^{3}\left\|v_{1}\right\|^{2}_{H^{2}(\mathcal{M})};$ (A.39)

  2. 2.

let the function $v_{2}\in H^{2}(\partial\mathcal{M})$,
$d_{2}:\partial\mathcal{M}\times\partial\mathcal{M}\to\mathbb{R}$ be defined
as
$d_{2}({x},{y})=v_{2}({x})-v_{2}({y})-({x}-{y})\cdot\nabla_{\partial\mathcal{M}}^{{y}}v_{2}({y}),$
then we have

$\int_{\partial\mathcal{M}}\int_{\partial\mathcal{M}}|d_{2}({x},{y})|^{2}{R}_{\delta}({x},{y})d\tau_{{x}}d\tau_{{y}}\leq
C\delta^{3}\left\|v_{2}\right\|^{2}_{H^{2}(\partial\mathcal{M})}.$ (A.40)

In both (A.39) and (A.40), ${R}_{\delta}({x},{y})$ can be replaced by the
other kernel functions $\bar{R}_{\delta}({x},{y}),\
\underline{R}_{\delta}({x},{y})$.

######  Proof A.4 (Proof of Lemma A.3).

We first study (A.40). For compact $\partial\mathcal{M}$, there exists a
$2\delta$-net,
$\mathcal{N}_{\delta}=\\{\bm{q}_{k}\in\partial\mathcal{M},k=1,2,...,N\\}$,
such that
$\partial\mathcal{M}\subset\bigcup_{k=1}^{N}\mathcal{B}_{\bm{q}_{k}}^{2\delta}$.

If in addition, the minimum reach of $\partial\mathcal{M}$ is no less than
$100\delta$, then according to page 8-9 of [38], for each integer $k\in[1,N]$,
there exists a $C^{2}$ smooth parametrization
$\bm{\phi}_{k}:\mathbb{R}^{m-1}\to\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{4\delta}$
in the neighborhood of $\bm{q}_{k}$ with width $4\delta$, with the following
local-small-deformation property:

$\frac{1}{2}|\theta_{1}-\theta_{2}|\leq|\bm{\phi}_{k}(\theta_{1})-\bm{\phi}_{k}(\theta_{2})|\leq
2|\theta_{1}-\theta_{2}|,\qquad\forall\
\theta_{1},\theta_{2}\in\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{4\delta}.$
(A.41)

For each fixed $k$ and any
${y}\in\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{2\delta},\
{x}\in\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{4\delta}$, we denote
$\bm{\alpha}=\bm{\phi}_{k}^{-1}({y}),\ \bm{\beta}=\bm{\phi}_{k}^{-1}({x}),\
\bm{\xi}=\bm{\beta}-\bm{\alpha}$. Using Newton-Leibniz formula, we get

$\begin{split}&d_{2}({x},{y})=v_{2}({x})-v_{2}({y})-({x}-{y})\cdot\nabla_{\partial\mathcal{M}}^{{y}}v_{2}({y})\\\
&=\int_{0}^{1}\nabla_{\partial\mathcal{M}}^{j}v_{2}(\bm{\phi}_{k}(\bm{\alpha}+t_{1}\bm{\xi}))\partial_{i}\phi_{k}^{j}(\bm{\alpha}+t_{1}\bm{\xi})\xi^{i}dt_{1}-\int_{0}^{1}\nabla_{\partial\mathcal{M}}^{j}v_{2}(\bm{\phi}_{k}(\bm{\alpha}))\partial_{i}\phi_{k}^{j}(\bm{\alpha}+t_{1}\bm{\xi})\xi^{i}dt_{1}\\\
&=\xi^{i}\xi^{i^{\prime}}\int_{0}^{1}\int_{0}^{1}\nabla_{\partial\mathcal{M}}^{j}\nabla_{\partial\mathcal{M}}^{j^{\prime}}v_{2}(\bm{\phi}_{k}(\bm{\alpha}+t_{1}t_{2}\bm{\xi}))t_{1}\partial_{i^{\prime}}\phi_{k}^{j^{\prime}}(\bm{\alpha}+t_{1}t_{2}\bm{\xi})\partial_{i}\phi_{k}^{j}(\bm{\alpha}+t_{1}\bm{\xi})dt_{2}dt_{1}\\\
&\leq
C\delta^{2}\int_{0}^{1}\int_{0}^{1}|\nabla_{\partial\mathcal{M}}\nabla_{\partial\mathcal{M}}v_{2}(\bm{\phi}_{k}(\bm{\alpha}+t_{1}t_{2}\bm{\xi}))|dt_{2}dt_{1}\\\
&\leq
C\delta^{2}(\int_{0}^{1}\int_{0}^{1}|\nabla_{\partial\mathcal{M}}\nabla_{\partial\mathcal{M}}v_{2}(\bm{\phi}_{k}(\bm{\alpha}+t_{1}t_{2}\bm{\xi}))|^{2}dt_{2}dt_{1})^{\frac{1}{2}},\end{split}$
(A.42)

this implies

$\begin{split}&\int_{\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{2\delta}}\int_{\partial\mathcal{M}}|d_{2}({x},{y})|^{2}{R}_{\delta}({x},{y})d\tau_{{x}}d\tau_{{y}}\\\
\leq&C\delta^{4}\int_{\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{2\delta}}\int_{\partial\mathcal{M}}\int_{0}^{1}\int_{0}^{1}|\nabla_{\partial\mathcal{M}}\nabla_{\partial\mathcal{M}}v_{2}(\bm{\phi}_{k}(\bm{\alpha}+t_{1}t_{2}\bm{\xi}))|^{2}dt_{2}dt_{1}{R}_{\delta}({x},{y})d\tau_{{x}}d\tau_{{y}}\\\
=&C\delta^{4}\int_{0}^{1}\int_{0}^{1}(\int_{\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{2\delta}}\int_{\partial\mathcal{M}\cap\mathcal{B}_{{y}}^{2\delta}}|\nabla_{\partial\mathcal{M}}\nabla_{\partial\mathcal{M}}v_{2}(\bm{\phi}_{k}(\bm{\alpha}+t_{1}t_{2}\bm{\xi}))|^{2}{R}_{\delta}({x},{y})d\tau_{{x}}d\tau_{{y}})dt_{2}dt_{1}\\\
\leq&C\delta^{4}\max\limits_{0\leq t\leq
1}\int_{\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{2\delta}}\int_{\partial\mathcal{M}\cap\mathcal{B}_{{y}}^{2\delta}}|\nabla_{\partial\mathcal{M}}\nabla_{\partial\mathcal{M}}v_{2}(\bm{\phi}_{k}(\bm{\alpha}+t\bm{\xi}))|^{2}{R}_{\delta}({x},{y})d\tau_{{x}}d\tau_{{y}}\\\
\leq&C\delta^{4}\max\limits_{0\leq t\leq
1}\int_{\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{2\delta}}\int_{\partial\mathcal{M}\cap\mathcal{B}_{{y}}^{2\delta}}|\nabla_{\partial\mathcal{M}}\nabla_{\partial\mathcal{M}}v_{2}({z})|^{2}{R}_{\delta}({y}+\frac{1}{4t}({z}-{y}),{y})(Cd\tau_{{z}})d\tau_{{y}}\\\
\leq&C\delta^{4}\int_{\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{4\delta}}(\max\limits_{0\leq
t\leq
1}\int_{\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{2\delta}}{R}_{\delta}({y}+\frac{1}{4t}({z}-{y}),{y})d\tau_{{y}})\
|\nabla_{\partial\mathcal{M}}\nabla_{\partial\mathcal{M}}v_{2}({z})|^{2}d\tau_{{z}}\\\
\leq&C\delta^{3}\int_{\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{4\delta}}|\nabla_{\partial\mathcal{M}}\nabla_{\partial\mathcal{M}}v_{2}({z})|^{2}d\tau_{{z}}\leq
C\delta^{3}\left\|v_{2}\right\|^{2}_{H^{2}(\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{4\delta})};\end{split}$
(A.43)

here the third last inequality is from the change of variable
${z}=\bm{\phi}_{k}(\bm{\alpha}+t\bm{\xi})$ and the assumption (A.41) of
$\bm{\phi}_{k}$ that $|{z}-{y}|\leq 2t|\bm{\xi}|\leq 4t|{x}-{y}|$. We thereby
conclude

$\begin{split}\int_{\partial\mathcal{M}}&\int_{\partial\mathcal{M}}|d_{2}({x},{y})|^{2}{R}_{\delta}({x},{y})d\tau_{{x}}d\tau_{{y}}\leq\sum\limits_{k=1}^{N}\int_{\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{2\delta}}\int_{\partial\mathcal{M}}|d_{2}({x},{y})|^{2}{R}_{\delta}({x},{y})d\tau_{{x}}d\tau_{{y}}\\\
&\leq
C\delta^{3}\sum\limits_{k=1}^{N}\left\|v_{2}\right\|^{2}_{H^{2}(\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{4\delta})}\leq
C\delta^{3}\left\|v_{2}\right\|^{2}_{H^{2}(\partial\mathcal{M})}.\end{split}$
(A.44)

The proof of (A.39) is similar. We will utilize the same $2\delta$-net
$\mathcal{N}_{\delta}=\\{\bm{q}_{k}\in\partial\mathcal{M},k=1,2,...,N\\}$ such
that
$\partial\mathcal{M}\subset\bigcup_{k=1}^{N}\mathcal{B}_{\bm{q}_{k}}^{2\delta}$.

If the minimum reach of $\mathcal{M}$ is no less than $100\delta$, then
according to page 8-9 of [38], for each integer $k\in[1,N]$, there exists a
$C^{2}$ smooth parametrization
$\bm{\Phi}_{k}:\mathbb{R}^{m}\to\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{4\delta}$
in the neighborhood of $\bm{q}_{k}$ with width $4\delta$, with the following
local-small-deformation property:

$\frac{1}{2}|\theta_{1}-\theta_{2}|\leq|\bm{\Phi}_{k}(\theta_{1})-\bm{\Phi}_{k}(\theta_{2})|\leq
2|\theta_{1}-\theta_{2}|,\qquad\forall\
\theta_{1},\theta_{2}\in\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{4\delta}.$
(A.45)

For each fixed $k$ and any
${y}\in\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{2\delta},\
{x}\in\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{4\delta}$, we denote
$\bm{\alpha}=\bm{\Phi}_{k}^{-1}({y}),\ \bm{\beta}=\bm{\Phi}_{k}^{-1}({x}),\
\bm{\xi}=\bm{\beta}-\bm{\alpha}$. By the similar argument as (A.42) and
(A.43), we derive

$\begin{split}d_{1}({x},{y})\leq
C\delta^{2}(\int_{0}^{1}\int_{0}^{1}|\nabla_{\mathcal{M}}\nabla_{\mathcal{M}}v_{1}(\bm{\Phi}_{k}(\bm{\alpha}+t_{1}t_{2}\bm{\xi}))|^{2}dt_{2}dt_{1})^{\frac{1}{2}},\end{split}$
(A.46)
$\begin{split}&\int_{\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{2\delta}}\int_{\mathcal{M}}|d_{1}({x},{y})|^{2}{R}_{\delta}({x},{y})d\mu_{{x}}d\tau_{{y}}\\\
\leq&C\delta^{4}\max\limits_{0\leq t\leq
1}\int_{\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{2\delta}}\int_{\mathcal{M}\cap\mathcal{B}_{{y}}^{2\delta}}|\nabla_{\mathcal{M}}\nabla_{\mathcal{M}}v_{1}(\bm{\Phi}_{k}(\bm{\alpha}+t\bm{\xi}))|^{2}{R}_{\delta}({x},{y})d\mu_{{x}}d\tau_{{y}}\\\
\leq&C\delta^{4}\max\limits_{0\leq t\leq
1}\int_{\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{2\delta}}\int_{\mathcal{M}\cap\mathcal{B}_{{y}}^{2\delta}}|\nabla_{\mathcal{M}}\nabla_{\mathcal{M}}v_{1}({z})|^{2}{R}_{\delta}({y}+\frac{1}{4t}({z}-{y}),{y})(Cd\mu_{{z}})d\tau_{{y}}\\\
\leq&C\delta^{4}\int_{\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{4\delta}}(\max\limits_{0\leq
t\leq
1}\int_{\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{2\delta}}{R}_{\delta}({y}+\frac{1}{4t}({z}-{y}),{y})d\tau_{{y}})\
|\nabla_{\mathcal{M}}\nabla_{\mathcal{M}}v_{1}({z})|^{2}d\mu_{{z}}\\\
\leq&C\delta^{3}\int_{\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{4\delta}}|\nabla_{\mathcal{M}}\nabla_{\mathcal{M}}v_{1}({z})|^{2}d\mu_{{z}}\leq
C\delta^{3}\left\|v_{1}\right\|^{2}_{H^{2}(\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{4\delta})};\end{split}$
(A.47)

this implies

$\begin{split}&\int_{\partial\mathcal{M}}\int_{\mathcal{M}}|d_{1}({x},{y})|^{2}{R}_{\delta}({x},{y})d\mu_{{x}}d\tau_{{y}}\leq\sum\limits_{k=1}^{N}\int_{\partial\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{2\delta}}\int_{\mathcal{M}}|d_{1}({x},{y})|^{2}{R}_{\delta}({x},{y})d\mu_{{x}}d\tau_{{y}}\\\
&\leq
C\delta^{3}\sum\limits_{k=1}^{N}\left\|v_{1}\right\|^{2}_{H^{2}(\mathcal{M}\cap\mathcal{B}_{\bm{q}_{k}}^{4\delta})}\leq
C\delta^{3}\left\|v_{1}\right\|^{2}_{H^{2}(\mathcal{M})}.\end{split}$ (A.48)

We then complete the proof of Lemma A.3.

No we move back to (A.38). In fact, by utilizing lemma A.3, we obtain

$\begin{split}&\int_{\partial\mathcal{M}}(\int_{\mathcal{M}}\delta|v_{1}({x})-v_{1}({y})|\tilde{R}_{\delta}({x},{y})d\mu_{{x}})^{2}d\tau_{{y}}\\\
\leq&C\delta^{2}\int_{\partial\mathcal{M}}(\int_{\mathcal{M}}|d_{1}({x},{y})|\tilde{R}_{\delta}({x},{y})d\mu_{{x}})^{2}d\tau_{{y}}+C\delta^{4}\int_{\partial\mathcal{M}}(\int_{\mathcal{M}}|\nabla_{\mathcal{M}}v_{1}({y})|\tilde{R}_{\delta}({x},{y})d\mu_{{x}})^{2}d\tau_{{y}}\\\
\leq&C\delta^{2}\int_{\partial\mathcal{M}}(\int_{\mathcal{M}}|d_{1}({x},{y})|^{2}\tilde{R}_{\delta}({x},{y})d\mu_{{x}}\int_{\mathcal{M}}\tilde{R}_{\delta}({x},{y})d\mu_{{x}})d\tau_{{y}}\\\
&+C\delta^{4}\int_{\partial\mathcal{M}}(\int_{\mathcal{M}}|\nabla_{\mathcal{M}}v_{1}({y})|^{2}\tilde{R}_{\delta}({x},{y})d\mu_{{x}}\int_{\mathcal{M}}\tilde{R}_{\delta}({x},{y})d\mu_{{x}})d\tau_{{y}}\\\
\leq&C\delta^{2}\int_{\partial\mathcal{M}}\int_{\mathcal{M}}|d_{1}({x},{y})|^{2}\tilde{R}_{\delta}({x},{y})d\mu_{{x}}d\tau_{{y}}+C\delta^{4}\int_{\partial\mathcal{M}}|\nabla_{\mathcal{M}}v_{1}({y})|^{2}d\tau_{{y}}\\\
\leq&C\delta^{6}\left\|v_{1}\right\|_{H^{2}(\mathcal{M})}^{2}+C\delta^{4}\left\|v_{1}\right\|_{H^{1}(\partial\mathcal{M})}^{2}\leq
C\delta^{4}\left\|v_{1}\right\|_{H^{2}(\mathcal{M})}^{2}\leq
C\delta^{4}\left\|u\right\|_{H^{4}(\mathcal{M})}^{2};\end{split}$ (A.49)

and similarly

$\begin{split}&\int_{\partial\mathcal{M}}(\int_{\partial\mathcal{M}}\delta^{2}|v_{2}({x})-v_{2}({y})|\tilde{R}_{\delta}({x},{y})d\tau_{{x}})^{2}d\tau_{{y}}\\\
\leq&C\delta^{4}\int_{\partial\mathcal{M}}(\int_{\partial\mathcal{M}}|d_{2}({x},{y})|\tilde{R}_{\delta}({x},{y})d\tau_{{x}})^{2}d\tau_{{y}}+C\delta^{6}\int_{\partial\mathcal{M}}(\int_{\partial\mathcal{M}}|\nabla_{\partial\mathcal{M}}v_{2}({y})|\tilde{R}_{\delta}({x},{y})d\tau_{{x}})^{2}d\tau_{{y}}\\\
\leq&C\delta^{4}\int_{\partial\mathcal{M}}(\int_{\partial\mathcal{M}}|d_{2}({x},{y})|^{2}\tilde{R}_{\delta}({x},{y})d\mu_{{x}}\int_{\partial\mathcal{M}}\tilde{R}_{\delta}({x},{y})d\mu_{{x}})d\tau_{{y}}\\\
&+C\delta^{6}\int_{\partial\mathcal{M}}(\int_{\partial\mathcal{M}}|\nabla_{\partial\mathcal{M}}v_{2}({y})|^{2}\tilde{R}_{\delta}({x},{y})d\mu_{{x}}\int_{\partial\mathcal{M}}\tilde{R}_{\delta}({x},{y})d\mu_{{x}})d\tau_{{y}}\\\
\leq&C\delta^{3}\int_{\partial\mathcal{M}}\int_{\partial\mathcal{M}}|d_{2}({x},{y})|^{2}\tilde{R}_{\delta}({x},{y})d\tau_{{x}}d\tau_{{y}}+C\delta^{4}\int_{\partial\mathcal{M}}|\nabla_{\partial\mathcal{M}}v_{2}({y})|^{2}d\tau_{{y}}\\\
\leq&C\delta^{6}\left\|v_{2}\right\|_{H^{2}(\partial\mathcal{M})}^{2}+C\delta^{4}\left\|v_{2}\right\|_{H^{1}(\partial\mathcal{M})}^{2}\leq
C\delta^{4}\left\|v_{2}\right\|_{H^{2}(\partial\mathcal{M})}^{2}\leq
C\delta^{4}\left\|u\right\|_{H^{4}(\mathcal{M})}^{2}.\end{split}$ (A.50)

Therefore we have complete the $L_{{y}}^{2}(\partial\mathcal{M})$ control of
the right hand side of (A.38) through (A.49) and (A.50). We then apply (A.38)
to conclude that

$\left\|\nabla_{\partial\mathcal{M}}s\right\|_{L^{2}(\partial\mathcal{M})}+\left\|\Delta_{\partial\mathcal{M}}s\right\|_{L^{2}(\partial\mathcal{M})}\leq
C\delta^{2}\left\|u\right\|_{H^{4}(\mathcal{M})}.$ (A.51)

## References

  * [1] B. Alali and M. Gunzburger.  Peridynamics and material interfaces.  Journal of Elasticity, 120, 2010.
  * [2] F. Andreu, J. M. Mazon, J. D. Rossi, and J. Toledo.  Nonlocal diffusion problems.  Math. Surveys Monogr. 165, AMS, Providence, RI, 2010.
  * [3] E. Askari, F. Bobaru, R. Lehoucq, M. Parks, S. Silling, and O. Weckner.  Peridynamics for multiscale materials modeling.  J. Physics.: Conf. Ser, 125, 2008.
  * [4] G. Barles, E. Chasseigne, C. Georgelin, and E. Jakobsen.  On neumann type problems for nonlocal equations set in a half space.  Trans. Am. Math. Soc., 366, 2014a.
  * [5] G. Barles, C. Georgelin, and E. Jakobsen.  On neumann and oblique derivatives boundary conditions for nonlocal elliptic equations.  J. Differ. Equ, 256, 2014b.
  * [6] R. Barreira, C. Elliott, and A. Madzvamuse.  The surface finite element method for pattern formation on evolving biological surfaces.  J. Math. Biol., 63:1095–1119, 2011.
  * [7] M. Belkin and P. Niyogi.  Laplacian eigenmaps for dimensionality reduction and data representation.  Neural Computation, 15(6):1373–1396, 2003.
  * [8] F. Bobaru, M. Yang, L. Alves, S. Silling, E. Askari, and J. Xu.  Convergence, adaptive refinement, and scaling in 1d peridynamics.  Int. J. Numer. Methods Eng., 77, 2009.
  * [9] J. W. Cahn, P. Fife, and O. Penrose.  A phase-field model for diffusion-induced grain-boundary motion.  Ann. Statist., 36(2):555–586, 2008.
  * [10] P. T. Choi, K. C. Lam, and L. M. Lui.  Flash: fast landmark aligned spherical harmonic parameterization for genus-0 closed brain surfaces.  SIAM Journal on Imaging Sciences, 8:67–94, 2015.
  * [11] R. R. Coifman, S. Lafon, A. B. Lee, M. Maggioni, F. Warner, and S. Zucker.  Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps.  In Proceedings of the National Academy of Sciences, pages 7426–7431, 2005.
  * [12] C. Cortazar, M. Elgueta, J. Rossi, and N. Wolanski.  How to approximate the heat equation with neumann boundary conditions by nonlocal diffusion problems.  Archive Ration. Mech. Anal., 187, 2008.
  * [13] K. Dayal and K. Bhattacharya.  Kinetics of phase transformations in the peridynamic formulation of continuum mechanics.  J. Mech. Pays. Solids, 54, 2006.
  * [14] S. Dipierro, X. Ros-Oton, and E. Valdinoci.  Nonlocal problems with neumann boundary conditions.  Rev. Mat. Iberoam, 2017.
  * [15] Q. Du, M. Gunzburger, R. B. Lehoucq, and K. Zhou.  Analysis and approximation of nonlocal diffusion problems with volume constraints.  SIAM Review, 54:667–696, 2012.
  * [16] Q. Du and Z. Shi.  A nonlocal stokes system with volume constraints.  Numerical Mathematics: Theory, Methods and Applications, 15, 2022\.
  * [17] Q. Du, X. Tian, and Z. Zhou.  Nonlocal diffusion models with consistent local and fractional limits.  arXiv:2203.00167v3, 2022.
  * [18] Q. Du, H. Xie, and X. Yin.  On the convergence to local limit of nonlocal models with approximated interaction neighborhoods.  SIAM Journal of Numerical Analysis, 60(4), 2022.
  * [19] C. Eilks and C. M. Elliott.  Numerical simulation of dealloying by surface dissolution via the evolving surface finite element method.  J. Comput. Phys., 227:9727–9741, 2008.
  * [20] C. M. Elliott and B. Stinner.  Modeling and computation of two phase geometric biomembranes using surface finite elements.  J. Comput. Phys., 229:6585–6612, 2010.
  * [21] S. Ganesan and L. Tobiska.  A coupled arbitrary lagrangian eulerian and lagrangian method for computation of free-surface flows with insoluble surfactants.  J. Comput. Phys., 228:2859–2873, 2009.
  * [22] X. Gu, Y. Wang, T. F. Chan, P. M. Thompson, and S.-T. Yau.  Genus zero surface conformal mapping and its application to brain surface mapping.  IEEE TMI, 23:949–958, 2004.
  * [23] A. J. James and J. Lowengrub.  A surfactant-conserving volume-of-fluid method for interfacial flows with insoluble surfactant.  J. Comput. Phys., 201:685–722, 2004.
  * [24] C.-Y. Kao, R. Lai, and B. Osting.  Maximization of laplace-beltrami eigenvalues on closed riemannian surfaces.  ESAIM: Control, Optimization and Calculus of Variations, 23:685–720, 2017.
  * [25] R. Lai, Z. Wen, W. Yin, X. Gu, and L. Lui.  Folding-free global conformal mapping for genus-0 surfaces by harmonic energy minimization.  Journal of Scientific Computing, 58:705–725, 2014.
  * [26] R. Lai and H. Zhao.  Multi-scale non-rigid point cloud registration using robust sliced-wasserstein distance via laplace-beltrami eigenmap.  SIAM Journal on Imaging Sciences, 10:449–483, 2017.
  * [27] H. Lee and Q. Du.  Second order accurate dirichlet boundary conditions for linear nonlocal diffusion problems.  Commun. Math. Sci., 20, 2022.
  * [28] Z. Li, Z. Shi, and J. Sun.  Point integral method for solving poisson-type equations on manifolds from point clouds with convergence guarantees.  Communications in Computational Physics, 22(1):228–258, 2017.
  * [29] T. W. Meng, P. T. Choi, and L. M. Lui.  Tempo: feature-endowed teichmuller extremal mappings of point clouds.  SIAM Journal on Imaging Sciences, 9:1582–1618, 2016.
  * [30] T. Mengesha and Q. Du.  Characterization of function space of vector fields and an application in nonlinear peridynamics.  Nonlinear Anal.:Theory Methods Appl., 140, 2016.
  * [31] M. P. Neilson, J. A. Mackenzie, S. D. Webb, and R. H. Insall.  Modelling cell movement and chemotaxis using pseudopod-based feedback.  SIAM J. Sci. Comput., 33:1035–1057, 2011.
  * [32] S. Osher, Z. Shi, and W. Zhu.  Low dimensional manifold model for image processing.  SIAM Journal on Imaging Sciences, 10(4), 2017.
  * [33] E. Oterkus and E. Madenci.  Peridynamic analysis of fiber-reinforced composed materials.  J. Mech. Mater. Struct., 7, 2012.
  * [34] G. Peyré.  Manifold models for signals and images.  Computer Vision and Image Understanding, 113:248–260, 2009.
  * [35] M. Reuter, F.-E. Wolter, and N. Peinecke.  Laplace-beltrami spectra as ”shape-dna” of surfaces and solids.  Computer-Aided Design, 38(4):342–366, 2006.
  * [36] Z. Shi.  Enforce the dirichlet boundary condition by volume constraint in point integral method.  Commun. Math. Sci, 15(6), 2017.
  * [37] Z. Shi.  Nonlocal approximation of elliptic operators with anisotropic coefficients on manifold.  Commun. Math. Sci., 17(3), 2019.
  * [38] Z. Shi and J. Sun.  Convergence of the point integral method for poisson equation on point cloud.  Research in the Mathematical Sciences, 4(1), 2017.
  * [39] Z. Shi and B. Wang.  Convergence of the weighted nonlocal laplacian on random point cloud.  J. Comput. Phys., 39(6), 2021.
  * [40] S. Silling, O. Weckner, E. Askari, and F. Bobaru.  Crack nucleation in a peridnamic solid.  Int. J. Fract., 162, 2010.
  * [41] Y. Tao, X. Tian, and Q. Du.  Nonlocal diffusion and peridynamic models with neumann type constraints and their numerical approximations.  Applied Mathematics and Application, 305, 2017.
  * [42] M. Taylor and D. Steigmann.  A two-dimensional peridynamic model for thin plates.  Math. Mech. Solids, 20, 2015.
  * [43] T. Wang and Z. Shi.  A nonlocal diffusion model with h1 convergence for dirichlet boundary.  arXiv:2302.03441v1, 2023.
  * [44] T. W. Wong, L. M. Lui, X. Gu, P. Thompson, T. Chan, and S.-T. Yau.  Instrinic feature extraction and hippocampal surface registration using harmonic eigenmap.  Technical Report, UCLA CAM Report 11-65, 2011.
  * [45] H. You, X. Lu, N. Trask, and Y. Yu.  A neumann-type boundary condition for nonlocal problems.  Mathematical Models and Methods in Applied Sciences, 2018.
  * [46] Y. Zhang and Z. Shi.  A nonlocal model of elliptic equation with jump coefficients on manifold.  Commun. Math. Sci., 19(7), 2021.
  * [47] Y. Zhang and Z. Shi.  A second-order nonlocal approximation for manifold poisson model with dirichlet boundary.  Researches in the Mathematical Sciences, 10(36), 2023.
  * [48] Y. Zhang and Z. Shi.  Truncation error analysis for an accurate nonlocal approximation to manifold poisson models with dirichlet boundary.  To appear in Communications in Mathematical Sciences, 2024.
  * [49] K. Zhou and Q. Du.  Mathematical and numerical analysis of linear peridynamic models with nonlocal boundary conditions.  SIAM J. Numer. Anal., 48:1759–1780, 2010.