-- Verify formulas for Bhabha 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)) v11 = (p1[4], p1[2] + i p1[3], E + m, 0) / sqrt(E + m) v12 = (p1[2] - i p1[3], -p1[4], 0, E + m) / sqrt(E + m) u21 = (E + m, 0, p2[4], p2[2] + i p2[3]) / sqrt(E + m) u22 = (0, E + m, p2[2] - i p2[3], -p2[4]) / sqrt(E + m) v31 = (p3[4], p3[2] + i p3[3], E + m, 0) / sqrt(E + m) v32 = (p3[2] - i p3[3], -p3[4], 0, E + m) / sqrt(E + m) u41 = (E + m, 0, p4[4], p4[2] + i p4[3]) / sqrt(E + m) u42 = (0, E + m, p4[2] - i p4[3], -p4[4]) / sqrt(E + m) 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)) "Verify Casimir trick" v1 = (v11,v12) u2 = (u21,u22) v3 = (v31,v32) u4 = (u41,u42) v1bar = dot(conj(v1),gamma0) -- adjoint of v1 u4bar = dot(conj(u4),gamma0) -- adjoint of u4 M1(a,b,c,d) = -e^2 dot(dot(v1bar[a],gammaT,v3[c]),dot(u4bar[d],gammaL,u2[b])) M2(a,b,c,d) = e^2 dot(dot(v1bar[a],gammaT,u2[b]),dot(u4bar[d],gammaL,v3[c])) M11 = sum(a,1,2,sum(b,1,2,sum(c,1,2,sum(d,1,2, M1(a,b,c,d) conj(M1(a,b,c,d)) )))) M12 = sum(a,1,2,sum(b,1,2,sum(c,1,2,sum(d,1,2, M1(a,b,c,d) conj(M2(a,b,c,d)) )))) M21 = sum(a,1,2,sum(b,1,2,sum(c,1,2,sum(d,1,2, M2(a,b,c,d) conj(M1(a,b,c,d)) )))) M22 = sum(a,1,2,sum(b,1,2,sum(c,1,2,sum(d,1,2, M2(a,b,c,d) conj(M2(a,b,c,d)) )))) check(M12 == M21) 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(X1,gammaT,X3,gammaT),1,4) T2 = contract(dot(X4,gammaL,X2,gammaL),1,4) f11 = contract(dot(T1,transpose(T2))) T = contract(dot(X1,gammaT,X2,gammaT,X4,gammaL,X3,gammaL),1,6) f12 = contract(T,1,3,1,2) T1 = contract(dot(X1,gammaT,X2,gammaT),1,4) T2 = contract(dot(X4,gammaL,X3,gammaL),1,4) f22 = contract(dot(T1,transpose(T2))) check(e^4 f11 == M11) check(-2 e^4 f12 == M12 + M21) check(e^4 f22 == M22) "ok" "Verify probability density" 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 p14^2 - 32 m^2 p12 + 32 m^2 p13 - 32 m^2 p14 - 32 m^4) check(f22 == 32 p13^2 + 32 p14^2 + 64 m^2 p12 + 64 m^4) s = dot(p1 + p2, gmunu, p1 + p2) t = dot(p1 - p3, gmunu, p1 - p3) u = dot(p1 - p4, gmunu, p1 - p4) check(f11 == 8 u^2 + 8 s^2 - 64 u m^2 - 64 s m^2 + 192 m^4) check(f12 == -8 u^2 + 64 u m^2 - 96 m^4) check(f22 == 8 u^2 + 8 t^2 - 64 u m^2 - 64 t m^2 + 192 m^4) d11 = t^2 d12 = s t d22 = s^2 f = e^4 (f11/d11 - 2 f12/d12 + f22/d22) / 4 -- high energy approximation m = 0 check(f11 == 8 u^2 + 8 s^2) check(f12 == -8 u^2) check(f22 == 8 u^2 + 8 t^2) check(s == 4 E^2) check(t == -2 E^2 (1 - cos(theta))) check(u == -2 E^2 (1 + cos(theta))) A = (2 (1 + cos(theta))^2 + 8) / (1 - cos(theta))^2 B = -2 (1 + cos(theta))^2 / (1 - cos(theta)) C = 1 + cos(theta)^2 check(f11 / t^2 == 4 A) check(-2 f12 / (s t) == 4 B) check(f22 / s^2 == 4 C) check(f == e^4 (A + B + C)) check(f == e^4 ((cos(theta)^2 + 3) / (cos(theta) - 1))^2) m = quote(m) -- verify integral f = (cos(theta)^2 + 3)^2 / (cos(theta) - 1)^2 I = 16 / (cos(theta) - 1) - 1/3 cos(theta)^3 - cos(theta)^2 - 9 cos(theta) - 16 log(1 - cos(theta)) check(f sin(theta) == d(I,theta)) "ok" Run