/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 19.05.7 ] */
/* [wxMaxima: title   start ]
Helicity Amplitudes and Differential Scattering 
Cross Section for W boson Decay
    W(-, p1,λ)  --> e(-,p3,σ3) nu-bar(p4,σ4 = 1)
   [wxMaxima: title   end   ] */


/* [wxMaxima: comment start ]
W-decay.wxm,  July 26, 2019
Edwin (Ted) Woollett, Maxima by Example,  
ch. 12,   Dirac Algebra and Quantum Electrodynamics, ver. 3, "Dirac3"
 http://web.csulb.edu/~woollett/ 
 We use Heaviside-Lorentz electromagnetic units.
In general, we follow the conventions in Peskin & Schroeder: Quantum Field Theory.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Reference: Mark Thomson, Modern Particle Physics, Sec. 15.1.1

From Thomson (15.2), for given W^{-} polarization λ, electron helicity σ3, and 
electron anti-neutrino helicity σ4, the invariant amplitude is:

Mfi = (gw/sqrt(2)) ε(p1,λ)_{μ} ubar (p3, σ3) γ^{μ} (1/2) ( 1 - γ^5 ) v (p4, σ4)

The electron neutrino is left-handed and the electron anti-neutrino is
right handed (ignoring any possible non-zero mass for the neutrino),
so σ4 = +1.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
load(dirac3);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: section start ]
Particle 4-momentum Components in the  W Boson Rest Frame
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
In the rest frame of the initial W^{-}, (the center of momentum frame) we define the 
components of the particle 4-momenta using comp_def (pa (Ea, pax, pay, paz)) 
syntax (see next inputs).

With M standing for the W mass and m standing for the electron mass,
we use the fact that M >> m to ignore the electron mass. Then both the
final electron and antineutrino have the same 3-momentum magnitude p
and the same energy E = p.

The initial energy of the system is M and the final energy of the 
system is 2 E = M, by conservation of energy,  so p =  E = M/2.

p1 = W boson 4-momentum, p3 = electron 4-momentum, p4 = antineutrino 4-momentum.
The electron emerges in the "forward" direction, making an angle θ with the positive
z axis, and the antineutrino has  3-momentum in the opposite direction.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
assume (E > 0,  th > 0, th < %pi)$
comp_def ( p1( M,0,0,0),             
             p3 (E,E*sin(th),0,E*cos(th)),
             p4 (E,-E*sin(th),0,-E*cos(th)) )$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
listarray (p4);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
p4[0];
/* [wxMaxima: input   end   ] */


/* [wxMaxima: subsect start ]
Polarization 4-vectors for the W Boson at Rest
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
Polarization 4-vectors for the W boson *at rest*  (See Thomson ref. above)
See our file massive-spin1-pol.wxmx for derivations.

ep1RH is the polarization 4-vector describing a right-handed circularly 
polarized spin 1 massive particle with Sz = +1,
ep1LH  is the polarization 4-vector describing a left-handed circularly 
polarized spin 1 massive particle, with Sz = -1, 
and  ep1L  is the "longitudinal" polarization 4-vector describing a spin 1 
massive particle,  with Sz = 0. 

These three polarization vectors respectively stand for the 
states |S, Sz> = |1, 1>, |1, -1>, |1, 0>.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
comp_def ( ep1RH (0, -1/sqrt(2),  -%i/sqrt(2), 0 ),
                    ep1LH (0,  1/sqrt(2),  -%i/sqrt(2),  0),
                    ep1L (0, 0, 0, 1))$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
listarray (ep1LH);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
ep1LH[0];
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Define the function Ep1( lam, nu ) where lam is short for "lambda",
which takes on values 1, 0, -1, and describes the value of Sz of the W(-),
and the Lorentz index nu takes on values [0, 1, 2, 3].
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Ep1(lam, nu) :=  if lam = 1 then ep1RH[nu] 
                else if lam = -1 then ep1LH[nu]
                else if lam = 0 then ep1L[nu]
              else (print ("error"), done) $
/* [wxMaxima: input   end   ] */


