How to vary a second order function with respect to the metric tensor? Can anybody help me to prove this relation, how is it is valid ?
\begin{equation}
\frac{\delta}{\delta g^{\mu\nu}}\nabla_{\sigma}\Bigr(\alpha(x^{\beta})\,\frac{\nabla^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Bigl)=\frac{1}{2}g_{\mu\nu}\nabla_{\sigma}\Bigr(\alpha(x^{\beta})\,\frac{\nabla^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Bigl)
\end{equation}
Here $x^{\beta}$ is just space time coordinate.  $\alpha(x^{\beta})$ and  ${\phi(x^\beta)}$ are scalar functions and Einstein summation  convention is used. The term that I take variation with respect to metric tensor $g_{\mu\nu}$ is just a scalar function.
In addition, how can I find 
\begin{equation}
\frac{\delta}{\delta g^{\mu\nu}}(\Box{\phi})=?
\end{equation}
Also
\begin{equation}
\frac{\delta}{\delta g^{\mu\nu}}(\Box^n{\phi})=?
\end{equation}
Regarding the first question, I tried the following way,
\begin{equation}
\frac{\delta}{\delta g^{\mu\nu}}\nabla_{\sigma}\Bigr(\alpha(x^{\beta})\,\frac{\nabla^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Bigl)=\frac{\delta}{\delta g^{\mu\nu}}\frac{1} {\sqrt{-g}}\partial_{\sigma}\Biggr({\sqrt{-g}}\alpha(x^{\beta})\,\frac{\partial^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Biggl)\\ 
\end{equation}
\begin{equation}
\frac{\delta}{\delta g^{\mu\nu}}\nabla_{\sigma}\Bigr(\alpha(x^{\beta})\,\frac{\nabla^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Bigl)=\frac{\delta}{\delta g^{\mu\nu}}\Bigr(\frac{1} {\sqrt{-g}}\Bigl)\partial_{\sigma}\Biggr({\sqrt{-g}}\alpha(x^{\beta})\,\frac{\partial^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Biggl)+\Bigr(\frac{1} {\sqrt{-g}}\Bigl)\frac{\delta}{\delta g^{\mu\nu}}\partial_{\sigma}\Biggr({\sqrt{-g}}\alpha(x^{\beta})\,\frac{\partial^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Biggl)\\ 
\end{equation}
Using the result
\begin{equation}
\frac{\delta}{\delta g^{\mu\nu}}\Bigr(\frac{1} {\sqrt{-g}}\Bigl)=(\frac{-1}{-g})\frac{-1}{2}{\sqrt{-g}}\, g_{\mu\nu}=\frac{1}{2}g_{\mu\nu}\frac{1}{\sqrt{-g}}
\end{equation}
I get
\begin{equation}
\frac{\delta}{\delta g^{\mu\nu}}\nabla_{\sigma}\Bigr(\alpha(x^{\beta})\,\frac{\nabla^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Bigl)=\frac{1}{2} g_{\mu\nu}\frac{1} {\sqrt{-g}}\partial_{\sigma}\Biggr({\sqrt{-g}}\alpha(x^{\beta})\,\frac{\partial^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Biggl)+\Bigr(\frac{1} {\sqrt{-g}}\Bigl)\frac{\delta}{\delta g^{\mu\nu}}\partial_{\sigma}\Biggr({\sqrt{-g}}\alpha(x^{\beta})\,\frac{\partial^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Biggl)\\  
\end{equation}
\begin{equation}
\frac{\delta}{\delta g^{\mu\nu}}\nabla_{\sigma}\Bigr(\alpha(x^{\beta})\,\frac{\nabla^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Bigl)=\frac{1}{2}g_{\mu\nu}\nabla_{\sigma}\Bigr(\alpha(x^{\beta})\,\frac{\nabla^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Bigl)+\Bigr(\frac{1} {\sqrt{-g}}\Bigl)\frac{\delta}{\delta g^{\mu\nu}}\partial_{\sigma}\Biggr({\sqrt{-g}}\alpha(x^{\beta})\,\frac{\partial^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Biggl)\\  
\end{equation}
My trouble is, I was not able to prove 
\begin{equation}
\Bigr(\frac{1} {\sqrt{-g}}\Bigl)\frac{\delta}{\delta g^{\mu\nu}}\partial_{\sigma}\Biggr({\sqrt{-g}}\alpha(x^{\beta})\,\frac{\partial^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Biggl)=0\\
\end{equation}
The difficulties is, I don't know how to evaluate the variation of the derivative of the determinant of the metric tensor.
 A: Thank you for reading my post and for the comments.
I found a way to prove the last term vanishes. This was trivial.
I need to prove,
\begin{equation}
\Bigr(\frac{1} {\sqrt{-g}}\Bigl)\frac{\delta}{\delta g^{\mu\nu}}\partial_{\sigma}\Biggr({\sqrt{-g}}\alpha(x^{\beta})\,\frac{\partial^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Biggl)=0\\
\end{equation}
This is the part of the action
\begin{equation}
A=\int{d^4x \sqrt{-g}\frac{1} {\sqrt{-g}}\frac{\delta}{\delta g^{\mu\nu}}\partial_{\sigma}\Biggr({\sqrt{-g}}\alpha(x^{\beta})\,\frac{\partial^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Biggl)}=\int{d^4x \frac{\delta}{\delta g^{\mu\nu}}\partial_{\sigma}\Biggr({\sqrt{-g}}\alpha(x^{\beta})\,\frac{\partial^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Biggl)}
\end{equation}
\begin{equation}
A=\int{d^4x \frac{\delta}{\delta g^{\mu\nu}}\partial_{\sigma}\Biggr({\sqrt{-g}}\alpha(x^{\beta})\,\frac{\partial^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Biggl)}=\frac{\delta}{\delta g^{\mu\nu}}\int{d^4x \partial_{\sigma}\Biggr({\sqrt{-g}}\alpha(x^{\beta})\,\frac{\partial^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Biggl)}
\end{equation}
\begin{equation}
A=\frac{\delta}{\delta g^{\mu\nu}}\int{d^4x \partial_{\sigma}\Biggr({\sqrt{-g}}\alpha(x^{\beta})\,\frac{\partial^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Biggl)}=\frac{\delta}{\delta g^{\mu\nu}}\Biggr({\sqrt{-g}}\alpha(x^{\beta})\,\frac{\partial^{\sigma}{\phi(x^\beta)}}{\phi(x^\beta)}\Biggl)|_{on\,the\,boundary}=0
\end{equation}
This is true if the variation of the boundary is zero, then this surface term vanishes.
