I'm trying to work through a derivation of the stability operator from minimal surface theory. Suppose $\Sigma^k$ is a minimal submanifold of $\mathbb{R}^n$, and suppose $X$ is a normal vector field on $\Sigma$ with compact support vanishing on the boundary.
Part of the derivation involves the integration by parts $\int_{\Sigma} \langle \nabla_{\Sigma}^N X, \nabla_{\Sigma}^N X \rangle = - \int_{\Sigma} \langle X, \Delta^N_{\Sigma} X \rangle$, where $\nabla_{\Sigma}^N$ is the normal projection of the covariant derivative on $\Sigma$, and $\Delta^N_{\Sigma}$ is the normal Laplacian defined by $\Sigma_{i = 1}^k \nabla_{E_i}^N \nabla_{E_i}^N X - \nabla^N_{\left(\nabla_{E^i} E_i\right)^T} X$, with $E_i$ being an orthonormal frame for $\Sigma$. Why does this hold? I'm aware of Green's identity holding for the Laplacian and gradient of scalar functions, but the operators involved here are normal projections, and I want to know why the identity still holds in this case.
I have tried to compute $\nabla_{\Sigma}^N (\nabla_{\Sigma}^N X)$ in hopes that its inner product with $X$ is the same as that with $\Delta_{\Sigma}^N X$, up to some terms which vanish under the integral like a divergence, but I haven't gotten anywhere.
The normal bundle along the submanifold in the Riemannian background is a special case of the following situation.
Let $M$ be a suitably good manifold of dimension $n$ with a Riemannian metric $g$. It has the Levi-Civita connection $\nabla = \nabla^g$ compatible with $g$, that is $\nabla g = 0$. Taking care with the orientation, there is also the Riemannian volume form $\mathrm{d}M$ on $M$.
Furthermore, let $E \to M$ be a (finite-dimesional) vector bundle over $M$ and forgive me if I sloppily denote the space of its sections by $E$ as well. Let $h: E \otimes E \to \mathbb{R}$ be a vector-bundle metric. Let $\nabla^E$ be some metric connection in $E$, that is $\nabla^E h = 0$.
There is a whole family of tensor products of the tangent, cotangent, and the given vector bundle $E$ on the manifold $E$, which are endowed with the corresponding tensor-product connection. More precisely, on each of such bundles, the connection acts by extending the Levi-Civita connection and the connection $\nabla^E$ using the Leibniz rule and linearity. Each of those connections we are going to denote simply by $\nabla$ and refer to them collectively as the coupled Levi-Civita connection.
The vector-bundle connection Laplacian $\Delta^E$ is defined on any section $X$ of the bundles, described above, by the formula $$ \Delta^E X := g^{a b} \nabla_a \nabla_b X $$ where $\nabla$ is the coupled Levi-Civita connection. The formula for the normal Laplacian in the OP is just a particular case of this.
Now I claim that $$ \big( \nabla^a h(X, \nabla_a X) \big) \mathrm{d}M = \mathrm{d}\omega $$ for some $(n - 1)$-form $\omega$ on $M$. This allows us to use the general Stokes theorem $$ \int_{M} \mathrm{d}\omega = \int_{\partial M} \omega $$ and, if we take into account the OP restrictions on the sections $X$ with regards to the boundary, we can say that $$ \int_{M} \big( \nabla^a h(X, \nabla_a X) \big) \mathrm{d}M = 0 $$
Using the Leibniz rule, we calculate $$ \nabla^a h(X, \nabla_a X) = h(\nabla^a X, \nabla_a X) + h(X, \nabla^a \nabla_a X) = |\nabla X|^2 + h(X, \Delta^E X) $$
From the last display, integrating and using the Stokes theorem, we obtain $$ \int_{M} |\nabla X|^2 \mathrm{d}M = - \int_M h(X, \Delta^E X) \mathrm{d}M $$
Proof of the claim. This is going to be a calculation using an abstract index notation technique for dealing with differential forms.
A $(n)$-form $\xi$ using the index notation can be presented as $\xi_{a_1 \dots a_n}$, where the sequentially numbered indexes are assumed to be antisymmetrized: $$ \xi_{a_1 \dots a_n} = \xi_{[a_1 \dots a_n]} $$
The exterior derivative $\mathrm{d}$ on the differential form $\xi$ as above is defined by the formula: $$ (\mathrm{d} \xi)_{a_0 a_1 \dots a_n} := (n + 1) \nabla_{a_0} \xi_{a_1 \dots a_n} $$
Recall that the volume form $\mathrm{d} M$ is a ${n}$-form, so we can write it as $(\mathrm{d} M)_{a_1 \dots a_n}$ to emphasize that.
Now we want $\xi$ to be a $1$-form (a covector field). We can introduce a notation, just for the convenience sake, for the interior product of the covector $\xi$ with the volume form $\mathrm{d} M$ as so: $$ (\iota_{\xi} \mathrm{d} M)_{a_2 \dots a_n} := g^{a a_1} \xi_a (\mathrm{d} M)_{a_1 \dots a_n} $$ Clearly, $\iota_{\xi} \mathrm{d} M$ is a $(n-1)$-form.
Let us compute the exterior derivative of $\iota_{\xi} \mathrm{d} M$ to see what happens. $$ \big( \mathrm{d} (\iota_{\xi} \mathrm{d} M) \big)_{a_0 a_2 \dots a_n} = n \nabla_{a_0 } (\iota_{\xi} \mathrm{d} M)_{a_2 \dots a_n} = n \nabla_{a_0 } g^{a a_1} \xi_a (\mathrm{d} M)_{a_1 \dots a_n} = \\ = n g^{a a_1} \big( \nabla_{a_0 } \xi_a (\mathrm{d} M)_{a_1 \dots a_n} \big) = n g^{a a_1} \big( (\nabla_{a_0 } \xi_a) (\mathrm{d} M)_{a_1 \dots a_n} + \xi_a \nabla_{a_0 } (\mathrm{d} M)_{a_1 \dots a_n} \big) = \\ = n g^{a a_1} \big( (\nabla_{a_0 } \xi_a) (\mathrm{d} M)_{a_1 \dots a_n} $$
The last expression in the display above is a $(n)$-form on $M$, so it must be proportional to $\mathrm{d} M$: $$ g^{a a_1} \big( (\nabla_{a_0 } \xi_a) (\mathrm{d} M)_{a_1 \dots a_n} = f (\mathrm{d} M)_{a_0 a_2 \dots a_n} $$
Contracting with $\mathrm{d} M$ and using its properties (see R.Wald, General Relativity, p.433), we can find that $$ \nabla^b \xi_b = n f $$
We see that for any covector $\xi$ the so-called divergence term $(\nabla^a \xi_a) \mathrm{d} M$ is an exact differential form: $$ (\nabla^a \xi_a) \mathrm{d} M = \mathrm{d}(\iota_{\xi} \mathrm{d} M) $$
Taking $\xi_a = h(X, \nabla_a X)$ we obtain the claimed property.