# Quasi-linear fractional-order operators in Lipschitz domains

Juan Pablo Borthagaray  Instituto de Matemática y Estadística “Rafael
Laguardia”, Facultad de Ingeniería, Universidad de la República, Montevideo,
Uruguay<EMAIL_ADDRESS>,  Wenbo Li  Institute of Computational
Mathematics and Scientific/Engineering Computing of the Chinese Academy of
Sciences, Beijing 100190 China<EMAIL_ADDRESS>and  Ricardo H.
Nochetto  Department of Mathematics and Institute for Physical Science and
Technology, University of Maryland, College Park, MD 20742, USA
<EMAIL_ADDRESS>
###### Abstract.

We prove Besov boundary regularity for solutions of the homogeneous Dirichlet
problem for fractional-order quasi-linear operators with variable coefficients
on Lipschitz domains $\Omega$ of $\mathbb{R}^{d}$. Our estimates are
consistent with the boundary behavior of solutions on smooth domains and apply
to fractional $p$-Laplacians and operators with finite horizon. The proof
exploits the underlying variational structure and uses a new and flexible
local translation operator. We further apply these regularity estimates to
derive novel error estimates for finite element approximations of fractional
$p$-Laplacians and present several simulations that reveal the boundary
behavior of solutions.

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

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

##  1\. Introduction

In recent years, fractional-order and, more generally, nonlocal operators have
received a great deal of attention in applied sciences and engineering. This
is mainly because such operators arise in jump processes modeling the
ubiquitous phenomenon of anomalous diffusion [41]. In this vein, the
fractional Laplacian, an outstanding nonlocal operator, arises as a limit of a
long-jump random walk [46]. Among other applications of nonlocal operators, we
mention finance [21, 38], ground-water solute transport [8], and biological
systems with binding, crowding, or trapping, such as electrodiffusion of ions
within nerve cells [36, 37].

For problems with a variational structure, finite element methods provide the
best approximation in the energy norm, and are amenable to an analysis with
low regularity conditions. In our setting, the latter is fundamental because
solutions of fractional-order problems generically develop algebraic boundary
layers. Solution regularity estimates in the Sobolev scale are a key
ingredient to prove a priori convergence rates for the finite element
discretization of such problems.

However, most progress in that direction and most computational studies have
been limited to either linear or semi-linear problems. This paper deals with
fractional-order quasi-linear operators. We prove elliptic regularity
estimates up to the boundary of the domain, which is only assumed to be
bounded and Lipschitz. The model operator we consider is the so-called
$(p,s)$-fractional Laplacian ($s\in(0,1)$, $p\in(1,\infty)$), but our theory
is also valid for a broader class of operators, including operators with
finite horizon. In this regard, we remark that our regularity estimates for
finite-horizon operators are even new for linear problems. As an application
of our regularity estimates, we consider direct finite element discretization
of the problems under study and prove convergence rates in the energy norm.

Let us make precise the problem setting in this paper. Let
$\Omega\subset{\mathbb{R}^{d}}$ ($d\geq 1$) be a bounded, Lipschitz domain,
$s\in(0,1)$, and $p\in(1,\infty)$. We consider energy functionals with domain
the fractional-order Sobolev space $\widetilde{W}^{s}_{p}(\Omega)$, namely
functions in $W^{s}_{p}(\mathbb{R}^{d})$ that vanish in
$\Omega^{c}:=\mathbb{R}^{d}\setminus\overline{\Omega}$. More precisely, for a
given function
$G\colon\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}\to(0,\infty)$, with
$(x,y,\rho)\mapsto G(x,y,\rho)$, and
$f\in(\widetilde{W}^{s}_{p}(\Omega))^{\prime}$, we are interested in
minimizers of the energy

(1.1)
$\mathcal{F}(u):=\iint_{Q_{\Omega}}G\left(x,y,\frac{u(x)-u(y)}{|x-y|^{s}}\right)\frac{1}{|x-y|^{d}}\,dydx-\langle
f,u\rangle.$

Above, $\langle\cdot,\cdot\rangle$ stands for the duality pairing between
$(\widetilde{W}^{s}_{p}(\Omega))^{\prime}$ and $\widetilde{W}^{s}_{p}(\Omega)$
and

$Q_{\Omega}:=({\mathbb{R}^{d}}\times{\mathbb{R}^{d}})\setminus(\Omega^{c}\times\Omega^{c}).$

Specific requirements on $G$ are listed in Hypothesis 2.1 below. The Gateaux
differential of $\mathcal{F}$ at $u$ is given by
$\mathcal{A}\colon\widetilde{W}^{s}_{p}(\Omega)\to(\widetilde{W}^{s}_{p}(\Omega))^{\prime}$,

(1.2)
$\mathcal{A}u(x):=\int_{{\mathbb{R}^{d}}}\left[G_{\rho}\left(x,y,\frac{u(x)-u(y)}{|x-y|^{s}}\right)-G_{\rho}\left(y,x,\frac{u(y)-u(x)}{|x-y|^{s}}\right)\right]\frac{1}{|x-y|^{d+s}}\,dy.$

For the moment, let us assume that $G$ satisfies the relation
$G(x,y,\rho)=G(y,x,-\rho)$ for a.e. $x,y,\rho$. While this assumption allows
us to write the minimization problem in a strong form in a concise fashion, it
is not necessary for our theoretical results. Under this additional condition,
we have $G_{\rho}(x,y,\rho)=-G_{\rho}(y,x,-\rho)$ for a.e. $x,y,\rho$ and we
can write

(1.3)
$\mathcal{A}u(x):=2\int_{{\mathbb{R}^{d}}}G_{\rho}\left(x,y,\frac{u(x)-u(y)}{|x-y|^{s}}\right)\frac{1}{|x-y|^{d+s}}\,dy.$

Minimizers of (1.1) are weak solutions of the homogeneous Dirichlet problem
for the operator $\mathcal{A}$:

(1.4) $\left\\{\begin{array}[]{rl}\mathcal{A}u=f&\mbox{in }\Omega,\\\
u=0&\mbox{in }\Omega^{c}.\end{array}\right.$

We assume standard hypotheses on $G$ in order to apply the direct method in
the calculus of variations. As a prototypical example, we consider
$G(x,y,\rho)=\frac{C_{d,s,p}}{2p}|\rho|^{p}$ with $C_{d,s,p}$ defined below.
Then, $G_{\rho}(x,y,\rho)=\frac{C_{d,s,p}}{2}|\rho|^{p-2}\rho$, and

(1.5)
$\mathcal{A}u(x)=(-\Delta)^{s}_{p}u(x):=C_{d,s,p}\int_{\mathbb{R}^{d}}\frac{|u(x)-u(y)|^{p-2}(u(x)-u(y))}{|x-y|^{d+sp}}\,dy$

is the so-called fractional $(p,s)$-Laplacian (or fractional $p$-Laplacian of
order $s$). We define the normalizing constant $C_{d,s,p}$ as

(1.6)
$C_{d,s,p}=\frac{s(1-s)p\;\Gamma(\frac{ps+d}{2})\;2^{2s-2}}{\pi^{\frac{d-1}{2}}\Gamma(\frac{(p-2)s+3}{2})\Gamma(2-s)}.$

This choice is somewhat arbitrary, but for $p=2$ it allows us to recover the
integral fractional Laplacian, which is the pseudodifferential operator with
symbol $|\xi|^{2s}$. Moreover, for every smooth function $v\in
C_{c}^{\infty}({\mathbb{R}^{d}})$ we have the asymptotic behaviors [14, 39,
22]

(1.7) $\lim_{s\to 0^{+}}(-\Delta)^{s}_{p}v=|v|^{p-2}v,\quad\lim_{s\to
1^{-}}(-\Delta)^{s}_{p}v=-\nabla\cdot(|\nabla v|^{p-2}\nabla v).$

We also point out that the integral in (1.5) needs to be understood in the
principal value sense if $s\geq 1-\frac{1}{p}$.

Another way to write the operator in (1.5) is

(1.8)
$(-\Delta)^{s}_{p}u(x)=2\int_{\mathbb{R}^{d}}\left(\frac{|u(x)-u(y)|}{|x-y|^{s}}\right)^{p-2}\frac{(u(x)-u(y))}{|x-y|^{d+2s}}\,dy,$

which suggests that, heuristically, one can understand the fractional
$(p,s)$-Laplacian as a weighted fractional Laplacian of order $s$, with a
weight $\left(\frac{|u(x)-u(y)|}{|x-y|^{s}}\right)^{p-2}$. This is analogous
to the local case, for which the $p$-Laplacian
$(-\Delta)_{p}u:=\mbox{div}(|\nabla u|^{p-2}\nabla u)$ can be regarded as a
Laplacian with weight $|\nabla u|^{p-2}$. The Dirichlet problem for the local
$p$-Laplacian arises in a number of models of physical processes, including
non-Newtonian fluids [5], turbulent flows in porous media [25], and global
climate modeling [26]. We refer to [7] for a historical account and other
applications of this operator, and to [6, 18, 31] for its numerical treatment.

The representation (1.8) also shows that the operator (1.5) corresponds to a
degenerate diffusion if $p>2$ and to a singular one if $p<2$. We refer to [17]
for several motivations for considering nonlinear operators like (1.5), to
[42] for a thorough discussion about existence and regularity results for
problems driven by the fractional $(p,s)$-Laplacian, and to [23] for a
monotone finite difference scheme with consistency error estimates for $C^{4}$
functions and applications to the Cauchy problem for such an operator.

Depending on whether the resulting operator $\mathcal{A}$ in (1.2) is
degenerate or singular, our regularity estimates are somewhat different from
one another. To make the point clear, let us focus on the case of the
$(p,s)$-Laplacian, although we emphasize that our estimates are valid for more
general operators (cf. Theorem 3.1 below). We derive shift estimates in Besov
norms of the form

(1.9)
$\displaystyle\|u\|_{\dot{B}^{s+\frac{1}{p}}_{p,\infty}(\Omega)}\lesssim\|f\|_{B^{-s+\frac{1}{p^{\prime}}}_{p^{\prime},1}(\Omega)}^{\frac{1}{p-1}},$
$\displaystyle\mbox{ if }p\geq 2,$ (1.10)
$\displaystyle\|u\|_{\dot{B}^{s+\frac{1}{2}}_{p,\infty}(\Omega)}\lesssim\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}^{\frac{2-p}{p-1}}\|f\|_{B^{-s+\frac{1}{2}}_{p^{\prime},1}(\Omega)},$
$\displaystyle\mbox{ if }p\leq 2,$

where $p^{\prime}=\frac{p}{p-1}$; these Besov estimates extend classical ones
[43, Theorems 2 and 2’] to the fractional setting. To check optimality, we
consider the prototypical $1d$-function $v(x)=x_{+}^{s}$ which mimics the
boundary behavior of solutions of (1.4) for the operator (1.5). A simple
calculation using second differences shows that $v\in
B^{s+1/p}_{p,\infty}(\Omega)$ for all $p\in(1,\infty)$, which revals that
(1.9) is optimal while (1.10) is suboptimal. Moreover, by a simple embedding
argument, (1.9) and (1.10) give rise to Sobolev regularity estimates. Estimate
(1.9) for $p=2$ turns out to be consistent with well-known optimal regularity
for solutions to the Dirichlet problem for the integral fractional Laplacian
on smooth domains, cf. [1, 33, 47]. Importantly, our estimates are valid for
Lipschitz domains and in that sense generalize the ones derived in [4, 11] to
a quasi-linear setting. Additionally, our estimates are valid under general
conditions on the function $G$. In this vein, we point out to [30] where, for
a class of nonlinear operators related to the ones in this work, analysis is
performed in fractional-order Orlicz-Sobolev spaces and Hölder regularity
estimates are derived for Dirichlet problems on bounded $C^{1,1}$ domains.

The paper is organized as follows. Section 2 collects preliminary material
about function spaces and Lipschitz domains, introduces a flexible local
translation operator that plays an instrumental role in our derivation of
regularity estimates, specifies the assumptions we require on the energy, and
discusses the use of localized translations in the proof of regularity of
energy minimizers. Section 3 contains the core of the paper, and studies the
regularity of solutions through the derivation of suitable energy bounds. It
also discusses the extension of the technique to operators with finite horizon
and truncated Laplacians in the linear setting. Section 4 proposes and
analyzes a finite element discretization of problems of the form (1.4), and
exploits (1.9) and (1.10) to prove error bounds for all $p\in(1,\infty)$.
Finally, Section 5 exhibits some numerical experiments that explore the
accuracy of this approach and the boundary behavior of solutions to the
$(p,s)$-Laplacian (1.5) and linear truncated Laplacians.

##  2\. Notation and assumptions

This section establishes the notation and collects some preliminary results.
We provide some discussion on function spaces and Lipschitz domains. We
analyze function space characterizations by means of translation operators,
discuss the relation between these translations and the regularity of
minimizers of (1.1), and introduce a suitable localized translation operator
to derive regularity estimates. Finally, we make explicit assumptions on the
energy, discuss some of their consequences, and comment on how they apply to
the model operator (1.5).

###  2.1. Sobolev and Besov spaces

Here, we briefly review some important facts about Sobolev and Besov spaces.
We follow the notation from [11] and refer to that work for further details.

Given $\sigma\in(0,1)$ and $p\in[1,\infty)$, we consider the zero-extension
Sobolev space

$\widetilde{W}^{\sigma}_{p}(\Omega):=\big{\\{}v\in
W^{\sigma}_{p}(\mathbb{R}^{d})\colon\textrm{supp~{}}v\subset\overline{\Omega}\big{\\}};$

this is a Banach space furnished with the norm

$\|v\|_{\widetilde{W}^{\sigma}_{p}(\Omega)}:=|v|_{W^{\sigma}_{p}({\mathbb{R}^{d}})}=\left(\frac{C_{d,s,p}}{2}\iint_{{\mathbb{R}^{d}}\times{\mathbb{R}^{d}}}\frac{|v(x)-v(y)|^{p}}{|x-y|^{d+\sigma
p}}\,dx\,dy\right)^{1/p}.$

Because functions in $\widetilde{W}^{\sigma}_{p}(\Omega)$ vanish in
$\Omega^{c}$, the integrand above vanishes on $\Omega^{c}\times\Omega^{c}$ and
one can effectively compute the integral over
$Q_{\Omega}=({\mathbb{R}^{d}}\times{\mathbb{R}^{d}})\setminus(\Omega^{c}\times\Omega^{c})$.

We define Besov spaces by real interpolation. Given a pair of compatible
Banach spaces $(X_{0},X_{1})$, $u\in X_{0}+X_{1}$, $t>0$, and
$p\in[1,\infty)$, we consider the $K$-functional

(2.1) $K(t,u):=\inf\left\\{\|u_{0}\|_{X_{0}}+t\|u_{1}\|_{X_{1}}\colon
u=u_{0}+u_{1},\ u_{0}\in X_{0},\ u_{1}\in X_{1}\right\\}.$

For $\theta\in(0,1)$ and $q\in[1,\infty]$, we define the interpolation spaces

$\big{[}X_{0},X_{1}\big{]}_{\theta,q}:=\\{u\in
X_{0}+X_{1}\colon\|u\|_{(X_{0},X_{1})_{\theta,q}}<\infty\\},$

where

(2.2)
$\|u\|_{[X_{0},X_{1}]_{\theta,q}}:=\begin{cases}\Big{[}q\theta(1-\theta)\int_{0}^{\infty}t^{-(1+\theta
q)}|K(t,u)|^{q}\,dt\Big{]}^{1/q}&\mbox{if }1\leq q<\infty,\\\ \sup_{t>0}\
t^{-\theta}|K(t,u)|&\mbox{if }q=\infty.\end{cases}$

The normalization factor $q\theta(1-\theta)$ in (2.2) guarantees the correct
scalings in the limits $\theta\to 0$, $\theta\to 1$ and $q\to\infty$ for norm
continuity; see [40, Appendix B] for a detailed proof in the case $p=2$.
Because we are interested in spaces with differentiability order between zero
and two, we let $X_{0}:=L^{p}(\Omega)$, $X_{1}:=W^{2}_{p}(\Omega)$,
$\sigma\in(0,2)$ and $q\in[1,\infty]$ to define the Besov spaces

$B^{\sigma}_{p,q}(\Omega):=\big{[}L^{p}(\Omega),W^{2}_{p}(\Omega)\big{]}_{\sigma/2,q},\quad\dot{B}^{\sigma}_{p,q}(\Omega):=\\{v\in
B^{\sigma}_{p,q}(\Omega):\textrm{supp~{}}v\subset\overline{\Omega}\\}.$

By reiteration, we have the following result regarding interpolation of Besov
spaces (cf. [9, Theorem 6.4.5]): given $\sigma_{0}\neq\sigma_{1}$, $1\leq
p,q_{0},q_{1},r\leq\infty$ and $0<\theta<1$,

(2.3)
$\left(B^{\sigma_{0}}_{p,q_{0}}(\Omega),B^{\sigma_{1}}_{p,q_{1}}(\Omega)\right)_{\theta,r}=B^{\sigma}_{p,r}(\Omega),\quad\mbox{where
}\sigma=(1-\theta)\sigma_{0}+\theta\sigma_{1}.$

Importantly, we have $B^{\sigma}_{p,p}(\Omega)=W^{\sigma}_{p}(\Omega)$ for all
$p\in[1,\infty)$, $\sigma\in(0,2)\setminus\\{1\\}$. In the case $\sigma=1$, we
only have the equality $B^{1}_{2,2}(\Omega)=H^{\sigma}(\Omega)$, while
$B^{1}_{p,p}(\Omega)\subset W^{1}_{p}(\Omega)$ if $p<2$ and
$B^{1}_{p,p}(\Omega)\supset W^{1}_{p}(\Omega)$ if $p>2$, cf. [3, §7.67].
Moreover, we have the following inclusions between Besov spaces on bounded
Lipschitz domains [45, §3.2.4, §3.3.1]:

$\begin{array}[]{llll}B^{\sigma}_{p,q_{0}}(\Omega)\subset
B^{\sigma}_{p,q_{1}}(\Omega),&\mbox{ if }\sigma>0,&1\leq p\leq\infty,&1\leq
q_{0}\leq q_{1}\leq\infty;\\\ B^{\sigma_{1}}_{p,q_{1}}(\Omega)\subset
B^{\sigma_{0}}_{p,q_{0}}(\Omega)&\mbox{ if }0<\sigma_{0}<\sigma_{1},&1\leq
p\leq\infty,&1\leq q_{0},q_{1}\leq\infty.\end{array}$

We are interested in making precise the statement about inclusion of a higher-
order Besov space with integrability index $p\in[1,\infty)$ and second
parameter $q=\infty$ into a lower-order Sobolev space with the same
integrability index. Concretely, the next lemma shows the scaling of the
continuity constant.

######  Lemma 2.1 (embedding).

