We extend the classical Bernstein technique to the setting of integro-differential operators. As a consequence, we provide first and one-sided second derivative estimates for solutions to fractional equations, including some convex fully nonlinear equations of order smaller than two —for which we prove uniform estimates as their order approaches two. Our method is robust enough to be
applied to some Pucci-type extremal equations and to obstacle problems for fractional operators, although several of the results are new even in the linear case.
We also raise some intriguing open questions, one of them concerning the “pure” linear fractional Laplacian, another one
being the validity of one-sided second derivative estimates for Pucci-type convex equations associated to linear operators with general kernels.

[2010]35R11, 47G20, 35B65, 35B50

§ INTRODUCTION

The Bernstein technique is a powerful tool to establish derivative estimates, through the use of auxiliary functions and the maximum principle, for solutions of elliptic equations. The goal of this paper is to extend this method to the setting of fractional equations. Up to our knowledge, this is done in the current work for the first time, even in the linear case.
The technique will allow us to establish first and one-sided second derivative estimates for a large class of integro-differential equations, including some fully nonlinear equations of order smaller than two.

Fractional diffusions arise in classical models, notably in the description of water waves, of atom dislocations in crystals, and in the displacement of
an elastic membrane on a thin obstacle. These problems can be efficiently attacked
by transforming them into a fractional setting on a lower dimensional (or boundary) object.
More recent models include energy transfer in nanotubes,
plasma physics, price oscillations in stock markets, and biological dispersals
in sparse environments —see e.g. <cit.> and references

The Bernstein technique —as introduced in the local case by Bernstein
himself <cit.>—
relies on considering some auxiliary functions
which involve the solution itself, its derivatives, and suitable cutoff functions. In view of certain equations (inequalities, rather) satisfied by the auxiliary functions and thanks to the maximum principle, they allow to estimate higher derivatives of the solution in terms of lower order ones,
by paying a price in the size of the reference domain.

Let us recall this procedure with the simplest example in local equations, the Laplace operator. This will serve as a preparation for the nonlocal framework.
Given two functions $u$ and $\eta$ ($u$ must be thought as the solution of an equation, while $\eta$ will be a cutoff), both smooth enough, consider the auxiliary function

\begin{equation}\label{FACIL0}
\varphi:= \eta^2(\partial_e u)^2+\sigma u^2,
\end{equation}

where $e\in\R^n$ is a unit vector, $|e|=1$, and $\sigma>0$ is a constant. We have that

\begin{equation}\label{PREFA}\begin{split}
-\Delta\varphi \;&=\;
-2|\nabla \eta|^2 (\partial_e u)^2-
2\eta \Delta\eta (\partial_e u)^2-
8\eta \nabla\eta \cdot \nabla (\partial_e u)\partial_e u\\
2\eta^2|\nabla\partial_e u|^2-
2\eta^2\partial_e u\, \partial_e \Delta u
-2 \sigma |\nabla u|^2-
2 \sigma u \Delta u.\end{split}\end{equation}

Since, by the Cauchy-Schwarz inequality

$$ \big|8\eta \nabla\eta \cdot \nabla(\partial_e u)\partial_e u\big|\le
2\eta^2\,|\nabla\partial_e u|^2
+8|\nabla\eta|^2 (\partial_e u)^2,

(<ref>) yields

\begin{equation*}\begin{split}
-\Delta\varphi \;&\le\;
6|\nabla \eta|^2 (\partial_e u)^2-
2\eta \Delta\eta (\partial_e u)^2-2 \sigma |\nabla u|^2
+2\eta^2\partial_e u\, (-\Delta )\partial_eu
+2 \sigma u (-\Delta u).
\end{split}\end{equation*}

In particular, by choosing $\sigma\ge C_n \|\eta\|^2_{C^2(\R^n)}$ for an appropriate constant $C_n$ depending on $n$ (more specifically, on the precise way in which the $C^2$-norm of a function in $\R^{n}$ is defined),
we obtain that

\begin{equation}\label{FACIL}\begin{split}
-\Delta\left( \eta^2(\partial_e u)^2+\sigma u^2\right) \;\le\;
2\eta^2\partial_e u\, (-\Delta )\partial_eu
+2 \sigma u \,(-\Delta u) \qquad \text{ if } \sigma\ge C_n \|\eta\|^2_{C^2(\R^n)}.
\end{split}\end{equation}

This is a “clean”, key inequality satisfied by any function $u$ (not necessarily a solution of an equation).

Now, if we assume the function $u$
to be harmonic, say $-\Delta u=0$ in $B_1\subset\R^n$, we deduce from (<ref>)
that $-\Delta\varphi\le0$ in $B_1$. If, in addition, we take the function $\eta\in C^\infty_c(B_1)$ to have compact support in $B_1$ and to satisfy $\eta=1$ in $B_{1/2}$,
the maximum principle for subharmonic functions ensures that $\varphi$
attains its maximum along $\partial B_1$. As a consequence,

\sup_{B_{1/2}} (\partial_e u)^2\le
\sup_{B_{1/2}}\varphi\le
\sup_{B_1}\varphi=\sup_{\partial B_1}\varphi= \sigma\sup_{\partial B_1}u^2\le \sigma \|u\|_{L^\infty(B_1)}^2,$$

thus yielding an explicit interior gradient estimate for the solution $u$. As we will see, simple variations
of this method, in which higher derivatives are taken into account
within the auxiliary function, lead to higher order estimates as well.

In spite of its rather elementary flavor, the Bernstein method is a powerful nonvariational tool that finds applications in several contexts and
for a large number of equations. The quadratic auxiliary function above (which is the one that we will consider within the nonlocal setting) finds applications even for second order fully nonlinear uniformly elliptic equations; see the monograph <cit.>.
More sophisticated auxiliary functions (with other nonlinear dependences on $u$ and $\partial_{e}u$) lead to gradient estimates for the prescribed mean curvature equation; see e. g. <cit.>
and references therein.

Instead, to the best of our knowledge, the Bernstein method has not been yet considered in relation with
fractional and integro-differential equations. We undertake this task here for quadratic auxiliary functions as in (<ref>).
In accordance with our proofs and results, we must merge
the operators that we treat into two
categories. The first one consists of
equations which admit a local extension in one more
dimension, as it is the case of the fractional Laplacian (the generator of the simplest $n$-dimensional Lévy jump process).
This class will also include several types of nonlocal convex fully nonlinear equations built from linear fractional operators. They arise
in optimization and in game theory, when players wish
to strategically optimize some diffusion parameters.
Finally, the class will also cover superpositions of fractional operators of different orders,
describing phenomena in which more than a single diffusion regime takes place.

For this first category of equations, we establish first and one-sided second derivative estimates.
While first derivative estimates have already been obtained by other methods in <cit.> and <cit.> (see a more detailed description below), our one-sided second derivative bounds (also called semiconcavity bounds) are new and somehow surprising. Indeed, since the order of the operators will be, in general, smaller than two, one should not expect a regularity theory up to the second order. Still, we will be able to prove one-sided second derivative estimates for some convex fully nonlinear fractional equations and, furthermore,
with estimates which will be uniform as the order of the operator approaches two.

Previously, one-sided second derivative bounds for fractional problems had been proved only for the thin obstacle problem, by Athanasopoulos and Caffarelli <cit.>.
We will address their result in Corollary <ref> below. In <cit.>
the Bernstein technique was already used, but with a different, less flexible, auxiliary function than in the local theory or in the current work. Their auxiliary function is linear in the second derivatives, while ours is quadratic and thus, as we will see, it has already allowed for applications to more general situations in the thin obstacle problem.

Our second category of equations consists of linear integro-differential
operators with general kernels, as well as fully nonlinear operators
built from them. In this case, the kernels are not rotationally invariant, in
general, and thus this setting
takes into account also anisotropic environments. In this
situation, a local extension problem will not be available and the
Bernstein technique becomes significantly more involved. Indeed, the
nonlocal contributions prevent the easy computational simplifications
that occur in the local case[Probably these additional
difficulties arising in the nonlocal setting
are not only a technical issue: as a matter of fact, in striking
contrast with the classical case,
the data at infinity can produce arbitrary modifications on the solutions
of a nonlocal equation; see <cit.>.] —essentially those
shown above for the Laplace operator— and the exact analogue of
inequality (<ref>) will remain, for general kernels, as an open
question. However, one of the main contributions of this paper is to
establish a similar inequality (Theorem <ref>) that,
despite including an error term, will still lead to first derivative
estimates. Instead, the validity of one-sided second derivative bounds for
general kernels remains as an open question.

We will raise some intriguing open problems in these respects. They will refer, thus, to the inequalities satisfied by the auxiliary functions in the
Bernstein technique, and will concern either the fractional Laplacian or operators with general kernels.
Some of them will deal with sharp inequalities that we are able to prove only using the local extension —and, hence, inequalities that we know to be true. We will ask for “downstairs” proofs of them, with the hope of being able to extend them to general kernels.

As already mentioned, first derivative bounds have already been proved for large classes of fully nonlinear integro-differential equations, using different methods than ours, by Caffarelli and Silvestre <cit.> and by Barles, Chasseigne, Ciomaga, and Imbert <cit.>.
The seminal work <cit.> established a $C^{1+\alpha}$ bound for a large class of fully nonlinear integro-differential equations that includes Isaacs-type equations made from uniformly elliptic linear operators with general kernels in the class ${\mathcal{L}}_{1}$. Their proof relies on ABP-type and Harnack inequalities, and thus it is an extension of the Krylov-Safonov local theory. Instead, the work <cit.> relies on the Ishii-Lions method (where an auxiliary function with doubled variables is used), leads to a Lipschitz bound, and requires Hölder continuous coefficients, but allows for weaker ellipticity assumptions. Thus, we are presenting here a third approach that applies to some new equations but, at the same time, not to all of theirs.

First derivative estimates do not require the convexity of the fully nonlinear operator, as shown in <cit.>. However, to simplify the exposition and for the reasons explained next in relation with a forthcoming viscosity theory, in this paper we restrict ourselves to convex equations.

The results in both <cit.> and <cit.> apply to viscosity solutions. Instead, to keep the technical points as simple as possible, in this first work we will restrict ourselves to smooth solutions. Indeed, the extension to the viscosity framework will require the use of incremental quotients within the auxiliary functions, as well as introducing a new averaged zeroth order term. We have not found this setting implemented in the literature, even in the local case,[Recall that in the monograph <cit.>, for instance, the Bernstein technique is carried out only for smooth solutions.] and thus it will be the content of a forthcoming article <cit.> on the viscosity framework.
As a byproduct of it, the use of incremental quotients within the auxiliary functions will allow us to include more general Bellman and Isaacs-type fractional equations —as in the optimal setting of <cit.>. We discuss these issues in Section <ref> below, while
Section <ref> presents known results on the existence and regularity of solutions to our equations.

In the next subsections, we describe in detail the
framework and results of our work.

§.§ Pucci-type equations in the presence of extensions

We start dealing with the case of Pucci-type equations associated to affine transformations of the fractional Laplacian with elliptic matrices.
To built them, given constants $0<\lambda\leq \Lambda$ we let

\begin{equation}\label{ELLMA}
{\mathcal{A}}={\mathcal{A}}_{\lambda,\Lambda}\text{ be the set of $n$-dimensional symmetric
matrices with eigenvalues in }[\lambda,\Lambda].
\end{equation}

Now, for a given $s\in(0,1)$, we define the operator

\begin{equation} \label{TUTTIS}\begin{split}
{\mathcal{L}}_A u(x)\,&:= c_{n,s} \, {\rm P.V.}\,
\int_{\R^n}\frac{u(x)-u(y)}{|A(x-y)|^{n+2s}}\, dy\\&
:=c_{n,s} \, \lim_{\e\searrow0}\,
\int_{\R^n\setminus B_\e(x)}\frac{u(x)-u(y)}{|A(x-y)|^{n+2s}}\, dy,
\end{split}\end{equation}

where $c_{n,s}>0$ is a suitable normalizing
constant which makes that, when $A={\rm Id}$ is the identity, ${\mathcal{L}}_{\rm Id}$ becomes a fraction of the classical Laplacian, that is,

{\mathcal{L}}_{\rm Id} = (-\Delta)^{s}.

The above limit is well defined, and finite, whenever $u$ is a $C^{1,1}$ function (locally) which is bounded in all of $\R^n$.
To ease the notation, the principal value ${\rm P.V.}$ will be omitted from now on.

We can now consider the maximal operator

\begin{equation}\label{TUTTIS2}
{\mathcal{M}}_{{\mathcal{A}}} u(x):=
\sup_{A\in {\mathcal{A}}}\Big( {\mathcal{L}}_A u(x)-g_A(x)\Big),\end{equation}

where $g_A$
is a given Lipschitz continuous function in a ball $B_R\subset\R^n$,
for every $A\in {{\mathcal{A}}}$.
We will also assume continuity of $g_A$
with respect to the set of parameters $A\in{{\mathcal{A}}}$:

\begin{equation}\label{Gcon}
{\mbox{if~$A_k\to A$ as~$k\to+\infty$, then }}\lim_{k\to+\infty}
g_{A_k}(x)= g_A(x) {\mbox{ for all }}x\in B_R.
\end{equation}

Some existence and regularity results for (<ref>) will be mentioned in Section <ref>.

By developing a Bernstein technique in this framework, we
establish first and one-sided second derivative bounds
for solutions of (<ref>).
Our estimates are uniform as the order of the operators converges to two. In this respect,
note that the operators ${\mathcal{L}}_A$
recover, in the limit $s\nearrow1$, every
second order linear elliptic operator in nondivergence form
(see Remark <ref> below or Section 6 in <cit.>).
Hence, the Bernstein method that we develop here
provides a unified theory up to $s=1$, with uniform constants in the bounds.

To prove the estimates, we first need to extend the Bernstein technique to the simplest linear operator: the fractional Laplacian $(-\Delta)^s$ with $0<s<1$, as defined above. In the case of first derivative estimates, the computations (<ref>)-(<ref>) for the classical Laplacian will easily carry over the extension operator for the fractional Laplacian, leading to the following analogue of (<ref>). Note that the result is uniform as $s$ tends to 1. To guarantee that the fractional Laplacian is well defined when acting on a smooth function $u$ and also on the auxiliary functions built from its first derivatives, we will assume that both $u$ and $\nabla u$ are bounded in all of $\R^n$, that is, $u\in W^{1,\infty}(\R^n)$.

Let $s\in (0,1)$, $u \in C^\infty (\R^n)\cap W^{1,\infty} (\R^n)$,
$\eta \in C^\infty (\R^n)\cap W^{2,\infty} (\R^n)$,
$\sigma>0$, and $e\in\R^{n}$ with $|e|=1$.
Then, we have

\begin{equation}\label{VE:FRA-new}\begin{split}
\eta^2(\partial_e u)^2+\sigma u^2\right)\;\le\;
2\eta^2\partial_e u\,(- \Delta)^s \partial_eu
+2 \sigma u\, (- \Delta)^s u\qquad \text{ if } \sigma\ge \sigma_0,
\end{split}\end{equation}

everywhere  in all of $\R^n$, for some constant $\sigma_0$ depending only on $n$ and $\Vert\eta\Vert_{C^2(\R^n)}$ —and, in particular, independent of $s$.

This proposition will be proved in Section <ref>,
where other related results (involving second derivatives of $u$)
will also be presented.

In relation with (<ref>), we will later comment on the
related inequality that we prove in order to obtain one-sided second
derivative bounds; see (<ref>) and
Proposition <ref> (as well as Proposition <ref>). Instead, for general kernels, the corresponding inequality is unkown; see Open problem <ref> .

Given a Lipschitz function $f$, we consider the equation

\begin{equation}\label{MUCCI:EQ}
{\mathcal{M}}_{{\mathcal{A}}} u(x):=
\sup_{A\in {\mathcal{A}}}\Big( {\mathcal{L}}_A u(x)-g_A(x)\Big)=
f(x)\quad {\mbox{for all }}x\in B_R.
\end{equation}

To state our precise estimates, we define

\begin{equation}
\label{991199}
G_0:=\sup_{{x\in B_R}\atop{A\in{\mathcal{A}}}} | (f+g_{A})_-(x)|, \qquad
G_1:=\sup_{{x\in B_R}\atop{A\in{\mathcal{A}}}} |\nabla (f+g_{A})(x)|,
\end{equation}

