-- Verify formulas for Moller scattering E = sqrt(p^2 + m^2) p1 = (E, 0, 0, p) p2 = (E, 0, 0, -p) p3 = (E, p sin(theta) cos(phi), p sin(theta) sin(phi), p cos(theta)) p4 = (E, -p sin(theta) cos(phi), -p sin(theta) sin(phi), -p cos(theta)) -- spinors N = (E + m)^4 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) 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 w31 = zero(2,2,4) w42 = zero(2,2,4) w41 = zero(2,2,4) w32 = zero(2,2,4) for(a,1,2, for(b,1,2, w31[a,b] = dot(u3bar[a], gammaT, u1[b]), w42[a,b] = dot(u4bar[a], gammaL, u2[b]), w41[a,b] = dot(u4bar[a], gammaT, u1[b]), w32[a,b] = dot(u3bar[a], gammaL, u2[b]) )) -- matrix elements a1 = zero(2,2,2,2) a2 = zero(2,2,2,2) for(a,1,2, for(b,1,2, for(c,1,2, for(d,1,2, a1[a,b,c,d] = dot(w31[c,a], w42[d,b]), a2[a,b,c,d] = dot(w41[d,a], w32[c,b]) )))) "Sum over spin states" a11 = sum(conj(a1) a1) a12 = sum(conj(a2) a1) a22 = sum(conj(a2) a2) S = a11 - a12 - conj(a12) + a22 S "Verify 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))) check(f11 == a11 / N) check(f12 == a12 / N) check(f22 == a22 / N) "ok" "Verify probability density" s = dot(p1 + p2, gmunu, p1 + p2) t = dot(p1 - p3, gmunu, p1 - p3) u = dot(p1 - p4, gmunu, p1 - p4) d11 = t^2 d12 = t u d22 = u^2 f = (f11/d11 - 2 f12/d12 + f22/d22) / 4 check(f11 == 8 s^2 + 8 u^2 - 64 s m^2 - 64 u m^2 + 192 m^4) check(f12 == -8 s^2 + 64 s m^2 - 96 m^4) check(f22 == 8 s^2 + 8 t^2 - 64 s m^2 - 64 t m^2 + 192 m^4) p12 = dot(p1, gmunu, p2) p13 = dot(p1, gmunu, p3) p14 = dot(p1, gmunu, p4) check(f11 == 32 p12^2 + 32 p14^2 - 64 m^2 p13 + 64 m^4) check(f12 == -32 p12^2 + 32 m^2 p12 + 32 m^2 p13 + 32 m^2 p14 - 32 m^4) check(f22 == 32 p12^2 + 32 p13^2 - 64 m^2 p14 + 64 m^4) -- verify for m = 0 m = 0 check(s == 4 E^2) check(t == 2 E^2 (cos(theta) - 1)) check(u == -2 E^2 (cos(theta) + 1)) check(f == 4 (cos(theta)^2 + 3)^2 / sin(theta)^4) -- Zee formula g = 2 ((1 + cos(theta/2)^4) / sin(theta/2)^4 + 2 / (sin(theta/2)^2 cos(theta/2)^2) + (1 + sin(theta/2)^4) / cos(theta/2)^4) check(f == g) m = quote(m) -- undo m = 0 -- verify integral f = (cos(theta)^2 + 3)^2 / sin(theta)^4 I = -8 cos(theta) / sin(theta)^2 - cos(theta) check(f sin(theta) == d(I,theta)) "ok" Run