Differential operator in Mathematica

I am fairly new to Mathematica and I wanted to know if I can define the following operator

$E:=(x^2sqrt{1-K(x^2+y^2+z^2}+y^2+z^2)gamma^{1}partial_{1}+(y^2sqrt{1-K(x^2+y^2+z^2}+x^2+z^2)gamma^{2}partial_{2}$

where the gamma’s are actually matrices, but that isn’t pertinent to the calculation; what is, however, is the position in which they are since they are three possible values that they can assume when multiplied together:

$mu=nu=0 rightarrow -2$

$mu=nuin {1,2,3} rightarrow 2$

$muneq nu in{0,1,2,3} rightarrow 0 $

since they satisfy the anticommutation relation, ${gamma^{mu},gamma^{nu}}=gamma^{mu}gamma^{nu}+gamma^{nu}gamma^{mu}=2eta^{munu}I_{4}$, such that I can then square the operator and apply it to some function $psi(t,x,y,z)$, $E^2psi$. I’ve been trying to use Nest and have even been looking at other posts for some ideas, but I can’t seem to wrap my finger on how I can get this to work. Any help????

https://mathematica.stackexchange.com/questions/72433/polynomial-expansion-of-operatorfunct