\begin{equation}
\label{9911994}
G_{e,1}:=\sup_{{x\in B_R}\atop{A\in{\mathcal{A}}}} |\partial_e (f+g_{A})(x)|, \quad\text{ and }\quad
G_{e,2}:=\sup_{{x\in B_R}\atop{A\in{\mathcal{A}}}} \big|\big(\partial_e^2 (
\end{equation}

Here and throughout the paper, we use the notation
$r_+:=\max\{r,0\}$ and $r_-:=\max\{-r,0\}$ to denote the positive
and negative parts, respectively. Notice that $G_{e,2}$ is finite only when $f+g_A$
is a semiconcave function in $B_R$ for all $A\in{\mathcal{A}}$. We will need to make this assumption for the validity of the second-derivative estimates —but not for the first derivative bounds. We recall that $f+g_A$ is semiconcave in $B_R$ when $(f+g_A)(x)-C_A|x|^2$ is a concave function
of $x\in B_R$ for some constant $C_A$. In such case, $f+g_A$ will be pointwise
twice differentiable at almost every point in $B_R$, and thus the quantity $G_{e,2}$ defined above
makes sense.

The following is our result. Here we need the smooth function $u$ to belong to $W^{2,\infty}(\R^n)$, since, within the proofs, the fractional operators will act on
derivatives of $u$ up to order two; in this way, second derivatives will be smooth functions bounded in all space.

Given $s\in (0,1)$, $0<\lambda\leq\Lambda$, and
Lipschitz functions $f$ and $g_{A}$ in $B_R$, assume (<ref>) and let $u \in C^\infty (\R^n)\cap W^{2,\infty} (\R^n)$ be a solution of (<ref>).

\begin{equation}\label{9ixsjXVaBAM-1}
\sup_{B_{R/2}} |\nabla u|
\le \frac{C}{R}\,\Big( \| u\|_{L^\infty(\R^n)}+
\,\| u\|_{L^\infty(B_R)}^{1/2}+R^{1+2s}\,G_1
\Big)
\end{equation}

for some constant $C$ depending only
on $n$, $\lambda$, and $\Lambda$.

If $f+g_A$ is, in addition,  semiconcave in $B_R$
for all $A\in{\mathcal{A}}$, then we have

\begin{equation}\label{9ixsjXVaBAM-2}\begin{split}&
\sup_{B_{R/2}} \partial^2_e u\le \frac{C}{R^2}\,
\Big(
\| u\|_{L^\infty(\R^n)}+R^{s}\,G_0^{1/2}
\,\| u\|_{L^\infty(B_R)}^{1/2}+R^{1+2s}\,G_{e,1}
\Big)
\end{split}\end{equation}

for every $e\in\R^n$ with $|e|=1$,
where $C$ is as before.

We stress that the one-sided second derivative bound (<ref>)
is somehow surprising, since the order $2s$ of the operator
is smaller than $2$. Hence, full bounds on the second
derivatives are not be expected to hold. The best
regularity theory available for this equation arrives at $C^{1+\varepsilon}\cup C^{2s+\varepsilon}$,
with $\max\{1+\varepsilon,2s+\varepsilon\}<2$; see <cit.> and <cit.>, respectively, and Section <ref> below.

Since solutions of equation (<ref>) will not be, in general, smooth (see Section <ref>), the full power of our techniques will be accomplished by their extension to the viscosity framework, carried out in our forthcoming article <cit.>. In any case, notice that the current setting for the gradient estimate of Theorem <ref> applies to a large number of equations. Indeed, given any smooth function $u\in W^{2,\infty} (\R^n)$ we may define $f$ by (<ref>), a function that, taking into account the presence of the supremum in this definition, will satisfy our Lipschitz assumption on the right-hand side. Instead, the one-sided second derivative estimate may not apply, since the function $f$ defined in this way may fail to be semiconcave.

Once estimate (<ref>) is extended to viscosity solutions $u$, it will establish that they are semiconcave functions in $B_{R/2}$. From this
and the  Alexandrov-Buselman-Feller theorem, we will have that every viscosity solution $u$ is pointwise twice differentiable at almost every point (in the sense
of Chapter 1 in <cit.>). That is, it
can be approximated by a polynomial of degree two at almost every point.

Recall that our estimates are uniform as $s\nearrow 1$ and thus provide a uniform theory up to the second order case.
The proof of the first derivative estimate in the theorem
will follow from inequality (<ref>), by
choosing an appropriate cutoff function $\eta$ and after a change of variables to replace $(-\Delta)^s$ by the operators ${\mathcal{L}}_A$. As mentioned, (<ref>) will be easily proved using the extension property of the fractional Laplacian.

Establishing the one-sided second derivative bounds will be similar, but involves new ingredients that will be especially relevant when dealing with operators having general kernels, discussed in next subsection. Indeed, in (<ref>) we first need to replace $u$ by $v=\partial_e u$, but since we only expect a one-sided second derivative bound from above, we must consider instead the auxiliary function

\begin{equation}\label{eqn-v-pos}
\eta^2(\partial_{e} v)_{+}^2+\sigma v^2
\end{equation}

involving a positive part,[In the local case there is no need to consider positive parts; see <cit.>. It is enough to apply the maximum principle in a ball intersected with
the set where $\partial_{e} v=\partial_{e}^{2} u$ is positive, and then check that
the auxiliary function is controlled on the boundary of such set. This approach does not work in the nonlocal framework due to the influence of the exterior datum.] where $v=\partial_{e}u$.
Also, since the maximum principle involves the exterior data and $v=\partial_{e}u$ will not be controlled (by our first derivative estimates) in all of $\R^{n}$ but only where the equation is satisfied
(unless our solution is global), we will need
to consider an auxiliary function involving $u$, $\partial_{e}u$, and $(\partial_{e}^{2} u)_{+}$ all together, as well as two cutoffs;
see Theorem <ref>
and the comments preceding it. In
the next subsection we will comment on the analogue of (<ref>) for the auxiliary function (<ref>);
see Open problem <ref>.

A similar result to Theorem <ref> but dealing with linear and convex fully nonlinear operators with no definite order will be presented in Subsection <ref>.

§.§ Pucci-type equations for general integro-differential operators

In this paper we also take into account the case of operators with more
general kernels, that is, kernels which are not pure powers, neither
rotationally invariant. Here no local extension problem will be
available and, for this reason, the computations in the Bernstein
technique become purely nonlocal and very delicate. Still, our
main contributions, Theorems <ref> and <ref> below, establish a key inequality for the auxiliary
function that includes a new error term but still leads to first
derivative bounds.

Let $K:\R^n\to (0,+\infty]$ satisfy

\begin{equation}\label{EVEN}
K(z)=K(-z)\qquad{\mbox{ for all }} z\in\R^n\setminus\{0\}
\end{equation}

and, for some $s\in(0,1)$,

\begin{equation}\label{KLIM}
\frac{\CDUEDUE\,s(1-s)}{|z|^{n+2s}}\le
K(z)\le \frac{\CDUE \,s(1-s)}{|z|^{n+2s}}\qquad{\mbox{ for all }}z\in\R^n\setminus\{0\},\end{equation}

where $0<\CDUEDUE\leq\CDUE$ are given constants. In our main results we will also assume $K$ to be $C^{2}$ in $\R^n\setminus\{0\}$ and to satisfy

\begin{equation}\label{KC1}
|z|\,| \nabla K(z)|
+ |z|^2\,|D^2 K(z)|
\le \CUNO \, K(z)\qquad{\mbox{ for all }}\;z\in\R^n\setminus\{0\},\end{equation}

for some constant $\CUNO>0$. This is the class ${\mathcal{L}}_{2}$ of kernels introduced in <cit.>.
We consider the linear operator

\begin{equation} \label{OP:L:BIS}
{\mathcal{L}}_K u(x)\,:=\int_{\R^n}\big( u(x)-u(y)\big)\,K(x-y)\,dy
\end{equation}

defined, as before, in the principal value sense. The operator is well defined on $C^{1,1}$ functions which are bounded in all of $\R^n$. When $K(z)=c_{n,s} |z|^{-n-2s}$, it is the fractional Laplacian.

The assumptions in (<ref>)-(<ref>)
are satisfied by the class of general stable symmetric
operators, where the kernels are
given by

\begin{equation}\label{SETA} K(z) :=\frac1{2 \,|z|^{n+2s}}\left(
a\left(-\frac{z}{|z|}\right) \right),\end{equation}

under appropriate hypotheses on the positive function $a$. See, for instance, (1.3) in <cit.>.

Integro-differential operators with such kernels naturally arise in the Lévy-Khintchine probabilistic formula, to take into
account Poisson processes with jumps; see e.g. Section 2.2 in <cit.>.
They possess applications in several fields; see e.g. Sections 1.2 and 1.3
in <cit.>.
In spite of many similarities with the case of the fractional Laplacian, they also
present some important differences, in terms of regularity results, with respect to
the fractional Laplacian.[See for instance the counterexample
to regularity in <cit.>, which shows that general integral kernels
with singular spectral measures
may propagate the boundary behavior thus producing interior singularities,
in contrast with the case of the fractional Laplacian.
More precisely, in <cit.>, one considers operators of the form

$$ Lu(x):=\int_0^{+\infty}\left(\int_{\partial B_1}\frac{2u(x)-u(x+\rho\omega)-u(x-\rho\omega)}{\rho^{1+2s}}
d\widetilde a(\omega)
\right)\,d\rho,$$

where $\widetilde a$ is a nonnegative measure on $\partial B_1$
(called in jargon the “spectral measure”) satisfying suitable ellipticity assumptions,
see (1.3) in <cit.>. In the specific case of absolutely continuous measures
of the form $ d\widetilde a(\omega)=(a(\omega)+a(-\omega))/4$,
one recovers the kernel $K$ in (<ref>).]

For the operators (<ref>), with general kernels, no extension technique is available. Therefore, the
Bernstein method must rely solely on integral computations made
“downstairs”, that is, in $\R^{n}$. This turns out to be a very
delicate issue. In fact, the validity of the key inequality
(<ref>) for the fractional Laplacian remains unknown in the
case of the operator ${\mathcal{L}}_{K}$ (see Open problem <ref>
below). Our main contribution is to prove the inequality with an error
term $E$ which will be absorbable (by scaling properties) at the end of
the proof of first derivative estimates.

The analogue for ${\mathcal{L}}_{K}$ of the pivotal inequality (<ref>), which will include an error term (generically named $E$ in the following statement), will be proved using the next proposition.

Let $K$ satisfy (<ref>)
and (<ref>), and let ${\mathcal{L}}_K$ be defined by (<ref>). Given a function $u\in C^\infty(\R^{n})\cap W^{1,\infty}(\R^n)$,
$\eta\in C^\infty(\R^n)\cap L^{\infty}(\R^n)$,
$e\in\R^{n}$ with $|e|=1$, $\sigma>0$, and $\CALERRE\in\R$, consider

\begin{equation}\label{COM0129pq} \varphi:=\eta^2 (\partial_e u)^2 +\sigma u^2.\end{equation}

Then, the inequality

\begin{equation}\label{MAY18-1}
{\mathcal{L}}_K \varphi\le 2\eta^2\,\partial_e u\,{\mathcal{L}}_K \partial_e u+2\sigma u\,{\mathcal{L}}_K u+\CALERRE
\end{equation}

holds at a point $x\in\R^n$ if and only if

\begin{equation}\label{MAY18-2}
\begin{split}&
2\int_{\R^n}\eta(x)\,\big( \eta(x)-\eta(y)\big)\,\partial_e u(x)\,\partial_e u(y)\,K(x-y)\,dy
\\ &\qquad\le\,\int_{\R^n}\big|\eta(x)\,\partial_e u(x)-\eta(y)\,\partial_e u(y)\big|^2\,K(x-y)\,dy
\\&\qquad\qquad+\sigma\int_{\R^n}\big|u(x)- u(y)\big|^2\,K(x-y)\,dy

It is simple to check that all integrals in (<ref>) are well defined in the principal value sense; see the comments following Proposition <ref>. A useful aspect of (<ref>) is that all terms in its right-hand side are nonnegative, which is not necessarily the case in (<ref>).

Using the criterium stated in the previous proposition, we will establish the following key inequality for the operator ${\mathcal{L}}_K$. It will lead to first derivative estimates. Note that the inequality differs from a possible optimal one by a “small error”; see the open problems below.
Its proof will contain several quite surprising weighted integral cancellations.

Let $K$
satisfy (<ref>),
(<ref>), and (<ref>), and let
${\mathcal{L}}_K$ be defined by (<ref>). Let $u\in C^\infty(\R^n)\cap W^{1,\infty}(\R^n)$ and $\eta\in C^\infty(\R^n)\cap W^{2,\infty}(\R^n)$.

Then, for every $\e>0$ there exists a constant $\sigma_\e>0$
depending only on $\e$, $\|\eta\|_{C^2(\R^n)}$, and
the structural constants $n$, $s$, $\CDUEDUE$,
$\CDUE$, and $\CUNO$ in (<ref>) and (<ref>), such that

\begin{equation}\label{T6}
{\mathcal{L}}_K \big( \eta^2(\partial_e u)^2+\sigma_\e u^2\big)\le
2\eta^2\partial_e u \, {\mathcal{L}}_K\partial_e u+
2\sigma_\e \,u{\mathcal{L}}_Ku+\e^2\,\|\partial_e u\|^2_{L^\infty(B_3)}
\quad{\mbox{everywhere in }}B_2,
\end{equation}

for every $e\in\R^n$ with $|e|=1$.

Note that the error depends on the $L^{\infty}$-norm of the first derivative in a larger ball than the ball where the inequality is claimed.

Recall that our proof of (<ref>)
heavily relies on the extension method.
It would be very interesting to prove inequality (<ref>)
without using the extension. We only know how to do this when the function $u$ is assumed to be $s$-harmonic; see Lemma <ref> and its proof.[Note that our proof of Lemma <ref> (which does not use the extension) is uniform as $s$ tends to 1. This is also the case for the extension proof of (<ref>). Instead, Theorem <ref> (and as a consequence, Theorem <ref> below) are not. In this respect, it would be very interesting to find a proof of Theorem <ref> which is uniform as $s\nearrow 1$.]
In view that we know the criterium (<ref>) to be true with $E=0$ when $K(z)=|z|^{-n-2s}$, we still find intriguing not to be able to prove it directly in $\R^{n}$ (for this kernel and with $E=0$) without using the extension.
Finding such a proof could shed light into the following question.

Given a kernel $K$ as above, does Theorem <ref>
hold true without the additional small remainder? That is,
does (<ref>) hold true with $\e=0$?

Theorem <ref>
will serve as the cornerstone to prove
first derivative estimates for maximal type operators
involving general fractional kernels.
To state the result, we consider a compact set of indexes ${\mathcal{B}}$, as well as
kernels $K_B$ and Lipschitz functions $g_B:B_R\subset\R^n\to\R$, for every $B\in
{\mathcal{B}}$, satisfying

\begin{equation}\label{Gcon-BIS}{\mbox{if~$B_j\to B$ as~$j\to+\infty$, then }}\lim_{j\to+\infty}g_{B_j}(x)= g_B(x) {\mbox{ for all }}x\in B_R.\end{equation}

The following is one of our main results.

Let ${\mathcal{B}}$ be a compact set and $\{ K_B\}_{B\in{\mathcal{B}}}$ be kernels satisfying (<ref>),
(<ref>), and (<ref>) $($all with the same structural constants $s$, $\CDUEDUE$,
$\CDUE$, and $\CUNO)$. For $B\in{\mathcal{B}}$, let
$g _B$ be a Lipschitz function in $B_R$
and assume that (<ref>)
is satisfied. Let also $f$ be Lipschitz in $B_R$ and set

\begin{eqnarray*}
\| ( f+g_B)_-\|_{L^\infty(B_R)} \quad\text{ and }\quad
\|\nabla (f+g_B)\|_{L^\infty(B_R)}.
\end{eqnarray*}

Let $u\in C^\infty(\R^{n})\cap W^{1,\infty}(\R^n)$ be a
solution of

\begin{equation}\label{1.23BIS} \sup_{B\in{\mathcal{B}}}\Big(
= f(x)\quad {\mbox{for all }}x\in B_R.\end{equation}

\begin{equation*}
\sup_{B_{R/2}} |\nabla u|\le
\frac{C}{R}\,\Big( \| u\|_{L^\infty(\R^n)}+R^{s}
H_0^{1/2}\| u\|_{L^\infty(B_R)}^{1/2}
+R^{1+2s} H_1
\Big)
\end{equation*}

for some constant $C$ depending only on
the structural constants $n$, $s$, $\CDUEDUE$,
$\CDUE$, and $\CUNO$ in (<ref>)
and (<ref>).

Recall that this result applies to a large number of equations, even if it assumes $u$ to be smooth,  by the discussion following Theorem <ref>.

We do not know whether one-sided
second derivative bounds hold true for general kernels, in the setting of the previous result.
Recall that they do hold, by Theorem <ref>, for the Bellman operator built from affine transformations of the fractional Laplacian.

To establish such a result, one would need to prove an inequality similar to (<ref>), but with $u$ and $\partial_{e} u$ replaced by $v$ and $(\partial_{e} v)_{+}$, as explained in (<ref>) (here the arbitrary function $v$ plays the role of $\partial_{e} u$). Recall that the positive part comes from the fact that we only expect one-sided estimates for the second derivatives. For the fractional Laplacian we know that such inequality holds, and without an error term, by the results of Section <ref>. For general kernels, the analogue of the key inequality (<ref>) (but now for the auxiliary function involving the positive part) is stated in Proposition <ref>. We do not know, however, how to prove it, even with error terms.

§.§ Other fully nonlinear equations
in the presence of extensions

This subsection presents our results
for fully nonlinear equations with no definite order. They are built as superposition of linear operators
which can be studied through local extensions.

Given a probability measure $\mu$ on $[0,1]$, i.e.,

\begin{equation}\label{1.0} \mu\ge0\qquad{\mbox{and}}\qquad\mu([0,1])=\mu(\R)=1,\end{equation}

we define

\begin{equation} \label{OP:L}
{\mathcal{L}}_\mu u(x):=\int_0^1 (-\Delta)^s u(x)\,d\mu(s).\end{equation}

In case of $\mu$ being a Dirac's delta
at some $s\in[0,1]$,
$ {\mathcal{L}}_\mu$ reduces to the fractional Laplacian $(-\Delta)^s$. For $s=0$, $(-\Delta)^{0}$ denotes the Identity, $(-\Delta)^{0}u=u$.
The interest of including $s=0$ is to allow a unified treatment
of fully nonlinear equations and obstacle problems; see Corollary <ref> below.

The operators ${\mathcal{L}}_\mu$ have been studied in <cit.>, in relation with Allen-Cahn type equations, through local extension methods.

Let $J\in\N$, $\mu_1,\dots,\mu_J$ be probability measures on $[0,1]$,
and $F\in C(\R^J)$.
We assume that there exist functions[Even though the functions $
\alpha_j$ are not required to be continuous, they are assumed to be defined
everywhere, and not almost everywhere. This is consistent with the general setting of
the article, in which the equations are supposed to be satisfied
everywhere in a given domain —indeed, our terminology “for all” has to be taken literally, and
not in the meaning of “for almost all”. This will important when proving the maximum principle, since it will require to evaluate the equation at a maximum point. This framework coincides, for instance,
with the one of Chapter 3 in the monograph <cit.>.]
$\alpha_1,\dots,\alpha_J$ defined in $\R^{J}$ and a constant $\vartheta_0>0$
such that

\begin{equation}\label{CONV:POSITIVA}
\alpha_j(p)\ge0\qquad{\mbox{for every $p\in
\R^J$ and $j\in\{1,\dots,J\}$,}}\end{equation}

\begin{equation}\label{MON:F}
\sum_{j=1}^J\alpha_j(p)\ge\vartheta_0\qquad{\mbox{for every $p\in
\R^J$}},
\end{equation}

\begin{equation}\label{CONV:F}
\qquad{\mbox{for every~$q$ and $p$ in $\R^J$.}}
\end{equation}

Now, given $R>0$, let $f$ and $g_1,\dots,g_J$ be continuous functions in $B_R$.
We consider solutions $u:\R^n\to\R$ of

\begin{equation}\label{EQ}
F\big( {\mathcal{L}}_{\mu_1}u(x)- g _1(x), \;\cdots,\; {\mathcal{L}}_{\mu_J}u(x)- g _J(x) \big)=f(x)\qquad{\mbox{for all }}x\in B_R.
\end{equation}

By (<ref>), $F$ is convex. Note that the hypotheses on $F$ represent, all three together, convexity and a quantification of ellipticity.
They are satisfied, for instance, by Bellman-type operators, for which $F(p)=\max\{p_{1},\ldots,p_{J}\}$.
Indeed, it suffices to define $\alpha_{j}(p)=1$ if $j$ is the smallest index for which $p_{j}=F(p)$ and $\alpha_{j}(p)=0$ for all other indexes. Since the operators $ {\mathcal{L}}_{\mu_j}$ are of indefinite order, our class of Bellman-type equations includes the obstacle problem for the fractional Laplacian, which we treat in Corollary <ref> below.

The three assumptions on $F$ are also satisfied if $F$ is $C^1$, convex,
in each of its coordinate variables, and satisfies
$ \Sigma_{j=1}^J\partial_{p_j}F(p)\ge\vartheta_0$ for every $p\in \R^J$ —here we take
$\alpha_j(p):=\partial_{p_j}F (p)$.

Note that we are dealing with a very general class of equations of indefinite order, which includes the model equation

\begin{equation}\label{FORM83} F\big( (-\Delta)^{s_1} u,
\;\dots,\; (-\Delta)^{s_J} u \big)=f(x),\end{equation}

with $s_j\in[0,1]$.

A more general framework consists of making affine changes of variables for each index $j\in J$,
see[In this way one includes here those equations of Subsection <ref> built from a finite number $J$ of linear operators. In fact, since our estimates will be independent of $J$, it is possible to
deduce Theorem <ref> in that subsection from the results of the current one,
by an approximation and limiting argument.]
Remark <ref>.
In such generality, if we assumed $F$ to be uniformly elliptic but not necessarily convex, equation (<ref>) would include the classical inf-sup Isaacs equations built from a finite number of second order linear operators; see <cit.>.
Therefore, since we will establish one-sided second derivative bounds, the convexity assumption on $F$ in the following theorem cannot be dropped, in view of the classical counterexamples to $C^{1,1}$ regularity by
Nadirashvili and Vlăduţ <cit.>
for nonconvex fully nonlinear equations of second order.[Two comments are in order here. First, once a one-sided second derivative bound for a second order fully nonlinear uniformly elliptic equation is established, it automatically leads to full second derivative estimates (by using the equation itself; see the
Bernstein technique described in Chapter 9 of <cit.>,
and in particular inequality (9.5) combined with Lemma 6.4
in <cit.>). Second, recall that Isaacs equations cover all possible fully
nonlinear elliptic equations of second order; see Remark 1.5 in <cit.>.]

Among the operators ${\mathcal{L}}_\mu$ in (<ref>)–(<ref>), a special class
is provided by the ones for
which $\mu(\{0\})=0$, since they involve only nontrivial
fractional (and perhaps also classical) derivatives
(namely, they do not involve the identity operator).
In our following result, which deals with the operators ${\mathcal{L}}_{\mu_1},\dots,
{\mathcal{L}}_{\mu_J}$ in (<ref>), the notation

\[ \omega_0:=\begin{dcases}
0 & {\mbox{ if~$\mu_j(\{0\})=0$ for all~$j\in\{1,\dots,J\}$,}}\\
1 & {\mbox{otherwise,}}\end{dcases}\]

will be useful.
In addition, given functions $f$ and $g_{j}$ as in its statement, we define the quantities

\begin{equation*}
\gamma_1 :=
\|\nabla f\|_{L^\infty(B_R)}+\sup_{p\in\R^J}
\sum_{j=1}^J \alpha_j(p)\,\| \nabla g_j\|_{L^\infty(B_R)},
\end{equation*}

\begin{equation*}
\gamma_2 :=\|(D^2f)_+\|_{L^\infty(B_R)}+\sup_{p\in\R^J}
\sum_{j=1}^J \alpha_j(p)\,\| (D^2 g_j)_+ \|_{L^\infty(B_R)}.
\end{equation*}

Our result establishes first and one-sided second derivative bounds for solutions of (<ref>). The estimates are uniform in the number $J$ of operators.

Given $\mu$ and $F$ satisfying (<ref>)-(<ref>), let
$f$ and $g _j$ be Lipschitz functions
in $B_R$
for $j=1,\ldots,J$, and $u \in C^{\infty}(\R^{n})\cap W^{2,\infty}(\R^{n})$ be a
solution of (<ref>) in $B_R$. Set

\begin{equation}\label{EFFEDE} g(x):=-
F\big(-g _1(x),\dots,-g _J(x)\big)\qquad
{\mbox{ for all }}x\in B_R.\end{equation}

\begin{equation}\label{EST-MAIN-1}\begin{split}
\sup_{B_{R/2}} |\nabla u|\le\,& C\,\Big(
( \omega_0+R^{-1})\,\| u\|_{L^\infty(\R^n)}\\&\quad
\|(f+g)_-\|_{L^\infty(B_R)}^{1/2}\|u\|_{L^\infty(B_R)}^{1/2}+

for some constant $C$ depending only
on $n$ and $\vartheta_0$.

If $f$ and $g _j$ are, in addition, semiconcave
in $B_R$
for $j=1,\ldots,J$, then we have

\begin{equation}\label{EST-MAIN-2}\begin{split}
\sup_{B_{R/2}} \partial^2_e u\le\,& C\, \Big( ( \omega_0+R^{-1}) R^{-1}\,\| u\|_{L^\infty(\R^n)}
\\&\qquad+(1+R^2)\,(R^{-1}\gamma_1+\gamma_2)\Big)\end{split}\end{equation}

for every $e\in\R^n$ with $|e|=1$, where $C$ is as before.

It is interesting to point out that, when all functions $g_{j}\equiv 0$, we do not need to assume the nonlinearity $F$ to be globally Lipschitz, since this is not required by our hypotheses (<ref>)-(<ref>).
Instead, when
the functions $g_j$ are not identically zero, the quantities $\gamma_{1}$ and $\gamma_{2}$ will be finite if one assumes $F$ to be globally Lipschitz.

We also remark that the estimates of
Theorem <ref> are new, to the best of our
knowledge, even in the case when $F$ is linear,
even for ${\mathcal{L}}_{\mu_j}=(-\Delta)^{s_j}$,
and even when all the functions $g_{j}$ are taken to be zero.

The unusual dependence on the radius $R$
of the right-hand side
of (<ref>) and (<ref>) is due
to the fact that our operator is not scale invariant,
since it is of indefinite order. In this respect, the term $1+R^2$ takes into account
all possible scaling exponents induced by the operator,
in which the fractional order ranges from $0$ (corresponding to $R^0=1$)
to $2$ (corresponding to $R^2$).

As in Subsection <ref>, the one-sided second derivative estimate (<ref>) is somehow surprising since, for operators which could be of order smaller than two,
second derivative estimates are not expected to hold.

Theorem <ref> can be easily extended by changing
the definition (<ref>) of ${\mathcal{L}}_\mu$
to involve more general operators of the form

\begin{equation}\label{STxWi} \int_0^1
\left( -\sum_{i,j=1}^n M_{ij}(s)\, \partial_{x_i x_j}^2 u\right)^s\,d\mu(s),
\end{equation}

containing fractions of second order elliptic operators
with constant coefficients instead of fractions of the Laplacian. Here, for every $s\in[0,1]$ we are given a matrix $M(s)\in {\mathcal{A}}_{\Lambda^{-2},\lambda^{-2}}$, where the classes $ {\mathcal{A}}$ were defined in (<ref>).
The proofs would remain almost unchanged by performing, for every $s\in [0,1]$,
the change of variables $\overline x=Ax$ (and then consider the function $\overline u (\overline x)= u(A^{-1}\overline x)$), where $A=A(s)\in
{\mathcal{A}}_{\lambda,\Lambda}$ is such that $A^{-2}=M:=M(s )$ —as done later in (<ref>) within the proof of Theorem <ref>. In this way, one can see that

\left( -\sum_{i,j=1}^n M_{ij}(s)\, \partial_{x_i x_j}^2 u\right)^s (x) =
c_{n,s}\, \det A(s) \int_{\R^n}\frac{u(x)-u(y)}{|A(s)\,(x-y)|^{n+2s}} dy,

which coincides with the operator ${\mathcal{L}}_{A}$ previously considered in (<ref>), up to a multiplicative constant. These fractional operators thus possess the same type of
extension properties as the fractional Laplacian.

The previous equality can be alternatively checked through Fourier symbols since

\begin{eqnarray*}&&
\left(\sum_{i,j=1}^n M_{ij}\xi_i\xi_j \right)^s=
\left(\sum_{i,j,k=1}^n (A^{-1})_{ik}(A^{-1})_{kj}\, \xi_i\xi_j \right)^s=|\overline{\xi}|^{2s}
=\frac{c_{n,s}}{2}\int_{\R^n} \frac{2-e^{i\overline{z}\cdot\overline{\xi}}-e^{-i\overline{z}\cdot\overline{\xi}}}{|\overline{z}|^{n+2s}}\,d\overline{z}\\
=\frac{c_{n,s}}{2}\int_{\R^n} \frac{2-e^{i(A^{-1}\overline{z})\cdot\xi}-e^{-i(A^{-1}\overline{z})\cdot\xi}}{|\overline{z}|^{n+2s}}\,d\overline{z}=
\frac{c_{n,s}\,\det A}{2}\,
\int_{\R^n} \frac{2-e^{iz\cdot\xi}-e^{-iz\cdot\xi}}{|Az|^{n+2s}}\,dz,
\end{eqnarray*}

where $\overline\xi:=A^{-1}\xi$.

Theorem <ref> includes, as a particular case, the obstacle problem for the fractional
also called “thin obstacle problem” or “Signorini problem”. For this, we take the measures to be $\mu_{0}=\delta_{s}$ for some $s\in [0,1]$ and $\mu_{1}=\delta_{0}$, and $F$ to be the $\max$ operator. Within the estimates of Theorem <ref>,
in this particular case we will be able to control the $L^\infty$-norm of $(f+g)_-$ by that of $u$.

Let $s\in[0,1]$, $f$ and $g$ be Lipschitz functions in $B_1$,
and $u\in C^{\infty}(\R^{n})\cap W^{2,\infty}(\R^{n})$ be a
solution of

\begin{equation}\label{934rtiI234}
\max\Big\{
(-\Delta)^{s}u, \;u- g
\Big\}=f \quad {\mbox{everywhere in }} B_1.
\end{equation}

\begin{equation}\label{KSM-P8SKNDN}
\sup_{B_{1/2}} |\nabla u|\le C\,\Big(\| u\|_{L^\infty(\R^n)}
+\| \nabla f\|_{L^\infty(B_1)}+\| \nabla g \|_{L^\infty(B_1)}
\Big),\end{equation}

for some constant $C$ depending only
on $n$.

If $f$ and $g$ are, in addition, semiconcave in $B_1$, then we have

\begin{equation}\label{SEMI}\begin{split}&
\sup_{B_{1/2}} \partial^2_e u\le C\, \Big(\| u\|_{L^\infty(\R^n)}
+\| \partial_e f\|_{L^\infty(B_1)}
+\| (\partial_e^2 f)_+\|_{L^\infty(B_1)}\\&\qquad\qquad\qquad\qquad
+\|\partial_e g \|_{L^\infty(B_1)}
+\|(\partial_e^2 g )_+\|_{L^\infty(B_1)
\end{split}\end{equation}

for every $e\in\R^n$ with $|e|=1$,
where $C$ is as before.

In the forthcoming work <cit.>, the use of incremental quotients, as sketched in Section <ref> below, will allow us
to remove the unnecessary hypothesis $u\in C^\infty(\R^{n})\cap W^{2,\infty} (\R^{n})$ in the above result.
Still,  by the discussion following Theorem <ref>,
the gradient estimate of Corollary <ref> applies to a large number of equations.[However, it does not apply when the right-hand side $f$ is identically $0$, unless only one operator  acts, by chance, in equation (<ref>) (that is, unless $u$ does not touch the obstacle $g$ or $u\equiv g$).]

The bound (<ref>) recovers the semiconcavity
estimate for the thin obstacle problem, first proved by Athanasopoulos and Caffarelli <cit.>
and later extended by Fernández-Real <cit.> to the fully nonlinear thin obstacle problem.
In these papers the Bernstein technique was already used, but with a less flexible
auxiliary function than
in the current work: their auxiliary function is linear in
the second derivatives, while ours is quadratic. The quadratic structure has already
allowed further applications in obstacle problems. Indeed, in private communication
to the authors of <cit.> (an article that cites ours),
Fernández-Real and Jhaveri have used our method in a situation where a
polynomial solving the thin obstacle problem is subtracted to the solution and
have gotten, in this way, estimates independent of the polynomial. This required
the use of the quadratic auxiliary function, as well as the use of incremental
quotients —as we introduce in Section <ref> below
and develop in our forthcoming work <cit.>.

§.§ Other fully nonlinear equations involving
integro-differential operators

Here we present regularity results
for other fully nonlinear equations built from linear integro-differential operators
with general kernels.

As in the previous subsection, we let $J\in\N$ and $F\in C(\R^{J})$
satisfy (<ref>), (<ref>), and (<ref>)
for some functions $\alpha_j$ and a constant $\vartheta_0>0$. Given $s\in(0,1)$,
we consider kernels $K_1,\dots,K_J$
satisfying (<ref>), (<ref>), and (<ref>)
and the corresponding operators ${\mathcal{L}}_{K_1},\dots,
{\mathcal{L}}_{K_J}$ defined by (<ref>).

We deal with
solutions $u:\R^n\to\R$ of

\begin{equation}\label{EQ-onthetrace}
F\big( {\mathcal{L}}_{K_1}u-g_1, \;\cdots,\; {\mathcal{L}}_{K_J}u-g_J\big)=f
\quad{\mbox{ everywhere in }}B_R\subset\R^n,
\end{equation}

where $f$ and $g _1,\dots, g _J$
are Lipschitz functions in $B_R$.
We prove the following first derivative bound.

Let $u\in C^{\infty}(\R^{n})\cap W^{1,\infty}(\R^{n})$ be a
solution of (<ref>) everywhere
in $B_R$.

\begin{equation}\label{EST-MAIN-1-on the trace}
\sup_{B_{R/2}} |\nabla u|\le C\,R^{-1}\,
\|u\|_{L^\infty(\R^n)}
\|u\|_{L^\infty(B_{R})}^{1/2}

for some constant $C$ depending only
on the structural constants $n$, $s$, $\CDUEDUE$, $\CDUE$,
and $\CUNO$
in (<ref>), (<ref>), and (<ref>), and
on $\vartheta_0$ in (<ref>). Here above,
$g$ is defined by (<ref>) while $
\gamma_1$
is the constant defined right before Theorem <ref>.

io = [trapezium, trapezium left angle=70,
trapezium right angle=110, minimum width=2cm, minimum height=1cm,
text centered, draw=black, top color =violet!60, drop shadow]
decision = [ellipse, draw, top color=red!10, bottom color=red!50,
text width=11.5em, text badly centered, node distance=3cm, inner sep=0pt, drop shadow]
cloud = [draw, trapezium,
minimum height=2em,  text width=11.5em,
fill=yellow!60, node distance=3cm, text centered, drop shadow]
bust = [draw, tape,
minimum height=2em,  text width=13em,
fill=green!60, node distance=3cm, text centered, drop shadow]
block = [rectangle, draw, top color =cyan!60,
text width=11em, text centered, rounded corners, minimum height=4em, drop shadow]
block0 = [semicircle, draw, top color =violet!60,
text width=11em, text centered, minimum height=4em, drop shadow]
line = [draw, -]

[node distance = 2cm, auto, remember picture]
[decision] (A1) Operators
with extension;
[decision, right of=A1, node distance=7.6cm] (B1) Operators
without extension;
[cloud, below of=A1, node distance=2.5cm, top color=yellow!20,
bottom color=yellow!60] (A2) Gradient and
semiconcavity estimates;
[cloud, below of=B1, node distance=2.5cm, top color=yellow!20,
bottom color=yellow!60] (B2) Only gradient estimates;
[bust, below of=A2, node distance=2.5cm, top color=green!20,
bottom color=green!60] (AA3) Key inequalities

Propositions <ref> and <ref>;
[bust, below of=B2, node distance=2.5cm, top color=green!20,
bottom color=green!60] (B3) Key inequalities

Propositions <ref> and <ref>.

Theorem <ref>;
[block, below of=B3, node distance=3cm] (B4) Maximal operators

Theorem <ref> ;
[block, left of=B4, node distance=7.6cm] (A4) Maximal operators

Theorem <ref>;
[block, below of=A4, node distance=3cm] (A5) Convex fully nonlinear

Theorem <ref>;
[block, below of=B4, node distance=3cm] (B5) Convex fully nonlinear

Theorem <ref>;
[block, below of=A5, node distance=3cm] (A6) Obstacle problems

Corollary <ref>;
(00)[text width=4cm, above of=A1] ;
[io, right of=00, node distance=3.8cm] (11) Integro-differential operators;

[line] (11) – (A1);
[line] (11) – (B1);
[line] (A1) – (A2);
[line] (B1) – (B2);
[line] (A2) – (AA3);
[line] (AA3) – (A4);
[line] (B2) – (B3);
[line] (B3) – (B4);
[line] (A4) – (A5);
[line] (B4) – (B5);
[line] (A5) – (A6);

A road map of this article

§.§ Organization of the paper

A road map of this article is depicted in Figure <ref>.

The rest of this paper is organized as follows.

Section <ref> discusses
regularity and existence issues. In
Section <ref> we sketch an
extension of our method to the viscosity
solutions setting, whose full details will be provided
in our forthcoming paper <cit.>.

In Section <ref>
we present the necessary material on
linearized operators and the
maximum principles needed for the proofs of our main results.

Section <ref>
we state and prove a general statement
(namely, Theorem <ref>) which will be pivotal to obtain
the main results of this paper.

Section <ref> contains
the proofs of the results presented in
Subsections <ref>
and <ref>, dealing with
operators “with extensions”. More precisely,
in Subsection <ref> we
discuss Proposition <ref>
and its variants needed for the proof of the main results,
Subsection <ref> contains the proofs of
Theorem <ref>, Theorem <ref>,
and Corollary <ref>.
In the proofs of
Section <ref> we also exploit a maximum principle
in the halfspace for weighted elliptic equations, contained in
Appendix <ref>.

Section <ref> is devoted to the arguments needed
to treat the operators defined
“downstairs” in Subsections <ref>
and <ref>. Subsection <ref> contains
a key inequality for auxiliary functions, Proposition <ref>, which will complement Proposition <ref>.
In Subsection <ref>
we prove Theorem <ref>. Then,
in Subsection <ref>,
by suitable scaled estimates (which are recalled in Appendix <ref>
for the sake of completeness),
we will be able to “reabsorb” the remainder term
and complete the proof of Theorem <ref>.
Subsection <ref> contains the
proof of Theorem <ref>. Finally, in Subsection <ref>
we will give a proof of the key inequality
of Proposition <ref> without using the extension,
but with the additional assumption that $u$
is $s$-harmonic; this will rely on a modification
of the previous key inequality whose precise statement and proof
is contained in Appendix <ref>.

§ EXISTENCE AND REGULARITY OF SOLUTIONS

In order to keep our arguments as simple as possible, in this paper
the main computations are performed
assuming that the solution $u$ is smooth. Note that our auxiliary functions depend on the first or second derivatives of the solution $u$, and that we need the operator to act on them. Thus, in this article we need the solution to be at least $C^{1+2s}$ or $C^{2+2s}$, respectively. In the next section, we will sketch a forthcoming extension to embrace viscosity solutions. Here we discuss known existence and regularity results for the equations that we cover.

Concerning the existence results
for concave nonlocal fully nonlinear equations (possibly including
also the case of rough kernels), Joaquim Serra proved in
Theorem 1.3
of <cit.>
that if ${\mathcal{B}}$ is a family of indexes
and for all $B\in{\mathcal{B}}$ the kernel $K_B$
satisfies the evenness and ellipticity conditions (<ref>)
and (<ref>),
and $g_B$
and $u_0$ are Hölder continuous, with $u_0$ bounded in $\R^n\setminus B_1$, then
the problem

$$ \begin{cases}\displaystyle
\sup_{B\in{\mathcal{B}}}\Big(
{\mathcal{L}}_{K_B}u-g_B\Big)=0 & {\mbox{ in }}B_1,\\
u=u_0 & {\mbox{ in }}\R^n\setminus B_1,
\end{cases}$$

admits a unique viscosity solution which is continuous everywhere
and $C^{2s+\beta}_{\rm loc}(B_1)$ for some $\beta\in(0,1)$.
This setting applies to equations (<ref>) and (<ref>)
presented in Theorems <ref>
and <ref> of this paper.

See also <cit.>
for other existence and regularity results in related (but quite different) fractional settings.

A general approach to regularity is provided by the notion of
elliptic operators in nonlinear integro-differential equations arising from Lévy processes,
as given by Caffarelli and Silvestre in Definition 3.1 of <cit.>. When our source terms $g_A$, $g_B$, $g_1,\dots,g_J$ are constants,[The constancy of the source terms is needed to make the operator translation invariant,
as requested on page 603 of <cit.>.]
this setting
our Pucci-type equations (<ref>) and (<ref>),
as well as
the fully nonlinear framework presented here in (<ref>) (when all the operators have
the same order)
and (<ref>)
when $F:\R^J\to\R$ is differentiable
and $\partial_{p_j} F\in [C^{-1},C]$ for some $C\ge1$
and all $j\in\{1,\dots,J\}$.
Under the assumptions
that the source terms $g_A$, $g_B$, $g_1,\dots,g_J$ are constants,
that the operators have all the same order and satisfy (<ref>),
(<ref>), and a bound for $|\nabla K|$ as in (<ref>), in light of Theorem 13.1 in <cit.>
we know that the solutions in Theorems <ref>,
and <ref>
are locally $C^{1+\beta}$.
This also applies to the case dealt with
in Theorem <ref> when all
the nonlocal operators ${\mathcal{L}}_{\mu_j}$
reduce to the fractional Laplacian of some order $s$ (with the same order $s$
for all the operators). See also
Theorem 27 of <cit.> for related results
with uniform estimates.

Similarly, in the
case of vanishing source terms $f$, $g_A$, and $g_B$,
Theorem 1.1 of <cit.> leads to the local $C^{2s+\beta}$ regularity for equations (<ref>)
and (<ref>) of definite order $2s$ —for
the latter, assuming that all the kernels satisfy
assumptions (<ref>),
(<ref>), and (<ref>).
In particular, since $2s+\beta>2s$,
the corresponding equation is satisfied also in the pointwise sense
(see e.g.
the comment before Theorem 1.1 in <cit.>,
and Proposition 2.1.4 of <cit.>
for a detailed proof of this fact).

In any case, the regularity of type $C^{1+\beta}$
or $C^{2s+\beta}$ is not sufficient for the techniques discussed
in this paper, and this is the reason for which we are taking additional
regularity assumptions in our main results; see Section <ref>
for a discussion on how these additional
regularity hypotheses will be weakened in the forthcoming article <cit.>.

General results dealing with higher regularity for nonlocal
fully nonlinear equations in
domains need to address
two difficulties, namely the nonlinearity of the equation and the possible
singularity that external data may induce due to the nonlocal structure of the problem.
For linear equations satisfied in the whole of $\R^n$
the regularity theory
is well understood, and in this case solutions are $C^\infty$.
Indeed, in this situation
one can conclude via a bootstrap argument that
solutions are smooth even in case of integro-differential operators with
rough kernels,
since one can differentiate the equation without introducing errors that come from rough exterior data
(e.g., one could proceed as in Section 4 of <cit.>
without having to introduce additional cutoff functions).

To understand the effect that being a global solution has
on regularity, one can consider the toy model given by

\begin{equation}\label{TOY1}
F\big( (-\Delta)^s u(x), u(x)\big)=0 \qquad{\mbox{ for all }}x\in\R^n,
\end{equation}

with $F$ smooth and with first derivatives bounded and bounded away from zero. Assume first
that $s\in(1/2,1)$.
In this case, once $u$ belongs to $C^{2s+\beta}(B_1(x_0))$ for some $\beta\in(0,1)$ —a regularity which is known from <cit.>—
for all $x_0\in\R^n$, given a direction $e$,
one can differentiate (<ref>) and obtain

$$ \partial_{p_1} F\big( (-\Delta)^s u(x), u(x)\big)\,(-\Delta)^s u_e(x)
+\partial_{p_2} F\big( (-\Delta)^s u(x), u(x)\big)\, u_e(x)=0\qquad{\mbox{ for all }}x\in\R^n.$$

$$ \,(-\Delta)^s u_e(x)=f(x)\qquad{\mbox{ for all }}x\in\R^n,$$

$$ G(p_1,p_2):= -
\frac{ \partial_{p_2} F(p_1,p_2) }{
\partial_{p_1} F(p_1,p_2) }$$

$$ f(x):=G\big( (-\Delta)^s u(x), u(x)\big)\, u_e(x).$$

Thus, one can apply the local regularity theory for the fractional Laplacian in $B_{1/2}(x_0)$
(see e.g. Proposition 2.1.11 in <cit.>) and obtain that

\begin{equation}\label{yhgbgfhfgjuF} \| u_e \|_{C^{1+\alpha}(B_{1/2}(x_0))}\le
C\,\Big( \|u_e\|_{L^\infty(\R^n)}+\|f\|_{L^\infty(B_1(x_0))}\Big)\le

for some $\alpha\in(0,1)$,
up to renaming $C$.
But, since $s>1/2$,
for every $x_0\in\R^n$,

$$ |u_e(x_0)|\le \| u\|_{C^1(B_1(x_0))}\le \| u\|_{C^{2s+\beta}(B_1(x_0))}.$$

Hence, we have a uniform bound on $\|u\|_{C^{2+\alpha}(B_{1/2}(x_0))}$,
for all $x_0\in\R^n$, and in particular a bound on $\|D^2 u\|_{L^\infty(\R^n)}$.
Since this is a global bound, one can iterate the procedure and obtain that $u\in C^\infty(\R^n)$.
When $s\le1/2$, this method needs
to be modified, by iterating a Hölder regularity result
on the incremental quotients; see e.g.
pages 634-635 in <cit.>.
In both cases, we stress that this procedure only works for global solutions,
since (<ref>)
requires $x_0$ to be an arbitrary point in $\R^n$.

The bootstrap method
can also be applied to fully nonlinear operators
when the equation is satisfied in a bounded domain with good exterior data;
see Theorem 1.5 in <cit.>
and Theorem 6 in <cit.> to implement the bootstrap
regularity of Schauder type.

For equations satisfied on bounded domains,
when the operator is built, roughly speaking, by the sum of fractional Laplacians of different orders
that include the classical Laplacian, then the solutions
are typically $C^\infty$, in view of the regularizing effect of the
higher order operator. In this setting
one can include also some nonlinear terms. A simple example consists of $u\in C^{2s+\beta}(B_1)$
being a solution of the equation

$$ -\Delta u(x)+\widetilde F\big( (-\Delta)^s u(x)\big)=0\qquad{\mbox{for all }}x\in B_1,$$

with $\widetilde F\in C^\infty(\R)$.
We then have that the map $ x\mapsto f(x):=-\widetilde F\big( (-\Delta)^s u(x)\big)$
is locally Hölder continuous
and hence we can apply
the classical Schauder theory to
obtain that second derivatives of $u$ are
Hölder continuous in $B_{1/2}$.
Then, by
bootstrapping, we
conclude that $u$ has as many derivatives as we wish.
We remark that this setting is a particular case of that in (<ref>)
by choosing $F(p_1,p_2):=p_1+\widetilde F(p_2)$ and $s_1:=1$.

There are other special situations in which $C^\infty$
solutions in a bounded domain can be constructed, as established in
Theorem 1.1 of <cit.>.
This paper provides cases of nonlocal fully nonlinear
equations (rather concrete ones) whose Dirichlet problems
possess a
unique and
smooth solution. More specifically,
as detailed in Definitions 2.3 and 2.10 in <cit.>,
one can consider “nice weights”
which make the bootstrap regularity compatible with the convex
structure of the equation. This setting provides
smooth solutions for
concave elliptic operators acting
on integral expressions of the form

$$ \int_{\R^n}\big(u(x+y)+u(x-y)-2u(x)\big)\,\frac{y_i\,y_j\,\rho(x,y)}{|y|^{n+2s+2}}\,dy,$$

where $\rho$ is smooth, bounded, and
bounded away from zero, with derivatives of order $j$
in the variable $y$ bounded by $C_j|y|^{-j}$ for some constant $C_j$.
In particular, Theorem 1.1 in <cit.>
gives the existence of a $C^\infty$ solution of

$$ \begin{cases}\displaystyle
\int_{\R^n}\big(u(x + y) + u(x - y) - 2u(x)\big)\frac{
\rho(x,y)
}{|y|^{n+2s}}\,dy=f(x) & {\mbox{ for all }}x\in B_1,
\\u=g&{\mbox{ in }}\R^n\setminus B_1,
\end{cases}$$

provided that $f$ is $C^\infty$,
$g$ is bounded and uniformly continuous, and $\rho$ is a nice weight
as above.

In general, however, the nonlocal setting is not expected
to always provide $C^\infty$ solutions to general
fully nonlinear equations in bounded domains. This is due to the fact that
the linearized equation exhibits coefficients which
depend on the global data and are in general not better
than Hölder continuous (no matter how regular the solution
is in the interior of the domain).
As a consequence, the Schauder theory cannot be applied to
bootstrap regularity. As a matter of fact, the higher regularity
theory for fully nonlinear elliptic equations of nonlocal type in bounded domains
is, at the moment,
a field of research still under investigation.
It would be desirable to understand natural assumptions guaranteeing bootstrap regularity of $ C^\infty$ type.

§ TOWARDS A BERNSTEIN TECHNIQUE FOR VISCOSITY SOLUTIONS

Recall that the gradient bounds for integro-differential equations proved in
<cit.> and <cit.>, through the Krylov-Safonov and the Ishii-Lions doubling variables methods respectively, are carried out, in those papers, within the viscosity solutions framework. In the local second order theory, both methods had also been extended to the viscosity setting: see <cit.> for instance, and see <cit.>, respectively. In these last two papers, Barles developed a very general viscosity framework to get derivative estimates using the Ishii-Lions method of doubling variables. It is called there
“weak Bernstein's method”; “weak” since it applies to weak solutions in the viscosity sense, and “Bernstein's” since the auxiliary functions with doubled variables lead to some structural assumptions on the equation that are needed to obtain estimates and are closely related to those in
Bernstein's technique.

Since the Bernstein method towards gradient estimates is extended
to the nonlocal framework in the current paper for the first time,
we chose to simplify the exposition working with smooth solutions.
In this section, however,
we sketch some ingredients from our
forthcoming work <cit.>
in which we will extend the theory to get
regularity and derivative bounds for viscosity solutions.
To accomplish such a task using the Bernstein technique, one must replace the use of derivatives by incremental quotients within the auxiliary functions. Even in the local second order case, we could not find such an extension within the literature. As we will see next, the use of incremental quotients requires to introduce a new, averaged, zeroth order term within the auxiliary function.

As a byproduct of the extension to viscosity solutions,
the use of incremental quotients within the auxiliary functions will allow us to include the optimal class for Bellman and Isaacs-type fractional equations, as in the settings of <cit.>.
In doing this, an important ingredient
will come from Section 5 of <cit.>, which developed
the theory of viscosity solutions in a general
nonlocal and nonlinear framework. This work showed, in particular, that the
difference of two viscosity
solutions belongs to the appropriate
class of subsolutions (this is needed to
apply the incremental quotients method, by comparing
a solution with a small translation of itself) and that
a convenient maximum principle holds for such subsolutions.

In <cit.>, given a continuous function $u$, a unit vector $e\in\R^{n}$,
a smooth function $\eta$, $h\in\left(-\frac12,\frac12\right)$
and a positive constant $\sigma$, we will
consider the auxiliary function

\begin{equation}\label{PHICOME3S}
\varphi_{h,e}(x):= \eta^2(x) \,
\big(\delta_{h,e} u(x)\big)^2 + \sigma \int_0^1 u^2(x+the) \, dt ,\end{equation}

$$\delta_{h,e} u(x) :=\frac{ u(x+he) -u(x) }h$$

is an incremental quotient. This is in analogy with the classical setting in (<ref>), where here we replaced the derivative of $u$ by an incremental quotient and, as the most important novelty, we averaged the values of $u^{2}$ in a segment on the given direction.

\begin{eqnarray*}
(\delta_{h,e}u(x))^2=\left|\int_0^1 \partial_e u(x+the)\,dt\right|^2 \le
\int_0^1 |\nabla u(x+the)|^2\,dt,
\end{eqnarray*}

the computations in (<ref>) and (<ref>)
for a harmonic function $u$ in $B_1$,
are now replaced by the following ones:

\begin{equation}\label{FACIL-N}\begin{split}
\Delta\varphi_{h,e} \;&=\;
2|\nabla \eta|^2 (\delta_{h,e} u)^2+2\eta \Delta\eta (\delta_{h,e} u)^2+
8\eta \delta_{h,e} u\nabla \eta \cdot \nabla \delta_{h,e} u\\
&\qquad+2\eta^2|\nabla \delta_{h,e} u|^2+
2\eta^2\delta_{h,e} u \delta_{h,e} \Delta u
+2 \sigma \int_0^1 \big(|\nabla u|^2+ u \Delta u\big)(x+the)\,dt\\&=\;
2|\nabla \eta|^2 (\delta_{h,e} u)^2+2\eta \Delta\eta (\delta_{h,e} u)^2+
8\eta \delta_{h,e} u\nabla \eta \cdot \nabla \delta_{h,e} u\\
&\qquad+2\eta^2|\nabla \delta_{h,e} u|^2
+2 \sigma \int_0^1 |\nabla u(x+the)|^2\,dt
\\ &\ge\;
-6|\nabla \eta|^2 (\delta_{h,e} u)^2+
2\eta \Delta\eta (\delta_{h,e} u)^2+2 \sigma \int_0^1 |\nabla u(x+the)|^2\,dt\\

in $B_{1/2}$, as long as $|h|<1/2$
and $\sigma \geq
C_{n}\|\eta\|^2_{C^2(\R^{n})}$ for some positive constant $C_{n}$ depending only on $n$ (in fact, on the way the $C^{2}$-norm of a function in $\R^{n}$ is defined).

We observe that not only the computations in (<ref>)
extend the classical ones in (<ref>) and (<ref>) to deal with incremental quotients in a way which appears to be new with respect to the existing literature, but it is also flexible in regard to the positive part. That is, if we replace the function in (<ref>) by

\begin{equation*}
\psi_{h,e}(x):= \eta^2(x) \,
\big(\delta_{h,e} v(x)\big)^2_+ + \sigma \int_0^1 v^2(x+the) \, dt ,\end{equation*}

then, instead of (<ref>), for $v$ harmonic in $B_1$ we have

\begin{equation*}\begin{split}
\Delta\psi_{h,e} \;&=\;
2|\nabla \eta|^2 (\delta_{h,e} v)_+^2+2\eta \Delta\eta (\delta_{h,e} v)_+^2+
8\eta (\delta_{h,e} v)_+\,\nabla \eta \cdot \nabla \delta_{h,e} v\\
&\quad+2\eta^2\chi_{\{ \delta_{h,e} v>0\}}\,|\nabla \delta_{h,e} v|^2+
2\eta^2(\delta_{h,e} v )_+\,\delta_{h,e} \Delta v
+2 \sigma \int_0^1 \big(|\nabla v|^2+ v \Delta v\big)(x+the)\,dt\\&=\;
2|\nabla \eta|^2 (\delta_{h,e} v)_+^2+2\eta \Delta\eta (\delta_{h,e} v)_+^2+
8\eta (\delta_{h,e} v)_+\,\nabla \eta \cdot \nabla \delta_{h,e} v\\
&\quad+2\eta^2\chi_{\{ \delta_{h,e} v>0\}}\,|\nabla \delta_{h,e} v|^2
+2 \sigma \int_0^1 |\nabla v(x+the)|^2\,dt
\\ &\ge\;
-6|\nabla \eta|^2 (\delta_{h,e} v)_+^2+
2\eta \Delta\eta (\delta_{h,e} v)_+^2+2 \sigma \int_0^1 |\nabla v(x+the)|^2\,dt\\

again as long as $\sigma \geq
and $|h|<1/2$.

These computations will work similarly for the $s$-harmonic
extension of a function, thus extending
the theory developed
in the current paper to the setting of viscosity solutions to some
fractional equations.

Instead, additional arguments will be needed to deal with
integro-differential operators not having a local extension problem.
For instance, integration by parts will be replaced by the identity

$$ \int_{\R^n} \delta_{h,e} F(x)\,G(x)\,dx=-
\int_{\R^n} F(x)\, \delta_{-h,e}G(x)\,dx.$$

More importantly, the computations involved in the proof of the equivalent inequalities
stated in
Propositions <ref> and <ref>
need to be modified in this case. That is, defining

$$ u_{h,e}(x):=\left(
\int_0^1 u^2(x+the) \, dt\right)^{\frac12},$$

one checks that

\begin{eqnarray*}
&& {\mathcal{L}}_K \Big( \eta^2(x)\,
\big(\delta_{h,e} u(x)\big)^2+\sigma
\,u_{h,e}^2(x)\Big)
-2\eta^2(x)\,\delta_{h,e} u(x)\,{\mathcal{L}}_K
\delta_{h,e} u(x)
\\&&\qquad\qquad-2\sigma\,u_{h,e}(x)\,{\mathcal{L}}_K u_{h,e}(x)
&&\qquad= \int_{\R^n}
\Big( \eta^2(x)\,\big(\delta_{h,e} u(x)\big)^2
-\eta^2(y)\,\big(\delta_{h,e} u(y)\big)^2\Big)\,K(x-y)\,dy\\&&\qquad\qquad
+\sigma\int_{\R^n} \big(u_{h,e}^2(x)-u_{h,e}^2(y)\big)\,K(x-y)\,dy\\
&&\qquad\qquad-2\eta^2(x)\, \delta_{h,e} u(x)\,
\int_{\R^n}\big(\delta_{h,e} u(x)
-\delta_{h,e} u(y)\big)
\,K(x-y)\,dy\\
&&\qquad\qquad-2\sigma\,u_{h,e}(x)\,\int_{\R^n} \big(u_{h,e}(x)-u_{h,e}(y)\big)\,K(x-y)\,dy
\\&&\qquad\qquad
\\&&\qquad= \int_{\R^n}
\Big( \eta^2(x)\,\big(\delta_{h,e} u(x)\big)^2
-\eta^2(y)\,\big(\delta_{h,e} u(y)\big)^2\Big)\,K(x-y)\,dy
\\&&\qquad\qquad-2\eta^2(x)\, \delta_{h,e} u(x)\,
\int_{\R^n}\big(\delta_{h,e} u(x)-
\delta_{h,e} u(y)\big)
\,K(x-y)\,dy\\
&&\qquad= \int_{\R^n}\Big( 2\eta^2(x)\,
\delta_{h,e} u(x)\, \delta_{h,e} u(y)
-\eta^2(x)\, \big(\delta_{h,e} u(x)\big)^2-
\eta^2(y)\, \big(\delta_{h,e} u(y)\big)^2
\Big)\,K(x-y)\,dy
\\ &&\qquad=2\int_{\R^n} \eta(x)\,\big(\eta(x)-\eta(y)\big)\,
\delta_{h,e} u(x)\, \delta_{h,e} u(y)
\,K(x-y)\,dy
\\&&\qquad\qquad-
\int_{\R^n}\big|
\eta(x)\, \delta_{h,e} u(x)-
\eta(y)\,\delta_{h,e} u(y)
\big|^2\,K(x-y)\,dy.
\end{eqnarray*}

Given $\CALERRE\in\R$, this shows the equivalence between the pointwise

\begin{equation*}
{\mathcal{L}}_K \varphi_{h,e}\le 2\eta^2\,\delta_{h,e} u\,{\mathcal{L}}_K (\delta_{h,e} u)+2\sigma u_{h,e}\,{\mathcal{L}}_K u_{h,e}+\CALERRE
\end{equation*}

\begin{equation*}
\begin{split}&
2\int_{\R^n}\eta(x)\,\big( \eta(x)-\eta(y)\big)\,\delta_{h,e} u(x)\,\delta_{h,e} u(y)\,K(x-y)\,dy
\\ &\qquad\le\,\int_{\R^n}\big|\eta(x)\,\delta_{h,e} u(x)-\eta(y)\,\delta_{h,e} u(y)\big|^2\,K(x-y)\,dy
\\&\qquad\qquad+\sigma\int_{\R^n}\big|u_{h,e}(x)- u_{h,e}(y)\big|^2\,K(x-y)\,dy

§ LINEARIZED OPERATOR AND A MAXIMUM ESTIMATE

§.§ The linearized operator

Here we study the linearized equation
associated with the nonlinear problems (<ref>) and
(<ref>). To this end,
we denote by ${\mathcal{L}}_1,\dots,{\mathcal{L}}_J$ the linear operators ${\mathcal{L}}_{\mu_1},\dots,
{\mathcal{L}}_{\mu_J}$, respectively ${\mathcal{L}}_{K_1},\dots,
{\mathcal{L}}_{K_J}$, in these equations.
We will always assume that the convexity and ellipticity conditions (<ref>),
and (<ref>) are satisfied,
for some constant $\vartheta_0>0$.

We set $g$ as in (<ref>), and, given a function $u$,
we use the short notation

\begin{equation}\label{defalphak} \alpha_j:=
\alpha_j
\big({\mathcal{L}}_1u(x)-g_1(x),\dots,{\mathcal{L}}_Ju(x)-g_J(x)\big)\end{equation}

we consider the operator

\begin{equation}\label{L:EQ:OP}
Lv := \sum_{j=1}^J \alpha_j \,{\mathcal{L}}_{j}v.\end{equation}

It will be important that the functions $\alpha_{j}$ are defined at all points of our domain —and not only almost everywhere— since we will need to evaluate them at a maximum point of an auxiliary function. This will always be
possible since ${\mathcal{L}}_j u$ and $g_{j}$ will be finite and well defined at all points, by the regularity assumed on $u$ and since  $g_{j}$ are continuous functions.

The relevance of the linearized operator $L$
is given by the fact that the solution and its derivatives
satisfy suitable inequalities with respect to $L$, as stated in
the following result. Here, we remark that since $g_j$ is not better than semiconcave, its first and second derivatives only exist almost everywhere.

Let $f$ and $g_j$
be Lipschitz functions in $B_R$ for $j=1,\dots,J$,
$u\in C^\infty(B_R)\cap W^{1,\infty}(\R^n)$ be a solution of (<ref>) or of (<ref>) everywhere
in $B_R$, and $e\in\R^n$ satisfy $|e|=1$.
Let also $\gamma_1$ and $\gamma_2$ be defined right before the statement of Theorem <ref>.

\begin{eqnarray*}
&& Lu\ge f+g
\quad{\mbox{ everywhere in }}B_R\end{eqnarray*}

\begin{eqnarray*}
&& L\,\partial_{e} u
= \partial_e f+
\sum_{j=1}^J \alpha_j\partial_e g_j
\quad{\mbox{ almost everywhere in }}B_R

$$ | L\,\partial_{e} u|\le\gamma_1
\quad{\mbox{ everywhere in }}B_R.$$

If in addition $f$ and $g_j$
are semiconcave in $B_R$ and $u\in W^{2,\infty}(\R^n)$, then

\[ L\,\partial_e^2 u \le \partial_e^2f+
\sum_{j=1}^J \alpha_j\partial_e^2 g_j
\quad{\mbox{ almost everywhere in }}B_R.
\]

$$ L\,\partial_e^2 u\le\gamma_2
\quad{\mbox{ everywhere in }}B_R.$$

We prove the first claim
by using the equation in
either (<ref>) or (<ref>), in combination
with the convexity assumption (<ref>), exploited here
with $q:=(-g_1,\dots,-g_J)$ and $p:=({\mathcal{L}}_1 {u}-g_1,
\dots,{\mathcal{L}}_J {u}-g_J)$. In this way,
recalling the definition of $g$ in (<ref>), we have that, everywhere in $B_R$,

\begin{eqnarray*} -g
&=& F(-g_1,\dots,-g_J)-F({\mathcal{L}}_1 {u}-g_1,
\dots,{\mathcal{L}}_J {u}-g_J)+f
\\&\ge& f-\sum_{j=1}^J\alpha_j({\mathcal{L}}_1 {u}-g_1
,\dots,{\mathcal{L}}_J {u}-g_J)\,
\end{eqnarray*}

as desired.

Now we prove the second claim. To this end,
we let $\e>0$ and exploit the convexity assumption (<ref>)
at $x\in B_R$,

\begin{equation}\label{qep:1} q:=
({\mathcal{L}}_1u(x\pm \e e)-g_1(x\pm \e e)
,\dots,{\mathcal{L}}_Ju(x\pm \e e)-g_J(x\pm \e e))
\end{equation}

\begin{equation}\label{qep:2}

Using either (<ref>) or (<ref>), this leads to

\begin{eqnarray*}
\frac1\e\big(f(x\pm \e e)-f(x)
\big)
\Big\{
F\big({\mathcal{L}}_1u(x\pm \e e)
-g_1(x\pm \e e),\dots,{\mathcal{L}}_Ju(x\pm \e e)-g_J(x\pm \e e)\big)
\\&&\qquad-F\big({\mathcal{L}}_1u(x)
\Big\}
\\
\alpha_j
\big({\mathcal{L}}_1u(x)-g_1(x),\dots,{\mathcal{L}}_Ju(x)-g_J(x)\big)\\
{\mathcal{L}}_j u(x\pm \e e)-{\mathcal{L}}_j u(x)
-g_j(x\pm \e e)+g_j(x)

Taking the limit as $\e\searrow0$, we thereby find that,
for almost
every $x\in B_R$,

\begin{eqnarray*}
\pm \partial_e f(x)&\ge&\pm\sum_{j=1}^J
\alpha_j
\big({\mathcal{L}}_1u(x)-g_1(x),\dots,{\mathcal{L}}_Ju(x)-g_J(x)\big)
\big({\mathcal{L}}_j\partial_e u(x)-\partial_eg_j(x)\big)\\&=&
\pm L\,\partial_eu(x)\mp\sum_{j=1}^J
\alpha_j \partial_eg_j (x).
\end{eqnarray*}

From this, we deduce the second claim.

Finally, we prove the third one.
For this,
exploiting again (<ref>)
with $q$ and $p$ as above,
we have that

\begin{eqnarray*}
&&\frac1{\e^2}\big(f(x+\e e)+f(x-\e e)-2f(x)
\big)\\&&\qquad=\frac{1}{\e^2}\,\Bigg\{\Big(
F\big({\mathcal{L}}_1u(x+ \e e)-g_1(x+ \e e),\dots,{\mathcal{L}}_Ju
(x+ \e e)-g_J(x+ \e e)\big)\\&&\qquad\qquad
\Big(
F\big({\mathcal{L}}_1u(x-\e e)-g_1(x-\e e),\dots,{\mathcal{L}}_J
u(x- \e e)-g_J(x-\e e)\big)\\&&\qquad\qquad
\\
\alpha_j
\big({\mathcal{L}}_1u(x)-g_1(x),\dots,{\mathcal{L}}_Ju(x)-g_J(x)\big)\\
&&\qquad\qquad\cdot\Bigg\{ \frac{
\big(
{\mathcal{L}}_ju(x+\e e)-{\mathcal{L}}_ju(x)\big)+
\big( {\mathcal{L}}_ju(x-\e e)-{\mathcal{L}}_ju(x)\big)
&&\qquad\qquad\qquad-\frac{\big(g_j(x+\e e)-g_j(x)\big)+\big(g_j(x-\e e)-g_j(x)\big)}{\e^2}\Bigg\}
\\
\alpha_j
\;\Bigg\{
\frac{
{\mathcal{L}}_ju(x+\e e)+
{\mathcal{L}}_ju(x-\e e)-2{\mathcal{L}}_ju(x)
\frac{
g_j(x+\e e)+
g_j(x-\e e)-2g_j(x)
\Bigg\}.
\end{eqnarray*}

Sending $\e\searrow0$, we conclude that $\partial_e^2f\ge\sum_{j=1}^J
\alpha_j
\,\Big( {\mathcal{L}}_j\partial^2_e u(x)-\partial^2_e g_j(x)\Big)$
for almost
every $x\in B_R$,
as desired.

The last two claims are obtained by the previous
computations by recalling the definitions of $\gamma_1$
and $\gamma_2$ given right before the statement of
Theorem <ref>.

§.§ A maximum estimate

Here we show that the linearized operator $L$ in (<ref>)
satisfies the maximum principle, as well as a quantitative maximum estimate in the case
of nonzero right-hand sides.
This is the content of the following result, which will be proved using a barrier function. As customary, we use the notation $C_{\rm b}(\R^n):=C(\R^n)\cap L^\infty(\R^n)$ to denote the space of bounded and continuous
functions over all $\R^n$.

Let $R>0$ and $\gamma_0\ge0$.
Let $\varphi\in C_{\rm b}(\R^n)\cap {{C}^{1,1}(B_R)}$ be a
function in $\R^n$ such that

\begin{equation}\label{98suidvyfzsdaq}
L\varphi\le\gamma_0\quad{\mbox{ everywhere in }}B_R,
\end{equation}

where $L$ is given by (<ref>)–(<ref>) for some $u\in C^\infty(B_R)\cap W^{1,\infty}(\R^n)$
and continuous functions $ g _1,\dots, g _J$ in $B_R$.

\begin{equation} \label{DClap}\sup_{{B_R}}\varphi \le
\sup_{\R^n\setminus{B_R}}\varphi +C_R\,\gamma_0.\end{equation}

Here, when ${\mathcal{L}}_{j}={\mathcal{L}}_{\mu_j}$
for $j\in\{1,\dots,J\}$,
the constant $C_R>0$
is of the form

\begin{equation}\label{CR pc-1}

for some constant $C$ which
depends only on $n$ and $\vartheta_0$.
Recall that $\vartheta_0$ is the constant in
hypothesis (<ref>).

Instead when ${\mathcal{L}}_{j}={\mathcal{L}}_{K_j}$ for $j\in\{1,\dots,J\}$,
with $K_1,\dots,K_J$ satisfying (<ref>) and (<ref>),
the constant $C_R>0$
is of the form

\begin{equation}\label{CR pc-2}

for some constant $C$ which
depends only on $n$, $\vartheta_0$, and
the structural constant $\CDUEDUE$ in (<ref>).

Concerning the statement (<ref>),
we stress that, in our applications, we will need to use Proposition <ref>
for $C_{\rm b}(\R^n)\cap {{C}^{1,1}(B_R)}$ functions $\varphi$
which are not $C^2$. For instance, in some cases the auxiliary function $\varphi$ will contain a term involving $(\partial_e^2u)_+^2$.
Note that this function is locally $C^{1,1}$ when $u$ is smooth, even if $(\partial_e^2u)_+$ is not $C^{1,1}$
in general.

For a $C_{\rm b}(\R^n)\cap {{C}^{1,1}(B_R)}$ function $\varphi$,
we now define
a precise meaning to (<ref>) at all points of the domain (here $B_R$),
obtaining a finite value for $L\varphi$.

First, we use the usual
principal value definition for all the integro-differential operators considered through the paper, since the regularity of $\varphi$ is sufficient
to compute the integrals (involved in $L\varphi$) in the principal value sense, obtaining a finite value.

Secondly, when $L\varphi$ involves computing the classical Laplacian —as it may be the case for ${\mathcal{L}}_\mu$ in (<ref>)—, we define $-\Delta\varphi$
through second incremental quotients as

\begin{equation}\label{CAV}
\end{equation}

which is finite since $\varphi\in C^{1,1}$.
As a matter of fact, in our applications, inequality (<ref>)
will be satisfied even when writing (<ref>) with a $\limsup$ instead of $\liminf$.

To establish Proposition <ref>, we start giving an auxiliary barrier function.

Let $s\in(0,1)$ and $K$ satisfy (<ref>)
and (<ref>), and let

\begin{equation}\label{BETA:DEF} \beta(x):= \left\{ \begin{matrix}
|x|^2 -2 & {\mbox{ if $x\in B_{10}$, }}\\
98 & {\mbox{ otherwise.}}
\end{matrix}
\right.\end{equation}

\begin{equation}\label{CASO:K}
{\mathcal{L}}_K \beta\le -c\quad{\mbox{everywhere in }}B_1,\end{equation}

for some constant $c>0$ depending only on $n$ and
the structural constant $\CDUEDUE$ in (<ref>).

Similarly, with ${\mathcal{L}}_\mu$ defined by (<ref>)–(<ref>),
we have that

\begin{equation}\label{CASO:MU}
{\mathcal{L}}_\mu \beta\le -c\quad{\mbox{everywhere in }}B_1,\end{equation}

for some constant $c>0$ depending only on $n$.

In particular, for every $s\in[0,1]$,

\begin{equation}\label{CASO:MU:LA}
(-\Delta)^s \beta\le -c\quad{\mbox{everywhere in }}B_1,\end{equation}

for some constant $c>0$ depending only on $n$.

Let us start by proving (<ref>).
Let $x\in B_1$.
We observe that if $z\in B_{9}$ then $|x\pm z|\le|x|+|z|<10$, and hence $
\beta(x\pm z)
=|x\pm z|^2-2
As a consequence, for every $x\in B_1$ and $z\in B_{9}$,

\begin{equation}\label{HneIJA-a}
\begin{split}& \beta(x+z)+\beta(x-z)-2\beta(x)
\\&\qquad\quad =\big(
|x|^2+|z|^2+2x\cdot z-2\big)+\big(
|x|^2+|z|^2-2x\cdot z-2\big)-2\big(
|x|^2-2\big)= 2|z|^2.\end{split}\end{equation}

\begin{equation}\label{891-234}
\int_{B_{9}} \Big( \beta(x+z)+\beta(x-z)-2\beta(x)\Big)\,K(z)\,dz=
2\int_{B_{9}} |z|^2\,K(z)\,dz.
\end{equation}

On the other hand, if $z\in \R^n\setminus B_{9}$,
we have that $|x\pm z|\ge|z|-|x|\ge8$.
Therefore, for every $z\in \R^n\setminus B_{9}$,
we have that $\beta(x\pm z)\ge62$.
Since $\beta(x)\le-1$, we find in this case that

\begin{equation*}
\int_{\R^n\setminus B_9} \Big( \beta(x+z)+\beta(x-z)-2\beta(x)\Big)\,K(z)\,dz
\ge 126\int_{\R^n\setminus B_9} K(z)\,dz.
\end{equation*}

This inequality
and (<ref>),
together with the symmetry (<ref>) of the kernel
and the lower bound (<ref>),
lead to

\begin{eqnarray*}
-{\mathcal{L}}_K \beta(x)&=&
\frac12\,\int_{\R^n}\big( \beta(x+z)+\beta(x-z)-2\beta(x)\big)\,K(z)\,dz\\
&\ge& \int_{B_{9}} |z|^2\,K(z)\,dz
+\int_{\R^n\setminus B_9} K(z)\,dz\\&\ge&
\CDUEDUE\;s\,(1-s)\,\left(
\int_{B_{9}} \frac{dz}{|z|^{n+2s-2}}
+\int_{\R^n\setminus B_9} \frac{dz}{|z|^{n+2s}}
\right)\\&\ge&
\widetilde{C}\,\CDUEDUE,\end{eqnarray*}

for some constant $\widetilde C$ only depending on $n$. This gives the desired result (<ref>).

Now we prove (<ref>).
For this, we let $s\in(0,1)$
and apply (<ref>) with $K(z):=c_{n,s}\,|z|^{-n-2s}$, obtaining (<ref>)
with $c>0$ depending only on $n$
(since it is well known that for this kernel one may take $C_1$ to depend only on $n$).
We also notice that $-\Delta \beta =-2n$
and $\beta\le -1$ in $ B_1$, and therefore (<ref>)
is satisfied also for $s=0$
and $s=1$, up to renaming $c$. Consequently, the claim (<ref>) follows
the definition (<ref>) of ${\mathcal{L}}_\mu$.

Note that claim (<ref>),
that we have already proved,
can also be regarded as a particular case of (<ref>) by taking $\mu$ to be a Dirac's delta.

The following refinement of
Lemma <ref> deals with
rescaled versions of the function $\beta$ in (<ref>).
While the scale invariance of the fractional Laplacian
makes the results for
this kind of rescaled versions obvious,
some care is needed for the operators of indefinite order.

Let $s\in(0,1)$,
$K$ satisfy (<ref>)
and (<ref>), and let $\beta$ be as in (<ref>).
Given $R>0$, let

\begin{equation}\label{LAMujwns-A}
{\mathcal{L}}_K \beta_R\le -\frac{c}{R^{2s}}\quad{\mbox{everywhere in }}B_R,\end{equation}

for some constant $c>0$ depending only on $n$ and
the structural constant $\CDUEDUE$ in (<ref>).

Instead, with ${\mathcal{L}}_\mu$ defined by (<ref>)–(<ref>),
we have that

\begin{equation}\label{LAM:SE}
{\mathcal{L}}_\mu \beta_R\le -\frac{c}{1+R^2}\quad{\mbox{everywhere in }}B_R,\end{equation}

for some constant $c>0$ depending only on $n$.

We observe that claim (<ref>)
follows directly from (<ref>) applied to the rescaled kernel

\begin{equation}\label{AGOSTO}
K^{[ R ]}(z):=

Hence we focus on the proof of (<ref>). To this end,
we make use of (<ref>), according to which, for every $s\in[0,1]$,

$$ (-\Delta)^s \beta_R(x)= \frac{1}{R^{2s}}\,(-\Delta)^s \beta\left(\frac{x}{R}\right)
\le -\frac{c}{R^{2s}}$$

for all $x\in B_R$, where $c>0$ depends only on $n$.
We also point out that, by Young's inequality with exponents $1/s$ and $1/(1-s)$, we have
$R^{2s}\le sR^2+(1-s)
\le R^2+1$.
From these observations, we conclude that

\begin{eqnarray*}
{\mathcal{L}}_\mu \beta_R(x) =\int_0^1 (-\Delta)^s \beta_R(x)\,d\mu(s)
\le -c\,\int_0^1 \frac{d\mu(s)}{R^{2s}}\le -c\,\int_0^1 \frac{d\mu(s)}{1+R^{2}}
\end{eqnarray*}

which establishes (<ref>).

With the above barriers at hand, we are in the position of proving the maximum principle
with estimate which
is suitable for our goals.

We give a unified proof of (<ref>)
for both the cases in which ${\mathcal{L}}_{j}={\mathcal{L}}_{\mu_j}$
and ${\mathcal{L}}_{j}={\mathcal{L}}_{K_j}$. To this end,
it is convenient to define

\begin{equation}\label{STRU} C_R:=
\begin{cases}
C\,(1+R^2), & {\mbox{ if }}
C\,R^{2s}, &
{\mbox{ if }} {\mathcal{L}}_{j}={\mathcal{L}}_{K_j},
\end{cases}
\end{equation}

with $C$ to be chosen conveniently large.

Take $\beta$ as in Lemma <ref>
and define

\begin{equation}\label{09p29} \varphi_*(x):=\varphi(x)+ \frac{
\end{equation}

We claim that

\begin{equation}\label{8u8g8fia82}
\sup_{ \R^n} \varphi_*
=\sup_{ \R^n\setminus B_R} \varphi_*\,.
\end{equation}

Once this is proved, from
the fact that $-2\le\beta\le98$, we deduce

$$ \sup_{ B_R} \varphi\le\sup_{ B_R} \varphi_*
\sup_{ \R^n\setminus B_R} \varphi_*
\sup_{ \R^n\setminus B_R} \varphi +\frac{98\,C_R\,\gamma_0}{100}

and conclude (<ref>).

To prove claim (<ref>), it is enough to establish that

\begin{equation}\label{8u8g8fia82-PRE}
\sup_{ \R^n} \varphi_*
=\sup_{ \R^n\setminus B_{R'}} \varphi_*\,
\end{equation}

for every $R'\in(0,R)$, since $\varphi_*$ is continuous in $\R^n$.
For this,
we argue by contradiction and
assume that
there exists $x_*\in\overline{ B}_{R'}\subset B_R$ such that $ \varphi_*(x_*)=\sup_{\R^n}\varphi_*$.
This leads to

\begin{equation}\label{k:cas1-1}
\end{equation}

for all $j\in\{1,\dots,J\}$ and, in particular, to $(-\Delta)^s\varphi_*(x_*)\ge0$
for all $s\in(0,1)$.
This inequality holds also for $s=1$,
in view of (<ref>)
and that $x_*\in B_R$,
and for $

$$ \varphi_*(x_*)=\sup_{\R^n}\varphi_*\ge\sup_{x\in\R^n\setminus
+ \frac{
C_R\,\gamma_0}{100}\,\beta\left(\frac{x}{R}\right)\right) \geq0.$$

From all
this, we obtain that

\begin{equation}\label{k:cas1-2} {\mathcal{L}}_{\mu_j}
\varphi(x_*)\ge0\end{equation}

for all $j\in\{1,\dots,J\}$.

Now, collecting the inequalities (<ref>)
and (<ref>), and recalling
that $\alpha_j\ge0$ by (<ref>)
and the definition (<ref>)
of $L$,
we deduce

\begin{equation}\label{k:cas1-3} L\varphi(x_*)\ge0.\end{equation}

Let $\beta_R(x):=\beta(x/R)$ be as
in Lemma <ref>.
In light of (<ref>), (<ref>), and (<ref>),
we have that, for all $j\in\{1,\dots,J\}$,

$$ {\mathcal{L}}_{K_j} \beta_R\le -\frac{c\,C}{C_R}
\qquad{\mbox{ and }}\qquad
{\mathcal{L}}_{\mu_j} \beta_R\le -\frac{c\,C}{C_R}\qquad{\mbox{everywhere in }}B_R,$$

and consequently, recalling
the lower bound (<ref>) on $\alpha_j$,

$$ L\beta_R
\le
\le-\frac{c\,C\,\vartheta_0}{C_R}\qquad{\mbox{everywhere in }}B_R.$$

Here, $c$ is the constant in Lemma <ref>.
Making use of this inequality, (<ref>), and (<ref>), and observing that $x_*\in\overline{B}_{R'}\subset B_R$,
we conclude

$$ L\varphi(x_*)=L\varphi_*(x_*)- \frac{
C_R\,\gamma_0}{100}\,L\beta_R(x_*)\ge \frac{c\,C\,\vartheta_0\,\gamma_0}{100}.

This, combined with (<ref>), gives that $
which provides a contradiction if $C$ is taken large enough, depending only on $n$
and $\vartheta_0$, and additionally on
the structural constant $\CDUEDUE$ in (<ref>)
if ${\mathcal{L}}_{j}={\mathcal{L}}_{K_j}$,
in view of the dependences of $c$ stated in Lemma <ref>.
This argument establishes
claim (<ref>) and finishes the proof.

§ A UNIFIED APPROACH TOWARDS DERIVATIVE ESTIMATES

Most of our theorems
will follow from our next result, Theorem <ref>.
It deals with solutions of “linear operators
possibly varying from point to point and satisfying a maximum principle
with estimate”.
To state our result precisely, given $R>0$, we consider a family of
linear operators $\{L^{(x)}\}_{x\in B_R}$,
of any of the types that we have considered previously in the paper,
and we set the following terminology.

We say that the family $\{
L^{(x)}\}_{x\in B_R}$ satisfies the maximum principle with
estimate in $B_R$ if
there exists a constant $C_R$
such that,
for every $\gamma_0\ge0$ and every nonnegative function $
\varphi\in C_{\rm b}(\R^n)\cap {{C}^{1,1}(B_R)}$,
the following statement holds true:

\begin{equation}\label{MAXPART1}
\inf_{y\in B_R}\{ L^{(y)}\varphi(x)\}\le\gamma_0
\qquad{{\mbox{for all }}}\;x\in B_R,
\end{equation}

\begin{equation}\label{MAXPART}
\sup_{{B_R}}\varphi \le
\sup_{\R^n\setminus{B_R}}\varphi +C_R\gamma_0.
\end{equation}

We recall that the computation of $L^{(y)}\varphi$
with $\varphi$ only of class $C_{\rm b}(\R^n)\cap {{C}^{1,1}(B_R)}$ is intended in the light of Remark <ref>
(in particular, when $L^{(y)}$ involves classical second order operators,
the incremental quotient setting
in (<ref>) must be adopted).
This comment concerns (<ref>) and also the left-hand side
of (<ref>) below. Note, instead, that
the right-hand side of (<ref>)
is well defined and finite since the
functions $\partial^2_eu$, $\partial_eu$, and $u$
are assumed in the theorem to be
smooth and globally bounded.

We now
provide a result that will
establish first and one-sided second derivative
bounds in the case of operators possessing a local extension. It will serve both for the Pucci-type operators as well as for our other fully nonlinear equations
defined through the function $F$. In addition, the same method, properly modified, will be used for equations with no extension property.

We introduce the auxiliary function

\begin{equation}\label{78AJJQoeAAiak0IIA}
\varphi(x) :=\overline\eta_R^2(x)\,\big(\partial_e^2u(x)\big)^2_+
+ \tau R^{-2}\,{{\eta}}_R^2(x)\big(\partial_e u(x)\big)^2+
\sigma R^{-4}\,\left(u(x)-\sup_{B_R} u\right)^2,\end{equation}

where $\overline\eta_R$ and $\eta_R
$ are smooth functions with compact support in $B_{R/2}$
and $B_R$, respectively. In addition, we will need that $\eta_R=1$
in $B_{R/2}$ in order to be able to verify (in the next section)
the key inequality (<ref>) below.
As we will see in the proof of the following theorem, the first
and second derivative bounds (<ref>)
and (<ref>) will follow from making
two different appropriate choices of the
pair of cutoff functions $(\overline\eta_R,{\eta}_R)$.

Let us explain why we are forced,
in this nonlocal theory, to consider
the auxiliary function (<ref>)
involving both $\partial_e u$ and $
\partial_e^2u$ at the same time, as well as two different cutoff functions. This
has not been considered in the local theory. Indeed,
for local equations one works first with $\eta^2(\partial_e u)^2+\sigma
u^2$ as in (<ref>), obtaining first derivative
estimates in a ball. One considers next $\overline\eta^2(\partial_e v)_+^2+\sigma
v^2$ with $v=\partial_e u$,
as in (<ref>). This allows to control $\partial_e^2 u$
by above in a smaller ball by the previous control of $\partial_e u$
in the larger ball. However, since in the nonlocal setting the maximum
principle involves the whole exterior datum
(and not only the boundary datum), this second
choice would require to control the first derivatives $v=\partial_e u$
in all space, which cannot be achieved for an equation posed
in a ball. We solve this trouble by introducing the new
auxiliary function (<ref>).
Note anyway that the simpler choice (<ref>)
would still work if our nonlocal equation were posed in all of $\R^n$.
This is why we have chosen to refer to this simpler test function
in the Introduction, especially in order to make the Open problem <ref>
as simple as possible.

Let $\sigma\ge1$, $\tau\ge1$, and $R>0$.
Assume that $\{
L^{(x)}\}_{x\in B_R}$ is a family of linear operators
satisfying the maximum principle with estimate in $B_R$,
with constant $C_R$, according to Definition <ref>.

Given any $e\in\R^n$ with $|e|=1$,
$\overline\eta_R\in C^\infty_c(B_{R/2})$,
${\eta}_R\in C^\infty_c(B_R)$
with $\eta_R=1$ in $B_{R/2}$,
and $u\in C^\infty(\R^n)\cap W^{2,\infty}(\R^n)$,
consider the function $\varphi$ in (<ref>).

Assume that, for every such choice of $\sigma$, $\tau$, $R$, $e$, $\overline\eta_R$, $\eta_R$, and $u$,
we have

\begin{equation} \label{L-EQAU-B1}\begin{split}
\le\,&
2\overline\eta_R^2(x)\,(\partial_e^2 u(x))_+\,L^{(x)}\partial_e^2 u(x)\\&\quad
+ 2\tau R^{-2}\,{\eta}_R^2(x)\,\partial_e u(x)\,
L^{(x)}\partial_e u(x)+2\sigma R^{-4}\,\left(u(x)-\sup_{B_R}u\right)\,
L^{(x)} {\left(u-\sup_{B_R}u\right)}(x)
\end{split}\end{equation}

every $x\in B_R$.

We then have,
for every $u\in C^\infty(\R^n)\cap W^{2,\infty}(\R^n)$,

\begin{equation}\label{LAPRIMAGENE}
\sup_{B_{R/2}} |\partial_e u|\le
C\,\Big(C_Ra_1+R^{-1} \big(C_R a_0\|u\|_{L^\infty(B_R)}\big)^{1/2}+
\Big)\end{equation}

\begin{equation}\label{LASECOGENE}
\sup_{B_{R/4}} \partial^2_e u\le C\,\Big(
C_R a_2+C_R R^{-1} a_1
+R^{-2}\big(C_R a_0 \| u\|_{L^\infty(B_R)}\big)^{1/2}+R^{-2}\| u\|_{L^\infty(\R^n)}
\Big),\end{equation}

\begin{equation}\label{1NOTA a 012}
a_0:=\sup_{x\in B_R} \left(L^{(x)}{\left(u-\sup_{B_R}u\right)}(x)\right)_- \, ,
\end{equation}

\begin{equation}\label{2NOTA a 012}
a_1:=\sup_{x\in B_R} \big|L^{(x)}\partial_e u( x)\big|,\end{equation}

\begin{equation}\label{3NOTA a 012}
a_2:=\sup_{x\in B_R}\big(
L^{(x)}\partial_e^2 u( x)\big)_+\,,
\end{equation}

and $C$ is a constant depending only on $n$,
$\sigma$, and $\tau$.

Note that the functions $\varphi$ in (<ref>)
belong to $C_{\rm b}(\R^n)\cap {{C}^{1,1}(B_R)}$, as required in the previous section,
by basic properties
of the positive part and the square power
in $\big(\partial_e^2u\big)^2_+$, and since $\partial_e^2u\in C^\infty(\R^n)$.

We let

\[
\Phi(x):=a_2\,
\overline\eta_R^2(x)\,(\partial_e^2 u(x))_+
+a_1R^{-2} {\eta}_R^2(x)\,|\partial_e u(x)|
+a_0R^{-4} \| u\|_{L^\infty(B_R)}.\]

We claim that

\begin{equation}\label{9OA-CLAI}
\sup_{{B_R}}\varphi
\le C_\sharp\left(C_R
\sup_{B_R} \Phi+R^{-4}\| u\|_{L^\infty(\R^n)}^2
\right)\end{equation}

for some constant $C_\sharp\ge1$ depending only on $n$,
$\sigma$, and $\tau$.
To check this, we use (<ref>) to find that, for all $x\in B_R$,

\begin{eqnarray*}
\inf_{y\in B_R}\{ L^{(y)}\varphi(x)\}&\le&
\\&\le&
2\overline\eta_R^2(x)\,(\partial_e^2 u( x))_+
\,L^{( x)}\partial_e^2 u( x)
+ 2\tau R^{-2}\, {\eta}_R^2( x)\,\partial_e u( x)\,
L^{( x)} \partial_e u( x)\\&&\qquad+2\sigma R^{-4}\,\left(u( x)-
\sup_{B_R}u\right)L^{( x)} {\left(u-\sup_{B_R}u\right)}(x)\\&\le&
\sup_{B_R}\Big(
2a_2\overline\eta_R^2\,(\partial_e^2 u)_+
+ 2a_1\,\tau R^{-2}\,{\eta}_R^2\,|\partial_e u|\Big)
+4a_0\,\sigma R^{-4}\,\|u\|_{L^\infty(B_R)}
\\ &\le&C\sup_{B_R}\Phi,\end{eqnarray*}

for some constant $C$ depending only on $n$, $\sigma$, and $\tau$.
Hence, (<ref>) is satisfied with $
\gamma_0:=
C\sup_{B_R}\Phi$, and thus (<ref>) yields that $
\sup_{{B_R}}\varphi \le
\sup_{\R^n\setminus{B_R}}\varphi +C_R\gamma_0$.
As a consequence, recalling (<ref>),
we have that $
\sup_{{B_R}}\varphi \le 4\sigma\,R^{-4}\|u\|_{L^\infty(\R^n)}^2
This establishes (<ref>), as desired.

We stress that the constant $C_\sharp$ in (<ref>)
does not depend on $\eta_R$ and $\overline\eta_R$, and hence
we can now take appropriate choices for these two functions.
Note first that, from (<ref>),

\begin{equation*}\begin{split}
&\sup_{{B_R}} \Big( \overline\eta_R^2(\partial_e^2 u)_+^2+
\tau R^{-2}\,\eta_R^2|\partial_e u|^2\Big)
\\
\sup_{{B_R}}\varphi\\&\qquad
\le\, C_\sharp\left(C_R
\sup_{B_R} \Phi+R^{-4}\| u\|_{L^\infty(\R^n)}^2
\right)\\ &\qquad
\le\, C_\sharp\left(C_R
\Big( \sup_{B_R}\big( a_2\overline\eta_R^2 \,(\partial_e^2 u)_+
+R^{-2}a_1{\eta}_R^2\,|\partial_e u|\big)
+R^{-4}a_0 \| u\|_{L^\infty(B_R)}\Big)+R^{-4}\| u\|_{L^\infty(\R^n)}^2
\right)
\\ &\qquad
\le\, C_\sharp\Bigg(\frac1{2C_\sharp}\sup_{B_R}\Big(\overline\eta_R^2 \,(\partial_e^2 u)_+^2 +
\tau R^{-2} {\eta}_R^2\,|\partial_e u|^2\Big)
\\&\qquad\qquad+C_\sharp C_R^2 a_2^2\sup_{B_R}\overline\eta_R^2
+C_\sharp C_R^2 R^{-2}a_1^2\sup_{B_R}\eta_R^2
+C_RR^{-4}a_0 \| u\|_{L^\infty(B_R)}+R^{-4}\| u\|_{L^\infty(\R^n)}^2
\Bigg),
\end{split}\end{equation*}

where we have used a
Cauchy-Schwarz inequality in the last step,
and therefore

\begin{eqnarray*}&& \frac12\,
\sup_{{B_R}} \Big( \overline\eta_R^2(\partial_e^2 u)_+^2+\tau R^{-2}\,
\eta_R^2|\partial_e u|^2\Big)\\&&\qquad\le
C_R^2 a_2^2\sup_{B_R} \overline\eta_R^2
+C_R^2 R^{-2}a_1^2\sup_{B_R}\eta_R^2
+C_RR^{-4}a_0 \| u\|_{L^\infty(B_R)}+R^{-4}\| u\|_{L^\infty(\R^n)}^2
\right).\end{eqnarray*}

We now make our two choices of cutoff
functions. First, choosing $\overline\eta_R:=0$
and $\eta_R$ with $|\eta_R|\le1$ and $\eta_R=1$
in $B_{R/2}$, from the last inequality
we obtain (<ref>).
Next, with the same choice of $\eta_R$, we now choose $\overline\eta_R\in C^\infty_c(B_{R/2})$
satisfying $\overline\eta_R=1$ in $B_{R/4}$. In this way, we infer (<ref>).

§ OPERATORS HAVING LOCAL EXTENSION

In this section, we prove the main results of
Subsections <ref> and <ref>,
Theorem <ref>, Theorem <ref>, and
Corollary <ref>.
To this end,
we will exploit the key inequality of Proposition <ref>
in a more general version, as given
in the forthcoming
Proposition <ref>. To state it,
we consider a direction $e\in\R^n$, with $|e|=1$,
and two cutoff functions $\eta$ and $\overline\eta$ such that

\begin{equation}\label{118}
\eta \text{ and } \overline{\eta}\text{ belong to } C^\infty_c (B_{1/2})
\text{ and }
{\mbox{$\eta(x)=1$ for all~$x\in B_{1/4}$.}}\end{equation}

We then define

\begin{equation}\label{Rq2sed234}
{{\eta}}_R(x):=\eta\left( \frac{x}{2R}\right)
\qquad{\mbox{ and }}\qquad
\overline\eta_R(x):=\overline{{{\eta}}}\left( \frac{x}{R}\right)

Note that this setting covers the situation required in Theorem <ref>.

Recall the comments after Definition <ref>
for the precise meaning of the left-hand side of (<ref>) below (when $s=1$),
and note that the proposition is uniform in $s$ (in fact, we include the cases $s=0$
and $s=1$).

Let $\eta_R$ and $\overline\eta_R$
be as in (<ref>)-(<ref>),
$s\in[0,1]$, and $\kappa
\in\R$. Let $u\in C^\infty(\R^n)\cap W^{2,\infty}(\R^n)$.
For every $x\in\R^n$, let

$$ \varphi(x) :=\overline\eta_R^2(x)\,\big( \partial_e^2u(x)\big)^2_+
+ \tau\,R^{-2}\,\eta_R^2(x)\big(\partial_e u(x)\big)^2+
\sigma\,R^{-4}\,\big(u(x)-\kappa\big)^2.$$

Then, there exist positive constants $\sigma_0$ and $\tau_0$, depending only on $n$,
and $\|\overline{{{\eta}}}\|_{C^2(\R^n)}$, such that

\begin{equation}\label{65GHA:98a124gH}
\begin{split}
& (-\Delta)^{s} \varphi \le
2\overline\eta_R^2\,(\partial_e^2 u)_+\,
(-\Delta)^s \partial_e^2 u+ 2\tau\,R^{-2}\,{\eta}_R^2\,\partial_e u\,(-\Delta)^s\partial_e u\\&\qquad\qquad\qquad
+2\sigma \,R^{-4}\,(u-\kappa)\,(-\Delta)^s {(u-\kappa)}\qquad\qquad\qquad
{\mbox{if~$\tau\ge \tau_0$ and~$\sigma\ge \sigma_0 \tau$,}}\end{split}\end{equation}

everywhere in all of $\R^n$.

The proof of Proposition <ref> will be given after Corollary <ref>
and will rely on the forthcoming auxiliary calculations.

§.§ Computations in the extended space

Throughout this section,
we consider $s\in(0,1)$ and $a:=1-2s\in(-1,1)$.
Given $u\in C^\infty(\R^n)\cap W^{2,\infty}(\R^n)$,
we take $U$ to be the $s$-harmonic extension of $u$
in $\R^{n+1}_+:=
\R^n\times(0,+\infty)=\R^n\times\R_+$, i.e.,
the unique
bounded solution of the problem

\begin{equation}\label{EQ:EXT:U}
\left\{
\begin{matrix}
{\rm div}\, \big(y^a \nabla U(x,y)\big)=0 & {\mbox{ for all }}(x,y)\in\R^{n+1}_+,\\
U(x,0)=u(x) & {\mbox{ for all }}x\in\R^n.
\end{matrix}
\right.\end{equation}

Such solution can be found, for instance, by convolving $u$ against
the fractional Poisson kernel; see <cit.>.
As a result, since $u$ is continuous,
we have that $U$ is continuous
in $\overline{\R^{n+1}_+}$.
It is also clear that $U$ is smooth in $\R^{n+1}_+$. Note that the solution $U$
can also be obtained by direct minimization of
an energy functional; see Section 3
in <cit.>.

The uniqueness of the bounded extension $U$ follows from Lemma <ref>, which is a more
general result that we will need later in this section. See
Corollary 3.5 in <cit.> for an alternative proof of uniqueness.

For notational convenience, we also set

\begin{equation}\label{defla}
L_a U :=- {\rm div}\, \big(y^a \nabla U\big).
\end{equation}

Of course, no confusion should arise with the
linearized operator $L$
introduced in (<ref>). For further reference, we point out that

\begin{equation}\label{AUS:fg}
L_a (VW) = (L_a V)\,W+V\,(L_a W) - 2y^a \nabla V\cdot \nabla W.
\end{equation}

With a slight abuse of notation, we identify the direction $e\in\R^n$
with the same
direction in $\R^{n+1}$
when we consider the directional derivatives along $e$.
In this setting, we have the
following result.

Let $\eta_R$ be as in (<ref>)-(<ref>) and $\kappa\in\R$.
For $(x,y)\in\R^{n+1}_+$, let

\big( U(x,y)-\kappa\big)^2$$

$$ \Psi_1(x,y):={{\eta}}_R^2(x)\,
\big( \partial_e U(x,y)\big)^2.$$

\begin{equation}
\label{LE:AUS:2} L_a \Psi_0 =- 2y^a \,|\nabla U|^2
\end{equation}

\begin{equation} \label{87t58-0dst}
L_a \Psi_1 \le- y^a\,{{\eta}}^2_R\,\big| \nabla\partial_e U\big|^2
+C\chi_{B_R}\,R^{-2}\,y^a\big( \partial_e U\big)^2
\end{equation}

in $\R^{n+1}_+$,
for some constant $C$
depending only on $n$ and $\|\eta\|_{C^2(\R^n)}$.

The statement and proof of Lemma <ref>
concern a
general constant $\kappa\in\R$. Later it will be convenient
to choose $\kappa$ to be the supremum of the solution.

By (<ref>)
we have that $ L_a \Psi_0=2(L_a U)\,(U-\kappa)-2y^a \,|\nabla U|^2$,
whence (<ref>) plainly follows.

Now, to prove (<ref>),
since $L_a U=0$, and thus $L_a \partial_eU=0$,
we infer that

\begin{equation} \label{LE:AUS:2:ZZ}
L_a( \partial_e U)^2 = -2y^a \,|\nabla \partial_e U|^2.
\end{equation}

Furthermore, by (<ref>),

\begin{equation}\label{HG:A1}
L_a {{\eta}}_R^2 = 2{{\eta}}_R\, (L_a
{{\eta}}_R )- 2y^a |\nabla{{\eta}}_R|^2 =-
2y^a {{\eta}}_R\,\Delta_x{{\eta}}_R- 2y^a |\nabla_x{{\eta}}_R|^2.
\end{equation}

Now, we notice that $\Psi_1={{\eta}}_R^2( \partial_e U)^2 $.
Accordingly, we
use (<ref>), (<ref>), and (<ref>)
to see that

\begin{eqnarray*}
L_a\Psi_1 &=&(
L_a {{\eta}}_R^2)\,( \partial_e U)^2 +{{\eta}}_R^2\,(L_a ( \partial_e U)^2  )-
2y^a \nabla_x {{\eta}}_R^2\cdot \nabla_x( \partial_e U)^2  \\
&=& -2y^a
+|\nabla_x{{\eta}}_R|^2)\,\big( \partial_e U\big)^2
-2y^a \,{{\eta}}_R^2\,|\nabla \partial_e U|^2
\\&&\qquad- 8y^a \,{{\eta}}_R\, \partial_e U\,\nabla_x{{\eta}}_R
\cdot \nabla_x \partial_e U.
\end{eqnarray*}

Also, we observe that, by the Cauchy-Schwarz inequality,

$$ 8y^a\, \big| {{\eta}}_R\,\partial_e U\,\nabla_x{{\eta}}_R
\cdot \nabla_x \partial_e U\big|\le
y^a \, {{\eta}}_R^2\,|\nabla\partial_e U|^2+
16\, y^a\,|\nabla_x{{\eta}}_R|^2\, (\partial_e U)^2.$$

Consequently, we obtain that

$$ L_a\Psi_1 \le- y^a \,{{\eta}}_R^2\, |\nabla\partial_e U|^2
+7\,|\nabla_x{{\eta}}_R|^2)\,\big( \partial_e U\big)^2.$$

Hence, the desired inequality (<ref>) follows,
since ${\eta}_R$ is supported in $B_R$,
$|\nabla_x{{\eta}}_R|\le \|\eta\|_{C^1(\R)}\,R^{-1}$,
and $|\Delta_x{{\eta}}_R|\le n\,\|\eta\|_{C^2(\R)}\,R^{-2}$.

We now start estimating the operator $L_a$
acting on the second derivatives of $U$.
The following (nonoptimal)
inequality is all what we will need subsequently.

Let $\overline\eta_R$ be as in (<ref>)-(<ref>).
For $(x,y)\in\R^{n+1}_+$,
$$\Psi_2(x,y):=\overline\eta_R^2 (x)\,\big( \partial^2_e U(x,y)\big)^2_+ .$$

Then, there exists a constant $C$,
depending only on $n$
and $\|\overline{{{\eta}}}\|_{C^2(\R^n)}$, such that

\begin{equation}\label{DR:ZZ}
C\,\chi_{B_{R/2}}\,R^{-2}\,y^a \,(\partial^2_e U)^2
\qquad{\mbox{ everywhere in }}\R^{n+1}_+.
\end{equation}

Concerning claim (<ref>),
we point out that,
since $U\in C^\infty(\R^{n+1}_+)$,
we have that $\partial_e^2U\in C^\infty(\R^{n+1}_+)$
and accordingly

\begin{equation}\label{ksd-we}
\big( \partial^2_e U\big)^2_+\in C^{1,1}_{\rm loc}(\R^{n+1}_+).\end{equation}

This is sufficient to compute the left-hand side of (<ref>)
everywhere by writing $L_a$
in nondivergence form and using
the incremental quotient limit definition in (<ref>). As an alternative,
one could also interpret (<ref>) in  the weak sense,
since (<ref>) yields that $\big( \partial^2_e U \big)^2_+\in H^1_{\rm loc}(\R^{n+1}_+)$.

Recalling (<ref>),
we have that

\begin{equation}\label{HG:A0:ZZ}
L_a\Big(\overline\eta_R^2 \,\big( \partial^2_e U\big)^2\Big)=
L_a \overline\eta_R^2\,\big( \partial^2_e U\big)^2+
\overline\eta_R^2\,L_a \big( \partial^2_e U\big)^2-
8y^a \overline\eta_R\,\partial^2_e U\,
\nabla \overline\eta_R\cdot\nabla \partial^2_e U
\end{equation}

and, since $L_a U=0$,

\begin{equation}\label{HG:A1:11}
L_a\big((\partial^2_e U)^2\big) = 2\partial^2_e U
\, L_a\partial^2_e U
- 2y^a |\nabla \partial^2_e U|^2=- 2y^a |\nabla \partial^2_e U|^2.
\end{equation}

$ L_a {{\overline\eta_R^2}}
2y^a {{\overline\eta_R}}\,\Delta_x{{\overline\eta_R}}-
2y^a |\nabla_x{{\overline\eta_R}}|^2$.
By inserting this and (<ref>) into (<ref>), we conclude that

\begin{equation}\label{9:AH:ZZ}
\begin{split}
L_a \Big(\overline\eta_R^2 \,\big( \partial^2_e U\big)^2\Big) \,&=\, -
2y^a \overline\eta_R\,\Delta_x\overline\eta_R\,(\partial^2_e U)^2
- 2y^a |\nabla_x\overline\eta_R|^2\,(\partial^2_e U)^2
- 2y^a \overline\eta_R^2\,|\nabla \partial^2_e U|^2\\&\qquad\qquad
-8 y^a \overline\eta_R \,\partial_e^2 U\,
\nabla_x \overline\eta_R\cdot\nabla_x\partial^2_e U\,.\end{split}\end{equation}

Now we use the
Cauchy-Schwarz inequality, finding that

$$ -8 y^a \overline\eta_R \,\partial_e^2 U
\nabla_{x} \overline\eta_R\cdot\nabla_x \partial^2_e U \le
y^a \,\overline\eta_R^2\, |\nabla \partial^2_e U|^2
+16\, y^a\, |\nabla_x\overline\eta_R|^2 \,(\partial_e^2 U)^2.$$

Since $\overline\eta_R$ is supported in $B_{R/2}$,
$|\nabla_x\overline\eta_R|\le \|\overline{{{\eta}}}\|_{C^1(\R)}\,R^{-1}$,
and $|\Delta_x\overline\eta_R|\le n\,\|\overline{{{\eta}}}\|_{C^2(\R)}\,R^{-2}$,
the latter estimate and (<ref>) give

\begin{equation}\label{CH964}
L_a \Big(\overline\eta_R^2 \,\big( \partial^2_e U\big)^2\Big)\le
C\,\chi_{B_{R/2}}\,R^{-2}\,y^a \,(\partial^2_e U)^2
- y^a \overline\eta_R^2\,|\nabla \partial^2_e U|^2\le
C\,\chi_{B_{R/2}}\,R^{-2}\,y^a \,(\partial^2_e U)^2\end{equation}

everywhere in $\R^{n+1}_+$, for a suitable constant $C$
depending only on $n$
and $\|\overline{{{\eta}}}\|_{C^2(\R^n)}$.

Now, we use this inequality to prove (<ref>).
For this, let $(x,y)\in\R^{n+1}_+$ and
we observe that if $\Psi_2(x,y)>0$, then $\Psi_2=\overline\eta_R^2 \,\big( \partial^2_e U\big)^2$
in a small neighborhood of $(x,y)$ and thus (<ref>)
follows in this case directly from (<ref>).

If instead $(x,y)\in\R^{n+1}_+$ is such that $\Psi_2(x,y)=0$,
the second incremental quotient
definition (<ref>) gives, since $\Psi_2\ge0$, that $
Also, since $\Psi_2\in C^1_{\rm loc}(\R^{n+1}_+)$ by (<ref>),
the fact that $\Psi_2\ge0=\Psi_2(x,y)$
leads to $\nabla\Psi_2(x,y)=0$
as a result,

\begin{eqnarray*}
\end{eqnarray*}

In particular, also
in this case $L_a\Psi_2(x,y)\le
C\,\chi_{B_{R/2}}\,R^{-2}\,y^a \,(\partial^2_e U)^2$. The proof
of (<ref>) is thereby complete.

By combining Lemmata <ref> and <ref> we obtain:

Let $\eta_R$ and $\overline\eta_R$
be as in (<ref>)-(<ref>) and $
\sigma$, $\tau$, $\kappa\in\R$.
For $(x,y)\in\R^{n+1}_+$, let

\begin{equation}\label{DEF:FI}
\Phi(x,y):=\overline\eta_R^2 (x)\,\big(\partial^2_e U(x,y) \big)_+^2 +
\tau\,R^{-2}\,{{\eta}}_R^2(x)\big( \partial_e U(x,y)\big)^2 +
\sigma\,R^{-4} \big(U(x,y)-\kappa\big)^2.\end{equation}

Then, there exist
positive constants $\sigma_0$ and $\tau_0$, depending only
on $n$, $\|\eta\|_{C^2(\R^n)}$,
and $\|\overline{{{\eta}}}\|_{C^2(\R^n)}$, such that $L_a \Phi\le0$
everywhere in $\R^{n+1}_+$
if $\tau\ge\tau_0$ and $\sigma\ge \sigma_0 \tau$.

In the notation of
Lemmata <ref> and <ref>
we have that $
\Phi= \Psi_2+\tau R^{-2}\Psi_1+\sigma R^{-4}\Psi_0$, and consequently

\begin{equation*}
\begin{split}
C\,\chi_{B_{R/2}}\,R^{-2}\,y^a \,(\partial^2_e U)^2
-\tau\, R^{-2}\,y^a\,{{\eta}}_R^2\,\big| \nabla\partial_e U\big|^2\\ &\qquad\qquad\qquad
+C\,\tau\,\chi_{B_R}\,R^{-4}\,y^a\big( \partial_e U\big)^2
-2\sigma\,R^{-4}\,y^a \,|\nabla U|^2.\end{split}\end{equation*}

Now, we notice that $ |\partial^2_e U|
\le |\nabla\partial_e U|$ and that, by (<ref>), $
\chi_{B_{R/2}}\le {{\eta}}_R^2$.
Therefore, if $\tau$ is sufficiently large
as stated in the corollary, we deduce

\begin{equation*}
L_a\Phi\;\le\; C\tau\,R^{-4}\,y^a\big( \partial_e U\big)^2
-2\sigma\,R^{-4}\,y^a \,|\nabla U|^2.\end{equation*}

Thus, if $\sigma$ is sufficiently large,
we conclude $L_a\Phi\le0$,
which proves the desired result.

The previous computations allow us to prove
the main inequality of this section.

We first
consider the case $s=0$. We have that

\begin{eqnarray*}
&& \varphi -
2\overline\eta_R^2\,(\partial_e^2 u)_+^2-
2\tau\,R^{-2}\,{\eta}_R^2\,(\partial_e u)^2
-2\sigma \,R^{-4}\,(u-\kappa)^2\\
-\overline\eta_R^2\,( \partial_e^2u)^2_+
- \tau\,R^{-2}\,\eta_R^2 (\partial_e u )^2-
\sigma\,R^{-4}\,(u-\kappa)^2\le0,
\end{eqnarray*}

which establishes the desired result in this case.

Now, we focus on the case $s\in(0,1)$. From this,
the case $s=1$ can be obtained in the limit[It is interesting to point out that while
we can
obtain (<ref>) with $s=1$ by a limit argument as $s\nearrow1$,
when $s=0$ we needed to perform a direct —though simple— computation
and we could not argue by sending $s\searrow0$. Indeed, the limit as $s\searrow0$ of the fractional Laplacian
involves a term of the type $s\int_{\R^n\setminus B_1}u(y)\,|y|^{-n-2s}\,dy$
which goes to zero if $u$ is compactly supported or decreases
sufficiently fast at infinity, but not in general.
In particular, it is not always true that $ \lim_{s\searrow0}(-\Delta)^su(x)=u(x)$,
since the left-hand side is invariant if one replaces $u$ by $u-\kappa$,
for $\kappa\in\R$, while the right-hand side is not;
see Proposition 4.4 in <cit.>
for additional details.]
(also, one could proceed
direct computations, as in <cit.>).
For $s\in(0,1)$,
we let $U$ be the $s$-harmonic extension of $u$ in $\R^{n+1}_+$,
and we recall that $a= 1-2s$. By <cit.>, up to a positive multiplicative constant depending
only on $s$ —that we do not write since it plays no role to establish (<ref>)— we have

\begin{equation}\label{A1.13p:ZZ}
-\lim_{y\searrow0} y^{a} \partial_y \partial_e U
= (-\Delta)^{s} \partial_e u
\end{equation}

\begin{equation}\label{B1.13p:ZZ}
-\lim_{y\searrow0} y^{a} \partial_y \partial^2_{e} U=
(-\Delta)^{s} \partial_{e}^2 u.\end{equation}

Moreover, if $\Phi^*$ is
the $s$-harmonic extension of $\varphi$ in $\R^{n+1}_+$,
we also have that

\begin{equation}\label{h67io:ZZ}
-\lim_{y\searrow0} y^a \partial_y \Phi^*
= (-\Delta)^{s} \varphi\end{equation}

(recall that $\varphi$ is locally $C^{1,1}$ and bounded in $\R^n$).

The key point is to consider the function

$$ \Phi(x,y):=
\overline\eta_R^2(x)\,\big(\partial_e^2 U(x,y)\big)^2_+
+ \tau\,R^{-2}\,\eta_R^2(x)\big(\partial_e U(x,y)\big)^2+
\sigma\,R^{-4}\,\big(U(x,y)-\kappa\big)^2.$$

Now, by Corollary <ref>, we know that $
L_{a} \Phi\le0$ almost everywhere in $\R^{n+1}_+$,
as long as $\sigma$ and $\tau$ are taken
as in the statement
of Proposition <ref>.
Thus, since $L_{a}\Phi^*=0$,

\begin{equation*}
L_{a} (\Phi-\Phi^*)\le0\quad{\mbox{ everywhere in }}\R^{n+1}_+.\end{equation*}

Notice also that $\Phi-\Phi^*$
vanishes in $\R^n \times\{0\}$,
and hence the maximum principle in Lemma <ref>
gives that $\Phi-\Phi^*\le0$ in $\R^{n+1}_+$.
As a consequence, since this function vanishes in $\R^n\times\{0\}$,
we find that $\lim_{y\searrow0} y^{a}\, \partial_y (\Phi-\Phi^*)\le0$.
Therefore, by (<ref>),
the definition of $\Phi$, and recalling again that $(\partial^2_e U)_+^2$ is a $C^{1,1}_{\rm loc}$ function,

\begin{equation*}
\begin{split}
(-\Delta)^{s} \varphi \;&\le\;-
\lim_{y\searrow0} y^{a} \partial_y \Phi\\
2y^{a}\,\overline\eta_R^2\,(\partial_e^2 U)_+\,
\partial_y\partial_e^2 U
+ 2y^{a}\,\tau\,R^{-2}\,\eta_R^2\,\partial_e U\,\partial_y\partial_e U
\\&\qquad\qquad+2y^{a}\,\sigma\,R^{-4}\,(U-\kappa)\,\partial_y
\end{split}\end{equation*}

We conclude, since $\partial_e U(x,0)=\partial_e u(x)$ and $\partial_e^2U(x,0)=
\partial_e^2u(x)$, that

\begin{equation}\label{g68:BIS:ZZ}\begin{split}
(-\Delta)^{s} \varphi\; \le\;&-
\lim_{y\searrow0}\Big(
2y^{a}\,\overline\eta_R^2\,(\partial_e^2 u)_+\,
\partial_y \partial_e^2 U
+ 2y^{a}\,\tau\,R^{-2}\,{\eta}_R^2\,\partial_e u\,\partial_y \partial_e U
\\&\qquad+2y^{a}\,\sigma\,R^{-4}\,(u-\kappa)\,\partial_y (U-\kappa)\Big)
\\ =\;&2\overline\eta_R^2\,(\partial_e^2 u)_+
(-\Delta)^s \partial_e^2 u
+ 2\tau\,R^{-2}\,\eta_R^2\,\partial_e u\,(-\Delta)^s\partial_e u
\\&\qquad
+2\sigma\,R^{-4}\,(u-\kappa)\,(-\Delta)^s {(u-\kappa)},\end{split}\end{equation}

where (<ref>) and (<ref>)
were used in the last step.

The following result is a simple, but
useful, improvement of Proposition <ref>,
in which we obtain that $\varphi$ is a suitable subsolution
with respect to the linearized operator (<ref>)
(in our current setting, the linear operators $
in (<ref>) boil down to the linear operators $

Let $\eta_R$ and $\overline\eta_R$
be as in (<ref>)-(<ref>).
Let $\kappa\in\R$ and $u\in C^\infty(\R^n)\cap W^{2,\infty}(\R^n)$.
For every $x\in\R^n$, let

\begin{equation}\label{9203-139pq} \varphi(x) :=\overline\eta_R^2(x)\,\big(
\partial_e^2u(x)\big)^2_+
+ \tau\,R^{-2}\,\eta_R^2(x)\big(\partial_e u(x)\big)^2+
\sigma\,R^{-4}\,\big(u(x)-\kappa\big)^2.\end{equation}

Then, there exist positive constants $\sigma_0$ and $\tau_0$, depending only on $n$,
and $\|\overline{{{\eta}}}\|_{C^2(\R^n)}$, such that,
if $\tau\ge \tau_0$ and $\sigma\ge \sigma_0 \tau$, then
in all of $\R^n$ we have

\begin{equation}\label{9203-139pq-BISB} L\varphi \le
2\overline\eta_R^2\,(\partial_e^2 u)_+
\,
L\,\partial_e^2 u
+ 2\tau\,R^{-2}\,{\eta}_R^2\,\partial_e u\,
L\,\partial_e u+2\sigma\,R^{-4}\,(u-\kappa)\,L{(u-\kappa)}.\end{equation}

We write (<ref>)
for every $s\in[0,1]$,
and we integrate with respect to the measure $\mu_j$,
for every given $j\in\{1,\dots,J\}$.
We find that

\begin{equation*}
{\mathcal{L}}_{\mu_j} \varphi\le
2\overline\eta_R^2\,(\partial_e^2 u)_+
\,{\mathcal{L}}_{\mu_j}(\partial_e^2 u)
+ 2\tau\,R^{-2}\,\eta_R^2\,\partial_e u\,{\mathcal{L}}_{\mu_j}
\partial_e u
+2\sigma\,R^{-4}\,(u-\kappa)\,{\mathcal{L}}_{\mu_j} {(u-\kappa)}

Now, we
multiply the above expression by $\alpha_j$, as defined
in (<ref>), using (<ref>),
and then we sum the inequality over $j\in\{1,\dots,J\}$. In this way,
recalling also the definition of $L$ in (<ref>), we conclude the proof.

We remark that Proposition <ref>
is the simplified version of Proposition <ref>
which already leads to first derivative estimates. On the other hand,
Proposition <ref> is stated for a more general cutoff function than the one in
Proposition <ref>. The proof of Proposition <ref> follows the same lines
as that of Proposition <ref>, and only requires Lemma <ref>.

Though not explicitly used in this article, we next state the simplest inequality for an auxiliary function
which leads to one-sided derivative estimates (at least for global solutions). For more details see
the comments before and after (<ref>), in Open problem <ref>,
and before Theorem <ref>. Its proof also follows the same lines as that
of Proposition <ref>.

Let $s\in(0,1)$,
$v\in C^\infty(\R^n)\cap W^{1,\infty}(\R^n)$, $\eta\in C^\infty(\R^n)\cap W^{2,\infty}(\R^n)$,
$e\in \R^n$ with $|e|=1$, and $\sigma>0$.

$$ \psi:=\eta^2 \big(
\partial_e v\big)^2_{+} +\sigma v^2.$$

Then, in all of $\R^n$ we have

\begin{equation*}
(-\Delta)^s \psi\le 2\eta^2\,(\partial_e v)_{+}\,
(-\Delta)^s\partial_e v
+2\sigma v\,(-\Delta)^s v\quad{\mbox{ if }}\sigma\ge\sigma_0,
\end{equation*}

for some constant $\sigma_0$ depending only on $n$ and $\|\eta\|_{C^2(\R^n)}$ —and,
in particular, independent of $s$.

§.§ Proofs of Theorem <ref>, Theorem <ref>, and
Corollary <ref>

We can now prove these results using Theorem <ref>, via an appropriate choice
of the linear operators $\{L^{(x)}\}_{x\in B_R}$.

We may assume that $R:=1$ by considering the rescaled
function $\widetilde u(x):=u(Rx)$ in $B_1$ and the same operator $\mathcal{L}_A$,
but now in $B_1$.
Note that,
since ${\mathcal{A}}$
is compact and we assume the continuity hypothesis (<ref>), observing also that ${\mathcal{L}}_Au(x)$
is continuous with respect to $A$ by the dominated convergence theorem,
we have that, given $x\in B_1$,
${\mathcal{M}}_{{\mathcal{A}}} u(x)={\mathcal{L}}_{A_x}u(x)-g_{A_x}(x)$
for some $A_x\in{\mathcal{A}}$.

We start verifying the hypotheses of Theorem <ref>
for the family of operators $L^{(x)}={\mathcal{L}}_{A_x}$.
We first check that
the maximum principle
with estimate, as considered in Definition <ref>,
is satisfied in this case. Indeed,
assume that (<ref>) holds true with $R=1$. By (<ref>)
applied to the kernel $K(z)=c_{n,s}|Az|^{-n-2s}$, we know that
there exists
a function $\beta$ with $|\beta|\le100$ in $\R^n$
and ${\mathcal{L}}_A\beta\le-c$
in $B_1$
for every $A\in{\mathcal{A}}$.
Here $c>0$ depends only on $n$, $\lambda$, and $\Lambda$.

Define $\widetilde\varphi:=\varphi+\gamma_0\beta/c$.
By (<ref>), we have that $
\inf_{y\in B_1}\{ {\mathcal{L}}_{A_y}\widetilde\varphi(x)\}\le 0$
for every $x\in B_1$.
From this (and again the fact that the previous infimum will be
equal to ${\mathcal{L}}_{A^x}\widetilde\varphi(x)$ for some $A^x\in{\mathcal{A}}$ —$A^x$
perhaps different from $A_x$), one sees that $\sup_{\R^n}\widetilde\varphi$
cannot be achieved in $B_1$ (except when $\widetilde\varphi$ is constant),
and therefore in any case

\begin{equation*}
\max_{\overline{B}_1}\widetilde\varphi\le \sup_{\R^n\setminus B_1}\varphi+\frac{100\,
\gamma_0}c.
\end{equation*}

Consequently, we find that

$$ \sup_{B_1}\varphi\le \sup_{\R^n\setminus B_1}\varphi+\frac{200\,
\gamma_0}c,$$

and hence the maximum principle with estimate, as considered in Definition <ref>,
is satisfied.

Now we check that (<ref>) holds true in the current situation
if we take $\tau$ and $\sigma$ appropriately.
by the definition (<ref>) of $\mathcal{L}_A$, we know that

\begin{equation}\label{TRANSFO} {\mathcal{L}}_A u(x)=\frac1{\det A}(-\Delta)^s
where~$u_A( \overline{x}):=u(A^{-1}\overline{x})$.}}\end{equation}

Also $\partial_{e_A} u_A(\overline{x})=\partial_{e} u(A^{-1}\overline{x})$
and $\partial_{e_A}^2 u_A(\overline{x})=\partial_{e}^2 u(A^{-1}\overline{x})$, with $e_A:=Ae$.
That is, setting $e_A':=e_A/|e_A|$, we have that $ \partial_{e} u(A^{-1}\overline{x})=
|e_A|\,\partial_{e_A'} u_A(\overline{x})$ and $
\partial_{e}^2 u(A^{-1}\overline{x})=|e_A|^2\,
\partial_{e_A'}^2 u_A(\overline{x})$.
Thus, we have that

\begin{equation}\label{TRANS-90sccjer}
(-\Delta)^s \partial_{e'_A} u_A(\overline{x})=
\frac{\det A}{|e_A|}\,{\mathcal{L}}_A\partial_{e} u(x)\quad{\mbox{
and }}\quad
(-\Delta)^s \partial_{e'_A}^2u_A(\overline{x})=\frac{\det A}{|e_A|^2}\,{\mathcal{L}}_A\partial_{e}^2 u(x).\end{equation}

At the same time, if we let $\varphi$ be as in (<ref>) (dropping the index $R$,
since here $R=1$), $\sigma_A:=|e_A|^{-4}\sigma$,
$\tau_A:=|e_A|^{-2}\tau$, $\eta_A(\overline{x}):=\eta(A^{-1}\overline{x})$,
\overline\eta(A^{-1}\overline{x})$,

$$ \overline\varphi_A(\overline{x}):=
\overline\eta^2_A(\overline{x})\,\big(\partial_{e'_A}^2u_A(\overline{x})\big)^2_+
+\tau_A \,\eta^2_A(\overline{x})
\big(\partial_{e'_A} u_A(\overline{x})\big)^2+\sigma_A\,
\left( u_A(\overline{x})-\sup_{A(B_1)} u_A\right)^2,$$

we find that, if $x\in B_1$, and thus $\overline{x}=Ax\in A(B_1)$,

\begin{align*}
=\frac{|e_A|^4}{\det A}\,(-\Delta)^s\overline\varphi_A(\overline{x}).
\end{align*}

As a result,
in view of[Strictly speaking, we are applying here a small variation
of Proposition <ref>,
since the condition (<ref>)
on the cutoff functions reads here $
\eta_A$, $\overline{\eta}_A\in C^\infty_c (A(B_{1/2}))$ and $\eta_A=1$
in $A( B_{1/4})$. It is easy to verify that
this fact does not alter the conclusion (<ref>).
Notice also that $\tau_A$ and $\sigma_A$ are comparable
to $\tau$ and $\sigma$, respectively, with constants depending
only on $n$, $\lambda$, and $\Lambda$.]
Proposition <ref>,
if $\tau$ and $\sigma$ are sufficiently large, we conclude that,
for every $x\in B_1$
(and thus $\overline{x}=Ax\in A(B_1)$),

\begin{align}\label{0o9o08ddf9}
\frac{|e_A|^4}{\det A}\,\Bigg\{
\big(\partial_{e'_A}^2u_A(\overline{x})\big)_+\,
(-\Delta)^s \partial_{e'_A}^2u_A(\overline{x})
\\&\qquad+ 2\tau_A\,\eta^2_A(\overline{x})\,\partial_{e'_A} u_A
(\overline{x})\,(-\Delta)^s\partial_{e'_A} u_A(\overline{x})\\&\qquad
+2\sigma_A \,
\left( u_A(\overline{x})-\sup_{A(B_1)} u_A\right)
\,(-\Delta)^s \left( u_A-\sup_{A(B_1)} u_A\right)(\overline{x})
\Bigg\}.
\end{align}

this and (<ref>), we get

\begin{align*}
\Bigg\{
\big(\partial_{e}^2u(x)\big)_+\,
{\mathcal{L}}_A \partial_{e}^2u(x)
+  2\tau\,\eta^2_A(\overline{x})\,
\partial_{e} u(x)\,{\mathcal{L}}_A\,\partial_{e} u(x)\\&\qquad
+ 2\sigma \,
\left( u(x)-\sup_{B_1} u\right)
\,{\mathcal{L}}_A \left( u-\sup_{B_1} u\right)(x)
\Bigg\}
\end{align*}

for all $x\in B_1$.
Hence, inequality (<ref>) holds true with $
L^{(x)}={\mathcal{L}}_{A_x}$, after choosing, given $x$, $A=A_x$.

With this, we are in the position of applying Theorem <ref>. To this end, we estimate the quantities $a_0$,
$a_1$, and $a_2$ given in its statement. The arguments here are similar to those of Subsection <ref>.
First, by equation (<ref>),

\begin{equation}\label{3.37BIS}
{\mathcal{L}}_{A_x} u(x)=f(x)+g_{A_x}(x)
\end{equation}

for every $x\in B_1$.
Also, using again equation (<ref>),
for every $x$, $y\in B_1$ we have that

\begin{equation}\label{X560199339999:1}\begin{split}&
{\mathcal{M}}_{{\mathcal{A}}} u(y)=
f(y)=f(x)+(f(y)-f(x))\\&\qquad={\mathcal{M}}_{{\mathcal{A}}} u(x) +f(y)-f(x)
\end{equation}

Therefore, if $e\in\R^n$ and $|e|=1$, given $h\in(0, 1-|x|)$,
we see that

\begin{eqnarray*}
{\mathcal{L}}_{A_x} \left(\frac{u(x\pm he)-u(x)}{h}\right)&=&
\frac{{\mathcal{L}}_{A_x}u(x\pm he)-{\mathcal{L}}_{A_x}u(x)}{h}\\
&\le& \frac{f(x\pm he)-f(x)}{h}+\frac{
g_{A_x}(x\pm he)-g_{A_x}(x)
\end{eqnarray*}

and, as a consequence, $
\pm{\mathcal{L}}_{A_x} \partial_e u(x)\le\pm
\big(\partial_e (f+g_{A_x})(x)
\big)$
for almost every $x\in B_1$.
Thanks to the possible sign choice in this inequality,
and to the fact that $f$ and $g_A$ are Lipschitz functions, and thus pointwise differentiable almost everywhere,
we thereby deduce that

\begin{equation}\label{X560199339999:2}
{\mathcal{L}}_{A_x} \partial_e u(x)=
\partial_e (f+g_{A_x})(x)
\end{equation}

for almost every $x\in B_1$.
In addition, exploiting (<ref>) once again,

\begin{equation*} \begin{split}&
{\mathcal{L}}_{A_x} \left(\frac{u(x+ he)+u(x- he)-2u(x)}{h^2}\right)\\&\qquad=
\frac{{\mathcal{L}}_{A_x}u(x+ he)+{\mathcal{L}}_{A_x}u(x- he)-

Thus, if $f+g_{A_x}$ is semiconcave, and hence pointwise twice
differentiable almost everywhere, we have that,
for almost every $x\in B_1$,

\begin{equation}\label{DESDF23} {\mathcal{L}}_{A_x}\partial_e^2 u (x)\le\partial_e^2 (f+ g_{A_x})(x)

In the notation of Theorem <ref>
(with $L^{(x)}$
there corresponding to ${\mathcal{L}}_{A_x}$ here),
(<ref>), (<ref>), and (<ref>)
lead to

\begin{equation}\label{IN939ik3me8834yg}
a_0=\sup_{x\in B_1}\big(f(x)+ g_{A_x}(x)\big)_-
\le \sup_{{x\in B_1}\atop{A\in{\mathcal{A}}}} \big(f(x)+
\big)_-
=\| (f+g_A )_-\|_{L^\infty({\mathcal{A}}\times B_1)}

\begin{equation}\label{IN939ik3me8834yg2}
\sup_{x\in B_1}\big|\partial_e (f+ g_{A_x})(x)\big|
\le\sup_{{x\in B_1}\atop{A\in{\mathcal{A}}}}\big|
\partial_e(f+ g_{A})(x)
\big|=
\|\partial_e(f+ g_{A})  \|_{L^\infty({\mathcal{A}}\times B_1)},
\end{equation}

\begin{equation}\label{IN939ik3me8834yg3}
a_2\le\sup_{x\in B_1}\big(\partial_e^2(f+ g_{A_x})(x)
\big)_+
\le\sup_{{x\in B_1}\atop{A\in{\mathcal{A}}}}\big(
\partial_e^2(f+ g_{A})(x)\big)_+=
\| (\partial_e^2(f+g_{A}) )_+ \|_{L^\infty({\mathcal{A}}\times B_1)}
\,.
\end{equation}

These considerations and Theorem <ref> lead to

\begin{equation}\label{NBCshnr:1}
\sup_{B_{1/2}} |\partial_e u|
\le C\,\Big( \|\partial_e(f+ g_{A}) \|_{L^\infty({\mathcal{A}}\times B_1)}
\| (f+ g_{A} )_-\|_{L^\infty({\mathcal{A}}\times B_1)}
\,\| u\|_{L^\infty(B_1)}
\big)^{1/2}+\| u\|_{L^\infty(\R^n)}
\Big)\end{equation}

\begin{equation}\label{NBCshnr:2}\begin{split}
\sup_{B_{1/4}} \partial^2_e u\le\,& C\,
\Big(
\| (\partial_e^2(f+ g_{A}) )_+ \|_{L^\infty({\mathcal{A}}\times B_1)}
+\|\partial_e(f+g_{A})\|_{L^\infty({\mathcal{A}}\times B_1)}\\&\qquad+\big(
\| (  f+g_{A})_-\|_{L^\infty({\mathcal{A}}\times B_1)}
\,\| u\|_{L^\infty(B_1)}
\big)^{1/2}+\| u\|_{L^\infty(\R^n)}\Big)

for some constant $C$ depending only on $n$, $\lambda$, and $\Lambda$.
After scaling, this gives (<ref>),
and also (<ref>)
with $B_{R/2}$ replaced by $B_{R/4}$. From this result,
it is easy to deduce the second derivative bound
as stated in Theorem <ref> (that is, in $B_{R/2}$)
by a covering argument. To this end,
we simply use the result that we have obtained
in a family of balls $\{B_i\}$ covering $B_{R/2}$,
with radii comparable with $R$, and such that $4B_i\subset B_R$.

We use again
Theorem <ref>, but now
taking, for all $x\in B_R$, $L^{(x)}=L$ to be
the linearized operator
introduced in (<ref>). Let us check that
the hypotheses of Theorem <ref> are fulfilled.
First of all, by
Proposition <ref>,
if $\varphi\in C_{\rm b}(\R^n)\cap {{C}^{1,1}(B_R)}$ is nonnegative
and satisfies $L\varphi\le\gamma_0$ in $B_R$, for some $\gamma_0\ge0$,
then $\sup_{{B_R}}\varphi \le
\sup_{\R^n\setminus{B_R}}\varphi +C_R\,\gamma_0$,

\begin{equation}\label{CIEweRRAsdfJ}
C_R:=\hat C\,(1+R^2),\end{equation}
where $\hat C$
depends only on $n$ and $\vartheta_0$.
This says that $L$ satisfies the maximum principle in $B_R$ with constant $C_R$, according to Definition <ref>.
Moreover, choosing $
\kappa:=\sup_{B_R}u$
in Corollary <ref>, we deduce that (<ref>) holds for $\sigma$ and $\tau$
large enough depending only[Here we use that the proof
of Theorem <ref> requires only two choices of pairs $(\overline\eta_R,\eta_R)$
for cutoff functions.]
on $n$.

Therefore, from
Theorem <ref> we obtain

\begin{equation}\label{7hb2qujjtui501-1}
\sup_{B_{R/2}} |\partial_e u|\le
C\,\Big(C_Ra_1+R^{-1}\big(C_R a_0\|u\|_{L^\infty(B_R)}\big)^{1/2}+
\Big)\end{equation}

\begin{equation}\label{7hb2qujjtui501-2}
\sup_{B_{R/4}} \partial^2_e u\le C\,\Big(
C_R a_2+C_R R^{-1} a_1
+R^{-2}\big(C_R a_0 \| u\|_{L^\infty(B_R)}\big)^{1/2}+R^{-2}\| u\|_{L^\infty(\R^n)}
\Big),\end{equation}

where $C$ is a constant depending only on $n$.

The quantities $a_0$, $a_1$, and $a_2$ in the statement of Theorem <ref> are given by

\begin{equation*}
a_0:=\sup_{x\in B_R} \Big(L \big(u-\sup_{B_R}u\big)( x)\Big)_- \, ,\quad
a_1:=\sup_{x\in B_R} \big|L\,\partial_e u( x)\big|,\quad
a_2:=\sup_{x\in B_R}\big(
L\,\partial_e^2 u( x)\big)_+\,.
\end{equation*}

Now, by
Lemma <ref>, recalling the definition
of $\omega_0$ given before the statement of
Theorem <ref>,
and taking into account that $L$ could contain
the operator $(-\Delta)^0={\rm Id}$,
in $B_R$ we have that

\begin{eqnarray*}&&
\Big(L\big(u-\sup_{B_R}u\big)\Big)_-=\max\Big\{-Lu+L\sup_{B_R}u,0\Big\}
\le\max\{ -f-g+\omega_0\|u\|_{L^\infty(B_R)} ,0\}\\&&\qquad\qquad\le
\| (f+g)_-\|_{L^\infty(B_R)}+\omega_0\| u\|_{L^\infty(B_R)},
\end{eqnarray*}

$\big| L\,\partial_{e} u \big| \le
\gamma_1$,
and $
\big(L\,\partial_e^2 u\big)_+\le \gamma_2$.
Here, we are using the notation introduced before Theorem <ref>
for $\omega_0$ (that is, $\omega_0=0$ if the operator does not involve the identity and $\omega_0=1$ otherwise).
Hence, we conclude that

$$ a_0\le\,\|(f+g)_- \|_{L^\infty(B_R)}+\omega_0\| u\|_{L^\infty(B_R)},\qquad
a_1\le\gamma_1,\quad{\mbox{ and }}\quad

This and the first derivative bound (<ref>) lead to

\begin{equation}\label{7hb2qujjtui501-1BOS}\begin{split}
\sup_{B_{R/2}} |\partial_e u|\,\le\,&
C\,\Big(C_R\,\gamma_1+R^{-1}\big(C_R \,(\| (f+g)_-\|_{L^\infty(B_R)}+\omega_0\|u\|_{L^\infty(B_R)})
\|u\|_{L^\infty(B_R)}\big)^{1/2}\\&\qquad+
\Big).\end{split}\end{equation}

Similarly, in light of the second derivative bound (<ref>),
we see that

\begin{equation}\label{7hb2qujjtui501-2BOS}\begin{split}
\sup_{B_{R/4}} \partial^2_e u\le\,& C\,\Big(
C_R \gamma_2+C_R R^{-1} \gamma_1
\\&\qquad+R^{-2}\big(C_R \, (\| (f+g)_-\|_{L^\infty(B_R)}+\omega_0\|u\|_{L^\infty(B_R)})
\| u\|_{L^\infty(B_R)}\big)^{1/2}+R^{-2}\| u\|_{L^\infty(\R^n)}
\Big).\end{split}\end{equation}

These observations complete
the proof of
Theorem <ref>, after
using a covering argument as at the end of the proof of
Theorem <ref>.

Corollary <ref> follows immediately from Theorem <ref>, as we show
next. Alternatively, the corollary could also be established
following the lines of the proof of Theorem <ref> —hence, using the $\max$
structure of the operator instead of introducing the linearized operator $L$
as in Theorem <ref>.

We are in the setting of Theorem <ref>, with $J=2$,
$\mathcal{L}_{\mu_1}=(-\Delta)^s$, $g_1=0$,
$\mathcal{L}_{\mu_2}=(-\Delta)^0$ (the identity), $g_2=g$, and $F$
being the $\max$ operator —which satisfies (<ref>)-(<ref>)
as pointed out
right after (<ref>). In particular, we have $\vartheta_0=1$ and $|\alpha_j|\le1$ everywhere.