[proofplan]
We differentiate the Ricci tensor through its expression as the contraction of the curvature of the Levi-Civita connection. Under Ricci flow, the time derivative of the Christoffel symbols is a first covariant derivative of the Ricci tensor. Substituting this variation into the Ricci tensor variation gives second covariant derivatives of $\operatorname{Ric}$; the contracted Bianchi identity converts the trace-divergence terms into scalar-curvature Hessians, and the remaining commutators of covariant derivatives produce precisely the two quadratic curvature terms.
[/proofplan]
[step:Differentiate the Levi-Civita connection under Ricci flow]
Fix $t \in I$. All tensors, covariant derivatives, contractions, and index raising in this proof are taken with respect to the metric $g(t)$. Let $\nabla$ denote the Levi-Civita connection of $g(t)$, and let $\Gamma_{ij}^k$ denote its Christoffel symbols in the chosen coordinate chart, so that
\begin{align*}
\nabla_{\partial_{x_i}}\partial_{x_j}=\sum_{k=1}^{n}\Gamma_{ij}^k\partial_{x_k}.
\end{align*}
Define the metric variation tensor
\begin{align*}
h_{ij} := \partial_t g_{ij}.
\end{align*}
For a general metric variation, differentiating the identities $\nabla_i g_{j\ell}=0$ and the torsion-free symmetry $\Gamma_{ij}^k=\Gamma_{ji}^k$ gives
\begin{align*}
\partial_t\Gamma_{ij}^k
=
\frac{1}{2}g^{k\ell}
\left(
\nabla_i h_{j\ell}
+
\nabla_j h_{i\ell}
-
\nabla_\ell h_{ij}
\right).
\end{align*}
Since the Ricci flow equation gives $h_{ij}=-2\operatorname{Ric}_{ij}$, this becomes
\begin{align*}
\partial_t\Gamma_{ij}^k
=
-g^{k\ell}
\left(
\nabla_i\operatorname{Ric}_{j\ell}
+
\nabla_j\operatorname{Ric}_{i\ell}
-
\nabla_\ell\operatorname{Ric}_{ij}
\right).
\end{align*}
[guided]
We first need the infinitesimal change of the Levi-Civita connection. Let
\begin{align*}
h_{ij} := \partial_t g_{ij}
\end{align*}
be the metric variation tensor. Because $\nabla$ is the Levi-Civita connection of the fixed metric $g(t)$, it satisfies $\nabla_i g_{j\ell}=0$ and has symmetric Christoffel symbols $\Gamma_{ij}^k=\Gamma_{ji}^k$.
The standard variation formula for the Levi-Civita connection is obtained by differentiating metric compatibility and solving for $\partial_t\Gamma_{ij}^k$:
\begin{align*}
\partial_t\Gamma_{ij}^k
=
\frac{1}{2}g^{k\ell}
\left(
\nabla_i h_{j\ell}
+
\nabla_j h_{i\ell}
-
\nabla_\ell h_{ij}
\right).
\end{align*}
Here $g^{k\ell}$ denotes the components of the inverse metric $g(t)^{-1}$, and the covariant derivatives are taken using $\nabla$ at the fixed time $t$.
For Ricci flow, the variation tensor is not arbitrary: the equation $\partial_t g_{ij}=-2\operatorname{Ric}_{ij}$ says
\begin{align*}
h_{ij}=-2\operatorname{Ric}_{ij}.
\end{align*}
Substituting this into the general variation formula gives
\begin{align*}
\partial_t\Gamma_{ij}^k
=
\frac{1}{2}g^{k\ell}
\left(
-2\nabla_i\operatorname{Ric}_{j\ell}
-2\nabla_j\operatorname{Ric}_{i\ell}
+2\nabla_\ell\operatorname{Ric}_{ij}
\right).
\end{align*}
Factoring out the coefficient $-2$ and cancelling it with $1/2$ yields
\begin{align*}
\partial_t\Gamma_{ij}^k
=
-g^{k\ell}
\left(
\nabla_i\operatorname{Ric}_{j\ell}
+
\nabla_j\operatorname{Ric}_{i\ell}
-
\nabla_\ell\operatorname{Ric}_{ij}
\right).
\end{align*}
This is the key input: the time derivative of the connection is expressed entirely in terms of first covariant derivatives of the Ricci tensor.
[/guided]
[/step]
[step:Convert the connection variation into second derivatives of the Ricci tensor]
The Ricci tensor is the contraction $\operatorname{Ric}_{ij}=R^{k}_{\ ikj}$. Since the difference of two connections is tensorial, the variation of Ricci is
\begin{align*}
\partial_t\operatorname{Ric}_{ij}
=
\nabla_k(\partial_t\Gamma_{ij}^k)
-
\nabla_j(\partial_t\Gamma_{ik}^k).
\end{align*}
Using $\nabla g^{-1}=0$ and the formula from the previous step,
\begin{align*}
\nabla_k(\partial_t\Gamma_{ij}^k)
=
-\nabla^k\nabla_i\operatorname{Ric}_{jk}
-\nabla^k\nabla_j\operatorname{Ric}_{ik}
+
\Delta\operatorname{Ric}_{ij}.
\end{align*}
Also,
\begin{align*}
\partial_t\Gamma_{ik}^k
=
-g^{k\ell}
\left(
\nabla_i\operatorname{Ric}_{k\ell}
+
\nabla_k\operatorname{Ric}_{i\ell}
-
\nabla_\ell\operatorname{Ric}_{ik}
\right).
\end{align*}
The first traced term is $-\nabla_i R$, where $R=g^{k\ell}\operatorname{Ric}_{k\ell}$ is the scalar curvature, and the last two traced terms cancel after relabelling the contracted indices. Hence
\begin{align*}
\partial_t\Gamma_{ik}^k=-\nabla_i R.
\end{align*} Therefore
\begin{align*}
\partial_t\operatorname{Ric}_{ij}
=
\Delta\operatorname{Ric}_{ij}
-\nabla^k\nabla_i\operatorname{Ric}_{jk}
-\nabla^k\nabla_j\operatorname{Ric}_{ik}
+
\nabla_j\nabla_i R.
\end{align*}
[/step]
[step:Commute derivatives and use the contracted Bianchi identity]
The contracted Bianchi identity gives
\begin{align*}
\nabla^k\operatorname{Ric}_{jk}
=
\frac{1}{2}\nabla_j R.
\end{align*}
For the covariant $2$-tensor $\operatorname{Ric}$, the Ricci commutation identity gives
\begin{align*}
\nabla^k\nabla_i\operatorname{Ric}_{jk}
=
\nabla_i\nabla^k\operatorname{Ric}_{jk}
-
R_{ikj\ell}\operatorname{Ric}^{k\ell}
+
\operatorname{Ric}_{ik}\operatorname{Ric}^{k}_{\ j}.
\end{align*}
Using the contracted Bianchi identity in the first term on the right-hand side,
\begin{align*}
\nabla^k\nabla_i\operatorname{Ric}_{jk}
=
\frac{1}{2}\nabla_i\nabla_j R
-
R_{ikj\ell}\operatorname{Ric}^{k\ell}
+
\operatorname{Ric}_{ik}\operatorname{Ric}^{k}_{\ j}.
\end{align*}
Interchanging $i$ and $j$ gives the same formula:
\begin{align*}
\nabla^k\nabla_j\operatorname{Ric}_{ik}
=
\frac{1}{2}\nabla_j\nabla_i R
-
R_{jki\ell}\operatorname{Ric}^{k\ell}
+
\operatorname{Ric}_{jk}\operatorname{Ric}^{k}_{\ i}.
\end{align*}
Using the curvature symmetry $R_{jki\ell}=R_{ikj\ell}$ after contraction against the symmetric tensor $\operatorname{Ric}^{k\ell}$, and using $\operatorname{Ric}_{jk}\operatorname{Ric}^{k}_{\ i}=\operatorname{Ric}_{ik}\operatorname{Ric}^{k}_{\ j}$, we obtain
\begin{align*}
\nabla^k\nabla_i\operatorname{Ric}_{jk}
+
\nabla^k\nabla_j\operatorname{Ric}_{ik}
=
\nabla_i\nabla_j R
-
2R_{ikj\ell}\operatorname{Ric}^{k\ell}
+
2\operatorname{Ric}_{ik}\operatorname{Ric}^{k}_{\ j}.
\end{align*}
[/step]
[step:Substitute the commuted expression and cancel the scalar Hessian]
Substitute the identity from the previous step into
\begin{align*}
\partial_t\operatorname{Ric}_{ij}
=
\Delta\operatorname{Ric}_{ij}
-\nabla^k\nabla_i\operatorname{Ric}_{jk}
-\nabla^k\nabla_j\operatorname{Ric}_{ik}
+
\nabla_j\nabla_i R.
\end{align*}
Since $R$ is a scalar function, its Hessian is symmetric:
\begin{align*}
\nabla_j\nabla_i R=\nabla_i\nabla_j R.
\end{align*}
Therefore
\begin{align*}
\partial_t\operatorname{Ric}_{ij}
=
\Delta\operatorname{Ric}_{ij}
-
\left(
\nabla_i\nabla_j R
-
2R_{ikj\ell}\operatorname{Ric}^{k\ell}
+
2\operatorname{Ric}_{ik}\operatorname{Ric}^{k}_{\ j}
\right)
+
\nabla_i\nabla_j R.
\end{align*}
Cancelling the scalar Hessian terms gives
\begin{align*}
\partial_t\operatorname{Ric}_{ij}
=
\Delta\operatorname{Ric}_{ij}
+
2R_{ikj\ell}\operatorname{Ric}^{k\ell}
-
2\operatorname{Ric}_{ik}\operatorname{Ric}^{k}_{\ j}.
\end{align*}
This is the asserted coordinate evolution equation for the Ricci tensor under Ricci flow.
[/step]
