-- Verify formulas for Compton scattering -- Run time about 30 seconds -- Safari recommended, click through alerts on other browsers E = sqrt(omega^2 + m^2) p1 = (omega, 0, 0, omega) p2 = (E, 0, 0, -omega) p3 = (omega, omega sin(theta) cos(phi), omega sin(theta) sin(phi), omega cos(theta)) p4 = (E, -omega sin(theta) cos(phi), -omega sin(theta) sin(phi), -omega cos(theta)) -- spinors N = (E + m)^2 u21 = (p2[1] + m, 0, p2[4], p2[2] + i p2[3]) u22 = (0, p2[1] + m, p2[2] - i p2[3], -p2[4]) u41 = (p4[1] + m, 0, p4[4], p4[2] + i p4[3]) u42 = (0, p4[1] + m, p4[2] - i p4[3], -p4[4]) u2 = (u21,u22) 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)) u4bar = dot(conj(u4),gamma0) -- adjoint of u4 q1 = p1 + p2 q2 = p4 - p1 qslash1 = dot(q1,gmunu,gamma) qslash2 = dot(q2,gmunu,gamma) a1 = zero(2,2,4,4) a2 = zero(2,2,4,4) for(a,1,2, for(b,1,2, a1[a,b] = -dot(u4bar[b], gammaT, qslash1 + m I, gammaT, u2[a]), a2[a,b] = -dot(u4bar[b], gammaT, qslash2 + m I, gammaT, u2[a]) )) "Sum over spin states" a11 = sum(a,1,2, sum(b,1,2, contract(dot(a1[a,b], gmunu, transpose(conj(a1[a,b])), gmunu)) )) a12 = sum(a,1,2, sum(b,1,2, contract(dot(a1[a,b], gmunu, conj(a2[a,b]), gmunu)) )) a22 = sum(a,1,2, sum(b,1,2, contract(dot(a2[a,b], gmunu, transpose(conj(a2[a,b])), gmunu)) )) S = a11 + a12 + conj(a12) + a22 S "Verify Casimir trick" pslash2 = dot(p2,gmunu,gamma) pslash4 = dot(p4,gmunu,gamma) P2 = pslash2 + m I P4 = pslash4 + m I Q1 = qslash1 + m I Q2 = qslash2 + m I T = dot(P2,gammaT,Q1,gammaT,P4,gammaL,Q1,gammaL) f11 = contract(contract(contract(T,3,4),2,3)) T = dot(P2,gammaT,Q2,gammaT,P4,gammaL,Q1,gammaL) f12 = contract(contract(contract(T,3,5),2,3)) T = dot(P2,gammaT,Q2,gammaT,P4,gammaL,Q2,gammaL) f22 = contract(contract(contract(T,3,4),2,3)) 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 = (s - m^2)^2 d12 = (s - m^2) (u - m^2) d22 = (u - m^2)^2 f = (f11/d11 + 2 f12/d12 + f22/d22) / 4 check(f11 == -8 s u + 24 s m^2 + 8 u m^2 + 8 m^4) check(f12 == 8 s m^2 + 8 u m^2 + 16 m^4) check(f22 == -8 s u + 8 s m^2 + 24 u m^2 + 8 m^4) p12 = dot(p1,gmunu,p2) p13 = dot(p1,gmunu,p3) p14 = dot(p1,gmunu,p4) check(f11 == 32 p12 p14 + 64 m^2 p12 - 32 m^2 p13 - 32 m^2 p14 + 32 m^4) check(f12 == 16 m^2 p12 - 16 m^2 p14 + 32 m^4) check(f22 == 32 p12 p14 + 32 m^2 p12 - 32 m^2 p13 - 64 m^2 p14 + 32 m^4) -- verify high energy approximation m = 0 m = 0 check(s == 4 omega^2) check(u == -2 omega^2 (cos(theta) + 1)) check(f == 2 ((cos(theta) + 1) / 2 + 2 / (cos(theta) + 1))) m = quote(m) -- undo m = 0 "ok" "Verify lab frame" -- rotate to lab frame (electron at rest) Lambda = ((E/m,0,0,omega/m),(0,1,0,0),(0,0,1,0),(omega/m,0,0,E/m)) p1 = dot(Lambda,p1) p2 = dot(Lambda,p2) p3 = dot(Lambda,p3) p4 = dot(Lambda,p4) check(s == dot(p1 + p2, gmunu, p1 + p2)) check(t == dot(p1 - p3, gmunu, p1 - p3)) check(u == dot(p1 - p4, gmunu, p1 - p4)) p12 = dot(p1,gmunu,p2) p13 = dot(p1,gmunu,p3) p14 = dot(p1,gmunu,p4) check(f11 == 32 p12 p14 + 64 m^2 p12 - 32 m^2 p13 - 32 m^2 p14 + 32 m^4) check(f12 == 16 m^2 p12 - 16 m^2 p14 + 32 m^4) check(f22 == 32 p12 p14 + 32 m^2 p12 - 32 m^2 p13 - 64 m^2 p14 + 32 m^4) -- verify s, t, and u in the lab frame omegaL = dot(p1, (1,0,0,0)) omegaLp = dot(p3, (1,0,0,0)) check(omegaL == omega^2 / m + omega E / m) check(omegaLp == omega^2 cos(theta) / m + omega E / m) check(s == m^2 + 2 m omegaL) check(t == 2 m (omegaLp - omegaL)) check(u == m^2 - 2 m omegaLp) check(f == 2 (omegaL/omegaLp + omegaLp/omegaL + (m/omegaL - m/omegaLp + 1)^2 - 1)) -- verify integral R = hbar omega / (m c^2) omegap = omega / (1 + R (1 - cos(theta))) f = (omegap / omega)^2 (omega / omegap + omegap / omega - sin(theta)^2) I = -cos(theta) / R^2 + log(1 + R (1 - cos(theta))) (1/R - 2/R^2 - 2/R^3) - 1 / (2 R (1 + R (1 - cos(theta)))^2) + 1 / (1 + R (1 - cos(theta))) (-2/R^2 - 1/R^3) check(f sin(theta) == d(I,theta)) "ok" Run