[step:Bessel potential of $\mu$ lies in $W^{1, p}$ via finite Wolff and Riesz energy]Let $G_1$ denote the Bessel kernel of order $1$, the function on $\mathbb{R}^n$ whose Fourier transform is $(1 + 4\pi^2 |\xi|^2)^{-1/2}$. The Bessel-potential space $\{G_1 * f : f \in L^p(\mathbb{R}^n)\}$ coincides with the Sobolev space $W^{1, p}(\mathbb{R}^n)$ for $1 < p < \infty$ as Banach spaces, with equivalent norms (Adams-Hedberg, *Function Spaces and Potential Theory*, §1.2.4 and Thm 1.2.3); we will denote this Bessel-potential space by $\mathcal{B}^{1, p}(\mathbb{R}^n)$ to avoid confusion with Lorentz spaces. We define the Bessel potential of $\mu$ by
\begin{align*}
u: \mathbb{R}^n \to [0, \infty], \qquad u(x) = (G_1 * \mu)(x) = \int_{\mathbb{R}^n} G_1(x - y)\, d\mu(y).
\end{align*}
The kernel $G_1$ satisfies the standard pointwise bound (Adams-Hedberg §1.2.5; Maz'ya §11.3.1)
\begin{align*}
G_1(x) \le C(n) \cdot |x|^{-(n - 1)} \qquad \text{for } |x| \le 1,
\end{align*}
with exponential decay at infinity, so $G_1$ is comparable to the Riesz kernel $|x|^{-(n - 1)}$ on bounded regions.
The key analytic fact is the Hedberg-Wolff theorem (Hedberg and Wolff, *Thin sets in nonlinear potential theory*, Ann. Inst. Fourier 33 (1983); see Adams-Hedberg, Thm 4.5.4 and Cor 4.5.5): for $1 < p < n$ and any non-negative Radon measure $\mu$ of compact support,
\begin{align*}
\|G_1 * \mu\|_{W^{1, p}(\mathbb{R}^n)}^{p} \asymp \mathcal{W}_{1, p}(\mu) := \int_{\mathbb{R}^n} \int_0^\infty \left( \frac{\mu(B(x, r))}{r^{n - p}} \right)^{p'-1}\, \frac{dr}{r}\, d\mu(x),
\end{align*}
where $p' = p/(p - 1)$ is the Hölder conjugate and $\mathcal{W}_{1, p}(\mu)$ is the Wolff $(1, p)$-energy of $\mu$. (For $p = 2$ the Wolff energy reduces, after a Fubini computation, to the linear Riesz $(n-2)$-energy; for $p \ne 2$ it is genuinely nonlinear.)
The Riesz $(n-p)$-energy $I_{n-p}(\mu) := \int\int |x - y|^{-(n-p)}\, d\mu(y)\, d\mu(x)$ provides a sufficient linear surrogate: a standard dyadic computation (Adams-Hedberg, Prop 4.5.6) shows that the Frostman growth $\mu(B(x, r)) \le r^s$ with $s > n - p$ forces both $\mathcal{W}_{1, p}(\mu) < \infty$ and $I_{n - p}(\mu) < \infty$, and that
\begin{align*}
\mathcal{W}_{1, p}(\mu) \le C(n, p, s, D)\, I_{n-p}(\mu) \qquad \text{when } \mu(B(x, r)) \le r^s,\ s > n - p.
\end{align*}
Hence it suffices to bound $I_{n-p}(\mu)$. We claim the Frostman bound forces $I_{n - p}(\mu) < \infty$.
*Riesz energy from Frostman.* Fix $x \in \operatorname{supp}\mu$ and let $D := \operatorname{diam}(K) < \infty$. Decomposing the inner integral dyadically with $A_k := \{y : 2^{-k - 1} D < |x - y| \le 2^{-k} D\}$,
\begin{align*}
\int_{\mathbb{R}^n} \frac{d\mu(y)}{|x - y|^{n - p}} \le \sum_{k = 0}^\infty (2^{-k - 1} D)^{-(n - p)}\, \mu(B(x, 2^{-k} D)) \le 2^{n - p} D^{-(n - p)} \sum_{k = 0}^\infty 2^{k(n - p)}\, (2^{-k} D)^s,
\end{align*}
where the last step uses $\mu(B(x, r)) \le r^s$. Simplifying,
\begin{align*}
\int_{\mathbb{R}^n} \frac{d\mu(y)}{|x - y|^{n - p}} \le C(n, p)\, D^{s - (n - p)} \sum_{k = 0}^\infty 2^{-k(s - (n - p))} = \frac{C(n, p)\, D^{s - (n - p)}}{1 - 2^{-(s - (n - p))}},
\end{align*}
which is finite because $s - (n - p) > 0$. Integrating in $\mu$,
\begin{align*}
I_{n - p}(\mu) \le \frac{C(n, p)\, D^{s - (n - p)}}{1 - 2^{-(s - (n - p))}} \cdot M.
\end{align*}
Combining with $\mathcal{W}_{1, p}(\mu) \le C(n, p, s, D)\, I_{n-p}(\mu)$ and the Hedberg-Wolff identity $\|G_1 * \mu\|_{W^{1, p}}^p \asymp \mathcal{W}_{1, p}(\mu)$, we conclude $u = G_1 * \mu \in W^{1, p}(\mathbb{R}^n)$ with $\|u\|_{W^{1, p}}^p \le C(n, p, s, D)\, M$.[/step]