Let $\Omega\subset{\mathbb{R}^{d}}$ be a bounded Lipschitz domain,
$p\in[1,\infty)$, $\sigma\in(0,2)\setminus\\{1\\}$, and
$\varepsilon\in(0,2-\sigma)$. Then,
$B^{\sigma+\varepsilon}_{p,\infty}(\Omega)\subset W^{\sigma}_{p}(\Omega)$ with

(2.4)
$\|v\|_{W^{\sigma}_{p}(\Omega)}\lesssim\left(\frac{\sigma(2-\sigma)}{\varepsilon}\right)^{\frac{1}{p}}\|v\|_{B^{\sigma+\varepsilon}_{p,\infty}(\Omega)}\quad\forall\,v\in
B^{\sigma+\varepsilon}_{p,\infty}(\Omega).$

###### Proof.

We exploit the characterization of Besov and fractional-order Sobolev spaces
as interpolation spaces between integer-order Sobolev spaces. More precisely,
if the $K$-functional corresponds to interpolation between the spaces
$X_{0}=L^{p}(\Omega)$ and $X_{1}=W^{2}_{p}(\Omega)$, we recall the norm
definitions

$\|v\|_{B^{\sigma+\varepsilon}_{p,\infty}(\Omega)}=\sup_{t>0}\,\Big{(}t^{-\frac{\sigma+\varepsilon}{2}}|K(t,v)|\Big{)},$

and

$\|v\|_{W^{\sigma}_{p}(\Omega)}^{p}=\frac{p\sigma(2-\sigma)}{4}\int_{0}^{\infty}t^{-1-\frac{\sigma
p}{2}}|K(t,v)|^{p}dt$

for $\sigma\in(0,2)\setminus\\{1\\}$, according to the remark following (2.3).
We split the integral above as the sum of the integrals between $0$ and $N$
and between $N$ and $\infty$, with $N>0$ to be chosen. A straightforward
calculation gives

$\int_{0}^{N}t^{-1-\frac{\sigma
p}{2}}|K(t,v)|^{p}dt\leq\sup_{t>0}\,t^{-\frac{(\sigma+\varepsilon)p}{2}}|K(t,v)|^{p}\,\int_{0}^{N}t^{-1+\frac{\varepsilon
p}{2}}dt=\frac{2N^{\frac{\varepsilon p}{2}}}{\varepsilon
p}\|v\|_{B^{\sigma+\varepsilon}_{p,\infty}(\Omega)}^{p}.$

Additionally, for any $v\in B^{\sigma+\varepsilon}_{p,\infty}(\Omega)$ and
$t\geq 0$, we choose the trivial decomposition $v=v+0$ in (2.1) to obtain

$|K(t,v)|\leq\|v\|_{L^{p}(\Omega)}\leq\|v\|_{W^{\sigma}_{p}(\Omega)}.$

This gives rise to

$\int_{N}^{\infty}t^{-1-\frac{\sigma
p}{2}}|K(t,v)|^{p}dt\leq\|v\|_{W^{\sigma}_{p}(\Omega)}^{p}\int_{N}^{\infty}t^{-1-\frac{\sigma
p}{2}}dt=\frac{2N^{-\frac{\sigma
p}{2}}}{p\sigma}\|v\|_{W^{\sigma}_{p}(\Omega)}^{p},$

and thus

$\|v\|_{W^{\sigma}_{p}(\Omega)}^{p}\leq\frac{\sigma(2-\sigma)N^{\frac{\varepsilon
p}{2}}}{2\varepsilon}\|v\|_{B^{\sigma+\varepsilon}_{p,\infty}(\Omega)}^{p}+\frac{(2-\sigma)N^{-\frac{\sigma
p}{2}}}{2}\|v\|_{W^{\sigma}_{p}(\Omega)}^{p}.$

It now suffices to fix $N$ sufficiently large so that
$\frac{(2-\sigma)N^{-\frac{\sigma p}{2}}}{2}<1$ and kick back the last term in
the right hand side above to arrive to the desired estimate (2.4). ∎

For $\sigma\in(0,1)$ and $p,q\in[1,\infty]$, we define

$B^{-\sigma}_{p,q}(\Omega):=\big{(}L^{p}(\Omega),W^{-1}_{p}(\Omega)\big{)}_{\sigma,q},$

and point out that, if $p,q\in(1,\infty]$ and
$p^{\prime},q^{\prime}\in[1,\infty)$ are the conjugate exponents, we have the
duality [11]

$\dot{B}^{\sigma}_{p,q}(\Omega)=(B^{-\sigma}_{p^{\prime},q^{\prime}}(\Omega))^{\prime}.$

It is common practice to furnish Besov spaces with equivalent norms based on
$L^{p}$-norms of difference quotients, instead of the interpolation norm.
Given $\rho>0$,

$\Omega_{\rho}:=\\{x\in\Omega:\textrm{dist}(x,\partial\Omega)<\rho\\},\quad\Omega^{\rho}:=\\{x\in\mathbb{R}^{d}:\textrm{dist}(x,\partial\Omega)>\rho\\},$

and a set of admissible directions $D\subset{\mathbb{R}^{d}}$, typically a
ball, we denote

$|v|_{B^{\sigma}_{p,q}(\Omega;D)}:=\Big{(}q\sigma(2-\sigma)\int_{D}\frac{\|v_{h}-2v+v_{-h}\|_{L^{p}(\Omega_{|h|})}^{q}}{|h|^{d+q\sigma}}dh\Big{)}^{1/q}$

for $p,q\in[1,\infty)$ while for $q=\infty$ we let

$|v|_{B^{\sigma}_{p,\infty}(\Omega;D)}:=\sup_{h\in
D}\frac{\|v_{h}-2v+v_{-h}\|_{L^{p}(\Omega_{|h|})}}{|h|^{\sigma}},$

where $v_{h}(x):=v(x+h)$ is the translation with vector
$h\in{\mathbb{R}^{d}}$. It is well-known that, if $D$ is a ball, then the norm
$\|\cdot\|_{L^{p}(\Omega)}+|\cdot|_{B^{\sigma}_{p,q}(\Omega;D)}$ is equivalent
to the Besov norm $\|\cdot\|_{B^{\sigma}_{p,q}(\Omega)}$ defined through
interpolation [3, Theorem 7.47]. Moreover, [11, Proposition 2.2] shows that
balls $D$ can be replaced by suitable convex cones in the definition of Besov
seminorms for $q=\infty$. More precisely, let us assume
$D\subset{\mathbb{R}^{d}}$ is bounded and star-shaped with respect to the
origin. We say that $D$ generates ${\mathbb{R}^{d}}$ if there exists
$\rho_{0}(D)>0$ such that for every $\rho\leq\rho_{0}(D)$ and every $h\in
D_{\rho}(0)$, the ball of radius $\rho$ and center $0$, there exists
$\\{h_{j}\\}_{j=1}^{d}\subset D\cup(-D)$ satisfying

$h=\sum_{j=1}^{d}h_{j},\quad\sum_{j=1}^{d}|h_{j}|\leq c|h|$

with a constant $c>0$ only dependent on $D$. The following equivalence is
proved in [11, Proposition 2.2].

######  Proposition 2.1 (Besov seminorms using cones).

Let $D$ be a convex cone generating ${\mathbb{R}^{d}}$ and let
$B\subset{\mathbb{R}^{d}}$ be a ball. If $\sigma\in(0,2)$ and
$p\in[1,\infty)$, then for every function
$v\colon{\mathbb{R}^{d}}\to\mathbb{R}$ we have
$|v|_{B^{\sigma}_{p,\infty}(\Omega;D)}\simeq|v|_{B^{\sigma}_{p,\infty}(\Omega;B)}$.

We will decompose $\Omega$ into overlapping balls and apply Proposition 2.1 to
subdomains $\omega$ made of intersections of such balls with $\Omega$. The
convex cone $D$ will depend on $\omega$ and will be dictated by the Lipschitz
property of $\Omega$. However, we will omit writing $D$ in
$B^{\sigma}_{p,q}(\omega)$ for simplicity of notation but without compromising
clarity.

We can estimate higher-order Besov norms, possibly of order higher than one,
in terms of difference quotients of Besov norms of order less than one. We
express this instrumental reiteration property as follows and refer to [11,
Proposition 2.1].

######  Proposition 2.2 (reiteration of Besov seminorms).

Let $\omega\subset{\mathbb{R}^{d}}$ be a bounded Lipschitz domain,
$s\in(0,1)$, $p,q\in[1,\infty]$, $\sigma\in(0,1]$ and let $D$ be a set
generating ${\mathbb{R}^{d}}$ and star-shaped with respect to the origin.
Then,

$\begin{split}&|v|_{B^{s+\sigma}_{p,q}(\omega)}\lesssim\left(\int_{D}\frac{|v-v_{h}|^{q}_{W^{s}_{p}(\omega)}}{|h|^{d+q\sigma}}\,dh\right)^{1/q},\quad
q\in[1,\infty),\\\ &|v|_{B^{s+\sigma}_{p,\infty}(\omega)}\lesssim\sup_{h\in
D}\frac{1}{|h|^{\sigma}}|v-v_{h}|_{W^{s}_{p}(\omega)}.\end{split}$

###  2.2. Lipschitz domains

We next briefly state a few well-known but relevant results regarding
Lipschitz domains in ${\mathbb{R}^{d}}$.

######  Definition 2.1 (admissible outward vectors).

For every $x_{0}\in{\mathbb{R}^{d}}$ and $\rho\in(0,1]$, we define the set of
admissible outward vectors

$\mathcal{O}_{\rho}(x_{0})=\\{h\in{\mathbb{R}^{d}}\colon|h|\leq\rho,(B_{2\rho}(x_{0})\setminus\Omega)+th\subset\Omega^{c},\
\forall t\in[0,1]\\}.$

An important fact about bounded Lipschitz domains is that they satisfy a
uniform cone property. This can be stated in the following fashion [32,
§1.2.2].

######  Proposition 2.3 (uniform cone property).

If $\Omega$ is a bounded Lipschitz domain, then there exist $\rho\in(0,1]$,
$\theta\in(0,\pi]$ and a map $\textbf{n}\colon{\mathbb{R}^{d}}\to S^{d-1}$
such that, for every $x\in{\mathbb{R}^{d}}$,

$\mathcal{C}_{\rho}(\textbf{n}(x),\theta):=\\{h\in{\mathbb{R}^{d}}\colon|h|\leq\rho,\
h\cdot\textbf{n}\geq|h|\cos\theta\\}\subset\mathcal{O}_{\rho}(x).$

Besov seminorms can be equivalently written as sums of norms over partitions,
as long as the partitions have some overlap. We refer to [11, Lemma 2.6].

######  Lemma 2.2 (localization).

Let $p,q\in[1,\infty]$ and $\sigma\in(0,2)$. Let $\\{D_{j}\\}_{j=1}^{J}$ be a
finite covering of $\Omega$ by balls of radius $\rho$,
$D_{j}=D_{\rho}(x_{j})$. Then, $v\in B^{\sigma}_{p,q}(\Omega)$ if and only if
$v\big{|}_{\Omega\cap B_{j}}\in B^{\sigma}_{p,q}(\Omega\cap D_{j})$ for all
$j=1,\ldots,J$, and

(2.5)
$\|v\|_{B^{\sigma}_{p,q}(\Omega)}^{p}\simeq\sum_{j=1}^{J}\|v\|_{B^{\sigma}_{p,q}(\Omega\cap
D_{j})}^{p}.$

Moreover, for $\delta\geq\rho$, let $\\{D_{j}\\}_{j=1}^{J}$ be a finite
covering of $\Omega^{\delta}$ and let $v\colon{\mathbb{R}^{d}}\to\mathbb{R}$
be such that $\mbox{supp}(v)\subset\overline{\Omega}$. Then,
$v\in\dot{B}^{\sigma}_{p,q}(\Omega)$ if and only if $v\big{|}_{D_{j}}\in
B^{\sigma}_{p,q}(D_{j})$ for all $j=1,\ldots,J$, and

(2.6)
$\|v\|_{\dot{B}^{\sigma}_{p,q}(\Omega)}^{p}\simeq|v|_{\dot{B}^{\sigma}_{p,q}(\Omega)}^{p}\simeq\sum_{j=1}^{J}|v|_{B^{\sigma}_{p,q}(D_{j})}^{p}.$

The equivalence constants above depend on $s,p,q,\Omega$ and the covering
chosen.

###  2.3. Localized translation operator

Our next goal is to construct a smooth operator that resembles a translation
around a certain given point $x_{0}\in{\mathbb{R}^{d}}$, while coincides with
the identity away from $x_{0}$. Such localized translation operator plays an
instrumental role in our derivation of regularity estimates.

Given $x_{0}$ and $\rho$, we fix a cut-off function $\phi$ such that
$0\leq\phi\leq 1$, $\phi\equiv 1$ on the ball $D_{\rho}(x_{0})$ of radius
$\rho$ centered at $x_{0}$, $\mbox{supp}(\phi)\subset B_{2\rho}(x_{0})$. Given
$h\in{\mathbb{R}^{d}}$, we define

(2.7) $T_{h}v(x):=v\big{(}x+h\phi(x)\big{)}=\big{(}v\circ S_{h}\big{)}(x),$

where the map $S_{h}:=I+h\phi$ is defined from ${\mathbb{R}^{d}}$ to
${\mathbb{R}^{d}}$. We restrict our consideration to $|h|$ small enough such
that the Jacobian of $S_{h}$ satisfies

$\frac{1}{2}I\preccurlyeq\nabla S_{h}=I+h\otimes\nabla\phi\preccurlyeq 2I.$

and thus $S_{h}$ is a one-to-one mapping from $D_{2\rho}(x_{0})$ to
$D_{2\rho}(x_{0})$. It is also one-to-one from ${\mathbb{R}^{d}}$ to
${\mathbb{R}^{d}}$ and coincides with the identity in $D_{2\rho}(x_{0})^{c}$.

In [43, 11], the localized translation operator $\widetilde{T}_{h}v:=\phi
v_{h}+(1-\phi)v$ was employed instead. The translation operator $T_{h}$ in
(2.7) is somewhat more flexible than $\widetilde{T}_{h}$, in the sense that it
gives rise to cleaner regularity estimates in which a priori one gains one
full derivative; compare the right-hand side in (3.5) below with the one in
[11, formula (3.3)]. This leads to a simpler bootstrapping argument than in
[11]. We return to this point in Remark 9 below.

######  Remark 1 (properties of $S_{h}$ and $S_{h}^{-1}$).

Some important properties of the transformations $S_{h}$ and $S_{h}^{-1}$
follow immediately from their definitions. We have the inequalities

(2.8) $\displaystyle|S_{h}(x)-x|\lesssim|h|,\quad\forall
x\in{\mathbb{R}^{d}},$ (2.9) $\displaystyle\left|\det(\nabla
S_{h}(x))-1\right|\lesssim|h|,\quad\forall x\in{\mathbb{R}^{d}},$ (2.10)
$\displaystyle\left|\frac{|x-y|}{|S_{h}(x)-S_{h}(y)|}-1\right|\lesssim|h|,\quad\forall
x,y\in{\mathbb{R}^{d}}.$

Analogous properties also hold for $S^{-1}_{h}$.

######  Remark 2 (boundedness of translations).

Clearly, the operator $T_{h}$ in (2.7) is bounded from
$L^{p}({\mathbb{R}^{d}})$ to $L^{p}({\mathbb{R}^{d}})$ and, more in general,
from $W^{k}_{p}({\mathbb{R}^{d}})$ to $W^{k}_{p}({\mathbb{R}^{d}})$ for every
$k\in\mathbb{N}$. Therefore, it is also bounded from
$B^{\sigma}_{p,q}({\mathbb{R}^{d}})$ to $B^{\sigma}_{p,q}({\mathbb{R}^{d}})$
for any non-integer $\sigma>0$. Moreover, if $h\in\mathcal{O}_{\rho}(x_{0})$
is an admissible outward vector (cf. Definition 2.1) and
$v\in\dot{B}^{\sigma}_{p,q}(\Omega)$, we have

$\displaystyle x\in\Omega^{c}\cap D_{2\rho}(x_{0}):\quad 0\leq\phi(x)\leq
1\Rightarrow S_{h}(x)=x+h\phi(x)\in\Omega^{c}\Rightarrow T_{h}v(x)=0,$
$\displaystyle x\in\Omega^{c}\setminus
D_{2\rho}(x_{0}):\quad\phi(x)=0\Rightarrow T_{h}v(x)=v(x)=0.$

Therefore, $T_{h}$ is also a bounded operator from
$\dot{B}^{\sigma}_{p,q}(\Omega)$ to $\dot{B}^{\sigma}_{p,q}(\Omega)$.

######  Lemma 2.3 (moduli of continuity).

Let $p,q\in[1,\infty]$. For all $\sigma\in(0,1)$ we have

(2.11)
$\|v-T_{h}v\|_{L^{p}(D_{2\rho}(x_{0}))}\lesssim|h|^{\sigma}\|v\|_{B^{\sigma}_{p,q}(D_{2\rho}(x_{0}))}\quad\forall
v\in B^{\sigma}_{p,q}(D_{2\rho}(x_{0})).$

Moreover, for all $r>0$ and $\sigma\in[0,1]$, we have

(2.12)
$\|v-T_{h}v\|_{B^{r}_{p,q}(D_{2\rho}(x_{0}))}\lesssim|h|^{\sigma}\|v\|_{B^{r+\sigma}_{p,q}(D_{2\rho}(x_{0}))}\quad\forall
v\in B^{r+\sigma}_{p,q}(D_{2\rho}(x_{0})).$

###### Proof.

In first place, we remark that because $T_{h}v=v\circ S_{h}$, $T_{h}v-v$
vanishes on $D_{2\rho}(x_{0})^{c}$. The boundedness of $T_{h}$ on
$L^{p}({\mathbb{R}^{d}})$ yields

$\|v-T_{h}v\|_{L^{p}(D_{2\rho}(x_{0}))}\lesssim\|v\|_{L^{p}(D_{2\rho}(x_{0}))},$

while a standard calculation gives

(2.13) $\|v-T_{h}v\|_{L^{p}(D_{2\rho}(x_{0}))}\lesssim|h|\,\|\nabla
v\|_{L^{p}(D_{2\rho}(x_{0}))}.$

Estimate (2.11) follows by interpolation.

We next consider higher-order derivatives. Given a positive integer $k$ and
$\alpha$ with $|\alpha|=k$, a direct calculation gives for any sufficiently
smooth $v$ and $x\in D_{2\rho}(x_{0})$,

$\left|\left(D^{\alpha}v\right)\circ S_{h}(x)-D^{\alpha}\left(v\circ
S_{h}(x)\right)\right|\lesssim|h|\|\phi\|_{C^{k}({\mathbb{R}^{d}})}\sum_{0<|\alpha^{\prime}|\leq
k}\big{|}\big{(}D^{\alpha^{\prime}}v\big{)}\circ S_{h}(x)\big{|}.$

Consequently, applying (2.13) to $\left(D^{\alpha}v\right)\circ S_{h}$, we
deduce

