How to show $(\delta S)Z = \delta (S(\cdot, Z))+\frac{1}{2}\langle S, \mathcal L_Zg\rangle$? I want  to show
$$
(\delta S)Z = \delta (S(\cdot, Z))+\frac{1}{2}\langle S, \mathcal L_Zg\rangle  
\tag{1}
$$
where $\delta= -tr_{12}\nabla$ is divergence operator, $\mathcal L$ is Lie derivative, $S\in\Gamma(Sym^2 T^*M)$
What I try: I use abstract index notation. First
$$
(\delta S)Z =-(\nabla^a S_{aj})Z =  -(\nabla^a S_{aj})Z^j
$$
Besides,
$$
\delta (S(\cdot,Z)) =-\nabla^a[S(\cdot, Z)]_a
$$
since $S(\cdot, Z)= S_{ij}Z^j$, I think $[S(\cdot, Z)]_a=S_{aj}Z^j$. Therefore
$$
\delta(S(\cdot, Z)) =-\nabla^a(S_{aj}Z^j)=-(\nabla ^a S_{aj})Z^j - S_{aj}(\nabla^a Z^j)
=(\delta S)Z- S_{aj}(\nabla^a Z^j)
\tag{2}
$$
On the other hand, when $\delta$ is restricted on $\Gamma (Sym^2 T^*M)$, its formal adjoint is
$\omega\rightarrow \frac{1}{2}\mathcal L_{\omega^\sharp}g$. Therefore, I have
$$
\frac{1}{2}\langle S, \mathcal L_Z g \rangle = \langle \delta(S), Z^\flat \rangle
=
\langle -\nabla^a S_{aj} , Z^i {g_{ik}} \rangle
=
-(\nabla^a S_{aj}) Z^j
\tag{3}
$$
However, I can't get (1) from (2) and (3).
If the answer is complex, a picture of handwritten draft is enough for me. Thanks very much.
 A: It's good that you are trying to use the abstract index notation, but you need one more piece for you puzzle, the formula for the Lie derivative of the metric:
$$
(\mathcal L_Z g)_{a b} = \nabla_a Z_b + \nabla_b Z_a = 2 \nabla_{(a} Z_{b)}
$$
where $Z^a$ is a vector field, and $Z_a = g_{a b} Z^b$ is the corresponding $1$-form. Of course, $\nabla$ is the Levi-Civita connection, corresponding to the Riemannian metric $g_{a b}$, and we use the latter along with its inverse $g^{a b}$ to raise and lower the indices without mention (the Ricci calculus).
With that in mind, the calculation will be straightforward:
$$
\delta (S(\cdot, Z)) = - \nabla^a (S_{a b} Z^b) = - (\nabla^a S_{a b})Z^b - S_{a b} \nabla^a Z^b = \\
 (\delta S) Z - S_{a b} \nabla^{(a} Z^{b)} = 
(\delta S) Z - \tfrac{1}{2} S_{a b} (\mathcal L_Z g)^{a b} = (\delta S) Z - \tfrac{1}{2} \langle  S, \mathcal L_Z g \rangle
$$
Here I use the notation for the symmetric part $t_{(a b)} := \tfrac{1}{2}(t_{a b} + t_{b a}) $ of a $2$-tensor $t_{a b}$ and the fact that if a tensor $s_{a b}$ is symmetric ($s_{a b} = s_{b a} = s_{(a b)} $), then for any tensor $t_{a b}$ we have the identity:
$$
s_{a b} t^{a b} = s_{(a b)} t^{a b} = s_{a b} t^{(a b)}
$$
