/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 17.10.0 ] */
/* [wxMaxima: title   start ]
ee-mumu-HE.wxm   High Energy Limit of

    e(-,p1,σ1) + e(+,p2,σ2) --> μ(-,p3,σ3) + μ(+,p4,σ4)

   [wxMaxima: title   end   ] */


/* [wxMaxima: comment start ]
ee-mumu-HE.wxm,  Oct. 18, 2018
Edwin (Ted) Woollett, Maxima by Example,  
ch. 12,   Dirac Algebra and Quantum Electrodynamics, ver. 3, "Dirac3"
 http://web.csulb.edu/~woollett/ 

"High Energy"  means center of momentum frame particle
energy E  >>  M  >>  m so we can "neglect" the particle masses.

 We use Heaviside-Lorentz electromagnetic units.
In general, we follow the conventions in Peskin & Schroeder: Quantum Field Theory.
   [wxMaxima: comment end   ] */


/* [wxMaxima: section start ]
HE UNPOLARIZED CROSS SECTION USING TRACES AND CONTRACTIONS
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
First load the dirac3 package by loading dirac3.mac:
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
Note from the last two printouts, that invar_flag = true and stu_flag = false (our defaults).

The fact that invar_flag = true means that the symbolic trace function tr
will automatically make use of the list invarR (produced by the function set_invarR) which defines
the replacements for D(pa,pb) in terms of the Mandelstam variables s, t, u.
(After the package loads, invarR has the default value [].)

The fact that stu_flag = false means that the function noncov will not
automatically replace the Mandelstam variables s, t, u, by s_th, t_th, and u_th.
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
UNPOLARIZED CROSS SECTION Using Mandelstam Variables s, t, and u
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
This first section works out the unpolarized differential
  cross section in the high energy (HE) limit, first in an arbitrary frame in terms
  of the Mandelstam variables s, t, and u, and then  in the CMS frame, using symbolic methods.

These symbolic methods use our dirac3 package functions set_invarR, tr, Con, pullfac,
 comp_def, VP, sub_stu, to_ao2, fr_ao2

The amplitude for specific helicities σ1, σ2, σ3, σ4, used in Sec. 2 which employs
explicit Dirac spinors, is

M(σ1, σ2, σ3, σ4) = -e^2 Mr = (v2bar γ^μ u1) ( u3bar γ_μ v4)/ s,
summed over the Lorentz index μ , where s = (p1 + p2)^2 = (p3 + p4)^2.
We call Mr the "reduced" amplitude.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
To introduce Mandelstam variables s, t, and u into the calculation,
first POPULATE THE LIST invarR OF symbolic D(pa,pb) 4-VEC DOT PRODUCT VALUES,
  using  set_invarR( D(p1,p1) = r11, D(p1,p2) = r12, ...)
using 
    p1 + p2 = p3 + p4, 
    s = (p1+p2)^2 = (p3+p4)^2,
     t = (p1-p3)^2 = (p2-p4)^2, and  
     u = (p1-p4)^2 = (p2-p3)^2 

These definitions are independent of coordinate frame.
For later work, m is the electron (positron) mass, and M is the mass of the muon (antimuon).

D(pa,pb)  is a symbolic representation of the 4-vector dot product (scalar product)
of 4-vectors pa and pb.

In this high-energy limit, E >>  M >>  m so effectively  M --> 0, m--> 0 and
 D(p1,p1) approx 0  instead of m^2, etc.

The dimensions of s, t and u are all the same
[s]  = [t]  = [u]  = E^2
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
set_invarR (D(p1,p1) = 0,
         D(p1,p2) = s/2,
         D(p1,p3) =  - t/2,
         D(p1,p4) =  - u/2,
         D(p2,p2) = 0,
         D(p2,p3) =  - u/2,
         D(p2,p4) =  - t/2,
         D(p3,p3) = 0,
         D(p3,p4) = s/2,
         D(p4,p4) = 0)$
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
There is only one diagram in lowest order, with a virtual photon in the
intermediate state.

omit the factor -e^2 from Mfi for now, and let MfiSQ be the sum over all four
helicities of the absolute value squared of the (reduced) amplitude,
and write MfiSQ = Mn/s^2 where s = (p1 + p2)^2 = E(cm)^2

   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Mn : factor (Con (tr(p2,mu,p1,nu) * tr(p3,mu,p4,nu), mu,nu) );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
MfiSQ : Mn/s^2;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
We have absorbed e^4 into a factor we call A, with e^2 = 4πα.

In general 2 --> 2 scattering,  A = (α^2/(4 s)) (pf/pi)  but in 
the high energy limit we neglect masses so  |p-vec|  
approx = particle energy, same for each lepton here,
so pf = pi, and A =  α^2/(4 s)
  
 Averaging over initial spins means we need to divide A by 4
  to get the unpolarized differential cross section in an arbitrary
  coordinate frame but in high energy limit.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
A : alpha^2/(4*s);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
dsigdo_unpol_HE : (A/4)*MfiSQ;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: subsect start ]
CONVERSION TO EXPLICIT FUNCTION OF SCATTERING ANGLE θ 
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
   CONVERSION TO EXPLICIT FUNCTION OF SCATTERING ANGLE θ 
(in center of momentum frame); th stands for the scattering angle.
Replace s, t, and u with explicit functions s_th, t_th, u_th of th,
using our package function sub_stu(expr).
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
The syntax comp_def ( pa(Ea, pax, pay, paz) ) assigns components to the array pa,
in which pa[0] = Ea, pa[1] = pax, etc.
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsubsect start ]
Invocation of comp_def(...) to define CMS 4-momentum components
   [wxMaxima: subsubsect end   ] */


