[proofplan]
We compare the minimizer $u$ with a competitor that agrees with $u$ near $\partial B(x_0,R)$ and is equal to the constant $c$ on a smaller ball. The growth assumptions convert minimality into an estimate for the energy on an inner ball by the energy on an annulus plus an $L^p$ oscillation term. A hole-filling iteration over nested radii absorbs the annular gradient term and yields the stated estimate at radius $r$.
[/proofplan]
[step:Construct a cutoff competitor equal to $c$ on the inner ball]
Set $B_\rho:=B(x_0,\rho)$ for $0<\rho\le R$, and in particular $B_R:=B(x_0,R)$. Fix radii $r<s<t<R$ and a vector $c\in\mathbb R^m$. By the standard smooth cutoff construction for concentric Euclidean balls, choose $\eta\in C_c^\infty(B_t)$ such that $0\le \eta\le 1$, $\eta=1$ on $B_s$, and $|\nabla\eta(x)|\le 2/(t-s)$ for every $x\in B_t$.
Define $v:B_R\to\mathbb R^m$ by $v(x)=u(x)-\eta(x)(u(x)-c)$ for $x\in B_t$ and $v(x)=u(x)$ for $x\in B_R\setminus B_t$. Since $\eta$ has compact support in $B_t\subset B_R$, the difference $v-u=-\eta(u-c)$ belongs to $W^{1,p}_0(B_R;\mathbb R^m)$ by the zero-trace characterization of compactly supported Sobolev perturbations. Hence $\operatorname{Tr}v=\operatorname{Tr}u$ on $\partial B_R$, so $v\in\mathcal A_u$.
The weak product rule for multiplication by a smooth bounded function with bounded gradient gives
\begin{align*}
\nabla v=(1-\eta)\nabla u-(u-c)\otimes\nabla\eta
\end{align*}
on $B_t$. On $B_s$, we have $\eta=1$ and $\nabla\eta=0$ a.e., so $v=c$ and therefore $\nabla v=0$ a.e. on $B_s$.
[guided]
The competitor must be defined on the whole ball $B_R$, because the admissible class $\mathcal A_u$ is imposed on $\partial B_R$, not on $\partial B_t$. We therefore modify $u$ only inside $B_t$ and keep it unchanged outside $B_t$.
For $0<\rho\le R$, define $B_\rho:=B(x_0,\rho)$. Let $C_c^\infty(B_t)$ denote the space of smooth real-valued functions with compact support in $B_t$. Fix $r<s<t<R$ and $c\in\mathbb R^m$. We choose a smooth cutoff function $\eta\in C_c^\infty(B_t)$ with $0\le \eta\le 1$, $\eta=1$ on $B_s$, and
\begin{align*}
|\nabla\eta(x)|\le \frac{2}{t-s}
\end{align*}
for every $x\in B_t$. This is the standard Euclidean cutoff construction between two concentric balls.
Define $v:B_R\to\mathbb R^m$ by setting $v(x)=u(x)-\eta(x)(u(x)-c)$ for $x\in B_t$ and $v(x)=u(x)$ for $x\in B_R\setminus B_t$. Equivalently, on $B_t$ we have $v=c+(1-\eta)(u-c)$. Let $W^{1,p}_0(B_R;\mathbb R^m)$ denote the closure of $C_c^\infty(B_R;\mathbb R^m)$ in the $W^{1,p}(B_R;\mathbb R^m)$ norm. Since $\eta$ has compact support in $B_t$, the difference $v-u=-\eta(u-c)$ has zero trace on $\partial B_R$ and belongs to $W^{1,p}_0(B_R;\mathbb R^m)$. Thus $\operatorname{Tr}v=\operatorname{Tr}u$ on $\partial B_R$, so $v\in\mathcal A_u$.
The weak product rule applies because $\eta$ is smooth and bounded with bounded gradient, while $u-c\in W^{1,p}(B_t;\mathbb R^m)$. Hence on $B_t$,
\begin{align*}
\nabla v=(1-\eta)\nabla u-(u-c)\otimes\nabla\eta.
\end{align*}
On $B_s$, the cutoff is identically $1$, so $\nabla\eta=0$ a.e. on $B_s$ and $v=c$ there. Therefore $\nabla v=0$ a.e. on $B_s$. This is the point of the construction: the competitor deletes the gradient of $u$ in the inner ball and pays for that deletion only in the annulus $B_t\setminus B_s$.
[/guided]
[/step]
[step:Compare the energies and isolate the inner gradient term]
Because $F$ is finite and convex on the finite-dimensional space $\mathbb R^{m\times n}$, it is continuous and therefore Borel measurable. The maps $\nabla u$ and $\nabla v$ are measurable, so $F(\nabla u)$ and $F(\nabla v)$ are measurable. The growth bound $F(A)\le \Lambda(1+|A|^p)$ and the fact that $u,v\in W^{1,p}(B_R;\mathbb R^m)$ imply $F(\nabla u),F(\nabla v)\in L^1(B_R)$; the lower bound $F(A)\ge \lambda|A|^p-\Lambda$ controls the negative part.
Since $u$ minimizes $I$ over $\mathcal A_u$ and $v\in\mathcal A_u$,
\begin{align*}
\int_{B_R}F(\nabla u)\,d\mathcal L^n \le \int_{B_R}F(\nabla v)\,d\mathcal L^n
\end{align*}
Because $v=u$ a.e. on $B_R\setminus B_t$, subtracting the common outer contribution yields
\begin{align*}
\int_{B_t}F(\nabla u)\,d\mathcal L^n \le \int_{B_t}F(\nabla v)\,d\mathcal L^n
\end{align*}
Let $E_{s,t}:=B_t\setminus B_s$. Since $\nabla v=0$ a.e. on $B_s$,
\begin{align*}
\int_{B_s}F(\nabla u)\,d\mathcal L^n+\int_{E_{s,t}}F(\nabla u)\,d\mathcal L^n \le F(0)\mathcal L^n(B_s)+\int_{E_{s,t}}F(\nabla v)\,d\mathcal L^n
\end{align*}
The growth bounds give $F(0)\le \Lambda$ and $F(A)\ge \lambda|A|^p-\Lambda$, hence
\begin{align*}
\lambda\int_{B_s}|\nabla u|^p\,d\mathcal L^n \le 2\Lambda\mathcal L^n(B_s)+\Lambda\mathcal L^n(E_{s,t})+\int_{E_{s,t}}F(\nabla v)\,d\mathcal L^n
\end{align*}
Applying the upper growth bound on the last term and using $\mathcal L^n(B_s)\le \mathcal L^n(B_t)$ and $\mathcal L^n(E_{s,t})\le \mathcal L^n(B_t)$ gives
\begin{align*}
\lambda\int_{B_s}|\nabla u|^p\,d\mathcal L^n \le 3\Lambda\mathcal L^n(B_t)+\Lambda\int_{E_{s,t}}|\nabla v|^p\,d\mathcal L^n
\end{align*}
[/step]
[step:Estimate the annular gradient of the competitor]
For a.e. $x\in E_{s,t}$,
\begin{align*}
|\nabla v(x)|\le |\nabla u(x)|+|u(x)-c|\,|\nabla\eta(x)|
\end{align*}
Using $|a+b|^p\le 2^{p-1}(|a|^p+|b|^p)$ in $\mathbb R^{m\times n}$ and $|\nabla\eta|\le 2/(t-s)$, we obtain
\begin{align*}
|\nabla v(x)|^p\le 2^{p-1}|\nabla u(x)|^p+2^{2p-1}(t-s)^{-p}|u(x)-c|^p
\end{align*}
for a.e. $x\in E_{s,t}$. Integrating and combining with the previous step yields
\begin{align*}
\int_{B_s}|\nabla u|^p\,d\mathcal L^n \le C_0\int_{E_{s,t}}|\nabla u|^p\,d\mathcal L^n+\frac{C_0}{(t-s)^p}\int_{B_t}|u-c|^p\,d\mathcal L^n+C_0t^n
\end{align*}
where $C_0>0$ depends only on $n,p,\lambda,\Lambda$.
[guided]
The competitor gradient on the annulus contains two terms: the old gradient of $u$ and the transition cost caused by the cutoff. For a.e. $x\in E_{s,t}$, the identity
\begin{align*}
\nabla v=(1-\eta)\nabla u-(u-c)\otimes\nabla\eta
\end{align*}
and the bounds $0\le \eta\le 1$ and $|\nabla\eta(x)|\le 2/(t-s)$ give
\begin{align*}
|\nabla v(x)|\le |\nabla u(x)|+|u(x)-c|\,|\nabla\eta(x)|
\end{align*}
Applying the elementary convexity inequality $|a+b|^p\le 2^{p-1}(|a|^p+|b|^p)$ in the finite-dimensional normed space $\mathbb R^{m\times n}$ gives
\begin{align*}
|\nabla v(x)|^p\le 2^{p-1}|\nabla u(x)|^p+2^{2p-1}(t-s)^{-p}|u(x)-c|^p
\end{align*}
for a.e. $x\in E_{s,t}$. Substituting this estimate into the energy comparison from the previous step and dividing by $\lambda$ yields
\begin{align*}
\int_{B_s}|\nabla u|^p\,d\mathcal L^n\le C_0\int_{E_{s,t}}|\nabla u|^p\,d\mathcal L^n+\frac{C_0}{(t-s)^p}\int_{B_t}|u-c|^p\,d\mathcal L^n+C_0t^n
\end{align*}
Here the constant $C_0>0$ depends only on $n,p,\lambda,\Lambda$: the dependence on $n$ enters through the bound $\mathcal L^n(B_t)=\mathcal L^n(B(0,1))t^n$, and the remaining constants come from the growth bounds and the elementary $p$-power inequality.
[/guided]
[/step]
[step:Convert the annular estimate into a hole-filling inequality]
Define $\Phi:(0,R)\to[0,\infty)$ by
\begin{align*}
\Phi(\rho)=\int_{B_\rho}|\nabla u|^p\,d\mathcal L^n
\end{align*}
Since $\int_{E_{s,t}}|\nabla u|^p\,d\mathcal L^n=\Phi(t)-\Phi(s)$, the previous estimate gives
\begin{align*}
\Phi(s)\le C_0(\Phi(t)-\Phi(s))+\frac{C_0}{(t-s)^p}\int_{B_t}|u-c|^p\,d\mathcal L^n+C_0t^n
\end{align*}
After adding $C_0\Phi(s)$ to both sides and dividing by $1+C_0$, we obtain
\begin{align*}
\Phi(s)\le \theta\Phi(t)+\frac{C_1}{(t-s)^p}\int_{B_t}|u-c|^p\,d\mathcal L^n+C_1t^n
\end{align*}
where $\theta:=C_0/(1+C_0)\in(0,1)$ and $C_1>0$ is a structural constant depending only on $n,p,\lambda,\Lambda$.
[guided]
The annular term is the only obstruction to the final estimate. We rewrite it in a form where it can be absorbed by iteration.
Define $\Phi:(0,R)\to[0,\infty)$ by
\begin{align*}
\Phi(\rho)=\int_{B_\rho}|\nabla u|^p\,d\mathcal L^n
\end{align*}
Because $B_s\subset B_t$ and $E_{s,t}=B_t\setminus B_s$, additivity of the [Lebesgue integral](/page/Lebesgue%20Integral) over disjoint measurable sets gives
\begin{align*}
\int_{E_{s,t}}|\nabla u|^p\,d\mathcal L^n=\Phi(t)-\Phi(s)
\end{align*}
Substituting this identity into the annular estimate gives
\begin{align*}
\Phi(s)\le C_0(\Phi(t)-\Phi(s))+\frac{C_0}{(t-s)^p}\int_{B_t}|u-c|^p\,d\mathcal L^n+C_0t^n
\end{align*}
Now the term $-C_0\Phi(s)$ on the right is useful: it lets us move part of the left-hand energy to the left side. Adding $C_0\Phi(s)$ to both sides yields
\begin{align*}
(1+C_0)\Phi(s)\le C_0\Phi(t)+\frac{C_0}{(t-s)^p}\int_{B_t}|u-c|^p\,d\mathcal L^n+C_0t^n
\end{align*}
Dividing by $1+C_0$ gives
\begin{align*}
\Phi(s)\le \theta\Phi(t)+\frac{C_1}{(t-s)^p}\int_{B_t}|u-c|^p\,d\mathcal L^n+C_1t^n
\end{align*}
where $\theta:=C_0/(1+C_0)\in(0,1)$.
The strict inequality $\theta<1$ is the decisive point. It is what permits the hole-filling iteration to absorb the larger-radius energy into the left-hand side after summing over nested radii. The constant $C_1$ records only the structural dependence inherited from $C_0$, so subsequent enlargements of $C_1$ do not introduce dependence on $u,r,R,s,t$, or $c$.
[/guided]
[/step]
[step:Iterate between $r$ and $R$ to absorb the larger-radius energy]
Set
\begin{align*}
A:=C_1\int_{B_R}|u-c|^p\,d\mathcal L^n+C_1R^n(R-r)^p.
\end{align*}
For every $r\le s<t<R$, the inclusions $B_t\subset B_R$ and the bound $t^n\le R^n$ give
\begin{align*}
\Phi(s)\le \theta\Phi(t)+\frac{A}{(t-s)^p}.
\end{align*}
Choose $q\in(\theta^{1/p},1)$, so that $\theta q^{-p}<1$, and define $\rho_k:=r+(R-r)(1-q^k)$ for $k\in\mathbb N$. Applying the preceding inequality with $s=\rho_k$ and $t=\rho_{k+1}$ gives
\begin{align*}
\Phi(\rho_k)\le \theta\Phi(\rho_{k+1})+A(R-r)^{-p}(1-q)^{-p}q^{-kp}.
\end{align*}
For each $k$, multiply this inequality by $\theta^k$. Summing the resulting inequalities from $k=1$ to $N$ yields
\begin{align*}
\sum_{k=1}^N\theta^k\Phi(\rho_k)\le \sum_{k=1}^N\theta^{k+1}\Phi(\rho_{k+1})+A(R-r)^{-p}(1-q)^{-p}\sum_{k=1}^N(\theta q^{-p})^k.
\end{align*}
After cancellation of the telescoping terms,
\begin{align*}
\theta\Phi(\rho_1)\le \theta^{N+1}\Phi(\rho_{N+1})+A(R-r)^{-p}(1-q)^{-p}\sum_{k=1}^N(\theta q^{-p})^k.
\end{align*}
Because $u\in W^{1,p}(B_R;\mathbb R^m)$, we have
\begin{align*}
\Phi(\rho_{N+1})\le \int_{B_R}|\nabla u|^p\,d\mathcal L^n<\infty.
\end{align*}
Since $0<\theta<1$, it follows that $\theta^{N+1}\Phi(\rho_{N+1})\to 0$ as $N\to\infty$. Hence
\begin{align*}
\Phi(\rho_1)\le C_2A(R-r)^{-p}
\end{align*}
for $C_2:=\theta^{-1}(1-q)^{-p}\sum_{k=1}^{\infty}(\theta q^{-p})^k$, which depends only on $p$ and $\theta$. Since $r<\rho_1$ and $B_r\subset B_{\rho_1}$,
\begin{align*}
\int_{B_r}|\nabla u|^p\,d\mathcal L^n\le \Phi(\rho_1)\le \frac{C_2C_1}{(R-r)^p}\int_{B_R}|u-c|^p\,d\mathcal L^n+C_2C_1R^n.
\end{align*}
Renaming $C:=C_1C_2$, and recalling that $C_1$ depends only on $n,p,\lambda,\Lambda$ while $C_2$ depends only on $p$ and $\theta$, we obtain a constant depending only on $n,p,\lambda,\Lambda$, hence also only on $n,m,p,\lambda,\Lambda$. This is the desired Caccioppoli estimate.
[guided]
The recursive inequality from the previous step has the form needed for hole filling: a smaller-radius energy is bounded by a strict fraction of a larger-radius energy plus an explicit error. Define
\begin{align*}
A:=C_1\int_{B_R}|u-c|^p\,d\mathcal L^n+C_1R^n(R-r)^p.
\end{align*}
For every $r\le s<t<R$, the inclusions $B_t\subset B_R$ and the estimate $t^n\le R^n$ give
\begin{align*}
\Phi(s)\le \theta\Phi(t)+\frac{A}{(t-s)^p}.
\end{align*}
Choose $q\in(\theta^{1/p},1)$, so $\theta q^{-p}<1$, and define
\begin{align*}
\rho_k:=r+(R-r)(1-q^k)
\end{align*}
for $k\in\mathbb N$. Then $r<\rho_k<\rho_{k+1}<R$ and
\begin{align*}
\rho_{k+1}-\rho_k=(R-r)(1-q)q^k.
\end{align*}
Applying the recursive estimate with $s=\rho_k$ and $t=\rho_{k+1}$ gives
\begin{align*}
\Phi(\rho_k)\le \theta\Phi(\rho_{k+1})+A(R-r)^{-p}(1-q)^{-p}q^{-kp}.
\end{align*}
Multiplying by $\theta^k$ and summing from $k=1$ to $N$ gives a telescoping inequality:
\begin{align*}
\sum_{k=1}^N\theta^k\Phi(\rho_k)\le \sum_{k=1}^N\theta^{k+1}\Phi(\rho_{k+1})+A(R-r)^{-p}(1-q)^{-p}\sum_{k=1}^N(\theta q^{-p})^k.
\end{align*}
After cancelling the common terms, we obtain
\begin{align*}
\theta\Phi(\rho_1)\le \theta^{N+1}\Phi(\rho_{N+1})+A(R-r)^{-p}(1-q)^{-p}\sum_{k=1}^N(\theta q^{-p})^k.
\end{align*}
Because $u\in W^{1,p}(B_R;\mathbb R^m)$, the quantity $\int_{B_R}|\nabla u|^p\,d\mathcal L^n$ is finite, and therefore
\begin{align*}
\Phi(\rho_{N+1})\le \int_{B_R}|\nabla u|^p\,d\mathcal L^n<\infty.
\end{align*}
Since $0<\theta<1$, the terminal term $\theta^{N+1}\Phi(\rho_{N+1})$ tends to $0$. Since $\theta q^{-p}<1$, the geometric series converges. This is an ordinary limit of non-negative [real numbers](/page/Real%20Numbers) in the finite telescoping inequality, so letting $N\to\infty$ yields
\begin{align*}
\Phi(\rho_1)\le C_2A(R-r)^{-p},
\end{align*}
where
\begin{align*}
C_2:=\theta^{-1}(1-q)^{-p}\sum_{k=1}^{\infty}(\theta q^{-p})^k.
\end{align*}
Finally $B_r\subset B_{\rho_1}$, so
\begin{align*}
\int_{B_r}|\nabla u|^p\,d\mathcal L^n\le \Phi(\rho_1)\le \frac{C_2C_1}{(R-r)^p}\int_{B_R}|u-c|^p\,d\mathcal L^n+C_2C_1R^n.
\end{align*}
Renaming $C:=C_1C_2$ gives the asserted estimate with a constant depending only on $n,m,p,\lambda,\Lambda$.
[/guided]
[/step]