-- Verify formulas for Compton scattering 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)) 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) 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)) q1 = p1 + p2 q2 = p4 - p1 qslash1 = dot(q1,gmunu,gamma) qslash2 = dot(q2,gmunu,gamma) "Verify Casimir trick" u41bar = dot(conj(u41),gamma0) -- adjoint of u41 u42bar = dot(conj(u42),gamma0) -- adjoint of u42 M111 = dot(u41bar, -i e gammaT, qslash1 + m I, -i e gammaT, u21) M112 = dot(u41bar, -i e gammaT, qslash1 + m I, -i e gammaT, u22) M121 = dot(u42bar, -i e gammaT, qslash1 + m I, -i e gammaT, u21) M122 = dot(u42bar, -i e gammaT, qslash1 + m I, -i e gammaT, u22) M211 = dot(u41bar, -i e gammaT, qslash2 + m I, -i e gammaT, u21) M212 = dot(u41bar, -i e gammaT, qslash2 + m I, -i e gammaT, u22) M221 = dot(u42bar, -i e gammaT, qslash2 + m I, -i e gammaT, u21) M222 = dot(u42bar, -i e gammaT, qslash2 + m I, -i e gammaT, u22) M11 = contract(dot(M111, gmunu, transpose(conj(M111)), gmunu)) + contract(dot(M112, gmunu, transpose(conj(M112)), gmunu)) + contract(dot(M121, gmunu, transpose(conj(M121)), gmunu)) + contract(dot(M122, gmunu, transpose(conj(M122)), gmunu)) M12 = contract(dot(M111, gmunu, conj(M211), gmunu)) + contract(dot(M112, gmunu, conj(M212), gmunu)) + contract(dot(M121, gmunu, conj(M221), gmunu)) + contract(dot(M122, gmunu, conj(M222), gmunu)) M21 = contract(dot(M211, gmunu, conj(M111), gmunu)) + contract(dot(M212, gmunu, conj(M112), gmunu)) + contract(dot(M221, gmunu, conj(M121), gmunu)) + contract(dot(M222, gmunu, conj(M122), gmunu)) M22 = contract(dot(M211, gmunu, transpose(conj(M211)), gmunu)) + contract(dot(M212, gmunu, transpose(conj(M212)), gmunu)) + contract(dot(M221, gmunu, transpose(conj(M221)), gmunu)) + contract(dot(M222, gmunu, transpose(conj(M222)), gmunu)) check(M12 == M21) 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(T,3,4,2,3,1,2) T = dot(P2,gammaT,Q2,gammaT,P4,gammaL,Q1,gammaL) f12 = contract(T,3,5,2,3,1,2) T = dot(P2,gammaT,Q2,gammaT,P4,gammaL,Q2,gammaL) f22 = contract(T,3,4,2,3,1,2) 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 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) s = dot(p1 + p2, gmunu, p1 + p2) t = dot(p1 - p3, gmunu, p1 - p3) u = dot(p1 - p4, gmunu, p1 - p4) 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) d11 = (s - m^2)^2 d12 = (s - m^2) (u - m^2) d22 = (u - m^2)^2 f = e^4 (f11/d11 + 2 f12/d12 + f22/d22) / 4 -- high energy approximation m = 0 check(s == 4 omega^2) check(u == -2 omega^2 (cos(theta) + 1)) check(f == 2 e^4 ((cos(theta) + 1) / 2 + 2 / (cos(theta) + 1))) m = quote(m) "ok" "Verify lab frame" 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 e^4 (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