-- Verify formulas for electron-positron annihilation

-- Run time about 30 seconds

-- Safari recommended, click through alerts on other browsers

E = sqrt(p^2 + m^2)

p1 = (E, 0, 0, p)
p2 = (E, 0, 0, -p)

p3 = (E,
      E sin(theta) cos(phi),
      E sin(theta) sin(phi),
      E cos(theta))

p4 = (E,
      -E sin(theta) cos(phi),
      -E sin(theta) sin(phi),
      -E cos(theta))

-- spinors

N = (E + m)^2

u11 = (E + m, 0, p1[4], p1[2] + i p1[3])
u12 = (0, E + m, p1[2] - i p1[3], -p1[4])

v21 = (p2[4], p2[2] + i p2[3], E + m, 0)
v22 = (p2[2] - i p2[3], -p2[4], 0, E + m)

u1 = (u11,u12)
v2 = (v21,v22)

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))

v2bar = dot(conj(v2),gamma0) -- adjoint of v2

q1 = p1 - p3
q2 = p1 - p4

qslash1 = dot(q1,gmunu,gamma)
qslash2 = dot(q2,gmunu,gamma)

t = dot(q1,gmunu,q1)
u = dot(q2,gmunu,q2)

a1 = zero(2,2,4,4)
a2 = zero(2,2,4,4)

for(a,1,2, for(b,1,2,
 a1[a,b] = dot(v2bar[b], gammaT, qslash1 + m I, gammaT, u1[a]),
 a2[a,b] = dot(v2bar[b], gammaT, qslash2 + m I, gammaT, u1[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"

pslash1 = dot(p1,gmunu,gamma)
pslash2 = dot(p2,gmunu,gamma)

P1 = pslash1 + m I
P2 = pslash2 - m I

Q1 = qslash1 + m I
Q2 = qslash2 + m I

T = dot(P1,gammaT,Q1,gammaT,P2,gammaL,Q1,gammaL)
f11 = contract(contract(contract(T,3,4),2,3))

T = dot(P1,gammaT,Q2,gammaT,P2,gammaL,Q1,gammaL)
f12 = contract(contract(contract(T,3,5),2,3))

T = dot(P1,gammaT,Q2,gammaT,P2,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 = (t - m^2)^2
d12 = (t - m^2) (u - m^2)
d22 = (u - m^2)^2

f = (f11/d11 + 2 f12/d12 + f22/d22) / 4

check(f11 == 8 t u - 24 t m^2 - 8 u m^2 - 8 m^4)
check(f12 == 8 s m^2 - 32 m^4)
check(f22 == 8 t u - 8 t 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 p13 p14 - 32 m^2 p12 + 64 m^2 p13 + 32 m^2 p14 - 64 m^4)
check(f12 == 16 m^2 p13 + 16 m^2 p14 - 32 m^4)
check(f22 == 32 p13 p14 - 32 m^2 p12 + 32 m^2 p13 + 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))

a = 1 + cos(theta)
b = 1 - cos(theta)

check(f == 2 (a/b + b/a))

m = quote(m) -- undo m = 0

-- verify integral

a = 1 + cos(theta)
b = 1 - cos(theta)
f = a/b + b/a
I = 2 cos(theta) + 2 log(1 - cos(theta)) - 2 log(1 + cos(theta))
check(f sin(theta) == d(I,theta))

"ok"