Hay una muy agradable derivación de esta identidad el uso de formas diferenciales, y se evita completamente todo el desorden de los símbolos de Christoffel.
La cosa agradable sobre formas diferenciales es que el exterior derivada puede calcularse mediante cualquier operador de la derivada, lo que nos permite comparar las expresiones obtenemos con la derivada covariante a la expresión que se obtendría utilizando derivadas parciales.
El tratamiento de la $\phi$ como una 0-forma, vamos a calcular $*d*d\phi$. Hacerlo con covariante derivados de reproducir las $\nabla_a \nabla^a \phi$, mientras que en derivadas parciales que vamos a conseguir la fórmula que involucra $\sqrt{-g}$:
\begin{align}
*d*d\phi &= \frac{1}{4!}\epsilon_{abcd}\cdot4\nabla^{[a}\epsilon^{|e|bcd]}\nabla_e\phi \\
&=\frac{4}{4!}\epsilon_{abcd}\epsilon^{ebcd}\nabla^a\nabla_e\phi\\
&=\frac{4\cdot 3!}{4!}(-\delta^e_a)\nabla^a\nabla_e\phi \\
&=-\nabla_a\nabla^a\phi
\end{align}
Hemos utilizado que la derivada covariante de la $\epsilon$ tensor es de cero, y también el hecho de que cuando usted contrata a los índices de dos $\epsilon$ tensores de obtener una generalizada delta de Kronecker. El signo menos viene del hecho de que $\epsilon$ está normalizado con un factor de $\sqrt{-g}$, mientras que cuando se levanta todos los índices de la épsilon del tensor con una relación inversa entre la métrica de obtener un total factor de $g^{-1}$, lo $\dfrac{\sqrt{-g}}{g} = -\dfrac{1}{\sqrt{-g}}$.
Ahora hacemos lo mismo con derivadas parciales, y el uso de la normalización de la $\epsilon$ tensor, $\epsilon_{abcd} = \sqrt{-g}\tilde{\epsilon}_{abcd}$ donde $\tilde{\epsilon}_{abcd}$ es la alternancia de símbolo, con valores de $+1$, $-1$ o $0$ (y, por tanto, ha de fuga derivadas parciales).
\begin{align}
*d*d\phi &=\frac{1}{4!}\epsilon^{abcd}\cdot4\partial_{[a}(\epsilon_{|e|bcd]}\partial^e\phi)\\
&=\frac{4}{4!}\frac{-1}{\sqrt{-g}}\tilde{\epsilon}^{abcd}\partial_a(\sqrt{-g}\tilde{\epsilon}_{ebcd}g^{ef}\partial_f\phi)\\
&=-\frac{4\cdot3!}{4!}\delta^e_a\frac{1}{\sqrt{-g}}\partial_a(\sqrt{-g})g^{ef}\partial_f\phi)\\
&=-\frac{1}{\sqrt{-g}}\partial_a(\sqrt{-g}g^{af}\partial_f\phi)
\end{align}
Por lo tanto, vemos por qué las dos expresiones son iguales: está íntimamente relacionado con las formas diferenciales descripción de este operador.
Para formas diferenciales de rango $1$ o mayor, el operador $\nabla_a\nabla^a$ está relacionado con el de Laplace-Beltrami operador $\triangle = *d*d+d*d*$, pero no es exactamente igual a este; la diferencia de los dos operadores es proporcional a la curvatura del tensor. Sin embargo, el análogo de la derivación anterior puede ser derivado a un diferencial de $p$forma $\alpha$, con el operador $*d*$, es decir,
$$\nabla_{a_1} \alpha^{a_1\ldots a_p}= -*d*\alpha = \frac{1}{\sqrt{-g}}\partial_{a_1}\left(\sqrt{-g} \alpha^{a_1\ldots a_p}\right)$$
Y, como otros han mencionado, la derivación no se basa en las ecuaciones de movimiento, y por lo tanto tiene off-shell.