/* [wxMaxima: section start ]
W^{-} Decay Helicity Amplitudes
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
See our expression for the decay amplitude Mfi at the top of this worksheet.

The vertex factor includes the so-called "chiral projection operator"
S(-1) = (1/2) ( 1 - γ^5 ), which has the explicit matrix representation  P(-1)
in our package.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
I4;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
(I4 - Gam[5])/2;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
P(-1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
P(-1) is "idempotent":  P(-1) . P(-1) = P(-1)
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
P(-1) . P(-1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
The explicit matrix P(s) is represented by  S(s) when constructing a symbolic trace or contraction.

To evaluate the decay amplitudes, we first need to construct the appropriate Dirac spinors.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
GENERAL COMMENTS: DIRAC SPINORS FOR SPIN 1/2 PARTICLES
                           AND ANTI-PARTICLES

  If upa is a Dirac spinor (column vector) corresponding to a
    lepton with 4-momentum pa, then  sbar(upa) is the barred Dirac 
    spinor (row vector) upa^{+} γ^0,   where (something)^{+} indicates 
     Hermitian conjugate.

   If vpa is a Dirac spinor (column vector) corresponding to an
     anti-lepton with 4-momentum pa, then sbar(vpa) is the barred 
     Dirac spinor (row vector) vpa^{+} γ^0,  where (something)^{+} 
     indicates the Hermitian conjugate.
    
   For leptons and anti-leptons with mass m,  our *normalization* is 
   upa_bar . upa = 2m and  vpa_bar . vpa = - 2m . These right hand sides
   are then equal to zero for neutrinos and anti-neutrinos (m = 0).
    
    UU(E,p,θ,φ,σ) is the Dirac spinor (4-element column vector) corresponding to a lepton with 
      relativistic energy  E = sqrt(p^2 + m^2) (if the lepton mass is m),  3-momentum magnitude p, 
      with direction of    3-momentum vector defined by the spherical polar angles θ and φ, and 
      helicity σ equal to   +1 for R  and  -1 for L.

   VV(E,p,θ,φ,σ) is the Dirac spinor (4-element column vector) corresponding to an anti-lepton 
     with  relativistic energy  E = sqrt(p^2 + m^2) (if the anti-lepton mass is m), 3-momentum magnitude p, 
     with direction of   3-momentum vector defined by the spherical polar angles θ and φ, and 
     helicity σ  equal to    +1 for R and  -1 for L. 
    
    The parameters E, p, θ, φ for each lepton or anti-lepton should agree with those specified in 
     the comp_def statement above.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
As a reminder, we repeat:

From Thomson (15.2), for given W^{-} polarization λ, electron helicity σ3, and 
electron anti-neutrino helicity σ4, the invariant amplitude is:

Mfi = (gw/sqrt(2)) ε(p1,λ)_{μ} ubar (p3, σ3) γ^{μ} (1/2) ( 1 - γ^5 ) v (p4, σ4)

The electron neutrino is left-handed and the electron anti-neutrino is
right handed (ignoring any possible non-zero mass for the neutrino),
so σ4 = +1.

We can write the amplitude 

Mfi = (gw/sqrt(2) * Mr 

in terms of the "reduced amplitude" Mr:
 
Mr =  ε(p1,λ)_{μ}   j^{μ}, 

where the leptonic weak charged current 4-vector is 

j^{μ} =  ubar (p3, σ3)  γ^{μ} S (-1) v (p4, σ4 = +1),    and where 

S(σ) =  (1/2) ( 1 + σ  γ^5 )  is the "chiral projection operator",  represented
by P(σ) with explicit matrix expressions. (See Thomson p. 141)
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
For the right-handed electron anti-neutrino with zero mass,
the spherical polar coordinate angles are
th4 = %pi - th, and phi4 = %pi when defining the spinor.

We could call this vp4_RH (for right-handed) since we are taking
σ4 = +1 for the (right-handed) anti-neutrino in the final state, but for
simplicity of notation we just call it vp4.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
vp4 : VV (E,E, %pi - th, %pi, 1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
The weak leptonic current 4-vector is proportional to S (-1) v (p4, σ4 = +1) and
this factor reduces to  v (p4, σ4 = +1)
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
P(-1) . vp4;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
so we can write the weak leptonic current 4-vector in the
simplified form:

  j^{μ} =  ubar (p3, σ3)  γ^{μ}  v (p4, σ4 = +1).

Because of the S(-1) operator, a hypothetical *left-handed* anti-neutrino gives
 zero contribution to the current 4-vector and thus zero contribution to the
decay process.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
vp4_LH : VV (E,E, %pi - th, %pi, -1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
P(-1) . vp4_LH;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
For use later when we use trace methods, let's establish the identity
(in matrix form here)

vp4 . sbar (vp4)  =   P(-1) . sL(p4)

which is a special case of the relation for any massless antileption:

v(p,σ) . sbar ( v(p,σ) ) = P(- σ) . sL (p)
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
is (equal ( vp4 . sbar (vp4), P(-1) . sL(p4) ) ) ;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
A : vp4 . sbar(vp4);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
B : P(-1) . sL(p4);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
A : fr_ao2 (A,th);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
is (equal (A, B) );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
 We write Mfi =  (gw/sqrt(2))  Mr, where the "reduced amplitude"
  can thus be written in the simplified form

Mr =  ε(p1,λ)_{μ} ubar (p3, σ3) γ^{μ}  v (p4, 1)

and we define a function Mr (Wlam, sv3) where Wlam stands for
the W-boson polarization parameter  λ, and sv3 stands for the electron helicity σ3.

vp4, defined above,  is a globally known matrix.

Wlam takes on the three values 1, 0, -1
sv3  takes on the two values 1, -1.

Here is our first version for the reduced amplitude Mr
function definition. We use mcon (expr, index) for the
needed contraction on mu.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Mr (Wlam, sv3) := 
  block ([up3b],
  up3b : sbar ( UU(E,E, th, 0, sv3)),
  mcon ( Ep1(Wlam, mu) * (up3b . Gam[mu] .  vp4),  mu) )$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
Mr (1,1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
Mr (1, -1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Our package function fr_ao2 (expr, angle) converts expr from a function of angle/2 to
a function of angle.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
factor (  fr_ao2 (%,th) );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Recall that M is the mass of the W-boson.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
  %, E = M/2;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Quoting Thomson, p. 409, "In the ultra-relativistic limit, where the helicity states are the same as the chiral
states, only left-handed helicity particle states and right-handed helicity antiparticle states contribute to
the weak interaction."

Hence we got zero amplitude for the case σ3 = +1 describing a right-handed electron helicity state ( h = 1/2).
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Now let's try to improve this function by incorporating
the desired simplifications:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Mr (Wlam, sv3) := 
  block ([up3b],
  up3b : sbar ( UU(E,E, th, 0, sv3)),
  mcon ( Ep1(Wlam, mu) * (up3b . Gam[mu] . vp4),  mu),
   factor (fr_ao2 (%%, th)),
   subst(E = M/2, %%))$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
We can make a table of reduced amplitudes using do loops and print.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
sL : [1, -1]$
wL : [1, 0, -1]$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
block ([ l1, s3],
    print ("  "),
    print ("  l1    s3        Mr "),
    for l1 in wL do
        for s3 in sL do  (
            print ("  "),
            print ("  ",l1,"   ",s3, "     ", Mr (l1, s3)) ),
    print ("  ") )$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: section start ]
Unpolarized Decay Rate from Helicity Amplitudes
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
 Use these amplitudes to calculate the unpolarized decay rate
of the W boson in its rest frame.

Recall that Mfi =  (gw/sqrt(2))  Mr, and we define  (averaging over the
three possible values of the initial helicity of the W boson)

mssq  =  Σ(sW)  Σ(se) |Mr (sW, se) |^2

MfiSQ =  (1/3)  Σ(sW)  Σ(se) |Mfi (sW, se) |^2  = (gw^2 / 6) mssq

   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
mssq : 0$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
block([l1, s3, temp],
    print ("  "),
    print ("  l1    s3        Mr "),
    for l1 in wL do
        for s3 in sL do (
                     temp :  Mr (l1, s3),                     
                     mssq : mssq + Avsq(temp),
                     print ("  "),
                     print ("  ",l1,"   ",s3, "     ", temp ) ),
    print ("   "))$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
mssq;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
mssq : trigsimp(mssq);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
MfiSQ : mssq *gw^2/6;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
This expression for the "spin-averaged matrix element squared" 
MfiSQ = < |Mfi|^2 > agrees with Thomson, Eq (15.6), p. 411.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Use this result with the two-body decay rate written down in the Sec. 2.6 of 
mbe12dirac3.pdf. 

Since we neglect the mass of the electron, p* in that expression is E = M/2.

MfiSQ is independent of angle, so the integral over the solid angle gives a factor
of  4 π steradians.  We then get for the decay rate in natural units

   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Gamma : gw^2 * M / (48 * %pi);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
 which agrees with Thomson (15.7).
   [wxMaxima: comment end   ] */