$\displaystyle\|v-T_{h}v\|_{W^{k}_{p}(D_{2\rho}(x_{0}))}=\|v-v\circ
S_{h}\|_{W^{k}_{p}(D_{2\rho}(x_{0}))}$
$\displaystyle\leq\\!\\!\\!\\!\sum_{0\leq|\alpha|\leq
k}\\!\\!\\!\\!\|\\!\left(D^{\alpha}v\right)\\!\circ\\!S_{h}\\!-\\!D^{\alpha}\left(v\\!\circ\\!S_{h}\right)\\!\|_{L^{p}(D_{2\rho}(x_{0}))}\\!+\\!\|\\!\left(D^{\alpha}v\right)\\!\circ\\!S_{h}\\!-\\!D^{\alpha}v\|_{L^{p}(D_{2\rho}(x_{0}))}$
$\displaystyle\lesssim|h|\\!\\!\sum_{0<|\alpha^{\prime}|\leq
k+1}\|D^{\alpha^{\prime}}v\|_{L^{p}(D_{2\rho}(x_{0}))}\leq|h|\|v\|_{W^{k+1}_{p}(D_{2\rho}(x_{0}))}.$

By using the boundedness of $T_{h}$ on $W^{k}_{p}({\mathbb{R}^{d}})$ and
interpolation, we deduce that

$\|v-T_{h}v\|_{W^{k}_{p}(D_{2\rho}(x_{0}))}\lesssim|h|^{\sigma}\|v\|_{B^{k+\sigma}_{p,q}(D_{2\rho}(x_{0}))}\quad\forall
v\in B^{k+\sigma}_{p,q}(D_{2\rho}(x_{0}))$

for $\sigma\in(0,1)$. Finally, we obtain (2.12) by combining this estimate
with (2.3). ∎

######  Remark 3.

The only properties of $S_{h}$ we exploited in the previous lemma are the fact
that it maps $D_{2\rho}(x_{0})$ onto $D_{2\rho}(x_{0})$, and that the
resulting translation operator $T_{h}v=v\circ S_{h}$ is stable in
$W^{k}_{p}(D_{2\rho}(x_{0}))$ and satisfies properties like (2.13). Thus, the
same arguments can be applied to the translation operator $v\circ S_{h}^{-1}$.
In particular, we have

$\|v-v\circ
S_{h}^{-1}\|_{W^{k}_{p}(D_{2\rho}(x_{0}))}\lesssim|h|^{\sigma}|v|_{B^{k+\sigma}_{p,q}(D_{2\rho}(x_{0}))}\quad\forall
v\in B^{k+\sigma}_{p,q}(D_{2\rho}(x_{0}))$

for $\sigma\in(0,1)$, and

$\|v-v\circ
S_{h}^{-1}\|_{B^{r}_{p,q}(D_{2\rho}(x_{0}))}\lesssim|h|^{\sigma}\|v\|_{B^{r+\sigma}_{p,q}(D_{2\rho}(x_{0}))}\quad\forall
v\in B^{r+\sigma}_{p,q}(D_{2\rho}(x_{0}))$

for $r>0,\sigma\in[0,1]$.

###  2.4. Assumptions on the energy

We recall the energy (1.1),

$\mathcal{F}(u)=\iint_{Q_{\Omega}}G\left(x,y,\frac{u(x)-u(y)}{|x-y|^{s}}\right)\frac{1}{|x-y|^{d}}\,dydx-\langle
f,u\rangle.$

Here, we list the conditions we require on this functional, discuss some
consequences of these conditions, and how they apply to the problems we are
interested in.

######  Hypothesis 2.1.

The function
$G\colon\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}\to(0,\infty)$
satisfies the following conditions.

1\. Convexity. The function $\rho\mapsto G(x,y,\rho)$ is uniformly convex for
each $x,y\in{\mathbb{R}^{d}}$.

2\. Space continuity. There exists $\beta\in(0,1]$ such that $G$ is
$\beta$-Hölder continuous with respect to the space variables: for all
$x,y,x^{\prime},y^{\prime}\in{\mathbb{R}^{d}}$ and $\rho\in\mathbb{R},$

(2.14) $|G(x,y,\rho)-G(x^{\prime},y^{\prime},\rho)|\leq
C(|x-x^{\prime}|^{\beta}+|y-y^{\prime}|^{\beta})|\rho|^{p}.$

3\. $p$-growth. There exists some $p\in(1,\infty)$ such that, for all
$x,y\in{\mathbb{R}^{d}}$, $G(x,y,\rho)$ is differentiable with respect to
$\rho$ with

(2.15) $|G(x,y,\rho)|\leq C|\rho|^{p},\qquad|G_{\rho}(x,y,\rho)|\leq
C|\rho|^{p-1}.$

4\. Monotonicity. There exists $\alpha>0$ such that, if $2\leq p<\infty$ in
the $p$-growth condition, then for all $x,y,\in{\mathbb{R}^{d}}$ and
$\rho,\rho^{\prime}\in\mathbb{R},$

$(G_{\rho}(x,y,\rho)-G_{\rho}(x,y,\rho^{\prime}))(\rho-\rho^{\prime})\geq\alpha|\rho-\rho^{\prime}|^{p},$

while if $1<p<2$, then for all $x,y,\in{\mathbb{R}^{d}}$ and
$\rho,\rho^{\prime}\in\mathbb{R},$

$(G_{\rho}(x,y,\rho)-G_{\rho}(x,y,\rho^{\prime}))(\rho-\rho^{\prime})\geq\alpha|\rho-\rho^{\prime}|^{2}\big{|}|\rho|+|\rho^{\prime}|\big{|}^{p-2}.$

5\. Continuity. There exists $c>0$ such that, if $2\leq p<\infty$ in the
$p$-growth condition, then

$|G_{\rho}(x,y,\rho)-G_{\rho}(x,y,\rho^{\prime})|\leq
c|\rho-\rho^{\prime}|\,\big{|}|\rho|+|\rho^{\prime}|\big{|}^{p-2}\quad\forall
x,y,\in{\mathbb{R}^{d}},\ \rho,\rho^{\prime}\in\mathbb{R},$

while if $1<p<2$, then

$|G_{\rho}(x,y,\rho)-G_{\rho}(x,y,\rho^{\prime})|\leq
c|\rho-\rho^{\prime}|^{p-1}\quad\forall x,y,\in{\mathbb{R}^{d}},\
\rho,\rho^{\prime}\in\mathbb{R}.$

######  Remark 4 (symmetry).

While not strictly needed for our purposes, the following assumption is
practical for the analysis and applies to a general class of operators:

6\. Symmetry. The function $G$ is symmetric with respect to the space
variables and with respect to $\rho$,

$G(x,y,\rho)=G(y,x,\rho),\quad G(x,y,\rho)=G(x,y,-\rho),\quad\forall
x,y\in{\mathbb{R}^{d}},\rho\in\mathbb{R}.$

Under this symmetry assumption, the operator $\mathcal{A}$ associated with the
energy minimization problem becomes (1.3); otherwise, it takes the form (1.2).

######  Remark 5 (monotonicity).

The monotonicity hypothesis above implies the following estimates for the
operator $\mathcal{A}$ in (1.2). If $2\leq p<\infty$, there exists $\alpha>0$
such that, for all $u,v\in\widetilde{W}^{s}_{p}(\Omega),$

(2.16)
$\langle\mathcal{A}u-\mathcal{A}v,u-v\rangle\geq\alpha\|u-v\|_{\widetilde{W}^{s}_{p}(\Omega)}^{p};$

hence $\mathcal{A}$ is $p$-coercive in $\widetilde{W}^{s}_{p}(\Omega)$.
Instead, if $1<p<2$, we have for all $u,v\in\widetilde{W}^{s}_{p}(\Omega)$

(2.17)
$\langle\mathcal{A}u-\mathcal{A}v,u-v\rangle\geq\alpha\|u-v\|_{\widetilde{W}^{s}_{p}(\Omega)}^{2}\left(\|u\|_{\widetilde{W}^{s}_{p}(\Omega)}+\|v\|_{\widetilde{W}^{s}_{p}(\Omega)}\right)^{p-2},$

whence $\mathcal{A}$ is $2$-coercive on bounded sets in
$\widetilde{W}^{s}_{p}(\Omega)$:

$\langle\mathcal{A}u-\mathcal{A}v,u-v\rangle\geq
C(R)\|u-v\|_{\widetilde{W}^{s}_{p}(\Omega)}^{2}$

provided
$\|u\|_{\widetilde{W}^{s}_{p}(\Omega)},\|v\|_{\widetilde{W}^{s}_{p}(\Omega)}\leq
R$. It is worth realizing that (2.16) cannot hold for $p<2$ and $G$ smooth and
convex [43, Remark 2.1]. This fact is responsible for the dichotomy between
(1.9) and (1.10) and reveals that our variational approach, which hinges on
(2.17) for $p<2$, cannot improve upon (1.10).

######  Remark 6 (continuity).

The continuity hypothesis implies that the operator $\mathcal{A}$ satisfies
the following bounds: if $p\in(1,2]$ then for all
$u,v\in\widetilde{W}^{s}_{p}(\Omega)$,

(2.18) $\|\mathcal{A}u-\mathcal{A}v\|_{W^{-s}_{p^{\prime}}(\Omega)}\leq
C\|u-v\|_{\widetilde{W}^{s}_{p}(\Omega)}^{p-1},$

while if $p\in[2,\infty)$ then for all $u,v\in\widetilde{W}^{s,p}(\Omega)$,

(2.19) $\|\mathcal{A}u-\mathcal{A}v\|_{W^{-s}_{p^{\prime}}(\Omega)}\leq
C\left(\|u\|_{\widetilde{W}^{s}_{p}(\Omega)}+\|v\|_{\widetilde{W}^{s}_{p}(\Omega)}\right)^{p-2}\|u-v\|_{\widetilde{W}^{s}_{p}(\Omega)}.$

######  Remark 7 (solution operator).

The uniform convexity of $G$ yields existence and uniqueness of weak
solutions: given $f\in W^{-s}_{p^{\prime}}(\Omega)$, the problem

(2.20)
$u_{f}\in\widetilde{W}^{s}_{p}(\Omega)\quad\langle\mathcal{A}u_{f},v\rangle=\langle
f,v\rangle\quad\forall v\in\widetilde{W}^{s}_{p}(\Omega),$

admits a unique solution, where

$\langle\mathcal{A}u,v\rangle=\iint_{{\mathbb{R}^{d}}\times{\mathbb{R}^{d}}}\widetilde{G}\left(x,y,\frac{u(x)-u(y)}{|x-y|^{s}}\right)\frac{(u(x)-u(y))(v(x)-v(y))}{|x-y|^{d+2s}}\,dx\,dy$

and $\widetilde{G}(x,y,\rho):=G_{\rho}(x,y,\rho)/\rho$. Testing (2.20) with
$v=u_{f}$ (or, equivalently, setting $v\equiv 0$ in (2.16)–(2.17)), we
immediately reach the stability estimate

(2.21)
$\|u_{f}\|_{\widetilde{W}^{s}_{p}(\Omega)}\leq\frac{1}{\alpha^{\frac{1}{p-1}}}\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}^{\frac{1}{p-1}}.$

We now assess the continuity properties of the solution operator $f\mapsto
u_{f}$. If we assume $\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}\leq K$, then $u_{f}$
satisfies
$\|u_{f}\|_{\widetilde{W}^{s}_{p}(\Omega)}\leq\left(\frac{K}{\alpha}\right)^{\frac{1}{p-1}}$
in view of (2.21). This shows that, denoting

$\overline{B}_{K}:=\\{f\in
W^{-s}_{p^{\prime}}(\Omega)\colon\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}\leq
K\\},$

and using (2.17), then the solution operator defined on $\overline{B}_{K}$ is
Lipschitz continuous,

(2.22) $\|u_{f}-u_{g}\|_{\widetilde{W}^{s}_{p}(\Omega)}\leq
c(K)\|f-g\|_{W^{-s}_{p^{\prime}}(\Omega)}$

for $p\in(1,2)$. In contrast, if $p\in[2,\infty)$, then the solution map is
Hölder continuous on $W^{-s}_{p^{\prime}}(\Omega)$ because of (2.16),

(2.23)
$\|u_{f}-u_{g}\|_{\widetilde{W}^{s}_{p}(\Omega)}\leq\frac{1}{\alpha^{\frac{1}{p-1}}}\|f-g\|_{W^{-s}_{p^{\prime}}(\Omega)}^{\frac{1}{p-1}}.$

Fractional $(p,s)$-Laplacians. For $p\in(1,\infty)$ we consider
$G(x,y,\rho)=\frac{C_{d,s,p}}{2p}|\rho|^{p}$ in (1.1), which gives rise to
(1.5). It is clear that the parameter $p$ in this definition corresponds to
the parameter $p$ in Hypothesis 2.1, and therefore the convexity and
$p$-growth conditions (2.15) hold. Moreover, because $G$ is independent of the
space variables, it is trivially symmetric with respect to $x,y$, and we can
take $\beta=1$ in (2.14). The monotonicity and continuity assumptions are
satisfied because of the following auxiliary identities [31, §5]: for all
$a,b\in\mathbb{R}$, we have

