[proofplan]
We rescale to the unit ball, because the normalized Dirichlet energy is invariant under the dilation $x=x_0+ry$. The small-energy hypothesis gives a Campanato decay estimate for energy-minimizing manifold-valued maps; iterating that estimate yields Morrey growth for the energy on every smaller ball. Morrey's Dirichlet growth theorem converts this decay into Hölder continuity, and the extrinsic harmonic map equation then gives smoothness and the stated interior gradient estimate after scaling back.
[/proofplan]
[step:Rescale the map to the unit ball]
Define the rescaled Sobolev map $v: B(0,2)\to N$ by
\begin{align*}
v(y):=u(x_0+ry).
\end{align*}
Since $u\in W^{1,2}(B(x_0,2r);N)$, the Sobolev chain rule for affine changes of variables gives $v\in W^{1,2}(B(0,2);N)$ and $dv(y)=r\,du(x_0+ry)$ for $\mathcal L^m$-a.e. $y\in B(0,2)$. Let $B(a,\rho)\subset B(0,2)$ be a ball, and let $\widetilde v\in W^{1,2}(B(a,\rho);N)$ be a competitor for $v$ with the same trace as $v$ on $\partial B(a,\rho)$. Define the corresponding competitor $\widetilde u: B(x_0+ra,r\rho)\to N$ by
\begin{align*}
\widetilde u(x):=\widetilde v\left(\frac{x-x_0}{r}\right).
\end{align*}
This affine pullback identifies competitors for $v$ on $B(a,\rho)$ with competitors for $u$ on $B(x_0+ra,r\rho)$, and the same change-of-variables computation rescales both Dirichlet energies by the positive factor $r^{m-2}$. Since $u$ is energy-minimizing, $v$ is energy-minimizing on $B(0,2)$.
Using the change of variables $x=x_0+ry$, so that $d\mathcal L^m(x)=r^m\,d\mathcal L^m(y)$, we first have
\begin{align*}
\int_{B(0,2)} |dv(y)|^2\,d\mathcal L^m(y)=\int_{B(0,2)} r^2 |du(x_0+ry)|^2\,d\mathcal L^m(y).
\end{align*}
The same substitution sends $B(0,2)$ to $B(x_0,2r)$, hence
\begin{align*}
\int_{B(0,2)} |dv(y)|^2\,d\mathcal L^m(y)=r^{2-m}\int_{B(x_0,2r)} |du(x)|^2\,d\mathcal L^m(x)<\varepsilon_0.
\end{align*}
[/step]
[step:Apply small-energy decay to obtain Morrey growth]
Let $E_v(a,\rho)$ denote the normalized energy
\begin{align*}
E_v(a,\rho):=\rho^{2-m}\int_{B(a,\rho)} |dv(y)|^2\,d\mathcal L^m(y)
\end{align*}
for $a\in B(0,3/2)$ and $0<\rho\le 1/2$. Choose constants $\varepsilon_*>0$, $\theta\in(0,1/4)$, and $\gamma\in(0,1)$ from the small-energy decay theorem for energy-minimizing harmonic maps into the compact embedded manifold $N$: if $z\in B(0,2)$, $B(z,s)\subset B(0,2)$, and $E_v(z,s)\le \varepsilon_*$, then
\begin{align*}
E_v(z,\theta s)\le \theta^{2\gamma}E_v(z,s).
\end{align*}
We choose the theorem's constant $\varepsilon_0$ small enough below; in particular require $\varepsilon_0\le 2^{2-m}\varepsilon_*$. For every $a\in B(0,3/2)$, the inclusion $B(a,1/2)\subset B(0,2)$ gives
\begin{align*}
E_v(a,1/2)=2^{m-2}\int_{B(a,1/2)} |dv(y)|^2\,d\mathcal L^m(y)\le 2^{m-2}\int_{B(0,2)} |dv(y)|^2\,d\mathcal L^m(y)<\varepsilon_*.
\end{align*}
Thus the decay theorem applies first on $B(a,1/2)$ and then, by induction, on each smaller concentric ball $B(a,\theta^j/2)$.
Iterating the decay inequality gives, for every integer $k\ge0$,
\begin{align*}
E_v(a,\theta^k/2)\le \theta^{2\gamma k}E_v(a,1/2).
\end{align*}
If $0<\rho\le1/2$, choose $k\in\mathbb N\cup\{0\}$ such that $\theta^{k+1}/2<\rho\le\theta^k/2$. Since $B(a,\rho)\subset B(a,\theta^k/2)$, domain monotonicity for nonnegative integrals gives
\begin{align*}
\int_{B(a,\rho)} |dv|^2\,d\mathcal L^m\le \int_{B(a,\theta^k/2)} |dv|^2\,d\mathcal L^m.
\end{align*}
By the definition of normalized energy and the iterated decay estimate,
\begin{align*}
\int_{B(a,\theta^k/2)} |dv|^2\,d\mathcal L^m\le (\theta^k/2)^{m-2+2\gamma}E_v(a,1/2).
\end{align*}
Because $\theta^{k+1}/2<\rho$, we have $\theta^k/2<\theta^{-1}\rho$, and therefore
\begin{align*}
\int_{B(a,\rho)} |dv|^2\,d\mathcal L^m\le \theta^{-(m-2+2\gamma)}\rho^{m-2+2\gamma}E_v(a,1/2).
\end{align*}
Using $E_v(a,1/2)\le 2^{m-2}\int_{B(0,2)} |dv|^2\,d\mathcal L^m$, defining $C_1:=2^{m-2}\theta^{-(m-2+2\gamma)}$ yields
\begin{align*}
\int_{B(a,\rho)} |dv|^2\,d\mathcal L^m\le C_1\rho^{m-2+2\gamma}\int_{B(0,2)} |dv|^2\,d\mathcal L^m
\end{align*}
for all $a\in B(0,3/2)$ and $0<\rho\le1/2$.
We also require $\varepsilon_0\le \eta_0/C_1$, where $\eta_0=\eta_0(m,N,\gamma)$ is the small Morrey-norm threshold in the interior bootstrap lemma used below.
[guided]
The purpose of this step is to turn small total energy into energy decay at every smaller scale on a domain that still contains a neighbourhood of $\overline{B(0,1)}$. Define
\begin{align*}
E_v(a,\rho):=\rho^{2-m}\int_{B(a,\rho)} |dv(y)|^2\,d\mathcal L^m(y),
\end{align*}
where $a\in B(0,3/2)$ and $0<\rho\le1/2$. This is the scale-invariant quantity for the Dirichlet energy in dimension $m$.
We use the small-energy decay theorem for energy-minimizing harmonic maps. Its hypotheses are satisfied on $B(a,1/2)$ because $v\in W^{1,2}(B(0,2);N)$ is energy-minimizing, $B(a,1/2)\subset B(0,2)$ for $a\in B(0,3/2)$, and
\begin{align*}
E_v(a,1/2)=2^{m-2}\int_{B(a,1/2)} |dv(y)|^2\,d\mathcal L^m(y)\le2^{m-2}\int_{B(0,2)} |dv(y)|^2\,d\mathcal L^m(y)<\varepsilon_*.
\end{align*}
This is why we choose $\varepsilon_0:=2^{2-m}\varepsilon_*$: it compensates for the factor in the normalized energy at radius $1/2$.
The theorem gives
\begin{align*}
E_v(a,\theta/2)\le\theta^{2\gamma}E_v(a,1/2).
\end{align*}
Since the right-hand side is at most $E_v(a,1/2)<\varepsilon_*$, the smallness hypothesis is available again on $B(a,\theta/2)$. Repeating this induction gives
\begin{align*}
E_v(a,\theta^k/2)\le\theta^{2\gamma k}E_v(a,1/2)
\end{align*}
for every integer $k\ge0$.
Now let $0<\rho\le1/2$ and choose $k\in\mathbb N\cup\{0\}$ such that $\theta^{k+1}/2<\rho\le\theta^k/2$. The inclusion $B(a,\rho)\subset B(a,\theta^k/2)$ and nonnegativity of $|dv|^2$ give
\begin{align*}
\int_{B(a,\rho)} |dv|^2\,d\mathcal L^m\le\int_{B(a,\theta^k/2)} |dv|^2\,d\mathcal L^m.
\end{align*}
By the definition of $E_v(a,\theta^k/2)$ and the iterated decay estimate,
\begin{align*}
\int_{B(a,\theta^k/2)} |dv|^2\,d\mathcal L^m\le(\theta^k/2)^{m-2+2\gamma}E_v(a,1/2).
\end{align*}
Because $\rho>\theta^{k+1}/2$, we have $\theta^k/2<\theta^{-1}\rho$, so
\begin{align*}
(\theta^k/2)^{m-2+2\gamma}\le\theta^{-(m-2+2\gamma)}\rho^{m-2+2\gamma}.
\end{align*}
Finally $E_v(a,1/2)\le2^{m-2}\int_{B(0,2)} |dv|^2\,d\mathcal L^m$, hence
\begin{align*}
\int_{B(a,\rho)} |dv|^2\,d\mathcal L^m\le C_1\rho^{m-2+2\gamma}\int_{B(0,2)} |dv|^2\,d\mathcal L^m,
\end{align*}
where $C_1:=2^{m-2}\theta^{-(m-2+2\gamma)}$ depends only on $m$ and $N$ through the constants in the decay theorem.
[/guided]
[/step]
[step:Convert Morrey growth into Hölder continuity]
We apply Morrey's Dirichlet growth theorem on the open ball $B(0,3/2)$. Its local hypothesis requires constants $M>0$ and $\gamma\in(0,1)$ such that, for every $a\in B(0,3/2)$ and every radius $\rho$ with $0<\rho\le \min\{1/2,\operatorname{dist}(a,\partial B(0,3/2))\}$, one has
\begin{align*}
\int_{B(a,\rho)} |dv|^2\,d\mathcal L^m\le M\rho^{m-2+2\gamma}.
\end{align*}
The previous step verifies this hypothesis with
\begin{align*}
M:=C_1\int_{B(0,2)} |dv|^2\,d\mathcal L^m.
\end{align*}
Therefore Morrey's theorem gives a representative, still denoted $v$, and a constant $C_2=C_2(m,N)>0$ such that
\begin{align*}
|v(y)-v(z)|\le C_2 |y-z|^\gamma\left(\int_{B(0,2)} |dv|^2\,d\mathcal L^m\right)^{1/2}
\end{align*}
for all $y,z\in B(0,3/2)$. In particular $v$ is Hölder continuous on $B(0,3/2)$.
[/step]
[step:Bootstrap the extrinsic harmonic map equation]
Let $A: N\to \operatorname{Bil}(\mathbb R^q\times\mathbb R^q,\mathbb R^q)$ denote the second fundamental form of the embedded submanifold $N\subset\mathbb R^q$, so $A(p)=A_p$ for each $p\in N$.
Because $v$ is energy-minimizing, it is a weak harmonic map and satisfies the extrinsic Euler-Lagrange system
\begin{align*}
\Delta v=-A(v)(\nabla v,\nabla v)
\end{align*}
in the sense of distributions on $B(0,1)$. Let
\begin{align*}
E:=\int_{B(0,2)} |dv|^2\,d\mathcal L^m.
\end{align*}
The Morrey growth estimate from the previous step implies
\begin{align*}
\rho^{2-m}\int_{B(a,\rho)}|dv|^2\,d\mathcal L^m \le C_1\rho^{2\gamma}E\le C_1\varepsilon_0\le \eta_0
\end{align*}
for every $a\in B(0,3/2)$ and every $0<\rho\le 1/2$ with $B(a,\rho)\subset B(0,3/2)$.
We now use the following standard interior bootstrap lemma of Hildebrandt-Kaul-Widman type for weak harmonic maps. For fixed $m$, compact embedded target $N$, and exponent $\gamma\in(0,1)$, there are constants $\eta_0>0$ and $C_3<\infty$ such that if $w\in W^{1,2}(B(0,3/2);N)$ is a weak solution of $\Delta w=-A(w)(\nabla w,\nabla w)$, is Hölder continuous with exponent $\gamma$, and satisfies the preceding Morrey bound with norm at most $\eta_0$, then $w$ is smooth on $B(0,1)$ and
\begin{align*}
\sup_{B(0,1)} |dw|^2\le C_3\int_{B(0,3/2)} |dw|^2\,d\mathcal L^m.
\end{align*}
This is an elliptic [regularity theorem](/theorems/2750) for the semilinear system after Morrey decay and Hölder continuity have already been established; it does not assert that small total energy alone gives regularity. Its hypotheses hold for $w=v$ by energy minimizing, the displayed Morrey bound, and Morrey continuity. Since $B(0,3/2)\subset B(0,2)$, the estimate gives
\begin{align*}
\sup_{B(0,1)} |dv|^2\le C_3\int_{B(0,2)} |dv|^2\,d\mathcal L^m,
\end{align*}
after increasing $C_3$ if necessary.
[/step]
[step:Scale the smoothness and gradient bound back to $B(x_0,r)$]
For $x\in B(x_0,r)$ write $y=(x-x_0)/r\in B(0,1)$. Since $u(x)=v(y)$ and $du(x)=r^{-1}dv(y)$, taking suprema over corresponding points gives
\begin{align*}
\sup_{B(x_0,r)} |du|^2=r^{-2}\sup_{B(0,1)} |dv|^2.
\end{align*}
Using the scale-one gradient estimate for $v$ and then the change of variables $x=x_0+ry$ gives
\begin{align*}
\sup_{B(x_0,r)} |du|^2\le C_3 r^{-2}\int_{B(0,2)} |dv|^2\,d\mathcal L^m.
\end{align*}
The rescaling identity for the Dirichlet energy is
\begin{align*}
\int_{B(0,2)} |dv(y)|^2\,d\mathcal L^m(y)=r^{2-m}\int_{B(x_0,2r)} |du(x)|^2\,d\mathcal L^m(x).
\end{align*}
Therefore
\begin{align*}
\sup_{B(x_0,r)} |du|^2\le C_3 r^{-m}\int_{B(x_0,2r)} |du|^2\,d\mathcal L^m.
\end{align*}
Set $C:=C_3$. Since $v$ is smooth on $B(0,1)$ and the dilation $y\mapsto x_0+ry$ is smooth with smooth inverse, $u$ is smooth on $B(x_0,r)$. This proves the theorem.
[/step]