**Step 1** [claim] For any fixed point $x \in \mathbb{R}^n$ and radius $r > 0$, the integration of the function differences over the open ball $B(x,r) \subseteq \mathbb{R}^n$ with respect to the $n$-dimensional Lebesgue measure $\Lambda^n$ satisfies the bound: \begin{align*} \frac{1}{\int_{B(x,r)} 1 \, d\Lambda^n(y)} \int_{B(x,r)} |u(x) - u(y)| \, d\Lambda^n(y) \le \frac{1}{n \alpha_n} \int_{B(x,r)} \frac{|\nabla u(z)|}{|x-z|^{n-1}} \, d\Lambda^n(z) \end{align*} where $\alpha_n$ is the volume of the unit ball in $\mathbb{R}^n$. [/claim] [proof] Fix a direction vector $w \in \partial B(0,1)$. Using the Fundamental Theorem of Calculus along the line segment from $x$ to $x+sw$, we compute the difference using the one-dimensional Lebesgue measure $\Lambda^1$: \begin{align*} u(x+sw) - u(x) = \int_0^s \partial_t u(x+tw) \, d\Lambda^1(t) = \int_0^s \nabla u(x+tw) \cdot w \, d\Lambda^1(t) \end{align*} Taking the absolute value and applying the Cauchy-Schwarz inequality with $|w| = 1$, we obtain: \begin{align*} |u(x+sw) - u(x)| \le \int_0^s |\nabla u(x+tw)| \, d\Lambda^1(t) \end{align*} We integrate this inequality over the unit sphere $\partial B(0,1)$ with respect to the Hausdorff surface measure $\sigma$: \begin{align*} \int_{\partial B(0,1)} |u(x+sw) - u(x)| \, d\sigma(w) \le \int_{\partial B(0,1)} \int_0^s |\nabla u(x+tw)| \, d\Lambda^1(t) \, d\sigma(w) \end{align*} We multiply both sides by $s^{n-1}$ and integrate with respect to the variable $s$ over the interval $[0, r]$ using the one-dimensional Lebesgue measure $\Lambda^1$: \begin{align*} \int_0^r \int_{\partial B(0,1)} |u(x+sw) - u(x)| s^{n-1} \, d\sigma(w) \, d\Lambda^1(s) \le \int_0^r s^{n-1} \int_{\partial B(0,1)} \int_0^s |\nabla u(x+tw)| \, d\Lambda^1(t) \, d\sigma(w) \, d\Lambda^1(s) \end{align*} We apply the polar coordinates formula for Lebesgue measure, which states that the $n$-dimensional Lebesgue measure $\Lambda^n$ decomposes geometrically as $d\Lambda^n(y) = s^{n-1} \, d\sigma(w) \, d\Lambda^1(s)$ under the substitution $y = x+sw$. The left-hand side transforms exactly into the volume integral over the ball $B(x,r)$: \begin{align*} \int_0^r \int_{\partial B(0,1)} |u(x+sw) - u(x)| s^{n-1} \, d\sigma(w) \, d\Lambda^1(s) = \int_{B(x,r)} |u(x) - u(y)| \, d\Lambda^n(y) \end{align*} For the right-hand side, since $s \le r$ and the integrand is strictly non-negative, we enlarge the upper limit of the inner integral with respect to $t$ from $s$ to $r$: \begin{align*} \int_0^s |\nabla u(x+tw)| \, d\Lambda^1(t) \le \int_0^r |\nabla u(x+tw)| \, d\Lambda^1(t) \end{align*} This uniform upper limit allows us to factor the integrals via Fubini's Theorem: \begin{align*} \int_0^r s^{n-1} \left( \int_{\partial B(0,1)} \int_0^r |\nabla u(x+tw)| \, d\Lambda^1(t) \, d\sigma(w) \right) \, d\Lambda^1(s) = \left( \int_0^r s^{n-1} \, d\Lambda^1(s) \right) \left( \int_{\partial B(0,1)} \int_0^r |\nabla u(x+tw)| \, d\Lambda^1(t) \, d\sigma(w) \right) \end{align*} Evaluating the first integral gives $\int_0^r s^{n-1} \, d\Lambda^1(s) = \frac{r^n}{n}$. For the second integral, we apply the polar coordinates formula in reverse using the change of variables $z = x+tw$, where the Lebesgue measure $\Lambda^n$ decomposes as $d\Lambda^n(z) = t^{n-1} \, d\sigma(w) \, d\Lambda^1(t)$. We multiply and divide the integrand by $t^{n-1}$: \begin{align*} \int_{\partial B(0,1)} \int_0^r \frac{|\nabla u(x+tw)|}{t^{n-1}} t^{n-1} \, d\Lambda^1(t) \, d\sigma(w) = \int_{B(x,r)} \frac{|\nabla u(z)|}{|x-z|^{n-1}} \, d\Lambda^n(z) \end{align*} Combining these evaluations bounds the integral over $B(x,r)$: \begin{align*} \int_{B(x,r)} |u(x) - u(y)| \, d\Lambda^n(y) \le \frac{r^n}{n} \int_{B(x,r)} \frac{|\nabla u(z)|}{|x-z|^{n-1}} \, d\Lambda^n(z) \end{align*} We divide by the measure of the ball, given by $\int_{B(x,r)} 1 \, d\Lambda^n(y) = \alpha_n r^n$, to obtain the final inequality: \begin{align*} \frac{1}{\int_{B(x,r)} 1 \, d\Lambda^n(y)} \int_{B(x,r)} |u(x) - u(y)| \, d\Lambda^n(y) \le \frac{1}{n \alpha_n} \int_{B(x,r)} \frac{|\nabla u(z)|}{|x-z|^{n-1}} \, d\Lambda^n(z) \end{align*} [/proof]

