[proofplan] The proof compares the variational equation against tangent test fields with the full distributional equation in $\mathbb{R}^q$. The key geometric input is that the Euclidean second derivative of a map into $N$ splits into tangent and normal parts, and the normal part is exactly the second fundamental form applied to first derivatives. One direction follows by testing the distributional equation against tangent fields, because $A_u(du,du)$ is normal. For the converse, arbitrary smooth ambient variations are projected back to $N$ by the nearest-point projection; differentiating the energy of the projected variation at the variation parameter, not differentiating $u$ twice, produces the distributional identity. [/proofplan] [step:Fix the projected weak formulation and the distributional notation] We use the extrinsic weak formulation by projected ambient variations: $u$ is weakly harmonic if, for every $\phi\in C_c^\infty(U;\mathbb{R}^q)$ and every sufficiently small projected variation $u_t:U\to N$ of the form \begin{align*} u_t(x)=\pi(u(x)+t\phi(x)), \end{align*} the first variation of the Dirichlet energy satisfies \begin{align*} \left.\frac{d}{dt}\right|_{t=0}\frac{1}{2}\int_U g^{\alpha\beta}\,\partial_{x_\alpha}u_t\cdot \partial_{x_\beta}u_t\,d\operatorname{vol}_g(x)=0. \end{align*} Here $\pi$ denotes the nearest-point projection on a tubular neighbourhood of $N$, constructed below. This avoids testing distributions against a field of the form $d\pi_u\phi$, which need not be a smooth $\mathbb{R}^q$-valued [test function](/page/Test%20Function) when $u$ is only Sobolev. The distributional Laplacian $\Delta_g u\in\mathcal{D}'(U;\mathbb{R}^q)$ is defined by \begin{align*} \Delta_g u(\phi) := -\int_U g^{\alpha\beta}\,\partial_{x_\alpha}u\cdot \partial_{x_\beta}\phi\,d\operatorname{vol}_g(x) \end{align*} for every $\phi\in C_c^\infty(U;\mathbb{R}^q)$. In local coordinates $(x_1,\dots,x_m)$ on $U$, define the $L^1_{\mathrm{loc}}(U;\mathbb{R}^q)$ vector field $A_u(du,du):U\to\mathbb{R}^q$ by \begin{align*} A_u(du,du)(x) := g^{\alpha\beta}(x)A_{u(x)}\bigl(\partial_{x_\alpha}u(x),\partial_{x_\beta}u(x)\bigr) \end{align*} for almost every $x\in U$. This is locally integrable because $A$ is bounded on the compact submanifold $N$ and $du\in L^2_{\mathrm{loc}}(U;T^*U\otimes\mathbb{R}^q)$. [/step] [step:Test the distributional equation only against smooth ambient fields] Assume \begin{align*} \Delta_g u + A_u(du,du)=0 \end{align*} in $\mathcal{D}'(U;\mathbb{R}^q)$, and let $\phi\in C_c^\infty(U;\mathbb{R}^q)$ be arbitrary. Using the nearest-point projection $\pi$ from the tubular neighbourhood construction below, define $u_t:U\to N$ by \begin{align*} u_t(x)=\pi(u(x)+t\phi(x)) \end{align*} for sufficiently small $|t|$. The differentiation formula justified in the next step gives \begin{align*} \left.\frac{d}{dt}\right|_{t=0}\frac{1}{2}\int_U g^{\alpha\beta}\,\partial_{x_\alpha}u_t\cdot\partial_{x_\beta}u_t\,d\operatorname{vol}_g(x) = \int_U g^{\alpha\beta}\,\partial_{x_\alpha}u\cdot\partial_{x_\beta}\phi\,d\operatorname{vol}_g(x) - \int_U A_u(du,du)\cdot\phi\,d\operatorname{vol}_g(x). \end{align*} The assumed distributional equation tested against the smooth ambient field $\phi$ says \begin{align*} 0 = -\int_U g^{\alpha\beta}\,\partial_{x_\alpha}u\cdot\partial_{x_\beta}\phi\,d\operatorname{vol}_g(x) + \int_U A_u(du,du)\cdot\phi\,d\operatorname{vol}_g(x). \end{align*} Hence the first variation above vanishes for every $\phi\in C_c^\infty(U;\mathbb{R}^q)$, so $u$ is weakly harmonic in the projected-variation sense. [/step] [step:Project smooth ambient variations back to $N$ and differentiate the energy] Assume now that $u$ is weakly harmonic, and let $\phi\in C_c^\infty(U;\mathbb{R}^q)$ be arbitrary. Since $N$ is compact and embedded, the [Tubular Neighbourhood Theorem](/theorems/???) gives an open 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 every $p\in N$, the differential $d\pi_p:\mathbb{R}^q\to T_pN$ is the [orthogonal projection](/theorems/437) onto $T_pN$. Since $u(U)\subset N$ almost everywhere and $\operatorname{supp}\phi$ is compact, there exists $\varepsilon>0$ such that $u(x)+t\phi(x)\in\mathcal{T}$ for almost every $x\in U$ and every $t\in(-\varepsilon,\varepsilon)$. For $t\in(-\varepsilon,\varepsilon)$ define the projected variation $u_t:U\to N$ by \begin{align*} u_t(x)=\pi(u(x)+t\phi(x)). \end{align*} The Sobolev chain rule applies because $\pi$ is smooth with bounded first and second derivatives on the compact set swept out by $u(x)+t\phi(x)$. Hence, for almost every $x\in U$, \begin{align*} \left.\frac{d}{dt}\right|_{t=0}u_t(x)=d\pi_{u(x)}(\phi(x)) \end{align*} and \begin{align*} \left.\frac{d}{dt}\right|_{t=0}\partial_{x_\alpha}u_t(x) = d\pi_{u(x)}(\partial_{x_\alpha}\phi(x)) + d^2\pi_{u(x)}(\partial_{x_\alpha}u(x),\phi(x)). \end{align*} Let $K\subset U$ be a compact set containing $\operatorname{supp}\phi$, and let $M_1,M_2>0$ be bounds for $d\pi$ and $d^2\pi$ on the compact subset of $\mathcal{T}$ reached by $u(x)+t\phi(x)$. The difference quotients of $\partial_{x_\alpha}u_t$ at $t=0$ are bounded in $L^2(K)$ by a constant multiple of $|\partial_{x_\alpha}\phi|+|\partial_{x_\alpha}u|\,|\phi|$. Therefore the first-variation integrands are dominated in $L^1(K,d\operatorname{vol}_g)$ by a constant multiple of $|du|\,|d\phi|+|du|^2|\phi|$, which is integrable because $du\in L^2_{\mathrm{loc}}$ and $\phi,d\phi$ are bounded with compact support. The [dominated convergence theorem](/theorems/4) justifies differentiating the energy under the integral sign. Weak harmonicity of $u$ gives the vanishing of the first variation of the Dirichlet energy along $u_t$: \begin{align*} 0 = \int_U g^{\alpha\beta}\,\partial_{x_\alpha}u\cdot\left(d\pi_u(\partial_{x_\beta}\phi)+d^2\pi_u(\partial_{x_\beta}u,\phi)\right)\,d\operatorname{vol}_g. \end{align*} Because $\partial_{x_\alpha}u(x)\in T_{u(x)}N$ almost everywhere and $d\pi_{u(x)}$ is the orthogonal projection onto $T_{u(x)}N$, self-adjointness of this projection gives \begin{align*} \partial_{x_\alpha}u\cdot d\pi_u(\partial_{x_\beta}\phi) = \partial_{x_\alpha}u\cdot \partial_{x_\beta}\phi. \end{align*} For the second term, fix $p\in N$, tangent vectors $X,Z\in T_pN$, and an ambient vector $V\in\mathbb{R}^q$. Differentiating the identity $d\pi_p(Z)=Z$ in the ambient direction $V$ gives the standard Weingarten relation \begin{align*} Z\cdot d^2\pi_p(X,V)=-A_p(Z,X)\cdot V. \end{align*} Indeed, the tangent component of $V$ contributes no normal pairing, while the normal component differentiates the tangent projection by the shape operator associated to that normal vector. Applying this with $p=u(x)$, $Z=\partial_{x_\alpha}u(x)$, $X=\partial_{x_\beta}u(x)$, and $V=\phi(x)$ yields \begin{align*} \partial_{x_\alpha}u\cdot d^2\pi_u(\partial_{x_\beta}u,\phi) = -A_u(\partial_{x_\alpha}u,\partial_{x_\beta}u)\cdot\phi. \end{align*} Substituting these two identities into the [first variation formula](/theorems/2728) yields \begin{align*} 0 = \int_U g^{\alpha\beta}\,\partial_{x_\alpha}u\cdot\partial_{x_\beta}\phi\,d\operatorname{vol}_g - \int_U A_u(du,du)\cdot\phi\,d\operatorname{vol}_g. \end{align*} Therefore \begin{align*} \Delta_g u(\phi) = -\int_U A_u(du,du)\cdot\phi\,d\operatorname{vol}_g. \end{align*} [guided] The obstruction in the weak setting is that the pointwise tangent projection $d\pi_{u(x)}(\phi(x))$ need not be a smooth tangent field along the Sobolev map $u$. We avoid testing directly against that projected field. Instead, we build an actual admissible variation through maps into $N$ and differentiate the energy with respect to the scalar variation parameter. Since $N$ is compact and embedded, the [Tubular Neighbourhood Theorem](/theorems/???) gives an open 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 every $p\in N$, the differential $d\pi_p:\mathbb{R}^q\to T_pN$ is the orthogonal projection onto $T_pN$. Because $\phi$ has compact support and $N$ is compact, the set of points $u(x)+t\phi(x)$ remains in $\mathcal{T}$ for all sufficiently small $|t|$; choose $\varepsilon>0$ with this property. For $t\in(-\varepsilon,\varepsilon)$ define \begin{align*} u_t:U&\to N \end{align*} by \begin{align*} u_t(x)&=\pi(u(x)+t\phi(x)). \end{align*} The map $u_t$ is a Sobolev map into $N$ because $u\in W^{1,2}(U;N)$, $\phi$ is smooth and compactly supported, and composition with the smooth map $\pi$ is valid on the compact subset of $\mathcal{T}$ reached by the variation. The Sobolev chain rule gives, for almost every $x\in U$, \begin{align*} \left.\frac{d}{dt}\right|_{t=0}u_t(x)=d\pi_{u(x)}(\phi(x)) \end{align*} and \begin{align*} \left.\frac{d}{dt}\right|_{t=0}\partial_{x_\alpha}u_t(x) = d\pi_{u(x)}(\partial_{x_\alpha}\phi(x)) + d^2\pi_{u(x)}(\partial_{x_\alpha}u(x),\phi(x)). \end{align*} This is the replacement for the invalid computation with second derivatives of $u$: all derivatives fall either on the smooth variation field $\phi$ or on the smooth projection $\pi$. We must also justify moving the derivative in $t$ through the energy integral. Let $K\subset U$ be a compact set containing $\operatorname{supp}\phi$, and choose constants $M_1,M_2>0$ bounding $d\pi$ and $d^2\pi$ on the compact subset reached by $u(x)+t\phi(x)$. The chain-rule formula bounds the difference quotients of $\partial_{x_\alpha}u_t$ by a constant multiple of $|\partial_{x_\alpha}\phi|+|\partial_{x_\alpha}u|\,|\phi|$. After pairing with $\partial_{x_\alpha}u$, the first-variation integrand is dominated by a constant multiple of $|du|\,|d\phi|+|du|^2|\phi|$. This function lies in $L^1(K,d\operatorname{vol}_g)$ because $du\in L^2_{\mathrm{loc}}$ and $\phi,d\phi$ are bounded with compact support. Therefore the dominated convergence theorem permits differentiating the Dirichlet energy under the integral sign. Weak harmonicity says that the first derivative at $t=0$ of the Dirichlet energy of this projected variation vanishes. Thus \begin{align*} 0 = \int_U g^{\alpha\beta}\,\partial_{x_\alpha}u\cdot\left(d\pi_u(\partial_{x_\beta}\phi)+d^2\pi_u(\partial_{x_\beta}u,\phi)\right)\,d\operatorname{vol}_g. \end{align*} Since $\partial_{x_\alpha}u(x)$ is tangent to $N$ at $u(x)$ for almost every $x$, and since $d\pi_{u(x)}$ is the orthogonal projection onto $T_{u(x)}N$, we may remove the projection from the first term after pairing with $\partial_{x_\alpha}u$: \begin{align*} \partial_{x_\alpha}u\cdot d\pi_u(\partial_{x_\beta}\phi) = \partial_{x_\alpha}u\cdot \partial_{x_\beta}\phi. \end{align*} The second derivative of the nearest-point projection records the curvature of the embedded submanifold. With the convention $A_p(X,Y)=(D_XY)^\perp$, fix $p\in N$, $X,Z\in T_pN$, and $V\in\mathbb{R}^q$. Differentiating the projection identity $d\pi_p(Z)=Z$ in the ambient direction $V$ and using the Weingarten relation for the normal component of $V$ gives \begin{align*} Z\cdot d^2\pi_p(X,V)=-A_p(Z,X)\cdot V. \end{align*} The tangent component of $V$ gives zero on the right because $A_p(Z,X)$ is normal, and it gives zero on the left after pairing with $Z$ because $d^2\pi_p(X,V^\top)$ is normal. Applying the identity with $p=u(x)$, $Z=\partial_{x_\alpha}u(x)$, $X=\partial_{x_\beta}u(x)$, and $V=\phi(x)$ gives \begin{align*} \partial_{x_\alpha}u\cdot d^2\pi_u(\partial_{x_\beta}u,\phi) = -A_u(\partial_{x_\alpha}u,\partial_{x_\beta}u)\cdot\phi. \end{align*} Substitution gives \begin{align*} 0 = \int_U g^{\alpha\beta}\,\partial_{x_\alpha}u\cdot\partial_{x_\beta}\phi\,d\operatorname{vol}_g - \int_U A_u(du,du)\cdot\phi\,d\operatorname{vol}_g. \end{align*} Using the definition \begin{align*} \Delta_g u(\phi) = -\int_U g^{\alpha\beta}\,\partial_{x_\alpha}u\cdot\partial_{x_\beta}\phi\,d\operatorname{vol}_g, \end{align*} we obtain \begin{align*} \Delta_g u(\phi) = -\int_U A_u(du,du)\cdot\phi\,d\operatorname{vol}_g. \end{align*} This is exactly the desired distributional identity when the term $A_u(du,du)$ is interpreted as the [regular distribution](/page/Regular%20Distribution) defined by integration against $d\operatorname{vol}_g$. [/guided] [/step] [step:Conclude the distributional equation] The previous step proves that every $\phi\in C_c^\infty(U;\mathbb{R}^q)$ satisfies \begin{align*} \Delta_g u(\phi) + \int_U A_u(du,du)\cdot\phi\,d\operatorname{vol}_g =0. \end{align*} Equivalently, \begin{align*} \bigl(\Delta_g u + A_u(du,du)\bigr)(\phi)=0, \end{align*} where the regular distribution associated to $A_u(du,du)$ acts by \begin{align*} A_u(du,du)(\phi) = \int_U A_u(du,du)\cdot\phi\,d\operatorname{vol}_g. \end{align*} Because $\phi\in C_c^\infty(U;\mathbb{R}^q)$ was arbitrary, the identity holds in $\mathcal{D}'(U;\mathbb{R}^q)$. [/step]