/* [wxMaxima: comment start ]
The syntax comp_def ( pa(Ea, pax, pay, paz) ) assigns components to the array pa,
in which pa[0] = Ea, pa[1] = pax, etc.

momRS is a list maintained by the dirac package of 4-vectors 
which have been assigned components using comp_def.

Each such 4-vector corresponds to a zero-based array.
   [wxMaxima: comment end   ] */


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


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


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


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


/* [wxMaxima: subsubsect start ]
Use of 4-vector product function VP(...) to define s_th, t_th, and u_th
   [wxMaxima: subsubsect end   ] */


/* [wxMaxima: comment start ]
Let s_th, t_th, u_th be frame dependent values of the Mandelstam variables
s, t, and u. (Here CMS).
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
s_th : VP (p1+p2,p1+p2);
  
  t_th : factor (VP (p1-p3,p1-p3));
  
  u_th : factor (VP (p1-p4,p1-p4));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: subsubsect start ]
Use sub_stu(expr) to replace s by s_th, t by t_th, and u by u_th
   [wxMaxima: subsubsect end   ] */


/* [wxMaxima: input   start ] */
MfiSQ_th  :  factor( sub_stu(MfiSQ)  );
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
dsigdo : (A/4)*MfiSQ_th;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Since s = E(cm)^2 = (2*E)^2 = 4*E^2, this agrees with Peskin/Schroeder  Eq (5.14), page 137.

To get a Greek θ in the result, we need to replace th with theta:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
dsigdo, th = theta;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: subsect start ]
UNPOLARIZED CROSS SECTION  WITHOUT Using Mandelstam Variables
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
reset the list invarR of invariant replacement rules to avoid the Mandelstam variables s, t, and u,
using set_invarR.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
set_invarR ( D(p1,p1) = 0,D(p2,p2) = 0,D(p3,p3) = 0,D(p4,p4) = 0 )$
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
We calculate the traces symbolicaly using tr(...) and then calculate the contraction
of the Lorentz indices mu and nu using Con or scon.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Mn : factor (Con (tr(p2,mu,p1,nu) * tr(p3,mu,p4,nu), mu,nu) );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: subsubsect start ]
Use noncov(expr) To Make Use of Frame Dependent Quantities Defined by comp_def(...)
   [wxMaxima: subsubsect end   ] */


