Solution (source code)

= Solution

Varying the coordinate expression for the curvature and replacing partial derivatives by covariant ones gives
$$
\delta R_{ab}=\nabla_c\delta\Gamma^c{}_{ab}
-\nabla_b\delta\Gamma^c{}_{ac}.
$$
Therefore
$$
\delta R=R_{ab}\delta g^{ab}+g^{ab}\delta R_{ab}
=-R^{ab}h_{ab}+\nabla_aX^a,
$$
where
$$
\boxed{X^a=\nabla_bh^{ab}-\nabla^ah}.
$$
This is the <metric variation of scalar curvature>; its divergence contributes only a boundary term to the variation of the <Einstein-Hilbert action>.