**Step 2** [claim] There exists a constant $C_1 > 0$ depending only on $n$ and $p$ such that for any $x \in \mathbb{R}^n$ and $r > 0$: \begin{align*} \frac{1}{\int_{B(x,r)} 1 \, d\Lambda^n(y)} \int_{B(x,r)} |u(x) - u(y)| \, d\Lambda^n(y) \le C_1 r^{1 - \frac{n}{p}} \left( \int_{B(x,r)} |\nabla u(z)|^p \, d\Lambda^n(z) \right)^{\frac{1}{p}} \end{align*} [/claim] [proof] We apply Hölder's inequality to the integral on the right-hand side of the previous claim. Let $q = \frac{p}{p-1}$ be the conjugate exponent to $p$. We split the integrand into the product of $|\nabla u(z)|$ and $|x-z|^{1-n}$: \begin{align*} \int_{B(x,r)} \frac{|\nabla u(z)|}{|x-z|^{n-1}} \, d\Lambda^n(z) \le \left( \int_{B(x,r)} |\nabla u(z)|^p \, d\Lambda^n(z) \right)^{\frac{1}{p}} \left( \int_{B(x,r)} |x-z|^{-(n-1)q} \, d\Lambda^n(z) \right)^{\frac{1}{q}} \end{align*} We compute the second integral using the polar coordinates formula centered at $x$, where the radius is denoted $\rho = |x-z|$ and the Lebesgue measure $\Lambda^n$ decomposes into $\rho^{n-1} \, d\sigma(w) \, d\Lambda^1(\rho)$. The total surface area $\int_{\partial B(0,1)} 1 \, d\sigma(w)$ evaluates to $n \alpha_n$: \begin{align*} \int_{B(x,r)} |x-z|^{-(n-1)q} \, d\Lambda^n(z) &= \int_0^r \int_{\partial B(0,1)} \rho^{-(n-1)q} \rho^{n-1} \, d\sigma(w) \, d\Lambda^1(\rho) \\ &= n \alpha_n \int_0^r \rho^{n - 1 - (n-1)q} \, d\Lambda^1(\rho) \end{align*} Substituting $q = \frac{p}{p-1}$, the exponent of $\rho$ evaluates exactly to $\frac{1-n}{p-1}$. Because $p > n$, this exponent strictly exceeds $-1$, ensuring integrability at the origin. We perform the one-dimensional Lebesgue integration over $\rho$: \begin{align*} n \alpha_n \int_0^r \rho^{\frac{1-n}{p-1}} \, d\Lambda^1(\rho) = n \alpha_n \left[ \frac{\rho^{\frac{p-n}{p-1}}}{\frac{p-n}{p-1}} \right]_0^r = n \alpha_n \frac{p-1}{p-n} r^{\frac{p-n}{p-1}} \end{align*} We raise this result to the power of $\frac{1}{q} = \frac{p-1}{p}$: \begin{align*} \left( \int_{B(x,r)} |x-z|^{-(n-1)q} \, d\Lambda^n(z) \right)^{\frac{1}{q}} = \left( n \alpha_n \frac{p-1}{p-n} \right)^{\frac{p-1}{p}} r^{1 - \frac{n}{p}} \end{align*} We define the geometric measure constant $C_1 = \frac{1}{n \alpha_n} \left( n \alpha_n \frac{p-1}{p-n} \right)^{\frac{p-1}{p}}$ and substitute this bounding expression back into the initial inequality from Step 1 to secure the claim: \begin{align*} \frac{1}{\int_{B(x,r)} 1 \, d\Lambda^n(y)} \int_{B(x,r)} |u(x) - u(y)| \, d\Lambda^n(y) \le C_1 r^{1 - \frac{n}{p}} \left( \int_{B(x,r)} |\nabla u(z)|^p \, d\Lambda^n(z) \right)^{\frac{1}{p}} \end{align*} [/proof]