/* [wxMaxima: section start ]
Unpolarized Decay Rate Calculated Using Symbolic Trace Function tr
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
Recall that Mfi =  (gw/sqrt(2))  Mr, 

Mr =  ε(p1,λ)_{μ} ubar (p3, σ3) γ^{μ}  v (p4, σ4 = 1)

and we define  (averaging over the
three possible values of the initial helicity of the W boson)

MfiSQ =  (1/3)  Σ(sW)  Σ(se) |Mfi (sW, se) |^2  = (gw^2 / 6) mssq

mssq  =  Σ(sW)  Σ(se) |Mr (sW, se) |^2

In computing mssq we use the completeness relation for the W-boson
polarization 4-vectors, derived in massive-spin1-pol.wxm:

Σ(λ) ε^{μ} (p1, λ)   conjugate ( ε^{ν} (p1, λ) ) = - g^{μ, ν}  + (1/M^2)  p1^{μ}  p1^{ν} 

We also use:

Σ(σ3)   ( u3bar Γ u3 ) = trace (sL(p3) Γ )  =  trace ( Γ sL(p3)  )

The anti-neutrino is right handed (and massless),  there is no
sum over the anti-neutrino helicity, and  we can replace
v4  v4bar = P(-1) . sL(p4) (using the matrix representation).
(We showed this property above )

mssq then  amounts to the sum of two terms, which we write
mssq = B_term + A_term.

Since tr expects symbolic arguments P(-1) is represented by the symbol S(-1).
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
A : tr (mu,S(-1),p4,nu,p3);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Since S(-1) involves γ^5, the symbolic Eps factor appears from taking the trace, and
we will only be able to proceed with a final needed contraction on μ and ν if we
convert this to the non-covariant expression using our package function noncov.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
A_nc : noncov(A);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
ignoring the prefactor (1/M^2), 
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
A_term : mcon (p1[mu]*p1[nu]*A_nc, mu,nu);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
so the A_term doesn't contribute.  The B_term involves a contraction of the metric
tensor with the quantity A_nc; we must use mcon(expr, mu,nu) with explicit non-covariant
expressions.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
B_term : -mcon (gmet[mu,nu]*A_nc, mu,nu);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
mssq : trigsimp ( B_term), E = M/2;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
MfiSQ_tr : (gw^2/6) * mssq;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
is (equal (MfiSQ_tr, MfiSQ));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
so we get the same final result for MfiSQ = <|Mfi|^2>.
   [wxMaxima: comment end   ] */