$\left||a|^{p-2}a-|b|^{p-2}b\right|\leq\left\\{\begin{array}[]{rl}C|a-b|^{p-1}&\mbox{if
}1<p\leq 2,\\\ C|a-b|(|a|+|b|)^{p-2}&\mbox{if }2\leq
p<\infty,\end{array}\right.$

and

$\left(|a|^{p-2}a-|b|^{p-2}b\right)(a-b)\geq\left\\{\begin{array}[]{rl}\alpha|a-b|^{2}(|a|+|b|)^{p-2}&\mbox{if
}1<p\leq 2,\\\ \alpha|a-b|^{p}&\mbox{if }2\leq p<\infty.\end{array}\right.$

Therefore, the Hypothesis 2.1 covers fractional $(p,s)$-Laplace operators
(1.5).

###  2.5. Regularity of functionals

Inspired by [43], we introduce a notion of regularity of functionals that
measures their sensitivity with respect to a family of perturbations.

######  Definition 2.2 ($(T,D,\sigma)$-regularity).

Let $V$ be a Banach space, $K\subset V$ and $\sigma>0$. Given a family of maps
$T_{h}\colon K\to K$, with $h$ varying on a given set
$D\subset{\mathbb{R}^{d}},$ we say that a functional $\mathcal{F}$ is
$(T,D,\sigma)$-regular on $K$ if, for all $v\in K$,

$\omega(v)=\omega(v;\mathcal{F},T,D,\sigma):=\sup_{h\in
D}\frac{\mathcal{F}(T_{h}v)-\mathcal{F}(v)}{|h|^{\sigma}}<\infty.$

######  Remark 8 (subadditivity).

The modulus $\omega$ of $(T,D,\sigma)$-regularity is subadditive with respect
to the $\mathcal{F}$-argument:

(2.24)
$\omega(v;\mathcal{F}_{1}+\mathcal{F}_{2},T,D,\sigma)\leq\omega(v;\mathcal{F}_{1},T,D,\sigma)+\omega(v;\mathcal{F}_{2},T,D,\sigma).$

Thus, in order to prove the $(T,D,\sigma)$-regularity of
$\mathcal{F}_{1}+\mathcal{F}_{2}$, it suffices to show the regularity of each
of the two functionals separately.

A key consequence of the monotonicity assumption in Hypothesis 2.1 is the
following estimate [43, Theorem 1 and Corollary 1].

######  Lemma 2.4 (regularity and minimizers).

Let $x_{0}\in{\mathbb{R}^{d}}$, $\rho>0$ and $h\in\mathcal{O}_{\rho}(x_{0})$.
Consider translation operators
$T_{h}\colon\widetilde{W}^{s}_{p}(\Omega)\to\widetilde{W}^{s}_{p}(\Omega)$ as
in (2.7). If $u$ solves (1.4), the functional $\mathcal{F}$ defined in (1.1)
satisfies Hypotheses 2.1 and it is $(T,D,\sigma)$-regular on
$\widetilde{W}^{s}_{p}(\Omega)$ for some $\sigma>0$, then the following hold:

  * •

If $p\geq 2$, then

(2.25) $\alpha\|u-T_{h}u\|_{\widetilde{W}^{s}_{p}(\Omega)}^{p}\leq
p\omega(u)|h|^{\sigma}.$

  * •

If $1<p<2$, then

(2.26) $\alpha\|u-T_{h}u\|_{\widetilde{W}^{s}_{p}(\Omega)}^{2}\leq
C(p)\left(\|u\|_{\widetilde{W}^{s}_{p}(\Omega)}+\|T_{h}u\|_{\widetilde{W}^{s}_{p}(\Omega)}\right)^{2-p}\omega(u)|h|^{\sigma}.$

Let us explain the crucial role of Lemma 2.4 in the proof of regularity of
solutions, and how localized translations come into play. Let
$D_{\rho}(x_{0})$ be a ball with center $x_{0}\in\Omega$ and radius $\rho$
satisfying Proposition 2.3 (uniform cone property), and let
$D=\mathcal{C}_{\rho}(\textbf{n}(x_{0}),\theta)$. For $p\in[2,\infty)$ we
combine Proposition 2.2 (reiteration of Besov seminorms) with (2.25) to obtain

(2.27)
$\begin{split}|u|^{p}_{B^{s+\sigma/p}_{p,\infty}(D_{\rho}(x_{0}))}&\lesssim\sup_{h\in
D}\frac{|u-u_{h}|^{p}_{W^{s}_{p}(D_{\rho}(x_{0}))}}{|h|^{\sigma}}\\\
&\lesssim\sup_{h\in
D}\frac{\|u-T_{h}u\|^{p}_{\widetilde{W}^{s}_{p}(\Omega)}}{|h|^{\sigma}}\lesssim\omega(u;\mathcal{F},T,D,\sigma),\end{split}$

for $\sigma\in(0,1]$. Instead for $p\in(1,2)$, combining Proposition 2.2 with
(2.26) yields

(2.28)
$\begin{split}|u|^{p}_{B^{s+\sigma/2}_{p,\infty}(D_{\rho}(x_{0}))}&\lesssim\sup_{h\in
D}\frac{\|u-T_{h}u\|^{2}_{\widetilde{W}^{s}_{p}(\Omega)}}{|h|^{\sigma}}\\\
&\lesssim\Big{(}\|u\|_{\widetilde{W}^{s}_{p}(\Omega)}+\|T_{h}u\|_{\widetilde{W}^{s}_{p}(\Omega)}\Big{)}^{2-p}\omega(u;\mathcal{F},T,D,\sigma).\end{split}$

This shows that proving that $\mathcal{F}$ satisfies Definition 2.2
($(T,D,\sigma)$-regularity) for some $\sigma\in(0,1]$ gives rise to local
regularity of solutions.

##  3\. Regularity

In this section we obtain regularity estimates for solutions to (2.20). For
that purpose, we analyze the regularity of the functional $\mathcal{F}$ in
(1.1) in the sense of Definition 2.2, and exploit the crucial property that
the operator $T_{h}$ (2.7) is locally a translation and
$T_{h}v\in\widetilde{W}^{s}_{p}(\Omega)$ for all
$v\in\widetilde{W}^{s}_{p}(\Omega)$.

We split the energy in (1.1) as $\mathcal{F}=\mathcal{F}_{G}-\mathcal{F}_{1}$,
with

(3.1)
$\begin{split}&\mathcal{F}_{G}(u):=\iint_{{\mathbb{R}^{d}}\times{\mathbb{R}^{d}}}G\left(x,y,\frac{u(x)-u(y)}{|x-y|^{s}}\right)\frac{1}{|x-y|^{d}}\,dydx,\\\
&\mathcal{F}_{1}(u):=\langle f,u\rangle\end{split}$

and recall that, by the subadditivity property (2.24), we can treat the two
terms separately.

###  3.1. Regularity of the functionals

As a first step towards deriving regularity estimates for minimizers of (1.1),
we prove the regularity of $\mathcal{F}_{1}$ and $\mathcal{F}_{G}$ with
respect to the family of admissible outward vectors, cf. Definition 2.1. We
begin with an estimate on the linear part of the functional.

######  Proposition 3.1 (regularity of $\mathcal{F}_{1}$).

Let $q\in(1,\infty]$, $p\in(1,\infty)$, $p^{\prime}=p/(p-1)$, $\sigma\in(0,1]$
and $t\in(-1,\sigma)$. If $f\in B^{t}_{p^{\prime},q^{\prime}}(\Omega)$, then
$\mathcal{F}_{1}$ is $(T,\mathcal{O}_{\rho}(x_{0}),\sigma)$-regular in
$\dot{B}^{\sigma-t}_{p,q}(\Omega)$ for all $x_{0}\in\Omega$, $\rho>0$. Namely,
for all $v\in\dot{B}^{\sigma-t}_{p,q}(\Omega)$

(3.2)
$\sup_{h\in\mathcal{O}_{\rho}(x_{0})}\frac{\mathcal{F}_{1}(T_{h}v)-\mathcal{F}_{1}(v)}{|h|^{\sigma}}\lesssim\|f\|_{B^{t}_{p^{\prime},q^{\prime}}(D_{2\rho}(x_{0})\cap\Omega)}\|v\|_{B^{\sigma-t}_{p,q}(D_{2\rho}(x_{0}))}.$

###### Proof.

We split the proof into three steps depending on the range of $\sigma$. Let
$\sigma\in(0,1]$, $r\geq 0$, and $v\in\dot{B}^{r+\sigma}_{p,q}(\Omega)$.

1\. Case $t\in(-1,0]$: We use (2.7) to write $T_{h}v=v\circ S_{h}$ with
$S_{h}=I+h\phi$, and recall that $\mbox{supp}(\phi)\subset D_{2\rho}(x_{0})$,
so that $v-T_{h}v\in\dot{B}^{r}_{p,q}(D_{2\rho}(x_{0}))$. We thus obtain

$\mathcal{F}_{1}(T_{h}v)-\mathcal{F}_{1}(v)=\langle
f,T_{h}v-v\rangle\leq\|f\|_{B^{-r}_{p^{\prime},q^{\prime}}(D_{2\rho}(x_{0})\cap\Omega)}\|T_{h}v-v\|_{\dot{B}^{r}_{p,q}(D_{2\rho}(x_{0}))}.$

Next, we resort to (2.12) to deduce

$\|T_{h}v-v\|_{B^{r}_{p,q}(B_{2\rho}(x_{0}))}\leq|h|^{\sigma}\|v\|_{B^{r+\sigma}_{p,q}(D_{2\rho}(x_{0}))},$

which implies for all $v\in\dot{B}^{r+\sigma}_{p,q}(\Omega)$

(3.3)
$\left|\mathcal{F}_{1}(T_{h}v)-\mathcal{F}_{1}(v)\right|\lesssim|h|^{\sigma}\|f\|_{B^{-r}_{p^{\prime},q^{\prime}}(D_{2\rho}(x_{0})\cap\Omega)}\|v\|_{B^{r+\sigma}_{p,q}(D_{2\rho}(x_{0}))}.$

This establishes (3.2) upon setting $r=-t$.

2\. Case $t=\sigma$: Let $f\in
B^{\sigma}_{p^{\prime},q^{\prime}}(D_{2\rho}(x_{0})\cap\Omega)$ and change
variables to write

$\int_{\Omega}f(x)\,v(S_{h}(x))\,dx=\int_{S_{h}(\Omega)}f(S_{h}^{-1}(x))v(x)|\det\nabla
S_{h}^{-1}(x)|\,dx,$

whence,

$\displaystyle\left|\mathcal{F}_{1}(T_{h}v)-\mathcal{F}_{1}(v)\right|$
$\displaystyle=\left|\int_{B_{2\rho}(x_{0})\cap\Omega}\left[f(S_{h}^{-1}(x))|\det\nabla
S_{h}^{-1}(x)|-f(x)\right]v(x)\,dx\right|$
$\displaystyle\lesssim\|v\|_{L^{p}(D_{2\rho}(x_{0}))}\|(1-|\det\nabla
S_{h}^{-1}|)(f\circ
S_{h}^{-1})\|_{L^{p^{\prime}}(D_{2\rho}(x_{0})\cap\Omega)}$
$\displaystyle+\|v\|_{L^{p}(D_{2\rho}(x_{0}))}\|f\circ
S_{h}^{-1}-f\|_{L^{p^{\prime}}(D_{2\rho}(x_{0})\cap\Omega)},$

because $(1-|\det\nabla S_{h}^{-1}|)(f\circ S_{h}^{-1})$ and $f\circ
S_{h}^{-1}-f$ vanish on $D_{2\rho}(x_{0})^{c}$. We use Remark 1 (properties of
$S_{h}$ and $S_{h}^{-1}$) and Lemma 2.3 (moduli of continuity) to deduce

(3.4)
$\left|\mathcal{F}_{1}(T_{h}v)-\mathcal{F}_{1}(v)\right|\lesssim|h|^{\sigma}\|f\|_{B^{\sigma}_{p^{\prime},q^{\prime}}(D_{2\rho}(x_{0})\cap\Omega)}\|v\|_{L^{p}(D_{2\rho}(x_{0}))}.$

3\. Case $t\in(0,\sigma)$: Since the mapping
$(f,v)\to\mathcal{F}_{1}(T_{h}v)-\mathcal{F}_{1}(v)$ is bilinear, we may
interpolate between (3.3) with $r=0$ and (3.4), with the same $p$ and $q$ in
both expressions, to infer that (3.2) holds as well in this case. ∎

Next, we prove the regularity of the non-linear term $\mathcal{F}_{G}$,
defined in (3.1).

######  Proposition 3.2 (regularity of $\mathcal{F}_{G}$).

Let $s\in(0,1)$ and assume that $G$ satisfies Hypothesis 2.1 for some
$p\in(1,\infty)$ and $\beta\in(0,1]$. Then, the functional
$\mathcal{F}_{G}\colon\widetilde{W}^{s}_{p}(\Omega)\to\mathbb{R}$ defined in
(3.1) is $(T,\mathcal{O}_{\rho}(x_{0}),\beta)$-regular in
$\widetilde{W}^{s}_{p}(\Omega)$ for all $x_{0}\in\Omega$, $\rho>0$. Namely,
for all $v\in\widetilde{W}^{s}_{p}(\Omega)$ it holds that

(3.5)
$\sup_{h\in\mathcal{O}_{\rho}(x_{0})}\frac{\mathcal{F}_{G}(T_{h}v)-\mathcal{F}_{G}(v)}{|h|^{\beta}}\lesssim\iint_{Q_{B_{2\rho}(x_{0})}}\frac{|v(x)-v(y)|^{p}}{|x-y|^{d+sp}}dydx,$

where
$Q_{B_{2\rho}(x_{0})}:=(B_{2\rho}(x_{0})\times{\mathbb{R}^{d}})\cup({\mathbb{R}^{d}}\times
B_{2\rho}(x_{0}))$.

###### Proof.

The change of variables
$(x,y)\mapsto(S_{h}^{-1}(x),S_{h}^{-1}(y))=:(x_{h},y_{h})$ leads to

$\displaystyle\mathcal{F}_{G}(T_{h}v)-\mathcal{F}_{G}(v)$
$\displaystyle=\iint_{{\mathbb{R}^{d}}\times{\mathbb{R}^{d}}}\frac{G\left(x,y,\frac{T_{h}v(x)-T_{h}v(y)}{|x-y|^{s}}\right)-G\left(x,y,\frac{v(x)-v(y)}{|x-y|^{s}}\right)}{|x-y|^{d}}\,dydx$
$\displaystyle=\iint_{{\mathbb{R}^{d}}\times{\mathbb{R}^{d}}}\frac{G\left(x_{h},y_{h},\frac{v(x)-v(y)}{|x_{h}-y_{h}|^{s}}\right)}{|x_{h}-y_{h}|^{d}}\,J(x,y)dydx$
$\displaystyle\quad-\iint_{{\mathbb{R}^{d}}\times{\mathbb{R}^{d}}}G\left(x,y,\frac{v(x)-v(y)}{|x-y|^{s}}\right)\frac{1}{|x-y|^{d}}\,dydx,$

where $J(x,y):=|\det(\nabla S^{-1}_{h}(x))\det(\nabla S^{-1}_{h}(y))|$ for
conciseness. We further split

$\mathcal{F}_{G}(T_{h}v)-\mathcal{F}_{G}(v)=I+II+III,$

with

$\displaystyle I$
$\displaystyle:=\iint_{{\mathbb{R}^{d}}\times{\mathbb{R}^{d}}}\frac{G\left(x_{h},y_{h},\frac{v(x)-v(y)}{|x_{h}-y_{h}|^{s}}\right)-G\left(x_{h},y_{h},\frac{v(x)-v(y)}{|x-y|^{s}}\right)}{|x_{h}-y_{h}|^{d}}\,J(x,y)\,dydx,$
$\displaystyle II$
$\displaystyle:=\iint_{{\mathbb{R}^{d}}\times{\mathbb{R}^{d}}}\frac{G\left(x_{h},y_{h},\frac{v(x)-v(y)}{|x-y|^{s}}\right)-G\left(x,y,\frac{v(x)-v(y)}{|x-y|^{s}}\right)}{|x_{h}-y_{h}|^{d}}\,J(x,y)\,dydx,$
$\displaystyle III$
$\displaystyle:=\iint_{{\mathbb{R}^{d}}\times{\mathbb{R}^{d}}}G\left(x,y,\frac{v(x)-v(y)}{|x-y|^{s}}\right)\left(\frac{J(x,y)}{|x_{h}-y_{h}|^{d}}-\frac{1}{|x-y|^{d}}\right)dydx.$

We observe that, because $S_{h}$ is a one-to-one mapping on ${\mathbb{R}^{d}}$
that coincides with the identity over $B_{2\rho}(x_{0})^{c}$, all the
integrals above need to be computed on the set
$Q_{B_{2\rho}(x_{0})}=(B_{2\rho}(x_{0})\times{\mathbb{R}^{d}})\cup({\mathbb{R}^{d}}\times
B_{2\rho}(x_{0}))$.

Applying the Mean Value Theorem and using (2.9), (2.10) and the growth
condition (2.15), it follows that

(3.6)
$\left|I\right|\lesssim|h|\iint_{Q_{B_{2\rho}(x_{0})}}\frac{|v(x)-v(y)|^{p}}{|x-y|^{d+sp}}dydx.$

Next, we resort to the $\beta$-Hölder continuity of $G$ with respect to the
space variables (2.14), use (2.8) and (2.9), and obtain

(3.7) $\displaystyle\left|II\right|$
$\displaystyle\lesssim\iint_{Q_{B_{2\rho}(x_{0})}}\left(|x_{h}-x|^{\beta}+|y_{h}-y|^{\beta}\right)\frac{|v(x)-v(y)|^{p}}{|x-y|^{d+sp}}dydx$
$\displaystyle\lesssim|h|^{\beta}\iint_{Q_{B_{2\rho}(x_{0})}}\frac{|v(x)-v(y)|^{p}}{|x-y|^{d+sp}}dydx.$

Finally, we use (2.9) and (2.10) to derive

$\left|\frac{J(x,y)}{|S^{-1}_{h}(x)-S^{-1}_{h}(y)|^{d}}-\frac{1}{|x-y|^{d}}\right|\leq\frac{|h|}{|x-y|^{d}}.$

Therefore, combining this bound with condition (2.15), we obtain

(3.8) $\displaystyle|III|$
$\displaystyle\lesssim|h|\iint_{Q_{B_{2\rho}(x_{0})}}\frac{|v(x)-v(y)|^{p}}{|x-y|^{d+sp}}dydx.$

Collecting (3.6), (3.7), and (3.8), we deduce

$|\mathcal{F}_{G}(T_{h}v)-\mathcal{F}_{G}(v)|\lesssim|h|^{\beta}\iint_{Q_{B_{2\rho}(x_{0})}}\frac{|v(x)-v(y)|^{p}}{|x-y|^{d+sp}}dydx$

and conclude the proof. ∎

######  Remark 9 (translation operator).

The estimate (3.5) depends crucially on the structure of the translation
operator (2.7). Indeed, localizing through a composition enabled us to perform
a simple argument based on change of variables. In contrast, if one aims to
exploit the convexity of the functional $\mathcal{F}_{G}$ by using a
translation operator of the form $\widetilde{T}_{h}=\phi v_{h}+(1-\phi)v$ for
a suitable cutoff $\phi$, as in [43, 11], then one obtains a less accurate
estimate where the right-hand side involves higher-order norms of $u$ [11,
Proposition 3.2].

###  3.2. Regularity of solutions

We are now in position to prove the regularity of minimizers of the energy
(1.1) under Hypothesis 2.1. We first prove estimates in the Besov scale, and
afterwards derive estimates in Sobolev norms by using embeddings. Finally, we
show how to accommodate the theory to include finite-horizon operators.

######  Theorem 3.1 (maximal Besov regularity).

Let $\Omega$ be a bounded Lipschitz domain, $s\in(0,1)$,
$G\colon\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}\to(0,\infty)$
satisfy Hypothesis 2.1, with $p\in(1,\infty)$, $\beta\in(0,1]$, and let
$p^{\prime}=\frac{p}{p-1}$.

If $p\geq 2$ and $f\in
B^{-s+\frac{\beta}{p^{\prime}}}_{p^{\prime},1}(\Omega)$, then
$u\in\dot{B}^{s+\frac{\beta}{p}}_{p,\infty}(\Omega)$ and

(3.9)
$\|u\|_{\dot{B}^{s+\frac{\beta}{p}}_{p,\infty}(\Omega)}\lesssim\|f\|_{B^{-s+\frac{\beta}{p^{\prime}}}_{p^{\prime},1}(\Omega)}^{\frac{1}{p-1}}.$

If $p<2$ and $f\in B^{-s+\frac{\beta}{2}}_{p^{\prime},1}(\Omega)$, then
$u\in\dot{B}^{s+\frac{\beta}{2}}_{p,\infty}(\Omega)$ and

(3.10)
$\|u\|_{\dot{B}^{s+\frac{\beta}{2}}_{p,\infty}(\Omega)}\lesssim\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}^{\frac{2-p}{p-1}}\|f\|_{B^{-s+\frac{\beta}{2}}_{p^{\prime},1}(\Omega)}.$

###### Proof.

Let $\gamma=\frac{\beta}{p^{\prime}}$ if $p\geq 2$ and
$\gamma=\frac{\beta}{2}$ if $p<2$. We first observe that problem (1.4) is
well-posed, because $f\in B^{-s+\gamma}_{p^{\prime},1}(\Omega)\subset
W^{-s}_{p^{\prime}}(\Omega)$ and, according to (2.21), we have

(3.11)
$\|u\|_{\widetilde{W}^{s}_{p}(\Omega)}\lesssim\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}^{\frac{1}{p-1}}\lesssim\|f\|_{B^{-s+\gamma}_{p^{\prime},1}(\Omega)}^{\frac{p^{\prime}}{p}}.$

We split the proof into three steps.

1\. Because $\Omega$ is a Lipschitz domain, by Proposition 2.3 (uniform cone
property) there exist $\rho$, $\varphi$ such that
$\mathcal{C}_{\rho}(\textbf{n}(x),\varphi)\subset\mathcal{O}_{\rho}(x)$ for
all $x\in{\mathbb{R}^{d}}$. We consider a finite covering of $\Omega^{\rho}$
by balls $D_{j}=D(x_{j},\rho)$ of radius $\rho$, $j=1,\ldots,M$. By the
localization estimate (2.6), it suffices to bound Besov seminorms over each of
the balls $D_{j}$.

We consider one of the balls $D_{j}$ in the covering and set
$\mathcal{C}_{j}=\mathcal{C}_{\rho}(\textbf{n}(x_{j}),\varphi)$. Importantly,
if $h\in\mathcal{C}_{j}$ we can guarantee that
$T_{h}u\in\widetilde{W}^{s}_{p}(\Omega)$.

Let $\sigma\in(0,\beta]$ and $t\in(-1,\sigma)$; our proof will use suitable
bounds on the $(T,\mathcal{C}_{j},\sigma)$-regularity modulus of
$\mathcal{F}$. We exploit the subadditivity of this functional, and
Propositions 3.1 (regularity of $\mathcal{F}_{1}$) with $q=\infty$ and 3.2
(regularity of $\mathcal{F}_{G})$):

(3.12)
$\displaystyle\omega(u;\mathcal{F},T,\mathcal{C}_{j},\sigma)\leq\omega(u;\mathcal{F}_{G},T,\mathcal{C}_{j},\sigma)+\omega(u;\mathcal{F}_{1},T,\mathcal{C}_{j},\sigma)$
$\displaystyle\lesssim\iint_{Q_{B_{2\rho}(x_{j})}}\frac{|u(x)-u(y)|^{p}}{|x-y|^{d+sp}}dydx+\|f\|_{B^{t}_{p^{\prime},1}(D_{2\rho}(x_{j})\cap\Omega)}\|u\|_{B^{\sigma-t}_{p,\infty}(D_{2\rho}(x_{j}))}.$

In view of (2.27) and (2.28), using (3.11) we deduce

(3.13)
$|u|_{B^{s+\sigma/q}_{p,\infty}(D_{j})}\lesssim\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}^{\frac{q-p}{p-1}}\omega(u;\mathcal{F},T,\mathcal{C}_{j},\sigma)$

provided $q:=\max\\{2,p\\}$. We next distinguish between $p\in[2,\infty)$ and
$p\in(1,2)$, and employ (3.12) and (3.13) with $t=-s+\beta/q^{\prime}$.

2\. Case $p\in[2,\infty)$. Upon choosing $q=p$ in (3.12), we deduce

