Existence and regularity for fractional Poisson-type equation According to Theorem 2.7 in the paper https://arxiv.org/pdf/1704.07560.pdf, we have the following classical results.
Let $s \in (0,1)$ and $1<p<\infty$. Then for any $F \in L^p(\mathbb{R}^n)$ there exists a unique weak solution $u$ of the equation $(-\Delta)^s u=F$ in $\mathbb{R}^n$, where $$(-\Delta)^s(x) = c_{n,s} P.V.\int_{\mathbb{R}^n} \frac{u(x)-u(y)}{|x-y|^{n+2s}}dy$$ denotes the fractional Laplacian. Moreover, $u$ belongs to the Bessel potential space $H^{2s,p}(\mathbb{R}^n) = \{u \in L^p(\mathbb{R}^n) \mid (-\Delta)^s u \in L^p(\mathbb{R}^n) \}$.
I have two questions regarding this probably well-known result.
First of all, I am unsure what exactly is meant by a weak solution in this context. Typically, in the case when $F \in L^2(\mathbb{R}^n)$ we say that $u \in H^s(\mathbb{R}^n)$ is a weak solution of the equation $(-\Delta)^s u=F$ in $\mathbb{R}^n$, if for any $\varphi \in H^s(\mathbb{R}^n)$ we have
$$ \frac{c_{n,s}}{2} \int_{\mathbb{R}^n} \int_{\mathbb{R}^n} \frac{(u(x)-u(y))(\varphi(x)-\varphi(y))}{|x-y|^{n+2s}}dydx = \int_{\mathbb{R}^n} F \varphi dx,$$
however I doubt that this notion of weak solutions is appropriate when $F \in L^p(\mathbb{R}^n)$ for $p \neq 2$ and the equation is posed on the whole space $\mathbb{R}^n$. Moreover, since the above paper contains no proof of the above result, I would be grateful for a brief outline of the proof of this result or a reference containing it.
 A: [This is an extension of my comment, and it does not really answer the original question].
The question includes the following statement: there is a unique solution to $(-\Delta)^s u = F$ for any $f \in L^p(\mathbb{R}^n)$. This is not quite correct: the solution need not exist, and if it exists, then it is unique.

Uniqueness is a simpler matter: it is straightforward to see that $(-\Delta)^s f(x) = 0$ for all $x \in \mathbb{R}^n$ if $f$ is a constant function. If $s > \tfrac{1}{2}$, then the same is true when $f$ is an affine function, that is, $f(x) = b \cdot x + c$ for some $b \in \mathbb{R}^n$ and $c \in \mathbb{R}$. It follows that the solution of $(-\Delta)^s u = F$ is never unique, unless we impose additional conditions on the behaviour of $u$ at infinity.
Interestingly, every $(-\Delta)^s$-harmonic function is an affine function. I believe this was part of the folklore until it has been rigorously proved by Mouhamed Moustapha Fall in [M. M. Fall, Entire $s$-harmonic functions are affine. Proc. Amer. Math. Soc. 144(6) (2016): 2587–2592] (link).

Existence is somewhat more complicated. If $p < \tfrac{n}{2 s}$, then $(-\Delta)^s$ is a one-to-one unbounded operator on $L^p(\mathbb{R}^n)$, and its inverse is the Riesz potential operator
$$ \mathscr{I}_{2s} u(x) = c_{n,-s} \int_{\mathbb{R}^n} \frac{u(y)}{|x - y|^{n - 2 s}} \, dy . $$
Furthermore, the Riesz potential operator is well-defined on all of $L^p(\mathbb{R}^n)$, and it defines a solution: for every $F \in L^p(\mathbb{R}^n)$ we have $(-\Delta)^s \mathscr{I}_{2s} F = F$. An excellent read on this subject is Stefan Samko's book [S. Samko, Hypersingular Integrals and Their Applications. CRC Press, London-New York (2001)]. The two results mentioned above, and a dozen of closely related statements, are given in Chapter 7 of this book.
A historical remark: The Riesz potential operator was introduced by Marcel Riesz in [M. Riesz, Intégrales de Riemann–Liouville et potentiels. Acta Sci. Math. Szeged 9 (1938), 1–42] (link), while $(-\Delta)^s$ and the relation between the two operators was first studied in [M. Riesz, L’intégrale de Riemann–Liouville et le problème de Cauchy. Acta Math. 81 (1949). 1–223] (link).

On the other hand, if $p \ge \tfrac{n}{2 s}$, then there are functions $F \in L^p(\mathbb{R}^n)$ such that there is no solution of $(-\Delta)^s u = F$. I do not have access to Samko's book at the moment, so I cannot see if that fact is given there; I will update this answer as soon as I check this.
There are at least two ways to convince oneself that this is true. For simplicity, let us consider $p = \infty$ and $F$ constant. Then the Fourier transform of $F$ is the Dirac delta at zero. On the other hand, the Fourier transform of $(-\Delta)^s u$ is $|\xi|^{2 s} \hat u(\xi)$, and there is no tempered distribution $u$ such that the above product is the Dirac delta at zero. This sounds simple, but lots of technicalities are hidden here: one needs to define the product of tempered distributions appropriately, show that the Fourier transform of $(-\Delta)^s u$ is indeed the product of $|\xi|^{2s}$ and $\hat u$ (which is well-known for Schwarz class functions, but far from obvious for tempered distributions), and finally prove that indeed the product of a tempered distribution and $|\xi|^{2 s}$ is never a Dirac delta. Details are partially given (in a different context) in two of my articles: Spectral analysis... and Ten equivalent definitions....
Another approach involves writing $u$ in a ball $B(0, r)$ in terms of the  Poisson integral and the Green function. Taking $r \to \infty$ one can apparently get a contradiction with integrability of $(1 + |x|)^{-n - \alpha} u(x)$; but I never really worked out the details.

Finally, let me stress that although I spent quite some time working on the fractional Laplace operator, I do not really know the PDE literature on that subject well, and it is quite likely that what I attempted to explain above can be found rigorously written in some paper.
