Commutator of Laplacian and a vector field on a Riemannian maniofld The following is from Taylor's PDE book, Vol. I, Chapter 2, Section 4, Equation 4.13. Also see exercises 1 and 2 of the same section.
At the end of that section, in proposition 4.2, it is stated that:

[In a Riemannian manifold $(M,g)$, equipped with the Levi-Civita connection,] a vector field $X$ commutes with [the Laplacian] $\Delta$ if and only if $X$ generates a group of isometries.

For a proof Taylor refers the reader to the following identity (in his notation $D_j=\partial/\partial x_j$, $X=X^jD_j$, and $\nabla_{D_j}X=X^{k}_{\:\:;j}D_k$ is the covariant derivative):
$$
[\Delta, X]u = (X^{j;k}+X^{k;j})u_{;j;k}+(X^{j;k}+X^{k;j})_{;j}u_{;k}
$$
where $u$ is an arbitrary scalar function. The identity is left to the reader, but I'm stomped proving it. At this point I'm not even certain this is true. In general, is there a source exploring the relation between Killing vector fields and Laplacian?
 A: For the next poor soul who spent a lot of time trying in vain to prove this statement, a statement which unfortunately is not true, I will leave the required modification and the sketch of the proof here.
I'll write the correct identity in a coordinate-free manner. Suppose $X$ is an arbitrary vector field. I will use the notation $\nabla_\mu T^{\cdots}_{\cdots}\equiv T^{\cdots}_{\cdots\:; \mu}$ for the covariant derivatives. Define the $(2,0)$ tensor $\mathrm{Def}(X)=\frac{1}{2}(\nabla^\mu X^{\nu}+\nabla^\nu X^{\mu})$. Then
$$
\frac{1}{2}[\Delta, X]u = \langle \mathrm{Def}(X), \nabla^2 u\rangle +
\langle \mathrm{div}\: \mathrm{Def}(X), du\rangle - \frac{1}{2}\langle\mathrm{grad}\:\mathrm{div}\:X, du\rangle
$$
The last term is what Taylor is missing (either that, or we need to add the volume-preserving assumption, i.e. $\mathrm{div}\: X=0$; which is reasonable given the context).
As for the proof, we first extend the definition of Laplacian to act on a tensor $T$ via $\Delta T = \mathrm{Tr}\nabla^2T=g^{ij}\nabla^2T(\partial_i,\partial_j)$. With this, it is straightfoward to verify
$$
[\Delta , X]u = 2\langle \mathrm{Def}(X), \nabla^2 u\rangle + \langle \Delta X, du\rangle + \langle X, \Delta du-d\Delta u)
$$
To go further, we calculate $\mathrm{div} \:\mathrm{Def}(X)$. This is obtained below
$$
2\mathrm{div} \:\mathrm{Def}(X)=
(\nabla_\mu \nabla^\mu X^\lambda +g^{\nu\lambda}\nabla_\mu\nabla_\nu X^\mu)\partial_\lambda=
\Delta X+ \mathrm{grad}\: \mathrm{div} X+g^{\nu\lambda}\mathrm{Ric}(X,\partial_\nu)\partial_\lambda
$$
Here, we used the commutator $[\nabla_\mu, \nabla_\nu]$ and use its relation to curvature tensor/Ricci tensor. Let $V=\mathrm{grad}\: u$. With this, one computes
$$
\langle \Delta X, du\rangle=
2\langle\mathrm{div} \mathrm{Def}(X),du\rangle- \langle\mathrm{grad}\: \mathrm{div} X, du\rangle-\mathrm{Ric}(X,V)
$$
Similarly, we prove $\langle X, \Delta du - d\Delta u\rangle = \mathrm{Ric}(V, X)$. Finally, as Ricci tensor is symmetric we are done.
A: As I commented I think that Taylor's claim is correct (and your argument shows that). Here's another, less computation-heavy argument: Let $(Φ_t)_t$ be the flow of $X$ (for simplicity lets pretend that it's defined on some open interval around $t = 0$). Then
$$
 [X, Δ_g] = \left.\frac{d}{d t}\right|_{t=0} Φ_t^* Δ_g
$$
so $[X, Δ] = 0$ is equivalent to $Φ_t^* Δ_g = Δ_g$ for all $t$. But the Laplacian is intrisically defined by the metric $g$, so $Φ^* Δ_g = Δ_{Φ^* g}$ for any diffeomorphism $Φ$. The claim now follows from the fact that the differential $Δ_g$ operator uniquely determines the metric $g$.
FAor that last fact is suffices to consider the case of $\mathbb{R}^n$ with some metric tensor given by $g_{ij}$. Recall that for a function $u$ with $u(0) = 0$ we have $u_{,i;j}(0) = u_{,ij}(0)$. Therefore $Δ_g u(0) = g^{ij} u_{,i;j}(0) = g^{ij} u_{,ij}(0)$. So if we plug in $u = x^i x^j$ we get
$$
Δ_g u(0) = g^{ij}(0)~~\text{for } i \neq j, \qquad Δ_g u(0) = 2 g^{ij}(0) ~~\text{for } i = j
$$
In any case we see that the differential operator $Δ_g$ determines the metric tensor at every point.