$|u|^{p}_{B^{s+\sigma/p}_{p,\infty}(D_{j})}\lesssim\iint_{Q_{B_{2\rho}(x_{j})}}\frac{|u(x)-u(y)|^{p}}{|x-y|^{d+sp}}dydx+\|f\|_{B^{t}_{p^{\prime},1}(D_{2\rho}(x_{j})\cap\Omega)}\|u\|_{B^{\sigma-t}_{p,\infty}(D_{2\rho}(x_{j}))}.$

Adding over $j$ for $1\leq J$, and using Lemma 2.2 (localization),
specifically (2.5) for the right-hand side and (2.6) for the left-hand side,
we obtain

$\|u\|^{p}_{\dot{B}^{s+\sigma/p}_{p,\infty}(\Omega)}\lesssim\|u\|_{\widetilde{W}^{s}_{p}(\Omega)}^{p}+\|f\|_{B^{t}_{p^{\prime},1}(\Omega)}\|u\|_{\dot{B}^{\sigma-t}_{p,\infty}(\Omega)}$

for all $\sigma\in(0,\beta],\ t\in(-1,\sigma)$, where the hidden constant
depends on $J$. We now replace the first term on the right-hand side via the
stability bound (3.11). This suggests the choices
$t=-s+\frac{\beta}{p^{\prime}}<\sigma\leq\beta$ and yields the bound

(3.14) $\|u\|^{p}_{\dot{B}^{s+\sigma/p}_{p,\infty}(\Omega)}\leq
C_{1}\|f\|_{B^{-s+\beta/p^{\prime}}_{p^{\prime},1}(\Omega)}^{p^{\prime}}+C_{2}\|f\|_{B^{-s+\beta/p^{\prime}}_{p^{\prime},1}(\Omega)}\|u\|_{\dot{B}^{\sigma+s-\beta/p^{\prime}}_{p,\infty}(\Omega)}$

for suitable constants $C_{1},C_{2}>0$ depending on $d,s,p,\Omega,J$. We
observe that the differentiability index $\sigma+s-\beta/p^{\prime}$ of $u$ on
the right-hand side gives a larger index $s+\sigma/p$ on the left-hand side
provided $\sigma\leq\beta$. We thus view (3.14) as a recursion relation and
set $\sigma_{0}=0$, and for $k\geq 0$ let $\sigma=\sigma_{k+1}\in(0,\beta]$
and

$s+\frac{\sigma_{k}}{p}=\sigma+s-\frac{\beta}{p^{\prime}}\quad\Rightarrow\quad\sigma_{k+1}:=\frac{\sigma_{k}}{p}+\frac{\beta}{p^{\prime}}.$

This implies $\sigma_{k}=\beta\left(1-\frac{1}{p^{k}}\right)$ and shows that
$\sigma_{k}\in(0,\beta)$ and $\sigma_{k}-t=\sigma_{k}+s-\beta/p^{\prime}\geq
s>0$ for all $k\geq 1$, as desired. We claim that

(3.15)
$\|u\|_{\dot{B}^{s+\sigma_{k}/p}_{p,\infty}(\Omega)}\leq\Lambda_{k}\|f\|_{B^{-s+\frac{\beta}{p^{\prime}}}_{p^{\prime},1}(\Omega)}^{\frac{1}{p-1}}$

for some uniformly bounded constants $\Lambda_{k}$ for $k\geq 0$. We argue by
induction. We first note that (3.11) and the continuity of the embedding
$\widetilde{W}^{s}_{p}(\Omega)\subset\dot{B}^{s}_{p,\infty}(\Omega)$ yield

$\|u\|_{\dot{B}^{s}_{p,\infty}(\Omega)}\leq\Lambda_{0}\|f\|_{B^{-s+\frac{\beta}{p^{\prime}}}_{p^{\prime},1}(\Omega)}^{\frac{1}{p-1}},$

for some $\Lambda_{0}:=\Lambda_{0}(d,s,p,\Omega)$; hence (3.15) is true for
$k=0$. We next set $\sigma=\sigma_{k+1}$ in (3.14) to arrive at

$\displaystyle\|u\|^{p}_{\dot{B}^{s+\sigma_{k+1}/p}_{p,\infty}(\Omega)}$
$\displaystyle\leq\left(C_{1}\|f\|_{B^{-s+\beta/p^{\prime}}_{p^{\prime},1}(\Omega)}^{\frac{1}{p-1}}+C_{2}\|u\|_{\dot{B}^{s+\sigma_{k}/p}_{p,\infty}(\Omega)}\right)\|f\|_{B^{-s+\beta/p^{\prime}}_{p^{\prime},1}(\Omega)}$
$\displaystyle\leq\left(C_{1}+C_{2}\Lambda_{k}\right)\|f\|_{B^{-s+\beta/p^{\prime}}_{p^{\prime},1}(\Omega)}^{p^{\prime}}.$

This shows that (3.15) holds for
$\Lambda_{k+1}:=(C_{1}+C_{2}\Lambda_{k})^{1/p}$; it remains to prove that such
a sequence is bounded. Let
$\Lambda:=\max\\{\Lambda_{0},(p^{\prime}C_{1}+C_{2}^{p^{\prime}})^{1/p}\\}$.
We obviously have $\Lambda_{0}\leq\Lambda$ and, if $\Lambda_{k}\leq\Lambda$,
applying Young’s inequality we obtain

$\begin{split}\Lambda_{k+1}^{p}&=C_{1}+C_{2}\Lambda_{k}\leq
C_{1}+C_{2}\Lambda\leq
C_{1}+\frac{C_{2}^{p^{\prime}}}{p^{\prime}}+\frac{\Lambda^{p}}{p}\\\
&\leq\frac{1}{p^{\prime}}\left(p^{\prime}C_{1}+C_{2}^{p^{\prime}}\right)+\frac{\Lambda^{p}}{p}\leq\Lambda^{p}.\end{split}$

Thus, replacing $\Lambda_{k}$ by $\Lambda$ and letting $k\to\infty$, we have
$\sigma_{k}\to\beta$ and deduce the desired estimate (3.9) for
$p\in[2,\infty)$.

3\. Case $p\in(1,2)$. We choose $q=2$ in (3.13) and $t=-s+\beta/2$ in (3.12).
After squaring, summing up over $j$ for $1\leq j\leq J$, and Lemma 2.2
(localization), we end up with

$\|u\|_{\dot{B}^{s+\sigma/2}_{p,\infty}(\Omega)}^{2}\lesssim\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}^{\frac{2-p}{p-1}}\left(|u|_{\widetilde{W}^{s}_{p}(\Omega)}^{p}+\|f\|_{B^{-s+\beta/2}_{p^{\prime},1}(\Omega)}\|u\|_{\dot{B}^{\sigma+s-\beta/2}_{p,\infty}(\Omega)}\right)$

for all $\sigma\in(0,\beta]$ such that $\sigma+s-\beta/2>0$; the hidden
constant depends on $J$. We next resort to (3.11) to obtain the following
counterpart of (3.14):

(3.16)
$\|u\|_{\dot{B}^{s+\sigma/2}_{p,\infty}(\Omega)}^{2}\leq\left(C_{1}\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}^{p^{\prime}}+C_{2}\|f\|_{B^{-s+\beta/2}_{p^{\prime},1}(\Omega)}\|u\|_{\dot{B}^{\sigma+s-\beta/2}_{p,\infty}(\Omega)}\right)\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}^{\frac{2-p}{p-1}}$

for all $\sigma\in(0,\beta]$ and with constants $C_{1}$ and $C_{2}$ that do
not depend on $u$ or $f$.

We now proceed as in Step 2 to exploit the improvement of the
differentiability index from the right-hand side to the left one in (3.16). To
this end, we set $\sigma_{0}=0,\sigma_{k+1}=\sigma\in(0,\beta]$ and rewrite

$\sigma+s-\frac{\beta}{2}=s+\frac{\sigma_{k}}{2}\quad\Rightarrow\quad\sigma_{k+1}=\frac{\beta+\sigma_{k}}{2}\quad\forall\,k\geq
0.$

This yields $\sigma_{k}=\beta(1-2^{-k})$, which satisfies the restrictions
$\sigma_{k}\in(0,\beta)$, $\sigma_{k}-t=\sigma_{k}+s-\frac{\beta}{2}\geq s>0$
for all $k\geq 1$, and $\sigma_{k}\to\beta$. We prove by induction that

(3.17)
$\|u\|_{\dot{B}^{s+\frac{\sigma_{k}}{2}}_{p,\infty}(\Omega)}\leq\Lambda_{k}\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}^{\frac{2-p}{p-1}+2^{-k}}\|f\|_{B^{-s+\frac{\beta}{2}}_{p^{\prime},1}(\Omega)}^{1-2^{-k}}\quad\forall
k\geq 0,$

with a uniformly bounded constant $\Lambda_{k}$. We first note that, according
to (3.11) and the continuity of the embedding
$\widetilde{W}^{s}_{p}(\Omega)\subset\dot{B}^{s}_{p,\infty}(\Omega)$, we have

$\|u\|_{\dot{B}^{s}_{p,\infty}(\Omega)}\leq\Lambda_{0}\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}^{\frac{1}{p-1}}$

for some $\Lambda_{0}:=\Lambda_{0}(d,s,p,\Omega)$, whence (3.17) is true for
$k=0$. Suppose that (3.17) holds for some $k$, and set $\sigma=\sigma_{k+1}$
in (3.16) to obtain

$\begin{split}\|u\|_{\dot{B}^{s+\sigma_{k+1}/2}_{p,\infty}(\Omega)}^{2}&\leq\left(C_{1}\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}^{p^{\prime}}+C_{2}\|f\|_{B^{-s+\beta/2}_{p^{\prime},1}(\Omega)}\|u\|_{\dot{B}^{s+\sigma_{k}/2}_{p,\infty}(\Omega)}\right)\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}^{\frac{2-p}{p-1}}\\\
&\leq\left(C_{1}\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}^{2-2^{-k}}+C_{2}\Lambda_{k}\|f\|_{B^{-s+\beta/2}_{p^{\prime},1}(\Omega)}^{2-2^{-k}}\right)\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}^{\frac{2(2-p)}{p-1}+2^{-k}}\\\
&\leq\left(C_{1}C_{3}^{2-2^{-k}}+C_{2}\Lambda_{k}\right)\|f\|_{B^{-s+\beta/2}_{p^{\prime},1}(\Omega)}^{2-2^{-k}}\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}^{\frac{2(2-p)}{p-1}+2^{-k}},\end{split}$

where $C_{3}:=C_{3}(\Omega,d,s,p,\beta)$ is the constant of the continuous
embedding $W^{-s}_{p^{\prime}}(\Omega)\subset
B^{-s+\beta/2}_{p^{\prime},1}(\Omega)$. Since $1\leq 2-2^{-k}\leq 2$, we have
$C_{3}^{2-2^{-k}}\leq\max\\{C_{3},C_{3}^{2}\\}=:C_{4}$ which gives rise to
(3.17) with constant

$\Lambda_{k+1}:=\left(C_{1}C_{4}+C_{2}\Lambda_{k}\right)^{1/2}.$

It only remains to show that $\Lambda_{k}\leq\Lambda$ for some $\Lambda>0$ and
all $k\geq 0$. Let
$\Lambda:=\max\\{\Lambda_{0},(2C_{1}C_{4}+C_{2}^{2})^{1/2}\\}$, then the same
argument as in Step 2 is also valid in this setting. Finally, estimate (3.10)
follows by letting $k\to\infty$. ∎

Theorem 3.1 gives the maximal Besov regularity one can expect via our
variational approach for solutions to problem (1.4). We recall that (3.9) is
optimal while (3.10) is suboptimal for the 1$d$ profile $(1-x^{2})_{+}^{s}$ in
$\Omega=(-1,1)$. This is due to the Hypothesis 2.1 (monotonicity) for
$p\in(1,2)$, which cannot be improved for a function $\rho\mapsto
G(\cdot,\cdot,\rho)$ convex and differentiable [43, Remark 2.1].

One can immediately deduce Sobolev regularity by combining Theorem 3.1 with
Lemma 2.1 (embedding). Even though such a lemma is stated for the spaces
$B^{\sigma+\varepsilon}_{p,\infty}(\Omega)$ and $W^{\sigma}_{p}(\Omega)$, the
proof for the zero-extension spaces
$\dot{B}^{\sigma+\varepsilon}_{p,\infty}(\Omega)$ and
$\widetilde{W}^{\sigma}_{p}(\Omega)$ follows by the same arguments.

######  Corollary 3.1 (maximal Sobolev regularity).

Let the assumptions of Theorem 3.1 be satisfies. For $p\in[2,\infty)$ and
$f\in B^{-s+\beta/p^{\prime}}_{p^{\prime},1}(\Omega)$ then

(3.18)
$\|u\|_{\widetilde{W}^{s+\frac{\beta}{p}-\varepsilon}_{p}(\Omega)}\lesssim\varepsilon^{-\frac{1}{p}}\|f\|_{B^{-s+\frac{\beta}{p^{\prime}}}_{p^{\prime},1}(\Omega)}^{\frac{1}{p-1}}$

is valid provided $\varepsilon\in(0,s+\beta/p)$. If $p\in(1,2)$ and $f\in
B^{-s+\beta/2}_{p^{\prime},1}(\Omega)$ then

(3.19)
$\|u\|_{\widetilde{W}^{s+\frac{\beta}{2}-\varepsilon}_{p}(\Omega)}\lesssim\varepsilon^{-\frac{1}{p}}\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}^{\frac{2-p}{p-1}}\|f\|_{B^{-s+\frac{\beta}{2}}_{p^{\prime},1}(\Omega)}$

holds provided $\varepsilon\in(0,s+\beta/2)$.

######  Remark 10 (interior Sobolev regularity).

Estimates (3.18) and (3.19) are valid up to the boundary of the domain. The
obstruction to higher regularity is due to boundary behavior. In the
superquadratic case $p\geq 2$, and for the $(p,s)$-fractional Laplacian (1.5)
–that gives rise to $\beta=1$ in Hypothesis 2.1– reference [15] derives the
following higher-order interior regularity provided $f\in
W^{s}_{p^{\prime}}(\Omega)$:

$\begin{split}\mbox{if }s\leq\frac{p-1}{p+1},&\quad\mbox{then
}u\in\bigcap_{\varepsilon>0}W^{s\frac{p+1}{p-1}-\varepsilon}_{p,loc}(\Omega),\\\
\mbox{if }s>\frac{p-1}{p+1},&\quad\mbox{then }u\in W^{1}_{p,loc}(\Omega)\mbox{
and }\nabla
u\in\bigcap_{\varepsilon>0}W^{s\frac{p+1}{p}-\frac{p-1}{p}-\varepsilon}_{p,loc}(\Omega).\end{split}$

We point out that if $s<\frac{1}{2}-\frac{1}{2p}$, estimate (3.18) is actually
stronger than the first statement above.

The technique in Theorem 3.1 also allows one to derive regularity estimates in
Besov spaces with a lower-order differentiability index whenever $f$ is less
regular than in such a theorem. We shall not consider this procedure, but
rather prove a lower regularity pickup by using interpolation theory. For that
purpose, we need the following nonlinear interpolation estimate (cf. [44,
Théorème I.1]).

######  Proposition 3.3 (nonlinear interpolation).

Let $A_{0}\subset A_{1}$, $B_{0}\subset B_{1}$ be Banach spaces, $U\subset
A_{1}$ a nonempty open set and $T\colon U\to B_{1}$ be a a function that maps
$A_{0}\cap U$ into $B_{0}$. Let us assume that there exist constants $c_{0}$,
$c_{1}$ such that

$\begin{split}&\|Tf\|_{B_{0}}\leq c_{0}\|f\|_{A_{0}}^{\alpha_{0}},\quad\forall
f\in A_{0}\cap U,\\\ &\|Tf-Tg\|_{B_{1}}\leq
c_{1}\|f-g\|_{A_{1}}^{\alpha_{1}},\quad\forall f,g\in U\end{split}$

for some $\alpha_{0}>0$, $\alpha_{1}\in(0,1]$. Then, if $\theta\in(0,1)$ and
$q\in[1,\infty]$, $T$ maps $(A_{0},A_{1})_{\theta,q}$ into
$(B_{0},B_{1})_{\eta,r}$, where
$\frac{1-\eta}{\eta}=\frac{\alpha_{1}}{\alpha_{0}}\frac{1-\theta}{\theta}$ and
$r=\max\\{1,\frac{q}{(1-\eta)\alpha_{0}+\eta\alpha_{1}}\\}$.

######  Corollary 3.2 (regularity pickup for rough data).

Let $\Omega$ be a bounded Lipschitz domain,
$G\colon\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}\to(0,\infty)$
satisfy Hypothesis 2.1, with $p\in(1,\infty)$, $\beta\in(0,1]$, and let
$p^{\prime}=\frac{p}{p-1}$, and $\theta\in(0,1)$. Then, the solution operator
$f\mapsto u$ is bounded between the following spaces

$\displaystyle\mbox{ if }p\geq 2\mbox{ and }f\in
W^{-s+\theta\frac{\beta}{p^{\prime}}}_{p^{\prime}}(\Omega)$ $\displaystyle\
\Rightarrow\ u\in\widetilde{W}^{s+\theta\frac{\beta}{p}}_{p}(\Omega);$
$\displaystyle\mbox{ if }1<p<2\mbox{ and }f\in
W^{-s+\theta\frac{\beta}{2}}_{p^{\prime}}(\Omega)$ $\displaystyle\
\Rightarrow\
u\in\widetilde{W}^{s+\theta\frac{\beta}{2}}_{p^{\prime}}(\Omega).$

###### Proof.

The proof follows by a direct application of Proposition 3.3. For $p\geq 2$,
we combine (3.9) and (2.23), while for $1<p<2$ we combine (3.10) and (2.22). ∎

###  3.3. Operators with finite horizon

Thus far, we have obtained regularity estimates under Hypothesis 2.1. In
particular, that hypothesis dictates the smoothness, growth and behavior at
$x=y$ of $G(x,y,\rho)$. However, such global behavior constraints can be
significantly relaxed to incorporate, for instance, finite-horizon operators.
For simplicity, we now assume $G$ takes the form

(3.20)
$G(x,y,\rho)=\frac{1}{2p}\,\psi\left(\frac{|x-y|}{\delta}\right)\,|\rho|^{p},$

where $1<p<\infty$ and $\psi:[0,\infty)\to[0,\infty)$ is a given function. In
case $\psi$ is supported in the unit interval $[0,1]$, the parameter
$\delta>0$ above is the horizon of the resulting operator

(3.21)
$\mathcal{A}_{\delta}u(x):=\int_{{\mathbb{R}^{d}}}\psi\left(\frac{|x-y|}{\delta}\right)\frac{|u(x)-u(y)|^{p-2}(u(x)-u(y))}{|x-y|^{d+sp}}\,dy.$

