This script can be pasted into the Eigenmath app.

 -- Verify Casimir trick for electron scattering E = sqrt(p^2 + m^2) p1 = (E, 0, 0, p) p2 = (E, 0, 0, -p) p3 = (E, p expsin(theta) expcos(phi), p expsin(theta) expsin(phi), p expcos(theta)) p4 = (E, -p expsin(theta) expcos(phi), -p expsin(theta) expsin(phi), -p expcos(theta)) u11 = (E + m, 0, p1[4], p1[2] + i p1[3]) u12 = (0, E + m, p1[2] - i p1[3], -p1[4]) u21 = (E + m, 0, p2[4], p2[2] + i p2[3]) u22 = (0, E + m, p2[2] - i p2[3], -p2[4]) u31 = (E + m, 0, p3[4], p3[2] + i p3[3]) u32 = (0, E + m, p3[2] - i p3[3], -p3[4]) u41 = (E + m, 0, p4[4], p4[2] + i p4[3]) u42 = (0, E + m, p4[2] - i p4[3], -p4[4]) u1 = (u11,u12) u2 = (u21,u22) u3 = (u31,u32) u4 = (u41,u42) N = (E + m)^4 I = ((1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)) gmunu = ((1,0,0,0),(0,-1,0,0),(0,0,-1,0),(0,0,0,-1)) gamma0 = ((1,0,0,0),(0,1,0,0),(0,0,-1,0),(0,0,0,-1)) gamma1 = ((0,0,0,1),(0,0,1,0),(0,-1,0,0),(-1,0,0,0)) gamma2 = ((0,0,0,-i),(0,0,i,0),(0,i,0,0),(-i,0,0,0)) gamma3 = ((0,0,1,0),(0,0,0,-1),(-1,0,0,0),(0,1,0,0)) gamma = (gamma0,gamma1,gamma2,gamma3) gammaT = transpose(gamma) gammaL = transpose(dot(gmunu,gamma)) u3bar = dot(conj(u3),gamma0) -- adjoint of u3 u4bar = dot(conj(u4),gamma0) -- adjoint of u4 "Sum over spin states" a11 = 0 a12 = 0 a22 = 0 for(a,1,2,for(b,1,2,for(c,1,2,for(d,1,2, a1 = dot(dot(u3bar[c],gammaT,u1[a]), dot(u4bar[d],gammaL,u2[b])), a2 = dot(dot(u4bar[d],gammaT,u1[a]), dot(u3bar[c],gammaL,u2[b])), a11 = a11 + a1 conj(a1), a12 = a12 + a1 conj(a2), a22 = a22 + a2 conj(a2) )))) a11 a12 a22 "Casimir trick" pslash1 = dot(p1,gmunu,gamma) pslash2 = dot(p2,gmunu,gamma) pslash3 = dot(p3,gmunu,gamma) pslash4 = dot(p4,gmunu,gamma) X1 = pslash1 + m I X2 = pslash2 + m I X3 = pslash3 + m I X4 = pslash4 + m I T1 = contract(dot(X3,gammaT,X1,gammaT),1,4) T2 = contract(dot(X4,gammaL,X2,gammaL),1,4) f11 = contract(dot(T1,transpose(T2))) T = contract(dot(X3,gammaT,X1,gammaT,X4,gammaL,X2,gammaL),1,6) f12 = contract(contract(T,1,3)) T1 = contract(dot(X4,gammaT,X1,gammaT),1,4) T2 = contract(dot(X3,gammaL,X2,gammaL),1,4) f22 = contract(dot(T1,transpose(T2))) f11 f12 f22 "Verfiy Casimir trick (1=ok)" f11 == 1/N a11 f12 == 1/N a12 f22 == 1/N a22