**Step 3** [claim] There exists a constant $C_2 > 0$ such that for all $x, y \in \mathbb{R}^n$, the function $u$ satisfies the uniform Hölder continuity bound: \begin{align*} |u(x) - u(y)| \le C_2 |x-y|^{1 - \frac{n}{p}} \|\nabla u\|_{L^p(\mathbb{R}^n)} \end{align*} [/claim] [proof] Fix $x, y \in \mathbb{R}^n$ such that $x \neq y$, and let $r = |x-y|$. Define the intersection of the two metric balls $W = B(x,r) \cap B(y,r) \subseteq \mathbb{R}^n$.

We define the average of $u$ over the set $W$ with respect to the Lebesgue measure $\Lambda^n$ as: \begin{align*} (u)_W = \frac{1}{\int_W 1 \, d\Lambda^n(z)} \int_W u(z) \, d\Lambda^n(z) \end{align*} We apply the triangle inequality: \begin{align*} |u(x) - u(y)| = |u(x) - (u)_W + (u)_W - u(y)| \le |u(x) - (u)_W| + |u(y) - (u)_W| \end{align*} We bound the first term by bringing the absolute value inside the measure integral: \begin{align*} |u(x) - (u)_W| = \left| \frac{1}{\int_W 1 \, d\Lambda^n(z)} \int_W (u(x) - u(z)) \, d\Lambda^n(z) \right| \le \frac{1}{\int_W 1 \, d\Lambda^n(z)} \int_W |u(x) - u(z)| \, d\Lambda^n(z) \end{align*} Since $W \subseteq B(x,r)$, we expand the domain of integration to the entire ball $B(x,r)$: \begin{align*} \frac{1}{\int_W 1 \, d\Lambda^n(z)} \int_W |u(x) - u(z)| \, d\Lambda^n(z) \le \frac{1}{\int_W 1 \, d\Lambda^n(z)} \int_{B(x,r)} |u(x) - u(z)| \, d\Lambda^n(z) \end{align*} The Lebesgue measure of the intersection $W$ scales proportionally to $r^n$, ensuring there exists a dimensional constant $\gamma_n > 0$ such that $\int_W 1 \, d\Lambda^n(z) = \gamma_n r^n$. Because $\int_{B(x,r)} 1 \, d\Lambda^n(z) = \alpha_n r^n$, we substitute the exact relationship $\int_W 1 \, d\Lambda^n(z) = \frac{\gamma_n}{\alpha_n} \int_{B(x,r)} 1 \, d\Lambda^n(z)$ into the inequality: \begin{align*} \frac{1}{\int_W 1 \, d\Lambda^n(z)} \int_{B(x,r)} |u(x) - u(z)| \, d\Lambda^n(z) = \frac{\alpha_n}{\gamma_n} \frac{1}{\int_{B(x,r)} 1 \, d\Lambda^n(z)} \int_{B(x,r)} |u(x) - u(z)| \, d\Lambda^n(z) \end{align*} We apply the geometric measure estimate from Step 2, recognizing that $y$ in Step 2 is merely a dummy integration variable, which we now write as $z$: \begin{align*} \frac{\alpha_n}{\gamma_n} \frac{1}{\int_{B(x,r)} 1 \, d\Lambda^n(z)} \int_{B(x,r)} |u(x) - u(z)| \, d\Lambda^n(z) \le \frac{\alpha_n C_1}{\gamma_n} r^{1 - \frac{n}{p}} \left( \int_{B(x,r)} |\nabla u(z)|^p \, d\Lambda^n(z) \right)^{\frac{1}{p}} \end{align*} We strictly enlarge the integration domain of the gradient norm to the global manifold $\mathbb{R}^n$: \begin{align*} |u(x) - (u)_W| \le \frac{\alpha_n C_1}{\gamma_n} r^{1 - \frac{n}{p}} \|\nabla u\|_{L^p(\mathbb{R}^n)} \end{align*} By symmetric reasoning centered strictly at $y$, the second term satisfies an identical upper bound: \begin{align*} |u(y) - (u)_W| \le \frac{\alpha_n C_1}{\gamma_n} r^{1 - \frac{n}{p}} \|\nabla u\|_{L^p(\mathbb{R}^n)} \end{align*} Summing these bounds and defining the compound constant $C_2 = 2 \frac{\alpha_n C_1}{\gamma_n}$, we establish: \begin{align*} |u(x) - u(y)| \le C_2 r^{1 - \frac{n}{p}} \|\nabla u\|_{L^p(\mathbb{R}^n)} = C_2 |x-y|^{1 - \frac{n}{p}} \|\nabla u\|_{L^p(\mathbb{R}^n)} \end{align*} [/proof]