Let us assume that $\psi\in L^{\infty}({\mathbb{R}^{d}})$. Then, the choice
(3.20) trivially fulfills conditions convexity, $p$-growth, and continuity;
moreover, we note it satisfies the symmetry condition. The only two missing
assumptions from Hypothesis 2.1 in this setting are space continuity and
monotonicity. At this point, we can define the energy norm induced by $G$,

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}:=\langle\mathcal{A}_{\delta}v,v\rangle^{\frac{1}{p}}=\left(\frac{1}{2}\iint_{{\mathbb{R}^{d}}\times{\mathbb{R}^{d}}}\psi\left(\frac{|x-y|}{\delta}\right)\frac{|v(x)-v(y)|^{p}}{|x-y|^{d+sp}}\,dydx\right)^{\frac{1}{p}}$

and realize it satisfies

(3.22)
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\leq\left(\frac{\|\psi\|_{L^{\infty}({\mathbb{R}^{d}})}}{C_{d,s,p}}\right)^{\frac{1}{p}}\,\|v\|_{\widetilde{W}^{s}_{p}(\Omega)}\quad\forall
v\in\widetilde{W}^{s}_{p}(\Omega),$

where $C_{d,s,p}$ is the constant given in (1.6).

Some form of non-degeneracy is needed in order to have a reverse inequality to
(3.22). We assume there exists some $r>0$ such that $\psi\geq\psi_{0}>0$ on
the interval $[0,r]$. This implies the following property:

4’. Local monotonicity. There exists $\alpha>0$ such that for all
$x,y\in{\mathbb{R}^{d}}$ with $|x-y|\leq r\delta$ and all
$\rho,\rho^{\prime}\in\mathbb{R}$,

$\displaystyle(G_{\rho}(x,y,\rho)-G_{\rho}(x,y,\rho^{\prime}))(\rho-\rho^{\prime})\geq\alpha|\rho-\rho^{\prime}|^{p}$
$\displaystyle\quad\mbox{ if }p\geq 2,$
$\displaystyle(G_{\rho}(x,y,\rho)-G_{\rho}(x,y,\rho^{\prime}))(\rho-\rho^{\prime})\geq\alpha|\rho-\rho^{\prime}|^{2}||\rho|+|\rho^{\prime}||^{p-2}$
$\displaystyle\quad\mbox{ if }1<p<2.$

######  Lemma 3.1 (non-degeneracy).

Let $\psi\geq\psi_{0}>0$ on the interval $[0,r]$. Then, there exists a
constant $C=C(d,p,s,\Omega,r,\psi_{0},\delta)$ such that

(3.23) $\|v\|_{\widetilde{W}^{s}_{p}(\Omega)}\leq
C\,{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\quad\forall
v\in\widetilde{W}^{s}_{p}(\Omega).$

###### Proof.

We invoke the localization estimate of [27, Lemma 7], which reads

$\displaystyle\iint_{B_{R}(0)\times B_{R}(0)}$
$\displaystyle\frac{|v(x)-v(y)|^{p}}{|x-y|^{d+sp}}\,dydx$
$\displaystyle\leq\left(\frac{3R}{r\delta}\right)^{p(1-s)}\iint_{B_{R}(0)\times
B_{R}(0)}\frac{|v(x)-v(y)|^{p}}{|x-y|^{d+sp}}\,\chi_{\\{|x-y|\leq
r\delta\\}}\,dydx,$

and is valid for all $R>r\delta>0$ and $s\in(0,1)$. This implies

$\iint_{B_{R}(0)\times
B_{R}(0)}\frac{|v(x)-v(y)|^{p}}{|x-y|^{d+sp}}\,dydx\leq\frac{2(3R)^{p(1-s)}}{(r\delta)^{p(1-s)}\psi_{0}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{p}.$

regardless of whether or not the support of $\psi$ is compact. Given
$v\in\widetilde{W}^{s}_{p}(\Omega)$, we use this bound with $R>0$ sufficiently
large so that $\Omega\subset B_{R}(0)$ and $\mbox{dist}(\Omega,\partial
B_{R}(0))\geq\frac{R}{2}$, exploit the fact that $|x-y|\geq\frac{R}{2}$ for
all $x\in\Omega$ and $y\in B_{R}(0)^{c}$ and integrate in polar coordinates to
get

$\displaystyle\|v\|_{\widetilde{W}^{s}_{p}(\Omega)}^{p}$
$\displaystyle=\frac{C_{d,s,p}}{2}\iint_{B_{R}(0)\times
B_{R}(0)}\frac{|v(x)-v(y)|^{p}}{|x-y|^{d+sp}}\,dydx$
$\displaystyle+C_{d,s,p}\iint_{B_{R}(0)\times
B_{R}(0)^{c}}\frac{|v(x)|^{p}}{|x-y|^{d+sp}}\,dydx$
$\displaystyle\leq\frac{(3R)^{p(1-s)}C_{d,s,p}}{(r\delta)^{p(1-s)}\psi_{0}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{p}+\frac{2^{sp}C_{d,s,p}\omega_{d-1}}{sp\,R^{sp}}\|v\|_{L^{p}(\Omega)}^{p},$

where $\omega_{d-1}=|\mathbb{S}^{d-1}|$ denotes the $(d-1)$-dimensional
measure of the unit sphere $\mathbb{S}^{d-1}=\partial B_{1}(0)$ in
${\mathbb{R}^{d}}$. We next use the well-known Poincaré inequality

$\|v\|_{L^{p}(\Omega)}\leq
C(\Omega,p)\|v\|_{\widetilde{W}^{s}_{p}(\Omega)}\qquad\forall
v\in\widetilde{W}^{s}_{p}(\Omega),$

and take $R>0$ sufficiently large such that

$\frac{2^{sp}C_{d,s,p}\omega_{d-1}}{sp\,R^{sp}}\|v\|_{L^{p}(\Omega)}^{p}\leq\frac{1}{2}\|v\|_{\widetilde{W}^{s}_{p}(\Omega)}^{p}\quad\forall
v\in\widetilde{W}^{s}_{p}(\Omega),$

to obtain a constant $C=C(d,p,s,\Omega,r,\psi_{0},\delta)$ such that (3.23)
holds. ∎

By combining (3.22) and (3.23), we deduce that the energy norm
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
is equivalent to the $\widetilde{W}^{s}_{p}(\Omega)$-norm. Consequently, the
Dirichlet problem for the operator $\mathcal{A}_{\delta}$ defined in (3.21) is
well-posed in $\widetilde{W}^{s}_{p}(\Omega)$ uniformly in $s$: if $f\in
W^{-s}_{p^{\prime}}(\Omega)$, there exists a unique
$u\in\widetilde{W}^{s}_{p}(\Omega)$ satisfying

$\langle\mathcal{A}_{\delta}u,v\rangle=\langle f,v\rangle\quad\forall
v\in\widetilde{W}^{s}_{p}(\Omega),$

and we have the stability bound
$\|u\|_{\widetilde{W}^{s}_{p}(\Omega)}\lesssim\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}$.

Another consequence of the equivalence between the energy norm
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
and the $\widetilde{W}^{s}_{p}(\Omega)$-norm is that the variational approach
of Sections 3.1 and (3.2) hinges on the regularity of $\mathcal{F}_{1}$ and
$\mathcal{F}_{G}$ and still applies in this context regardless of the support
of $\psi$. We state this next.

######  Corollary 3.3 (operator with Hölder continuous $\psi$).

Let $\psi$ in (3.20) be globally $\beta$-Hölder continuous with
$\beta\in(0,1]$ and satisfy $\psi\geq\psi_{0}>0$ on the interval $[0,r]$ for
$r>0$. Then the maximal regularity estimates (3.9) for $2\leq p<\infty$ and
(3.10) for $1<p<2$ are valid regardless of the support of $\psi$.

######  Remark 11 (tempered fractional Laplacians).

Besides being applicable to finite-horizon operators with $\beta$-Hölder
continuous kernel, the previous result is valid for a family of tempered
fractional $p$-Laplacians. Concretely, we let $\lambda=\frac{1}{\delta}>0$ and
$\psi(\rho)=e^{-\lambda\rho}$ to obtain

$\mathcal{A}_{\lambda}u(x):=\int_{\mathbb{R}^{d}}\frac{|u(x)-u(y)|^{p-2}(u(x)-u(y))}{e^{\lambda|x-y|}\,|x-y|^{d+sp}}\,dy.$

In the linear setting ($p=2$), this operator arises from the study of tempered
Lévy flights and has been investigated for example in [24]. For the
homogeneous Dirichlet problem associated to the operator
$\mathcal{A}_{\lambda}$ above, the regularity estimates in Theorem 3.1 are
valid with $\beta=1$.

Finite-horizon operators in practice usually involve a discontinuous function
$\psi$ such as the characteristic function of $[0,r]$. This does not fit
within the preceding theory but maximal Besov regularity is still valid, at
least for linear operators, provided $\psi$ is Hölder in a neighborhood of the
origin. We explore this next, but before we point out that even local
regularity seems excessive, an interesting question to investigate.

We need to make the following local regularity assumption to compensate for
the lack of space continuity hypothesis of $G$:

2’. Local regularity. There exists $r>0$ such that $\psi$ is of class
$C^{\beta}$ on the interval $[0,r]$ for some $\beta\in(0,1]$, and
$\psi\geq\psi_{0}>0$ on $[0,r]$.

######  Theorem 3.2 (linear finite-horizon operator with discontinuous
$\psi$).

Let $\psi$ satisfy the previous local regularity assumption with some
$\beta\in(0,1].$ If $G(x,y,\rho)$ is of the form (3.20) with $p=2$, then the
following maximal regularity holds

(3.24)
$\|u\|_{\dot{B}^{s+\frac{\beta}{2}}_{2,\infty}(\Omega)}\lesssim\|f\|_{B^{-s+\frac{\beta}{2}}_{2,1}(\Omega)}.$

###### Proof.

We resort to a perturbation argument. We proceed in several steps.

1\. Perturbation: Let $\widetilde{\psi}\in C^{\beta}[0,\infty)$ coincide with
$\psi$ on $[0,r]$ and rewrite the equation $\mathcal{A}_{\delta}u=f$ with
$u=0$ in $\Omega^{c}$ as

$\widetilde{\mathcal{A}}_{\delta}u=f+\big{(}\widetilde{\mathcal{A}}_{\delta}u-\mathcal{A}_{\delta}u\big{)}=\widetilde{f}.$

The operator
$\mathcal{B}_{\delta}:=\widetilde{\mathcal{A}}_{\delta}-\mathcal{A}_{\delta}$
is a convolution operator that reads as follows in terms of the function
$\varphi:=\widetilde{\psi}-\psi$, which vanishes on $[0,r]$:

$\mathcal{B}_{\delta}u(x)=\int_{\mathbb{R}^{d}}\varphi\Big{(}\frac{|z|}{\delta}\Big{)}\frac{u(x)-u(x-z)}{|z|^{d+2s}}dz=Ku(x)+k*u(x).$

Above, $k(z)=\varphi\big{(}\frac{|z|}{\delta}\big{)}|z|^{-d-2s}$ and
$K=\int_{\mathbb{R}^{d}}k(z)dz<\infty$ because $\psi\in
L^{\infty}(\mathbb{R}^{d})$.

2\. Properties of $\mathcal{B}_{\delta}$:
$\mathcal{B}_{\delta}:\dot{B}^{t}_{2,\infty}(\Omega)\to{B}^{t}_{2,\infty}({\mathbb{R}^{d}})$
is a linear bounded operator

$\|\mathcal{B}_{\delta}u\|_{\dot{B}^{t}_{2,\infty}(\Omega)}\lesssim\|u\|_{\dot{B}^{t}_{2,\infty}(\Omega)}\quad\forall\,0<t<1.$

It suffices to examine $k*u$ which, using its linear structure and Young’s
inequality with $K=\|k\|_{L^{1}(\mathbb{R}^{d})}$, yields

$\|k*u(\cdot+h)-k*u\|_{L^{2}(\mathbb{R}^{d})}\leq
K\|u(\cdot+h)-u\|_{L^{2}(\mathbb{R}^{d})}\quad\forall\,h\in\mathbb{R}^{d}.$

Consequently, we deduce the asserted estimate from

$\frac{\|\mathcal{B}_{\delta}u(\cdot+h)-\mathcal{B}u\|_{L^{2}(\mathbb{R}^{d})}}{|h|^{t}}\leq
K\frac{\|u(\cdot+h)-u\|_{L^{2}(\mathbb{R}^{d})}}{|h|^{t}}\quad\forall\,0<t<1.$

3\. Regularity of functionals: Combining the local estimate (3.2) of
Proposition 3.1 (regularity of $\mathcal{F}_{1}$) for $\widetilde{f}$, the
argument of Theorem 3.1 (maximal regularity) leading to (3.12) and (3.14)
yields the following estimate for any $\sigma\in(0,\beta]$ and
$t\in(-s,\sigma)$

$\|u\|_{\dot{B}^{s+\sigma/2}_{2,\infty}(\Omega)}^{2}\lesssim\|f\|_{B^{t}_{2,1}(\Omega)}^{2}+\|f\|_{B^{t}_{2,1}(\Omega)}\|u\|_{B^{\sigma-t}_{2,\infty}(\Omega)}+\|\mathcal{B}_{\delta}u\|_{B^{t}_{2,1}(\Omega)}\|u\|_{B^{\sigma-t}_{2,\infty}(\Omega)},$

where we have used the definition of $\widetilde{f}$ and Step 2 to obtain the
last term. To be able to iterate this estimate we observe that
$\|\mathcal{B}_{\delta}u\|_{B^{t}_{2,1}(\Omega)}\lesssim\|\mathcal{B}_{\delta}u\|_{B^{t+2\epsilon}_{2,\infty}(\Omega)}$
for any $\epsilon>0$ and to simplify the subsequent derivation we take
$\epsilon=\frac{s}{2}$ so that $t+2\epsilon=t+s>0$. We thus have
$\|\mathcal{B}_{\delta}u\|_{B^{t}_{2,1}(\Omega)}\lesssim\|u\|_{B^{t+2\epsilon}_{2,\infty}(\Omega)}$
and

(3.25)
$\|u\|_{\dot{B}^{s+\sigma/2}_{2,\infty}(\Omega)}^{2}\lesssim\|f\|_{B^{t}_{2,1}(\Omega)}^{2}+\|f\|_{B^{t}_{2,1}(\Omega)}\|u\|_{B^{\sigma-t}_{2,\infty}(\Omega)}+\|u\|_{B^{t+s}_{2,\infty}(\Omega)}\|u\|_{B^{\sigma-t}_{2,\infty}(\Omega)}.$

4\. Preliminary regularity: We now claim that
$u\in\dot{B}^{\frac{\beta}{2}}_{2,\infty}(\Omega)$ with

(3.26)
$\|u\|_{\dot{B}^{\frac{\beta}{2}}_{2,\infty}(\Omega)}\lesssim\|f\|_{B^{-s+\frac{\beta}{2}}_{2,1}(\Omega)}.$

If $s\geq\frac{\beta}{2}$ this is straightforward, because of the stability
bound $\|u\|_{\widetilde{H}^{s}(\Omega)}\lesssim\|f\|_{H^{-s}(\Omega)}$ and
the continuity of the embedding
$\widetilde{H}^{s}(\Omega)\subset\dot{B}^{\beta/2}_{2,\infty}(\Omega)$.

In the case $s<\frac{\beta}{2}$, we iterate (3.25). We set $\sigma_{0}:=0$ and

(3.27) $\sigma_{k+1}:=\sigma_{k}+\frac{s}{2}.$

In the first few iterations we could have
$0<\sigma_{k}\leq-s+\frac{\beta}{2}$, in which case we set
$\sigma:=2\sigma_{k+1}=s+2\sigma_{k}\in[s,\beta-s]\subset(0,\beta]$ and
$t:=\sigma_{k}=\frac{\sigma-s}{2}\in[0,\frac{\sigma}{2})\subset(-s,\sigma)$.
This choice of parameters yields $\sigma-t=t+s=s+\sigma_{k}$, whence (3.25)
reads

$\|u\|_{\dot{B}^{s+\sigma_{k+1}}_{2,\infty}(\Omega)}^{2}\lesssim\|f\|_{B^{t}_{2,1}(\Omega)}^{2}+\|f\|_{B^{t}_{2,1}(\Omega)}\|u\|_{B^{s+\sigma_{k}}_{2,\infty}(\Omega)}+\|u\|_{B^{s+\sigma_{k}}_{2,\infty}(\Omega)}^{2}.$

Additionally, continuity of the embedding
$B^{-s+\frac{\beta}{2}}_{2,1}(\Omega)\subset B^{t}_{2,1}(\Omega)$, gives

$\|u\|_{\dot{B}^{s+\sigma_{k+1}}_{2,\infty}(\Omega)}^{2}\lesssim\|f\|_{B^{-s+\frac{\beta}{2}}_{2,1}(\Omega)}^{2}+\|u\|_{B^{s+\sigma_{k}}_{2,\infty}(\Omega)}^{2},$

and the bound
$\|u\|_{\dot{B}^{s+\sigma_{k}}_{2,\infty}(\Omega)}\lesssim\|f\|_{B^{-s+\frac{\beta}{2}}_{2,1}(\Omega)}$
valid for $k=0$ implies

$\|u\|_{\dot{B}^{s+\sigma_{k+1}}_{2,\infty}(\Omega)}\lesssim\|f\|_{B^{-s+\frac{\beta}{2}}_{2,1}(\Omega)},$

for as long as $\sigma_{k}\leq-s+\frac{\beta}{2}$. Moreover, because
$\sigma_{k+1}=\sigma_{k}+\frac{s}{2},$ we have a regularity improvement by the
fixed amount $\frac{s}{2}$ in each iteration. Therefore, after a finite (but
$s$-dependent) number $k_{*}$ of iterations, we reach
$\sigma_{k_{*}}>-s+\frac{\beta}{2}$ and deduce the validity of the regularity
bound (3.26).

5\. Final regularity: We now assume (3.26) and define the new sequence

$\sigma_{0}:=0,\quad\sigma_{k+1}:=\frac{\beta}{4}+\frac{\sigma_{k}}{2}\quad\Rightarrow\quad\sigma_{k}=\frac{\beta}{2}\left(1-\frac{1}{2^{k}}\right)\to\frac{\beta}{2}.$

We fix $\sigma=\sigma_{k}+\frac{\beta}{2}\in(0,\beta]$ and
$t=-s+\frac{\beta}{2}\in(-s,\sigma)$ in (3.25), and note that
$s+\sigma/2=s+\sigma_{k+1},$ $\sigma-t=s+\sigma_{k},$ and
$t+s=\frac{\beta}{2}$, to arrive at

$\|u\|_{\dot{B}^{s+\sigma_{k+1}}_{2,\infty}(\Omega)}^{2}\lesssim\|f\|_{B^{-s+\frac{\beta}{2}}_{2,1}(\Omega)}^{2}+\|f\|_{B^{-s+\frac{\beta}{2}}_{2,1}(\Omega)}\|u\|_{B^{s+\sigma_{k}}_{2,\infty}(\Omega)}+\|u\|_{B^{\frac{\beta}{2}}_{2,\infty}(\Omega)}\|u\|_{B^{s+\sigma_{k}}_{2,\infty}(\Omega)}.$

Using (3.26), we get

$\|u\|_{\dot{B}^{s+\sigma_{k+1}}_{2,\infty}(\Omega)}^{2}\leq\left(C_{1}\|f\|_{B^{-s+\frac{\beta}{2}}_{2,1}(\Omega)}+C_{2}\|u\|_{B^{s+\sigma_{k}}_{2,\infty}(\Omega)}\right)\|f\|_{B^{-s+\frac{\beta}{2}}_{2,1}(\Omega)},$

with constants $C_{1},C_{2}$ depending on $\Omega,d$ and $s$. We finally
proceed as in the proof of Theorem 3.1 (maximal Besov regularity) to show

$\|u\|_{\dot{B}^{s+\sigma_{k+1}}_{2,\infty}(\Omega)}\leq\Lambda\|f\|_{B^{-s+\frac{\beta}{2}}_{2,1}(\Omega)}$

for a finite number $\Lambda>0$. We trivially have
$\|u\|_{\dot{B}^{s+\sigma_{0}}_{2,\infty}(\Omega)}\leq\Lambda_{0}\|f\|_{B^{-s+\beta/2}_{2,1}(\Omega)}$
for some $\Lambda_{0}(d,s,\Omega)>0$. We define
$\Lambda_{k+1}:=(C_{1}+C_{2}\Lambda_{k})^{1/2}$ and realize that

$\|u\|_{\dot{B}^{s+\sigma_{k+1}}_{2,\infty}(\Omega)}\leq\Lambda_{k+1}\|f\|_{B^{-s+\frac{\beta}{2}}_{2,1}(\Omega)}.$

Since $\Lambda_{k+1}\leq\max\\{\Lambda_{0},(2C_{1}+C_{2})^{1/2}\\}=:\Lambda$,
passing to the limit $k\to\infty$, the desired estimate (3.24) follows
immediately. ∎

######  Remark 12 (truncated fractional Laplacians).

Estimate (3.24) holds whenever $\psi$ is locally $\beta$-Hölder continuous at
the origin. Consequently, it applies with $\beta=1$ to (linear) truncated
fractional Laplacians [16]

$\mathcal{A}u(x):=C(d,s,\delta)\int_{B_{\delta}(x)}\frac{u(x)-u(y)}{|x-y|^{d+2s}}\,dy,$

for which $\psi(\rho)=\chi_{[0,1]}(\rho)$.

In the same fashion as Corollary 3.1, the following maximal Sobolev regularity
holds for operators of the form (3.21) and, in particular, for linear
truncated fractional Laplacians and tempered fractional $p$-Laplacians.

######  Corollary 3.4 (maximal Sobolev regularity).

Under the hypothesis of either Corollary 3.3 for any $p\in(1,\infty)$ or of
Theorem 3.2 for $p=2$, for all $\varepsilon>0$ sufficiently small and
$q=\max\\{2,p\\}$ there holds

$\|u\|_{\widetilde{W}^{s+\frac{\beta}{q}-\varepsilon}_{p}(\Omega)}\lesssim\varepsilon^{-\frac{1}{p}}\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}^{\frac{q-p}{p-1}}\|f\|_{B^{-s+\frac{\beta}{q^{\prime}}}_{p^{\prime},1}(\Omega)}^{\frac{1}{q-1}}.$

