[proofplan] The proof converts constrained stationarity into an unconstrained Euclidean distributional equation on the Riemannian domain $(U,g)$. Weak harmonicity is understood as stationarity of the Dirichlet energy under compactly supported projected ambient variations $u_t=\Pi(u+t\psi)$, where $\Pi$ is a nearest-point projection to $N$. The first variation of these projected variations yields a weak identity containing the Hessian of $\Pi$, and the identity $D^2\Pi_p(X,Y)=A_p(X,Y)$ for tangent vectors converts that term into the second fundamental form. With the convention $A_p(X,Y)=(D_XY)^\perp$ and $\Delta_g u(\psi)=-\int_U(du,d\psi)_g\,d\operatorname{vol}_g$, the resulting distributional equation is $\Delta_g u-A(u)(du,du)_g=0$, and the reverse implication is obtained by reading the same computation backward for projected ambient variations. [/proofplan]

[step:Fix the Riemannian domain and choose a tubular projection] Let $(U,g)$ denote the smooth Riemannian domain in the statement. Let $d\operatorname{vol}_g$ denote its Riemannian volume measure, and define the distributional Laplace-Beltrami operator on an ambient map $u\in W^{1,2}(U;\mathbb{R}^q)$ by \begin{align*} \Delta_g u(\psi):=-\int_U (du,d\psi)_g\,d\operatorname{vol}_g \end{align*} for every $\psi\in C_c^\infty(U;\mathbb{R}^q)$. In this theorem, weak harmonicity means stationarity of the Dirichlet energy under every compactly supported projected ambient variation $u_t=\Pi(u+t\psi)$ generated by a test field $\psi\in C_c^\infty(U;\mathbb{R}^q)$. By the [Tubular Neighbourhood Theorem](/page/Tubular%20Neighbourhood%20Theorem) applied to the compact smooth embedding $N\subset\mathbb{R}^q$, there is an open tubular neighbourhood $\mathcal{T}\subset\mathbb{R}^q$ of $N$ and a smooth nearest-point projection \begin{align*} \Pi:\mathcal{T}\to N . \end{align*} For each $p\in N$, the derivative \begin{align*} D\Pi_p:\mathbb{R}^q\to T_pN \end{align*} is the [orthogonal projection](/theorems/437) onto $T_pN$. [/step]

