Riemannian manifold: Integration by parts of Lie derivative How can we apply integration by parts to the Lie derivative?
Background: In the Hamiltonian formulation of general relativity, we have the momentum constraint (using abstract index notation, $D_a$ is the covariant derivative on the 3-dimensional spatial slice with positive definite metric, $N^a$ is a vector field)
$$C_{\vec{N}}=2\int\mathrm{d}^3x\,(D_aN_b)p^{ab}=\int\mathrm{d}^3x\,(\mathcal{L}_{\vec{N}}q_{ab})p^{ab}=-\int\mathrm{d}^3x\,(\mathcal{L}_{\vec{N}}p^{ab})q_{ab}$$
In this situation, $p^{ab}$ is a symmetric (2,0) tensor density of weight 1 and $q_{ab}$ the spatial metric which is just a symmetric (0,2) tensor field. Now I could prove the given equalities exploiting the fact that $p^{ab}$ is density 1 and $q_{ab}$ is the metric. This means the statement
$$\int\mathrm{d}^3x\,(\mathcal{L}_{\vec{N}}q_{ab})P^{ab}=-\int\mathrm{d}^3x\,(\mathcal{L}_{\vec{N}}P^{ab})q_{ab}$$
holds in general for $q_{ab}$ being the metric and $P^{ab}$ being symmetric and of density 1. To prove this, I wrote everything in covariant derivatives and showed equality.
However, I'm wondering if there is a more geometrical way to see that this is true (maybe directly from Stokes' theorem). Furthermore, I'm wondering if the statement is always true or really just in the restricted situation that I could prove.
 A: 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.  
