Which Sobolev spaces are function spaces? Let us consider the Sobolev space $W^{s,p}(\mathbb R^n)$, $s \in \mathbb R$ (I am most interested in the case $n=1$). If $s \geq 0$, it embeds to $L^p(\mathbb R^n)$, so its elements may be interpreted as functions (modulo redefinitions on sets of measure zero). If $s$ is sufficiently negative then not all elements of $W^{s,p}(\mathbb R^n)$ are functions, for example the Dirac's delta may belong to $W^{s,p}(\mathbb R^n)$. I would like to ask if there is some range of negative $s$ for which all elements in $W^{s,p}(\mathbb R^n)$ are genuine functions (e.g. from $L^1_{\mathrm{loc}}(\mathbb R^n))$.
 A: The answer to your question is no: for every $s<0$ and $p\in (1,\infty)$ there exists a $\psi\in W^{s,p}$ that cannot be represented by a function in $L^1_{\mathrm{loc}}$.
It is not so easy to provide a complete answer, but I want to stress the following points, which help identifying distributions in $W^{-s,p}$ which are not representable by functions, in a sort of minimal way.

*

*The problem is not one of locality but rather one of differentiability/continuity;


*The case of general $p$ should be addressed by resorting to (Bessel) capacities;


*The case $p=2$ is well exemplified by the theory of Gaussian fields.
By 1., I mean that there is little difference in considering $L^1$ rather than $L^1_{\mathrm{loc}}$. Indeed, the definition of $W^{1,p}$ by the OP (the one in the comments, via Bessel potentials) makes perfect sense (and is one standard definition) on closed manifolds (compact no boundary, with $\Delta=\Delta_g$ the Laplace–Beltrami operator) or on bounded domains (with the Neumann Laplacian). In these cases, the distinction between $L^1$ and $L^1_{\mathrm{loc}}$ is immaterial.
By 2., I mean that – morally – we are looking for distributions represented by measures that do not charge sets of vanishing Bessel $C_{s,p}$-capacity (definition below). In the same sense in which $L^p$ functions are defined up to sets of vanishing measure, $W^{s,p}$ functions are defined up to sets of vanishing $C_{s,p}$-capacity (henceforth: $C_{s,p}$-polar sets). This is important because it allows us to get rid of the representation problem and turns everything into a summability problem. Indeed let $\mathcal M_{s,p}$ be the class of (non-negative not necessarily $\sigma$-finite) measures not charging $C_{s,p}$-polar sets. Then $\mu(f):= \int f\mathrm{d}\mu$ is well-defined for every $\mu\in\mathcal M_{s,p}$ whenever $f\geq 0$ is defined up to $C_{s,p}$-polar sets (in particular, if $f$ is some representative of a function in $W^{s,p}$).
This shows that $f\mapsto \mu f$ is representable as an integral, and therefore we only need to ask when is this integral finite.
$C_{s,p}$-polarity is well-related to Hausdorff dimension, thus we are looking at measures not charging sets of too small Hausdorff dimension. On the other hand, representability by functions is related to the property of not charging measure $0$ sets. There is a large class of measures in between these two extrema. This justifies why, in dimension $1$ and for $p=2$, we look for measures charging (thin) Cantor(-like) sets.
By 3., I mean that if one is content with distributions that are not representable by functions in $L^2_{\mathrm{loc}}$ (or, equivalently, $L^2$), then many examples come from (derivatives of) the Brownian motion. Again morally, this is very simply understood by recalling that $B_t$ is almost surely $C^{0,1/2-}$ (i.e. $(1/2-\epsilon)$-Hölder continuous for all $\epsilon>0$) thus its derivative is a good starting point. The same holds for the Gaussian Free Field in dimension 2, and for log-correlated (fractional) Gaussian fields, all objects that are in $W^{0-,2}$ (and in fact not representable by functions at all, see [B] below).
Bessel kernel representation: Let $p_t(x,y):=(4\pi t)^{-n/2} e^{-|x-y|^2/4t}$ be the $n$-dimensional heat kernel. The $s$-power of the resolvent of $-\Delta$ at $1$ is represented by the Bessel kernel $G_s(x,y)$
$$(1-\Delta)^{-s}f=\int_{\mathbb R^n}\underbrace{\frac{1}{\Gamma(s)}\int_0^\infty e^{-t} t^{s-1} p_t(x,y) dt}_{=:G_s(x,y)} f(y) dy$$
where $\Gamma$ is Euler Gamma function.
(if $f$ is a distribution, then the rhs should be interpreted as a convolution of the integral kernel acting in the variable $y$ with $f$ also acting in the variable $y$, resulting in a distribution acting in the variable $x$)
$G_s(x,y)$ is the fundamental solution to the Poisson equation for $(1-\Delta)^{-s}$, i.e.
$$(1)\qquad\qquad (1-\Delta_y)^{-s} G_s(x,y)=\delta_x$$
Since $(1-\Delta)^{-s}(1-\Delta)^{-r}=(1-\Delta)^{-(s+r)}$, the kernels $(G_s)_{s\geq 0}$ satisfy the semigroup property w.r.t. convolution. For integer $k$, this reads $G_k=G_1^{*k}$, where ${}^{*k}$ is the $k$-th convolution power.
Comments: the definition with the heat kernel is useful to motivate that the discussion below also holds on (closed) manifolds and beyond.
Bessel capacities: (see e.g. [A, §2.6])
Define
$$C_{s,p}(E):= \inf \{\|f\|_p^p: G_s*f\geq 1 \text{ on } E, f\geq 0\}\qquad E\subset \mathbb R^n$$
([A, 2.6.2]).
Then ($\mathscr H^d$ is the $d$-dimensional Hausdorff measure):

