Clifford Algebras and Spinors III: Bochner identity

Let $M$ be a Riemannian manifold, and let $Cl(M)$ be its Clifford bundle. Let $E \to M$ be any vector bundle with connection, and assume that $C^\infty(M, E)$ is a $Cl(M)$-module. We can define a Dirac operator $\mathcal{D}$ acting on sections of $E$ via the formula
\[ \mathcal{D} \sigma = \sum_{i=1}^n e_i \cdot \nabla_i \sigma \]
for any orthonormal frame $\{e_1, \dots, e_n\}$ on $M$, and where $\cdot$ denotes the Clifford module action. We demand that the connection on $E$ is compatible with Clifford multiplication in the following sense:
\[ \nabla_j (e_i \cdot \sigma) = (\nabla_j e_i) \cdot \sigma + e_i \cdot \nabla_j \sigma. \]

Let $R$ denote the curvature of $E$, i.e. we have
\[ [\nabla_i, \nabla_j] \sigma =  R(e_i, e_j) \sigma+ \nabla_{[e_i, e_j]} \sigma \]

We can define an endomorphism $\mathcal{R}$ on $E$ by
\[ \mathcal{R} = \frac{1}{2} \sum_{ij} R(e_i, e_j). \]

Theorem. We have the identity $\mathcal{D}^2 = -\Delta + \mathcal{R}$.

Proof. We compute
\begin{align}
\mathcal{D}^2 \sigma &= \sum_{ij} e_i \nabla_i \left( e_j \nabla_j \sigma \right) \\
&= \sum_{ij} e_i e_j \nabla_i \nabla_j \sigma + e_i ( \nabla_i e_j ) \nabla_j \sigma \\
&= -\Delta \sigma + \frac{1}{2}\sum_{ij}e_i e_j [\nabla_i, \nabla_j] \sigma + \sum_{ij} e_i ( \nabla_i e_j) \nabla_j \sigma \\
&= -\Delta \sigma + \frac{1}{2}\sum_{ij}e_i e_j R(e_i, e_j) \sigma + \frac{1}{2}\sum_{ij} e_i e_j \nabla_{[e_i, e_j]} \sigma+ \sum_{ij} e_i ( \nabla_i e_j) \nabla_j \sigma \\
&= (-\Delta + \mathcal{R})\sigma + \frac{1}{2} \sum_{ij} \left( e_i e_j \nabla_{[e_i, e_j]}\sigma +  e_i (\nabla_i e_j) \nabla_j + e_j (\nabla_j e_i) \nabla_i \right)\sigma
\end{align}
We will be done provided we can show that the last term vanishes. Notice that it is fully tensorial, since it can be expressed as $\mathcal{D}^2 + \Delta - \mathcal{R}$. On the other hand, the terms $[e_i, e_j]$ and $\nabla_j e_i$ are (by definition!) proportional to Christoffel symbols. Since we can always choose a frame so that these vanish at a point, these terms must vanish identically. Hence we have $0 = \mathcal{D}^2 + \Delta - \mathcal{R}$, as desired.