[proofplan] The strategy mirrors the Sobolev case [$L^{p^*}$-Differentiability of Sobolev Functions](/theorems/3081) but with the structural complications of $BV$. We exploit the Sobolev-Poincaré inequality on balls in the $BV$ form: for $w \in BV(B(z, \rho))$, $(\fint_{B(z, \rho)} |w - w_{B(z, \rho)}|^{1^*} d\mathcal{L}^n)^{1/1^*} \le C_{\mathrm{SP}}\, |Dw|(B(z, \rho))/\rho^{n-1}$, where $1^* = n/(n-1)$. The argument splits the total variation $Du = \nabla u \cdot \mathcal{L}^n + D^s u$ into its absolutely continuous and singular parts, and works at points $x$ that are simultaneously: (i) Lebesgue points of the precise representative of $u$, (ii) $L^1$-Lebesgue points of $\nabla u$ in the strong sense $\fint_{B(x, \rho)} |\nabla u - \nabla u(x)| \, d\mathcal{L}^n \to 0$, and (iii) points at which $D^s u$ has zero density, $|D^s u|(B(x, \rho))/\rho^n \to 0$. The first two hold $\mathcal{L}^n$-a.e. by the Lebesgue Differentiation Theorem; the third holds $\mathcal{L}^n$-a.e. because $D^s u \perp \mathcal{L}^n$ implies the Radon-Nikodym derivative of $|D^s u|$ with respect to $\mathcal{L}^n$ vanishes a.e. As in the Sobolev case, we decompose $B(x, r)$ into dyadic shells, apply the $BV$ Sobolev-Poincaré on each dyadic ball to the Taylor remainder $v_x$ and to the auxiliary function $\tilde u_x$, and telescope the local means. The Hardy-Littlewood maximal function provides a summable dominator, and a parameter-version Dominated Convergence argument concludes $\Phi(x, r) \to 0$. Approximate differentiability follows by Markov's inequality. [/proofplan]