/* [wxMaxima: comment start ]
We now use noncov (expr) to evaluate this expression in the center of momentum
frame. This works because we have already used comp_def (...) to define the components
of the relevant 4-momenta above.
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: input   start ] */
MfiSQ_th : %/s_th^2;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
dsigdo : (A/4)*MfiSQ_th;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
which agrees with our results using the Mandelstam variable route.

To summarize, this symbolic method needs either an empty list invarR,
or at a list invarR which has no Mandelstam variable replacement
information regarding the dot products of the 4-momenta involved.

But this symbolic method needs the use of comp_def to define 4-momenta
components, and also,  the evaluation of the denominator(s) of the amplitude
contributions, in this example s^2  which requires s_th : VP(p1 + p2, p1 + p2).
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Recall from the top of this worksheet that since stu_flag = false, noncov does NOT
replace s by s_th, t by t_th, u by u_th.  You can still use sub_stu(expr) for this job
(as long as you have defined s_th, t_th, and u_th using VP(pa,pb) ).
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
sub_stu(s);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
sub_stu(Mn/s^2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: subsect start ]
UNPOLARIZED CROSS SECTION, WITHOUT  Mandelstam Variables, Using nc_tr(...)
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
nc_tr  is a special function which is not just tr followed by noncov, but
  rather as TR1 generates the trace of each term of the initial argument
  expansion, the function noncov is immediately applied to the trace result,
  and the output is accumulated in a sum which is eventually returned by
  nc_tr. For expressions with many terms, this is a more efficient method
of symbolic trace evaluation for the whole complicated expression.

We first make sure invarR is either an empty list or a list of replacement rules
which don't involve the Mandelstam variables s, t, and u.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: subsubsect start ]
For Contraction, Use  mcon (expr, index1,index2,...) with traces produced using nc_tr
   [wxMaxima: subsubsect end   ] */


/* [wxMaxima: comment start ]
The most efficient (fast) contraction method with traces produced by using nc_tr or m_tr
is mcon (expr,  contraction indices). Speed criteria become more important as you
deal with an expression with many terms.

   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Mn : mcon (nc_tr(p2,mu,p1,nu)*nc_tr(p3,mu,p4,nu),mu,nu);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
We find that trigsimp(expr) does not convert to a form
involving cos(th)^2  for this example
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
so we resort to using the dirac package trigsimp alternative:   ts(expr, var):
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
ts(Mn, th);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
we now get MfiSQ (to compare with our previous work)
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
%/s_th^2;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
and dsigdo:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
(A/4)*%;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
%, th = theta;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: subsect start ]
Explicit Matrix Traces and Contractions: the m_tr(...) Method
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
m_tr uses the same syntax as tr, but translates everything into explicit matrices
and doesn't call the symbolic tr function, so the nature of the invarR list does
not matter.  

The explicit matrix method using m_tr for traces and mcon for contraction
is the fastest general method.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
m_tr (p3,mu,p1,nu) is converted into the explicit matrix expression:

mat_trace (sL(p3) . Gam[mu] . sL(p1) . Gam[nu] )

where mat_trace is defined in the linearalgebra package and is
automatically loaded when we first try to use mat_trace.

In this example, p1[0] = E, p1[1] = 0, p1[2] = 0, p1[3] = E,  and
p1 is a "zero based array" holding the components of the
4-momentum vector p1. The components were defined
above using the comp_def function.

sL(p1) is the Feynman slash  matrix corresponding to 
the 4-momentum p1.

Gam[mu] is the explicit matrix representation of γ^μ, for example.
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsubsect start ]
Use mcon (expr, index1,index2,..) for traces produced using m_tr or mat_trace
   [wxMaxima: subsubsect end   ] */


