Is there an invariant way to derive the energy-momentum tensor? On a (pseudo-)Riemannian manifold $(M, g)$ I can define the following action for any $\phi \in C ^{\infty}(M)$:
$$
  \mathcal{S}(\phi) = \int_M g(\text{grad }\phi, \text{grad }\phi)
  \mathrm{d} V.
$$
According to the Wikipedia page, the energy-momentum tensor associated
with this action is given by the variation through a ''spacetime
translation''. I take this to mean that given a vector field $X$ on
$M$ with flow $\psi$, the composition $\Phi(t) :=  \phi \circ \psi_t$
gives a variation of $\phi$ as $t$ varies from some small $-\epsilon$
to $\epsilon$ (at least if $X$ is compactly supported). Then for any
$v \in T_p M$
\begin{align*}
g_p(\text{grad } \phi \circ \psi_t|_p, v) &= \mathrm{d}( \phi \circ
                                       \psi_t)|_p (v) =  v( \phi
                                       \circ \psi_t)
  \\
  &=
                                      \mathrm{d} \psi_t|_p(v) \phi =
                                   \mathrm{d} \phi|_{\psi_t(p)}(\mathrm{d} \psi_t
    |_p(v))
  \\
  &= g_{\psi_t(p)}(\text{grad } \phi |_{\psi_t(p)}, \mathrm{d}
    \psi_t|_p(v))
  \\
  &= g_p((\mathrm{d} \psi_t) ^{\top}|_{\psi_t(p)}\text{grad } \phi |_{\psi_t(p)}, v),
\end{align*}
where $\mathrm{d} \psi_t|_{\psi_t(p)} : T_{\psi_t(p)} M \to T_pM$ is
the adjoint of $\mathrm{d} \psi_t|_p : T_p M \to T_{\psi_t(p)} M$ with
respect to $g$. Thus, by diffeomorphism invariance of integrals:
\begin{align*}
\mathcal{S}(\Phi(t)) &= \int_M g_p(\text{grad }\phi \circ \psi_t|_p,
                       \text{grad }\phi \circ \psi_t|_p) \mathrm{d} V
  \\
  &= \int_M g_p((\mathrm{d} \psi_t) ^{\top}|_{\psi_t(p)}\text{grad }
    \phi |_{\psi_t(p)}, (\mathrm{d} \psi_t)
    ^{\top}|_{\psi_t(p)}\text{grad } \phi |_{\psi_t(p)}) \mathrm{d} V
  \\
&= \int_M g_{\psi_{-t}(p)}((\mathrm{d} \psi_t) ^{\top}\text{grad }
    \phi |_p, (\mathrm{d} \psi_t)
    ^{\top}\text{grad } \phi |_p) \psi_{-t} ^* \mathrm{d} V
\end{align*}
Upon taking the derivative with respect to $t$ and using the Leibniz
rule I see that
\begin{align*}
\frac{\mathrm{d}}{\mathrm{d} t} \Big|_{t = 0} \mathcal{S}(\Phi(t)) &=
                                                                     \int_M
                                                                     \frac{\mathrm{d}}{\mathrm{d}
                                                                     t}
                                                                     \Big|_{t
                                                                     =
                                                                     0}g_{\psi_{-t}(p)}((\mathrm{d} \psi_t) ^{\top}\text{grad }
    \phi |_p, (\mathrm{d} \psi_t)
    ^{\top}\text{grad } \phi |_p) \mathrm{d} V
  \\
  &- \int_M g_p(\text{grad } \phi|_p, \text{grad } \phi|_p)
    \text{div}(X) \mathrm{d} V.
\end{align*}
If I compare this with the the contraction of the actual energy-momentum
tensor against the deformation tensor $\pi ^{\alpha \beta} = \nabla
^{\alpha} X ^{\beta} + \nabla ^{\beta} X ^{\alpha}$ I get
\begin{align*}
T_{\alpha \beta} \pi ^{\alpha \beta} &= \nabla_{\alpha} \phi
                                       \nabla_{\beta} \phi 
                                       (\nabla ^{\alpha} X ^{\beta} +
                                       \nabla ^{\beta} X ^{\alpha}) -
                                       \frac{1}{2} g_{\alpha \beta}
                                       \nabla ^{\gamma} \phi
                                       \nabla_{\gamma} \phi (\nabla ^{\alpha} X ^{\beta} +
                                       \nabla ^{\beta} X ^{\alpha})
  \\
  &= 2g(\nabla_{\text{grad } \phi} X, \text{grad }\phi) - g(\text{grad
    }\phi, \text{grad }\phi) \text{div}(X).
\end{align*}
The divergence terms match, so all I have to do now is show that
\begin{align*}
  \frac{\mathrm{d}}{\mathrm{d} t} \Big|_{t =
    0}g_{\psi_{-t}(p)}((\mathrm{d} \psi_t) ^{\top}\text{grad } \phi
  |_p, (\mathrm{d} \psi_t) ^{\top}\text{grad } \phi |_p) &= 2
                                                           g(\nabla_{\text{grad }\phi} X, \text{grad }\phi)
  \\
  &= \mathcal{L}_X g(
  \text{grad }\phi, \text{grad }\phi).
