Dirac Operator in de Sitter background

Question

I am attempting to write the Dirac operator in a curved background, and eventually solve the equation as a second order PDE, since I am attempting to bring it into Klein-Gordon form. Essentially what I am trying to do is for the Dirac operator \begin{equation}D=\gamma^{a} e_{a}^{\mu}\left(\partial_{\mu}+\frac{1}{2} \sigma^{b c} \omega_{b c \mu}\right)=\gamma^{a} D_{a} \equiv \gamma^{\mu} D_{\mu}\end{equation} have the Lichnerowicz formula \begin{equation}D^{2}=\nabla^{2}-\frac{1}{4} R\end{equation} where \begin{equation}\nabla^{2} \equiv \frac{1}{\sqrt{g}} D_{\mu} \sqrt{g} g^{\mu \nu} D_{\nu}\end{equation} essentially, $D_{\mu}=\partial_{\mu}+\frac{1}{2} \sigma^{b c} \omega_{b c \mu}$, $\partial_{\mu}$ is the partial/directional derivative and $g=\det{g_{\mu \nu}}$. The operator would be applied to the spinor $\Psi$, such that I could essentially write the Dirac equation in Klein-Gordon form as \begin{equation}\left(D^{2} + m^{2}\right)\Psi=0\end{equation} So far, I have the following code

xu = {t, x, y, z};
n = Length[xu];
gd = {{-1, 0, 0, 
    0}, {0, (1 - K*(y^2 + z^2))/(1 - K*(x^2 + y^2 + z^2)), 
    K*x*y/(1 - K*(x^2 + y^2 + z^2)), 
    K*x*z/(1 - K*(x^2 + y^2 + z^2))}, {0, 
    K*x*y/(1 - K*(x^2 + y^2 + z^2)), (1 - K*(x^2 + z^2))/(1 - 
       K*(x^2 + y^2 + z^2)), K*y*z/(1 - K*(x^2 + y^2 + z^2))}, {0, 
    K*x*z/(1 - K*(x^2 + y^2 + z^2)), 
    K*y*z/(1 - K*(x^2 + y^2 + z^2)), (1 - K*(x^2 + y^2))/(1 - 
       K*(x^2 + y^2 + z^2))}}; (*Let this be the metric*)
gu = Simplify[Inverse[gd]];(*Let this be its inverse*)
\[CapitalGamma] = 1/2 Table[
    Sum[ gu[[\[Alpha], \[Beta]]] ( \!\(
\*SubscriptBox[\(\[PartialD]\), \(xu[[\[Mu]]]\)]\(gd[[\[Beta], \ \
\[Nu]]]\)\) +  \!\(
\*SubscriptBox[\(\[PartialD]\), \(xu[[\[Nu]]]\)]\(gd[[\[Beta], \ \
\[Mu]]]\)\) -  \!\(
\*SubscriptBox[\(\[PartialD]\), \(xu[[\[Beta]]]\)]\(gd[[\[Mu], \ \
\[Nu]]]\)\) ), {\[Beta], 1, n} ],
    {\[Alpha], 1, n}, {\[Mu], 1, n}, {\[Nu], 1, 
     n} ];(*Christoffel symbols/affine connection*)
Simplify[Flatten[\[CapitalGamma]]];
Eud = {{1, 0, 0, 
    0}, {0, (y^2 + z^2 + x^2 Sqrt[1/(1 - K (x^2 + y^2 + z^2))])/(x^2 +
        y^2 + z^2), (x y (-1 + 
         Sqrt[1/(1 - K (x^2 + y^2 + z^2))]))/(x^2 + y^2 + 
       z^2), (x z (-1 + Sqrt[1/(1 - K (x^2 + y^2 + z^2))]))/(x^2 + 
       y^2 + z^2)}, {0, (x y (-1 + 
         Sqrt[1/(1 - K (x^2 + y^2 + z^2))]))/(x^2 + y^2 + z^2), (x^2 +
        z^2 + y^2 Sqrt[1/(1 - K (x^2 + y^2 + z^2))])/(x^2 + y^2 + 
       z^2), (y z (-1 + Sqrt[1/(1 - K (x^2 + y^2 + z^2))]))/(x^2 + 
       y^2 + z^2)}, {0, (x z (-1 + 
         Sqrt[1/(1 - K (x^2 + y^2 + z^2))]))/(x^2 + y^2 + 
       z^2), (y z (-1 + Sqrt[1/(1 - K (x^2 + y^2 + z^2))]))/(x^2 + 
       y^2 + z^2), (x^2 + y^2 + 
       z^2 Sqrt[1/(1 - K (x^2 + y^2 + z^2))])/(x^2 + y^2 + 
       z^2)}} ;(*Let this be the tetrad/vierbein*)
Edu = Simplify[
   Inverse[Eud]]; (*The dual of the tetrad, i.e. the inverse vierbein*)
(\[Eta] = 
   Simplify[
    Table[Eud[[a]].gu.Eud[[b]], {a, 1, n}, {b, 1, n}]]) // MatrixForm
Simplify[Table[ Sum[Eud[[a, \[Mu]]] Edu[[b, \[Mu]]], {\[Mu], 1, n}], {a, 1, 
   n}, {b, 1, 
   n} ]] // MatrixForm (*This checks the duality condition of the \
tetrad and its dual*)
Omega := omega = Simplify[ Table[
      Sum[
      Edu[[a, \[Rho]]] Eud[[
        b, \[Nu]]] \[CapitalGamma][[\[Nu], \[Mu], \[Rho]]], {\[Rho], 
       1, n}, {\[Nu], 1, n}] 
     - Sum[Edu[[a, \[Nu]]] \!\(
\*SubscriptBox[\(\[PartialD]\), \(xu[[\[Mu]]]\)]\(Eud[[
         b, \ \[Nu]]]\)\), {\[Nu], 1, n}],
    {\[Mu], 1, n}, {b, 1, n}, {a, 1, 
     n} ]] (*spin connection of the above metric*)
listOmega := 
 Table[If[UnsameQ[Omega[[\[Mu], b, a]], 
    0], {ToString[\[CapitalOmega][\[Mu], b, a]], 
    Omega[[\[Mu], b, a]]}] , {\[Mu], 1, n}, {b, 1, n}, {a, 1, b}]
TableForm[Partition[DeleteCases[Flatten[listOmega], Null], 2], 
 TableSpacing -> {2, 2}]

I am a bit confused in writing the part involving the term $\sigma^{bc}=\frac{i}{2}[\gamma^{b},\gamma^{c}]$ in the spin connection. I am a bit at a loss in how to write the commutator, in order to define the operator which will inevitable act on the spinor. Must I use some sort of package such as FeynCalc or FeynRules? My ultimate goal is to define the Laplacian squared above, i.e. $\nabla^{2}$, which I will then use to derive the squared Dirac operator and esentially bring my Dirac equation to Klein-Gordon form. Any help would be greatly appreciated.