[step:Differentiate projected ambient variations to obtain the weak projected identity]Let \begin{align*} \psi:U\to\mathbb{R}^q \end{align*} be a compactly supported smooth Euclidean test field. Since $\operatorname{supp}\psi$ is compact and $\psi$ is bounded, there exists $\varepsilon>0$ such that $u(x)+t\psi(x)\in\mathcal{T}$ for almost every $x\in U$ whenever $|t|<\varepsilon$. Define the projected variation $u_t:U\to N$ by \begin{align*} u_t(x)=\Pi(u(x)+t\psi(x)). \end{align*} The [Sobolev Chain Rule](/page/Sobolev%20Chain%20Rule) for smooth maps with bounded derivatives on the compact set reached by $u+t\psi$ gives $u_t\in W^{1,2}(U;N)$ and \begin{align*} \left.\frac{d}{dt}\right|_{t=0}u_t=D\Pi_u\psi \end{align*} in $W^{1,2}_0(U;\mathbb{R}^q)$. Define the energy functional $E:W^{1,2}(U;N)\to\mathbb{R}$ by \begin{align*} E[v]=\frac{1}{2}\int_U |dv|_g^2\,d\operatorname{vol}_g . \end{align*} The [First Variation Formula](/page/First%20Variation%20Formula) for projected Sobolev variations applies to this curve because $u_t\in W^{1,2}(U;N)$, $u_t=u$ outside the compact set $\operatorname{supp}\psi$, and the preceding chain-rule bounds imply differentiability of $t\mapsto u_t$ and $t\mapsto du_t$ at $t=0$ in the relevant $W^{1,2}$ and $L^2$ norms. The chain rule gives \begin{align*} du_t=D\Pi_{u+t\psi}(du+t\,d\psi). \end{align*} For $|t|$ small, the image of $u+t\psi$ lies in a compact subset $K_{\mathcal T}\subset\mathcal T$ on which $D\Pi$ and $D^2\Pi$ are bounded. Define \begin{align*} M_1&:=\sup_{y\in K_{\mathcal T}}\|D\Pi_y\|_{\mathrm{op}},& M_2&:=\sup_{y\in K_{\mathcal T}}\|D^2\Pi_y\|_{\mathrm{op}}. \end{align*} The chain rule gives \begin{align*} du_t=D\Pi_{u+t\psi}(du+t\,d\psi), \end{align*} so $|du_t|_g\leq M_1(|du|_g+|d\psi|_g)$. The difference quotient for $du_t$ decomposes into the term from differentiating $D\Pi_{u+t\psi}$ in the direction $\psi$, bounded in $L^2$ by $M_2\|\psi\|_{L^\infty}|du+t\,d\psi|_g$, and the term from differentiating $du+t\,d\psi$, bounded in $L^2$ by $M_1|d\psi|_g$. Thus the difference quotients are dominated in $L^2$ by \begin{align*} M_2\|\psi\|_{L^\infty}(|du|_g+|d\psi|_g)+M_1|d\psi|_g. \end{align*} Since $du\in L^2_{\mathrm{loc}}(U)$ and $d\psi$ is compactly supported and bounded, dominated convergence justifies differentiating $d(\Pi(u+t\psi))$ in $L^2$ and gives \begin{align*} 0 =\left.\frac{d}{dt}\right|_{t=0}E[u_t] =\int_U \bigl(du,d(D\Pi_u\psi)\bigr)_g\,d\operatorname{vol}_g . \end{align*} Writing the differential of $D\Pi_u\psi$ by the Sobolev chain rule yields \begin{align*} d(D\Pi_u\psi)=D^2\Pi_u(du,\psi)+D\Pi_u(d\psi). \end{align*} Because $du$ is tangent to $N$ almost everywhere and $D\Pi_u$ is the orthogonal projection onto $T_{u(x)}N$, we have \begin{align*} (du,D\Pi_u(d\psi))_g=(du,d\psi)_g. \end{align*} Therefore stationarity is equivalent to \begin{align*} \int_U (du,d\psi)_g\,d\operatorname{vol}_g + \int_U \bigl(du,D^2\Pi_u(du,\psi)\bigr)_g\,d\operatorname{vol}_g =0 \end{align*} for every $\psi\in C_c^\infty(U;\mathbb{R}^q)$.[/step]