\end{align*}
But this does not seem right. The left hand side is almost
$-\mathcal{L}_X g(\text{grad }\phi, \text{grad }\phi)$, in fact if $X$
is Killing this is the case, so I am not sure if the variation I am
using is the correct one, or what Wikipedia means by ''the energy
tensor is the Noether current of spacetime translations''. I have seen
derivations by using coordinates, but is there a way to proceed
invariantly similar to what I am trying to do above?
 A: Computing
\begin{align*}
\mathrm{d}(\phi \circ \psi_t)|_p(Y) &= Y|_p \phi \circ \psi_t =
                                      \mathrm{d}\psi_t|_p(Y) \phi =
                                      \mathrm{d} \phi|_{\psi_t(p)}
                                      (\mathrm{d} \psi_t|_p(Y)) =
                                      (\psi_t ^* \mathrm{d} \phi)_p(Y)
\end{align*}
yields
\begin{align*}
g_p ^{-1}(\mathrm{d}(\phi \circ \psi_t), \mathrm{d}(\phi \circ
  \psi_t)) &= g_p ^{-1}(\psi_t ^* \mathrm{d} \phi, \psi_t ^*
                \mathrm{d} \phi).
\end{align*}
Using diffeomorphism invariance
\begin{align*}
\mathcal{S}(\Phi(t)) = \int_M g ^{-1}_p(\mathrm{d} (\phi \circ \psi_t), \mathrm{d} ( \phi
  \circ \psi_t)) \mathrm{d}V &= \int_M g
                               ^{-1}_{\psi_{-t}(p)}((\psi_{-t} ^{-1}) ^* \mathrm{d} \phi, (\psi_{-t} ^{-1}) ^*
                               \mathrm{d} \phi) \phi_{-t} ^*\mathrm{d}V
  \\
  &= \int_M (\psi_{-t} ^* g ^{-1})_p(\mathrm{d} \phi, \mathrm{d} \phi)
    \phi_{-t} ^* \mathrm{d}V
\end{align*}
and taking derivatives we have
\begin{align*}
  \frac{\mathrm{d}}{\mathrm{d} t} \Big|_{t = 0} \mathcal{S}(\Phi(t))
  &= \int_M \mathcal{L}_{-X} g ^{-1}(\mathrm{d} \phi ,\mathrm{d}\phi)
    \mathrm{d} V + \int_M g ^{-1}(\mathrm{d} \phi, \mathrm{d} \phi)
    \mathcal{L}_{-X} \mathrm{d} V
  \\
  &=                                                                   -\int_M
                                                                     (\mathcal{L}_X
                                                                     g
                                                                     ^{-1}(\mathrm{d}
                                                                     \phi,
                                                                     \mathrm{d}\phi)
                                                                     +
                                                                     g
                                                                     ^{-1}(\mathrm{d}\phi,
                                                                     \mathrm{d}
                                                                     \phi)
                                                                     \text{div}(X))
                                                                     \mathrm{d} V.
\end{align*}
But the Lie derivative of a musical isomorphism is
\begin{align*}
\mathcal{L}_X Y ^{\flat}(Z) &= X Y ^{\flat}(Z) - Y ^{\flat}([X, Z]) = X
                              g(Y, Z) - g(Y, [X, Z])
  \\
                            &= g(\nabla_X Y, Z) + g(Y, \nabla_X Z - [X, Z])
  \\
  &= g(\nabla_X Y, Z) + g(Y, \nabla_Z X) =: (\nabla_XY) ^{\flat}(Z) +
    \theta_{Y}(Z),
\end{align*}
whence
\begin{align*}
\mathcal{L}_X g ^{-1}(Y ^{\flat}, Z ^{\flat}) &= X g ^{-1}(Y ^{\flat},
                                                Z ^{\flat}) - g
                                                ^{-1}(\mathcal{L}_X Y
                                                ^{\flat} , Z ^{\flat}) - g
                                                ^{-1}(Y ^{\flat}, \mathcal{L}_X Z
                                                ^{\flat} )
  \\
  &= X g (Y ,
                                                Z) - g
                                                ^{-1}((\nabla_X Y)
                                                ^{\flat} , Z ^{\flat}) - g
                                                ^{-1}(Y ^{\flat},
    (\nabla_X Z) ^{\flat} )
  \\
                                              &- g ^{-1}(\theta_{Y}, Z ^{\flat}) - g ^{-1}(Y ^{\flat}, \theta_{Z})
  \\
                                              &= - g(\theta_{Y} ^{\sharp}, Z) - g (Y , \theta_{Z} ^{\sharp})
  \\
                                              &= - \theta_Y(Z) - \theta_Z(Y)
  \\
  &= -g(Y, \nabla_Z X) - g(Z, \nabla_YX) = - \mathcal{L}_X(Y, Z).
\end{align*}
We can then write
$$
  \delta S = \int_M (\mathcal{L}_X g (\text{grad }\phi, \text{grad
  }\phi) - g(\text{grad }\phi, \text{grad }\phi) \text{div}(X))
  \mathrm{d} V,
$$
as required.
