A less experimental version of squarerootOfFourSigmaDotP.nb -- uses the results obtained from that code without actually taking the roots using the Eigensystem solution.

Defines:
* DiracMatrix (Weyl representation)
* complex -- allow for less general complex numbers, instead of: a + I b, complex[a, b] always formed from real a, b.
* diracMatrix (using complex instead of Complex)
* pauliMatrix (using complex instead of Complex)
* p · σ
* p · (OverBar[σ])
* (Sqrt[p·σ])
* (Sqrt[p·(OverBar[σ])])
* u, v, (OverBar[u]), (OverBar[v]),

ClearAll[complex, normsq, norm, complexQ, notComplexQ, real, imag, conjugate, complexI]
complex/:complex[r1_,i1_]+complex[r2_,i2_]:=complex[r1+r2,i1+i2];
complex/:r1_+complex[r2_,i2_]:=complex[r1+r2,i2];

complex/:-complex[re_,im_]:=complex[-re,-im];

complex/:complex[re_]:=re;
complex/:complex[re_,0]:=re;

complex/:complex[r1_,i1_] complex[r2_,i2_]:=complex[r1 r2-i1 i2,r1 i2+r2 i1];


(* a simple Complex-like implementation that allows for explicitly encoding real and imaginary without using Assumptions. *)
complex/:Power[z_complex,2]:=complex[z] complex[z];
complex/:Power[z_complex,n_]:=Module[{r=norm[z]^(n/2),theta=n ArcTan[z//First,z//Last]},r complex[Cos[theta],Sin[theta]]];

normsq[z_complex]:=((z//First)^2+(z//Last)^2);
normsq[v_]:=v^2;
norm[v_] := Sqrt[normsq[v]];

complexQ[z_complex]:=True;
complexQ[_]:=False;
notComplexQ[v_]:=Not[complexQ[v]];

complex/:(v_?notComplexQ) complex[re_,im_]:=complex[v re,v im];

real[z_complex]:=(z//First);
imag[z_complex]:=(z//Last);
real[ex_]:=ex;
imag[ex_]:=0;
conjugate[z_complex]:=complex[z//First,-z//Last];
conjugate[ex_]:=ex;

complexI:=complex[0,1];

ClearAll[z2,DiracMatrix]

z2 =ConstantArray[0,{2,2}];

DiracMatrix::usage = "DiracMatrix[0-3].";
DiracMatrix[k_Integer] :=  ({{z2, PauliMatrix[k]}, {-PauliMatrix[k], z2}} // ArrayFlatten);
DiracMatrix[0] = ({{z2, PauliMatrix[0]}, {PauliMatrix[0], z2}} // ArrayFlatten);
DiracMatrix[5] = I DiracMatrix[0].DiracMatrix[1].DiracMatrix[2].DiracMatrix[3] ;

ClearAll[diracMatrix]
diracMatrix::usage = "Like DiracMatrix, but uses complex instead of Complex";
diracMatrix[2] = (DiracMatrix[2] /. {I -> complexI,-I -> -complexI});
diracMatrix[k_] := DiracMatrix[k]

ClearAll[pauliMatrix]
pauliMatrix::usage="Like PauliMatrix, but uses complex instead of Complex";
pauliMatrix[2]=(PauliMatrix[2]/.{I->complexI,-I->-complexI});
pauliMatrix[k_]:=PauliMatrix[k]

ClearAll[pSlashMinusM]
pSlashMinusM::usage = "(Subscript[p, μ])(γ^μ) : parameters are the (Subscript[p, μ])'s.";
pSlashMinusM[w_, px_, py_, pz_, m_] := w diracMatrix[0] + px diracMatrix[1] + py diracMatrix[2] + pz diracMatrix[3] - m IdentityMatrix[4];





ClearAll[diracMatrix]
diracMatrix::usage = "Like DiracMatrix, but uses complex instead of Complex";
diracMatrix[2] = (DiracMatrix[2] /. {I -> complexI,-I -> -complexI});
diracMatrix[k_] := DiracMatrix[k]

ClearAll[pauliMatrix]
pauliMatrix::usage="Like PauliMatrix, but uses complex instead of Complex";
pauliMatrix[2]=(PauliMatrix[2]/.{I->complexI,-I->-complexI});
pauliMatrix[k_]:=PauliMatrix[k]



ClearAll[pdotsigma, pdotosigma]
pdotsigma::usage = "p · σ = (Subscript[p, 0]) - (OverVector[p]) · (OverVector[σ]) ; p's are upper index.";
pdotsigma[p0_, p1_, p2_, p3_] := p0 pauliMatrix[0] -p1 pauliMatrix[1] - p2 pauliMatrix[2] - p3 pauliMatrix[3] ;

pdotosigma::usage = "p · (OverBar[σ]) = (Subscript[p, 0]) + (OverVector[p]) · (OverVector[σ]) ; p's are upper index.";
pdotosigma[p0_, p1_, p2_, p3_] := pdotsigma[p0, -p1, -p2, -p3]

ClearAll[commutator, anticommutator]
commutator::usage = "commutator [a,b] = a b - b a, of two matrices a,b.";
commutator[a_, b_, sign_:-1] := (a.b +sign b.a);
anticommutator::usage = "anticommutator {a,b} = a b + b a, of two matrices a,b.";
anticommutator[a_, b_] := commutator[a,b,1];

ClearAll[sqrtpdotsigma, sqrtpdotosigma]
sqrtpdotsigma::usage = "returns: (Sqrt[p·σ]) ; params: (ω, (OverVector[p]), m) ; (m^2)= (ω^2)- ((OverVector[p])^2) ; p's are upper index";
sqrtpdotsigma[w_,px_,py_,pz_, m_] :=
(1/(Sqrt[2 w + 2 m])) {{w - pz + m, -px + complexI py},{-px - py complexI, w +pz + m}};
sqrtpdotsigma::usage = "returns: \!\(\*SqrtBox[\(p\\ ·\\ \*OverscriptBox[\(σ\), \(_\)]\)]\) ; params: (ω, (OverVector[p]), m) ; (m^2)= (ω^2)- ((OverVector[p])^2) ; p's are upper index";
sqrtpdotosigma[w_,px_,py_,pz_, m_] := sqrtpdotsigma[w, -px, -py, -pz, m];

ClearAll[u,v, eta, zeta]
eta[1] = {1,0};
eta[2] = {0,1};
zeta[1] = {1,0};
zeta[2] = {0,1};
u[w_,px_,py_,pz_,m_, k_] := {sqrtpdotsigma[w,px,py,pz, m].eta[k], sqrtpdotosigma[w,px,py,pz, m].eta[k]} // Flatten;
v[w_,px_,py_,pz_,m_,k_] := {sqrtpdotsigma[w,px,py,pz, m].eta[k], -sqrtpdotosigma[w,px,py,pz, m].eta[k]} // Flatten;


ClearAll[ubar,vbar];
ubar[w_,px_, py_, pz_,m_,k_] := Map[conjugate,u[w,px,py,pz,m,k]].diracMatrix[0];
vbar[w_,px_, py_, pz_,m_,k_] := Map[conjugate,v[w,px,py,pz,m,k]].diracMatrix[0];

ClearAll[ps, pos, rootps, rootpos, format]
format = {omega -> ω, pu1 -> "(p^1)", pu2 -> "(p^2)", pu3 -> "(p^3)"};
ps = pdotsigma[omega, pu1,pu2,pu3] ;
pos = pdotosigma[omega, pu1,pu2,pu3] ;
rootps = sqrtpdotsigma[omega, pu1,pu2,pu3,m] ;
rootpos = sqrtpdotosigma[omega, pu1,pu2,pu3,m] ;
Map[(MatrixForm[# /. format]&),{ps, pos, rootps, rootops}]// Column

simplify[mprod_] := Module[{s1,s2,s3, s4},
s1 = ExpandAll[mprod];
s2 = s1 /. omega^2 -> m^2 + pu1^2 + pu2^2 + pu3^2;
s3 = ExpandAll[s2 // FullSimplify];
s4 = s3 /. (pu1^2) -> -(pu2^2)-(pu3^2) -(m^2)+omega^2;
s4 // Simplify
]

checkps = simplify[rootps.rootps];
(checkps /. format) // MatrixForm
( (checkps - ps)/. format)// MatrixForm

checkpos = simplify[rootpos.rootpos];
(checkpos /. format) // MatrixForm
( (checkpos - pos)/. format)// MatrixForm

(*expect 2m, 2m, 0, 0*)
{((ubar[omega,pu1,pu2,pu3,m,1].u[omega,pu1,pu2,pu3,m,1])//simplify)/.format,
((ubar[omega,pu1,pu2,pu3,m,2].u[omega,pu1,pu2,pu3,m,2])//simplify)/.format,
((ubar[omega,pu1,pu2,pu3,m,1].u[omega,pu1,pu2,pu3,m,2])//simplify)/.format,
((ubar[omega,pu1,pu2,pu3,m,2].u[omega,pu1,pu2,pu3,m,1])//simplify)/.format} // Column


(*expect -2m, -2m, 0, 0*)
{((vbar[omega,pu1,pu2,pu3,m,1].v[omega,pu1,pu2,pu3,m,1])//simplify)/.format,
((vbar[omega,pu1,pu2,pu3,m,2].v[omega,pu1,pu2,pu3,m,2])//simplify)/.format,
((vbar[omega,pu1,pu2,pu3,m,1].v[omega,pu1,pu2,pu3,m,2])//simplify)/.format,
((vbar[omega,pu1,pu2,pu3,m,2].v[omega,pu1,pu2,pu3,m,1])//simplify)/.format} // Column

ClearAll[pSlashMinusM]
pSlashMinusM::usage = "(Subscript[p, μ])(γ^μ) : parameters are the (Subscript[p, μ])'s.";
pSlashMinusM[w_, px_, py_, pz_, m_] := w diracMatrix[0] + px diracMatrix[1] + py diracMatrix[2] + pz diracMatrix[3] - m IdentityMatrix[4];

(*Test: verify that these solutions work.  I've used (Subscript[p, k ])as the spatial vector coordinates so flipped signs for the (Subscript[p, μ])'s below.*)
{{
((pSlashMinusM[omega,-pu1,-pu2,-pu3,m].u[omega,pu1,pu2,pu3,m,1]//simplify)/.format)//MatrixForm,((pSlashMinusM[omega,-pu1,-pu2,-pu3,m].u[omega,pu1,pu2,pu3,m,2]//simplify)/.format)//MatrixForm},{((pSlashMinusM[omega,-pu1,-pu2,-pu3,-m].v[omega,pu1,pu2,pu3,m,1]//simplify)/.format)//MatrixForm,((pSlashMinusM[omega,-pu1,-pu2,-pu3,-m].v[omega,pu1,pu2,pu3,m,2]//simplify)/.format)//MatrixForm}
}//Grid

Table[Row[{Superscript["(Overscript[v, -])",i],Superscript[γ,k],Superscript[v,j]," = ", (vbar[omega,pu1,pu2,pu3,m,i].diracMatrix[k].v[omega,pu1,pu2,pu3,m,j]// simplify)/. format}],{i,1,2},{j,1,2},{k,0,3}]// Flatten // Column
Table[Row[{Superscript["(Overscript[u, -])",i],Superscript[γ,k],Superscript[u,j]," = ", (ubar[omega,pu1,pu2,pu3,m,i].diracMatrix[k].u[omega,pu1,pu2,pu3,m,j]// simplify)/. format}],{i,1,2},{j,1,2},{k,0,3}]// Flatten // Column

PauliMatrix[2]// MatrixForm

(rootps/. format) // MatrixForm
(rootpos/. format) // MatrixForm
((rootps.eta[1])/. format) // MatrixForm
((rootpos.eta[1])/. format) // MatrixForm
(u[omega,pu1,pu2,pu3,m,1]Sqrt[2 m + 2 omega]/. format) // MatrixForm
(DiracMatrix[0].DiracMatrix[1].u[omega,pu1,pu2,pu3,m,1]Sqrt[2 m + 2 omega]/. format) // MatrixForm

DiracMatrix[5]// MatrixForm
((MatrixForm[DiracMatrix[#].DiracMatrix[5] + DiracMatrix[5].DiracMatrix[#]])&/@ (Range[4]-1) ) // Column

(I/2) DiracMatrix[1].DiracMatrix[2] // MatrixForm

(1/4) DiracMatrix[0].DiracMatrix[3].DiracMatrix[1].DiracMatrix[2] // MatrixForm

(I/4)commutator[DiracMatrix[0],DiracMatrix[3]] // MatrixForm

Simplify[( Cosh[a/2] PauliMatrix[0] + Sinh[a/2] PauliMatrix[3]).(Cosh[a/2] PauliMatrix[0] + Sinh[a/2] PauliMatrix[3]), a ∈ Reals] // TrigToExp