I'm reading a book named The Ricci Flow in Riemannian Geometry. In this book, the author defined a divergence operator:
$$ \begin{aligned} \delta_g: &\Gamma(\operatorname{Sym}^2T^*M)\mapsto \Gamma(T^*M)\\ (\delta_gh)_k&=-g^{ij}\nabla_ih_{jk}. \end{aligned} $$ where $g$ is the Riemannian metric.
Suppose $M$ is compact without boundary. Then they claimed that the adjoint of $\delta_g$ with respect to $L^2$-norm is $$ \begin{aligned} \delta^*_g: &\Gamma(T^*M)\mapsto \Gamma(\operatorname{Sym}^2T^*M)\\ (\delta^*_g\omega)_{jk}&=\frac{1}{2}(\nabla_j\omega_k+\nabla_k\omega_j)=\frac{1}{2}(\mathscr{L}_{\omega^{\sharp}}g)_{jk}. \end{aligned} $$
I try to prove this by verifying $\int_M(\delta_gh)^k\omega_k dvol_g=\int_Mh^{jk}(\delta^*\omega)_{jk}dvol_g$. I try to rewrite the expression under the intergal sign as a divergence of a vecor field then use Stoke's formula. But I failed. Can anyone help to explain the calculation? Thanks in advance.
Proof by brutal computations
Since $\nabla g=0$, it follows from Leibniz rule that \begin{align} (\delta h)^k\omega_k &= g^{kl}(\delta h)_l\omega_k\\ &=-g^{kl}g^{ij} \nabla_ih_{jl}\omega_k\\ &= -\nabla_i(g^{kl}g^{ij}h_{jl}\omega_k) + g^{kl}g^{ij}h_{jl}\nabla_i\omega_k\\ &= -\nabla_i(h^{ik}\omega_k) + h^{ik}\nabla_i\omega_k. \end{align} Since $M$ is closed, the divergence Theorem yields $\int_M \nabla_i(h^{ik}\omega_k)\mathrm{d}vol_g = 0$, so that $$ \langle \delta h,\omega\rangle_{L^2} = \int_Mh^{ik}\nabla_i\omega_k \mathrm{d}vol_g. $$ This does not allow to conclude, since $\nabla_i\omega_k$ is not symmetric. But here is the trick: since $h$ is symmetric, we have $$ h^{ik}\nabla_i\omega_k = \frac{1}{2}(h^{ik} + h^{ki})\nabla_i\omega_k = \frac{1}{2}h^{ik}\nabla_i\omega_k + \frac{1}{2}h^{ki}\nabla_i\omega_k, $$ and since we are dealing with dummy variables, we have that $$ h^{ik}\nabla_i\omega_k = \frac{1}{2}h^{ik}(\nabla_i\omega_k + \nabla_k\omega_i). $$ Finally, one has $$ \langle \delta h,\omega\rangle = \int_M h^{ik}\cdot \frac{1}{2}(\nabla_i\omega_k + \nabla_k\omega_i)\mathrm{d}vol_g, $$ which allows to conclude that $(\delta^*\omega)_{ik} = \frac{1}{2}(\nabla_i\omega_k + \nabla_k\omega_i)$.
Coordinate-free proof
Let $\sigma \colon \Gamma(T^*M\otimes T^*M) \to \Gamma(T^*M \otimes T^*M)$ be the symmetrisation endomorphism given by $$ \sigma(h)(X,Y) = \frac{1}{2}(h(X,Y) + h(Y,X)). $$ One readily checks that $\sigma$ is self-adjoint, that is $\langle \sigma(h),k\rangle_{L^2(T^*M\otimes T^*M)} = \langle h,\sigma(k)\rangle_{L^2(T^*M\otimes T^*M)}$ whenever $h$ and $k$ are $2$-tensors.
Now, recall that the divergence $\delta \colon \Gamma(T^*M\otimes T^*M) \to \Gamma(T^*M)$ is the adjoint operator of the covariant derivative $\nabla \colon \Gamma(T^*M) \to \Gamma(T^*M\otimes T^* M)$. This means that $\langle \delta h, \omega\rangle_{L^2(T^*M)} = \langle h, \nabla \omega \rangle_{L^2(T^*M\otimes T^* M)}$ whenever $h \in \Gamma(T^*M \otimes T^* M)$ and $\omega \in \Gamma(T^*M)$.
Ultimately, notice that if $h \in \Gamma(\mathrm{Sym}^2 T^*M)$, then $\sigma(h) = h$. It finally follows that if $h \in \Gamma(\mathrm{Sym}^2 T^*M)$ and $\omega \in \Gamma(T^*M)$, then \begin{align} \langle \delta h, \omega \rangle_{L^2(T^*M)} &= \langle h, \nabla \omega\rangle_{L^2(T^*M \otimes T^*M)} \\ &= \langle \sigma(h), \nabla \omega \rangle_{L^2(T^*M \otimes T^*M)} \\ &= \langle h, \sigma (\nabla \omega) \rangle_{L^2(T^*M \otimes T^*M)} \\ &= \langle h, (\sigma\circ \nabla) (\omega) \rangle_{L^2(T^*M \otimes T^*M)} \\ &= \langle h, (\sigma\circ \nabla) (\omega) \rangle_{L^2(\mathrm{Sym}^2T^*M)}. \end{align} Hence, $\sigma\circ \nabla \colon \Gamma(T^*M) \to \Gamma(\mathrm{Sym}^2T^*M)$ is the adjoint of $\delta \colon \Gamma(\mathrm{Sym}^2T^*M) \to \Gamma(T^*M)$.
Why is it equal to $\frac{1}{2}\mathcal{L}_{\omega^{\sharp}}g$?
Let $X = \omega^{\sharp}$, so that $\omega(Y) = g(X,Y)$ for all $Y$. Let $Y$ and $Z$ be two vector fields. Then \begin{align} (\delta^*\omega)(Y,Z) &= \frac{1}{2} \big((\nabla_Y\omega)(Z) + (\nabla_Z\omega)(Y)\big) \\ &= \frac{1}{2} \big(Y(\omega(Z)) - \omega(\nabla_YZ) + Z(\omega(Y)) - \omega(\nabla_ZY) \big) \\ &= \frac{1}{2}\big(Yg(X,Z) - g(X,\nabla_YZ) + Zg(X,Y) - g(X,\nabla_ZY) \big) \\ &= \frac{1}{2}\big(g(\nabla_YX,Z) + g(\nabla_ZX,Y) \big), \end{align} and one easily verifies that this last expression is equal to $(\frac{1}{2}\mathcal{L}_Xg)(Y,Z)$ (see, for instance, the first part of this old answer of mine).