Say that $X:M \to \mathbb R^3$ is an isometric immersion of an oriented riemannian surface (oriented $2$ dimensional riemannian manifold). I understand there holds a vector-valued equation, namely
$$ \Delta X = 2HN, $$
where $H$ is the mean curvature (half the trace of the second fundamental form) and $N$ is the outward unit normal to $X(M)$. I've been unable to find a reference, however.
I would appreciate a reference or an indication of how to compute this formula.
On
Based on my study of @Ernie060's answer, I believe I've worked out a solution.
Fix an isometric immersion $X: M \to \mathbb R^3$. For some arbitrary point $p \in X(M)$, consider an orthonormal basis $e_1, e_2$ of $T_p X(M)$. Using existence theory for systems of linear ODE, in a small $X(M)$-neighborhood of $p$ let us extend $e_1, e_2$ to an orthonormal frame $E_1, E_2$ satisfying $\nabla_{E_i} E_j = 0$ for $1 \leq i,j \leq 2$. Here the covariant derivative is defined using the formula $$D_{E_i} E_j = \nabla_{E_i} E_j + (D_{E_i}E_j \cdot n)n,$$ where $n$ is the outward normal to $X(M)$. It follows that $E_1, E_2, n$ is an orthonormal frame for $\mathbb R^3$ in a small ball around $p$. In this context it is appropriate to regard $X$ as the inclusion $X(M) \to \mathbb R^3$. Extend $X$ to a ball around $p$ via the rule $X(p + tn) = X(p) = p$, so that $D_n X = 0$. We wish to calculate $$ \Delta X,$$ where $\Delta X = \mathrm{div}(\mathrm{grad} X)$, the Euclidean vector laplacian. Here we write $$\mathrm{grad} X = (D_{E_1} X) \otimes dE_1 + (D_{E_2} X)\otimes dE_2 + (D_{n} X)\otimes dn,$$ where $dE_i$ and $dn$ are the one-forms dual to $E_i$ and $n$. A quick calculation shows that $D_{E_i} X = E_i$ so that $$\mathrm{grad} X= E_1\otimes dE_1 + E_2\otimes dE_2.$$ Then $$\Delta X =\mathrm{div} (\mathrm {grad} X) = \mathrm{trace}(D\mathrm (grad X)) = D_{E_1} E_1 + D_{E_2}E_2.$$ As $\nabla_{E_i} E_j = 0$ we see that $D_{E_i} E_i = (D_{E_i} E_i \cdot n) n.$ Differentiating the equation $E_i \cdot n = 0$ gives us $$D_{E_i}E_i\cdot n = -D_{E_i}n\cdot E_i = II(E_i, E_i),$$ where $II$ is the second fundamental form of $X$. Then $$ D_{E_1} E_1 + D_{E_2}E_2 = (II(E_1, E_1) + II(E_2,E_2))n = 2Hn,$$ as desired.
One reference is Pseudo-Riemannian Geometry and Delta Invariants by Bang-Yen Chen. On page 41, the so-called formula of Beltrami is derived for any immersion $\mathbf{x}\colon M \to E^m_s$ of a pseudo-Riemannian $n$-manifold into the pseudo-Euclidean space.
I will give the calculations for your specific situation.
Proof. Let $v$ be an arbitrary vector in $\mathbb{E}^3$ and $p\in M$. If $\{e_1,e_2\}$ is an orthonormal basis of $T_p M$, we can extend $e_1, e_2$ to an orthonormal frame $E_1,E_2$ such that $$ \nabla_{E_i} E_j = 0 \quad \text{at $p$ for $i,j=1,2$,} $$ where $\nabla$ is the Levi-Civita connection of $M$. Then at $p$ we have $$ \begin{align*} (\Delta \langle \mathbf{x},v\rangle)_p &= \sum_{i=1}^2 e_i\langle E_i,v\rangle = \sum_{i=1}^2 \langle D_{e_i}E_i,v\rangle \\ &= \sum_{i=1}^2 \langle h(e_i,e_i),v\rangle = 2 \langle H,v\rangle(p). \end{align*} $$ The $D$ stands for the Levi-Civita on $\mathbb{E}^3$. Since both $\Delta x$ and $H$ are independent of the choice of local basis, we have $\langle\Delta x, v\rangle = n \langle H,v \rangle$ for any $v$. Since $v$ was arbitrary and the inner product is non-degenerate, the formula of Beltrami follows.
I also want to mention that your question is related to this recently asked question.