[step:Set up: Sobolev conjugate exponent, precise representatives, and the $BV$ structure] Throughout, $\Omega \subseteq \mathbb{R}^n$ is open, $n \ge 2$, and $u \in BV(\Omega)$. Define the Sobolev conjugate exponent \begin{align*} 1^* := \frac{n}{n - 1}, \end{align*} which satisfies $1^* > 1$ for $n \ge 2$ and is the Sobolev conjugate of $p = 1$ (in the limiting sense). By the structure theorem for $BV$ functions, the distributional gradient $Du$ is an $\mathbb{R}^n$-valued Radon measure of finite total variation $|Du|(\Omega) < \infty$, and decomposes as \begin{align*} Du = D^a u + D^s u, \qquad D^a u = \nabla u \cdot \mathcal{L}^n, \end{align*} where $D^a u$ is the absolutely continuous part with respect to $\mathcal{L}^n$ and $D^s u \perp \mathcal{L}^n$ is the singular part. The Radon-Nikodym density $\nabla u \in L^1_{\mathrm{loc}}(\Omega; \mathbb{R}^n)$. Define the precise representative of $u$, \begin{align*} u^*: \Omega &\to \mathbb{R} \cup \{*\} \\ x &\mapsto \begin{cases} \displaystyle \lim_{\rho \to 0^+} \fint_{B(x, \rho)} u(y) \, d\mathcal{L}^n(y) & \text{if the limit exists} \\ 0 & \text{otherwise.} \end{cases} \end{align*} By the Lebesgue Differentiation Theorem applied to $u \in L^1_{\mathrm{loc}}(\Omega)$, $u^*(x) = u(x)$ for $\mathcal{L}^n$-a.e. $x$. For $\nabla u \in L^1_{\mathrm{loc}}(\Omega; \mathbb{R}^n)$ define analogously the componentwise precise representative $(\nabla u)^*$. Henceforth replace $u$ by $u^*$ and $\nabla u$ by $(\nabla u)^*$, writing $u(x)$, $\nabla u(x)$ for the pointwise values without further comment. Define the Hardy-Littlewood maximal function \begin{align*} M f: \Omega &\to [0, \infty] \\ x &\mapsto \sup_{r > 0 : B(x, r) \subseteq \Omega} \fint_{B(x, r)} |f(y)| \, d\mathcal{L}^n(y). \end{align*} Fix $x_0 \in \Omega$ and $r_0 > 0$ with $\overline{B}(x_0, 2 r_0) \subseteq \Omega$. We prove the conclusion at $\mathcal{L}^n$-a.e. point of $B(x_0, r_0)$; since $\Omega$ admits a countable cover by such balls, the result extends to $\Omega$. Set \begin{align*} \Omega_0 := \{ x \in B(x_0, r_0) : (i), (ii), (iii), (iv) \text{ hold} \}, \end{align*} where: - (i) $x$ is a Lebesgue point of $u$, so $\lim_{\rho \to 0} \fint_{B(x, \rho)} u \, d\mathcal{L}^n = u(x)$; - (ii) $x$ is an $L^1$-Lebesgue point of $\nabla u$, so $\lim_{\rho \to 0} \fint_{B(x, \rho)} |\nabla u - \nabla u(x)| \, d\mathcal{L}^n = 0$; - (iii) $x$ has zero singular density, $\lim_{\rho \to 0} |D^s u|(B(x, \rho)) / \rho^n = 0$; - (iv) $M(|\nabla u - \nabla u(x)| \mathbb{1}_{B(x_0, 2 r_0)})(x) < \infty$ and $M^*(|D^s u|)(x) := \sup_{0 < \rho \le r_0} |D^s u|(B(x, \rho))/(\omega_n \rho^n) < \infty$. Conditions (i), (ii) hold $\mathcal{L}^n$-a.e. on $B(x_0, r_0)$ by the Lebesgue Differentiation Theorem applied to $u \in L^1_{\mathrm{loc}}$ and to $|\nabla u - \nabla u(x)| \in L^1_{\mathrm{loc}}$ (pointwise in $x$). Condition (iii) holds $\mathcal{L}^n$-a.e. by the Besicovitch differentiation theorem: since $|D^s u|$ is a Radon measure singular with respect to $\mathcal{L}^n$, the Radon-Nikodym derivative $D(|D^s u|, \mathcal{L}^n)(x) = \lim_{\rho \to 0} |D^s u|(B(x, \rho)) / (\omega_n \rho^n) = 0$ for $\mathcal{L}^n$-a.e. $x$. Condition (iv) holds $\mathcal{L}^n$-a.e. by the weak (1,1) inequality applied to $|\nabla u - \nabla u(x)| \mathbb{1}_{B(x_0, 2 r_0)} \in L^1$ and to the finite Radon measure $|D^s u| \mathbb{1}_{B(x_0, 2 r_0)}$. Hence $\mathcal{L}^n(B(x_0, r_0) \setminus \Omega_0) = 0$. We work with $x \in \Omega_0$ and $0 < r \le r_0$ so that $\overline{B}(x, r) \subseteq B(x_0, 2 r_0) \subseteq \Omega$. [/step]

