What is the equation of motion for the the following two actions? The Lagrangian densities that I need equations of motion for are the following:
$$L=\sqrt{-g} R^{\mu\nu} R_{\mu\nu}$$
and
$$L=\sqrt{-g} R^{\mu\nu\rho\sigma} R_{\mu\nu\rho\sigma}$$
That is, I want to know how to find the equations of motion varying with respect to the metric for Lagrangians consisting of the volume element and the square of the ricci tensor (the first Lagrangian) and then the volume element and the square of the Riemann tensor (the second Lagrangian). 
 A: You'll need to use the following formulae
\begin{equation}
\begin{split}
\delta \Gamma^c_{ab} &=  \frac{1}{2}  \left( \nabla_a h_b{}^c + \nabla_b h_a{}^c - \nabla^c h_{ab}  \right)   + O(h^2) ~, \\
\delta R^a{}_{bcd} &= \frac{1}{2}  \nabla_c \nabla_d h_b{}^a +  \frac{1}{2}  \nabla_c  \nabla_b h_d{}^a -  \frac{1}{2}  \nabla_c  \nabla^a h_{db}   -\frac{1}{2}  \nabla_d \nabla_c h_b{}^a   -\frac{1}{2}  \nabla_d  \nabla_b h_c{}^a + \frac{1}{2}  \nabla_d  \nabla^a h_{cb}    ~, \\
\delta R_{ab} &= \frac{1}{2} \left( \nabla_c \nabla_a h_b{}^c + \nabla_c \nabla_b h_a{}^c - \nabla^2 h_{ab} - \nabla_a \nabla_b h \right)  + O(h^2) ~, \\
\delta R &= - R_{ab} h^{ab} + \nabla_a \nabla_b h^{ab} - \nabla^2 h  + O(h^2) ~, \\
\delta \det g &=  h \det g  + O(h^2) ~. \\
\end{split}
\end{equation}
where $h_{ab} = \delta g_{ab}$. You should now be able to vary the two actions and determine the equations of motion. 
EDIT - Let me work it out for you for the first action, which is
$$
S = \int d^d x \sqrt{-g} g^{ac} g^{bd} R_{ab} R_{cd} 
$$
Then,
\begin{align}
\delta S &= \int d^d x \left[ \delta \sqrt{-g} g^{ac} g^{bd} R_{ab} R_{cd}  + 2 \sqrt{-g} \delta  g^{ac} g^{bd} R_{ab} R_{cd}  + 2 \sqrt{-g} g^{ac} g^{bd}  \delta R_{ab} R_{cd}  \right] \\
&= \int d^d x   \sqrt{-g} \left[ \frac{1}{2}  h R_{ab} R^{ab} - 2 h^{ab}    R_{ac} R^c{}_b   +  R^{ab} \left( 2 \nabla_c \nabla_a h_b{}^c  - \nabla^2 h_{ab} - \nabla_a \nabla_b h \right) \right] \\
&= \int d^d x   \sqrt{-g} \left[ \frac{1}{2} g_{ab} R_{cd} R^{cd} - 2     R_{ac} R^c{}_b   +  \nabla^c \nabla_a  R_{bc}  +  \nabla^c \nabla_b  R_{ac}  -  \nabla^2 R_{ab}   - g_{ab} \nabla_c  \nabla_d R^{cd}   \right] h^{ab}
\end{align}
We can then read off the equations of motion as
$$
\frac{1}{2} g_{ab} R_{cd} R^{cd} - 2     R_{ac} R^c{}_b   +  \nabla^c \nabla_a  R_{bc}  +  \nabla^c \nabla_b  R_{ac}  -  \nabla^2 R_{ab}   - g_{ab} \nabla_c  \nabla_d R^{cd}   = 0~. 
$$
A: These two forms are equivalent. This is for the Weyl tensor the same as
$$
{\cal L}~=~\sqrt{-g} C^{\mu\nu\rho\sigma} C_{\mu\nu\rho\sigma}
$$
where $C^{\mu\nu\rho\sigma}$ is the Weyl tensor. The Euler-Lagrange equation gives the following differential equation called the Bach equation
$$
\nabla_\mu\nabla_\nu C^{\mu\alpha\beta\nu}~+~\frac{1}{2}R_{\alpha\beta}C^{\mu\alpha\beta\nu}~=~0.
$$
This is pretty interesting, for in conformally flat spacetimes the Ricci curvature is a sort of eigenvalued tenor. Now to generalize to the Riemann curvature we have to use
$$
R^{\mu\alpha\beta\nu}~=~C^{\mu\alpha\beta\nu}~-~\frac{1}{n-2}\left(R_{\mu\nu}g_{\alpha\beta}~-~R_{\mu\beta}g_{\alpha\nu}~+~R_{\alpha\beta}g_{\mu\nu}~-~R_{\alpha\nu}g_{\mu\beta}\right)$$$$~-~\frac{1}{(n-1)(n-2)}R(g_{\mu\beta}g_{\alpha\nu}~-~R_{\mu\nu}g_{\alpha\beta}),
$$
where $n$ is the dimension of the space taken as $n=4$. You can see what this Bach equation looks like.
This Lagrangian turns up in string theory. where the gravitons are described as $\alpha'R^{abcd}R_{abcd}$ for $\alpha'$ the string parameter or tension.
