Eigenmath Demo

-- Verify Casimir trick for electron positron annihilation

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

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

p3 = (E,
      E expsin(theta) expcos(phi),
      E expsin(theta) expsin(phi),
      E expcos(theta))

p4 = (E,
      -E expsin(theta) expcos(phi),
      -E expsin(theta) expsin(phi),
      -E expcos(theta))

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)

N = (E + m)^2

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)

C1 = t - m^2
C2 = u - m^2

"Sum over spin states"

S = 0

for(s1,1,2,for(s2,1,2,
  a1 = dot(v2bar[s2],gammaT,qslash1 + m I,gammaT,u1[s1]),
  a2 = dot(v2bar[s2],gammaT,qslash2 + m I,gammaT,u1[s1]),
  a11 = contract(dot(a1,gmunu,transpose(conj(a1)),gmunu)),
  a12 = contract(dot(a1,gmunu,conj(a2),gmunu)),
  a22 = contract(dot(a2,gmunu,transpose(conj(a2)),gmunu)),
  f = a11/C1^2 + a12/(C1 C2) + conj(a12)/(C1 C2) + a22/C2^2,
  S = S + f
))

S = S / N -- normalize
S

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

T = f11/C1^2 + f12/(C1 C2) + conj(f12)/(C1 C2) + f22/C2^2
T

"Does S = T?"

test(S == T, "yes", "no")