Writing the Klein-Gordon in terms of the gamma matrices Hello I have a quick question regarding the Dirac gamma matrices or Dirac equation and the Klein-Gordon equation. Recall that the Klein-Gordon is given by the following 
\begin{equation}
\left(\Box + m^2\right)\Psi=0
\end{equation}
where $\Box$ is the D'Alembert operator.
Then, we have the Dirac equation written in terms of the gamma matrices (or covariant form)
\begin{equation}
\left( i\gamma^{\mu}\partial_{\mu}-m\right)\Psi=0
\end{equation}
I multiplied both sides of the Dirac equation by the $\beta$ matrix, such that I may be able to reduce the gamma matrices, $\gamma^{\mu}$, to the Dirac matrices, $\alpha^{i}$, after which I $\textbf{square}$ both sides of the equation, from which we have the Klein-Gordon equation after imposing certain restrictions on the alpha and beta matrices, i.e. 
\begin{equation}\left(-i \alpha_{x} \frac{\partial}{\partial x}-i \alpha_{y} \frac{\partial}{\partial y}-i \alpha_{z} \frac{\partial}{\partial z}+\beta m\right)\left(-i \alpha_{x} \frac{\partial}{\partial x}-i \alpha_{y} \frac{\partial}{\partial y}-i \alpha_{z} \frac{\partial}{\partial z}+\beta m\right) \psi=-\frac{\partial^{2} \psi}{\partial t^{2}}\end{equation}
\begin{equation}\begin{aligned}
\implies -\frac{\partial^{2} \psi}{\partial t^{2}}=&-\alpha_{x}^{2} \frac{\partial^{2} \psi}{\partial x^{2}}-\alpha_{y}^{2} \frac{\partial^{2} \psi}{\partial y^{2}}-\alpha_{z}^{2} \frac{\partial^{2} \psi}{\partial z^{2}}+\beta^{2} m^{2} \psi \\
&-\left(\alpha_{x} \alpha_{y}+\alpha_{y} \alpha_{x}\right) \frac{\partial^{2} \psi}{\partial x \partial y}-\left(\alpha_{y} \alpha_{z}+\alpha_{z} \alpha_{y}\right) \frac{\partial^{2} \psi}{\partial y \partial z}-\left(\alpha_{z} \alpha_{x}+\alpha_{x} \alpha_{z}\right) \frac{\partial^{2} \psi}{\partial z \partial x} \\
&-\left(\alpha_{x} \beta+\beta \alpha_{x}\right) m \frac{\partial \psi}{\partial x}-\left(\alpha_{y} \beta+\beta \alpha_{y}\right) m \frac{\partial \psi}{\partial y}-\left(\alpha_{z} \beta+\beta \alpha_{z}\right) m \frac{\partial \psi}{\partial z}
\end{aligned}\end{equation}
such that,
\begin{equation}\begin{aligned}
\alpha_{x}^{2}=\alpha_{y}^{2}=\alpha_{z}^{2}=\beta^{2} &=1 \\
\alpha_{j} \beta+\beta \alpha_{j} &=0 \\
\alpha_{j} \alpha_{k}+\alpha_{k} \alpha_{j} &=0 \quad(j \neq k)
\end{aligned}\end{equation}
Now, this works swell for the normal Dirac equation, but what about the Dirac equation which may have additional terms, i.e. the Dirac equation in curved spacetime, in which we have the spin connection, as well as other terms, i.e. the Levi-Civita term. This proves extremely inefficient, as multiplying both sides by the $\beta$ matrix, then squaring both sides gives an enormous amount of terms. Is there any other way in which perhaps we can derive the Klein-Gordon using gamma matrices, or perhaps derive the Klein-Gordon from the Dirac? Thank you for your help 
 A: Simply hit the Dirac equation with the complex conjugate of $(i\gamma^\mu \partial_\mu-m)$. So
\begin{align} (i\gamma^\mu \partial_\mu-m)\Psi &= 0 \\
-(i\gamma^\nu \partial_\nu+m)(i\gamma^\mu \partial_\mu-m)\Psi &= 0 \\
(\gamma^\mu\gamma^\nu\partial_\mu\partial_\nu+m^2)\Psi &= 0
\end{align}
Since both $\mu$ and $\nu$ are dummy indices it doesn't hurt to do a dumb move and rename them, so you get
\begin{align}
\left(\frac{1}{2}\big(\gamma^\mu\gamma^\nu\partial_\mu\partial_\nu+\gamma^\nu\gamma^\mu\partial_\nu\partial_\mu\big)+m^2\right)\Psi &= 0.
\end{align}
But now observe that $\partial_\mu\partial_\nu = \partial_\nu\partial_\mu$ so
\begin{align}
\left(\frac{1}{2}\big(\gamma^\mu\gamma^\nu+\gamma^\nu\gamma^\mu\big)\partial_\mu\partial_\nu+m^2\right)\Psi &= 0.
\end{align}
Finally the gamma matrices obey the Clifford algebra $\{\gamma^\mu,\gamma^\nu\}=2\eta^{\mu\nu}$, so the last equation becomes simply
\begin{align}
\left(\frac{1}{2}2\;\eta^{\mu\nu}\partial_\mu\partial_\nu+m^2\right)\Psi &= 0 \\
\left(\Box+m^2\right)\Psi &= 0,
\end{align}
which is the Klein-Gordon that you wanted.
A: For the curved space Dirac operator 
$$
{D}=  \gamma^a e^\mu_a \left(\partial_\mu + \textstyle{\frac 12} \sigma^{bc}\,\omega_{bc\mu}\right)= \gamma^a D_a\equiv \gamma^\mu D_\mu. 
$$
we have the simple looking Lichnerowicz formula
$$
D^2
= \nabla^2 -  \frac 14 R,
$$
where 
$$
\nabla^2\equiv\frac 1{\sqrt{g}} D_\mu \sqrt{g}g^{\mu\nu}  D_\nu
$$
is the "rough" or spin-connection Laplacian on spinors and $R$is the scalar curvature. Although named for Andre Lichnerowicz who published it in 1963 (A. Lichnerowicz, Spineurs harmonique,  Compt. Rend. Acad. Sci. Paris, Ser. A 257,  (1963) 7-9),
the     formula      appears as the very last equation in a paper  by Erwin Schroedinger writen some  30 years earlier (E. Schroedinger, Diracsches Elektron im Schwerefeld I, Sitzungsber. Preuss.  Akad.  Wiss., Phys.  Math. Kl. {11},  (1932) 105-128.). As you say, there are a lot of terms and quite a bit of manipulation (23 pages in Schroedinger's paper) to get the Schroedinger-Lichnerowicz formula! 
(Actually it's not that bad. You can do in a couple pf pages using the Bianchi identities.)
