# Monotone two-scale methods for a class of integrodifferential operators and
applications

Juan Pablo Borthagaray  Instituto de Matematica y Estadística “Rafael
Laguardia”, Facultad de Ingeniería, Universidad de la República, Montevideo,
Uruguay<EMAIL_ADDRESS>https://sites.google.com/view/jpborthagaray/home ,  Ricardo H. Nochetto
Department of Mathematics and Institute for Physical Science and Technology,
University of Maryland, College Park, MD 20742, USA<EMAIL_ADDRESS>https://www.math.umd.edu/ rhn/ ,  Abner J. Salgado  Department of Mathematics,
University of Tennessee, Knoxville, TN 37996, USA<EMAIL_ADDRESS>https://math.utk.edu/people/abner-salgado/ and  Céline Torres  Department of
Mathematics, University of Maryland, College Park, MD 20742, USA.
<EMAIL_ADDRESS>
###### Abstract.

We develop a monotone, two-scale discretization for a class of
integrodifferential operators of order $2s$, $s\in(0,1)$. We apply it to
develop numerical schemes, and convergence rates, for linear and obstacle
problems governed by such operators. As applications of the monotonicity, we
provide error estimates for free boundaries and a convergent numerical scheme
for a concave fully nonlinear, nonlocal, problem.

###### Key words and phrases:

Integrodifferential operators; Monotonicity; Obstacle problems; Fully
nonlinear nonlocal equations; Rate of convergence.

###### 2020 Mathematics Subject Classification:

35S15, 35R09, 45K05, 65R20, 65N06, 65N12, 65N15, 35R35, 65K15, 49L25,

JPB has been supported in part by Fondo Clemente Estable grant 2022-172393.

RHN has been supported in part by NSF grant DMS-1908267.

AJS has been supported in part by NSF grant DMS-2111228.

CT has been supported in part by NSF grant DMS-1908267 and the Swiss National
Science Foundation grant number P500PT_203082

Part of this work was completed while RHN and AJS were in residence at the
Institute for Computational and Experimental Research in Mathematics (ICERM)
in Providence, RI, during the Numerical PDEs: Analysis, Algorithms, and Data
Challenges semester program. ICERM is supported by NSF grant DMS-1929284.

##  1\. Introduction

In recent times, nonlocal models have gained a lot of popularity in the pure
and applied sciences. The reason behind this boom is manifold. From the point
of view of applications, it is claimed that they are able to encode a wider
range of phenomena when compared to their local counterparts. In this regard,
for instance, one can refer to the creation of peridynamics [65, 66], nonlocal
diffusion reaction equations [70], fractional Cahn-Hilliard models [4, 3],
fractional porous media equations [31], the fractional Schrödinger equation
[68], fractional viscoelasticity [26], fractional Monge-Ampère [18] and many
more. The interested reader can refer to the many existing overviews [12, 48,
32] for further references and insight.

In our opinion, nonlocal models started to gain the attention of the
mathematical community after the seminal work [20], where the authors showed
that the so-called fractional Laplacian (in the whole space) can be realized
as a degenerate elliptic operator in one more dimension. With this, many of
the techniques that were used to analyze local problems could now be applied
to nonlocal ones. Later, purely nonlocal techniques for many problems were
also introduced; see, for instance, [19, 45].

Regarding the numerical treatment of nonlocal problems, some early attempts
could be found in [41, 42, 43, 34, 33]. However, references [56, 14, 13, 1]
deserve special attention. In [56] the extension technique of [20] was
exploited to develop a numerical method, and its analysis, for the so-called
_spectral_ fractional Laplacian. References [14, 13] for the spectral and
integral fractional Laplacian, respectively, develop a nonconforming
approximation on the basis of an integral representation, a quadrature
formula, and spatial discretization. Finally, closer to our discussion here is
[1] where the _integral fractional Laplacian_ is considered. On the basis of
weighted Hölder regularity results developed in [61], the authors of this work
provide weighted Sobolev regularity, and construct a direct discretization
over graded meshes that, in two dimensions, is optimal with respect to degrees
of freedom.

The error estimates of [1, 15] are in the energy and $L^{2}(\Omega)$–norms.
Pointwise error estimates for the integral fractional Laplacian were obtained
in [39], where a monotone two-scale method for this operator was developed.
The error estimates in this work were based on weighted Hölder regularity and
the construction of suitable barriers.

Two-scale methods, such as those developed in [58, 46, 52, 54] and [39, 40]
naturally inherit a discrete maximum principle from the continuous operator.
This is in general not true for finite element approximations. For instance,
the finite element approximation of the classical Laplacian possesses a
discrete maximum principle only under certain geometric constraints on the
mesh [29, Section III.20]. We remark that for nonlocal operators, the lack of
monotonicity for the finite element approximation is not expected to be fixed
by only a geometric constraint as for its local counterpart, especially when
$s\approx 0$, where the stiffness matrix approaches a standard mass matrix.

This brings us then to the main goals of the present work. We will consider an
integrodifferential operator of order $2s$, $s\in(0,1)$,

(1)
${\mathcal{L}}_{\eta}[w](x)=\mathrm{v.p.}\int_{\mathbb{R}^{d}}\frac{w(x)-w(y)}{|x-y|^{d+2s}}\eta\left(\frac{x-y}{|x-y|}\right)\mathop{}\\!\mathrm{d}y,$

where the so-called _kernel_ $\eta:{\mathbb{S}}^{d-1}\to{\mathbb{R}}$ is
assumed to verify several properties that will be specified below. Operators
of this form have been extensively studied in probability and finance, as they
represent the generator of a $2s$–stable Lévy processes, see e.g., [30, 7].
The simplest, and most important, example is the aforementioned integral
fractional Laplacian $(-\Delta)^{s}$, which is defined by

(2) $\eta(\theta)\equiv
C(d,s)=\frac{4^{s}s\Gamma(s+d/2)}{\pi^{d/2}\Gamma(1-s)}>0.$

In [39], the authors propose a monotone discretization for the fractional
Laplacian based on the splitting of the integral in the definition of the
operator into a singular and a tail parts; the former is approximated by a
scaled finite-difference Laplacian. This, in essence, amounts to truncating
the kernel in a neighborhood of the origin. Inspired by that work, for a
broader class of operators with kernel $\eta\in{\mathcal{C}}(\lambda,\Lambda)$
(cf. Definition 1), we develop a monotone, two-scale method for the
discretization of the operator (1) and show its consistency under weighted
Hölder regularity assumptions. Instead of truncation, we propose to
_regularize_ the kernel in (1) near the origin. From an analysis perspective,
this means that we are using a zero-order operator for the approximation of
the singular part of the integral. Our motivation is to obtain approximations
of (1) that are robust in the limit $s\uparrow 1$, which may require wide
stencils.

As applications of the proposed scheme, we then consider three applications of
increasing difficulty. First we obtain pointwise rates of convergence for a
linear problem where the operator is of the form (1). Next, we consider an
obstacle problem for (1) and obtain pointwise rates of convergence. These
pointwise rates of convergence, in turn, allow us to provide error estimates
for free boundaries. Finally, we develop a numerical scheme for a class of
concave fully nonlinear integrodifferential operators of order $2s$. Rates of
convergence for problems of this type, however, are at this stage beyond our
reach. We need to point out that the rates of convergence are obtained under
realistic regularity assumptions. Namely, up to the boundary $C^{s}$
regularity, and interior $C^{\beta}$ regularity, with $\beta\in(2s,4)$. In the
linear case, the interior regularity is solely dictated by the regularity of
the forcing function, whereas the obstacle problem exhibits a reduced
regularity of $\beta\leq 1+s$.

Our presentation is organized as follows. Notation and the functional
framework we shall operate in is presented in Section 2. There we also
introduce the class of integrodifferential operators we shall be interested
in, and their most elementary properties. The two-scale discretization of our
operators is introduced in Section 3. We describe the action of each one of
the scales, namely regularization and discretization, and study the
consistency of approximations. The first application of our two-scale
discretization is the content of Section 4, where we study a linear problem
and provide pointwise rates of convergence for its numerical scheme. These
rates of convergence are obtained under realistic regularity assumptions.
Next, in Section 5, we study an obstacle problem. We propose and analyze a
pointwise convergent scheme. The rates are, once again, obtained under
realistic regularity assumptions. We then continue, in Sections 5.2 and 5.3,
with the study of the approximation of free boundaries. We provide a rate of
convergence for discrete boundaries, both in the presence of regular and
singular points. As a final application of our scheme, Section 6 studies a
concave fully nonlinear integrodifferential equation. We propose a convergent
scheme, albeit without explicit rates of convergence. Finally, Appendix A
provides some intuition, justification, and consequences for our design
choices behind the regularization scale.

##  2\. Notation and preliminaries

Let us introduce some notation and terminology that will be used throughout
our discussion. For $A,B\in\mathbb{R}$ the relation $A\lesssim B$ means that,
for a nonessential constant $c$, we have $A\leq cB$. The value of this
constant may change in every occurrence. $A\gtrsim B$ means $B\lesssim A$. If
$A\lesssim B\lesssim A$ we abbreviate this by saying $A\approx B$. The Landau
symbols, big-O and little-o, respectively, are ${\mathcal{O}}$ and
${\mathfrak{o}}$.

For $r>0$ and $x\in{\mathbb{R}}^{d}$ we denote by $B_{r}(x)$ the (open)
Euclidean ball of radius $r$ centered at $x$. We set $B_{r}=B_{r}(0)$. The
unit sphere in $\mathbb{R}^{d}$ is ${\mathbb{S}}^{d-1}=\partial B_{1}$.

Throughout our discussion $\Omega\subset\mathbb{R}^{d}$, with
$d\in{\mathbb{N}}$, is a bounded Lipschitz domain which we assume satisfies an
exterior ball condition. We denote its boundary by $\partial\Omega$ and its
complement by $\Omega^{c}={\mathbb{R}}^{d}\setminus\overline{\Omega}$.

For $x\in\Omega$ we define $\delta(x)=\operatorname*{dist}(x,\partial\Omega)$
and if $x,y\in\Omega$, then we set
$\delta(x,y)=\min\\{\delta(x),\delta(y)\\}$.

###  2.1. Function spaces and their norms

We will adhere to standard notation for the Lebesgue spaces, and their norms,
when they are defined either over the whole space or some domain. Since we are
concerned with pointwise estimates, we must work with functions that are at
least continuous. The space of functions $w:\bar{\Omega}\to{\mathbb{R}}$ that
are continuous is denoted by $C(\bar{\Omega})$. We recall that this space
endowed with the norm

$\|w\|_{C(\bar{\Omega})}=\sup_{x\in\bar{\Omega}}|w(x)|,$

is a Banach space. We also need spaces of continuously differentiable
functions. For $k\in{\mathbb{N}}$ we set

$C^{k}(\bar{\Omega})=\left\\{w:\bar{\Omega}\to{\mathbb{R}}\ \middle|\
D^{\beta}w\in C(\bar{\Omega}),|\beta|\leq k\right\\}.$

We also set $C^{0}(\bar{\Omega})=C(\bar{\Omega})$. The norm on these spaces is

$\|w\|_{C^{k}(\bar{\Omega})}=\sup_{|\beta|\leq
k}\|D^{\beta}w\|_{C(\bar{\Omega})}.$

To provide a more refined characterization of smoothness, for
$k\in{\mathbb{N}}_{0}$, and $\alpha\in(0,1]$ we define the Hölder spaces via

$C^{k,\alpha}(\bar{\Omega})=\left\\{w\in C^{k}(\bar{\Omega})\ \middle|\
[D^{\beta}w]_{C^{0,\alpha}(\bar{\Omega})}<\infty,|\beta|=k\right\\},$

where

$[w]_{C^{0,\alpha}(\bar{\Omega})}=\sup_{x,y\in\bar{\Omega}:x\neq
y}\frac{|w(x)-w(y)|}{|x-y|^{\alpha}}.$

The norm in the Hölder spaces is

$\|w\|_{C^{k,\alpha}(\bar{\Omega})}=\max\left\\{\|w\|_{C^{k}(\bar{\Omega})},\max_{|\beta|=k}[D^{\beta}w]_{C^{0,\alpha}(\bar{\Omega})}\right\\}.$

These spaces are complete.

We set $C^{k,0}(\bar{\Omega})=C^{k}(\bar{\Omega})$. In addition, whenever
$k\in{\mathbb{N}}_{0}$ and $\alpha\in(0,1]$ with
$k+\alpha\notin{\mathbb{N}}_{0}$, we may denote
$C^{k+\alpha}(\bar{\Omega})=C^{k,\alpha}(\bar{\Omega})$. Finally, for
$k\in{\mathbb{N}}_{0}$ and $\alpha\in(0,1]$ we will say that $w\in
C^{k,\alpha}(\Omega)$ if $w_{|\bar{U}}\in C^{k,\alpha}(\bar{U})$ for all
$U\Subset\Omega$. A further refinement of these spaces will be detailed when
needed.

We will also deal with fractional Sobolev spaces. For $r\in(0,1)$ we set

$H^{r}({\mathbb{R}}^{d})=\left\\{w\in L^{2}({\mathbb{R}}^{d})\ \middle|\
\|w\|_{H^{r}({\mathbb{R}}^{d})}<\infty\right\\},$

with

$\displaystyle\|w\|_{H^{r}({\mathbb{R}}^{d})}^{2}$
$\displaystyle=\|w\|_{L^{2}({\mathbb{R}}^{d})}^{2}+[w]_{H^{r}(\mathbb{R}^{d})}^{2},$
$\displaystyle[w]_{H^{r}(\mathbb{R}^{d})}^{2}$
$\displaystyle=\frac{C(d,r)}{2}\iint_{{\mathbb{R}}^{d}\times{\mathbb{R}}^{d}}\frac{|w(x)-w(y)|^{2}}{|x-y|^{d+2s}}\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}y,$

for which it is a Hilbert space.

We also set, again for $r\in(0,1)$,

${\widetilde{H}}^{r}(\Omega)=\left\\{w\in{L^{2}(\Omega)}\ \middle|\
{\widetilde{w}}\in H^{r}(\mathbb{R}^{d})\right\\},$

where by ${\widetilde{w}}$ we denote the extension by zero onto $\Omega^{c}$.
Owing to the fractional Poincaré inequality,

$\|w\|_{L^{2}(\Omega)}\lesssim[{\widetilde{w}}]_{H^{r}(\mathbb{R}^{d})},\quad\forall
w\in{\widetilde{H}}^{r}(\Omega),$

the seminorm $[\cdot]_{H^{r}(\mathbb{R}^{d})}$ is actually an equivalent norm
in ${\widetilde{H}}^{r}(\Omega)$ which makes it Hilbert. The dual of
${\widetilde{H}}^{r}(\Omega)$ is denoted by $H^{-r}(\Omega)$, and the duality
pairing will be $\langle\cdot,\cdot\rangle$. In what follows, if confusion
does not arise, we shall suppress the explicit mention of zero extensions.

###  2.2. The integrodifferential operator

We shall consider integrodifferential operators of the form (1). Regarding the
kernel, we encode its assumptions in the following definition.

######  Definition 1 (class ${\mathcal{C}}(\lambda,\Lambda)$).

Let $\lambda,\Lambda\in(0,\infty)$ with $\lambda\leq\Lambda$. We will say that
the kernel $\eta:{\mathbb{S}}^{d-1}\to{\mathbb{R}}$ belongs to the class
${\mathcal{C}}(\lambda,\Lambda)$ if it is:

  1. 1.

Symmetric, i.e., $\eta(\theta)=\eta(-\theta)$ for all
$\theta\in{\mathbb{S}}^{d-1}$.

  2. 2.

Elliptic, i.e., we have

$0<\lambda\leq\eta(\theta)\leq\Lambda,\quad\text{a.e.}\quad\theta\in{\mathbb{S}}^{d-1}.$

Owing to these assumptions, if $\eta\in{\mathcal{C}}(\lambda,\Lambda)$, a
fractional integration by parts shows that

$[w]_{H^{s}({\mathbb{R}}^{d})}^{2}\approx\langle{\mathcal{L}}_{\eta}[w],w\rangle,\quad\forall
w\in{\widetilde{H}}^{s}(\Omega),$

with equivalence constants that depend on $\eta$ only through the constants
$\lambda$ and $\Lambda$. In other words, for every
$\eta\in{\mathcal{C}}(\lambda,\Lambda)$ the expression

$\|w\|_{\eta,s}^{2}=\iint_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\frac{|w(x)-w(y)|^{2}}{|x-y|^{d+2s}}\eta\left(\frac{x-y}{|x-y|}\right)\mathop{}\\!\mathrm{d}y\mathop{}\\!\mathrm{d}x,$

is an equivalent Hilbertian norm on ${\widetilde{H}}^{s}(\Omega)$.

Owing to the positivity of the kernel, the operator (1) satisfies a comparison
principle.

######  Proposition 2 (comparison).

Let $K\colon\mathbb{R}^{d}\to\mathbb{R}$ be a positive function. Define the
operator ${\mathcal{L}}_{K}$ via

${\mathcal{L}}_{K}[w](x)=\mathrm{v.p.}\int_{\mathbb{R}^{d}}(w(x)-w(y))\,K(x-y)\mathop{}\\!\mathrm{d}y.$

Let $w\colon\mathbb{R}^{d}\to\mathbb{R}$ be such that, in the weak sense,
${\mathcal{L}}_{K}[w]\geq 0$ almost everywhere in $\Omega$ and $w\geq 0$ in
$\Omega^{c}$. Then, $w\geq 0$ in $\Omega$.

###### Proof.

See [60, Proposition 4.1]. ∎

As a final piece of preliminary notation we must introduce a set of
parameters, relations between them, and a space of functions with a prescribed
behavior away from the boundary. We let $\beta\in{\mathbb{R}}$, and assume it
satisfies

(3) $\beta>2s,\qquad\beta,\beta-2s\notin{\mathbb{N}}_{0}.$

Given such $\beta$ we define

(4) ${\mathcal{S}}^{\beta}(\Omega)=\left\\{w\in C^{\beta}(\Omega)\cap
C^{0,s}(\bar{\Omega})\ \middle|\ \|w\|_{C^{\beta}(\left\\{x\in\Omega\
\middle|\ \delta(x)\geq\rho\right\\})}\lesssim\rho^{s-\beta}\right\\}.$

The usefulness of this function class shall become clear once we perform the
error analysis. Indeed, as we shall show (cf. Lemma 11 below), solutions of
problems involving the operator ${\mathcal{L}}_{\eta}$ typically belong to
such class.

##  3\. Two-scale discretization

We now begin the discretization of our integrodifferential operator (1).
Besides consistency, a fundamental necessity, we wish to preserve its
comparison property detailed in Proposition 2.

These two requirements stand at odds of each other. For instance, when dealing
with (local) elliptic second order differential operators and their finite
element discretization, very stringent mesh requirements must be imposed to
retain a comparison principle, see [29, Section III.20] and [49, Section 3.5].
If, on the contrary, we discretize via finite differences, it is known that
wide stencils must be employed, see [49, Section 3.2]. To fulfill these two
conditions then, we will employ a two-scale discretization.

Two-scale discretizations have become a popular choice to develop schemes that
preserve the comparison principle for (local) elliptic second order
differential operators. As an example, the reader is referred to [58, 46, 52,
54]. Regarding the two-scale discretization of nonlocal problems we mention
[39, 40], where the authors proposed a two-scale discretization for the
fractional Laplacian.

###  3.1. The regularization scale

The first step in the approximation of (1) is regularization. The idea is
that, for some $\varepsilon>0$, we split the integral that defines the
operator in two parts: ${\mathbb{R}}^{d}\setminus B_{\varepsilon}$ and
$B_{\varepsilon}$. The integral in the small ball $B_{\varepsilon}$ is then
regularized by introducing a non singular kernel with suitable approximation
properties.