[step:Define the Taylor remainder $v_x$ and reduce $\Phi(x, r)$ to dyadic-ball averages] For $x \in \Omega_0$ and $0 < r \le r_0$, define the Taylor remainder \begin{align*} v_x: B(x, r) &\to \mathbb{R} \\ y &\mapsto u(y) - u(x) - \nabla u(x) \cdot (y - x). \end{align*} Since the affine map $y \mapsto u(x) + \nabla u(x) \cdot (y - x)$ is smooth, $v_x \in BV(B(x, r))$ with distributional gradient \begin{align*} D v_x = D u - \nabla u(x) \cdot \mathcal{L}^n = (\nabla u - \nabla u(x)) \cdot \mathcal{L}^n + D^s u \end{align*} on $B(x, r)$. The total variation of $D v_x$ on a Borel set $E \subseteq B(x, r)$ satisfies \begin{align*} |D v_x|(E) \le \int_E |\nabla u(y) - \nabla u(x)| \, d\mathcal{L}^n(y) + |D^s u|(E). \end{align*} Write $\omega_n := \mathcal{L}^n(B(0, 1))$ so $\mathcal{L}^n(B(x, \rho)) = \omega_n \rho^n$. For $0 < \mathcal{L}^n(E) < \infty$ and measurable $f$ on $E$, write $f_E := \fint_E f \, d\mathcal{L}^n$. The goal is: \begin{align*} \Phi(x, r) := \fint_{B(x, r)} \left( \frac{|v_x(y)|}{|y - x|} \right)^{1^*} d\mathcal{L}^n(y) \to 0 \qquad (r \to 0^+) \text{ for every } x \in \Omega_0. \end{align*} For $k \in \mathbb{N}_0 = \{0, 1, 2, \dots\}$, let $A_k := B(x, 2^{-k} r) \setminus B(x, 2^{-k-1} r)$ be the dyadic shell. Then $B(x, r) \setminus \{x\} = \bigsqcup_{k \ge 0} A_k$ (disjoint), and $\mathcal{L}^n(\{x\}) = 0$. For $y \in A_k$ we have $|y - x| \ge 2^{-k-1} r$, so $|y - x|^{-1^*} \le 2^{(k+1) 1^*} r^{-1^*}$. Splitting $\Phi(x, r)$ over shells: \begin{align*} \Phi(x, r) &= \frac{1}{\omega_n r^n} \sum_{k = 0}^\infty \int_{A_k} \frac{|v_x(y)|^{1^*}}{|y - x|^{1^*}} \, d\mathcal{L}^n(y) \\ &\le \frac{1}{\omega_n r^n} \sum_{k = 0}^\infty \frac{2^{(k+1) 1^*}}{r^{1^*}} \int_{A_k} |v_x|^{1^*} d\mathcal{L}^n. \end{align*} Since $A_k \subseteq B(x, 2^{-k} r)$ and $\mathcal{L}^n(B(x, 2^{-k} r)) = \omega_n (2^{-k} r)^n$, \begin{align*} \int_{A_k} |v_x|^{1^*} d\mathcal{L}^n \le \omega_n (2^{-k} r)^n \fint_{B(x, 2^{-k} r)} |v_x|^{1^*} d\mathcal{L}^n. \end{align*} Combining, \begin{align*} \Phi(x, r) \le \frac{2^{1^*}}{r^{1^*}} \sum_{k = 0}^\infty 2^{k(1^* - n)} \fint_{B(x, 2^{-k} r)} |v_x|^{1^*} d\mathcal{L}^n. \end{align*} This reduces the problem to bounding the dyadic-ball averages $\fint_{B(x, 2^{-k} r)} |v_x|^{1^*} \, d\mathcal{L}^n$. We will need the linear-mean-zero observation. By the change of variables $y = x + z$ on $B(x, \rho)$, and using that the Lebesgue measure on $B(0, \rho)$ is invariant under $z \mapsto -z$, \begin{align*} \fint_{B(x, \rho)} \nabla u(x) \cdot (y - x) \, d\mathcal{L}^n(y) = \nabla u(x) \cdot \fint_{B(0, \rho)} z \, d\mathcal{L}^n(z) = 0. \end{align*} Hence for the auxiliary function \begin{align*} \tilde u_x: B(x, r) &\to \mathbb{R} \\ y &\mapsto u(y) - \nabla u(x) \cdot (y - x), \end{align*} we have $\tilde u_x \in BV(B(x, r))$ with distributional gradient $D \tilde u_x = (\nabla u - \nabla u(x)) \cdot \mathcal{L}^n + D^s u$, and the linear-mean-zero observation gives \begin{align*} (\tilde u_x)_{B(x, \rho)} = u_{B(x, \rho)} - \nabla u(x) \cdot \fint_{B(x, \rho)}(y - x) \, d\mathcal{L}^n(y) = u_{B(x, \rho)}. \end{align*} [/step]

