Let $\left(M^n, g\right)$ be an $n$ dimensional closed manifold, $T$ a tensor on this manifold, we define the energy functional to be $$ \mathcal{F}_T=\int_M|T|^2 d vol(g). $$ Let $g$ be a family of metrics defined as $g=g_0+t h$, where $h$ is a symmetric 2 tensor. The gradient of this energy functional will be denoted by grad $\mathcal{F}_T$ and defined as $$ \frac{\partial}{\partial t} \mathcal{F}_T=\int_M\left\langle\operatorname{grad} \mathcal{F}_T, h\right\rangle d \ vol(g). $$
Now, we want to compute the gradient of the following functional: $$ \mathcal{F}_{\text {Ric }}=\int_M|R i c|^2 d \ vol(g). $$
We want to show that the gradient of the above functional is the following:
$$ \operatorname{gradF}_{\text {Ric }}=-\Delta R_{i j}-\frac{1}{2} \Delta R g_{i j}+\nabla_i \nabla_j R-2 R^{p q} R_{i p j q}+\frac{1}{2}|R i c|^2 g_{i j}. $$
I have computed the variation of the Ricci tensor, metric tensor, and volume element. Could you help me understand the following computation? How to show the equality of (1)?
$$ \begin{aligned} \frac{\partial}{\partial t} \int_M|R i c|^2 d \mu=& \int_M\left\langle R_{i j}, \frac{\partial}{\partial t} R_{i j}\right\rangle d \mu+\int_M\left\langle\frac{1}{2}|R i c|^2 g_{i j}-2 R_{i p} R_j^p, h_{i j}\right\rangle dvol(g) ....... (1) \\ =& \int_M\left\langle R_{i j},-\Delta h_{i j}\right\rangle d vol(g) +\int_M\left\langle R_{i j},-\nabla_i \nabla_j T r_g(h)\right\rangle dvol(g) \\ &+\int_M\left\langle R_{i j}, \nabla^p \nabla_i h_{j p}\right\rangle dvol(g)+\int_M\left\langle R_{i j}, \nabla^p \nabla_j h_{i p}\right\rangle dvol(g) \\ &+\frac{1}{2} \int_M\left\langle|\operatorname{Ric}|^2 g_{i j}, h_{i j}\right\rangle dvol(g) \end{aligned} $$ Finally, how can I proceed from here? Somehow I can use the divergence theorem? Thanks so much.
Edit: I think I got an answer (not precise though). At the very end of the calculation, he uses the fact that the divergence and gradient are "dual" to each other.