/* [wxMaxima: section start ]
Unpolarized Decay Rate Using Noncov-Symbolic Trace Function nc_tr
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
Using nc_tr instead of just tr causes noncov to be immediately called
on the result returned by the symbolic trace function tr, which automatically
contracts on repeated indices inside the tr(...) function.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
A : nc_tr (mu,S(-1),p4,nu,p3);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
A_term : mcon (p1[mu]*p1[nu]*A, mu,nu)/ (M^2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
B_term : -mcon (gmet[mu,nu]*A, mu,nu);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
mssq : trigsimp ( B_term), E = M/2;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: section start ]
Unpolarized Decay Rate Using the Matrix Trace Function m_tr
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
Both nc_tr and m_tr have the advantage of using the same formal arguments
that are used with the symbolic trace function tr. m_tr (expr) converts all arguments
to explicit 4 x 4 matrices.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
A : m_tr (mu,S(-1),p4,nu,p3);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
A_term : mcon (p1[mu]*p1[nu]*A, mu,nu)/ (M^2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
trigsimp(A_term);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
B_term : -mcon (gmet[mu,nu]*A, mu,nu);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
mssq : trigsimp ( B_term), E = M/2;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: section start ]
Unpolarized Decay Rate Calculated Using mat_trace and Explicit Matrices
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
Also note all the arguments to mat_trace are explicit matrices, with matrix P(-1) replacing the symbol S(-1), for
example.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
A : mat_trace (Gam[mu] . P(-1) . sL(p4) . Gam[nu] . sL(p3)) ;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
A_term : mcon (p1[mu]*p1[nu]*A, mu,nu)/ (M^2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
trigsimp (A_term);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
B_term : -mcon (gmet[mu,nu]*A, mu,nu);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
mssq : trigsimp ( B_term), E = M/2;
/* [wxMaxima: input   end   ] */



/* Old versions of Maxima abort on loading files that end in a comment. */
"Created with wxMaxima 19.05.7"$