[step:Key lemma: telescoping bound on the local mean $m_k(r)$ in terms of gradient and singular averages] Define \begin{align*} m_k(r) &:= (v_x)_{B(x, 2^{-k} r)} = \fint_{B(x, 2^{-k} r)} v_x \, d\mathcal{L}^n, \\ G_j(x, r) &:= \fint_{B(x, 2^{-j} r)} |\nabla u(y) - \nabla u(x)| \, d\mathcal{L}^n(y), \\ H_j(x, r) &:= \frac{|D^s u|(B(x, 2^{-j} r))}{\omega_n (2^{-j} r)^n}. \end{align*} By the linear-mean-zero observation of Step 2, $m_k(r) = u_{B(x, 2^{-k} r)} - u(x)$. [claim:Telescoping bound on $|m_k(r)|$] There exists a constant $C = C(n)$ such that for every $x \in \Omega_0$, every $0 < r \le r_0$, and every $k \in \mathbb{N}_0$, \begin{align*} |m_k(r)| \le C \cdot r \sum_{j = k}^\infty 2^{-j} \big( G_j(x, r) + H_j(x, r) \big). \end{align*} [/claim] *Proof of claim.* Recall the Sobolev-Poincaré inequality on balls for $BV$ (a standard result; for $w \in BV(B(z, \rho))$, \begin{align*} \left( \fint_{B(z, \rho)} |w - w_{B(z, \rho)}|^{1^*} d\mathcal{L}^n \right)^{1/1^*} \le C_{\mathrm{SP}} \cdot \frac{|Dw|(B(z, \rho))}{\rho^{n-1}}, \end{align*} where $C_{\mathrm{SP}} = C_{\mathrm{SP}}(n)$ is a dimensional constant. Apply this to $\tilde u_x \in BV(B(x, r))$ on the ball $B(x, 2^{-j} r)$, using $|D \tilde u_x|(B(x, 2^{-j} r)) \le \int_{B(x, 2^{-j} r)} |\nabla u - \nabla u(x)| \, d\mathcal{L}^n + |D^s u|(B(x, 2^{-j} r))$ and $(\tilde u_x)_{B(x, 2^{-j} r)} = u_{B(x, 2^{-j} r)}$: \begin{align*} &\left( \fint_{B(x, 2^{-j} r)} \big| \tilde u_x(y) - u_{B(x, 2^{-j} r)} \big|^{1^*} d\mathcal{L}^n(y) \right)^{1/1^*} \\ &\quad \le C_{\mathrm{SP}} \cdot \frac{|D \tilde u_x|(B(x, 2^{-j} r))}{(2^{-j} r)^{n-1}} \\ &\quad \le C_{\mathrm{SP}} \cdot \frac{1}{(2^{-j} r)^{n-1}} \big( \omega_n (2^{-j} r)^n \cdot G_j(x, r) + \omega_n (2^{-j} r)^n \cdot H_j(x, r) \big) \\ &\quad = C_{\mathrm{SP}} \omega_n \cdot 2^{-j} r \cdot \big( G_j(x, r) + H_j(x, r) \big). \end{align*} For each $j \ge 0$, since $B(x, 2^{-j-1} r) \subseteq B(x, 2^{-j} r)$ with volume ratio $2^{-n}$, by Jensen's inequality applied to the conjugate-exponent direction (which gives $|\fint f| \le \fint |f| \le (\fint |f|^{1^*})^{1/1^*}$ for $1^* \ge 1$), \begin{align*} \big| u_{B(x, 2^{-j-1} r)} - u_{B(x, 2^{-j} r)} \big| &= \big| (\tilde u_x)_{B(x, 2^{-j-1} r)} - (\tilde u_x)_{B(x, 2^{-j} r)} \big| \\ &\le \fint_{B(x, 2^{-j-1} r)} \big| \tilde u_x - (\tilde u_x)_{B(x, 2^{-j} r)} \big| \, d\mathcal{L}^n \\ &\le 2^n \fint_{B(x, 2^{-j} r)} \big| \tilde u_x - (\tilde u_x)_{B(x, 2^{-j} r)} \big| \, d\mathcal{L}^n \\ &\le 2^n \left( \fint_{B(x, 2^{-j} r)} \big| \tilde u_x - (\tilde u_x)_{B(x, 2^{-j} r)} \big|^{1^*} d\mathcal{L}^n \right)^{1/1^*} \\ &\le 2^n C_{\mathrm{SP}} \omega_n \cdot 2^{-j} r \cdot \big( G_j(x, r) + H_j(x, r) \big). \end{align*} The first equality used $(\tilde u_x)_{B(x, \rho)} = u_{B(x, \rho)}$ for both balls. Now apply condition (i) at the Lebesgue point $x$ of $u$: $u_{B(x, 2^{-j} r)} \to u(x)$ as $j \to \infty$. Hence the telescoping series \begin{align*} u(x) - u_{B(x, 2^{-k} r)} = \sum_{j = k}^\infty \big( u_{B(x, 2^{-j-1} r)} - u_{B(x, 2^{-j} r)} \big) \end{align*} converges absolutely (granted the bound just derived; we verify summability below), and the triangle inequality gives \begin{align*} |m_k(r)| = \big| u_{B(x, 2^{-k} r)} - u(x) \big| \le \sum_{j = k}^\infty \big| u_{B(x, 2^{-j-1} r)} - u_{B(x, 2^{-j} r)} \big| \le 2^n C_{\mathrm{SP}} \omega_n \cdot r \sum_{j = k}^\infty 2^{-j} \big( G_j(x, r) + H_j(x, r) \big). \end{align*} Absolute convergence of the right side follows from $G_j(x, r) \le M_*(x)$ and $H_j(x, r) \le M^*(|D^s u|)(x)$ (both finite by (iv)), together with $\sum_j 2^{-j} < \infty$. Setting $C := 2^n C_{\mathrm{SP}} \omega_n$ proves the claim. $\square$ [/step]

[step:Apply Sobolev-Poincaré on each dyadic ball and combine with the telescoping bound] Apply the $BV$ Sobolev-Poincaré inequality to $v_x \in BV(B(x, r))$ on the ball $B(x, 2^{-k} r)$, using $|D v_x|(B(x, 2^{-k} r)) \le \omega_n (2^{-k} r)^n (G_k(x, r) + H_k(x, r))$: \begin{align*} \left( \fint_{B(x, 2^{-k} r)} |v_x - m_k(r)|^{1^*} d\mathcal{L}^n \right)^{1/1^*} \le C_{\mathrm{SP}} \omega_n \cdot 2^{-k} r \cdot \big( G_k(x, r) + H_k(x, r) \big). \end{align*} By the triangle inequality in $L^{1^*}(B(x, 2^{-k} r))$ for the decomposition $v_x = (v_x - m_k(r)) + m_k(r)$, \begin{align*} \left( \fint_{B(x, 2^{-k} r)} |v_x|^{1^*} d\mathcal{L}^n \right)^{1/1^*} \le C_{\mathrm{SP}} \omega_n \cdot 2^{-k} r \cdot \big( G_k(x, r) + H_k(x, r) \big) + |m_k(r)|. \end{align*} Raising to the $1^*$-th power (a pure algebraic manipulation; not Hölder) and using $(a + b)^{1^*} \le 2^{1^*}(a^{1^*} + b^{1^*})$ for $a, b \ge 0$: \begin{align*} \fint_{B(x, 2^{-k} r)} |v_x|^{1^*} d\mathcal{L}^n \le 2^{1^*} \left( (C_{\mathrm{SP}} \omega_n)^{1^*} (2^{-k} r)^{1^*} \big( G_k(x, r) + H_k(x, r) \big)^{1^*} + |m_k(r)|^{1^*} \right). \end{align*} [/step]

[step:Sum the dyadic shell estimates and conclude $\Phi(x, r) \to 0$ via dominated convergence] Substitute the Step 4 dyadic-ball bound into the Step 2 reduction: \begin{align*} \Phi(x, r) &\le \frac{2^{1^*}}{r^{1^*}} \sum_{k = 0}^\infty 2^{k(1^* - n)} \fint_{B(x, 2^{-k} r)} |v_x|^{1^*} d\mathcal{L}^n \\ &\le \frac{2^{1^*}}{r^{1^*}} \sum_{k = 0}^\infty 2^{k(1^* - n)} \cdot 2^{1^*} \left( (C_{\mathrm{SP}} \omega_n)^{1^*} (2^{-k} r)^{1^*} (G_k + H_k)^{1^*} + |m_k(r)|^{1^*} \right) \\ &= 2^{2 \cdot 1^*} (C_{\mathrm{SP}} \omega_n)^{1^*} \underbrace{\sum_{k = 0}^\infty 2^{-k n} (G_k(x, r) + H_k(x, r))^{1^*}}_{=: S_I(r)} + \underbrace{\frac{2^{2 \cdot 1^*}}{r^{1^*}} \sum_{k = 0}^\infty 2^{k(1^* - n)} |m_k(r)|^{1^*}}_{=: S_{II}(r)}, \end{align*} where in the first sum we used $2^{k(1^* - n)} (2^{-k} r)^{1^*} / r^{1^*} = 2^{-k n}$. Note $1^* = n/(n-1)$, so $1^* - n = (n - n(n-1))/(n-1) = (n - n^2 + n)/(n-1) = n(2-n)/(n-1)$, which is non-positive for $n \ge 2$. We do not need the sign explicitly; the relevant identity is the cancellation $2^{k(1^* - n)} (2^{-k})^{1^*} = 2^{-k n}$, which we have already used. We treat $S_I(r)$ and $S_{II}(r)$ separately. Throughout this step the constants depend on $x$ via $M_*(x) := M(|\nabla u - \nabla u(x)| \mathbb{1}_{B(x_0, 2r_0)})(x)$ and $M^*(|D^s u|)(x)$, both finite by condition (iv); we are proving pointwise (in $x$) convergence, not a uniform bound. **Sum I.** For each fixed $k$, condition (ii) gives $G_k(x, r) \to 0$ as $r \to 0^+$, and condition (iii) gives $H_k(x, r) \to 0$ as $r \to 0^+$ (since $|D^s u|(B(x, 2^{-k} r)) / (2^{-k} r)^n \to 0$, the rescaling factor $2^{-k}$ being constant for fixed $k$). Hence $(G_k + H_k)^{1^*} \to 0$. The maximal-function dominators give, for $2^{-k} r \le r_0$, \begin{align*} G_k(x, r) \le M_*(x), \quad H_k(x, r) \le M^*(|D^s u|)(x), \end{align*} so $(G_k + H_k)^{1^*} \le (M_*(x) + M^*(|D^s u|)(x))^{1^*}$, a constant independent of $k$ and $r$. Therefore $2^{-k n} (G_k + H_k)^{1^*} \le 2^{-k n} (M_*(x) + M^*(|D^s u|)(x))^{1^*}$, summable in $k$. By the Dominated Convergence Theorem applied to the counting measure on $\mathbb{N}_0$ weighted by $2^{-k n}$ with parameter $r \to 0^+$ (the parameter version: for every sequence $r_\ell \to 0^+$, $(G_k(x, r_\ell) + H_k(x, r_\ell))^{1^*} \to 0$ pointwise in $k$ and is dominated by a summable sequence, so the weighted sum converges to $0$), \begin{align*} S_I(r) \to 0 \qquad (r \to 0^+). \end{align*} **Sum II.** Apply the telescoping bound from the Step 3 claim: \begin{align*} |m_k(r)|^{1^*} \le (2^n C_{\mathrm{SP}} \omega_n)^{1^*} r^{1^*} \left( \sum_{j = k}^\infty 2^{-j} (G_j(x, r) + H_j(x, r)) \right)^{1^*}. \end{align*} Using $G_j(x, r) + H_j(x, r) \le M_*(x) + M^*(|D^s u|)(x) =: K_*(x)$, the inner sum is bounded uniformly by $\sum_{j \ge k} 2^{-j} K_*(x) = 2^{1 - k} K_*(x)$, so \begin{align*} |m_k(r)|^{1^*} \le (2^n C_{\mathrm{SP}} \omega_n)^{1^*} r^{1^*} \cdot 2^{(1 - k) 1^*} K_*(x)^{1^*}. \end{align*} Multiplying by $2^{k(1^* - n)} / r^{1^*}$ and using $2^{k(1^* - n)} \cdot 2^{(1 - k) 1^*} = 2^{1^*} \cdot 2^{-k n}$: \begin{align*} 2^{k(1^* - n)} \frac{|m_k(r)|^{1^*}}{r^{1^*}} \le (2^n C_{\mathrm{SP}} \omega_n)^{1^*} \cdot 2^{1^*} \cdot 2^{-k n} K_*(x)^{1^*}. \end{align*} This is the summable dominator (in $k$) for $S_{II}(r)$, independent of $r$. For pointwise convergence in $k$ as $r \to 0^+$: for each fixed $k$ and each fixed $j \ge k$, condition (ii) gives $G_j(x, r) \to 0$ and condition (iii) gives $H_j(x, r) \to 0$ as $r \to 0^+$; the dominator $2^{-j} K_*(x)$ is summable in $j$. By DCT (parameter version, on the counting measure on $\{j : j \ge k\}$), \begin{align*} \sum_{j = k}^\infty 2^{-j} (G_j(x, r) + H_j(x, r)) \to 0 \qquad (r \to 0^+) \text{ for each fixed } k. \end{align*} Hence $|m_k(r)|^{1^*} / r^{1^*} \to 0$ as $r \to 0^+$ for each fixed $k$, and consequently $2^{k(1^* - n)} |m_k(r)|^{1^*} / r^{1^*} \to 0$. By DCT once more (now on the counting measure on $\mathbb{N}_0$ with the dominator $(2^n C_{\mathrm{SP}} \omega_n)^{1^*} \cdot 2^{1^*} \cdot 2^{-k n} K_*(x)^{1^*}$), \begin{align*} S_{II}(r) \to 0 \qquad (r \to 0^+). \end{align*} **Conclusion.** $\Phi(x, r) \le 2^{2 \cdot 1^*} (C_{\mathrm{SP}} \omega_n)^{1^*} S_I(r) + S_{II}(r) \to 0$ as $r \to 0^+$ for every $x \in \Omega_0$. This is the $L^{1^*}$-differentiability statement: for $\mathcal{L}^n$-a.e. $x \in \Omega$, \begin{align*} \left( \fint_{B(x, r)} \left( \frac{|u(y) - u(x) - \nabla u(x) \cdot (y - x)|}{|y - x|} \right)^{1^*} d\mathcal{L}^n(y) \right)^{1/1^*} \to 0 \qquad (r \to 0^+). \end{align*} [/step]

[step:Approximate differentiability via Markov's inequality] For $\varepsilon > 0$ and $r > 0$, define \begin{align*} E_\varepsilon(r) := \left\{ y \in B(x, r) \setminus \{x\} : \frac{|u(y) - u(x) - \nabla u(x) \cdot (y - x)|}{|y - x|} > \varepsilon \right\}. \end{align*} By Markov's inequality (the elementary bound $\mu(\{|f| > \lambda\}) \le \lambda^{-q} \int |f|^q \, d\mu$ for $q \ge 1$, $\lambda > 0$, valid because $\mathbb{1}_{\{|f| > \lambda\}} \le |f|^q / \lambda^q$ pointwise) with $q = 1^*$, $\lambda = \varepsilon^{1^*}$, $\mu = \mathcal{L}^n$, applied to $y \mapsto |v_x(y)|^{1^*} / |y - x|^{1^*}$ on $B(x, r)$: \begin{align*} \mathcal{L}^n(E_\varepsilon(r)) \le \frac{1}{\varepsilon^{1^*}} \int_{B(x, r)} \frac{|v_x(y)|^{1^*}}{|y - x|^{1^*}} \, d\mathcal{L}^n(y). \end{align*} Dividing by $\mathcal{L}^n(B(x, r)) = \omega_n r^n$: \begin{align*} \frac{\mathcal{L}^n(E_\varepsilon(r))}{\mathcal{L}^n(B(x, r))} \le \frac{\Phi(x, r)}{\varepsilon^{1^*}}. \end{align*} By Step 5, $\Phi(x, r) \to 0$ as $r \to 0^+$ for every $x \in \Omega_0$, so \begin{align*} \lim_{r \to 0^+} \frac{\mathcal{L}^n(E_\varepsilon(r))}{\mathcal{L}^n(B(x, r))} = 0 \qquad \text{for every } \varepsilon > 0. \end{align*} This is the definition of approximate differentiability of $u$ at $x$ with approximate derivative $\nabla u(x)$. The conclusion holds for $\mathcal{L}^n$-a.e. $x \in \Omega$, completing the proof. (This is the $BV$ counterpart of [$L^{p^*}$-Differentiability of Sobolev Functions](/theorems/3081). The new ingredient relative to the $W^{1,1}$ case is condition (iii) — the vanishing density of the singular part $D^s u$ at $\mathcal{L}^n$-a.e. $x$ — which controls the contribution of $D^s u$ to the dyadic Sobolev-Poincaré estimate.) [/step]