**Step 4** [claim] There exists a constant $C_3 > 0$ such that the function $u$ satisfies the uniform spatial bound $\|u\|_{L^\infty(\mathbb{R}^n)} \le C_3 \|u\|_{W^{1,p}(\mathbb{R}^n)}$. [/claim] [proof] Fix an arbitrary point $x \in \mathbb{R}^n$. We bound $|u(x)|$ by adding and subtracting the Lebesgue average over the open unit ball $B(x,1)$ and applying the triangle inequality: \begin{align*} |u(x)| \le |u(x) - (u)_{B(x,1)}| + |(u)_{B(x,1)}| \end{align*} For the first term, we set $r = 1$ in the fundamental inequality established in Step 2: \begin{align*} |u(x) - (u)_{B(x,1)}| \le C_1 \left( \int_{B(x,1)} |\nabla u(z)|^p \, d\Lambda^n(z) \right)^{\frac{1}{p}} \le C_1 \|\nabla u\|_{L^p(\mathbb{R}^n)} \end{align*} For the second term, we apply Hölder's inequality over the strict domain $B(x,1)$ with functions $u$ and the constant function $1$, integrating with respect to the Lebesgue measure $\Lambda^n$: \begin{align*} |(u)_{B(x,1)}| \le \frac{1}{\alpha_n} \int_{B(x,1)} |u(z)| \cdot 1 \, d\Lambda^n(z) \le \frac{1}{\alpha_n} \left( \int_{B(x,1)} |u(z)|^p \, d\Lambda^n(z) \right)^{\frac{1}{p}} \left( \int_{B(x,1)} 1^q \, d\Lambda^n(z) \right)^{\frac{1}{q}} \end{align*} Since the exact Lebesgue measure integral of $1$ evaluates to the volume $\alpha_n$, the rightmost factor scales to $\alpha_n^{\frac{1}{q}} = \alpha_n^{\frac{p-1}{p}}$. Thus: \begin{align*} |(u)_{B(x,1)}| \le \frac{1}{\alpha_n} \alpha_n^{\frac{p-1}{p}} \left( \int_{B(x,1)} |u(z)|^p \, d\Lambda^n(z) \right)^{\frac{1}{p}} = \alpha_n^{-\frac{1}{p}} \|u\|_{L^p(B(x,1))} \le \alpha_n^{-\frac{1}{p}} \|u\|_{L^p(\mathbb{R}^n)} \end{align*} Combining these bounds and taking the global supremum over all $x \in \mathbb{R}^n$, we isolate the uniform bound by defining $C_3 = \max\{C_1, \alpha_n^{-\frac{1}{p}}\}$: \begin{align*} \|u\|_{L^\infty(\mathbb{R}^n)} \le C_1 \|\nabla u\|_{L^p(\mathbb{R}^n)} + \alpha_n^{-\frac{1}{p}} \|u\|_{L^p(\mathbb{R}^n)} \le C_3 \left( \|\nabla u\|_{L^p(\mathbb{R}^n)} + \|u\|_{L^p(\mathbb{R}^n)} \right) = C_3 \|u\|_{W^{1,p}(\mathbb{R}^n)} \end{align*} [/proof]

**Step 5** [claim] The global Hölder norm of the function $u$ satisfies $\|u\|_{C^{0, 1-\frac{n}{p}}(\mathbb{R}^n)} \le C \|u\|_{W^{1,p}(\mathbb{R}^n)}$ for a universal constant $C > 0$. [/claim] [proof] The Hölder norm is uniquely defined as the direct sum of the supremum norm and the Hölder seminorm: \begin{align*} \|u\|_{C^{0, 1-\frac{n}{p}}(\mathbb{R}^n)} = \|u\|_{L^\infty(\mathbb{R}^n)} + \sup_{x \neq y} \frac{|u(x) - u(y)|}{|x-y|^{1-\frac{n}{p}}} \end{align*} We substitute the bounds explicitly derived in Step 4 and Step 3 into this metric definition: \begin{align*} \|u\|_{C^{0, 1-\frac{n}{p}}(\mathbb{R}^n)} \le C_3 \|u\|_{W^{1,p}(\mathbb{R}^n)} + C_2 \|\nabla u\|_{L^p(\mathbb{R}^n)} \end{align*} Because the $L^p$ norm of the gradient is bounded by the full Sobolev norm $\|\nabla u\|_{L^p(\mathbb{R}^n)} \le \|u\|_{W^{1,p}(\mathbb{R}^n)}$, we linearly combine the bounding terms: \begin{align*} \|u\|_{C^{0, 1-\frac{n}{p}}(\mathbb{R}^n)} \le (C_3 + C_2) \|u\|_{W^{1,p}(\mathbb{R}^n)} \end{align*} Defining the comprehensive constant $C = C_3 + C_2$ formally concludes the proof. [/proof]