"Gordon decomposition" p1 = (E1, p1x, p1y, p1z) p2 = (E2, p2x, p2y, p2z) u11 = (E1 + m1, 0, p1z, p1x + i p1y) / sqrt(E1 + m1) -- spin up u12 = (0, E1 + m1, p1x - i p1y, -p1z) / sqrt(E1 + m1) -- spin down u21 = (E2 + m2, 0, p2z, p2x + i p2y) / sqrt(E2 + m2) -- spin up u22 = (0, E2 + m2, p2x - i p2y, -p2z) / sqrt(E2 + m2) -- spin down E1 = sqrt(p1x^2 + p1y^2 + p1z^2 + m1^2) E2 = sqrt(p2x^2 + p2y^2 + p2z^2 + m2^2) "Verify spinor normalization" check(dot(conj(u11),u11) == 2 E1) check(dot(conj(u12),u12) == 2 E1) check(dot(conj(u21),u21) == 2 E2) check(dot(conj(u22),u22) == 2 E2) "ok" 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) -- transpose first two indices to make compatible with dot function gammaT = transpose(gamma) T = dot(gamma,gammaT) sigmamunu = i/2 (T - transpose(T,1,3)) -- transpose mu and nu sigmamunu = transpose(sigmamunu,3,4) -- transpose nu and beta G = (outer(p1 + p2, I) + i dot(sigmamunu, gmunu, p2 - p1)) / (m1 + m2) -- transpose first two indices to make compatible with dot function G = transpose(G) "Verify G for spin up-up" u1 = u11 u2 = u21 u2bar = dot(conj(u2),gamma0) A = dot(u2bar, gammaT, u1) B = dot(u2bar, G, u1) check(A == B) "ok" "Verify G for spin up-down" u1 = u11 u2 = u22 u2bar = dot(conj(u2),gamma0) A = dot(u2bar, gammaT, u1) B = dot(u2bar, G, u1) check(A == B) "ok" "Verify G for spin down-up" u1 = u12 u2 = u21 u2bar = dot(conj(u2),gamma0) A = dot(u2bar, gammaT, u1) B = dot(u2bar, G, u1) check(A == B) "ok" "Verify G for spin down-down" u1 = u12 u2 = u22 u2bar = dot(conj(u2),gamma0) A = dot(u2bar, gammaT, u1) B = dot(u2bar, G, u1) check(A == B) "ok" Run