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.