Einstein field equations from the Hilbert action and Bianchi identity

Step 1 (decompose the variation of the gravitational integrand). Write $R=g^{\mu\nu}R_{\mu\nu}$. By the product rule, $$\delta\big((R-2\Lambda)\sqrt{-g}\big) = R_{\mu\nu}\,\sqrt{-g}\,\delta g^{\mu\nu} + g^{\mu\nu}(\delta R_{\mu\nu})\sqrt{-g} + (R-2\Lambda)\,\delta\sqrt{-g}.$$

**Step 2 (Jacobi's formula for $\delta\sqrt{-g}$ and the Palatini identity).** Using Jacobi's identity $\delta(\det g)=(\det g)\,g^{\mu\nu}\delta g_{\mu\nu}=-(\det g)\,g_{\mu\nu}\delta g^{\mu\nu}$, we get $\delta\sqrt{-g}=-\tfrac12\sqrt{-g}\,g_{\mu\nu}\delta g^{\mu\nu}$. Meanwhile, the Palatini identity expresses the Ricci variation as a covariant divergence, $g^{\mu\nu}\delta R_{\mu\nu}=\nabla_\alpha\big(g^{\mu\nu}\delta\Gamma^\alpha_{\mu\nu}-g^{\alpha\nu}\delta\Gamma^\mu_{\mu\nu}\big)$, which integrates to zero on $M$ by Stokes's theorem for variations $\delta g^{\mu\nu}$ supported away from the boundary.

**Step 3 (assemble $\delta S=0$).** Combining Steps 1 and 2 and the definition $T_{\mu\nu}\equiv-\dfrac{2}{\sqrt{-g}}\dfrac{\delta S_{\mathrm{matter}}}{\delta g^{\mu\nu}}$, $$\delta S = \int_M \left[\frac{c^4}{16\pi G}\left(R_{\mu\nu}-\tfrac12 Rg_{\mu\nu}+\Lambda g_{\mu\nu}\right) - \frac12 T_{\mu\nu}\right]\delta g^{\mu\nu}\,\sqrt{-g}\,d^4x = 0.$$ Since $\delta g^{\mu\nu}$ is an arbitrary symmetric tensor variation, the bracket must vanish pointwise, giving $R_{\mu\nu} - \tfrac{1}{2}R g_{\mu\nu} + \Lambda g_{\mu\nu} = \dfrac{8\pi G}{c^4} T_{\mu\nu}$.

Step 4 (contracted Bianchi identity and conservation). The second Bianchi identity $\nabla_\lambda R_{\rho\sigma\mu\nu}+\nabla_\rho R_{\sigma\lambda\mu\nu}+\nabla_\sigma R_{\lambda\rho\mu\nu}=0$, contracted twice with $g^{\lambda\mu}g^{\rho\nu}$, yields $\nabla^\mu R_{\mu\nu}-\tfrac12\nabla_\nu R=0$, i.e. $\nabla^\mu G_{\mu\nu} = 0$ for $G_{\mu\nu}=R_{\mu\nu}-\tfrac12 Rg_{\mu\nu}$. Because $\nabla^\mu g_{\mu\nu}=0$ (metric compatibility), taking $\nabla^\mu$ of $R_{\mu\nu} - \tfrac{1}{2}R g_{\mu\nu} + \Lambda g_{\mu\nu} = \dfrac{8\pi G}{c^4} T_{\mu\nu}$ forces $\nabla^\mu T_{\mu\nu}=0$.