We begin by observing that, since the kernel $\eta$ is symmetric, we have

(5)
${\mathcal{L}}_{\eta}[w](x)=\frac{1}{2}\int_{{\mathbb{R}}^{d}}\left(2w(x)-w(x-y)-w(x+y)\right)\eta\left(\frac{y}{|y|}\right){\mathcal{K}}(|y|)\mathop{}\\!\mathrm{d}y,$

where ${\mathcal{K}}(r)=\tfrac{1}{r^{d+2s}}$. Let now $\varepsilon>0$. We
define

(6)
${\mathcal{L}}_{\eta,\varepsilon}[w](x)=\frac{1}{2}\int_{{\mathbb{R}}^{d}}\left(2w(x)-w(x-y)-w(x+y)\right)\eta\left(\frac{y}{|y|}\right){\mathcal{K}}_{\varepsilon}(|y|)\mathop{}\\!\mathrm{d}y,$

where the radial, nonsingular kernel ${\mathcal{K}}_{\varepsilon}$ is defined
as

(7)
${\mathcal{K}}_{\varepsilon}(r)=\begin{dcases}\frac{1}{r^{d+2s}},&r\geq\varepsilon,\\\
\frac{1}{\varepsilon^{d+2s}}+\gamma(r^{2}-\varepsilon^{2})+\nu(r^{3}-\varepsilon^{3}),&r<\varepsilon.\end{dcases}$

We point out that (7) gives

$\lim_{r\uparrow\varepsilon}{\mathcal{K}}_{\varepsilon}(\varepsilon)=\varepsilon^{-d-2s},\qquad\lim_{t\downarrow
0}\frac{{\mathcal{K}}_{\varepsilon}(t)-{\mathcal{K}}_{\varepsilon}(0)}{t}=0.$

The constants $\gamma$ and $\nu$ are chosen so that
${\mathcal{K}}_{\varepsilon}\in C^{1}([0,\infty))$ and

${\mathcal{L}}_{\eta}[q]={\mathcal{L}}_{\eta,\varepsilon}[q],\quad\forall
q\in{\mathbb{P}}_{2}.$

The smoothness requirement implies that

$2\gamma\varepsilon+3\nu\varepsilon^{2}=-(d+2s)\varepsilon^{-d-2s-1}.$

On the other hand, by the symmetry of $\eta$ we will have exactness provided
that

$\int_{0}^{\varepsilon}\frac{r^{2}}{r^{d+2s}}r^{d-1}\mathop{}\\!\mathrm{d}r=\int_{0}^{\varepsilon}r^{2}{\mathcal{K}}_{\varepsilon}(r)r^{d-1}\mathop{}\\!\mathrm{d}r.$

In short, the parameters $\gamma$ and $\nu$ have the values

(8) $\displaystyle\gamma$
$\displaystyle=-\frac{(4+d)(d+2s)(3+d+2s)}{4(1-s)}\varepsilon^{-d-2s-2}<0,$
$\displaystyle\nu$
$\displaystyle=\frac{(5+d)(d+2s)(2+d+2s)}{6(1-s)}\varepsilon^{-d-2s-3}>0.$

The following result studies the interior consistency of this regularization.
By interior we mean that we consider points $x\in\Omega$, such that
$B_{\varepsilon}(x)\subset\Omega$, i.e., all the points where the
regularization of the kernel takes place are contained in $\Omega$.

######  Theorem 3 (interior consistency).

Let $\Omega$ be a bounded Lipschitz domain that satisfies the exterior ball
condition. Assume that $\beta\in(2s,4]$ and
$\eta\in{\mathcal{C}}(\lambda,\Lambda)$. Let
$w\in{\mathcal{S}}^{\beta}(\Omega)$, where this class is defined in (4). Let
$\alpha_{0}>1$, $\varepsilon>0$, and $x\in\Omega$ be such that
$\delta(x)\geq\alpha_{0}\varepsilon$. Then, it holds that

$\left|{\mathcal{L}}_{\eta}[w](x)-{\mathcal{L}}_{\eta,\varepsilon}[w](x)\right|\lesssim\varepsilon^{\beta-2s}\delta(x)^{s-\beta}.$

The implicit constant depends on $d$, $s$, $\alpha_{0}$, $\beta$, $\lambda$,
and $\Lambda$.

###### Proof.

We need to estimate

$\left|{\mathcal{L}}_{\eta}[w](x)-{\mathcal{L}}_{\eta,\varepsilon}[w](x)\right|=\\\
\frac{1}{2}\left|\int_{B_{\varepsilon}}\left(2w(x)-w(x-y)-w(x+y)\right)\eta\left(\frac{y}{|y|}\right)\left(\frac{1}{|y|^{d+2s}}-{\mathcal{K}}_{\varepsilon}(|y|)\right)\mathop{}\\!\mathrm{d}y\right|\lesssim\\\
\int_{0}^{\varepsilon}\|w\|_{C^{\beta}(B_{r}(x))}r^{\beta}\left|\frac{1}{r^{d+2s}}-{\mathcal{K}}_{\varepsilon}(r)\right|r^{d-1}\mathop{}\\!\mathrm{d}r\\\
\lesssim\int_{0}^{\varepsilon}\left(\delta(x)-r\right)^{s-\beta}\left|\frac{1}{r^{d+2s}}-{\mathcal{K}}_{\varepsilon}(r)\right|r^{\beta+d-1}\mathop{}\\!\mathrm{d}r,$

where, in the last step, we used the interior Hölder estimate that is assumed
of $w$. Because $\delta(x)\geq\alpha_{0}\varepsilon$,

$\int_{0}^{\varepsilon}\left(\delta(x)-r\right)^{s-\beta}\left|\frac{1}{r^{d+2s}}-{\mathcal{K}}_{\varepsilon}(r)\right|r^{\beta+d-1}\mathop{}\\!\mathrm{d}r\lesssim\\\
\delta(x)^{s-\beta}\int_{0}^{\varepsilon}r^{\beta+d-1}\left(\frac{1}{r^{d+2s}}+{\mathcal{K}}_{\varepsilon}(r)\right)\mathop{}\\!\mathrm{d}r\lesssim\varepsilon^{\beta-2s}\delta(x)^{s-\beta}.\qed$

As it can be seen from the proof of Theorem 3, the interior consistency of our
regularization depends on how close we are to the boundary. In particular, we
need to have $B_{\varepsilon}(x)\Subset\Omega$. For this reason, given
$\varepsilon>0$ we let $\boldsymbol{\varepsilon}:\Omega\to(0,\infty)$ be
sufficiently smooth and satisfy

$\boldsymbol{\varepsilon}(x)\leq\min\left\\{\frac{\delta(x)}{2},\varepsilon\right\\},\quad\forall
x\in\Omega.$

With this function at hand, we then define

${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[w](x)={\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}(x)}[w](x).$

In order to properly leverage the boundary regularity of solutions, the choice
of $\boldsymbol{\varepsilon}$ is made more precise below; depending on the
problem under consideration.

Owing to the positivity of ${\mathcal{K}}_{\varepsilon}$, the comparison
principle from Proposition 2 holds for
${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}$. Thus, pointwise consistency
estimates for this operator can be achieved by combining a comparison
principle with a suitable barrier function. We thus construct a barrier.

######  Lemma 4 (barrier).

Let $\Omega$ be a bounded, Lipschitz domain. Define

$b(x)=\chi_{\Omega}(x),\quad\forall x\in\Omega.$

Then,