/* [wxMaxima: input   start ] */
Mn : mcon (m_tr(p2,mu,p1,nu) * m_tr(p3,mu,p4,nu), mu,nu);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
%/s_th^2;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
(A/4)*%;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
%, th = theta;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: subsect start ]
The Explicit Matrix Method Using mat_trace plus mcon Directly
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
Instead of using the convenient m_tr syntax (which follows the symbolic tr and nc_tr syntax), 
we can work directly with explicit matrices to get the same final results.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Mn : mcon (mat_trace (sL(p2) . Gam[mu] . sL(p1) . Gam[nu]) * mat_trace (sL(p3) . Gam[mu] . sL(p4) . Gam[nu]), mu,nu);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
factor(Mn) / s_th^2;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
(A/4) * %;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
%,  th = theta;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: section start ]
HE POLARIZED AMPLITUDES USING EXPLICIT DIRAC SPINORS
e(-,p1, σ1) +  e(+,p2, σ2)  -->  μ(-,p3, σ3) + μ(+,p4, σ4)
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
Using explicit Dirac matrices, in this section we use the package
    functions UU, VV, sbar, Gam, mcon, sub_stu, Avsq, ts, to_ao2,  fr_ao2, Esq_to_psq, Ep_to_m
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Since this worksheet deals with the high energy (HE) limit, 
we should make use of what Halzen & Martin call "helicity
conservation at high energies,"  (see p. 126).
They show that for lepton pair annihilation at high energies,
the incoming electron and positron have opposite helicities,
and likewise for the outgoing muon and antimuon.

A lengthier discussion can be found in Mark Thomson, Modern
Particle Physics, sections 6.2 through 6.4.
 
GENERAL COMMENTS: DIRAC SPINORS FOR 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.
    
 Our normalization is  upa_bar . upa = 2m and  vpa_bar . vpa = - 2m
    
    UU(E,p,θ,φ,σ) is the Dirac spinor (4-element column vector) corresponding to a lepton with 
      relativistic energy  E = sqrt(p^2 + m^2), 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), 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.
    
    In the center of momentum frame, the incident electron e(-) and positron e(+)have equal and opposite
    3-momentum vectors, and same for the outgoing muon μ(-) and anti-muon μ(+).
    
    These parameters should be compatible with those specified in the comp_def statement above.
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
e(-,p1, R) +  e(+,p2, R)  -->  μ(-,p3, R) + μ(+,p4, R)
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
We start with RR --> RR case, just to confirm
that we get an amplitude which is zero.

case  RR --> RR, ie, (+1,+1) --> (+1,+1) polarized amplitude.
  e(-,p1,R) + e(+,p2,R) --->  μ(-, p3, R)  + μ(+, p4, R)
  [ note:  R => σ = +1,  L =>  σ = -1 ]
      
  First define the needed Dirac spinors and barred spinors.

   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
