Do Lie derivatives commute with divergence: $\mathcal{L}_\xi \nabla_\mu V^\mu = \nabla_\mu \mathcal{L}_\xi V^{\mu}$? So my question is: given a vector field $V^\mu$ on some manifold with metric tensor $g_{\mu\nu}$ is it necessarily true that
\begin{equation}
\mathcal{L}_\xi \nabla_\mu V^\mu = \nabla_\mu \mathcal{L}_\xi V^{\mu} \ ?
\end{equation} 
where $\nabla_\mu$ denote the covariant derivative satisfying the metric compatibility condition $\nabla_\lambda g_{\mu\nu} = 0$ and $\mathcal{L}_\xi$ is the Lie derivative. If so how do I show that?
If not, then is there any condition we could impose on $\xi$ such that this is true? i.e. Is it true if $\xi$ is a Killing vector or something like that? 
 A: I think have got an answer but I'm not sure if it's right or not so I'm hoping that someone here could point out any potential mistake.
I think the conclusion is the identity is true if $\xi^\mu$ is a Killing vector. That means $\nabla_\nu \xi_\mu + \nabla_\mu \xi_\nu = 0$ and therefore $\nabla_\mu \xi^\mu = 0$. We have
\begin{align*}
\nabla_\mu \mathcal{L}_\xi V^\mu &= \nabla_\mu(\xi^\nu\nabla_\nu V^\mu - V^\nu \nabla_\nu \xi^\mu) \\
&= \xi^\nu \nabla_{\mu}\nabla_\nu V^\mu - V^\nu \nabla_\mu\nabla_\nu \xi^\mu \\
&= \xi^\nu\nabla_\nu\nabla_\mu V^\mu + \xi^\nu[\nabla_\mu, \nabla_\nu]V^\mu - V^\nu\nabla_\nu\nabla_\mu \xi^\mu - V^\nu[\nabla_\mu,\nabla_\nu]\xi^\mu \\
&= \xi^\nu\nabla_\nu\nabla_\mu V^\mu + \xi^\nu[\nabla_\mu, \nabla_\nu]V^\mu - V^\nu[\nabla_\mu,\nabla_\nu]\xi^\mu
\end{align*}
Now, using the identity
\begin{equation}
[\nabla_\mu, \nabla_\nu]v^\rho = R^\rho_{\ \sigma\mu\nu}v^\sigma
\end{equation}
where $R^\rho_{\ \sigma\mu\nu}$ is the Riemann curvature tensor. From this, we find that
\begin{equation}
\xi^\nu[\nabla_\mu,\nabla_\nu]V^\mu = R^\mu_{\ \sigma\mu\nu}\xi^\nu V^\sigma = R_{\sigma\nu}\xi^\nu V^\sigma
\end{equation}
and
\begin{equation}
V^\nu[\nabla_\mu,\nabla_\nu]\xi^\mu = R^{\mu}_{\ \sigma\mu\nu}\xi^\sigma V^\nu = R_{\sigma\nu}\xi^\sigma V^\nu = R_{\nu\sigma}\xi^\sigma V^\nu = \xi^\nu[\nabla_\mu,\nabla_\nu]V^\mu
\end{equation}
where $R_{\sigma\nu} = R^\mu_{\ \sigma\mu\nu}$ is the Ricci tensor satisfying a well-known symmetry property $R_{\mu\nu} = R_{\nu\mu}$. Therefore, the last two terms in the expression for $\nabla_\mu \mathcal{L}_\xi V^\mu$ will cancel. We are left with
\begin{equation}
\nabla_\mu \mathcal{L}_\xi V^\mu = \xi^\nu\nabla_\nu\nabla_\mu V^\mu = \xi^\nu\partial_\nu\nabla_\mu V^\mu = \mathcal{L}_\xi\nabla_\mu V^\mu
\end{equation} 
as required.
A: I was able to get $[L_\xi,\text{div}]X=0$ using coordinate invariant notation when $\xi$ is a Killing vector field.
Let $\epsilon$ be the volume form on the manifold. Then $L_X\epsilon = (\nabla\cdot X) \epsilon$. Then:
$$\nabla \cdot (L_\xi X) \epsilon = L_{L_\xi X} \epsilon = (L_\xi L_X - L_X L_\xi) \epsilon
\\
= L_\xi (\nabla \cdot X \epsilon) - L_X(\nabla \cdot \xi \epsilon)
\\
=L_\xi(\nabla \cdot X) \epsilon + \nabla \cdot X \nabla \cdot \xi \epsilon - L_X(\nabla \cdot \xi) \epsilon - \nabla \cdot \xi \nabla \cdot X \epsilon
$$
Since $\xi$ is Killing, its divergence vanishes everywhere, meaning everything but the first term on the last line vanishes. The volume form is nonzero everywhere, so $\nabla \cdot (L_\xi X)=L_\xi(\nabla \cdot X)$.
