Solution (source code)

= Solution

Put
$$
Y=\nabla_\rho\phi\nabla^\rho\phi.
$$
Using metric compatibility and the <symmetric Hessian of a scalar field>, the divergence of the proposed tensor is
$$
\begin{aligned}
2\nabla^\mu T_{\mu\nu}
&=(\nabla^\mu\nabla_\mu\phi)\nabla_\nu\phi
+\nabla^\mu\phi\,\nabla_\mu\nabla_\nu\phi
+A\nabla_\nu Y+2B\phi\nabla_\nu\phi\\
&=m^2\phi\nabla_\nu\phi
+\frac12\nabla_\nu Y
+A\nabla_\nu Y+2B\phi\nabla_\nu\phi\\
&=\left(A+\frac12\right)\nabla_\nu Y
+(m^2+2B)\phi\nabla_\nu\phi.
\end{aligned}
$$
For this to vanish for arbitrary solutions and local field data, the two coefficients must vanish. Hence
$$
\boxed{A=-\frac12,\qquad B=-\frac{m^2}{2}}.
$$
The resulting tensor is the <Conserved Klein-Gordon scalar stress-energy ansatz>.