$\delta(x)^{-2s}\lesssim{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[b](x).$

###### Proof.

Since, by construction, we have $\boldsymbol{\varepsilon}(x)\leq\delta(x)$ for
all $x\in\Omega$ we may write

${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[b](x)=\int_{{\mathbb{R}}^{d}}\left(1-b(y)\right)\eta\left(\frac{x-y}{|x-y|}\right){\mathcal{K}}_{\boldsymbol{\varepsilon}(x)}(|x-y|)\mathop{}\\!\mathrm{d}y\geq\lambda\int_{\Omega^{c}}\frac{1}{|x-y|^{d+2s}}\mathop{}\\!\mathrm{d}y.$

Since

$\delta(x)^{-2s}\lesssim\int_{\Omega^{c}}\frac{1}{|x-y|^{d+2s}}\mathop{}\\!\mathrm{d}y,$

with a hidden constant that depends on $d,$ $s,$ and $\Omega$, the result
follows. ∎

We are now in position to estimate the consistency of our regularized
operator.

######  Theorem 5 (consistency).

In the setting of Theorem 3 assume, in addition, that
$w,w_{\varepsilon}\in{\mathcal{S}}^{\beta}(\Omega)$ verify

${\mathcal{L}}_{\eta}[w]={\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[w_{\varepsilon}]\quad\text{
in }\Omega,\qquad w=w_{\varepsilon}=0,\quad\text{ in }\Omega^{c}.$

Then we have

$\|w-w_{\varepsilon}\|_{L^{\infty}(\Omega)}\lesssim\max\\{\varepsilon^{s},\varepsilon^{\beta-2s}\\}.$

###### Proof.

Let $x\in\Omega$. Since
$\boldsymbol{\varepsilon}(x)\leq\tfrac{1}{2}\delta(x)$, arguing as in the
proof of Theorem 3 we get

$\displaystyle{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[w-w_{\varepsilon}](x)$
$\displaystyle={\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[w](x)-{\mathcal{L}}_{\eta}[w](x)\lesssim\|w\|_{C^{\beta}(B_{\boldsymbol{\varepsilon}(x)}(x))}\boldsymbol{\varepsilon}(x)^{\beta-2s}$
$\displaystyle\lesssim(\delta(x)-\boldsymbol{\varepsilon}(x))^{s-\beta}\boldsymbol{\varepsilon}(x)^{\beta-2s}\lesssim
2^{s-\beta}\delta(x)^{s-\beta}\boldsymbol{\varepsilon}(x)^{\beta-2s}$
$\displaystyle\lesssim\delta(x)^{3s-\beta}\boldsymbol{\varepsilon}(x)^{\beta-2s}\delta(x)^{-2s}.$

If $3s-\beta\geq 0$ continue the estimate as

${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[w-w_{\varepsilon}](x)\lesssim\varepsilon^{\beta-2s}\delta(x)^{-2s}\lesssim\varepsilon^{\beta-2s}{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[b](x),$

where we used the barrier function of Lemma 4. If, on the other hand,
$3s-\beta<0$ we use that

$\delta(x)^{3s-\beta}\boldsymbol{\varepsilon}(x)^{\beta-2s}=\left(\frac{\boldsymbol{\varepsilon}(x)}{\delta(x)}\right)^{\beta-3s}\boldsymbol{\varepsilon}(x)^{s}\lesssim\varepsilon^{s}.$

Gathering both cases we conclude that, for every $x\in\Omega$, we have

${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[w-w_{\varepsilon}](x)\lesssim\max\\{\varepsilon^{\beta-2s},\varepsilon^{s}\\}{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[b](x).$

An application of the comparison principle for the operator
${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}$ allows us to conclude. ∎

###  3.2. The discretization scale

We assume that $\Omega$ is a convex polytope and let ${\mathscr{T}}=\\{T\\}$
be a conforming and shape regular simplicial triangulation of $\Omega$. The
elements $T\in{\mathscr{T}}$ are assumed to be closed. We set, for
$T\in{\mathscr{T}}$, $h_{T}=\operatorname*{diam}(T)$. We denote by
$\mathscr{N}_{{\mathscr{T}}}$ the set of vertices of ${\mathscr{T}}$. The
interior and boundary vertices are, respectively,

$\mathscr{N}_{{\mathscr{T}}}^{i}=\mathscr{N}_{{\mathscr{T}}}\cap\Omega,\qquad\mathscr{N}_{{\mathscr{T}}}^{\partial}=\mathscr{N}_{{\mathscr{T}}}\cap\partial\Omega.$

For each interior vertex $z\in\mathscr{N}_{{\mathscr{T}}}^{i}$ we define its
patch to be

$\omega_{z}=\bigcup\left\\{T\in{\mathscr{T}}\ \middle|\ z\in T\right\\}.$

By $h_{z}$ we denote the radius of the ball of maximal radius, centered at
$z$, that can be inscribed in $\omega_{z}$. Over such a triangulation we
define the following spaces of functions

(9) $\displaystyle{\mathbb{V}}_{\mathscr{T}}$
$\displaystyle=\left\\{w_{\mathscr{T}}\in C(\bar{\Omega})\ \middle|\
w_{{\mathscr{T}}|T}\in{\mathbb{P}}_{1},\ \forall T\in{\mathscr{T}}\right\\},$
(10) $\displaystyle{\mathbb{V}}_{\mathscr{T}}^{0}$
$\displaystyle=\left\\{w_{\mathscr{T}}\in C(\bar{\Omega})\ \middle|\
w_{{\mathscr{T}}|\partial\Omega}=0\right\\},$

Notice that any function $w_{\mathscr{T}}\in{\mathbb{V}}_{\mathscr{T}}^{0}$
can be trivially extended to $\Omega^{c}$ by zero. When this causes no
confusion, we shall not make a distinction between a function and its
extension.

It is a general fact that solutions to problems involving integrodifferential
operators, like (1), exhibit an algebraically singular behavior near the
boundary, independently of the smoothness of the problem data. The problems
that we shall be interested in are no exception; see the regularity results of
Sections 4.1 and 5. To compensate this we will consider a mesh that is graded
towards the boundary as it was studied in [1]. We consider a mesh size $h>0$
and parameter $\mu\geq 1$. Our mesh ${\mathscr{T}}$ is assumed to satisfy

(11) $h_{T}\approx\begin{dcases}h^{\mu},&T\cap\partial\Omega\neq\emptyset,\\\
h\operatorname*{dist}(T,\partial\Omega)^{\tfrac{\mu-1}{\mu}}&T\cap\partial\Omega=\emptyset.\end{dcases}$

As shown in [17, Remark 4.14] we have that

(12)
$\dim{\mathbb{V}}_{\mathscr{T}}^{0}\approx\begin{dcases}h^{(1-d)\mu},&\mu\geq\frac{d}{d-1},\\\
h^{-d}|\log h|,&\mu=\frac{d}{d-1},\\\ h^{-d},&\mu<\frac{d}{d-1}.\end{dcases}$

Finally, we observe that, under the condition (11), we have

(13) $h_{z}\approx h\delta(z)^{1-1/\mu}\quad\forall
z\in\mathscr{N}_{{\mathscr{T}}}^{i}.$

At this point we impose that, at least,
$\frac{1}{2}h_{z}\leq\boldsymbol{\varepsilon}(z)$. We shall later refine this
choice; see Definition 12 and 20 for the linear and obstacle problems,
respectively.

####  3.2.1. Consistency of interpolation

Let $I_{h}:C(\bar{\Omega})\to{\mathbb{V}}_{\mathscr{T}}$ be the Lagrange
interpolation operator. For $z\in\mathscr{N}_{{\mathscr{T}}}^{i}$ we wish to
estimate the consistency error

${\mathcal{E}}[w,z]={\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[I_{h}w](z)-{\mathcal{L}}_{\eta}[w](z),$

provided the function $w$ possesses suitable, but realistic, smoothness. We
begin by rewriting this error as

${\mathcal{E}}[w,z]=\left({\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[w](z)-{\mathcal{L}}_{\eta}[w](z)\right)+{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[I_{h}w-w](z).$

The first term entails the regularization error, and it was estimated in
Theorem 3. Our immediate goal shall be to estimate the second term.

######  Lemma 6 (refined interpolation estimate).

Let $\beta>0$ satisfy (3), $\bar{\beta}=\min\\{\beta,2\\}$, and
$w\in{\mathcal{S}}^{\beta}(\Omega)$. Furthermore, assume the mesh
${\mathscr{T}}$ satisfies (11). Then, for all $T\in{\mathscr{T}}$, we have

$\|w-I_{h}w\|_{L^{\infty}(T)}\lesssim\begin{dcases}h^{\bar{\beta}}\operatorname*{dist}(T,\Omega)^{s-\bar{\beta}/\mu},&T\cap\partial\Omega=\emptyset,\\\
h^{\mu s},&T\cap\partial\Omega\neq\emptyset.\end{dcases}$

Consequently, we have the global interpolation estimate

(14) $\|w-I_{h}w\|_{L^{\infty}(\Omega)}\lesssim\max\\{h^{\mu
s},h^{\bar{\beta}}\\}.$

###### Proof.

If $w\in C^{\alpha}(T)$ for some $\alpha\in(0,2]$, we have

$\|w-I_{h}w\|_{L^{\infty}(T)}\lesssim\|w\|_{C^{\alpha}(T)}h_{T}^{\alpha}.$

Now, if $T\cap\partial\Omega=\emptyset$ the definition of the class
${\mathcal{S}}^{\beta}(\Omega)$, given in (4), and the mesh grading imply that

$\|w-I_{h}w\|_{L^{\infty}(T)}\lesssim\operatorname*{dist}(T,\partial\Omega)^{s-\bar{\beta}}h^{\bar{\beta}}\operatorname*{dist}(T,\partial\Omega)^{\bar{\beta}(\mu-1)/\mu}\lesssim
h^{\bar{\beta}}\operatorname*{dist}(T,\Omega)^{s-\bar{\beta}/\mu}.$

If, on the contrary, $T\cap\partial\Omega\neq\emptyset$, we estimate

$\|w-I_{h}w\|_{L^{\infty}(T)}\lesssim\|w\|_{C^{0,s}(T)}h_{T}^{s}\lesssim
h^{\mu s}.$

Estimate (14) follows from a closer inspection of the case
$T\cap\partial\Omega=\emptyset$. If $s\geq\bar{\beta}/\mu$ then there is
nothing to prove, so we assume $s<\bar{\beta}/\mu$. Let
$z\in\mathscr{N}_{{\mathscr{T}}}^{i}$ be a vertex of $T$. By shape regularity
and (13), we have $h_{T}\approx h_{z}\approx h\delta(z)^{1-1/\mu},$ and
$\operatorname*{dist}(T,\partial\Omega)\approx\delta(z)$. Therefore, we can
write

$\|w-I_{h}w\|_{L^{\infty}(T)}\lesssim h^{\mu s}h^{\bar{\beta}-\mu
s}\delta(z)^{s-\bar{\beta}/\mu}\approx h^{\mu
s}\left(\frac{h_{z}}{\delta(z)}\right)^{\bar{\beta}-\mu s}\lesssim h^{\mu
s}.\qed$

We can now estimate the remaining consistency term.

######  Proposition 7 (consistency of interpolation).

Let the function $\boldsymbol{\varepsilon}:\Omega\to{\mathbb{R}}$ be such that
$\tfrac{1}{2}h_{z}\leq\boldsymbol{\varepsilon}(z)\leq\tfrac{1}{2}\delta(z)$
for all $z\in\mathscr{N}_{{\mathscr{T}}}^{i}$. In the setting of Lemma 6 we
have, for all $z\in\mathscr{N}_{{\mathscr{T}}}^{i}$,

(15)
$\left|{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[I_{h}w-w](z)\right|\lesssim
h^{\bar{\beta}}\delta(z)^{s-\bar{\beta}/\mu}\boldsymbol{\varepsilon}(z)^{-2s}+\max\\{h^{\mu
s},h^{\bar{\beta}}\\}\delta(z)^{-2s}.$

###### Proof.

Let $z\in\mathscr{N}_{{\mathscr{T}}}^{i}$. We observe that $I_{h}w(z)=w(z)$
and that $I_{h}w\equiv w\equiv 0$ on $\Omega^{c}$. Therefore, we have

${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[I_{h}w-w](z)=\int_{\Omega}\left(w(y)-I_{h}w(y))\right){\mathcal{K}}_{\boldsymbol{\varepsilon}(z)}(|z-y|)\eta\left(\frac{z-y}{|z-y|}\right)\mathop{}\\!\mathrm{d}y$

We define

(16) $\Omega^{\partial}=\bigcup\left\\{T\in{\mathscr{T}}\ \middle|\
T\cap\partial\Omega\neq\emptyset\right\\},$

and partition the integration domain into

$\bar{\Omega}=\bigcup_{i=1}^{3}D_{i},$

where

$\displaystyle D_{1}$ $\displaystyle=\Omega^{\partial}\cap
B_{\delta(z)/2}(z)^{c},$ $\displaystyle D_{2}$
$\displaystyle=(\Omega\setminus\Omega^{\partial})\cap B_{\delta(z)/2}(z)^{c},$
$\displaystyle D_{3}$ $\displaystyle=B_{\delta(z)/2}(z),$

and estimate each term separately.

Estimate on $D_{1}$: Since every point $y\in D_{1}$ belongs to a boundary-
touching element, we use the first part of Lemma 6 to write

$\left|\int_{D_{1}}\left(w(y)-I_{h}w(y))\right){\mathcal{K}}_{\boldsymbol{\varepsilon}(z)}(|z-y|)\eta\left(\frac{z-y}{|z-y|}\right)\mathop{}\\!\mathrm{d}y\right|\lesssim
h^{\mu s}\int_{D_{1}}\frac{1}{|z-y|^{d+2s}}\mathop{}\\!\mathrm{d}y\lesssim\\\
h^{\mu
s}\int_{B_{\delta(z)/2}(z)^{c}}\frac{1}{|z-y|^{d+2s}}\mathop{}\\!\mathrm{d}y\lesssim
h^{\mu s}\delta(z)^{-2s}.$

Estimate on $D_{2}$: Notice now that every element $y\in D_{2}$ belongs to a
non-boundary-touching element. We can then invoke the other case in the first
part of Lemma 6 to estimate

$\left|\int_{D_{2}}\left(w(y)-I_{h}w(y))\right){\mathcal{K}}_{\boldsymbol{\varepsilon}(z)}(|z-y|)\eta\left(\frac{z-y}{|z-y|}\right)\mathop{}\\!\mathrm{d}y\right|\lesssim
h^{\bar{\beta}}\int_{D_{2}}\frac{\delta(y)^{s-\bar{\beta}/\mu}}{|z-y|^{d+2s}}\mathop{}\\!\mathrm{d}y.$

Now, if $s-\tfrac{\bar{\beta}}{\mu}\geq 0$ we simply estimate

$h^{\bar{\beta}}\int_{D_{2}}\frac{\delta(y)^{s-\bar{\beta}/\mu}}{|z-y|^{d+2s}}\mathop{}\\!\mathrm{d}y\leq
h^{\bar{\beta}}\int_{B_{\delta(z)/2}(z)^{c}}\frac{1}{|z-y|^{d+2s}}\mathop{}\\!\mathrm{d}y\lesssim
h^{\bar{\beta}}\delta(z)^{-2s}.$

If, instead, $s-\tfrac{\bar{\beta}}{\mu}<0$ then

$h^{\bar{\beta}}\int_{D_{2}}\frac{\delta(y)^{s-\bar{\beta}/\mu}}{|z-y|^{d+2s}}\mathop{}\\!\mathrm{d}y\lesssim
h^{\bar{\beta}}h^{\mu(s-\bar{\beta}/\mu)}\int_{D_{2}}\frac{1}{|z-y|^{d+2s}}\mathop{}\\!\mathrm{d}y\lesssim
h^{\mu s}\delta(z)^{-2s}.$

Estimate on $D_{3}$: Observe that $\delta(y)\geq\delta(z)-|z-y|$ and, if $y\in
D_{3}$, we also have that
$\tfrac{1}{2}\delta(z)\leq\delta(y)\leq\tfrac{3}{2}\delta(z)$. Now, since
$z\in\mathscr{N}_{{\mathscr{T}}}^{i}$ we have that $\delta(z)\gtrsim h^{\mu}$.
Consequently, if $y\in T\cap D_{3}$ for some $T$ such that
$T\cap\partial\Omega=\emptyset$,

$h_{T}\lesssim h\delta(y)^{1-1/\mu}\lesssim h\delta(z)^{1-1/\mu}.$

We, once again, invoke the estimates of Lemma 6 to obtain

$\left|\int_{D_{3}\cap(\Omega\setminus\Omega^{\partial})}\left(w(y)-I_{h}w(y))\right){\mathcal{K}}_{\boldsymbol{\varepsilon}(z)}(|z-y|)\eta\left(\frac{z-y}{|z-y|}\right)\mathop{}\\!\mathrm{d}y\right|\lesssim\\\
h^{\bar{\beta}}\delta(z)^{s-\bar{\beta}/\mu}\int_{D_{3}\cap(\Omega\setminus\Omega^{\partial})}{\mathcal{K}}_{\boldsymbol{\varepsilon}(z)}(|z-y|)\mathop{}\\!\mathrm{d}y\lesssim\\\
h^{\bar{\beta}}\delta(z)^{s-\bar{\beta}/\mu}\left(\int_{\boldsymbol{\varepsilon}(z)}^{\delta(z)/2}\frac{1}{r^{1+2s}}\mathop{}\\!\mathrm{d}r+\int_{B_{\boldsymbol{\varepsilon}(z)}(z)}{\mathcal{K}}_{\boldsymbol{\varepsilon}(z)}(|z-y|)\mathop{}\\!\mathrm{d}y\right)\lesssim\\\
h^{\bar{\beta}}\delta(z)^{s-\bar{\beta}/\mu}\boldsymbol{\varepsilon}(z)^{-2s}.$

On the other hand, if $D_{3}\cap\Omega^{\partial}\neq\emptyset$ it means that
$\operatorname*{dist}(B_{\delta(z)/2}(z),\partial\Omega)\lesssim h^{\mu}$, and
therefore
$\delta(z)=\operatorname*{dist}(z,\partial\Omega)\lesssim\delta(z)/2+h^{\mu}$,
namely $\delta(z)\approx h^{\mu}\approx h_{z}$. Therefore, it must be
$\boldsymbol{\varepsilon}(z)\approx\delta(z)$. In such a case, on
$D_{3}\cap\Omega^{\partial}$, we have

$\left|\int_{D_{3}\cap\Omega^{\partial}}\left(w(y)-I_{h}w(y))\right){\mathcal{K}}_{\boldsymbol{\varepsilon}(z)}(|z-y|)\eta\left(\frac{z-y}{|z-y|}\right)\mathop{}\\!\mathrm{d}y\right|\lesssim\\\
h^{\mu
s}\int_{D_{3}\cap\Omega^{\partial}}{\mathcal{K}}_{\boldsymbol{\varepsilon}(z)}(|z-y|)\mathop{}\\!\mathrm{d}y\lesssim\\\
h^{\mu
s}\left(\int_{\boldsymbol{\varepsilon}(z)}^{\delta(z)/2}\frac{1}{r^{1+2s}}\mathop{}\\!\mathrm{d}r+\int_{B_{\boldsymbol{\varepsilon}(z)}(z)}{\mathcal{K}}_{\boldsymbol{\varepsilon}(z)}(|z-y|)\mathop{}\\!\mathrm{d}y\right)\lesssim\\\
h^{\mu s}\boldsymbol{\varepsilon}(z)^{-2s}\approx h^{\mu s}\delta(z)^{-2s}.$

Gathering all the obtained estimates leads to the claim. ∎

##  4\. The linear problem

Having studied a consistent and monotone discretization of the operator
${\mathcal{L}}_{\eta}$, defined in (1), we proceed to use this discretization
to propose and analyze numerical methods for problems of increasing
complexity. The first one shall be a linear one.

We consider the following problem: Let $s\in(0,1)$ and
$\eta\in{\mathcal{C}}(\lambda,\Lambda)$ for some $\lambda,\Lambda>0$. Given
$f\in H^{-s}(\Omega)$, find $u\in{\widetilde{H}}^{s}(\Omega)$ such that

(17) ${\mathcal{L}}_{\eta}[u]=f,\qquad\text{ in }\Omega.$

Notice that, by virtue of the definition of the solution space, we are
implicitly providing the exterior condition $u=0$ on $\Omega^{c}$. Owing to
the fact that $\|\cdot\|_{\eta,s}$ is an equivalent norm on
${\widetilde{H}}^{s}(\Omega)$, existence and uniqueness of a weak solution
follows immediately. In addition, since $\eta$ is positive, a nonlocal maximum
principle holds; see Proposition 2.

###  4.1. Regularity

The regularity properties of $u$, solution of (17), are of utmost relevance
for its numerical approximation. In contrast to local elliptic operators, it
is well known that solutions to (17) possess limited regularity near the
boundary, regardless of the smoothness of $\Omega$ and $f$; see [69, 38, 16].
The following is an _optimal_ regularity result.

######  Lemma 8 (optimal Hölder regularity).

Let $s\in(0,1)$ and let $\Omega\subset\mathbb{R}^{d}$ be a bounded Lipschitz
domain that satisfies the exterior ball condition. Let $f\in
L^{\infty}(\Omega)$ and $u\in{\widetilde{H}}^{s}(\Omega)$ be the (weak)
solution to (17). Then $u\in C^{0,s}(\bar{\Omega})$. Moreover, we have

$\|u\|_{C^{0,s}(\bar{\Omega})}\leq C\|f\|_{L^{\infty}(\Omega)},$

where the implicit constant depends only on $d,\Omega$ and $s$.

###### Proof.

See [63, Proposition 4.5]. ∎

The previous regularity is optimal, as exemplified by existing explicit
solutions to the fractional Laplace problem over a ball (cf. [37]). The reason
for this limited regularity gain lies in the fact that there is an algebraic
boundary singularity present in the solution, which can be characterized via

$u(x)\approx\operatorname*{dist}(x,\partial\Omega)^{s},$

as $x$ approaches $\partial\Omega$. Higher order regularity estimates can be
obtained if one takes into account such boundary behavior. In [62] such
results were obtained for the fractional Laplacian over weighted Hölder
spaces, where the weights are given by powers of the distance to the boundary.

######  Definition 9 (weighted Hölder space).

Let $\sigma>-1$ and $\beta=k+\gamma>0$ with
$k\in\mathbb{N}_{0},\gamma\in(0,1]$. For $w\in C^{k}(\Omega)$ define the
seminorm

$|w|_{\beta,\Omega}^{(\sigma)}=\sup_{x,y\in\Omega:x\neq
y}\left(\delta(x,y)^{\beta+\sigma}\frac{\left|D^{k}w(x)-D^{k}w(y)\right|}{|x-y|^{\gamma}}\right),$

and the norm

$\|w\|_{\beta,\Omega}^{(\sigma)}=\sum_{\ell=1}^{k}\sup_{x\in\Omega}\left(\delta(x)^{\ell+\sigma}|D^{\ell}w(x)|\right)+|w|_{\beta,\Omega}^{(\sigma)}+\begin{dcases}\sup_{x\in\Omega}\left(\delta(x)^{\sigma}|v(x)|\right),&\sigma\geq
0,\\\ \|v\|_{C^{0,-\sigma}(\bar{\Omega})},&\sigma<0.\end{dcases}$

The methods used in [61] can be adapted to operators of the form (1) to obtain
the following weighted Hölder regularity.

######  Lemma 10 (weighted Hölder estimates).

Let $\Omega$ be a bounded domain that satisfies the exterior ball condition,
and let $\beta>0$ satisfy (3). Moreover, we have

$\|u\|_{\beta,\Omega}^{(-s)}\lesssim\|u\|_{C^{0,s}(\mathbb{R}^{d})}+\|f\|_{\beta-2s,\Omega}^{(s)},$

where the implicit constant only dependsLet $f\in C^{\beta-2s}(\Omega)$ be
such that $\|f\|_{\beta-2s,\Omega}^{(s)}<\infty$. If, for $k\in{\mathbb{N}}$,
$\beta\geq k$, then we additionally assume that $\eta\in
C^{k}({\mathbb{S}}^{d-1})$. In this setting, the solution of (17) satisfies
$u\in C^{\beta}(\Omega)$. on $d$, $\Omega$, $\lambda$, $\Lambda$ and $s$.

###### Proof.

The result can be obtained by using the same procedure as in [61, Proposition
1.4] for $\eta\equiv 1$. For completeness, we repeat the main parts of the
argument.

Let $x_{0}\in\Omega$ and $R=\frac{\delta(x_{0})}{K}$ for some
$K\in\mathbb{N}$. We define $\tilde{u}(x)=u(x_{0}+Rx)-u(x_{0})$ and, following
[61, Proposition 1.4], it is possible to show that

(18) $\displaystyle\|\tilde{u}\|_{C^{0,s}(B_{1})}$ $\displaystyle\leq
CR^{s}[u]_{C^{0,s}(B_{R}(x_{0}))}$ (19)
$\displaystyle\|(1+|\cdot|)^{-d-2s}\tilde{u}(\cdot)\|_{L^{1}(\mathbb{R}^{n})}$
$\displaystyle\leq CR^{s}[u]_{C^{0,s}(\mathbb{R}^{d})}.$

In addition, if $\beta\leq 1$,

$\displaystyle\sup_{x,y\in
B_{1}}\frac{|\mathcal{L}_{\eta}[\tilde{u}](x)-\mathcal{L}_{\eta}[\tilde{u}](y)|}{|x-y|^{\beta}}$
$\displaystyle=R^{2s+\beta}\sup_{x,y\in
B_{R}(x_{0}))}\frac{\left|\mathcal{L}_{\eta}[u](x)-\mathcal{L}_{\eta}[u](y)\right|}{|x-y|^{\beta}}.$

We also observe that, since the coefficient $\eta$ is independent of $x$, for
any $k\in{\mathbb{N}}_{0}$

$\displaystyle
D^{k}\mathcal{L}_{\eta}[\tilde{u}](x)=R^{k}\mathcal{L}_{\eta}[D^{k}u](x).$

We use this in the case $\beta>1$ to assert that, for any $\beta>0$ that
satisfies (3),

(20)
$\displaystyle\|\mathcal{L}_{\eta}[\tilde{u}]\|_{C^{\beta}(B_{1})}\lesssim
R^{2s+\beta}\|\mathcal{L}_{\eta}[u]\|_{C^{\beta}(B_{R}(x_{0}))}\leq
R^{s}\|\mathcal{L}_{\eta}[u]\|_{\beta,\Omega}^{(s)}.$

The estimates so far do not use the regularity of the coefficient $\eta$, but
only its boundedness and translation invariance. Now we can repeat the proof
[61, Corollary 2.4] using the additional regularity on $\eta$ together with
the estimates (18)—(20) to get

(21) $\displaystyle\|\tilde{u}\|_{C^{\beta}(B_{1/2})}$
$\displaystyle\lesssim\|(1+|\cdot|)^{-d-2s}\tilde{u}(\cdot)\|_{L^{1}(\mathbb{R}^{d})}+\|\tilde{u}\|_{C^{\beta-2s}(B_{2})}+\|\mathcal{L}_{\eta}[\tilde{u}]\|_{C^{\beta-2s}(B_{2})}$
$\displaystyle\lesssim
R^{s}\left([u]_{C^{0,s}(\mathbb{R}^{d})}+\|f\|_{\beta-2s,\Omega}^{(s)}\right)+\|\tilde{u}\|_{C^{\beta-2s}(B_{2})},$

where we assumed that $K\geq 2$ to bound the term
$\|\mathcal{L}_{\eta}[\tilde{u}]\|_{C^{\beta-2s}(B_{2})}$ appropriately. We
repeat that the regularity of the coefficient $\eta$ is only used in order to
be able to stress the arguments used in the proof of [61, Corollary 2.4].

Now, for $\beta-2s\leq s$, the claim follows from estimate (21). Indeed, using
that, if $y\in B_{R}(x_{0})$, we have

$\|D^{\ell}\tilde{u}\|_{L^{\infty}(B_{1/2})}=R^{\ell}\|D^{\ell}u\|_{L^{\infty}(B_{R/2}(x_{0}))},$

and $\delta(x,y)/|x-y|\leq K$ if $x$ and $y$ are far away from each other.

If, instead, we have that $\beta-2s>s$, we repeat the argument in (21) for
$\|\tilde{u}\|_{C^{\beta-2s}(B_{2})}$ in total $m$ times until $\beta-2sm\leq
s$. Choosing $K$ large enough, but finite (depending on $s$), finishes the
proof. ∎

As a consequence of Lemmas 8 and 10, we have the following regularity estimate
away from the boundary.

######  Lemma 11 (interior Hölder estimate).

Let $\Omega$ be a bounded Lipschitz domain that satisfies the exterior ball
condition and let $f$, $\beta$, and $\eta$ satisfy the same assumptions as in
Lemma 10. For every $\rho>0$ we have

$\|u\|_{C^{\beta}(\left\\{x\in\Omega:\delta(x)\geq\rho\right\\})}\lesssim\rho^{s-\beta},$

for a constant only on
$\|f\|_{\beta-2s,\Omega}^{(s)},\|f\|_{L^{\infty}(\Omega)}$, $s$, and $\Omega$.

###### Proof.

Repeat the proof of [39, Corollary 2.5]. ∎

In short, the previous results show that $u$, the solution of (17), satisfies
$u\in{\mathcal{S}}^{\beta}(\Omega)$ and $\|u\|_{\beta,\Omega}^{(-s)}<\infty$.

###  4.2. Pointwise error estimates

We have reached the point where we can propose a numerical scheme for (17) and
provide an error analysis for it. To begin, we need to make a precise choice
of the regularization scale $\boldsymbol{\varepsilon}$.

######  Definition 12 ($\boldsymbol{\varepsilon}$ for the linear problem).

Let the vertex $z\in\mathscr{N}_{{\mathscr{T}}}^{i}$. We set

$\boldsymbol{\varepsilon}(z)=\frac{1}{2}h_{z}^{1/2}\delta(z)^{1/2}.$

At this point we make some comments regarding the choice in Definition 12.
First, notice that close to the boundary we have

$\boldsymbol{\varepsilon}(z)\approx h_{z}\approx\delta(z),$

which ensures that $B_{\boldsymbol{\varepsilon}(z)}(z)\subset\Omega$. On the
other hand, in the interior of the domain we have $\delta(z)\approx 1$, so
that

$\boldsymbol{\varepsilon}(z)\approx h_{z}^{1/2}.$

This scaling bears resemblance to the scalings used in two-scale methods for
(local) linear second order elliptic problems; see; [49, 50, 64, 46, 58, 53,
54].

We now describe the scheme. For $\boldsymbol{\varepsilon}$ given according to
Definition 12, our scheme seeks for $u_{h}\in{\mathbb{V}}_{\mathscr{T}}^{0}$
such that

(22)
${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[u_{h}](z)=f(z),\quad\forall
z\in\mathscr{N}_{{\mathscr{T}}}^{i}.$

We comment that, since this is a finite dimensional problem, existence and
uniqueness of a solution are implied immediately by the following discrete
comparison principle for the operator
${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}$.

######  Proposition 13 (discrete comparison principle).

Let $\boldsymbol{\varepsilon}$ be such that, for every
$z\in\mathscr{N}_{{\mathscr{T}}}^{i}$, we have
$\boldsymbol{\varepsilon}(z)\geq\tfrac{1}{2}h_{z}$. Assume that
$v_{h},w_{h}\in{\mathbb{V}}_{\mathscr{T}}^{0}$ are such that

(23)
${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[v_{h}](z)\geq{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[w_{h}](z),\quad\forall
z\in\mathscr{N}_{{\mathscr{T}}}^{i}.$

Then,

$v_{h}(z)\geq w_{h}(z)\quad\forall z\in\mathscr{N}_{{\mathscr{T}}}^{i}.$

###### Proof.

The proof follows by a simple argument. Suppose the inequality (23) is strict,
and that the function $v_{h}-w_{h}$ attains a (non-positive) minimum at an
interior node $z\in\mathscr{N}_{{\mathscr{T}}}^{i}.$ Then, we have

$v_{h}(z)-w_{h}(z)\leq v_{h}(y)-w_{h}(y)\ \implies\ v_{h}(z)-v_{h}(y)\leq
w_{h}(z)-w_{h}(y)\quad\forall y\in{\mathbb{R}}^{d}$

and, consequently,

${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[v_{h}](z)\leq{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[w_{h}](z).$

This shows that, in case of strict inequality in (23), it must be that
$v_{h}>w_{h}$ in $\Omega$.

Assume next that the inequality (23) is not strict. Consider the discrete
barrier function $b_{h}=I_{h}b$, where $b=\chi_{\Omega}$ and $I_{h}$ denotes
the Lagrange interpolant. We have
${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[b_{h}](z)>0$ for all
$z\in\mathscr{N}_{{\mathscr{T}}}^{i}$. Therefore, for $\epsilon>0$, we have
the strict inequality

${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[v_{h}+\epsilon
b_{h}](z)>{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[w_{h}](z)\quad\forall
z\in\mathscr{N}_{{\mathscr{T}}}^{i},$

from which it follows that $v_{h}(z)+\epsilon b_{h}(z)>w_{h}(z)$ for all
$z\in\mathscr{N}_{{\mathscr{T}}}^{i}$. Letting $\epsilon\to 0$, we obtain the
desired result. ∎

We now proceed with the error analysis of our scheme. Using the choice of
$\boldsymbol{\varepsilon}$ given by Definition 12, we begin by making the
results from Proposition 7 more precise.

######  Proposition 14 (consistency of interpolation for the linear problem).

Let the regularization scale $\boldsymbol{\varepsilon}$ verify Definition 12,
$\beta$ satisfy (3), $\bar{\beta}=\min\\{\beta,2\\}$, and
$w\in{\mathcal{S}}^{\beta}(\Omega)$. We have

$\left|{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[I_{h}w-w](z)\right|\lesssim\max\\{h^{\mu
s},h^{\bar{\beta}-s}\\}\delta(z)^{-2s}\quad\forall
z\in\mathscr{N}_{{\mathscr{T}}}^{i}.$

###### Proof.

Let $z\in\mathscr{N}_{{\mathscr{T}}}^{i}$. By (15), it suffices to show that

$h^{\bar{\beta}}\delta(z)^{s-\bar{\beta}/\mu}\boldsymbol{\varepsilon}(z)^{-2s}\lesssim\max\\{h^{\mu
s},h^{\bar{\beta}-s}\\}\delta(z)^{-2s}.$

We consider the set $\Omega^{\partial}$ given by (16).

If $z\in\Omega^{\partial}$, then $\delta(z)\approx h_{z}\approx h^{\mu}$.
Combining this with
$\boldsymbol{\varepsilon}(z)=\frac{1}{2}h_{z}^{1/2}\delta(z)^{1/2}\approx\delta(z)$,
we obtain

$h^{\bar{\beta}}\delta(z)^{s-\bar{\beta}/\mu}\boldsymbol{\varepsilon}(z)^{-2s}\approx
h^{\mu s}\delta(z)^{-2s}.$

In contrast, if $z\notin\Omega^{\partial}$, then $z\in T$ for some $T$ with
$T\cap\partial\Omega=\emptyset$ and therefore $h_{z}\approx h_{T}\approx
h\delta(z)^{1-1/\mu}$ and
$\boldsymbol{\varepsilon}(z)=\frac{1}{2}h_{z}^{1/2}\delta(z)^{1/2}$ If
$s\geq(\bar{\beta}-s)/\mu$, then we have

$h^{\bar{\beta}}\delta(z)^{s-\bar{\beta}/\mu}\boldsymbol{\varepsilon}(z)^{-2s}\approx
h^{\bar{\beta}-s}\delta(z)^{s-(\bar{\beta}-s)/\mu}\delta(z)^{-2s}\lesssim
h^{\bar{\beta}-s}\delta(z)^{-2s}.$

Otherwise, we write $h\approx h_{z}\delta(z)^{-1+1/\mu}$ and

$h^{\bar{\beta}}\delta(z)^{s-\bar{\beta}/\mu}\boldsymbol{\varepsilon}(z)^{-2s}\approx
h_{z}^{\bar{\beta}-s}\delta(z)^{-\bar{\beta}+2s}\delta(z)^{-2s}\approx h^{\mu
s}\left(\frac{h_{z}}{\delta(z)}\right)^{\bar{\beta}-s-\mu s}\delta(z)^{-2s}.$

Because $h_{z}\leq\delta(z)$ and $\bar{\beta}-s-\mu s>0$, the second term in
the right hand side is uniformly bounded above. ∎

It remains then to obtain error estimates. This is the content of the
following result.

######  Theorem 15 (error estimate).

Let $\Omega$ be a convex polytope, $s\in(0,1)$, $\beta\leq 4$ is such that (3)
holds. Define $\bar{\beta}=\min\\{\beta,2\\}$. Let $f\in
C^{\beta-2s}(\Omega)\cap{L^{\infty}(\Omega)}$ be such that
$\|f\|_{\beta-2s,\Omega}^{(s)}<\infty$. Assume that
$\eta\in{\mathcal{C}}(\lambda,\Lambda)$ and that if, for some
$k\in{\mathbb{N}}$, $\beta>k$, then $\eta\in C^{k}({\mathbb{S}}^{d-1})$. Let
$u\in{\mathcal{S}}^{\beta}(\Omega)$ solve (17) and
$u_{h}\in{\mathbb{V}}_{\mathscr{T}}^{0}$ solve (22). If ${\mathscr{T}}$
satisfies (11) and $\boldsymbol{\varepsilon}$ is chosen as in Definition 12,
we have

$\|u-u_{h}\|_{L^{\infty}(\Omega)}\lesssim\max\\{h^{\mu
s},h^{\bar{\beta}-s},h^{\beta/2-s}\\}.$

###### Proof.

We consider $u-u_{h}=\left(u-I_{h}u\right)+\left(I_{h}u-u_{h}\right)$. By
(14), to prove the claim, it suffices to estimate the second term. We do so by
estimating the consistency

$\left|{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[I_{h}u-u_{h}](z)\right|\leq\left|{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[I_{h}u-u](z)\right|+\left|{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[u-u_{h}](z)\right|=\mathrm{I}+\mathrm{II},$

for $z\in\mathscr{N}_{{\mathscr{T}}}^{i}$, and then applying the discrete
comparison principle from Proposition 13.

Proposition 14 yields

$\mathrm{I}\lesssim\max\\{h^{\mu
s},h^{\bar{\beta}-s}\\}\delta(z)^{-2s}\approx\max\\{h^{\mu
s},h^{\bar{\beta}-s}\\}{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[b](z),$

where we recall that $b$ is the barrier function introduced in Lemma 4.
Additionally, the choice of $\boldsymbol{\varepsilon}$, identity (13), and
Theorem 3 with $\alpha_{0}=2$, yield

$\displaystyle\mathrm{II}=\left|{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[u](z)-{\mathcal{L}}_{\eta}[u](z)\right|$
$\displaystyle\lesssim\boldsymbol{\varepsilon}(z)^{\beta-2s}\delta(z)^{s-\beta}$
$\displaystyle\approx
h^{\beta/2-s}\delta(z)^{s-\frac{\beta/2-s}{\mu}}{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[b](z).$

Thus, if $s-\frac{\beta/2-s}{\mu}\geq 0$ we are done. In case
$s<\frac{\beta/2-s}{\mu}$, we write instead $h\approx
h_{z}\delta(z)^{1/\mu-1}$ and obtain

$\displaystyle\mathrm{II}$ $\displaystyle\lesssim
h^{\beta/2-s}\delta(z)^{s-\frac{\beta/2-s}{\mu}}{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[b](z)$
$\displaystyle\approx h^{\mu s}h_{z}^{\beta/2-\mu
s-s}\delta(z)^{(1/\mu-1)(\beta/2-\mu
s-s)}\delta(z)^{s-\frac{\beta/2-s}{\mu}}{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[b](z)$
$\displaystyle\approx h^{\mu
s}\left(\frac{h_{z}}{\delta(z)}\right)^{\beta/2-\mu
s-s}{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[b](z).$

This concludes the proof. ∎

######  Remark 16 (complexity estimate).

Let us try to interpret the estimates of Theorem 15 in terms of degrees of
freedom. To shorten the notation we set, only for this discussion,
$N=\dim{\mathbb{V}}_{\mathscr{T}}^{0}$. Furthermore, we will assume that the
right hand side $f$ is as smooth as possible, i.e., we let $\beta=4$ and
satisfying (3). This regularity assumption implies that
$\beta/2=\bar{\beta}=2$.

First, if $d=2$, we can choose $\mu=2$ so that, according to (12), $N\approx
h^{-2}|\log h|$. As a consequence,

$\|u-u_{h}\|_{L^{\infty}(\Omega)}\lesssim\max\left\\{N^{-s}|\log
N|^{s},N^{-1+s/2}|\log N|^{1-s/2}\right\\}.$

Therefore, with respect to the number of degrees of freedom, and up to
logarithmic factors, we obtain convergence with order $s$ for $s\leq 2/3$ and
with order $1-s/2$ for $s>2/3$. Additionally we observe that, for $s>2/3$, one
does not need to take $\mu=2$ in two dimensions, and that the maximal
convergence rate is attained whenever $\mu s=2-s$ i.e., $\mu=\frac{2-s}{s}$;
any extra mesh refinement does not reflect in an improvement of convergence
rates, although it affects the conditioning of the resulting system.

On the other hand, for $d=3$ we must choose $\mu=\tfrac{3}{2}$ if we wish to
maintain a near optimal number of degrees of freedom, i.e., $N\approx
h^{-3}|\log h|$. Thus,

$\|u-u_{h}\|_{L^{\infty}(\Omega)}\lesssim\max\left\\{N^{-s/2}|\log
N|^{s/2},N^{-2/3+s/3}|\log N|^{2/3-s/3}\right\\}.$

Again, we observe that, if $-s/2\leq-2/3+s/3$, namely if $s\geq 4/5$, the
maximal convergence order $s/2$ is attained for meshes graded with
$\mu=\frac{2-s}{s}$.

######  Remark 17 (relationship between $\boldsymbol{\varepsilon}$ and $f$).

Definition 12 is suitable for sufficiently smooth right hand sides, namely, it
formally delivers optimal convergence rates in case $f\in
C^{\beta-2s}(\Omega)\cap{L^{\infty}(\Omega)}$ be such that
$\|f\|_{\beta-2s,\Omega}^{(s)}<\infty$ with $\beta\geq 4$. We recall that,
since we are using a second difference formula, the interior consistency of
the regularized operator cannot exploit any regularity beyond
${\mathcal{S}}^{4}(\Omega)$, cf. Theorem 3. Let us briefly comment on what one
can obtain when $f$ satisfies the assumptions above but with $\beta\in(2s,4)$
and satisfying (3).

We let
$\boldsymbol{\varepsilon}(z)=\frac{1}{2}h_{z}^{\alpha}\delta(z)^{1-\alpha}$
with $\alpha\in[0,1]$, that clearly satisfies
$\frac{1}{2}h_{z}\leq\boldsymbol{\varepsilon}(z)\leq\frac{1}{2}\delta(z)$. By
doing the same calculations as in Proposition 14, we obtain

$\left|{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[I_{h}w-w](z)\right|\lesssim\max\\{h^{\mu
s},h^{\bar{\beta}-2\alpha s}\\}\delta(z)^{-2s}\quad\forall
z\in\mathscr{N}_{{\mathscr{T}}}^{i}.$

Arguing then as in the proof of Theorem 15, we obtain the error estimate

$\|u-u_{h}\|_{L^{\infty}(\Omega)}\lesssim\max\\{h^{\mu
s},h^{\bar{\beta}-2\alpha s},h^{\alpha(\beta-2s)}\\}.$

Now, if $\beta\in(2,4)$, we have $\bar{\beta}=2$ and

$2-2\alpha s=\alpha(\beta-2s)\ \Rightarrow\alpha=\frac{2}{\beta}.$

Therefore, choosing
$\boldsymbol{\varepsilon}(z)=\frac{1}{2}h_{z}^{\frac{2}{\beta}}\delta(z)^{1-\frac{2}{\beta}}$
yields the error estimate

$\|u-u_{h}\|_{L^{\infty}(\Omega)}\lesssim\max\\{h^{\mu
s},h^{2-\frac{4s}{\beta}}\\}.$

In contrast, if $\beta\in(2s,2)$, we observe $\bar{\beta}=\beta$ and we get
that

$\beta-2\alpha s=\alpha(\beta-2s)\quad\implies\quad\alpha=1.$

In this low regularity case, setting
$\boldsymbol{\varepsilon}(z)=\frac{1}{2}h_{z}$ gives rise to

$\|u-u_{h}\|_{L^{\infty}(\Omega)}\lesssim\max\\{h^{\mu s},h^{\beta-2s}\\}.$

The latter will be of interest in the approximation of the obstacle problem in
the next section, and justifies Definition 20 below.

##  5\. The obstacle problem

As the next application of our two-scale discretization, we will consider the
following nonlinear problem. In the setting of Section 4 we assume that, in
addition, we have $\psi:\bar{\Omega}\to{\mathbb{R}}$ that satisfies $\psi<0$
on $\partial\Omega$. We seek for $u\in{\widetilde{H}}^{s}(\Omega)$ that
satisfies

(24)
$\min\left\\{{\mathcal{L}}_{\eta}[u]-f,u-\psi\right\\}=0,\quad\text{a.e.}\
\Omega.$

While existence and uniqueness of a weak solution is classical, the regularity
of such solution is more delicate. Following [17] we introduce the classes

(25)
${\mathcal{F}}_{s}(\bar{\Omega})=C^{3-2s+\epsilon}(\bar{\Omega}),\qquad\Psi=\left\\{\psi\in
C(\bar{\Omega}):\psi_{|\partial\Omega}<0\right\\}\cap C^{2,1}(\Omega),$

where $\epsilon>0$ is sufficiently small, so that
$1-2s+\epsilon\notin{\mathbb{N}}$.

For future use we define the contact and non-contact sets as follows:

$\Omega^{0}=\left\\{x\in\Omega\ \middle|\
u(x)=\psi(x)\right\\},\qquad\Omega^{+}=\left\\{x\in\Omega\ \middle|\
u(x)>\psi(x)\right\\}.$

In order to obtain rates of convergence, we must understand the regularity of
the solution. This can be achieved by combining the arguments in [17], [19],
and the regularity of the linear problem presented in Section 4.1. Namely, one
first proves the result for a problem in the whole space $\mathbb{R}^{d}$, and
then use a localization argument; see [17] for details in the case $\eta\equiv
1$.

######  Proposition 18 (regularity).

Assume that $f\in{\mathcal{F}}_{s}$ and $\psi\in\Psi$. Then,
$u\in{\widetilde{H}}^{s}(\Omega)$, the solution of (24) satisfies $u\in
C^{1,s}(\Omega)$ and

$\|{\mathcal{L}}_{\eta}[u]\|_{1-s,\Omega}^{(s)}<\infty.$

######  Remark 19 (pointwise evaluation).

Notice that, since $\psi<0$ on $\partial\Omega$, the solution to (24) solves
the linear problem ${\mathcal{L}}_{\eta}[u]=f$ in a neighborhood of the
boundary. By the interior regularity of the previous result we additionally
have ${\mathcal{L}}_{\eta}[u]\in C^{0,1-s}(\Omega)$. As a consequence, we have
that ${\mathcal{L}}_{\eta}[u]\in C^{0,1-s}(\bar{\Omega})$, meaning that
pointwise evaluation of the operator is meaningful.

###  5.1. Pointwise error estimates

Let us now provide a numerical scheme for the obstacle problem (24) and
pointwise error estimates for it. We seek for
$u_{h}\in{\mathbb{V}}_{\mathscr{T}}^{0}$ such that

(26)
$\min\left\\{{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[u_{h}](z)-f(z),u_{h}(z)-\psi(z)\right\\}=0,\quad\forall
z\in\mathscr{N}_{{\mathscr{T}}}^{i}.$

For this problem, we shall make a different choice of
$\boldsymbol{\varepsilon}$ than for the linear problem. The reason behind this
is that even if the data is sufficiently smooth, the solution to the obstacle
problem possesses a limited interior regularity, compare Lemma 10 with
Proposition 18; see also the discussion of the case $\beta\in(2s,2)$ in Remark
17.

######  Definition 20 (choice of $\boldsymbol{\varepsilon}$ for the obstacle
problem).

Let $z\in\mathscr{N}_{{\mathscr{T}}}^{i}$. We set

$\boldsymbol{\varepsilon}(z)=\frac{1}{2}h_{z}.$

Existence and uniqueness of $u_{h}$ follow from the fact that we are in finite
dimensions and the comparison principle for
${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}$. Of interest here is the
derivation of pointwise error estimates. The technique that we will use is
rather classical and can be traced back to [5, 51], see also [55]. We begin by
introducing the notions of sub- and supersolutions to the obstacle problem.

######  Definition 21 (sub- and supersolution).

We say that $u_{h}^{+}\in{\mathbb{V}}_{\mathscr{T}}$ is a supersolution to
(26) if $u_{h}^{+}\geq 0$ in $\Omega^{c}$ and, for all
$z\in\mathscr{N}_{{\mathscr{T}}}^{i}$, we have

$u_{h}^{+}(z)\geq\psi(z),\qquad{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[u_{h}^{+}](z)\geq
f(z).$

On the other hand, we say that $u_{h}^{-}\in{\mathbb{V}}_{\mathscr{T}}$ is a
subsolution to (26) if $u_{h}^{-}\leq 0$ in $\Omega^{c}$ and, for every
$z\in\mathscr{N}_{{\mathscr{T}}}^{i}$, if $u_{h}^{-}(z)\geq\psi(z)$, then

${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[u_{h}^{-}](z)\leq f(z).$

The comparison principle of the operator
${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}$ gives a comparison for sub-
and supersolutions.

######  Lemma 22 (discrete comparison).

Let $u_{h}^{+},u_{h}^{-}\in{\mathbb{V}}_{\mathscr{T}}$ be super- and
subsolutions to (26), and $u_{h}\in{\mathbb{V}}_{\mathscr{T}}^{0}$ be the
solution to (26). Then, for every $z\in\mathscr{N}_{{\mathscr{T}}}$ we have

$u_{h}^{-}(z)\leq u_{h}(z)\leq u_{h}^{+}(z).$

###### Proof.

We consider each inequality separately. Let $u_{h}^{-}$ be a subsolution and
consider the set of nodes

$C_{-}=\left\\{z\in\mathscr{N}_{{\mathscr{T}}}\ \middle|\
u_{h}^{-}(z)\geq\psi(z)\right\\}.$

Now, if $z\in C_{-}$, we have

${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[u_{h}^{-}](z)\leq
f(z)\leq{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[u_{h}](z).$

If, on the other hand $z\notin C_{-}$, then

$u_{h}^{-}(z)<\psi(z)\leq u_{h}(z).$

In summary, the function $w_{h}=u_{h}^{-}-u_{h}\in{\mathbb{V}}_{\mathscr{T}}$
verifies

${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[w_{h}](z)\leq 0,\quad z\in
C_{-},\qquad w_{h}(z)\leq 0,\quad z\notin C_{-}.$

A discrete comparison principle then implies that $w_{h}\leq 0$.

Let now $u_{h}^{+}$ be a supersolution. Consider now the discrete contact set

$C_{+}=\left\\{z\in\mathscr{N}_{{\mathscr{T}}}\ \middle|\
u_{h}(z)=\psi(z)\right\\},$

and observe that, if $z\in C_{+}$,

$u_{h}(z)=\psi(z)\leq u_{h}^{+}(z).$

On the other hand, if $z\notin C_{+}$ we have

${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[u_{h}^{+}](z)\geq
f(z)={\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[u_{h}](z).$

In conclusion, the function
$w_{h}=u_{h}^{+}-u_{h}\in{\mathbb{V}}_{\mathscr{T}}$ satisfies

${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[w_{h}]\geq 0,\quad z\notin
C_{+},\qquad w_{h}(z)\geq 0,\quad z\in C_{+}.$

Once again, a discrete comparison principle yields that $w_{h}\geq 0$. ∎

Next we need to present a suitable discrete proxy for $u$. Namely, we consider
$R_{h}u\in{\mathbb{V}}_{\mathscr{T}}^{0}$ to be such that, for every
$z\in\mathscr{N}_{{\mathscr{T}}}^{i}$,

(27)
${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[R_{h}u](z)={\mathcal{L}}_{\eta}[u](z).$

Recall that, as detailed in Remark 19, the right hand side is meaningful. The
approximation power of $R_{h}u$ is the content of the following result.

######  Corollary 23 (projection error).

Let $u$ be the solution of (24) and $R_{h}u$ be defined in (27). If the
regularization scale $\boldsymbol{\varepsilon}$ is chosen according to
Definition 20, then we have

(28) $\|u-R_{h}u\|_{L^{\infty}(\Omega)}\lesssim{\mathfrak{e}}(h)=\max\\{h^{\mu
s},h^{1-s}\\},$

with an implicit constant that is independent of $h$.

###### Proof.

We begin by recalling that, as indicated by Proposition 18, we must set
$\beta=1+s$.

Observe that, owing to (27),

$\mathrm{I}=\left|{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[u-R_{h}u]\right|=\left|{\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[u]-{\mathcal{L}}_{\eta}[u]\right|\lesssim\boldsymbol{\varepsilon}(z)^{1-s}\delta(z)^{-1},$

where, in the last step, we used Theorem 3.

Using the the current choice of $\boldsymbol{\varepsilon}$ we then continue
this estimate as

$\displaystyle\mathrm{I}$
$\displaystyle\lesssim\left(h\delta(z)^{1-\tfrac{1}{\mu}}\right)^{1-s}\delta(z)^{-1}=h^{1-s}\delta(z)^{-2s}\delta(z)^{s+\tfrac{s-1}{\mu}}\lesssim
h^{1-s}\delta(z)^{-2s},$

provided that $s+\tfrac{s-1}{\mu}\geq 0$.

If, on the other hand, we have $s+\tfrac{s-1}{\mu}<0$ we use that
$\delta(z)\gtrsim h_{z}\gtrsim h^{\mu}$ to obtain that

$h^{1-s}\delta(z)^{s+\tfrac{s-1}{\mu}}\lesssim
h^{1-s}h^{\mu\left(s+\tfrac{s-1}{\mu}\right)}\approx h^{\mu s},$

so that, in all cases,

$\mathrm{I}\lesssim\max\\{h^{\mu s},h^{1-s}\\}\delta(z)^{-2s}.$

The remaining of the proof follows exactly as the one of Theorem 15. ∎

With this proxy of the solution at hand we are ready to obtain error
estimates.

######  Theorem 24 (error estimates).

In the setting of Proposition 18 let $u$ solve (24) and
$u_{h}\in{\mathbb{V}}_{\mathscr{T}}^{0}$ solve (26). If ${\mathscr{T}}$
satisfies (11) and $\boldsymbol{\varepsilon}$ is chosen as in Definition 20,
we have that

$\|u-u_{h}\|_{L^{\infty}(\Omega)}=\tilde{\mathfrak{e}}(h)\lesssim{\mathfrak{e}}(h),$

where the quantity ${\mathfrak{e}}(h)$ is defined in (28).

###### Proof.

We will construct suitable super- and subsolutions to the discrete obstacle
problem (26) and apply the comparison principle of Lemma 22 to conclude.
Notice that, owing to (28), there is a sufficiently large $C>0$ for which

$u-C{\mathfrak{e}}(h)\leq R_{h}u\leq u+C{\mathfrak{e}}(h).$

Supersolution: Let
$u_{h}^{+}=R_{h}u+C_{1}{\mathfrak{e}}(h)\in{\mathbb{V}}_{\mathscr{T}}$, where
the constant $C_{1}>0$ is to be chosen, and notice that $u_{h}^{+}\geq 0$ in
$\Omega^{c}$. Moreover, for $z\in\mathscr{N}_{{\mathscr{T}}}^{i}$, with the
barrier function $b=\chi_{\Omega}$ and Lemma 4,

${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[u_{h}^{+}](z)={\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[R_{h}u](z)+C_{1}{\mathfrak{e}}(h){\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[b](z)>{\mathcal{L}}_{\eta}[u](z)\geq
f(z).$

In addition, if $C_{1}\geq C$, we have, for any
$z\in\mathscr{N}_{{\mathscr{T}}}^{i}$,

$u_{h}^{+}(z)\geq
u(z)+(C_{1}-C){\mathfrak{e}}(h)\geq\psi(z)+(C_{1}-C){\mathfrak{e}}(h)\geq\psi(z).$

We have then shown that $u_{h}^{+}$ is a supersolution for the obstacle
problem. This implies, via Lemma 22, that for every $x\in\Omega$ we have

$u_{h}(x)\leq u_{h}^{+}(x)\leq u(x)+(C+C_{1}){\mathfrak{e}}(h),$

so that

$u_{h}(x)-u(x)\lesssim{\mathfrak{e}}(h).$

Subsolution: We now define
$u_{h}^{-}=R_{h}u-C_{2}{\mathfrak{e}}(h)\in{\mathbb{V}}_{\mathscr{T}}$ where
$C_{2}>0$ is to be chosen. Notice that $u_{h}^{-}\leq 0$ in $\Omega^{c}$. We
will show that it is a subsolution. To see this, let
$z\in\mathscr{N}_{{\mathscr{T}}}^{i}$ and assume that
$u_{h}^{-}(z)\geq\psi(z)$. Then,

$\psi(z)\leq u_{h}^{-}(z)\leq u(z)-(C_{2}-C){\mathfrak{e}}(h)<u(z),$

provided $C_{2}>C$. The fact that this inequality is strict shows that

${\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[u_{h}^{-}](z)={\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[R_{h}u](z)-C_{2}{\mathfrak{e}}(h){\mathcal{L}}_{\eta,\boldsymbol{\varepsilon}}[b](z)<{\mathcal{L}}_{\eta}[u](z)=f(z),$

so that indeed this is a subsolution. We invoke once again Lemma 22 to obtain
that

$u(x)-(C+C_{2}){\mathfrak{e}}(h)\leq u_{h}^{-}(x)\leq u_{h}(x),\quad\forall
x\in\Omega,$

as we needed to show. ∎

######  Remark 25 (optimality and complexity).

We comment that the rate of convergence of Theorem 24, as expressed by the
quantity ${\mathfrak{e}}(h)$ is optimal for our proof technique. To see this,
we recall that near the boundary the solution to the obstacle problem behaves
like that of the linear problem, i.e., $u(z)\approx\delta(z)^{s}$, so that the
rate of interpolation is at best $h^{\mu s}$. On the other hand, the interior
regularity of the solution is, at best, $C^{1,s}$. Since our operator
${\mathcal{L}}_{\eta}$ is of order $2s$, and our proof technique is based on
comparison principles, the rate in the interior can be at best
$h^{1+s-2s}=h^{1-s}$, as we have obtained.

Finally, with the notation and conventions of Remark 16 let us present the
following complexity estimate

$\|u-u_{h}\|_{L^{\infty}(\Omega)}\lesssim\begin{dcases}\max\left\\{N^{-s}|\log
N|^{s},N^{(s-1)/2}|\log N|^{(1-s)/2}\right\\},&d=2,\\\
\max\left\\{N^{-s/2}|\log N|^{s/2},N^{(s-1)/3}|\log
N|^{(1-s)/3}\right\\},&d=3.\end{dcases}$

###  5.2. Regularity of free boundaries

Of particular importance in applications is the so-called free boundary, which
is the boundary of the contact set

$\Gamma=\partial\Omega^{0}\cap\Omega.$

In this section we are concerned with the regularity of this set.

####  5.2.1. Regular points

We begin with the regularity of the free boundary near _regular points_ which,
roughly speaking, are those at which the function $u-\psi$ grows at a rate
$1+s$. Define $d(x)=\operatorname*{dist}(x,\Omega^{0})$ to be the distance
from $x\in\Omega$ to the contact set. The following result is proved in [19,
Theorem 1.1].

######  Theorem 26 (regular points I).

Let $\alpha\in(0,\min\\{s,1-s\\})$. Let $\psi\in\Psi$ and $x_{0}\in\Gamma$ be
a regular point. Then, there exist positive constants $a(x_{0})$ and
$r(x_{0})$ such that

$u(x)-\psi(x)=a(x_{0})d(x)^{1+s}+{\mathfrak{o}}(|x-x_{0}|^{1+s+\alpha})$

for all $x\in B_{r(x_{0})}(x_{0})\cap\Omega^{+}$. Moreover, the set of points
satisfying this property is an open set of $\Gamma$ and is locally a
$C^{1,\gamma}$ graph for all $\gamma\in(0,s)$. Finally,

$u\in C^{1,s}\left(B_{r(x_{0})}(x_{0})\right)$

######  Remark 27 (singular points).

The points of $\Gamma$ that are not regular are called _singular_. According
to [19, Theorem 1.1], these points satisfy

$u(x)-\psi(x)={\mathfrak{o}}(|x-x_{0}|^{1+s+\alpha}).$

Singular points do in fact occur and are characterized in [6] for the case
$\eta\equiv 1$.

######  Remark 28 ($C^{1,\gamma}$ smoothness).

Consider an interface $\Gamma$, a point $x_{0}\in\Gamma$ with normal vector
$\boldsymbol{\nu}$, and $x=x_{0}+r_{0}\boldsymbol{\nu}$ for a sufficiently
small $r_{0}>0$. If the interface $\Gamma$ was $C^{2}$, then the closest point
to $x$ in $\Gamma$ is $x_{0}$. However, if $\Gamma$ is of class $C^{1,\gamma}$
with $\gamma<1$ (as in the conclusion of Theorem 26 ), then this is not the
case anymore and the distance $d(x,\Gamma)$ may be realized at a point
different than $x_{0}$.

The following variant of Theorem 26, which avoids the use of the distance
function $d$ is stated in [36, Theorem 4.4.1].

######  Theorem 29 (regular points II).

Let $\alpha\in(0,\min\\{s,1-s\\})$, $\theta>\max\\{0,2s-1\\}$, and
$\gamma\in(0,s)$. Let the obstacle $\psi\in C_{0}^{2,\theta}(\mathbb{R}^{d})$.
Let $x_{0}\in\Gamma$ be a free boundary point. Then, there is $r_{0}>0$ such
that, for all $x\in B_{r_{0}}(x_{0})$, we have:

  1. (i)

either

(29)
$u(x)-\psi(x)=a_{0}\left((x-x_{0})\cdot\boldsymbol{\nu}\right)_{+}^{1+s}+{\mathcal{O}}(|x-x_{0}|^{1+s+\gamma}),$

for some $a_{0}>0$, $\boldsymbol{\nu}\in{\mathbb{S}}^{d-1}$, and $\gamma>0$,

  2. (ii)

or

(30) $u(x)-\psi(x)={\mathcal{O}}(|x-x_{0}|^{1+s+\alpha}).$

Moreover, the set of points satisfying (29) (regular points) is an open set of
$\Gamma$, and it is locally a $C^{1,\gamma}$ manifold.

The following result is presented in [36, Corollary 4.5.3] without proof.
Since this estimate will be useful in our constructions we present a proof.

######  Lemma 30 (Hölder continuity).

Let $x_{0}\in\Gamma$ be a regular point, namely, one that satisfies (29). In
the setting of Theorem 29, the vector and scalar

${{\mathbf{b}}}(z)=\lim_{x\to z}\frac{1}{d(x)^{s}}\nabla(u(x)-\psi(x)),\qquad
a(z)=\frac{|{{\mathbf{b}}}(z)|}{1+s}$

are of class $C^{0,\gamma}$ in $B_{r_{0}}(x_{0})\cap\Gamma$.

###### Proof.

Let $x_{1},x_{2}\in\Gamma\cap B_{r_{0}}(x_{0})$ be two arbitrary free boundary
points which we assume to be regular points and so to satisfy Theorem 29.
Since the set of regular points is open in $\Gamma$, this may require further
restricting $r_{0}>0$.

Denote $v=u-\psi$. It is possible to show that [36, Proposition 4.4.15]

$\left\|\frac{1}{d^{s}}\nabla v\right\|_{C^{0,\gamma}(B_{r_{0}}(x_{0}))}\leq
c_{0},$

whence

$\left|\frac{1}{d(x)^{s}}\nabla v(x)-\frac{1}{d(y)^{s}}\nabla v(y)\right|\leq
c_{0}|x-y|^{\gamma},\qquad\forall x,y\in B_{r_{0}}(x_{0}).$

Since $\Gamma$ is $C^{1,\gamma}$ within $B_{r_{0}}(x_{0})$ we let
$\boldsymbol{\nu}_{1},\boldsymbol{\nu}_{2}\in{\mathbb{S}}^{d-1}$ be the unit
normals to $\Gamma$ at $x_{1},~{}x_{2}$ pointing towards $\Omega^{+}$. Let
$x\to x_{1}$, $y\to x_{2}$ and use the fact that

$\lim_{x\to x_{1}}\frac{1}{d(x)^{s}}\nabla
v(x)={{\mathbf{b}}}(x_{1}),\qquad\lim_{y\to x_{2}}\frac{1}{d(y)^{s}}\nabla
v(y)={{\mathbf{b}}}(x_{2}),$

to deduce that

$\left|{{\mathbf{b}}}(x_{1})-{{\mathbf{b}}}(x_{2})\right|\leq
c_{0}|x_{1}-x_{2}|^{\gamma},$

whence ${{\mathbf{b}}}\in C^{0,\gamma}(B_{r_{0}}(x_{0})\cap\Gamma)$. Since
$a(x_{i})=\tfrac{1}{1+s}|{{\mathbf{b}}}(x_{i})|$, we infer that

$|a(x_{1})-a(x_{2})|\lesssim||{{\mathbf{b}}}(x_{1})|-|{{\mathbf{b}}}(x_{2})||\leq|{{\mathbf{b}}}(x_{1})-{{\mathbf{b}}}(x_{2})|\leq
c_{0}|x_{1}-x_{2}|^{\gamma}.$

This is the desired estimate for $a$. ∎

To exploit the previous result we make the following convenient, but
realistic, regularity assumption on $\Gamma$.

######  Assumption 31 (regular points).

The free boundary $\Gamma$ consists only of regular points, namely those that
satisfy Theorem 26 or (29).

Next we discuss the fundamental _nondegeneracy properties_ (NDP). Given
$\varepsilon>0$ we let ${\mathfrak{S}}(\Gamma,\varepsilon)$ denote a strip of
thickness $\varepsilon$ around the free boundary $\Gamma$, namely,

${\mathfrak{S}}(\Gamma,\varepsilon)=\left\\{x\in\Omega\ \middle|\
d(x)=\operatorname*{dist}(x,\Gamma)<\varepsilon\right\\},\qquad{\mathfrak{S}}^{+}(\Gamma,\varepsilon)={\mathfrak{S}}(\Gamma,\varepsilon)\cap\Omega^{+}.$

The following result is known as an NDP in distance. It prescribes a pointwise
behavior of $u-\psi$ (growth with rate at least $1+s$) if one is close to a
regular point of the free boundary $\Gamma$.

######  Corollary 32 (NDP in distance).

Let $K\Subset\Omega$ and set $\tilde{\Gamma}=\Gamma\cap K$. If Assumption 31
is valid, then there exists constants $a,\varepsilon_{0}>0$ such that

$u(x)-\psi(x)\geq ad(x)^{1+s},\qquad\forall
x\in{\mathfrak{S}}^{+}(\tilde{\Gamma},\varepsilon_{0}).$

###### Proof.

We proceed in several steps.

  1. 1.

Since, by Assumption 31, every point $x_{0}\in\tilde{\Gamma}$ is regular we
have that $a(x_{0})=\tfrac{|{{\mathbf{b}}}(x_{0})|}{1+s}>0$, and $a$ and
${{\mathbf{b}}}$ are of class $C^{0,\gamma}$ in view of Lemma 30. We then
deduce that

$a=\frac{1}{2}\min_{x_{0}\in\tilde{\Gamma}}a(x_{0})>0,$

because $\tilde{\Gamma}$ is compact.

  2. 2.

Let $v=u-\psi$, $x_{0}\in\tilde{\Gamma}$, and $r_{0}>0$ be such that

$\left\|\frac{1}{d^{s}}\nabla v\right\|_{C^{0,\gamma}(B_{r_{0}}(x_{0}))}\leq
C_{0}=C(x_{0}).$

Given $x\in B_{r_{0}/2}(x_{0})$ we let $x_{1}\in\tilde{\Gamma}\cap
B_{r_{0}}(x_{0})$ be a point at minimal distance, i.e.,

$d(x)=|x-x_{1}|=(x-x_{1})\cdot\boldsymbol{\nu}_{1}.$

Then, for $y$ close to $x_{1}$, we have

$\left|\frac{1}{d(x)^{s}}\nabla v(x)-\frac{1}{d(y)^{s}}\nabla v(y)\right|\leq
C_{0}|x-y|^{\gamma},$

whence, upon computing the limit as $y\to x_{1}$, we get

$\left|\frac{1}{d(x)^{s}}\nabla v(x)-{{\mathbf{b}}}_{1}\right|\leq
C_{0}|x-x_{1}|^{\gamma},$

This implies that

$\frac{1}{d(x)^{s}}\partial_{\boldsymbol{\nu}_{1}}v(x)\geq{{\mathbf{b}}}_{1}\cdot\boldsymbol{\nu}_{1}-C_{0}|x-x_{1}|^{\gamma}=|{{\mathbf{b}}}_{1}|-C_{0}|x-x_{1}|^{\gamma}\geq(1+s)a,$

provided $r_{0}^{\gamma}C_{0}\leq(1+s)a$, upon restricting $r_{0}$ if
necessary. Therefore, we deduce the nondegeneracy property for the normal
derivative

$\partial_{\boldsymbol{\nu}_{1}}v(x)\geq(1+s)ad(x)^{s}=(1+s)a\left((x-x_{1})\cdot\boldsymbol{\nu}_{1}\right)_{+}^{s}.$

  3. 3.

Let $x(t)=tx+(1-t)x_{1}$ denote any point in the segment joining $x_{1}$ and
$x$, $t\in[0,1]$. Then $x_{1}\in\tilde{\Gamma}$ is again a point in
$\tilde{\Gamma}$ at a minimal distance to $x(t)$. The previous point implies
then that

$\frac{\mathop{}\\!\mathrm{d}v(x(t))}{\mathop{}\\!\mathrm{d}t}=\partial_{\boldsymbol{\nu}_{1}}v(x(t))|x-x_{1}|\geq(1+s)at^{s}|x-x_{1}|^{1+s},$

or

$v(x(t))=v(x(1))-v(x(0))=\int_{0}^{1}\frac{\mathop{}\\!\mathrm{d}v(x(t))}{\mathop{}\\!\mathrm{d}t}\mathop{}\\!\mathrm{d}t\geq
a|x-x_{1}|^{1+s},$

for all $x\in B_{r_{0}/2}(x_{0})$. This is the desired local nondegeneracy
property.

  4. 4.

We cover $\tilde{\Gamma}$ with balls $B_{r_{0}/2}(x_{0})$ for every
$x_{0}\in\tilde{\Gamma}$. Since $\tilde{\Gamma}$ is compact, there is a finite
subcovering

$\tilde{\Gamma}\subset\bigcup_{m=1}^{M}B_{r_{m}/2}(x_{m}).$

Finally, let $\varepsilon_{0}>0$ be the distance from $\tilde{\Gamma}$ to the
complement of $\cup_{m=1}^{M}B_{r_{m}/2}(x_{m})$. Then every
$x\in{\mathfrak{S}}^{+}(\tilde{\Gamma},\varepsilon_{0})$ belongs to a ball
$B_{r_{m}/2}(x_{m})$ for which the previous step applies.

This concludes the proof. ∎

######  Remark 33 (NDP).

Observe that Corollary 32 implies the following weaker form of nondegeneracy:
For all $x_{0}\in\tilde{\Gamma}$ and $r\in(0,\varepsilon_{0}]$ we have

(31) $\sup_{x\in B_{r}(x_{0})}\left(u(x)-\psi(x)\right)\geq ar^{1+s}.$

This inequality is also valid for all
$\varepsilon_{0}<r\leq\operatorname*{diam}(\Omega)$ because

$\varepsilon_{0}=\frac{\varepsilon_{0}}{r}r\geq\frac{\varepsilon_{0}}{\operatorname*{diam}(\Omega)}r,$

which yields (31) with the constant $a$ replaced by

${\widetilde{a}}=a\left(\frac{\varepsilon_{0}}{\operatorname*{diam}(\Omega)}\right)^{1+s}.$

We now make an explicit assumption about the boundary behavior of $u-\psi$,
which is not only useful for the subsequent argument, but it also has been
used in deriving regularity via a localization argument; see Proposition 18.

######  Assumption 34 (boundary behavior).

The obstacle is strictly negative on $\partial\Omega$, i.e., there is
$c_{0}>0$ for which

$\psi(x)\leq-c_{0}<0,\quad\forall x\in\partial\Omega.$

The assumption above, by continuity, assumes that the free boundary $\Gamma$
is uniformly away from the boundary of the domain $\partial\Omega$, and thus
the problem is linear in a neighborhood of $\partial\Omega$; see [17] for
details in the case $\eta\equiv 1$. In addition, this implies that the NDP in
distance of Corollary 32 holds for all
$x\in{\mathfrak{S}}^{+}(\Gamma,\varepsilon_{0})$.

Let now $\varepsilon\in(0,\varepsilon_{0}]$. We define a
$\varepsilon$–neighborhood of the free boundary $\Gamma$

${\mathfrak{N}}(\Gamma,\varepsilon)=\left\\{x\in\Omega^{+}\ \middle|\
0<u(x)-\psi(x)<\varepsilon^{1+s}\right\\}.$

######  Corollary 35 (comparing ${\mathfrak{N}}(\Gamma,\varepsilon)$ and
${\mathfrak{S}}(\Gamma,\varepsilon)$).

If Assumptions 31 and 34 hold, then there is
$\varepsilon_{1}\in(0,\varepsilon_{0}]$, with $\varepsilon_{0}>0$ defined in
Corollary 32, such that

${\mathfrak{N}}\left(\Gamma,a^{\tfrac{1}{1+s}}\varepsilon\right)\subset{\mathfrak{S}}(\Gamma,\varepsilon),\qquad\forall\varepsilon\in(0,\varepsilon_{1}].$

###### Proof.

Consider the compact set $\omega$ between
${\mathfrak{S}}^{+}(\Gamma,\varepsilon_{0})$ and $\partial\Omega$, namely,

$\omega=\overline{\Omega\setminus\left({\mathfrak{S}}^{+}(\Gamma,\varepsilon_{0})\cap\Omega^{0}\right)}.$

Set again $v=u-\psi$. In view of Assumption 34 and the fact that $v\in
C^{0,s}(\bar{\Omega})$, we deduce

$v\geq c_{0},\text{ on }\partial\Omega,\qquad v\geq
a\varepsilon_{0}^{1+s}\text{ on
}\partial{\mathfrak{S}}^{+}(\Gamma,\varepsilon_{0})\cap\Omega^{+},$

and that there is a constant $c_{1}>0$ such that

$v(x)\geq c_{1},\qquad\forall x\in\omega.$

Let $\varepsilon_{1}\in(0,\varepsilon_{0}]$ be given by
$a\varepsilon_{1}^{1+s}=c_{1}$. The definition of
${\mathfrak{S}}(\Gamma,\varepsilon)$ in conjunction with the NDP in distance
of Corollary 32 guarantee that for all $\varepsilon\in(0,\varepsilon_{1}]$

$v(x)\geq a\varepsilon^{1+s},\qquad\forall
x\in{\mathfrak{S}}(\Gamma,\varepsilon)^{c}\cap\Omega^{+}.$

Therefore, for all $x\in{\mathfrak{N}}(\Gamma,a^{\tfrac{1}{1+s}}\varepsilon)$
we have

$0<v(x)<\left(a^{\tfrac{1}{1+s}}\varepsilon\right)^{1+s}=a\varepsilon^{1+s},$

whence $x\in{\mathfrak{S}}(\Gamma,\varepsilon)$ as asserted. This concludes
the proof. ∎

We comment that the inclusion asserted in the previous Corollary is the
typical assumption that is made for the error estimates in distance for the
classical obstacle problem. See, for instance, [55, Section 2(d)].

####  5.2.2. Singular points

We now examine nondegeneracy for singular points. We base our discussion on
[6], which requires that Theorem 26 holds. We recall that a singular point
$x_{0}\in\Gamma$ corresponds to

$a(x_{0})=\frac{|{{\mathbf{b}}}(x_{0})|}{1+s}=0$

The characterization of singular points given in [6] hinges on the following
structural assumption.

######  Assumption 36 (singular points).

There is $c_{0}>0$ and $\gamma>0$ such that $\psi\in C^{3,\gamma}(\Omega)$,
$f+\Delta\psi\leq-c_{0}<0$ in $\\{x\in\Omega\ |\ \psi(x)>0\\}$ and
$\emptyset\neq\\{x\in\Omega\ |\ \psi(x)>0\\}\Subset\Omega$.

For the rest of the discussion regarding singular points we restrict our
attention to the _fractional Laplacian_ , i.e., $\eta\equiv 1$. The following
result is proved in [6, Lemma 3.2].

######  Lemma 37 (general growth).

If Assumption 36 holds, then there are constants $a,r_{0}>0$ such that, for
all $x_{0}\in\Gamma$

$\sup_{x\in B_{r}(x_{0})}\left(u(x)-\psi(x)\right)\geq ar^{2},\qquad\forall
r\in(0,r_{0}).$

Notice that this is similar to (31) but with exponent $2$ instead of $1+s$.

Consider now the following class of homogeneous polynomials of degree two:

${\mathbb{P}}_{2}^{+}=\left\\{p_{2}(x)=\frac{1}{2}x^{\intercal}{{\mathbf{A}}}x\
\middle|\ {{\mathbf{A}}}\in\mathbb{R}^{d\times
d}\setminus\\{{{\mathbf{O}}}\\},\
{{\mathbf{A}}}^{\intercal}={{\mathbf{A}}},\sigma({{\mathbf{A}}})\subset[0,\infty)\right\\}.$

The following result, which can be found in [6, Proposition 7.2], is crucial
to characterize singular points but does not play an important role in our
discussion.

$\Omega^{0}$$\Omega^{0}$$\Omega^{+}$$\Omega^{+}$$\Omega^{0}$$\Omega^{+}$$\Omega^{+}$
Figure 1. Examples of singular free boundary points where the function
$u-\psi$ has strict quadratic growth. We expect quadratic growth in the
direction of the red arrows. The case on the right is more degenerate than the
case on the left, and worse for the approximation of the free boundary
$\Gamma$ (depicted in blue). In both cases we have that
$\dim\ker{{\mathbf{A}}}=1$. A discrete free boundary $\Gamma_{\mathscr{T}}$,
which is at a distance ${\mathcal{O}}(\delta(h)^{1/2})$ is depicted in dashed
brown.

######  Proposition 38 (growth at singular points).

If Assumption 36 holds, then there is a modulus of continuity
$\omega:\mathbb{R}^{+}\to\mathbb{R}^{+}$ such that for any $x_{0}\in\Gamma$
singular, we have

$u(x)-\psi(x)=p_{2}^{x_{0}}(x-x_{0})+\omega(|x-x_{0}|)|x-x_{0}|^{2},$

for some $p_{2}^{x_{0}}\in{\mathbb{P}}_{2}^{+}$.

######  Remark 39 (growth at singular points).

In the setting of Proposition 38 it is important to notice that, if
$\dim\ker{{\mathbf{A}}}=k$, with $k\in\\{0,\ldots,d-1\\}$, then $u-\psi$
exhibits strict quadratic growth in the directions orthogonal to
$\ker{{\mathbf{A}}}$. Since there is at least one such direction, we conclude
that Proposition 38 implies Lemma 37.

The first example of this scenario is when $\dim\ker{{\mathbf{A}}}=0$. This
corresponds to an isolated contact point $x_{0}$, with the function $u-\psi$
growing quadratically in all directions emanating from $x_{0}$.

As a second example we consider the situations given in Figure 1. Although the
geometry is quite distinct, in both cases we have $\dim\ker{{\mathbf{A}}}=1$.
We expect quadratic growth in the directions of the red arrows. The case on
the right is more degenerate than the case on the left, and worse for the
approximation of the free boundary $\Gamma$. This will be discussed further
below.

###  5.3. Error estimates for free boundaries

Let us now put the pointwise error estimates of the previous sections to use
and, on the basis of the discussions of Section 5.2, obtain approximation
properties for the free boundary $\Gamma$. To start we mention that there is
$\psi_{h}=I_{h}\psi\in{\mathbb{V}}_{\mathscr{T}}$ (the Lagrange interpolant)
such that $\psi_{h}(z)=\psi(z)$ for all $z\in\mathscr{N}_{{\mathscr{T}}}$ and,
more importantly,

(32) $\|\psi-\psi_{h}\|_{L^{\infty}(\Omega)}=\sigma(h),$

for some function $\sigma:\mathbb{R}^{+}\to\mathbb{R}^{+}$ with
$\sigma(h)\downarrow 0$ as $h\downarrow 0$. Typically, and this shall be the
case if $\psi\in\Psi$,

$\sigma(h)\lesssim|\psi|_{W^{2,\infty}(\Omega)}h^{2}.$

Notice that (26), and its analysis, remain unchanged if we replace $\psi$ by
$\psi_{h}$.

The next step is to define the discrete free boundary $\Gamma_{h}$. To do so,
instead of looking at the zero level set of $u_{h}-\psi_{h}$ we consider the
level set at height

$\delta(h)=\tilde{\mathfrak{e}}(h)+\sigma(h),$

where $\tilde{\mathfrak{e}}$ was defined in Theorem 24, see also (28). Thus we
define the discrete noncontact and contact sets, respectively, to be

$\Omega_{\mathscr{T}}^{+}=\left\\{x\in\Omega\ \middle|\
u_{h}(x)-\psi_{h}(x)>\delta(h)\right\\},\quad\Omega_{\mathscr{T}}^{0}=\left\\{x\in\Omega\
\middle|\ u_{h}(x)-\psi_{h}(x)\leq\delta(h)\right\\}.$

The free boundary is then,

$\Gamma_{\mathscr{T}}=\partial\Omega^{+}_{\mathscr{T}}\cap\Omega=\partial\Omega_{\mathscr{T}}^{0}\cap\Omega.$

We intend to prove that, in a sense, $\Gamma$ and $\Gamma_{\mathscr{T}}$ are
close. This will be quantified by means of the _Hausdorff distance_.

######  Definition 40 (Hausdorff distance).

Let $A,B\subset{\mathbb{R}}^{d}$. Their Hausdorff distance is

$d_{H}(A,B)=\max\left\\{\max_{x\in A}\operatorname*{dist}(x,B),\max_{y\in
B}\operatorname*{dist}(y,A)\right\\}.$

The error estimate for free boundaries reads as follows.

######  Theorem 41 (free boundary approximation).

Let $\Gamma$ satisfy Assumptions 31 and 34, and let $a>0$ be given in
Corollary 32 and (31). Then, there is $h_{0}>0$ such that

$d_{H}(\Gamma,\Gamma_{\mathscr{T}})\leq\left(\frac{2}{a}\delta(h)\right)^{\tfrac{1}{1+s}},\qquad\forall
h\in(0,h_{0}].$

###### Proof.

We proceed in three steps.

  1. 1.

We first show that $\Omega^{0}\subset\Omega_{\mathscr{T}}^{0}$. To achieve
this consider $x\in\Omega^{0}$, so that $u(x)=\psi(x)$. Thus, using (28) and
(32) we have

$u_{h}(x)-\psi_{h}(x)=\left(u_{h}(x)-u(x)\right)+\left(\psi(x)-\psi_{h}(x)\right)\leq\tilde{\mathfrak{e}}(h)+\sigma(h)=\delta(h).$

In other words, $x\in\Omega_{\mathscr{T}}^{0}$.

  2. 2.

Next we show that

$\Gamma_{\mathscr{T}}\subset{\mathfrak{S}}\left(\Gamma,\left(\frac{2}{a}\delta(h)\right)^{\tfrac{1}{1+s}}\right).$

Let then $x\in\Gamma_{\mathscr{T}}$ and $x_{0}\in\Gamma$ be the closest point
to it, i.e.,

$|x-x_{0}|=d(x)=\operatorname*{dist}(x,\Gamma).$

From the previous step we know that $x\in\Omega^{+}$. Since

$u_{h}(x)-\psi_{h}(x)=\delta(h),$

and, in view of Theorem 24, $u(x)-u_{h}(x)\leq\tilde{\mathfrak{e}}(h)$ we
realize that

$\displaystyle u(x)-\psi(x)$
$\displaystyle\leq\left(u(x)-u_{h}(x)\right)+\left(u_{h}(x)-\psi_{h}(x)\right)+\left(\psi_{h}(x)-\psi(x)\right)$
$\displaystyle\leq\tilde{\mathfrak{e}}(h)+\delta(h)+\sigma(h)=2\delta(h).$

Let now $h_{0}>0$ be sufficiently small so that
$\delta(h_{0})\leq\tfrac{a}{2}\varepsilon_{1}^{1+s}$, where
$\varepsilon_{1}>0$ is given in Corollary 35. This implies

$x\in{\mathfrak{N}}(\Gamma,a^{\tfrac{1}{1+s}}\varepsilon_{1})\subset{\mathfrak{S}}(\Gamma,\varepsilon_{1}).$

Therefore, Corollary 32 yields

$u(x)-\psi(x)\geq a\operatorname*{dist}(x,\Gamma)^{1+s},$

whence,

$a\operatorname*{dist}(x,\Gamma)^{1+s}\leq 2\delta(h),$

i.e.,

$\operatorname*{dist}(x,\Gamma)\leq\left(\frac{2}{a}\delta(h)\right)^{\tfrac{1}{1+s}},$

as desired.

  3. 3.

We show that

$\Gamma\subset{\mathfrak{S}}\left(\Gamma_{\mathscr{T}},\left(\frac{2}{a}\delta(h)\right)^{\tfrac{1}{1+s}}\right).$

Indeed, let $x_{0}\in\Gamma$ and assume, for the sake of contradiction, that

$R=\operatorname*{dist}(x,\Gamma_{\mathscr{T}})>\left(\frac{2}{a}\delta(h)\right)^{\tfrac{1}{1+s}}.$

Notice that the first step of the proof yields that if
$B_{R}(x_{0})\subset\Omega_{\mathscr{T}}^{0}$, then

$u_{h}(y)-\psi_{h}(y)\leq\delta(h),\qquad\forall y\in B_{R}(x_{0}).$

We now recall (31) which is a consequence of Corollary 32: for all
$r\leq\operatorname*{diam}(\Omega)$

$\sup_{y\in B_{r}(x_{0})}\left(u(y)-\psi(y)\right)\geq ar^{1+s}.$

In other words, by compactness and continuity, there is $y\in B_{R}(x_{0})$
such that

$u(y)-\psi(y)\geq aR^{1+s}.$

On the other hand,

$\displaystyle u_{h}(y)-\psi_{h}(y)$
$\displaystyle=\left(u_{h}(y)-u(y)\right)+\left(u(y)-\psi(y)\right)+\left(\psi(y)-\psi_{h}(y)\right)$
$\displaystyle\geq-\tilde{\mathfrak{e}}(h)+aR^{1+s}-\sigma(h)>-\tilde{\mathfrak{e}}(h)+2\delta(h)-\sigma(h)>\delta(h),$

which is a contradiction. This proves the assertion and concludes the proof. ∎

######  Remark 42 (stability).

We observe that the nondegeneracy constant $a>0$, defined in Corollary 32,
acts as a stability parameter in the estimate of Theorem 41.

######  Remark 43 (localized estimate).

Let $K\Subset\Omega$ be a compact so that $\Gamma\cap K$ is made only of
regular points. We thus allow $\Gamma$ to have singular points in $\Gamma\cap
K^{c}$. Since the set of regular points is relatively open in $\Gamma$, we
realize that Corollary 32 is valid in $K$, i.e.,

$u(x)-\psi(x)\geq a_{K}d(x)^{1+s},\quad\forall
x\in{\mathfrak{S}}^{+}(\Gamma,\varepsilon)\cap K,$

for some constant $a_{K}>0$ that, in particular, depends on the compact $K$.
We can thus repeat the proof of Theorem 41 locally to deduce

$d_{H}(\Gamma\cap K,\Gamma_{\mathscr{T}}\cap
K)\leq\left(\frac{2}{a_{K}}\delta(h)\right)^{\tfrac{1}{1+s}},\qquad\forall
h\in(0,h_{0}].$

Theorem 41 and Remark 43 assume that, at least locally, there are no singular
free boundary points. Let us conclude the discussion by presenting some
results about the general case, namely when $\Gamma$ contains singular points.
An inspection of the proof of Theorem 41 shows that the second step cannot be
carried out, but the first and third one remain valid. In fact, the first step
hinges on the definition of $\Gamma_{\mathscr{T}}$ and the last step relies on
the growth condition (31) which, in principle, could be replaced by the
general growth condition provided in Lemma 37. This brings about the following
result.

######  Theorem 44 (error estimates for singular points).

If Assumption 36 holds and $\eta\equiv 1$, then for every $x\in\Gamma$ we have
that $x\in\Omega_{\mathscr{T}}^{0}$ and

(33)
$\operatorname*{dist}(x,\Gamma_{\mathscr{T}})\leq\left(\frac{2}{a}\delta(h)\right)^{\tfrac{1}{2}}.$

######  Remark 45 (a posteriori error estimation).

Estimate (33) establishes an error of $\Gamma$ relative to
$\Gamma_{\mathscr{T}}$. This is the spirit of an _a posteriori_ error
estimate. We refer to [57] for similar estimates for the classical obstacle
problem for the Laplace operator.

######  Remark 46 (regularity).

Estimate (33) requires no regularity of the free boundary $\Gamma$ beyond the
nondegeneracy property of Lemma 37, which relies on Assumption 36. The free
boundary regularity stated in Proposition 38 is not needed to assert (33). It
is then natural to wonder how the quadratic growth in certain directions
established in Proposition 38 can be exploited for free boundary
approximation. For instance, one may use the approximations that are depicted
in Figure 1.

In the left panel of Figure 1 the discrete free boundary is uniformly close to
$\Gamma$ and a global estimate of the form

$\operatorname*{dist}(x,\Gamma)\lesssim\delta(h)^{1/2},\qquad\forall
x\in\Gamma_{\mathscr{T}}$

is expected.

On the other hand, the right panel of Figure 1 shows that points of
$\Gamma_{\mathscr{T}}$ may be far away from $\Gamma$ and, thus, the estimate
above cannot be obtained.

To conclude the discussion of approximation of free boundaries at singular
points we mention that the scenarios depicted in Figure 1 also illustrate the
relevance of the localized error estimates alluded to in Remark 43.

##  6\. A concave, fully nonlinear, nonlocal, problem

As a final application of our two-scale discretization, we shall consider a
nonlocal Hamilton Jacobi Bellman equation. Let $s\in(0,1)$, and
$\eta_{i}\in{\mathcal{C}}(\lambda,\Lambda)$ ($i\in\\{1,2\\}$) for some
$\lambda,\Lambda$. Given $f\in C(\bar{\Omega})$ we must find
$u:{\mathbb{R}}^{d}\to{\mathbb{R}}$ such that

(34)
$\min\left\\{{\mathcal{L}}_{\eta_{1}}[u],{\mathcal{L}}_{\eta_{2}}[u]\right\\}=f,\quad\text{
in }\Omega,\qquad u=0,\quad\text{ in }\Omega^{c}.$

Equations of this form have gathered a lot of attention in recent times. We
refer the reader to, for instance, [21, 23, 44] for details regarding the
existence, uniqueness, and regularity of solutions. Regarding numerics, in the
case $\Omega=\mathbb{R}^{d}$, see [8, 59, 25, 10, 9, 27, 28]. To our
knowledge, however, no reference addresses the numerical approximation in the
case of a bounded domain, as in (34). Our goal here will be to provide then
the first convergent method for such problem. To achieve this we will get
inspiration from [47, 35, 2, 11], and relate this problem to a sequence of
obstacle problems.

###  6.1. Existence and uniqueness

Since it will be useful for our numerical purposes, we begin by discussing the
existence of solutions. Below, by $\div$ we mean the remainder of integer
division.

######  Theorem 47 (existence and uniqueness).

Assume that $f\in C(\bar{\Omega})$ and that ${\mathcal{L}}_{\eta_{i}}$ have a
common supersolution, i.e., there is $U\in C^{2}({\mathbb{R}}^{d})$ for which

${\mathcal{L}}_{\eta_{i}}[U]\geq f\quad\text{ in }\Omega,\qquad U\geq
0,\quad\text{ in }\Omega^{c}.$

Then, problem (34) has a unique solution. Moreover, this solution can be
obtained as the uniform limit of the following sequence:
$u_{0}\in{\widetilde{H}}^{s}(\Omega)$ solves

${\mathcal{L}}_{\eta_{1}}[u_{0}]=f,\quad\text{ in }\Omega.$

For $k\in{\mathbb{N}}$, let $i=k\div 2+1$. The function
$u_{k}\in{\widetilde{H}}^{s}(\Omega)$ solves

(35)
$\min\left\\{{\mathcal{L}}_{\eta_{i}}[u_{k}]-f,u_{k}-u_{k-1}\right\\}=0,\quad\text{
in }\Omega.$

###### Proof.

We split the proof in several steps.

Convergence: Notice that the sequence $\\{u_{k}\\}_{k\in{\mathbb{N}}_{0}}$, as
solutions of an obstacle problem with common right hand side, lie uniformly in
$C^{0,s}(\bar{\Omega})$. Moreover, by construction, the sequence is
nondecreasing, i.e., $u_{k}\geq u_{k-1}$ and bounded above. Indeed, since

${\mathcal{L}}_{\eta_{1}}[U]\geq f={\mathcal{L}}_{\eta_{1}}[u_{0}]\quad\text{
in }\Omega,\qquad U-u_{0}\geq 0,\quad\text{ in }\Omega^{c},$

we must have, by comparison, that $U\geq u_{0}$. Assume next that, for some
$k\in{\mathbb{N}}$, we have $u_{k-1}\leq U$. Let $x\in\Omega$. We either have
$u_{k}(x)=u_{k-1}(x)\leq U(x)$ or $u_{k}(x)>u_{k-1}(x)$. However, at such
points we must have

${\mathcal{L}}_{\eta_{i}}[u_{k}](x)=f(x)\leq{\mathcal{L}}_{\eta_{i}}[U](x).$

Therefore, if we denote $\Omega_{k}^{+}=\\{x\in\Omega\ |\
u_{k}(x)>u_{k-1}(x)\\}$ we see that

${\mathcal{L}}_{\eta_{i}}[U-u_{k}]\geq 0,\quad\text{ in }\Omega_{k}^{+},\qquad
U-u_{k}\geq 0,\quad\text{ in }(\Omega_{k}^{+})^{c}.$

Consequently, by comparison, we must have that $u_{k}\leq U$ in
$\Omega_{k}^{+}$. By Dini’s theorem then, there is $u\in C(\bar{\Omega})$ such
that $u_{k}\rightrightarrows u$.

Solution: The uniform convergence also implies that, for $i\in\\{1,2\\}$,

${\mathcal{L}}_{\eta_{i}}[u]\geq f,\quad\text{ in }\Omega.$

On the other hand, if we show that for every $k\in{\mathbb{N}}_{0}$

(36)
$\min\left\\{{\mathcal{L}}_{\eta_{1}}[u_{k}],{\mathcal{L}}_{\eta_{2}}[u_{k}]\right\\}\leq
f,\quad\text{ in }\Omega,$

we may pass to the limit and obtain that

${\mathcal{L}}_{\eta_{1}}[u]\geq f,\quad{\mathcal{L}}_{\eta_{2}}[u]\geq
f,\quad\min\left\\{{\mathcal{L}}_{\eta_{1}}[u],{\mathcal{L}}_{\eta_{2}}[u]\right\\}\leq
f,$

and thus $u$ must be a solution.

Proof of (36): We argue by induction. By the way the iterative scheme is
initialized, (36) holds for $k=0$. Assume now that the inequality holds for
some $k\in{\mathbb{N}}$. Consider the noncoincidence set

$x_{n}\in\Omega_{k}^{+}=\left\\{x\in\Omega\ \middle|\
u_{k+1}(x)>u_{k}(x)\right\\}.$

By the complementarity conditions we must have, for some $i\in\\{1,2\\}$, that

${\mathcal{L}}_{\eta_{i}}[u_{k+1}](x_{n})=f(x_{n}).$

If, on the other hand, we consider the coincidence set

$x_{c}\in\Omega_{k}^{0}=\left\\{x\in\Omega\ \middle|\
u_{k+1}(x)=u_{k}(x)\right\\}\cap\Omega,$

we see that $u_{k+1}(x_{c})=u_{k}(x_{c})$ and $u_{k+1}(y)\geq u_{k}(y)$ for
all $y\in\mathbb{R}^{d}$. By the inductive hypothesis, there is
$i_{0}\in\\{1,2\\}$ for which

${\mathcal{L}}_{\eta_{i_{0}}}[u_{k}](x_{c})\leq f(x_{c}),$

then

$\displaystyle{\mathcal{L}}_{\eta_{i_{0}}}[u_{k+1}](x_{c})$
$\displaystyle=\mathrm{v.p.}\int_{{\mathbb{R}}^{d}}\left(u_{k+1}(x_{c})-u_{k+1}(y)\right)\frac{1}{|x-y|^{d+2s}}\eta\left(\frac{x_{c}-y}{|x_{c}-y|}\right)\mathop{}\\!\mathrm{d}y$
$\displaystyle=\mathrm{v.p.}\int_{{\mathbb{R}}^{d}}\left(u_{k}(x_{c})-u_{k+1}(y)\right)\frac{1}{|x-y|^{d+2s}}\eta\left(\frac{x_{c}-y}{|x_{c}-y|}\right)\mathop{}\\!\mathrm{d}y$
$\displaystyle\leq\mathrm{v.p.}\int_{{\mathbb{R}}^{d}}\left(u_{k}(x_{c})-u_{k}(y)\right)\frac{1}{|x-y|^{d+2s}}\eta\left(\frac{x_{c}-y}{|x_{c}-y|}\right)\mathop{}\\!\mathrm{d}y$
$\displaystyle={\mathcal{L}}_{\eta_{i_{0}}}[u_{k}](x_{c})\leq f(x_{c}),$

as claimed.

Uniqueness: Follows by comparison. ∎

######  Remark 48 (obstacle).

Notice that the iterative scheme presented in Theorem 47 requires, at every
step the solution of an obstacle problem like the one described in Section 5.

######  Remark 49 (supersolution).

Theorem 47 relies on the existence of a common supersolution for the
operators. A possible common supersolution is

$U(x)=\frac{A}{2}|x|^{2}+B.$

The constant $B$ can be chosen so that $U\geq 0$ in $\Omega^{c}$, whereas we
can choose $A$, depending only on $d$, $s$, $\lambda$, $\Lambda$, and
$-\|f\|_{L^{\infty}(\Omega)}$, to obtain a supersolution.

######  Remark 50 (rate of convergence).

It is not known to us whether a rate of convergence for the iteration of
Theorem 47 can be established.

###  6.2. A convergent scheme

As a final application of our constructions we present a scheme for problem
(34). Inspired by the proof of Theorem 47 we consider the following iterative
scheme: $u_{h,0}\in{\mathbb{V}}_{\mathscr{T}}^{0}$ is such that

${\mathcal{L}}_{\eta_{1},\boldsymbol{\varepsilon}}[u_{h,0}](z)=f(z),\quad\forall
z\in\mathscr{N}_{{\mathscr{T}}}^{i}.$

For $k\in{\mathbb{N}}$, let $i=k\div 2+1$. The function
$u_{h,k}\in{\mathbb{V}}_{\mathscr{T}}^{0}$ solves

(37)
$\min\left\\{{\mathcal{L}}_{\eta_{i},\boldsymbol{\varepsilon}}[u_{h,k}](z)-f(z),u_{h,k}(z)-u_{h,k-1}(z)\right\\}=0,\quad\forall
z\in\mathscr{N}_{{\mathscr{T}}}^{i}.$

We immediately observe that, for every $k$, (37) has a unique solution. We
would now like to obtain convergence of $u_{h,k}$ to $u$, the solution to
(34). In order to achieve this, we introduce, for $k\in{\mathbb{N}}$, the
function ${\widetilde{u}}_{h,k}\in{\mathbb{V}}_{\mathscr{T}}^{0}$ that is the
solution of

(38)
$\min\left\\{{\mathcal{L}}_{\eta_{i},\boldsymbol{\varepsilon}}[{\widetilde{u}}_{h,k}](z)-f(z),{\widetilde{u}}_{h,k}(z)-I_{h}u_{k-1}(z)\right\\}=0,\quad\forall
z\in\mathscr{N}_{{\mathscr{T}}}^{i}.$

Notice that ${\widetilde{u}}_{h,k}$ is nothing but an approximation to the
function $u_{k}$ from the scheme of Theorem 47. As such we expect that
${\widetilde{u}}_{h,k}\to u_{k}$ with a given rate. We quantify this by
introducing the following assumption.

######  Assumption 51 (approximation).

There is a continuous function $\sigma:{\mathbb{R}}^{+}\to{\mathbb{R}}^{+}$
such that $\sigma(h)\downarrow 0$ as $h\downarrow 0$ for which

$\sup_{k\in{\mathbb{N}}_{0}}\|u_{k}-{\widetilde{u}}_{h,k}\|_{L^{\infty}(\Omega)}\leq\sigma(h),$

where $\\{u_{k}\\}_{k\in{\mathbb{N}}_{0}}$ are defined in Theorem 47 and
$\\{{\widetilde{u}}_{h,k}\\}_{k\in{\mathbb{N}}_{0}}$ in (38).

######  Remark 52 (smoothness and compatibility).

While we would like to assert that the rate $\sigma(h)$ of Assumption 51 is
that given by Theorem 24 we are unable to prove this. The reason is that the
error estimates of this Theorem hinge on the regularity of Proposition 18
which need the obstacle to belong to the class $\Psi$. This means that, for
every $k\in{\mathbb{N}}$, we must be able to assert that:

  1. 1.

$u_{k}\in C^{2,1}(\Omega)$. While we are not able to verify this directly, we
comment that this smoothness assumption is taken from [17], which in turn
follows the arguments of [24, 67]. This may not be sharp.

  2. 2.

$u_{k|\partial\Omega}<0$. However, we have $u_{k|\partial\Omega}=0$.

Clearly

$\tilde{\mathfrak{e}}(h)\leq\sigma(h),$

where $\tilde{{\mathfrak{e}}}(h)$ was defined in Theorem 24; see also (28).
Let us now show convergence.

######  Theorem 53 (convergence).

In the setting of Theorem 47 and under Assumption 51, let $u$ solve (34), and
$\\{u_{h,k}\\}_{k\in{\mathbb{N}}_{0},h>0}$ be the solutions to (37). Assume
that, as $k\uparrow\infty$ and $h\downarrow 0$, we have that

$k\,\sigma(h)\to 0.$

Then, $u_{h,k}\rightrightarrows u$.

###### Proof.

As mentioned above, the family
$\\{{\widetilde{u}}_{h,k}\\}_{k\in{\mathbb{N}}_{0}}\subset{\mathbb{V}}_{\mathscr{T}}^{0}$,
defined in (38), is an approximation of the obstacle problem defined in (35).

Consider now the difference ${\widetilde{u}}_{h,k}-u_{h,k}$. We claim that

(39)
$\|{\widetilde{u}}_{h,k}-u_{h,k}\|_{L^{\infty}(\Omega)}\leq\|u_{h,k-1}-I_{h}u_{k-1}\|_{L^{\infty}(\Omega)}.$

Indeed, if we define
$w_{h}={\widetilde{u}}_{h,k}+\|u_{h,k-1}-I_{h}u_{k-1}\|_{L^{\infty}(\Omega)}$
we see that

$w_{h}\geq{\widetilde{u}}_{h,k}=0,\quad\text{ in }\Omega^{c}.$

Moreover, for $z\in\mathscr{N}_{{\mathscr{T}}}^{i}$,

$w_{h}(z)\geq{\widetilde{u}}_{h,k}(z)+u_{h,k-1}(z)-I_{h}u_{k-1}(z)\geq
u_{h,k-1}(z)$

and ${\mathcal{L}}_{\eta_{i},\boldsymbol{\varepsilon}}[w_{h}](z)\geq f(z)$. In
other words, the function $w_{h}$ is a supersolution to the obstacle problem
(37). Lemma 22 then implies that $u_{h,k}\leq w_{h}$, i.e.,

$u_{h,k}-{\widetilde{u}}_{h,k}\leq\|u_{h,k-1}-I_{h}u_{k-1}\|_{L^{\infty}(\Omega)}.$

A similar argument shows the lower bound.

Using Assumption 51 we now iterate (39) to obtain

$\displaystyle\|u_{k}-u_{h,k}\|_{L^{\infty}(\Omega)}$
$\displaystyle\leq\|u_{k}-{\widetilde{u}}_{h,k}\|_{L^{\infty}(\Omega)}+\|{\widetilde{u}}_{h,k}-u_{h,k}\|_{L^{\infty}(\Omega)}$
$\displaystyle\leq\sigma(h)+\|u_{h,k-1}-I_{h}u_{k-1}\|_{L^{\infty}(\Omega)}$
$\displaystyle\leq\sigma(h)+\|u_{k-1}-I_{h}u_{k-1}\|_{L^{\infty}(\Omega)}+\|u_{h,k-1}-u_{k-1}\|_{L^{\infty}(\Omega)}$
$\displaystyle\leq 2\sigma(h)+\|u_{h,k-1}-u_{k-1}\|_{L^{\infty}(\Omega)}$
$\displaystyle\leq(k+1)\sigma(h)+\|u_{h,0}-u_{0}\|_{L^{\infty}(\Omega)}\leq(k+2)\sigma(h),$

where, in the last step, we used that $u_{h,0}$ is an approximation to the
solution of the linear problem.

The triangle inequality, and the assumption on $k$ and $h$, yield the result.
∎

##  Appendix A Approximation of integrodifferential operators of order $2s$ by
second order differential operators

To gain some intuition about what properties of the operator we would like to
preserve after regularization; and, in addition, to highlight the differences
between the case of constant $\eta$ (i.e., a fractional Laplacian) and a
variable one, here we inspect the approximation of the integrodifferential
operator

${\mathcal{I}}_{\varepsilon}[w](x)=\int_{B_{\varepsilon}}\frac{w(x+y)-2w(x)+w(x-y)}{|y|^{d+2s}}\eta\left(\frac{y}{|y|}\right)\mathop{}\\!\mathrm{d}y$

by considering its action on a quadratic, i.e., $w\in{\mathbb{P}}_{2}$. This
may help in justifying our choice, showing the stark difference between our
approach and that of [39, 40], as well as providing some intuition into
existing works that operate in a reverse way. For instance, [22, 58, 64]
approximate (local) second order elliptic differential operators by integral
operators like ${\mathcal{I}}_{\varepsilon}$. We emphasize, however, that to
obtain consistent, monotone, finite-difference schemes for arbitrary (local)
second order elliptic operators, one requires the use of wide stencils [49].

We begin by observing that, for $w\in{\mathbb{P}}_{2}$,

$w(x+y)-2w(x)+w(x-y)=D^{2}w(x):y\otimes y,$

where $D^{2}w(x)$ is the Hessian of $w$ at $x$, $\otimes$ is the outer
product, and $:$ is the Frobenius inner product. This means that

$\displaystyle{\mathcal{I}}_{\varepsilon}[w](x)$
$\displaystyle=\int_{B_{\varepsilon}}D^{2}w(x):y\otimes
y\frac{1}{|y|^{d+2s}}\eta\left(\frac{y}{|y|}\right)\mathop{}\\!\mathrm{d}y$
$\displaystyle=D^{2}w(x):\int_{B_{\varepsilon}}y\otimes
y\frac{1}{|y|^{d+2s}}\eta\left(\frac{y}{|y|}\right)\mathop{}\\!\mathrm{d}y={{\mathbf{A}}}:D^{2}w(x),$

where

${{\mathbf{A}}}=\int_{B_{\varepsilon}}y\otimes
y\frac{1}{|y|^{d+2s}}\eta\left(\frac{y}{|y|}\right)\mathop{}\\!\mathrm{d}y\in\mathbb{R}^{d\times
d},\qquad{{\mathbf{A}}}^{\intercal}={{\mathbf{A}}}.$

Consider now the particular case of $\eta\equiv 1$. The change of variables
$y=\varepsilon z$ with $z\in B_{1}$ shows that

${{\mathbf{A}}}=\varepsilon^{2(1-s)}\int_{B_{1}}z\otimes
z\frac{1}{|z|^{d+2s}}\mathop{}\\!\mathrm{d}z=\frac{\varepsilon^{2(1-s)}\omega_{d}}{2d(1-s)}{{\mathbf{I}}},$

where ${{\mathbf{I}}}$ is the identity matrix, and
$\omega_{d}=|{\mathbb{S}}^{d-1}|$. Therefore,

${\mathcal{I}}_{\varepsilon}[w](x)=\frac{\varepsilon^{2(1-s)}\omega_{d}}{2d(1-s)}\Delta
w(x),\qquad\forall w\in{\mathbb{P}}_{2}.$

Assume now that the coefficient $\eta$ is not constant, but even, i.e.,
$\eta(z)=\eta(-z)$. Notice that this assumption is included in the class
${\mathcal{C}}(\lambda,\Lambda)$ of Definition 1. The change of variables
$y=\varepsilon z$ implies that

${{\mathbf{A}}}=\varepsilon^{2(1-s)}\int_{B_{1}}z\otimes
z\frac{1}{|z|^{d+2s}}\eta\left(\frac{z}{|z|}\right)\mathop{}\\!\mathrm{d}z.$

We see that ${{\mathbf{A}}}$ is symmetric positive definite. Let us use polar
coordinates to obtain that

$\displaystyle{{\mathbf{A}}}$
$\displaystyle=\varepsilon^{2(1-s)}\int_{{\mathbb{S}}^{d-1}}\theta\otimes\theta\eta(\theta)\int_{0}^{1}\frac{r^{2}}{r^{d+2s}}r^{d-1}\mathop{}\\!\mathrm{d}r\mathop{}\\!\mathrm{d}\theta$
$\displaystyle=\frac{\varepsilon^{2(1-s)}}{2(1-s)}\int_{{\mathbb{S}}^{d-1}}\theta\otimes\theta\eta(\theta)\mathop{}\\!\mathrm{d}\theta=\frac{\varepsilon^{2(1-s)}}{2(1-s)}{{\mathbf{A}}}_{0},$

where, since $\eta$ is even, ${{\mathbf{A}}}_{0}$ is diagonal but anisotropic.
In this case

(40)
${\mathcal{I}}_{\varepsilon}[w](x)=\frac{\varepsilon^{2(1-s)}}{2(1-s)}{{\mathbf{A}}}_{0}:D^{2}w(x)=\frac{\varepsilon^{2(1-s)}}{2(1-s)}\nabla{\cdot}({{\mathbf{A}}}_{0}\nabla
w(x)),\qquad\forall w\in{\mathbb{P}}_{2}.$

Let us, finally, consider the general case of a coefficient of the form
$\eta(x,\tfrac{z}{|z|})$. In this case

${{\mathbf{A}}}_{0}={{\mathbf{A}}}_{0}(x)=\int_{{\mathbb{S}}^{d-1}}\theta\otimes\theta\eta(x,\theta)\mathop{}\\!\mathrm{d}\theta,$

and

(41)
${\mathcal{I}}_{\varepsilon}[w](x)=\frac{\varepsilon^{2(1-s)}}{2(1-s)}{{\mathbf{A}}}_{0}(x):D^{2}w(x),\qquad\forall
w\in{\mathbb{P}}_{2}.$

######  Remark 54 (anisotropy).

As we have already mentioned, it is not easy to construct approximations of
second order anisotropic operators, like those in (40) and (41), that are
monotone. In particular, (41) requires wide stencils.

Let us mention now what are the implications of the previous considerations
for the approximation of the operator ${\mathcal{L}}_{\eta}$, introduced in
(1). First, owing to the symmetry of the coefficient $\eta$, which is part of
Definition 1, we have (5), with ${\mathcal{K}}(r)=r^{-d-2s}$. Next, if we
assume that $w\in{\mathbb{P}}_{2}$, we have

$\displaystyle{\mathcal{L}}_{\eta}[w](x)$
$\displaystyle=\frac{1}{2}\int_{\mathbb{R}^{d}\setminus
B_{\varepsilon}}\left(2w(x)-w(x+y)-w(x-y)\right)\eta\left(\frac{y}{|y|}\right){\mathcal{K}}_{\varepsilon}(|y|)\mathop{}\\!\mathrm{d}y$
$\displaystyle-\frac{1}{2}{\mathcal{I}}_{\varepsilon}[w](x),$

where the smooth kernel ${\mathcal{K}}_{\varepsilon}$ was introduced in (7),
and

${\mathcal{I}}_{\varepsilon}[w](x)=\int_{B_{\varepsilon}}\frac{w(x+y)-2w(x)+w(x-y)}{|y|^{d+2s}}\eta\left(\frac{y}{|y|}\right)\mathop{}\\!\mathrm{d}y.$

Using (40) we then conclude that, if $w\in{\mathbb{P}}_{2}$,

$\displaystyle{\mathcal{L}}_{\eta}[w](x)$
$\displaystyle=\frac{1}{2}\int_{\mathbb{R}^{d}\setminus
B_{\varepsilon}}\left(2w(x)-w(x+y)-w(x-y)\right)\eta\left(\frac{y}{|y|}\right){\mathcal{K}}_{\varepsilon}(|y|)\mathop{}\\!\mathrm{d}y$
$\displaystyle-\frac{\varepsilon^{2(1-s)}}{1-s}\nabla{\cdot}({{\mathbf{A}}}_{0}\nabla
w(x)).$

In summary, by regularizing our operator, we have replaced it by a (local)
second order differential operator in divergence form, with a small
coefficient, and an operator of order zero. The monotonicity properties can
then be driven by the zero order operator without imposing any mesh
restrictions.

Finally, expression (5) motivates the consistency for quadratics. The symmetry
of the kernel transforms the difference inside the integral into a second
difference. A quadratic is the highest order polynomial for which a second
difference gives exactly the value of the (directional) derivative. Notice
that the approach from [39, 40], that merely truncates the kernel, does not
preserve this property.

## References

  * [1] Gabriel Acosta and Juan Pablo Borthagaray.  A fractional Laplace equation: regularity of solutions and finite element approximations.  SIAM J. Numer. Anal., 55(2):472–495, 2017.
  * [2] J. P. Agnelli, U. Kaufmann, and J. D. Rossi.  Numerical approximation of equations involving minimal/maximal operators by successive solution of obstacle problems.  J. Comput. Appl. Math., 342:133–146, 2018.
  * [3] Mark Ainsworth and Zhiping Mao.  Analysis and approximation of a fractional Cahn-Hilliard equation.  SIAM J. Numer. Anal., 55(4):1689–1718, 2017.
  * [4] Goro Akagi, Giulio Schimperna, and Antonio Segatti.  Fractional Cahn-Hilliard, Allen-Cahn and porous medium equations.  J. Differential Equations, 261(6):2935–2985, 2016.
  * [5] C. Baiocchi.  Estimations d’erreur dans $L^{\infty}$ pour les inéquations à obstacle.  In Mathematical aspects of finite element methods (Proc. Conf., Consiglio Naz. delle Ricerche (C.N.R.), Rome, 1975),, Lecture Notes in Math., Vol. 606,, pages 27–34. ,, 1977.
  * [6] Begoña Barrios, Alessio Figalli, and Xavier Ros-Oton.  Global regularity for the free boundary in the obstacle problem for the fractional Laplacian.  Amer. J. Math., 140(2):415–447, 2018.
  * [7] Jean Bertoin.  Lévy processes.  Cambridge tracts in mathematics n$\,{}^{\circ}$121\. Cambridge University Press, 1996.
  * [8] Imran H. Biswas, Indranil Chowdhury, and Espen R. Jakobsen.  On the rate of convergence for monotone numerical schemes for nonlocal Isaacs equations.  SIAM J. Numer. Anal., 57(2):799–827, 2019.
  * [9] Imran H. Biswas, Espen R. Jakobsen, and Kenneth H. Karlsen.  Error estimates for finite difference-quadrature schemes for a class of nonlocal Bellman equations with variable diffusion.  In Stochastic analysis and partial differential equations, volume 429 of Contemp. Math., pages 19–31. Amer. Math. Soc., Providence, RI, 2007.
  * [10] Imran H. Biswas, Espen R. Jakobsen, and Kenneth H. Karlsen.  Error estimates for a class of finite difference-quadrature schemes for fully nonlinear degenerate parabolic integro-PDEs.  J. Hyperbolic Differ. Equ., 5(1):187–219, 2008.
  * [11] Pablo Blanc, Juan P. Pinasco, and Julio D. Rossi.  Obstacle problems and maximal operators.  Adv. Nonlinear Stud., 16(2):355–362, 2016.
  * [12] Andrea Bonito, Juan Pablo Borthagaray, Ricardo H. Nochetto, Enrique Otárola, and Abner J. Salgado.  Numerical methods for fractional diffusion.  Comput. Vis. Sci., 19(5-6):19–46, 2018.
  * [13] Andrea Bonito, Wenyu Lei, and Joseph E. Pasciak.  Numerical approximation of the integral fractional Laplacian.  Numer. Math., 142(2):235–278, 2019.
  * [14] Andrea Bonito and Joseph E. Pasciak.  Numerical approximation of fractional powers of elliptic operators.  Math. Comp., 84(295):2083–2110, 2015.
  * [15] J.P. Borthagaray, D. Leykekhman, and R.H. Nochetto.  Local energy estimates for the fractional Laplacian.  SIAM J. Numer. Anal., 59(4):1918–1947, 2021.
  * [16] Juan Pablo Borthagaray and Ricardo H Nochetto.  Besov regularity for the Dirichlet integral fractional Laplacian in Lipschitz domains.  Journal of Functional Analysis, 284(6):109829, 2023.
  * [17] Juan Pablo Borthagaray, Ricardo H. Nochetto, and Abner J. Salgado.  Weighted Sobolev regularity and rate of approximation of the obstacle problem for the integral fractional Laplacian.  Math. Models Methods Appl. Sci., 29(14):2679–2717, 2019.
  * [18] Luis Caffarelli and Fernando Charro.  On a fractional Monge-Ampère operator.  Ann. PDE, 1(1):Art. 4, 47, 2015.
  * [19] Luis Caffarelli, Xavier Ros-Oton, and Joaquim Serra.  Obstacle problems for integro-differential operators: regularity of solutions and free boundaries.  Invent. Math., 208(3):1155–1211, 2017.
  * [20] Luis Caffarelli and Luis Silvestre.  An extension problem related to the fractional Laplacian.  Comm. Partial Differential Equations, 32(7-9):1245–1260, 2007.
  * [21] Luis Caffarelli and Luis Silvestre.  Regularity theory for fully nonlinear integro-differential equations.  Comm. Pure Appl. Math., 62(5):597–638, 2009.
  * [22] Luis Caffarelli and Luis Silvestre.  Smooth approximations of solutions to nonconvex fully nonlinear elliptic equations.  In Nonlinear partial differential equations and related topics, volume 229 of Amer. Math. Soc. Transl. Ser. 2, pages 67–85. Amer. Math. Soc., Providence, RI, 2010.
  * [23] Luis Caffarelli and Luis Silvestre.  The Evans-Krylov theorem for nonlocal fully nonlinear equations.  Ann. of Math. (2), 174(2):1163–1187, 2011.
  * [24] Luis A. Caffarelli, Sandro Salsa, and Luis Silvestre.  Regularity estimates for the solution and the free boundary of the obstacle problem for the fractional Laplacian.  Invent. Math., 171(2):425–461, 2008.
  * [25] Fabio Camilli and Espen R. Jakobsen.  A finite element like scheme for integro-partial differential Hamilton-Jacobi-Bellman equations.  SIAM J. Numer. Anal., 47(4):2407–2431, 2009.
  * [26] B. Carmichael, H Babahosseini, S. N. Mahmoodi, and M. Agah.  The fractional viscoelastic response of human breast tissue cells.  Physical Biology, 12(4):046001, May 2015.
  * [27] Indranil Chowdhury and Espen R. Jakobsen.  Precise error bounds for numerical approximations of fractional HJB equations.  arXiv:2308.16434, 2023.
  * [28] Indranil Chowdhury, Espen Robstad Jakobsen, and Robin Østern Lien.  Discretization of fractional fully nonlinear equations by powers of discrete Laplacians.  arXiv:2401.09926, 2024.
  * [29] P. G. Ciarlet.  Basic error estimates for elliptic problems.  In Handbook of numerical analysis, Vol. II, volume II of Handb. Numer. Anal., pages 17–351. North-Holland, Amsterdam, 1991.
  * [30] Rama Cont and Peter Tankov.  Financial Modelling with Jump Processes, volume 101.  01 2006.
  * [31] Arturo de Pablo, Fernando Quirós, Ana Rodríguez, and Juan Luis Vázquez.  A general fractional porous medium equation.  Comm. Pure Appl. Math., 65(9):1242–1284, 2012.
  * [32] Marta D’Elia, Qiang Du, Christian Glusa, Max Gunzburger, Xiaochuan Tian, and Zhi Zhou.  Numerical methods for nonlocal and fractional models.  Acta Numer., 29:1–124, 2020.
  * [33] Vincent J. Ervin, Norbert Heuer, and John Paul Roop.  Numerical approximation of a time dependent, nonlinear, space-fractional diffusion equation.  SIAM J. Numer. Anal., 45(2):572–591, 2007.
  * [34] Vincent J. Ervin and John Paul Roop.  Variational formulation for the stationary fractional advection dispersion equation.  Numer. Methods Partial Differential Equations, 22(3):558–576, 2006\.
  * [35] Lawrence C. Evans and Avner Friedman.  Optimal stochastic switching and the Dirichlet problem for the Bellman equation.  Trans. Amer. Math. Soc., 253:365–389, 1979.
  * [36] Xavier Fernádez-Real and Xavier Ros-Oton.  Integro-Differential elliptic equations.  Progress in Mathematics. Birkhauser, 2024.
  * [37] R. K. Getoor.  First passage times for symmetric stable processes in space.  Trans. Amer. Math. Soc., 101:75–90, 1961.
  * [38] Gerd Grubb.  Fractional Laplacians on domains, a development of Hörmander’s theory of $\mu$-transmission pseudodifferential operators.  Adv. Math., 268:478–528, 2015.
  * [39] Rubing Han and Shuonan Wu.  A monotone discretization for integral fractional Laplacian on bounded Lipschitz domains: pointwise error estimates under Hölder regularity.  SIAM J. Numer. Anal., 60(6):3052–3077, 2022.
  * [40] Rubing Han, Shuonan Wu, and Hao Zhou.  A monotone discretization for the fractional obstacle problem and its improved policy iteration.  Fractal and Fractional, 7(8):634, 2023.
  * [41] Yanghong Huang and Adam Oberman.  Numerical methods for the fractional Laplacian: a finite difference–quadrature approach.  SIAM J. Numer. Anal., 52(6):3056–3084, 2014.
  * [42] M. Ilic, F. Liu, I. Turner, and V. Anh.  Numerical approximation of a fractional-in-space diffusion equation. I.  Fract. Calc. Appl. Anal., 8(3):323–341, 2005.
  * [43] M. Ilic, F. Liu, I. Turner, and V. Anh.  Numerical approximation of a fractional-in-space diffusion equation. II. With nonhomogeneous boundary conditions.  Fract. Calc. Appl. Anal., 9(4):333–349, 2006.
  * [44] Tianling Jin and Jingang Xiong.  Schauder estimates for nonlocal fully nonlinear equations.  Ann. Inst. H. Poincaré C Anal. Non Linéaire, 33(5):1375–1407, 2016.
  * [45] Janne Korvenpää, Tuomo Kuusi, and Giampiero Palatucci.  The obstacle problem for nonlinear integro-differential operators.  Calc. Var. Partial Differential Equations, 55(3):Art. 63, 29, 2016\.
  * [46] Wenbo Li and Ricardo H. Nochetto.  Optimal pointwise error estimates for two-scale methods for the Monge-Ampère equation.  SIAM J. Numer. Anal., 56(3):1915–1941, 2018.
  * [47] P.-L. Lions and B. Mercier.  Approximation numérique des équations de Hamilton-Jacobi-Bellman.  RAIRO Anal. Numér., 14(4):369–393, 1980.
  * [48] Anna Lischke, Guofei Pang, Mamikon Gulian, and et al.  What is the fractional Laplacian? A comparative review with new results.  J. Comput. Phys., 404:109009, 62, 2020.
  * [49] Michael Neilan, Abner J. Salgado, and Wujun Zhang.  Numerical analysis of strongly nonlinear PDEs.  Acta Numer., 26:137–303, 2017.
  * [50] Michael Neilan, Abner J. Salgado, and Wujun Zhang.  Chapter 2 — The Monge-Ampère equation.  In Andrea Bonito and Ricardo H. Nochetto, editors, Geometric Partial Differential Equations - Part I, volume 21 of Handbook of Numerical Analysis, pages 105–219. Elsevier, 2020.
  * [51] Joachim Nitsche.  $L_{\infty}$-convergence of finite element approximations.  In Mathematical aspects of finite element methods (Proc. Conf., Consiglio Naz. delle Ricerche (C.N.R.), Rome, 1975),, Lecture Notes in Math., Vol. 606,, pages 261–274. ,, 1977.
  * [52] R. H. Nochetto, D. Ntogkas, and W. Zhang.  Two-scale method for the Monge-Ampère equation: convergence to the viscosity solution.  Math. Comp., 88(316):637–664, 2019.
  * [53] R. H. Nochetto, D. Ntogkas, and W. Zhang.  Two-scale method for the Monge-Ampère equation: pointwise error estimates.  IMA J. Numer. Anal., 39(3):1085–1109, 2019.
  * [54] Ricardo H. Nochetto and Dimitrios Ntogkas.  Convergent two-scale filtered scheme for the Monge-Ampère equation.  SIAM J. Sci. Comput., 41(2):B295–B319, 2019.
  * [55] Ricardo H. Nochetto, Enrique Otárola, and Abner J. Salgado.  Convergence rates for the classical, thin and fractional elliptic obstacle problems.  Philos. Trans. Roy. Soc. A, 373(2050):20140449, 14, 2015.
  * [56] Ricardo H. Nochetto, Enrique Otárola, and Abner J. Salgado.  A PDE approach to fractional diffusion in general domains: a priori error analysis.  Found. Comput. Math., 15(3):733–791, 2015.
  * [57] Ricardo H. Nochetto, Kunibert G. Siebert, and Andreas Veeser.  Fully localized a posteriori error estimators and barrier sets for contact problems.  SIAM J. Numer. Anal., 42(5):2118–2135, 2005.
  * [58] Ricardo H. Nochetto and Wujun Zhang.  Discrete ABP estimate and convergence rates for linear elliptic equations in non-divergence form.  Found. Comput. Math., 18(3):537–593, 2018.
  * [59] Christoph Reisinger and Yufei Zhang.  A penalty scheme and policy iteration for nonlocal HJB variational inequalities with monotone nonlinearities.  Comput. Math. Appl., 93:199–213, 2021.
  * [60] Xavier Ros-Oton.  Nonlocal elliptic equations in bounded domains: a survey.  Publ. Mat., 60(1):3–26, 2016.
  * [61] Xavier Ros-Oton and Joaquim Serra.  The Dirichlet problem for the fractional Laplacian: regularity up to the boundary.  J. Math. Pures Appl. (9), 101(3):275–302, 2014.
  * [62] Xavier Ros-Oton and Joaquim Serra.  Boundary regularity for fully nonlinear integro-differential equations.  Duke Math. J., 165(11):2079–2154, 2016.
  * [63] Xavier Ros-Oton and Joaquim Serra.  Regularity theory for general stable operators.  J. Differential Equations, 260(12):8675–8715, 2016.
  * [64] Abner J. Salgado and Wujun Zhang.  Finite element approximation of the Isaacs equation.  ESAIM Math. Model. Numer. Anal., 53(2):351–374, 2019.
  * [65] S. A. Silling.  Reformulation of elasticity theory for discontinuities and long-range forces.  J. Mech. Phys. Solids, 48(1):175–209, 2000.
  * [66] S. A. Silling, M. Epton, O. Weckner, J. Xu, and E. Askari.  Peridynamic states and constitutive modeling.  J. Elasticity, 88(2):151–184, 2007.
  * [67] Luis Silvestre.  Regularity of the obstacle problem for a fractional power of the Laplace operator.  Comm. Pure Appl. Math., 60(1):67–112, 2007.
  * [68] Bradley E. Treeby and B. T. Cox.  Modeling power law absorption and dispersion for acoustic propagation using the fractional Laplacian.  The Journal of the Acoustical Society of America, 127(5):2741–2748, 2010.
  * [69] M. I. Višik and G. I. Èskin.  Convolution equations in a bounded region.  Uspehi Mat. Nauk, (no. 3 (123)):89–152, 1965.
  * [70] Masakazu Yamamoto.  Asymptotic expansion of solutions to the dissipative equation with fractional Laplacian.  SIAM J. Math. Anal., 44(6):3786–3805, 2012.