up1:UU(E,E,0,0,1);
up3b:sbar(UU(E,E,th,0,1));
vp2b:sbar(VV(E,E,%pi,0,1));
vp4:VV(E,E,%pi-th,%pi,1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
check normalization of vp4; in this high energy limit, vp4b . vp4 = - 2 M  --> 0
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
grind(%);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
vp4b . vp4;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
For helicity quantum numbers RR -> RR,  the amplitude Mfi
 (pulling out -e^2) has the form Mfi = Mn/s 
 where  s = (p1 + p2)^2.

For contraction over Lorentz indices involving explict Dirac matrices,
use mcon (expr, index,...) syntax.  

Gam[mu] is an explicit matrix representing γ^{μ}
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Mn : mcon ( (vp2b . Gam[mu] . up1)*(up3b . Gam[mu] . vp4), mu);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
So in the high energy limit, the amplitude for RR --> RR is zero, and dsigdo_RR_RR = 0.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
"Helicity conservation" for massless particles requires the incoming e(-) and e(+) to have opposite
helicities; the same is true for the outgoing μ(-) and μ(+) in the final state.
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
e(-,p1, R) +  e(+,p2, L)  -->  μ(-,p3, R) + μ(+,p4, L)
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
So look at the (reduced) amplitude for e(-,p1,R) + e(+,p2,L) --> μ(-,p3,R) + μ(+,p4,L)
This is case  RL --> RL
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
up1:UU(E,E,0,0,1);
up3b:sbar(UU(E,E,th,0,1));
vp2b:sbar(VV(E,E,%pi,0, -1));
vp4:VV(E,E,%pi-th,%pi, -1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: subsubsect start ]
Use mcon (expr, index1, index2,...) for contractions of explicit matrix expressions
   [wxMaxima: subsubsect end   ] */


/* [wxMaxima: input   start ] */
Mn : mcon ( (vp2b . Gam[mu] . up1)*(up3b . Gam[mu] . vp4), mu);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
Mfi : Mn / s_th;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
convert f(th/2) to g(th) using package function fr_ao2(expr, th)
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
so the (reduced) amplitude for the case RL -->  RL is ( 1 + cos(θ) ), and reminding ourselves
of the definition of A,
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
dsigdo_RL_RL : A * Mfi^2;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
%, th = theta;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Since s = E(cm)^2 = (2 E)^2 = 4 E^2, this agrees with Peskin/Schroeder p. 143, Eq (5.22).
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
General Helicity Cases using explicit Dirac spinors
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
In the following function definition,
s1v is the helicity of the particle with 4-momentum p1 (the incident electron), etc.

Define the function he_me(s1v,s2v,s3v,s4v) which returns the high energy limit of
the (reduced) amplitude for the case that the helicities have the specified values.

note: s1v = (1 or -1) for (R, or L) , ie, right-handed or left-handed.
The quantitiy s_th  was defined above and  globally known.
E and th are undefined parameters.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
he_me (s1v,s2v,s3v,s4v) :=
 block ([up1,up3b,vp4,vp2b],
      up1:UU(E,E,0,0,s1v),
      vp2b:sbar(VV(E,E,%pi,0, s2v)),
      up3b:sbar(UU(E,E,th,0,s3v)),
      vp4:VV(E,E,%pi-th,%pi, s4v),         
      mcon ( (vp2b . Gam[mu] . up1)*(up3b . Gam[mu] . vp4), mu),      
      trigsimp (%%) /s_th,
      fr_ao2 (%%, th))$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Let's check the case RR --> RR we worked out above.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
which agrees.
Next check case RL --> RL
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
which also agrees.
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
Polarized amplitudes --> unpolarized cross section
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
AN ALTERNATIVE PATH TO THE UNPOLARIZED CROSS SECTION
 IS TO CALCULATE THE INDIVIDUAL AMPLITUDES FOR SPECIFIC 
POLARIZATIONS  OF THE INITIAL AND FINAL PARTICLES (SEE 
SEC. 3 "SPIN SUMS" IN CH. 12) AND THEN FIND THE 
UNPOLARIZED CROSS SECTION BY  COMPUTING THE SUM OF 
THE ABSOLUTE VALUE SQUARED OF EACH OF THE POLARIZED 
AMPLITUDES, WHICH WE NOW DO.

 The function Avsq(expr) is defined in dgmatrix3.mac 
and computes the absolute value squared.
  (We could also use here abs(expr)^2 )

mssq stands here for < |M|^2 > = sum over spins of  |M(σ1,σ2,σ3,σ4)|^2

Using the function he_me defined in terms of the amplitude expressed
in terms of explicit Dirac spinors to calculate mssq takes about 0.6 - 0.8 sec.
(although this includes the time spent calling print)
on the author's machine, and is the fastest route to the cross section.

Of course this time does not include the time spent in setting up the
function he_me.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
block([s1,s2,s3,s4,temp],    
    print ("  "),
    print (" s1 s2  s3  s4     amplitude     "), 
    for s1 in sL do
        for s2 in sL do
            for s3 in sL do
                for s4 in sL do (
                     temp :  he_me (s1,s2,s3,s4),                     
                     mssq : mssq + Avsq(temp),
                     print ("  "),
                      print (s1,s2,s3,s4," ",temp ) ),
    print ("   "))$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
time(%);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
we computed MfiSQ_th using nc_tr plus mcon in Sec. 1.3.1.
This should be the same as mssq.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
which shows equality with our symbolic trace, contraction route to the unpolarized cross section.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
dsigdo : (A/4) * factor(mssq);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
%, th = theta;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: section start ]
Trace Methods and Helicity Projection Matrices for Square of Polarized Amplitudes
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
We can check the (square) of the polarized amplitude calculated above 
  (using Dirac spinors for given helicity assignments) by using either 
 symbolic nc_tr(..) plus mcon method, or the explicit matrix method:  m_tr plus mcon .

