Variation of modified Einstein Hilbert Action In general relativity one can derive the Einstein Field Equations by the principle of least action through variations with respect to the inverse of the metric tensor. In some modified theories of gravity, such as the Brans-Dicke Theory, a scalar field is added to the Einstein Hilbert Action and the gravitational constant is replaced by a function of the scalar field. I am not quite sure how to derive the field equations from this action, more specifically the part where the scalar field is attached to the Ricci scalar $\phi R$.
The Brans-Dicke Action is $$S_{BD} = \int d^4x \sqrt{-g}\left[ \frac{1}{16\pi}\left(\phi R - \frac{\omega }{\phi}g^{ab} \partial _a\phi\partial _b \phi \right) +L_M \right].$$
The resulting field equation is $$G_{ab} = \frac{8\pi}{\phi}T_{ab}+\frac{\omega}{\phi^2} (\partial_a\phi\partial_b\phi-\frac{1}{2}g_{ab}\partial_c\phi\partial^c\phi) +\frac{1}{\phi}(\nabla_a\nabla_b\phi-g_{ab}\Box\phi).$$
I also want to derive a new field equation for practice. So my questions are:


*

*How does one derive the equations of motion?

*How to perform the variation of the following action? 
$$S=\int d^4x \sqrt{-g} \left[ \frac{1}{16\pi G} R - \phi( \nabla_{\mu} g_{ab} \nabla_{\nu}g_{ab}) - 2\Lambda+L_M) \right] $$
The Ricci scalar, the cosmological constant, and the matter Lagrangian will variate simply like the Einstein Hilbert Action to: 
$$\delta S = \int d^4x  \sqrt{-g} \left[ \frac{1}{\kappa} \left( R_{ab}-\frac{1}{2}Rg_{ab}+\Lambda g_{ab} \right) -T_{ab} \right]\delta g^{ab}.$$ What about the extra term?
Would one simply variate with respect to $\phi$, or is the variation of the covariant derivative of the metric tensor also required? If the latter is true, then would the variation of this extra term be $$\frac{\partial L}{\partial g_{ab}}-\partial _\mu\frac{\partial L}{\partial (\nabla_{\mu}g_{ab})}=0.$$
 Any help would be appreciated. 
By the way, is $\nabla_{\mu}g_{ab}\nabla_{\nu}g_{ab}$ an expression that shows the rate of change (derivative) of the metric tensor with respect to a coordinate $(t, x, y, z)$? 
 A: Find the answer to question 1 below. Question 2. is weird since $\nabla_\mu g_{\alpha\beta} = 0$ (when the connection is metric compatible) as mentioned by @Trimok. In any case, the variation of the action can be derived using the method described below. 
