This is the first time I face this problem, I think that the answer is positive **provided the connection is torsion free** as, in fact, happens for the Levi-Civita connection associated with a metric. Here is my proof.

Let us indicate the Lie derivative with respect to the vector field $Z$ by ${\cal L}_Z$. The action on tensor fields is completely determined by the following requirements.

[0] ${\cal L}_Z(X) = [Z,X]$.

[1] It trivially acts on scalar fields: ${\cal L}_Z f = Z(f)$.

[2] It commutes with contractions: ${\cal L}_Z(\langle X, \omega \rangle) = \langle {\cal L}_Z(X), \omega \rangle + \langle X, {\cal L}_Z(\omega) \rangle$.

[3] It acts as a derivative with respect to the tensor product: ${\cal L}_Z (T\otimes U) = ({\cal L}_Z(T) )\otimes U + T \otimes {\cal L}_Z(U)$.

The requirements [0] and [2] completely define the actions on covariant vector fields in terms of the action of vector fields on scalar fields (components of covariant fields):
$${\cal L}_Z(\omega) = \langle \partial_{x^a}, {\cal L}_Z(\omega) \rangle dx^a = Z(\langle \partial_{x^a}, \omega\rangle ) dx^a - \langle Z(\partial_{x^a}), \omega\rangle dx^a = Z(\omega_a) dx^a - \langle [Z, \partial_{x^a}], \omega \rangle dx^a\:,$$
i.e.,
$${\cal L}_Z(\omega)_a = Z(\omega_a) - \langle [Z, \partial_{x^a}], dx^b \rangle \omega_b\:. \tag{1}$$

From the definition of torsion we have that, if $\nabla$ is a torsion-free connection, then
$$[Z,X] = \nabla_Z X - \nabla_X Z$$

since the torsion is nothing but the difference of the two sides.
From [0] we conclude that
$${\cal L}_Z(X) = \nabla_Z X - \nabla_X Z\:, \tag{2}$$
or, in components:
$$({\cal L}_Z(X))^a = (Z^b\nabla_b X)^a - (\nabla_b Z)^aX^b\:,$$
and this identity just says that, for torsion-free connctions, **we can replace derivatives for covariant ones when computing the Lie derivative of contravariant vector fields.**

Let us pass to covariant vector fields. From (1) and (2) we have that
$${\cal L}_Z(\omega)_a = \nabla_Z(\omega_a) - \langle \nabla_Z( \partial_{x^a}), dx^b \rangle \omega_b + \langle \nabla_{\partial_{x^a}}(Z), dx^b \rangle \omega_b\:,$$
that is
$${\cal L}_Z(\omega)_a = (Z^b\nabla_b\omega)_a + \langle \nabla_{a}(Z), dx^b \rangle \omega_b\:,$$
namely
$${\cal L}_Z(\omega)_a = (Z^b\nabla_b\omega)_a + (\nabla_{a}Z)^b\omega_b\:,$$
and this identity just says that, for torsion-free connctions, **we can replace derivatives for covariant ones when computing the Lie derivative of covariant vector fields.**

Since generic tensor fields can be constructed as a linear combination of tensor products of contravariant and covariant vector fields, from [3] the found results can easily be generalized to all types of tensor fields. In computing the Lie derivative of a tensor field we can always replace the standard derivative for the covariant one at every occurrence provided the connection is torsion free.

Nice question. One way to think about it is that given a metric $g$, the statement $\mathcal L_Xg = 0$ says something about the *metric*, whereas $\nabla_Xg = 0$ says something about the *connection*. Now what $\mathcal L_Xg = 0$ says, is that the flow of $X$, where defined, is an isometry for the metric, while $\nabla_Xg = 0$ says that $\nabla$ transports a pair of tangent vectors along the integral curves of $X$ in such a way that their inner product remains the same.

As an example, consider the upper half plane model of the hyperbolic plane. Its metric is $y^{-2}(dx^2 + dy^2)$, so clearly $\partial_x$ is a Killing vector field; its flow, horizontal translation, is an isometry. The fact that $\nabla_{\partial_x}g = 0$ doesn't say anything about $g$, but it does say that Euclidean parallel transport is compatible with this directional derivative of the connection.

Now consider $\partial_y$. This of course is not a Killing vector field, since vertical translation is not an isometry. The connection however can be made such (by the theorem of Levi-Civita) that a pair of tangent vectors can be parallel transported in such a way that the inner product is preserved.

**EDIT**

I think I have a more illustrative example: consider the sphere embedded in $\Bbb R^3$. Pick an axis and take the velocity vector field $\xi$ associated to rotation around the axis at some constant angular velocity. Also consider a second vector field $\zeta$ that is everywhere (in a neighbourhood of the equator, extend in any smooth way toward the poles) proportional to $\xi$, but that has constant velocity everywhere, something like in this image

(downloaded from this page).

Obviously $\xi$ is a Killing field, as it integrates to an isometry. An immediate way to see that $\zeta$ is not, is by noting that curves parallel to the equator remain parallel to the equator under the flow of $\zeta$, hence so do their tangent vectors. What happens to a curve whose tangent vector at the equator points toward a pole, is that the flow of $\zeta$ moves the point at the equator over a smaller angle than a point above the equator, so these two vectors don't remain perpendicular. For parallel transport on the other hand, two perpendicular tangent vectors to a point at the equator will remain perpendicular both under $\xi$ and in $\zeta$, since they only depend on the restriction to the vector fields to the equator, where they are equal. This doesn't say anything about the vector field generating an isometry, i.e. being a Killing vector field.

## Best Answer

This is true in general, and there is a very nice geometrical reason why.

First use that the Lie derivative satisfies the Leibniz rule, $$£_N(q_{ab} p^{ab})=(£_Nq_{ab})p^{ab}+q_{ab}£_Np^{ab} $$ to rewrite the integral as $$\int d^3x (£_N q_{ab})p^{ab}= \int d^3x\,£_N(q_{ab}p^{ab}) - \int d^3x\,(£_N p^{ab})q_{ab} $$ Now note that the first integrand on the RHS is the Lie derivative of a scalar density. Scalar densities of weight 1 are more naturally thought of as a differential form of rank $d$, where $d$ is the dimension of the surface you are integrating over, so $3$ in this case. The 3-form in this case is $$p^{ab}q_{ab}\tilde{\epsilon}_{ijk} $$ where $\tilde{\epsilon}_{ijk}$ is the alternating symbol. The above expression is a tensor because $\tilde{\epsilon}$ is a density of weight -1, and $p^{ab}$ is weight 1. Now we apply Cartan's magic formula for differential forms, which says $$ £_N = i_N d + d i_N,$$ where $d$ is the exterior derivative, and $i_N$ denotes contraction with the vector $N$. So when the Lie derivative acts on a $d$ form, the $i_N d$ term vanishes because the exterior derivative of a $d$ form is a $d+1$ form, which vanishes for a dimension $d$ manifold. So the integral we are focusing on becomes $$ \int_S di_N(p^{ab}q_{ab} \tilde{\epsilon}) = \int_{\partial S}i_N(p^{ab}q_{ab}\tilde{\epsilon}).$$ Here we have used Stoke's theorem to relate the integral of an exact form to an integral of a form on the boundary of the region of integration. If the surface $S$ has no boundary, or if you assume the fields die off at the boundary, your identity results.

It should be clear that this argument will apply for any kind of situation like this: the Leibniz rule allows you to integrate by parts, and then by viewing the extra integral as a differential form you can show that it is always a boundary contribution.