-- Verify muon decay formulas p1x = 0 p1y = 0 p1z = 0 E1 = sqrt(p1x^2 + p1y^2 + p1z^2 + m1^2) E2 = sqrt(p2x^2 + p2y^2 + p2z^2 + m2^2) E3 = sqrt(p3x^2 + p3y^2 + p3z^2 + m3^2) E4 = sqrt(p4x^2 + p4y^2 + p4z^2 + m4^2) p1 = (E1, p1x, p1y, p1z) -- muon p2 = (E2, p2x, p2y, p2z) -- muon neutrino p3 = (E3, p3x, p3y, p3z) -- electron antineutrino p4 = (E4, p4x, p4y, p4z) -- electron u11 = (E1 + m1, 0, p1z, p1x + i p1y) / sqrt(E1 + m1) u12 = (0, E1 + m1, p1x - i p1y, -p1z) / sqrt(E1 + m1) u21 = (E2 + m2, 0, p2z, p2x + i p2y) / sqrt(E2 + m2) u22 = (0, E2 + m2, p2x - i p2y, -p2z) / sqrt(E2 + m2) v31 = (p3z, p3x + i p3y, E3 + m3, 0) / sqrt(E3 + m3) v32 = (p3x - i p3y, -p3z, 0, E3 + m3) / sqrt(E3 + m3) u41 = (E4 + m4, 0, p4z, p4x + i p4y) / sqrt(E4 + m4) u42 = (0, E4 + m4, p4x - i p4y, -p4z) / sqrt(E4 + m4) 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)) gamma5 = i dot(gamma0,gamma1,gamma2,gamma3) gamma = (gamma0,gamma1,gamma2,gamma3) gammaT = transpose(gamma) gammaL = transpose(dot(gmunu,gamma)) "Verify Casimir trick" u1 = (u11,u12) -- muon u2 = (u21,u22) -- muon neutrino v3 = (v31,v32) -- electron antineutrino u4 = (u41,u42) -- electron u2bar = dot(conj(u2),gamma0) -- adjoint of u2 u4bar = dot(conj(u4),gamma0) -- adjoint of u4 M(a,b,c,d) = dot(dot(u4bar[d],gammaT,I-gamma5,v3[c]),dot(u2bar[b],gammaL,I-gamma5,u1[a])) MM = sum(a,1,2,sum(b,1,2,sum(c,1,2,sum(d,1,2, M(a,b,c,d) conj(M(a,b,c,d)) )))) pslash1 = dot(p1,gmunu,gamma) pslash2 = dot(p2,gmunu,gamma) pslash3 = dot(p3,gmunu,gamma) pslash4 = dot(p4,gmunu,gamma) T1 = dot(pslash4,gammaT,I-gamma5,pslash3,gammaT,I-gamma5) T2 = dot(pslash2,gammaL,I-gamma5,pslash1,gammaL,I-gamma5) T1 = contract(T1,1,4) -- Tr T2 = contract(T2,1,4) -- Tr f = contract(dot(T1,transpose(T2))) check(f == MM) "ok" "Verify probability" check(f/4 == 64 dot(p1,gmunu,p3) dot(p2,gmunu,p4)) "ok" -- Compute muon lifetime GF = 1.1663787 10^(-5) GeV^(-2) mmu = 1.883531627 10^(-28) kilogram -- muon mass h = 6.62607015 10^(-34) joule second c = 299792458 meter / second joule = kilogram meter^2 / second^2 GeV = 10^9 1.602176634 10^(-19) joule kilogram = "kilogram" meter = "meter" second = "second" pi = float(pi) "Muon lifetime" 96 pi^2 h / (GF^2 (mmu c^2)^5) Run