We start with the BD action
$$
S = \frac{1}{16\pi}\int d^4 x \sqrt{-g} \left[ \phi R  - \frac{\omega}{\phi} g^{\mu\nu} \partial_\mu \phi \partial_\nu \phi  \right] + S_M
$$
where $S_M$ is the matter action. To determine Einstein's field equations, we vary the action w.r.t to the metric. We will use the formulae (ref. wikipedia)
\begin{equation}
\begin{split}
\delta R = R_{\mu\nu} \delta g^{\mu\nu} + \nabla_\sigma \left( g^{\mu\nu} \delta \Gamma^\sigma_{\mu\nu} - g^{\mu\sigma} \delta \Gamma^\rho_{\rho\mu} \right)
\end{split}
\end{equation}
The variation of the Christoffel tensor is
\begin{equation}
\begin{split}
\delta \Gamma^\lambda_{\mu\nu} & = \delta g^{\lambda\rho} g_{\rho\alpha} \Gamma^\alpha_{\mu\nu} + \frac{1}{2} g^{\lambda\rho} \left( \partial_\mu \delta g_{\nu\rho} + \partial_\nu \delta g_{\mu\rho} - \partial_\rho \delta  g_{\mu\nu} \right) \\
& = \frac{1}{2} g^{\lambda\rho} \left( \nabla_\mu \delta g_{\nu\rho}  +  \nabla_\nu \delta g_{\mu\rho}  -  \nabla_\rho \delta g_{\mu\nu}   \right) \\
& = - \frac{1}{2}  \left( g_{\nu\alpha} \nabla_\mu \delta g^{\alpha\lambda}  + g_{\mu\alpha}  \nabla_\nu \delta g^{\alpha\lambda}  -  g_{\mu\alpha} g_{\nu\beta} \nabla^\lambda \delta g^{\alpha\beta}   \right)
\end{split}
\end{equation}
where we used $\delta g_{\mu\nu} = - g_{\mu\alpha} g_{\nu\beta} \delta g^{\alpha\beta}$. This implies
\begin{equation}
\begin{split}
g^{\mu\nu} \delta \Gamma^\sigma_{\mu\nu}  &= -\nabla_\alpha \delta g^{\alpha\sigma} + \frac{1}{2} g_{\alpha\beta}\nabla^\sigma\delta g^{\alpha\beta}  \\
g^{\mu\sigma} \delta \Gamma^\lambda_{\lambda\mu} &=  - \frac{1}{2}g_{\alpha\beta} \nabla^\sigma \delta g^{\alpha\beta}  
\end{split}
\end{equation}
which implies
\begin{equation}
\begin{split}
\delta R = R_{\mu\nu} \delta g^{\mu\nu} - \nabla_\mu \nabla_\nu \delta g^{\mu\nu} + g_{\mu\nu} \nabla^2 \delta g^{\mu\nu} 
\end{split}
\end{equation}
Finally, from 1, we also have
$$
\delta \sqrt{-g} = - \frac{1}{2} \sqrt{-g} g_{\mu\nu} \delta g^{\mu\nu} 
$$
Finally, we are ready to compute the variation of the action. We have
\begin{equation}
\begin{split}
\delta S &=  \frac{1}{16\pi}\int d^4 x \delta \sqrt{-g} \left[ \phi R  - \frac{\omega}{\phi} g^{\mu\nu} \partial_\mu \phi \partial_\nu \phi  \right] \\
&~~~~~~~~~~~~ + \frac{1}{16\pi}\int d^4 x  \sqrt{-g} \left[ \phi \delta R  - \frac{\omega}{\phi} \delta g^{\mu\nu} \partial_\mu \phi \partial_\nu \phi  \right]  + \delta S_M \\
&= - \frac{1}{32\pi}\int d^4 x  \sqrt{-g} g_{\mu\nu} \left[ \phi R  - \frac{\omega}{\phi} g^{\alpha\beta} \partial_\alpha \phi \partial_\beta \phi  \right] \delta g^{\mu\nu} + \int d^4 x \frac{ \delta S_M}{\delta g^{\mu\nu}} \delta g^{\mu\nu}   \\
&~~~~~~~~~~~~ + \frac{1}{16\pi}\int d^4 x  \sqrt{-g} \left[ \left( \phi  R_{\mu\nu}  - \nabla_\mu \nabla_\nu \phi   + g_{\mu\nu} \nabla^2  \phi   \right)  - \frac{\omega}{\phi} \partial_\mu \phi \partial_\nu \phi  \right] \delta g^{\mu\nu}  
\end{split}
\end{equation}
Requiring that the variation of the action vanish (to leading order in $\delta g^{\mu\nu}$) gives
\begin{equation}
\begin{split}
G_{\mu\nu}   &=   - \frac{16\pi}{\phi \sqrt{-g}} \frac{ \delta S_M}{\delta g^{\mu\nu}}     +  \frac{\omega}{\phi^2} \left[  \partial_\mu \phi \partial_\nu \phi   -  \frac{1}{2} g_{\mu\nu}  \partial_\alpha \phi \partial^\alpha \phi  \right]   +\frac{1}{\phi} \left[  \nabla_\mu \nabla_\nu \phi   - g_{\mu\nu} \nabla^2  \phi    \right]  
\end{split}
\end{equation}
Recall that the stress-tensor is defined as
\begin{equation}
\begin{split}
T_{\mu\nu} =  - \frac{2}{\sqrt{-g}} \frac{ \delta S_M}{\delta g^{\mu\nu}} 
\end{split}
\end{equation}
Thus
\begin{equation}
\begin{split}
G_{\mu\nu}   &=   \frac{8\pi}{\phi}  T_{\mu\nu}   +  \frac{\omega}{\phi^2} \left[  \partial_\mu \phi \partial_\nu \phi   -  \frac{1}{2} g_{\mu\nu}  \partial_\alpha \phi \partial^\alpha \phi  \right]   +\frac{1}{\phi} \left[  \nabla_\mu \nabla_\nu \phi   - g_{\mu\nu} \nabla^2  \phi    \right]  
\end{split}
\end{equation}
which is the Brans-Dicke equation.