[guided]The correct test fields for a Sobolev map are obtained from ambient smooth fields, not by assuming that there are many smooth fields tangent to the rough pullback bundle $u^{-1}TN$. Let \begin{align*} \psi:U\to\mathbb{R}^q \end{align*} be compactly supported and smooth. The projection $\Pi$ converts the ambient perturbation into an admissible constrained variation $u_t:U\to N$ defined by \begin{align*} u_t(x)=\Pi(u(x)+t\psi(x)). \end{align*} The compact support and boundedness of $\psi$ imply that $u+t\psi$ stays in the tubular neighbourhood $\mathcal T$ for all sufficiently small $|t|$, so $u_t$ is well-defined and takes values in $N$. Since $\Pi$ is smooth with bounded first and second derivatives on the relevant compact subset of $\mathcal T$, the [Sobolev Chain Rule](/page/Sobolev%20Chain%20Rule) gives $u_t\in W^{1,2}(U;N)$. More explicitly, \begin{align*} du_t=D\Pi_{u+t\psi}(du+t\,d\psi). \end{align*} Let $K_{\mathcal T}\subset\mathcal T$ be a compact subset containing the image of $u+t\psi$ for all sufficiently small $|t|$, and define \begin{align*} M_1&:=\sup_{y\in K_{\mathcal T}}\|D\Pi_y\|_{\mathrm{op}},& M_2&:=\sup_{y\in K_{\mathcal T}}\|D^2\Pi_y\|_{\mathrm{op}}. \end{align*} Then $|du_t|_g\leq M_1(|du|_g+|d\psi|_g)$. For the derivative in $t$, the difference quotient has two contributions: differentiating $D\Pi_{u+t\psi}$ in the direction $\psi$ gives a term controlled by $M_2\|\psi\|_{L^\infty}|du+t\,d\psi|_g$, and differentiating $du+t\,d\psi$ gives a term controlled by $M_1|d\psi|_g$. Hence the difference quotients are dominated in $L^2$ by \begin{align*} M_2\|\psi\|_{L^\infty}(|du|_g+|d\psi|_g)+M_1|d\psi|_g. \end{align*} This function belongs to $L^2(U)$ because $du\in L^2_{\mathrm{loc}}(U)$ and $d\psi$ is smooth with compact support. Therefore dominated convergence permits differentiation of $d(\Pi(u+t\psi))$ in $L^2$ at $t=0$. We use the energy functional $E:W^{1,2}(U;N)\to\mathbb{R}$ defined by \begin{align*} E[v]=\frac{1}{2}\int_U |dv|_g^2\,d\operatorname{vol}_g . \end{align*} Weak harmonicity, in the sense fixed in the theorem statement, means that $E$ is stationary under every compactly supported projected ambient variation through $N$. The present curve is admissible because $u_t\in W^{1,2}(U;N)$ and $u_t=u$ off $\operatorname{supp}\psi$. Applying the [First Variation Formula](/page/First%20Variation%20Formula) to this projected variation gives \begin{align*} 0 =\left.\frac{d}{dt}\right|_{t=0}E[u_t] =\int_U \bigl(du,d(D\Pi_u\psi)\bigr)_g\,d\operatorname{vol}_g . \end{align*} Now we expand the derivative of the tangent projection. The Sobolev chain rule gives \begin{align*} d(D\Pi_u\psi)=D^2\Pi_u(du,\psi)+D\Pi_u(d\psi). \end{align*} The second term simplifies because $du$ is tangent to $N$ almost everywhere and $D\Pi_u$ is the orthogonal projection onto $T_{u(x)}N$. Hence projecting the second factor does not change its [inner product](/page/Inner%20Product) with $du$: \begin{align*} (du,D\Pi_u(d\psi))_g=(du,d\psi)_g. \end{align*} Substitution gives the weak projected identity \begin{align*} \int_U (du,d\psi)_g\,d\operatorname{vol}_g + \int_U \bigl(du,D^2\Pi_u(du,\psi)\bigr)_g\,d\operatorname{vol}_g =0. \end{align*} This identity is the rigorous replacement for speaking about the tangential component of $\Delta_g u$ directly; it only pairs distributions with smooth ambient test fields.[/guided]