a. $C_{s,p}$ is a Choquet capacity ([A, 2.6.8])
b. if $p\in (1,\infty)$ and $sp\leq n$, then $C_{s,p}E=0$ if $\mathscr H^{n-sp}E<\infty$. Conversely, if $C_{s,p}E=0$, then $\mathscr H^{(n-sp)+}E=0$ ([A, 2.6.16])
c. $f\in W^{s,p}$ is defined up to $C_{s,p}$-polar sets (consequence of [A,2.6.4])

About representability in point 2.
$\mu\in \mathcal M_{s,p}\cap \mathcal M_b$ (finite measures) is a continuous linear functional on $W^{s,p}\cap L^\infty(P_{s,p})$, where $L^\infty(P_{s,p})$ is the $L^\infty$-space of the $\sigma$-ideal $P_{s,p}$ of $C_{s,p}$-polar sets.
This is a consequence of b. and c. above.
Summability should be addressed as in the question pointed out in the comments. In particular, the need to restrict to $L^\infty(P_{s,p})$ in the above discussion arises from the fact that we chose $\mu\in\mathcal M_{s,p}\cap \mathcal M_b$ (i.e. the sharp threshold). If we were to choose $\mu\in\mathcal M_{s-,p}\cap \mathcal M_b$, then the summability of $\mu f$ would be a consequence of the $\mu$-volume growth of balls, which is analyzed in [A, 2.6.13].
Comments: We have just simplified the problem from discussing the representability of a distribution (representing the measure $\mu$) to discussing summability of an integral.
This is no minor thing, as one can better understand looking at when $\delta_0\in W^{-s,p}$. We have
$$(2)\qquad\qquad \delta_x\in W^{-s,p'} , \qquad s>n-n/p'$$
Indeed we know that $W^{s,p}\hookrightarrow \mathcal C^0$ by Sobolev–Morrey's embedding for $s>n/p$ or $s=n$ if $p=1$, so that $\delta_x$ can be tested against all functions in $W^{s,p}$ in the above cases, and therefore $\delta_x\in W^{-s,p'}$ whenever $s>n/p=n-n/p'$ (multiply by $n$ the equality $1=1/p+1/p'$) and $\delta_x\in W^{-n,\infty}$.
In this sense, the representability of $\delta_0\in W^{-s,p'}$ is equivalent to the Sobolev–Morrey embedding $W^{s,p}\hookrightarrow \mathcal C^0$.
Comments: The Bessel capacity is equivalent to the Riesz capacity in the question ..., but the Riesz capacity is less well-behaved on bounded domains and on manifolds.
About $q$-representability.
We will need some numerology, so everywhere in the following $s>0$, $p\in (1,\infty)$, and $p'$ is the dual Hölder exponent of $p$.
The proof should work also for the case $p=\infty$, but the definition of $W^{s,\infty}$ via Bessel potential is somewhat less standard.
Say that a distribution is $q$-bad if it is not representable by a function in $L^q_{\mathrm{loc}}$.
Let us show that:

for every $s>0$ and $p\in (1,\infty)$ there exists a distribution in $W^{-s,p'}$ that is $q$-bad for all $q\geq p'$

while taking for granted the following result (which I would like to regard as somewhat standard):

subcritical Sobolev embeddings are sharp, i.e. $W^{s,p}\hookrightarrow\!\!\!\!\!\not \,\,\,\, \mathcal C^0$ for any $s<n/p$ and $p\in(1,\infty)$.

Recall that $W^{s,p}$ is a local space: if $\phi\in\mathcal C^\infty$ and $f\in W^{s,p}$ then $\phi\cdot f\in W^{s,p}$ (also for $s<0$, in which case $\cdot$ denotes the product of a smooth function times a distribution)
When $-s<-n+n/p'$ the conclusion is already shown by (2). (Since $\delta_x$ is obviously $q$-bad for every $q\in [1,\infty]$.)
By (2) and the semigroup property of $G_s$ we have that $\psi_p:=(1-\Delta)^{-\frac{n}{2}+\frac{n}{2p'}}\delta_0\in W^{-s,p'}$ for all $s<0$. It suffices to show that it is a $q$-bad distribution for all $q>p'$. Argue that $\psi_p$ is $q$-representable, i.e. $\psi_p\in L^{q}_{\mathrm{loc}}$.
This is equivalent to say that
$$g\longmapsto \int \psi_p(x) g(x) dx ,\qquad \qquad g\in L^{q'}(B)$$
is a continuous linear functional on $L^{q'}(B)$, with  $B:=B_1(0)\subset \mathbb R^n$. Regarding $L^{q'}(B)\subset L^1$ (extending by $0$ outside of $B$), we have
$$\int \psi_p(x) g(x) dx = \delta_0 \left[(1-\Delta)^{-\frac{n}{2}+\frac{n}{2p'}} g\right]$$
which is a continuous linear functional if and only if $(1-\Delta)^{-\frac{n}{2}+\frac{n}{2p'}} g$ has a continuous representative. Since $g$ is arbitrary, this implies that $W^{n-n/p',q'}(B):=(1-\Delta)^{-\frac{n}{2}+\frac{n}{2p'}} L^{q'}(B)\hookrightarrow \mathcal C^0$.
Since Sobolev spaces are local, $W^{n-n/p',q'}(B) \hookrightarrow \mathcal C^0$ if and only if $W^{n-n/p',q'} \hookrightarrow \mathcal C^0$ (here we are using the extension property, which we can do since $B$ is a nice domain).
Now we use the fact that subcritical Sobolev embeddings are sharp. In order for $W^{n-n/p',q'} \hookrightarrow \mathcal C^0$ to hold we need that either $$n-n/p'>n/q' \qquad \text{or} \qquad n-n/p'=n \,\text{ and }\, q'=1$$ (which are the thresholds for Sobolev–Morrey's embeddings: $s_0>n/p_0$ or $s_0=n$ and $p_0=1$ with $s_0:= n-n/p'$ and $p_0:=q'$).
We have
\begin{align*}
n-n/p'>& n/q' \\
n/p>& n/q' \\
q'>& p\\
q<& p'
\end{align*}
About fractional Gaussian fields:
Let me give an example arising from point 3. in the beginning.
The fractional Gaussian Field $h$ of order $s$ on $\mathbb R^n$ is the $\mathcal D'$-valued centered Gaussian field with covariance
$$\mathbf E[\langle h \mid \phi \rangle \langle h \mid \psi \rangle]=\iint \phi(x) G_s(x,y) \psi(y) \mathrm{d} x \, \mathrm{d} y$$
If $s=n/2$, then $h$ is log-correlated (because $G_s$ has logarithmic divergence $-\ln |x-y|$ as $x\to y$) and satisfies $h\in W^{0-,2}$. In dimension 2 and for $s=1$, $h$ is the Gaussian Free Field (see e.g. [B]).
[A] W. P. Ziemer, Weakly differentiable functions, Sobolev spaces and functions of bounded variation, 1989
[B] N. Berestycki, Introduction to the Gaussian Free Field
and Liouville Quantum Gravity
