Metric variation of Laplace-Beltrami I am trying to vary the laplace-Beltrami operator with respect to the metric. Using the following two rules
\begin{align}
\frac{\delta g^{\alpha \beta}}{\delta g^{\mu \nu}} &=\frac{1}{2} \left[\delta^\alpha_\mu \delta^\beta_\nu + \delta^\alpha_\nu \delta^\beta_\mu \right]\\
\frac{\delta g_{\alpha \beta}}{\delta g^{\mu \nu}} &=- \frac{1}{2} \left[g_{\mu \alpha} g_{\nu \beta} + g_{\mu \beta} g_{\nu \alpha} \right],
\end{align}
I have come across the following conundrum when taking the metric variation of the Laplace-Beltrami operator $\square_x = g^{\mu \nu} \partial^x_\mu \partial^x_\nu$
\begin{align}
\frac{\delta \square}{\delta g^{\mu \nu}} &= \frac{\delta}{\delta g^{\mu \nu}} \left[g_{\alpha \beta} \partial^\alpha \partial^\beta \right]=-\frac{1}{2}  \left[g_{\mu \alpha} g_{\nu \beta} + g_{\mu \beta} g_{\nu \alpha} \right] \partial^\alpha \partial^\beta \\
&=-\partial_\mu \partial_\nu\\
&=\frac{\delta}{\delta g^{\mu \nu}} \left[g^{\alpha \beta} \partial_\alpha \partial_\beta \right] = + \frac{1}{2}  \left[\delta^\alpha_\mu \delta^\beta_\nu + \delta^\alpha_\nu \delta^\beta_\mu \right] \partial_\alpha \partial_\beta\\
&=+ \partial_\mu \partial_\nu.
\end{align}
Where did I make a mistake ? I assume it has something to do with the gradient being a covariant object, but I am yet to spot the missing minus sign - where does it come into play ? Thanks in advance for your help.
Edit: How I derived the formulas for the metric variation: Starting from $g_{\mu \alpha} g^{\alpha \nu} = \delta^\nu_\mu$, I find that
\begin{align} 
\delta g_{\mu \alpha} g^{\alpha \nu} + g_{\mu \alpha} \delta g^{\alpha \nu} &=0 \\
\Rightarrow \quad g_{\lambda \nu} g^{\alpha \nu} \delta g_{\mu \alpha} = \delta g_{\mu \lambda} &= -g_{\lambda \nu} g_{\mu \alpha} \delta g^{\alpha \nu} \\
\Leftrightarrow \quad \frac{\delta g_{\mu \lambda}}{\delta g^{\alpha \nu}} &= - g_{\lambda\nu} g_{\mu \alpha},
\end{align}
which I subsequently symmetrized in the Lorentz indices. Then I used that $g^{\mu \nu} g_{\mu \nu}=d$ to derive that
\begin{align} 
\delta(g^{\mu \nu} g_{\mu \nu}) &=0 = g^{\mu \nu } \delta g_{\mu \nu} + \delta g^{\mu \nu} g_{\mu \nu} \\
\Rightarrow \quad g^{\mu \nu} \delta g_{\mu \nu} &=- \delta g^{\mu \nu} g_{\mu \nu},
\end{align}
such that the other rule follows analogously with a minus sign and raised indices. Is this correct ?
 A: The “raised-index derivative” $\partial^\mu$ depends implicitly on the inverse metric, since it stands for
$$
g^{\mu \nu} \frac{\partial}{\partial x^\nu}. $$
The reason for this is that mathematically, the coordinate derivative with a lowered index can be defined without reference to a metric, and so it is natural to take it to be constant when doing a variation like this.  By contrast, the definition of the “raised-index” coordinate derivative relies on the metric and so we can't “hold it fixed” while varying with respect to the metric.
Finally, note that there's a subtlety in your derivation that you haven't taken into account.  If the background metric isn't flat, then the Laplace-Beltrami operator isn't $g^{\mu \nu} \partial_\mu \partial_\nu$ but $g^{\mu \nu} \nabla_\mu \nabla_\nu$.  If this operates on some scalar field $f$ then we have
$$
g^{\mu \nu} \nabla_\mu \nabla_\nu f = g^{\mu \nu} \left( \partial_\mu \partial_\nu f - \Gamma^{\rho} {}_{\mu \nu} \partial_\rho f \right)
$$
for some scalar field $f$.  Since the Christoffels $\Gamma^{\rho} {}_{\mu \nu}$ depend on the metric, they must be varied as well.  The resulting terms will involve coordinate derivatives of $g^{\mu \nu}$ which must then be integrated by parts to find $\delta (\Box f)/\delta g^{\mu \nu}$.  If the Laplace-Beltrami operator is applied to vector or tensor fields, additional terms involving the Christoffels will be present that must be varied as well.
A: It is not true that $g_{ab} \partial^a \partial^b = g^{ab} \partial_a \partial_b$, because the left derivative in the left hand side acts on a metric too. The correct definition is $\square = g^{ab}\partial_a \partial_b$, from which we get
$$\begin{align}
g_{ab} \partial^a \partial^b &= g_{ab} g^{ac} \partial_c (g^{bd} \partial_d) \\
&= g_{ab} g^{ac} g^{bd} \partial_c \partial_d + g_{ab} g^{ac} (\partial_c g^{bd}) \partial_d \\
&= g^{cd} \partial_c \partial_d + (\partial_b g^{bd}) \partial_d
\end{align}$$