[step:Rewrite the projected identity as the extrinsic distributional equation] For $p\in N$ and tangent vectors $X,Y\in T_pN$, the nearest-point projection identity and the convention $A_p(X,Y)=(D_XY)^\perp$ give \begin{align*} D^2\Pi_p(X,Y)=A_p(X,Y). \end{align*} Indeed, choose a smooth tangent vector field $Y_0$ near $p$ with $Y_0(p)=Y$ and differentiate the identity $D\Pi(Y_0)=Y_0$ in the tangent direction $X$. The tangential part is $D\Pi_p(D_XY_0)$ and the normal part is $D_XY_0-D\Pi_p(D_XY_0)=A_p(X,Y)$. Since $du$ is tangent to $N$ almost everywhere, the preceding identity gives \begin{align*} D^2\Pi_u(du,du)_g=A(u)(du,du)_g. \end{align*} It remains to compute the sign in the mixed term. For $p\in N$, $X\in T_pN$, and $Z\in\mathbb{R}^q$, choose a smooth tangent vector field $X_0$ near $p$ with $X_0(p)=X$ and differentiate the scalar identity $X_0\cdot D\Pi Z=X_0\cdot Z$ in the tangent direction $X$. Since $D\Pi_p$ is self-adjoint and $D\Pi_pX=X$, this gives \begin{align*} X\cdot D^2\Pi_p(X,Z)=(D_XX_0-D\Pi_p(D_XX_0))\cdot Z=A_p(X,X)\cdot Z. \end{align*} Applying this pointwise with $p=u(x)$, $X=du(x)$, and $Z=\psi(x)$, then contracting with the metric $g$, yields \begin{align*} \bigl(du,D^2\Pi_u(du,\psi)\bigr)_g = A(u)(du,du)_g\cdot\psi \end{align*} for every $\psi\in C_c^\infty(U;\mathbb{R}^q)$. Substituting this into the weak projected identity gives \begin{align*} \int_U (du,d\psi)_g\,d\operatorname{vol}_g + \int_U A(u)(du,du)_g\cdot\psi\,d\operatorname{vol}_g =0. \end{align*} By the definition of the distributional Laplace-Beltrami operator, \begin{align*} \Delta_g u(\psi):=-\int_U (du,d\psi)_g\,d\operatorname{vol}_g, \end{align*} so the last identity is exactly \begin{align*} \Delta_g u-A(u)(du,du)_g=0 \end{align*} in $\mathcal{D}'(U;\mathbb{R}^q)$. [/step]

[step:Check that the nonlinear term defines a distribution] Since $A$ is smooth and $N$ is compact, there is a constant $C_A>0$ such that \begin{align*} |A_p(X,Y)|\leq C_A|X|\,|Y| \end{align*} for every $p\in N$ and every $X,Y\in T_pN$. Because $du\in L^2_{\mathrm{loc}}(U)$, this estimate gives $A(u)(du,du)_g\in L^1_{\mathrm{loc}}(U;\mathbb{R}^q)$. Therefore the distributional equation \begin{align*} \Delta_g u - A(u)(du,du)_g = 0 \end{align*} in $\mathcal{D}'(U;\mathbb{R}^q)$ is well-defined. [/step]

[step:Recover weak harmonicity from the distributional equation] Conversely, assume \begin{align*} \Delta_g u - A(u)(du,du)_g = 0 \end{align*} in $\mathcal{D}'(U;\mathbb{R}^q)$. Let \begin{align*} \psi:U\to\mathbb{R}^q \end{align*} be any compactly supported smooth Euclidean test field, and define $u_t=\Pi(u+t\psi)$ for sufficiently small $|t|$. Testing the distributional equation against $\psi$ gives \begin{align*} \int_U (du,d\psi)_g\,d\operatorname{vol}_g + \int_U A(u)(du,du)_g\cdot\psi\,d\operatorname{vol}_g =0. \end{align*} Using again $D^2\Pi_p(X,Y)=A_p(X,Y)$ for tangent $X,Y$, the mixed identity \begin{align*} \bigl(du,D^2\Pi_u(du,\psi)\bigr)_g = A(u)(du,du)_g\cdot\psi, \end{align*} and the [Sobolev Chain Rule](/page/Sobolev%20Chain%20Rule) expansion of $d(D\Pi_u\psi)$, this identity is equivalent to \begin{align*} \int_U \bigl(du,d(D\Pi_u\psi)\bigr)_g\,d\operatorname{vol}_g=0. \end{align*} The [First Variation Formula](/page/First%20Variation%20Formula) for the admissible projected variation $u_t=\Pi(u+t\psi)$ therefore gives \begin{align*} \left.\frac{d}{dt}\right|_{t=0}E[u_t]=0. \end{align*} Since $\psi\in C_c^\infty(U;\mathbb{R}^q)$ was arbitrary, the energy is stationary under every compactly supported projected ambient variation. By the definition of weak harmonicity used in the theorem statement, $u$ is weakly harmonic. This proves the equivalence. [/step]