##  4\. Approximation

As an application of our regularity estimates, we consider discretizations of
the problem (1.4) by means of the finite element method with piecewise linear
continuous functions. From now on, we assume that $G$ satisfies Hypothesis
2.1.

Let $h_{0}>0$; for $h\in(0,h_{0}]$, we let $\mathcal{T}_{h}$ denote a
triangulation of $\Omega$, i.e., $\mathcal{T}_{h}=\\{T\\}$ is a partition of
$\Omega$ into simplices $T$ of diameter $h_{T}$ and
$h=\max_{T\in{\mathcal{T}_{h}}}h_{T}$. We assume the family
$\\{{\mathcal{T}_{h}}\\}_{h>0}$ to be shape-regular, namely,

$\sigma:=\sup_{h>0}\max_{T\in{\mathcal{T}_{h}}}\frac{h_{T}}{\rho_{T}}<\infty,$

where $\rho_{T}$ is the diameter of the largest ball contained in $T$. We take
elements to be closed sets.

Let $\mathcal{N}_{h}$ be the set of interior vertices of ${\mathcal{T}_{h}}$,
$N$ be its cardinality and $\\{\varphi_{i}\\}_{i=1}^{N}$ be the standard
piecewise linear Lagrangian basis, with $\varphi_{i}$ associated to the node
$\texttt{x}_{i}\in\mathcal{N}_{h}$. With this notation, the set of discrete
functions is

$\widetilde{\mathbb{V}}_{h}:=\left\\{v:{\mathbb{R}^{d}}\to\mathbb{R}\colon
v\in C_{0}({\mathbb{R}^{d}}),\ v=\sum_{i=1}^{N}v_{i}\varphi_{i}\right\\},$

where $v$ is trivially extended by zero outside $\Omega$. It is clear that
$\widetilde{\mathbb{V}}_{h}\subset\widetilde{W}^{s}_{p}(\Omega)$ for all
$s\in(0,1),p\in(1,\infty)$. Therefore, we consider a direct finite element
discretization and seek $u_{h}\in\widetilde{\mathbb{V}}_{h}$ such that

(4.1)
$\iint_{{\mathbb{R}^{d}}\times{\mathbb{R}^{d}}}\widetilde{G}\left(x,y,\frac{u_{h}(x)-u_{h}(y)}{|x-y|^{s}}\right)\frac{(u_{h}(x)-u_{h}(y))(v_{h}(x)-v_{h}(y))}{|x-y|^{d+2s}}\,dx\,dy=\langle
f,v_{h}\rangle$

for all $v_{h}\in\widetilde{\mathbb{V}}_{h}$, where we recall that
$\widetilde{G}(x,y,\rho)=G_{\rho}(x,y,\rho)/\rho$. Clearly, $u_{h}$ solves
(4.1) if and only if it is the minimizer of the restriction of the convex
functional $\mathcal{F}$ from (1.1) over the linear space
$\widetilde{\mathbb{V}}_{h}$; existence of discrete solutions follows
immediately. Moreover, if we take $v_{h}=u_{h}$ in (4.1), then we immediately
obtain the discrete stability bound

(4.2)
$\|u_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}\lesssim\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}^{\frac{1}{p-1}}.$

###  4.1. Localization and interpolation estimates

The seminorm $|\cdot|_{W^{s}_{p}({\mathbb{R}^{d}})}$ is nonlocal and,
consequently, is non-additive with respect to domain partitions. To localize
it, we first define the star (or patch) of a set $A\in\Omega$ by

$S_{A}^{1}:=\bigcup\big{\\{}T\in{\mathcal{T}_{h}}\colon T\cap
A\neq\emptyset\big{\\}}.$

Given $T\in{\mathcal{T}_{h}}$, the star $S_{T}^{1}$ of $T$ is the first ring
of $T$ and the star $S_{T}^{2}$ of $S_{T}^{1}$ is the second ring of $T$. We
have the following localization estimate, that can be proved with the same
arguments as in [28, 29]:

(4.3)
$|v|_{W^{s}_{p}(\Omega)}^{p}\leq\sum_{T\in{\mathcal{T}_{h}}}\left[\int_{T}\int_{S_{T}^{1}}\frac{|v(x)-v(y)|^{p}}{|x-y|^{d+sp}}\;dy\;dx+C(d,\sigma)\frac{2^{p}}{sph_{T}^{sp}}\|v\|^{p}_{L^{p}(T)}\right]$

for all $v\in W^{s}_{p}(\Omega)$.

This localization of fractional-order seminorms implies that, in order to
prove global interpolation estimates in $W^{s}_{p}(\Omega)$, it suffices to
produce over the set of patches $\\{T\times
S_{T}^{1}\\}_{T\in{\mathcal{T}_{h}}}$ plus local, zero-order contributions.

We point out, however, that clearly if one wants to have a zero-extension norm
on the left hand side in (4.3), then interactions between $\Omega$ and
$\Omega^{c}$ must be accounted for in the right hand side. For that purpose,
following [12], we introduce the extended stars

$\widetilde{S}_{T}^{1}:=\left\\{\begin{array}[]{rl}S_{T}^{1}&\textrm{if
}T\cap\partial\Omega=\emptyset,\\\
B_{T}&\textrm{otherwise,}\end{array}\right.$

where $B_{T}=B(x_{T},Ch_{T})$ is the ball of center $x_{T}$ and radius
$Ch_{T}$, with $x_{T}$ being the barycenter of $T$, and $C=C(\sigma)$ a shape
regularity dependent constant such that $S_{T}^{1}\subset B_{T}$. The extended
second ring $\widetilde{S}_{T}^{2}$ of $T$ is given by

$\widetilde{S}_{T}^{2}:=\bigcup\big{\\{}\widetilde{S}_{T^{\prime}}^{1}\colon
T^{\prime}\in{\mathcal{T}_{h}},T^{\prime}\cap
S_{T}^{1}\neq\emptyset\big{\\}}.$

The localization of the $W^{s}_{p}({\mathbb{R}^{d}})$-seminorm reads [12,
Lemma 4.1]:

(4.4)
$\|v\|_{\widetilde{W}^{s}_{p}(\Omega)}^{p}=|v|_{W^{s}_{p}({\mathbb{R}^{d}})}^{p}\leq\sum_{T\in{\mathcal{T}_{h}}}\left[\int_{T}\int_{\widetilde{S}_{T}^{1}}\frac{|v(x)-v(y)|^{p}}{|x-y|^{d+sp}}\;dy\;dx+C(d,\sigma)\frac{2^{p}}{sph_{T}^{sp}}\|v\|^{p}_{L^{p}(T)}\right]$

for all $v\in\widetilde{W}^{s}_{p}(\Omega)$.

Our use of (4.4) will be restricted to $v$ being an interpolation error; in
such a case, $v$ has vanishing means over elements and thus we can estimate
the scaled $L^{p}$ contributions in terms of local $W^{s}_{p}$ seminorms by
using Poincaré inequalities. We consider a suitable (such as Clément or Scott-
Zhang) quasi-interpolation operator,

$\Pi_{h}:\widetilde{W}^{s}_{p}(\Omega)\to\mathbb{V}_{h},$

which is stable in $W^{s}_{p}(\Omega)$ and for which one can prove the
following local approximation estimates (see, for example, [2, 13, 20]):

(4.5) $\begin{split}&\|v-\Pi_{h}v\|_{L^{p}(T)}\leq
C\,h_{T}^{pr}|v|_{W^{r}_{p}(S_{T}^{1})}^{p},\\\
&\int_{T}\int_{\widetilde{S}_{T}^{1}}\frac{|(v-\Pi_{h}v)(x)-(v-\Pi_{h}v)(y)|^{p}}{|x-y|^{d+sp}}\,dy\,dx\leq
C\,h_{T}^{p(r-s)}|v|_{W^{r}_{p}(\widetilde{S}_{T}^{2})}^{p},\end{split}$

for all $T\in{\mathcal{T}_{h}}$, $s\in(0,1)$, $r\in(s,2]$, $v\in
W^{r}_{p}(\widetilde{S}_{T}^{2})$, where $C=C(\Omega,d,s,\sigma,right)$.

Combining (4.4) and (4.5), we deduce localized interpolation error estimates.

######  Proposition 4.1 (localized interpolation estimates).

Let $s\in(0,1)$, $p\in(1,\infty)$, $r\in(s,2]$, and
$\Pi_{h}:\widetilde{W}^{s}_{p}(\Omega)\to\mathbb{V}_{h}$ be a quasi-
interpolation operator as above. If $v\in\widetilde{W}^{r}_{p}(\Omega)$, then

(4.6) $\|v-\Pi_{h}v\|_{\widetilde{W}^{s}_{p}(\Omega)}^{p}\leq
C\left(\sum_{T\in{\mathcal{T}_{h}}}h_{T}^{p(r-s)}|v|_{W^{r}_{p}(\widetilde{S}_{T}^{2})}^{p}\right)^{\frac{1}{p}},$

where $C=C(\Omega,d,s,\sigma,right)$.

###  4.2. Error estimates in $\widetilde{W}^{s}_{p}(\Omega)$

We next derive some error estimates for the finite element solutions discussed
in §4. We borrow techniques from the finite element analysis of classical
(local) quasi-linear problems. The technique presented in [31] or [19, §5.3]
exploits the continuity and monotonicity of the operator, but not the fact
that $u$ and $u_{h}$ solve respective minimization problems. We obtain
enhanced rates by adapting an approach by Chow [18] for the classical
$p$-Laplacian.

######  Theorem 4.1 (error estimates).

Let $\Omega$ be a bounded Lipschitz domain, assume that $G$ satisfies
Hypothesis 2.1, let $p\in(1,\infty)$, $\beta\in(0,1]$ be as in such
assumptions and let $p^{\prime}=\frac{p}{p-1}$. Assume $f\in
B^{-s+\gamma}_{p^{\prime},1}(\Omega)$, where
$\gamma=\max\\{\beta/p^{\prime},\beta/2\\}$. Let $u$ and $u_{h}$ be the
respective solutions of (2.20) and (4.1). Then, if $p\in(1,2]$ it holds that

(4.7)
$\|u-u_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}\lesssim\inf_{v_{h}\in\mathbb{V}_{h}}\|u-v_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}^{p/2}\lesssim
h^{\frac{\beta p}{4}}|\log h|^{\frac{1}{2}}.$

On the other hand, if $p\in[2,\infty)$, we have the error bound

(4.8)
$\|u-u_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}\lesssim\inf_{v_{h}\in\mathbb{V}_{h}}\|u-v_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}^{2/p}\lesssim
h^{\frac{2\beta}{p^{2}}}|\log h|^{\frac{2}{p^{2}}}.$

###### Proof.

For any $v\in\widetilde{W}^{s}_{p}(\Omega)$, using either (2.17) or (2.16), we
write

$\displaystyle\mathcal{F}(v)-\mathcal{F}(u)$
$\displaystyle=\int_{0}^{1}\langle\mathcal{F}^{\prime}(u+t(v-u))-\mathcal{F}^{\prime}(u),v-u\rangle\,dt$
$\displaystyle=\int_{0}^{1}\langle\mathcal{A}(u+t(v-u))-\mathcal{A}u,t(v-u)\rangle\,\frac{dt}{t}$
$\displaystyle\geq\left\\{\begin{array}[]{rl}\frac{C}{p}\|u-v\|_{\widetilde{W}^{s}_{p}(\Omega)}^{2}\left(\|u\|_{\widetilde{W}^{s}_{p}(\Omega)}+\|u-v\|_{\widetilde{W}^{s}_{p}(\Omega)}\right)^{p-2},&\mbox{if
}p\in(1,2],\\\
\frac{\alpha}{p}\|u-v\|_{\widetilde{W}^{s}_{p}(\Omega)}^{p},&\mbox{if
}p\in[2,\infty).\end{array}\right.$

In addition, if we use either (2.18) or (2.19), we obtain

$\displaystyle\mathcal{F}(v)-\mathcal{F}(u)$
$\displaystyle=\int_{0}^{1}\langle\mathcal{A}(u+t(v-u))-\mathcal{A}u,t(v-u)\rangle\,\frac{dt}{t}$
$\displaystyle\leq\left\\{\begin{array}[]{rl}\frac{C}{p}\|u-v\|_{\widetilde{W}^{s}_{p}(\Omega)}^{p},&\mbox{if
}p\in(1,2],\\\
\frac{C}{p}(\|u\|_{\widetilde{W}^{s}_{p}(\Omega)}+\|u-v\|_{\widetilde{W}^{s}_{p}(\Omega)})^{p-2}\|u-v\|_{\widetilde{W}^{s}_{p}(\Omega)}^{2},&\mbox{if
}p\in[2,\infty).\end{array}\right..$

The proof now follows easily. Indeed, for any
$v_{h}\in\mathbb{V}_{h}\subset\widetilde{W}^{s}_{p}(\Omega)$, we have, for
$p\in(1,2]$,

$\displaystyle
c(\|u\|_{\widetilde{W}^{s}_{p}(\Omega)}+\|u-u_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)})^{p-2}\|u-u_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}^{2}$
$\displaystyle\leq\mathcal{F}(u_{h})-\mathcal{F}(u)$
$\displaystyle\leq\mathcal{F}(v_{h})-\mathcal{F}(u)\leq
c\|u-v_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}^{p}.$

By the stability estimates (2.21) and (4.2), we have

$\|u\|_{\widetilde{W}^{s}_{p}(\Omega)}+\|u-u_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}\lesssim\|f\|_{W^{-s}_{p^{\prime}}(\Omega)}^{\frac{1}{p-1}}\lesssim\|f\|_{B^{-s+\gamma}_{p^{\prime},1}(\Omega)}^{\frac{1}{p-1}}$

and therefore

$\|u-u_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}\lesssim\|u-v_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}^{p/2},\quad\forall
v_{h}\in\mathbb{V}_{h},$

which proves the first inequality in (4.7).

Similar considerations yield, for $p\in[2,\infty)$,

$\|u-u_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}^{p}\lesssim(\|u\|_{\widetilde{W}^{s}_{p}(\Omega)}+\|u-v_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)})^{p-2}\|u-v_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}^{2},\quad\forall
v_{h}\in\mathbb{V}_{h},$

and thus the first inequality in (4.8) holds.

We now set $v_{h}={\Pi_{h}}u$, use the stability of $\Pi_{h}$ in
$\widetilde{W}^{s}_{p}(\Omega)$, the quasi-interpolation estimate (4.6) and
the regularity bounds (3.18) or (3.19), to conclude:

$\|u-u_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}\lesssim\begin{cases}\|u-v_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}^{p/2}\lesssim
h^{\frac{\beta p}{4}-\frac{\varepsilon
p}{2}}\,\varepsilon^{-\frac{1}{2}},&\mbox{ if }p\in(1,2],\\\
\|u-v_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}^{2/p}\lesssim
h^{\frac{2\beta}{p^{2}}-\frac{2\varepsilon}{p}}\,\varepsilon^{-\frac{2}{p^{2}}},&\mbox{
if }p\in[2,\infty).\end{cases}$

for $\varepsilon>0$ sufficiently small. The result follows by setting
$\varepsilon=|\log h|^{-1}$. ∎

##  5\. Computational exploration

This section presents several numerical experiments for the Dirichlet problem
(1.4). We use the finite element discretization discussed in Section 4 on
either quasi-uniform or shape-regular graded meshes ${\mathcal{T}_{h}}$ with
grading parameter $\mu\geq 1$ satisfying

