Solution (source code)

= Solution

For a scalar, $\nabla_a\phi=\partial_a\phi$, while its second <covariant derivative> is $\nabla_a\nabla_b\phi=\partial_a\partial_b\phi-\Gamma^c_{ab}\partial_c\phi$. <Partial derivatives> commute and the <Levi-Civita connection> is torsion-free, $\Gamma^c_{ab}=\Gamma^c_{ba}$, proving the asserted symmetry.

Write the <Einstein equation> without <cosmological constant> as $G_{ab}=\kappa T_{ab}$, $\kappa=8\pi G/c^4$. Tracing the given Ricci equation gives $R=g^{ab}\partial_a\phi\partial_b\phi$. Hence
$$
\boxed{T_{ab}=\kappa^{-1}\left(\partial_a\phi\partial_b\phi-\frac12g_{ab}\partial_c\phi\partial^c\phi\right).}
$$
In units $\kappa=1$ the prefactor is absent. <Metric compatibility>, the product rule and the just-proved <symmetric Hessian of a scalar field> give
$$
\kappa\nabla^aT_{ab}=(\nabla^a\nabla_a\phi)\partial_b\phi+\partial^a\phi\nabla_a\partial_b\phi-\partial^a\phi\nabla_b\partial_a\phi=(\nabla^a\nabla_a\phi)\partial_b\phi.
$$
The contracted <Bianchi identity> makes the left side zero. Wherever the gradient covector is nonzero, at least one component is nonzero, forcing \b[$\nabla_a\nabla^a\phi=0$], even when that gradient is null.

Finally $\Gamma^a_{ab}=\partial_b\log\sqrt{-g}$, obtained by differentiating the metric <determinant>. Thus for any vector $V$, $\nabla_aV^a=(\sqrt{-g})^{-1}\partial_a(\sqrt{-g}V^a)$. Applying it to $V^a=g^{ab}\partial_b\phi$ proves
$$
\boxed{\partial_a(\sqrt{-g}\,g^{ab}\partial_b\phi)=0.}
$$