Solution (source code)

= Solution

Define the <stress-energy tensor> by <metric variation> with the scalar held fixed:
$$
\delta S=\frac12\int d^4x\,\sqrt{-g}\,T^{ab}\delta g_{ab}.
$$
Use compactly supported variations to remove boundary terms. The first <covariant derivative> of a <scalar field> is an ordinary derivative, so it has no connection variation. Varying the kinetic term and the <volume form> gives
$$
T^{\mathrm{kin}}_{ab}=\nabla_a\Phi\nabla_b\Phi-\frac12g_{ab}\nabla_c\Phi\nabla^c\Phi.
$$

For the curvature term set $f=\Phi^2$. Using part (a), its variation is
$$
\delta S_\xi=-\xi\int\sqrt{-g}\left[
\left(\frac12Rf\,g^{ab}-fR^{ab}\right)h_{ab}
-f\Box h+f\nabla^a\nabla^b h_{ab}\right]d^4x.
$$
Applying <integration by parts> twice moves each pair of derivatives from the metric variation onto $f$. The result is
$$
\delta S_\xi=\xi\int\sqrt{-g}
\left[fG^{ab}+g^{ab}\Box f-\nabla^a\nabla^b f\right]h_{ab}\,d^4x,
$$
where $G_{ab}=R_{ab}-Rg_{ab}/2$ is the <Einstein tensor> and $\Box=\nabla^c\nabla_c$ is the <d'Alembert operator>. Comparing with the defining variation gives the <stress-energy tensor of a nonminimally coupled scalar field>:
$$
\boxed{T_{ab}=\nabla_a\Phi\nabla_b\Phi
-\frac12g_{ab}(\nabla\Phi)^2
+2\xi\left[\Phi^2G_{ab}+g_{ab}\Box(\Phi^2)-\nabla_a\nabla_b(\Phi^2)\right].}
$$
The factor $2\xi$ follows from the action's coupling $-\xi R\Phi^2$; a convention using $-\xi R\Phi^2/2$ would instead give $\xi$ in that bracket. At $\xi=0$ this reduces to the massless <Klein-Gordon scalar stress-energy tensor>.