(5.1)
$h_{T}\approx\left\\{\begin{array}[]{ll}C(\sigma)h^{\mu},&\overline{T}\cap\partial\Omega\neq\emptyset\\\
C(\sigma)h\;\textrm{dist}(T,\partial\Omega)^{(\mu-1)/\mu},&\overline{T}\cap\partial\Omega=\emptyset\end{array}\right.$

for every $T\in{\mathcal{T}_{h}}$. We refer to [10] for further details on
this grading strategy and additional computational experiments.

###  5.1. Fractional $(p,s)$-Laplacians

Along this section, we consider the energy minimization problem (1.1) with
$G(x,y,\rho)=\frac{C_{d,s,p}}{2p}|\rho|^{p}$ for $p\in(1,\infty)$, that gives
rise to the fractional $(p,s)$-Laplace operator (1.5).

######  Example 5.1 (boundary behavior).

We let $\Omega=(-0.5,0.5)^{2}\setminus[0,0.5)\times(-0.5,0]$ be an $L$-shaped
domain and $f=1$, and investigate the boundary behavior of numerical
solutions. From the analytical results for the linear problem ($p=2$), we
expect the solution to have a boundary behavior of the type

$u(x)\approx\textrm{dist}(x,\partial\Omega)^{\alpha(s,p)}.$

We estimate the power $\alpha(s,p)$ near different points on $\partial\Omega$:
the mid-point of one of the edges $(0,0.5)$, a convex corner $(-0.5,0.5)$ and
the reentrant corner $(0,0)$. We compute the numerical solutions on the graded
mesh with $\mu=2$ and $9467$ free nodes, and fit the power $\alpha(s,p)$ using
mesh points near the boundary points mentioned above, specifically along the
normal direction near the edge mid-point and along the bisectors of the angles
near the corners. We report the results we obtain in Table 1 and Table 2 for
$p=3$ and $p=1.75$, respectively.

Despite the limited mesh resolution, near the edge mid-point, we observe
$\alpha(s,p)\approx s$ for both $p=3$ and $p=1.75$. This is consistent with
the behavior shown in [34, 35] for domains satisfying an exterior ball
condition. In addition, we notice that $\alpha(s,p)>s$ near the convex corner,
$\alpha(s,p)<s$ near the reentrant corner, and the deviation from $s$ is
larger when $p=1.75$ compared to $p=3$.

Value of $s$ | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9
---|---|---|---|---|---|---|---|---|---
Edge mid-point | 0.11 | 0.21 | 0.31 | 0.41 | 0.51 | 0.61 | 0.70 | 0.80 | 0.90
Convex corner | 0.09 | 0.22 | 0.35 | 0.49 | 0.64 | 0.79 | 0.94 | 1.11 | 1.29
Reentrant corner | 0.06 | 0.16 | 0.24 | 0.31 | 0.39 | 0.46 | 0.54 | 0.62 | 0.70

Table 1. Example (5.1): Exponents $\alpha=\alpha(s,p)$ of boundary asymptotics
$u(x)\approx\textrm{dist}(x,\partial\Omega)^{\alpha}$ for $p=3$, different
values of fractional order $s$, and three qualitatively distinct boundary
points.

Value of $s$ | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9
---|---|---|---|---|---|---|---|---|---
Edge mid-point | 0.10 | 0.21 | 0.31 | 0.41 | 0.51 | 0.60 | 0.70 | 0.79 | 0.89
Convex corner | 0.11 | 0.27 | 0.46 | 0.66 | 0.86 | 1.07 | 1.29 | 1.51 | 1.76
Reentrant corner | 0.07 | 0.12 | 0.17 | 0.23 | 0.29 | 0.35 | 0.42 | 0.49 | 0.57

Table 2. Example (5.1): Exponents $\alpha=\alpha(s,p)$ of boundary asymptotics
$u(x)\approx\textrm{dist}(x,\partial\Omega)^{\alpha}$ for $p=1.75$, different
values of fractional order $s$, and three qualitatively distinct boundary
points.

######  Example 5.2 (convergence rates).

Consider the square domain $\Omega=(-0.5,0.5)^{2}$ and $f=1$. Since $f$ is
smooth, Corollary 3.1 (maximal Sobolev regularity) gives
$u\in\widetilde{W}^{s+1/p-\varepsilon}_{p}(\Omega)$ for $p\geq 2$ and
$u\in\widetilde{W}^{s+1/2-\varepsilon}_{p}(\Omega)$ for $1<p\leq 2$. We
compute numerical solutions for $p=3$, $p=1.75$ and different values of $s$ on
quasi-uniform meshes, and examine convergence rates in the energy norm. Since
the exact solutions $u$ are unknown, we use
$\|u_{h}-u_{h/2}\|_{\widetilde{W}^{s}_{p}(\Omega)}$ as a proxy for
$\|u_{h}-u\|_{\widetilde{W}^{s}_{p}(\Omega)}$. Table 3 summarizes our
findings.

Value of $s$ | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9
---|---|---|---|---|---|---|---|---|---
$p=3$ | 0.326 | 0.335 | 0.333 | 0.329 | 0.328 | 0.329 | 0.335 | 0.357 | 0.494
$p=1.75$ | 0.558 | 0.552 | 0.555 | 0.561 | 0.569 | 0.583 | 0.607 | 0.658 | 0.790

Table 3. Example 5.2: Convergence rates on uniform meshes for $p=1.75,3$ and
different values of $s$. They indicate that the theoretical rates in Theorem
4.1 (error estimates) might be suboptimal.

The rates are approximately $1/p\approx 0.33$ for $p=3$ except for the case
$s=0.9$, where we believe the discrepancy is due to the proxy solution not
being sufficiently refined in comparison to the rest of the experiments.
Although the rate $1/p$ is larger than the theoretical rate $2/p^{2}$ of (4.8)
in Theorem 4.1 (error estimates), we point out that it is consistent with the
best approximation error

$\inf_{v_{h}\in\mathbb{V}_{h}}\|u-v_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}\lesssim
h^{1/p}|\log h|^{1/p}.$

This indicates that, instead of the regularity of $u$, the suboptimal rates in
the theory might be a consequence of the suboptimal power $2/p$ in the error
estimate of (4.8)

$\|u-u_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}\lesssim\inf_{v_{h}\in\mathbb{V}_{h}}\|u-v_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}^{2/p}.$

Similarly, for $p=1.75$, we observe that the rates are approximately
$1/p\approx 0.57$ except for $s=0.7,0.8,0.9$, where we believe the meshes are
not fine enough to deliver accurate rates. This indicates that instead of the
regularity $u\in\widetilde{W}^{s+1/2-\varepsilon}_{p}(\Omega)$ proved in
Corollary 3.1, the solution $u$ in this example might satisfy
$u\in\widetilde{W}^{s+1/p-\varepsilon}_{p}(\Omega)$. In addition, the power
$p/2$ in the error estimate (4.7) of Theorem 4.1,

$\|u-u_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}\lesssim\inf_{v_{h}\in\mathbb{V}_{h}}\|u-v_{h}\|_{\widetilde{W}^{s}_{p}(\Omega)}^{p/2},$

might not be optimal as well. These two reasons together lead to the
theoretical suboptimal rate $p/4$ of (4.7).

###  5.2. Linear operators with finite horizon

We consider the operator $\mathcal{A}_{\delta}u$ defined in (3.21) and let
$p=2$, that gives rise to the linear fractional Laplacian with finite horizon
and variable diffusivity.

######  Example 5.3 (truncated fractional Laplacian in $1$D).

Consider $\Omega=(-1,1)$, $p=2,\delta=0.2,f=\frac{\delta^{2-2s}}{1-s}$ and
$\psi=\chi_{[0,1]}$ where $\chi_{I}$ is the characteristic function of $I$.
The resulting operator is the truncated fractional Laplacian we discussed in
Remark 12. We first compute numerical solutions for different $s$ using a mesh
graded according to (5.1) with $h=2^{-12},$ $\mu=2$ to investigate the
boundary behavior of solutions. Figure 5.1 displays the solutions we obtained
for several values of $s$.

Figure 5.1. Example 5.3: Numerical solutions of linear operator with finite
horizon for different values of $s$.

Assuming the solutions have an algebraic boundary behavior

$u(x)\approx\textrm{dist}(x,\partial\Omega)^{\alpha(s)},$

we estimate $\alpha(s)$ numerically in Table 4. We clearly observe
$\alpha(s)\approx s$; this is the same as for the fractional Laplacian
$(-\Delta)^{s}$.

Value of $s$ | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9
---|---|---|---|---|---|---|---|---|---
$\alpha(s)$ | 0.10 | 0.19 | 0.29 | 0.39 | 0.49 | 0.60 | 0.70 | 0.80 | 0.90

Table 4. Example 5.3: Boundary exponents $\alpha(s)$ for different values of
$s$. They confirm $u(x)\approx\textrm{dist}(x,\partial\Omega)^{s}$.

Next, we measure convergence rates for different $s$ on uniform meshes for $h$
from $2^{-8}$ to $2^{-12}$. Since we do not know a closed formula for the
solution $u$, we use $\|u_{h}-u_{h/2}\|_{\widetilde{H}^{s}(\Omega)}$ as a
proxy for $\|u_{h}-u\|_{\widetilde{H}^{s}(\Omega)}$ and present the rates in
terms of $h$ in Table 5. We observe the convergence rates are about $0.5$ for
all $s$, in agreement with the regularity
$u\in\widetilde{H}^{s+\frac{1}{2}-\varepsilon}(\Omega)$ proved in Corollary
3.4 (maximal Sobolev regularity) and a standard best-approximation argument.

Value of $s$ | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9
---|---|---|---|---|---|---|---|---|---
Rate | 0.530 | 0.508 | 0.500 | 0.498 | 0.497 | 0.498 | 0.498 | 0.498 | 0.501

Table 5. Example 5.3: Optimal convergence rates on uniform meshes for
different values of $s$.

## References

  * [1] H. Abels and G. Grubb.  Fractional-order operators on nonsmooth domains.  Journal of the London Mathematical Society, 107(4):1297–1350, 2023\.
  * [2] G. Acosta and J. Borthagaray.  A fractional Laplace equation: regularity of solutions and finite element approximations.  SIAM J. Numer. Anal., 55(2):472–495, 2017.
  * [3] R. Adams and J. Fournier.  Sobolev spaces.  Elsevier, 2003.
  * [4] G. Akagi, G. Schimperna, A. Segatti, and L. V. Spinolo.  Quantitative estimates on localized finite differences for the fractional Poisson problem, and applications to regularity and spectral stability.  Commun. Math. Sci., 16(4):913–961, 2018.
  * [5] C. Atkinson and C. Jones.  Similarity solutions in some non-linear diffusion problems and in boundary-layer flow of a pseudo-plastic fluid.  Quart. J. Mech. Appl. Math., 27(2):193–211, 1974.
  * [6] J. Barrett and W. Liu.  Finite element approximation of the $p$-Laplacian.  Math. Comp., 61(204):523–537, 1993.
  * [7] J. Benedikt, P. Girg, L. Kotrla, and P. Takáč.  Origin of the $p$-Laplacian and A. Missbach.  Electron. J. Differ. Equ., 2018:17, 2018.  Id/No 16.
  * [8] D. A. Benson, S. W. Wheatcraft, and M. M. Meerschaert.  Application of a fractional advection-dispersion equation.  Water resources research, 36(6):1403–1412, 2000.
  * [9] J. Bergh and J. Löfström.  Interpolation spaces: an introduction.  Springer-Verlag, Berlin, 1976.
  * [10] J. Borthagaray, W. Li, and R. Nochetto.  Fractional elliptic problems on Lipschitz domains: Regularity and approximation.  A3N2M: Approximation, Applications, and Analysis of Nonlocal, Nonlinear Models Proceedings of the 50th John H. Barrett Memorial Lectures, 2023.
  * [11] J. Borthagaray and R. Nochetto.  Besov regularity for the Dirichlet integral fractional Laplacian in Lipschitz domains.  J. Funct. Anal., 284(6):109829, 2023.
  * [12] J. Borthagaray and R. Nochetto.  Constructive approximation on graded meshes for the integral fractional Laplacian.  Constr. Approx., 57:463–487, 2023.
  * [13] J. Borthagaray, R. Nochetto, and A. 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.
  * [14] J. Bourgain, H. Brezis, and P. Mironescu.  Another look at Sobolev spaces.  In Optimal Control and Partial Differential Equations, pages 439–455, 2001.
  * [15] L. Brasco and E. Lindgren.  Higher Sobolev regularity for the fractional $p$-Laplace equation in the superquadratic case.  Adv. Math., 304:300–354, 2017.
  * [16] O. Burkovska and M. Gunzburger.  Regularity analyses and approximation of nonlocal variational equality and inequality problems.  J. Math. Anal. Appl., 478(2):1027–1048, 2019.
  * [17] L. Caffarelli.  Non-local diffusions, drifts and games.  In Nonlinear partial differential equations. The Abel symposium 2010\. Proceedings of the Abel symposium, Oslo, Norway, September 28–October 2, 2010, pages 37–52. Berlin: Springer, 2012.
  * [18] S.-S. Chow.  Finite element error estimates for non-linear elliptic equations of monotone type.  Numer. Math., 54(4):373–393, 1989.
  * [19] P. Ciarlet.  The finite element method for elliptic problems, volume 40.  Siam, 2002.
  * [20] P. Ciarlet, Jr.  Analysis of the Scott-Zhang interpolation in the fractional order Sobolev spaces.  J. Numer. Math., 21(3):173–180, 2013.
  * [21] R. Cont and P. Tankov.  Financial modelling with jump processes. Chapman Hall/CRC Financ. Math. Ser. Boca Raton, FL: Chapman and Hall/CRC, 2004.
  * [22] F. del Teso, D. Gómez-Castro, and J. Vázquez.  Three representations of the fractional $p$-Laplacian: semigroup, extension and Balakrishnan formulas.  Fract. Calc. Appl. Anal., 24(4):966–1002, 2021.
  * [23] F. del Teso, M. Medina, and P. Ochoa.  Higher-order asymptotic expansions and finite difference schemes for the fractional $p$-laplacian.  arXiv preprint arXiv:2303.02502, 2023.
  * [24] W. Deng, B. Li, W. Tian, and P. Zhang.  Boundary problems for the fractional and tempered fractional operators.  Multiscale Model. Simul., 16(1):125–149, 2018.
  * [25] J. Diaz and F. De Thelin.  On a nonlinear parabolic problem arising in some models related to turbulent flows.  SIAM J. Math. Anal., 25(4):1085–1111, 1994.
  * [26] J. I. Díaz, G. Hetzer, and L. Tello.  An energy balance climate model with hysteresis.  Nonlinear Anal., Theory Methods Appl., Ser. A, Theory Methods, 64(9):2053–2074, 2006.
  * [27] B. Dyda and M. Kassmann.  On weighted Poincaré inequalities.  Ann. Acad. Sci. Fenn. Math, 38(2):721–726, 2013.
  * [28] B. Faermann.  Localization of the Aronszajn-Slobodeckij norm and application to adaptive boundary element methods. I. The two-dimensional case.  IMA J. Numer. Anal., 20(2):203–234, 2000.
  * [29] B. Faermann.  Localization of the Aronszajn-Slobodeckij norm and application to adaptive boundary element methods. II. The three-dimensional case.  Numer. Math., 92(3):467–499, 2002.
  * [30] J. Fernández Bonder, A. Salort, and H. Vivas.  Interior and up to the boundary regularity for the fractional g-laplacian: The convex case.  Nonlinear Analysis, 223:113060, 2022.
  * [31] R. Glowinski and A. Marroco.  Sur l’approximation, par éléments finis d’ordre un, et la résolution, par pénalisation-dualité d’une classe de problèmes de Dirichlet non linéaires.  ESAIM: Mathematical Modelling and Numerical Analysis-Modélisation Mathématique et Analyse Numérique, 9(R2):41–76, 1975.
  * [32] P. Grisvard.  Elliptic problems in nonsmooth domains, volume 24 of Monographs and Studies in Mathematics.  Pitman (Advanced Publishing Program), Boston, MA, 1985.
  * [33] G. Grubb.  Fractional Laplacians on domains, a development of Hörmander’s theory of $\mu$-transmission pseudodifferential operators.  Adv. Math., 268:478–528, 2015.
  * [34] A. Iannizzotto, S. Mosconi, and M. Squassina.  Global Hölder regularity for the fractional $p$-Laplacian.  Rev. Mat. Iberoam., 32(4):1353–1392, 2016.
  * [35] A. Iannizzotto, S. J. N. Mosconi, and M. Squassina.  Fine boundary regularity for the degenerate fractional $p$-Laplacian.  J. Funct. Anal., 279(8):108659, 54, 2020.
  * [36] T. A. M. Langlands, B. I. Henry, and S. L. Wearne.  Fractional cable equation models for anomalous electrodiffusion in nerve cells: Infinite domain solutions.  J. Math. Biol., 59(6):761–808, 2009.
  * [37] T. A. M. Langlands, B. I. Henry, and S. L. Wearne.  Fractional cable equation models for anomalous electrodiffusion in nerve cells: finite domain solutions.  SIAM J. Appl. Math., 71(4):1168–1203, 2011.
  * [38] S. Z. Levendorskiĭ.  Pricing of the American put under Lévy processes.  Int. J. Theor. Appl. Finance, 7(3):303–335, 2004.
  * [39] V. Maz’ya and T. Shaposhnikova.  On the Bourgain, Brezis, and Mironescu theorem concerning limiting embeddings of fractional Sobolev spaces.  Journal of Functional Analysis, 195(2):230 – 238, 2002.
  * [40] W. McLean.  Strongly elliptic systems and boundary integral equations.  Cambridge university press, 2000.
  * [41] R. Metzler and J. Klafter.  The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics.  Journal of Physics A: Mathematical and Theoretical, 37(31):R161–R208, 2004.
  * [42] S. Mosconi and M. Squassina.  Recent progresses in the theory of nonlinear nonlocal problems.  In Bruno Pini Mathematical Analysis Seminar 2016, volume 7 of Bruno Pini Math. Anal. Semin., pages 147–164. Univ. Bologna, Alma Mater Stud., Bologna, 2016.
  * [43] G. Savaré.  Regularity results for elliptic equations in Lipschitz domains.  J. Funct. Anal., 152(1):176–201, 1998.
  * [44] L. Tartar.  Interpolation non linéaire et régularité.  J. Functional Analysis, 9:469–489, 1972.
  * [45] H. Triebel.  Theory of Function Spaces.  Modern Birkhäuser Classics. Springer Basel, 2010.
  * [46] E. Valdinoci.  From the long jump random walk to the fractional Laplacian.  Bol. Soc. Esp. Mat. Apl., S$\vec{\text{e}}$MA, 49:33–44, 2009\.
  * [47] M. I. Višik and G. I. Èskin.  Convolution equations in a bounded region.  Uspehi Mat. Nauk, 20(3 (123)):89–152, 1965.  English translation in Russian Math. Surveys, 20:86-151, 1965.