You should NOT use tr(...) plus Con (or scon) for this calculation.

The matrix method using m_tr syntax is designed to have
  arguments which look the same as the purely symbolic case nc_tr(...)

In both the nc_tr and m_tr methods, one must use a spin projection operator
 symbolized by S(σ, Sp ) (for the massive lepton or massive anti-lepton case) 
 or S(σ) for the massless lepton case and  S(-σ)  for the massless anti-lepton case. 

We are ignoring the mass m of the electron and anti-electron and the mass M of the muon and
anti-muon, so we use either S(σ) or S(-σ) as the spin projection operators.

For example of massive particles, see ee-mumu-AE.wxm

   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
Use of nc_tr plus mcon  for squared amplitudes and cross section using spin projection operators
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
general amplitude squared for specified helicities using nc_tr method
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
meSQ(s1v,s2v,s3v,s4v) := factor(
      mcon ( nc_tr(S(-s2v),p2,mu,S(s1v),p1,nu) * nc_tr(S(s3v),p3,mu,S(-s4v),p4,nu), mu,nu))/s_th^2$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
case RR --> RR should be zero
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
case RL -->  R L
   [wxMaxima: comment end   ] */


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


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


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


/* [wxMaxima: comment start ]
which agrees with the explicit Dirac spinor calculation of the square of the
polarized amplitude A_RL_RL

We can use these squared polarized amplitudes, produced using nc_tr plus mcon
to find the unpolarized differential scattering cross section  dsigdo.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
block([s1,s2,s3,s4,temp],    
    print ("  "),
    print("  using nc_tr(...) method "),
    print (" s1 s2  s3  s4     |amplitude|^2     "), 
    for s1 in sL do
        for s2 in sL do
            for s3 in sL do
                for s4 in sL do (
                    temp : meSQ(s1,s2,s3,s4),
                     mssq : mssq + temp,
                     print ("  "),
                      print (s1,s2,s3,s4," ",temp ) ),
    print ("   "))$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
time(%);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
ts (mssq,th), factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
dsigdo : (A/4) * %;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
which agrees with our previous results for the unpolarized differential scattering
cross section using the polarized amplitudes via explicit Dirac spinor matrices.

   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
Use of m_tr plus mcon  for squared amplitudes and cross section using spin projection operators
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
We define a function meSQ which uses m_tr and mcon and which returns the absolute value
squared of the (reduced) amplitude for specified helicities for this process.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
meSQ(s1v,s2v,s3v,s4v) := factor(
      mcon ( m_tr(S(-s2v),p2,mu,S(s1v),p1,nu) * m_tr(S(s3v),p3,mu,S(-s4v),p4,nu), mu,nu))/s_th^2$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
check RR --> RR case, which should be zero
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
check RL --> RL case
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
which agrees with the explicit Dirac spinor calculation of the square of the
polarized amplitude A_RL_RL

We can use these squared polarized amplitudes, produced using m_tr plus mcon
to find the unpolarized differential scattering cross section  dsigdo.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
block([s1,s2,s3,s4,temp],    
    print ("  "),
    print("  using m_tr(...) method "),
    print (" s1 s2  s3  s4     |amplitude|^2     "), 
    for s1 in sL do
        for s2 in sL do
            for s3 in sL do
                for s4 in sL do (
                    temp : meSQ(s1,s2,s3,s4),
                     mssq : mssq + temp,
                     print ("  "),
                      print (s1,s2,s3,s4," ",temp ) ),
    print ("   "))$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
time(%);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
ratsimp(mssq), factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
dsigdo : (A/4) * %;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: subsect start ]
Use of mat_trace plus mcon  for squared amplitudes using spin projection operators
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
The explicit m_tr method converts the symbols in the call to explixit matrices,
with the contraction carried out by mat_trace.

Here is an example for the RR --> RR case.

The matrix P(σ)  replaces the symbol S(σ).
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: subsect start ]
Explicit Spin Projection Matrix P(sv) for the High Energy Limit Case
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
NOTE: THIS IS A PEDAGOGICAL DIGRESSION
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Explicit Spin Projection Matrix P(sv) for the High Energy limit case has the properties:

let p be 4-momentum, pa be the magnitude of the corresponding 3-momentum
E  =  sqrt (m^2 + pa^2) ---> pa  for pa >> m  ( ==>  E >>  m )

massless FERMIONS:

u(p,R) . sbar (u(p,R))  =  P(1) . sL(p)
u(p,L) . sbar (u(p,L))  =  P(-1) . sL(p)

where u(p, s) : UU(E,pa,th,phi,s), 

u(p,R) = u(p,1)  = UU(E,pa,th,phi,1)
 u(p,L) = u(p,-1)  = UU(E,pa,th,phi,-1)
-----------------------------------------------------

Massless ANTI-FERMIONS

v(p,R) . sbar (v(p,R))  =  P(-1) . sL(p)
v(p,L) . sbar (v(p,L))  =  P(1) . sL(p)

where v(p, s) : VV(E,pa,th,phi,s), 

v(p,R) = v(p,1) =  VV(E,pa,th,phi,1),
v(p,L) = v(p,-1) =  VV(E,pa,th,phi,-1),

recall that the components of 4-momentum p are
such that listarray(p) returns the list [E, px, py, pz], where
px = pa sin(th) cos(phi), py = pa sin(th) sin(phi),
and pz = pa cos(th), using spherical polar coordinates,
and then pa^2 = px^2 + py^2 + pz^2.

In the high energy limit E >> M >> m, pa = E.

   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
In general,  P(s) is the explicit matrix  (I4 + s*Gam[5]) / 2 

Gam[5] is the 4 by 4 matrix representing γ^5.

I4 is the 4 by 4 unit matrix.

 s is restricted to the values [1, -1] for our purposes.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
is(equal( P(s), (I4 + s*Gam[5])/2 ));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
The simplest 4-momentum vector we have defined (above, using comp_def)
is p1. Let's use this 4-momentum to illustrate the above assertions.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
here we let up1 correspond to a RIGHT handed helicity fermion
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
up1:UU(E,E,0,0,1);
up1b : sbar (up1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
is (equal (up1 . up1b, P(1) . sL(p1)));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
here we let up1 correspond to a LEFT handed helicity fermion
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
up1:UU(E,E,0,0,-1);
up1b : sbar (up1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
is (equal (up1 . up1b, P(-1) . sL(p1)));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
here we let vp1 correspond to a RIGHT handed helicity anti-fermion
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
vp1 : VV(E,E,0,0,1);
vp1b : sbar (vp1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
is (equal (vp1 . vp1b, P(-1) . sL(p1)));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
here we let vp1 correspond to a LEFT handed helicity anti-fermion
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
vp1 : VV(E,E,0,0,-1);
vp1b : sbar (vp1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
is (equal (vp1 . vp1b, P(1) . sL(p1)));
/* [wxMaxima: input   end   ] */



/* Maxima can't load/batch files which end with a comment! */
"Created with wxMaxima"$