# ee-mumu-AE.wxmx    Arbitrary Energy Case of

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

ee-mumu-AE.wxmx,   Oct. 19, 2018
Edwin (Ted) Woollett, Maxima by Example,  
ch. 12,   Dirac Algebra and Quantum Electrodynamics, ver. 3, "Dirac3"

http://web.csulb.edu/~woollett/ 
woollett@charter.net

"Arbitrary Energy"  means center of momentum frame particle
energy E  is not necessarily large compared to electron mass m.

 We use Heaviside-Lorentz electromagnetic units.
In general, we follow the conventions in Peskin & Schroeder: Quantum Field Theory.

## UNPOLARIZED CROSS SECTION USING TRACES AND CONTRACTIONS

First load the dirac3 package by loading dirac3.mac:

load(dirac3);

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.

### UNPOLARIZED CROSS SECTION Using Mandelstam Variables s, t, and u

This first sub-section works out the unpolarized differential
  cross section for arbitrary energy (AE), 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, ts, Esq_to_psq, psq_to_Esq.

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 
where 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.

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

The dimensions of s, t and u are all the same
[s]  = [t]  = [u]  = E^2

set_invarR (D(p1,p1) = m^2,
         D(p1,p2) =  s/2 - m^2,
         D(p1,p3) = (m^2 + M^2)/2 - t/2,
         D(p1,p4) = (m^2 + M^2)/2 - u/2,
         D(p2,p2) = m^2,
         D(p2,p3) = (m^2 + M^2)/2 - u/2,
         D(p2,p4) = (m^2 + M^2)/2 - t/2,
         D(p3,p3) = M^2,
         D(p3,p4) = s/2 - M^2,
         D(p4,p4) = M^2)$

invarR;

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

Mn : factor (Con (tr(p2 - m, mu, p1 + m, nu) * tr(p3 + M, mu, p4 - M, nu), mu,nu) );

MfiSQ : Mn/s^2;

We have absorbed e^4 into a prefactor A, with e^2 = 4πα.

For  2 --> 2 scattering,  A = (α^2/(4 s)) (pf/pi)  where
 s = (p1 + p2)^2 = E(cm)^2, pf is the magnitude of the final
muon's 3-momentum, which we will call k, and pi is the
magnitude of the initial electron's 3-momentum, which we
will call p.
  
 Averaging over initial spins of the electron and positron
means we need to divide A by 4   to get the unpolarized 
differential cross section in an arbitrary  coordinate frame 
but in the high energy limit.

So dsigdo = (A/4) * MfiSQ

A : alpha^2 * k/(4*s*p);

Then the unpolarized differential scattering cross section is

dsigdo : (A/4) * MfiSQ;

We can recover the high energy limit (m = 0, M = 0, k = p effectively)

dsigdo, m = 0, M = 0, k = p;

which agrees with result in ee-mumu-HE.wxmx.

### CONVERSION TO EXPLICIT FUNCTION OF SCATTERING ANGLE θ

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

### Invocation of comp_def(...) to define CMS 4-momentum components

We first need to use comp_def(...) to define the components of the 4-momenta
p1, p2, p3, p4.

The syntax comp_def ( p1(Ea, pax, pay, paz) ) assigns components to the array pa,
in which pa[0] = Ea, pa[1] = pax, etc. See ee-mumu-HE.wxmx for an example.

Initial electron and positron have p1vec = -p2vec, with |p1vec| = p = |p2vec|, and hence
equal particle energies E = sqrt(m^2 + p^2). The center of momentum frame
system energy E(cm) = 2 E, and conservation of system total energy means
that the muon and anti-muon energies in the final state add up to 2 E. Conservation of
system 3-momentum means p3vec = -p4vec, with |p3vec| = k = |p4vec|, and
equal energies E(mu) = sqrt(M^2 + k^2). Then 2 E(mu) = E(cm) = 2 E implies
that E(mu) = E.

We can solve for k and p in terms of E and the masses:

p = sqrt(E^2 - m^2),  k = sqrt(E^2 - M^2)

th is the angle the final state muon 3-momentum vector p3vec makes with the z axis
(the initial state electron 3-momentum vector p1vec is in the direction of the
positive z-axis.

assume ( E > 0, p > 0, k > 0,  th >= 0, th <= %pi )$
comp_def ( p1( E,0,0,p),
            p2( E,0,0,-p),
            p3 (E,k*sin(th),0,k*cos(th)),
            p4 (E,-k*sin(th),0,-k*cos(th)) )$

comp_def creates a separate (zero based)  array for each 4-vector.
Use listarray to see the components

listarray(p4);

p4[0];

### Use of 4-vector product function VP(...) to define s_th, t_th, and u_th

Let s_th, t_th, u_th be frame dependent values of the Mandelstam variables
s, t, and u. (Here CMS).

Use VP(pa,pb) to form the dot product of the pair of 4-vectors. (VP = "vector product").

s_th : VP (p1+p2,p1+p2);
  
  t_th : VP (p1-p3,p1-p3);
  
  u_th : VP (p1-p4,p1-p4);

### Use sub_stu(expr) to replace s by s_th, t by t_th, and u by u_th

MfiSQ_th : sub_stu(MfiSQ);

We use the package trigsimp alternative ts (expr, angle) to replace sin(th) in terms of cos(th):

MfiSQ_th : factor (ts (MfiSQ_th, th));

This is dimensionless, since E, k, p each have the dimension of mass in our units.

Look at the high energy limit:

MfiSQ_th, k = E, p = E, m = 0, M = 0;

factor(%);

which agrees with the high energy limit.

Returning to the arbitrary energy case:

A;

dsigdo : (A/4)*MfiSQ_th;

In this expression we can replace p and k in terms of E, M, and m if we wish.

 we first check on the high energy limit:

dsigdo, k = E, p = E, m = 0, M = 0;

factor(%);

which agrees with the high energy limit.

### UNPOLARIZED CROSS SECTION  WITHOUT Using Mandelstam Variables

Instead of using sub_stu, we use noncov in this sub-section.

reset the list invarR of invariant replacement rules to avoid the Mandelstam variables s, t, and u,
using set_invarR.

set_invarR ( D(p1,p1) = m^2,D(p2,p2) = m^2,D(p3,p3) = M^2, D(p4,p4) = M^2 )$

invarR;

We calculate the traces symbolicaly using tr(...) and then calculate the contraction
of the Lorentz indices mu and nu using Con or scon.

Mn_meth2 : factor (Con (tr(p2 - m, mu, p1 + m, nu) * tr(p3 + M, mu, p4 - M, nu), mu,nu) );

### Use noncov(expr) To Make Use of Frame Dependent Quantities Defined by comp_def(...)

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.

noncov (Mn_meth2);

Mn_meth2_nc : pullfac(%,32);

MfiSQ_th_meth2 : %/s_th^2;

to compare this with method 1, express k and p in terms of E, M, and m.

psq_to_Esq(MfiSQ_th_meth2, k,E,M), expand;

MfiSQ_th_meth2_E : psq_to_Esq(%,p,E,m), expand;

Remind ourselves what we got for MfiSQ_th in the Mandelstam variable route:

MfiSQ_th;

psq_to_Esq(MfiSQ_th,k,E,M), expand;

MfiSQ_th_E : psq_to_Esq (%,p,E,m), expand;

is (equal (MfiSQ_th_E, MfiSQ_th_meth2_E));

so  the expressions are shown to be equivalent when we take into account the
dependence of k and p on E.

Method 2 (non-Mandelstam route) produced a simpler expression in terms of
E, k, p, M, and m, which results in:

A;

dsigdo_meth2 : (A/4) * factor( MfiSQ_th_meth2);

check the high energy limit

dsigdo_meth2, k = E, p = E, M = 0, m = 0, factor;

which agrees.

To summarize, this non-Mandelstam 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 before using tr(...).

But this second 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).

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

### UNPOLARIZED CROSS SECTION, WITHOUT  Mandelstam Variables, Using nc_tr(...)

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.

invarR;

### For Contraction, Use  mcon (expr, index1,index2,...) with traces produced using nc_tr

The most efficient (fast) contraction method with traces produced by using nc_tr
is mcon (expr,  contraction indices). Speed criteria become more important as you
deal with an expression with many terms.

mcon (nc_tr(p2 - m,mu,p1 + m,nu)*nc_tr(p3 + M,mu,p4 - M,nu),mu,nu);

ts (%, th), expand;

Mn_meth3 : pullfac (%, 32);

MfiSQ_th_meth3 :   Mn_meth3/s_th^2;

compare with method 2

MfiSQ_th_meth2;

is (equal (MfiSQ_th_meth3, MfiSQ_th_meth2));

### Explicit Matrix Traces and Contractions: the m_tr(...) Method

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.

 m_tr(p2 - m, mu, p1 + m, nu) is converted into the explicit matrix expression:

mat_trace ( (sL(p2) - m*I4)  . Gam[mu] . ( sL(p1) + m*I4)  . 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, p2[0] = E, p2[1] = 0, p2[2] = 0, p2[3] = - E,  and
p2 is a "zero based array" holding the components of the
4-momentum vector p2. The components were defined
above using the comp_def function.

sL(p1) is the Feynman slash  matrix corresponding to 
the 4-momentum p1, and I4 is the explicit 4 x 4 unit matrix.

Gam[mu] is the explicit matrix representation of γ^μ, for example.

### Use mcon (expr, index1,index2,..) for traces produced using m_tr or mat_trace

Mn : mcon (m_tr(p2 - m, mu, p1 + m, nu) * m_tr(p3 + M, mu, p4 - M, nu), mu,nu);

replace sin(th)^2  using the Dirac3 package function ts (expr, th)

ts (Mn, th), expand;

Mn_meth4 : pullfac (%, 32);

MfiSQ_th_meth4 :   Mn_meth4 / s_th^2;

compare with method 3

is (equal (MfiSQ_th_meth4, MfiSQ_th_meth3));

A;

dsigdo_meth4 : (A/4) * MfiSQ_th_meth4;

compare differential scattering cross section with method 2

is (equal (dsigdo_meth4, dsigdo_meth2));

what is the high energy limit?

dsigdo_meth4, k = E, p = E, m = 0, M = 0;

factor(%);

which agrees with the high energy limit.

### The Explicit Matrix Method Using mat_trace plus mcon 

Instead of using the convenient m_tr syntax (which follows the symbolic tr syntax), 
we can work directly with explicit matrices to get the same final results.

tr1 : mat_trace ( (sL(p2) - m*I4) . Gam[mu] . (sL(p1) + m*I4)  . Gam[nu]);

tr2 : mat_trace ( (sL(p3) + M*I4) . Gam[mu] . (sL(p4) - M*I4) . Gam[nu]);

mcon (mat_trace ( (sL(p2) - m*I4) . Gam[mu] . (sL(p1) + m*I4)  . Gam[nu]) * 
                      mat_trace ( (sL(p3) + M*I4) . Gam[mu] . (sL(p4) - M*I4) . Gam[nu]), mu,nu);

ts(%, th),  expand;

Mn_meth5 : pullfac (%, 32);

is (equal (Mn_meth5, Mn_meth4));

## POLARIZED AMPLITUDES USING EXPLICIT DIRAC SPINORS
e(-,p1, σ1) +  e(+,p2, σ2)  -->  μ(-,p3, σ3) + μ(+,p4, σ4)

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

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 (for electrons and
positrons respectively; for muons and antimuons, m --> M).
    
    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.

 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 agree with those specified in the comp_def statement above.
    

Recall from Sec. 1.2.1 that in CMS all particles have energy E, and
p is the magnitude of the 3-momentum of the electron and positron
in the initial state (4-momenta  p1 and p2), while k is the magnitude of the
3-momentum of the muon and anti-muon in the final state
(4-momenta  p3 and p4).

### e(-,p1, R) +  e(+,p2, R)  -->  μ(-,p3, R) + μ(+,p4, R)

We start with RR --> RR case,  bearing in mind that in the
high energy case (considered in ee-mumu-HE.wxm)
we got an amplitude which was 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.
The magnitude of the electron and positron 3-momentum is p.
The magnitude of the muon and anti-muon 3-momentum is k.
 All of the involved particles (in the CMS) have energy E.

up1:UU(E,p,0,0,1);
up3b:sbar(UU(E,k,th,0,1));
vp2b:sbar(VV(E,p,%pi,0,1));
vp4:VV(E,k,%pi-th,%pi,1);

check normalization of vp4; we should have vp4b . vp4 = - 2 M

vp4b : sbar (vp4);

grind(%);

vnorm : vp4b . vp4;

We first simplify this "by hand."

rootscontract(%);

subst(sqrt(E^2 - k^2) = M, %);

trigsimp(%);

We next make use of the package function Ep_to_m (expr,p,E,m) which replaces sqrt(E - p) * sqrt(E + p) by m

Ep_to_m(vnorm, k, E, M);

trigsimp(%);

### Use mcon (expr, index1, index2,...) for contractions of explicit matrix expressions

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 γ^{μ}

This numerator now is for the case RR --> RR by construction of our spinors.

Mn : mcon ( (vp2b . Gam[mu] . up1)*(up3b . Gam[mu] . vp4), mu);

Ep_to_m(%, k, E, M);

Ep_to_m(%, p, E, m);

convert this expression to a function of th instead of the current function of th/2,
using fr_ao2(expr, th)

Mn :  fr_ao2 (%,th);

s_th;

Mfi : Mn/s_th;

note that Mfi is real, so |Mfi|^2 = Mfi*Mfi = Mfi^2

A;

dsigdo_RR_RR : A * Mfi^2;

%, th = theta;

This expression has dimension [dsigdo] = 1/[s] = 1/[E^2} = 1/mass^2.

In the high energy (HE) limit, E >> M >> m , this differential scattering cross section
for the case RR --> RR  becomes small.

### e(-,p1, R) +  e(+,p2, L)  -->  μ(-,p3, R) + μ(+,p4, L)

Next we look at the amplitude for e(-,p1,R) + e(+,p2,L) --> μ(-,p3,R) + μ(+,p4,L)
This is case  RL --> RL

up1:UU(E,p,0,0,1);
up3b:sbar(UU(E,k,th,0,1));
vp2b:sbar(VV(E,p,%pi,0, -1));
vp4:VV(E,k,%pi-th,%pi, -1);

### Use mcon (expr, index1, index2,...) for contractions of explicit matrix expressions

Mn : mcon ( (vp2b . Gam[mu] . up1)*(up3b . Gam[mu] . vp4), mu);

Mfi : Mn / s_th;

Mfi : fr_ao2(Mfi, th);

so the (reduced) amplitude for the case RL -->  RL is ( 1 + cos(θ) ), the same as in the HE case, 
and reminding ourselves of the definition of A,

A;

dsigdo_RL_RL : A * Mfi^2;

%, th = theta;

Even though the amplitude is the same as the HE case, the differential
scattering cross section in the AE case will be less since here (k/p) < 1.

We can express the ratio (k/p) in terms of E, M, and m as:

(k/p) = √(E^2 - M^2) / √(E^2 - m^2)  < 1

### General Helicity Cases

n the following function definition,
s1v is the helicity of the particle with 4-momentum p1 (the incident electron), etc.

Define the function AE_me(s1v,s2v,s3v,s4v) which returns Mfi = Mn/s
for the case that the particle helicities have the specified values.

note: s1v = (1 or -1) for (R, or L) , ie, right-handed or left-handed.
The quantities  t_th2 and s_th  are assumed globally known.
E and th are undefined parameters.

The details of these steps are suggested by our interactive step by
step cases worked out above.

AE_me (s1v,s2v,s3v,s4v) :=
 block ([up1,up3b,vp4,vp2b],
      up1:UU(E,p,0,0,s1v),
      vp2b:sbar(VV(E,p,%pi,0, s2v)),
      up3b:sbar(UU(E,k,th,0,s3v)),
      vp4:VV(E,k,%pi-th,%pi, s4v),         
      mcon ( (vp2b . Gam[mu] . up1)*(up3b . Gam[mu] . vp4), mu),
     Ep_to_m(%%, k, E, M),
     Ep_to_m(%%, p, E, m)/s_th,
      fr_ao2 (%%, th))$

Let's check the case RR --> RR we worked out above.

AE_me(1,1,1,1);

This agrees with Mfi worked out above for the RR --> RR case.

Next check the case RL --> RL

AE_me(1,-1,1,-1);

which again agrees.   

### Polarized amplitudes --> unpolarized cross section

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 AE_me defined in terms of the amplitude expressed
in terms of explicit Dirac spinors to calculate mssq takes about 0.7 - 0.9 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 AE_me.

mssq : 0$
sL : [1, -1]$

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 :  AE_me (s1,s2,s3,s4),                     
                     mssq : mssq + Avsq(temp),
                     print ("  "),
                      print (s1,s2,s3,s4," ",temp ) ),
    print ("   "))$

We see that none of the amplitudes is zero.

time(%);

save this as mssq_DS for future use.

mssq_DS : mssq;

define mssq_ts as an expression involving all cosines

mssq_ts : expand ( ts (mssq,th) );

Let's compare this to MfiSQ_th_E from Sect. 1.3.1

MfiSQ_th_E;

is (equal (mssq_ts, MfiSQ_th_E));

so this method produces an equivalent result.

## Trace Methods and Helicity Projection Matrices for Square of Polarized Amplitudes

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).
 
[see ee-mumu-HE.wxm, where we ignore the masses of the particles, and
must use S(σ) for the massless lepton case and  S(-σ)  for the massless anti-lepton case. ]

### Use of nc_tr plus mcon  for squared amplitudes and cross section using spin projection operators

The new feature for massive fermions and anti-fermions is the spin 4-vector  Sp defined by
the particle's 4-momentum vector.

For a lepton or anti-lepton with mass m > 0, energy E,   3-momentum pvec,
3-momentum magnitude |pvec|, and 3-momentum unit vector
phat = pvec / |pvec|,

Sp = [ |pvec|/m, (E/m) phat ],  or in more detail

Sp = [ |pvec|/m, (E/m) phat_x,  (E/m) phat_y,  (E/m) phat_z ]

and this can be defined via our package using the notation
with comp_def (...), and specializing to our specific case
as a concrete example:

comp_def ( p1(E,0,0,p), Sp1 (p/m),0,0,E/m),
                      p2 (E,0,0, -p),  Sp2 (p/m, 0,0,-E/m),
                      p3 (E, k*sin(th), 0, k*cos(th)), Sp3 (k/M, E*sin(th)/M, 0, E*cos(th)/M),
                      p4 (E, -k*sin(th), 0, - k*cos(th)), Sp4 (k/M, - E*sin(th)/M, 0, - E*cos(th)/M) )

Since we have already used comp_def to define the components of the particle
4-momentua p1, p2, p3, and p4, we only need to define the components of
the four vectors Sp1, Sp2, Sp3, and Sp4.

comp_def (  Sp1 (p/m,0,0,E/m),
                   Sp2 (p/m, 0,0,-E/m),
                  Sp3 (k/M, E*sin(th)/M, 0, E*cos(th)/M),
                  Sp4 (k/M, - E*sin(th)/M, 0, - E*cos(th)/M) )$

If we use the symbolic nc_tr(...) method,  we should make use of an invarR list of
4-momentum product replacements.

set_invarR(D(p1,p1) = m^2,
             D(p2,p2) = m^2,
             D(p3,p3) = M^2,
             D(p4,p4) = M^2,
             D(p1,Sp1) = 0,
             D(Sp1,Sp1) = -1,
             D(p2,Sp2) = 0,
             D(Sp2,Sp2) = -1,
             D(p3,Sp3) = 0,
             D(Sp3,Sp3) = -1,
             D(p4,Sp4) = 0,
             D(Sp4,Sp4) = -1)$

invarR;

case RR --> RR

MSQn : mcon (nc_tr (mu,S(1,Sp1),p1 + m,nu,S(1,Sp2),p2 - m) * 
                                      nc_tr (mu,S(1,Sp4),p4 - M,nu,S(1,Sp3),p3 + M), mu,nu);

time(%);

MSQn : ts(MSQn, th), expand;

MSQn :  psq_to_Esq(MSQn,p,E,m), expand;

MSQn :  psq_to_Esq(MSQn ,k,E,M), expand;

MSQn/ s_th^2;

this agrees with  the RR --> RR case
computed in Sec. 2.1 using explicit Dirac spinors.

general amplitude squared for specified helicities using nc_tr method,
constructed using our experience simplifying the RR --> RR case above.

meSQ(s1v,s2v,s3v,s4v) := 
( mcon (nc_tr (mu, S(s1v,Sp1),p1 + m, nu, S(s2v,Sp2), p2 - m) * 
                                                    nc_tr (mu,S(s4v,Sp4), p4 - M, nu, S(s3v,Sp3), p3 + M), mu, nu),
    expand (ts (%%, th)),   expand (psq_to_Esq (%%,p,E,m)), factor (expand (psq_to_Esq(%%, k,E,M))),
    %%/ s_th^2)$

case RR --> RR

meSQ(1,1,1,1);

time(%);

The m_tr method is more than eight times faster at this job: see below.

### Use of m_tr plus mcon  for squared amplitudes and cross section using spin projection operators

The m_tr method converts the input to explcit matrices, doesn't call tr or nc_tr, hence doesn't care
about the contents of the invarR list.  The syntax of the m_tr method is formally the same
as the nc_tr method above. In general, the m_tr method is faster than the nc_tr method.


case RR --> RR

MSQn : mcon (m_tr (mu,S(1,Sp1),p1 + m,nu,S(1,Sp2),p2 - m) * 
                                      m_tr (mu,S(1,Sp4),p4 - M,nu,S(1,Sp3),p3 + M), mu,nu);

time(%);

ts(MSQn, th), expand;

psq_to_Esq (%,p,E,m),  expand;

psq_to_Esq (%,k,E,M),  expand;

%/ s_th^2;

this agrees with  the reduced amplitude for the RR --> RR case
computed in Sec. 2.1 using explicit Dirac spinors.

mesq_mtr : meSQ_mtr(1,1,1,-1);

mesq_nctr : meSQ(1,1,1,-1);

mesq_DS : AE_me(1,1,1,-1)^2;

ratsimp(mesq_mtr);

ts(%,th);

general amplitude squared for specified helicities using m_tr method,
constructed using our experience simplifying the RR --> RR case above.

meSQ_mtr (s1v,s2v,s3v,s4v) := 
( mcon (m_tr (mu, S(s1v,Sp1),p1 + m, nu, S(s2v,Sp2), p2 - m) * 
                                                    m_tr (mu,S(s4v,Sp4), p4 - M, nu, S(s3v,Sp3), p3 + M), mu, nu),
      expand (psq_to_Esq (%%,p,E,m)), expand (psq_to_Esq(%%, k,E,M)), trigsimp(%%), factor(%%),
   /*  ratsimp (%%),  */   %%/ s_th^2)$

check its operation for RR --> RR case

meSQ_mtr (1,1,1,1);

time(%);

case RR --> RL

meSQ_mtr (1,1,1,-1);

case RL --> RL

meSQ_mtr (1,-1,1,-1);

print out all reduced squared amplitudes and comput the sum using m_tr method

mssq : 0$
sL : [1, -1]$

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_mtr (s1,s2,s3,s4),
                     mssq : mssq + temp,
                     print ("  "),
                      print (s1,s2,s3,s4," ",temp ) ),
    print ("   "))$

time(%);

mssq;

expand(mssq);

compare this with mssq_DS computed above using explicit Dirac spinors

is (equal (%, mssq_DS));

which shows that the m_tr method agrees with the explcit Dirac spinor method.

### Use of mat_trace plus mcon  for squared amplitudes using spin projection operators

### Using mat_trace to get the square of the RR --> RR amplitude

/*  MSQn : mcon (m_tr (mu,S(1,Sp1),p1 + m,nu,S(1,Sp2),p2 - m) * 
                                      m_tr (mu,S(1,Sp4),p4 - M,nu,S(1,Sp3),p3 + M), mu,nu); */

MSQn : mcon (mat_trace (Gam[mu] . P(1,Sp1) . (sL(p1) + m*I4) . Gam[nu] . P(1,Sp2) . (sL(p2) - m*I4)) * 
                             mat_trace (Gam[mu] . P(1,Sp4) . (sL(p4) - M*I4) . Gam[nu] . P(1,Sp3) . (sL(p3) + M*I4) ), mu,nu);

ts(%, th), expand;

psq_to_Esq (%,p,E,m),  expand;

psq_to_Esq (%,k,E,M),  expand;

% / s_th^2;

which agrees with the square of the reduced RR --> RR amplitude evaluated using
explicit Dirac spinors. We see that we get the advantage of the speed of explict Dirac matrix
traces with a simple notation if we stick to the m_tr + mcon method.

### Explicit Spin Projection Matrix P(sv, Sp) for the Arbitrary Energy Case

The following section is purely for pedagogical reasons and may be skipped.

In general,  P(s,  Sp) is the explicit spin projection matrix  (I4 + s*Gam[5]  .  sL(Sp)   ) / 2 
where s is restricted to the values [1, -1] for our purposes.

P(s, Sp);

is (equal (P(s,Sp), (I4 + s*Gam[5]  .  sL(Sp)   ) / 2 ));

We repeat here comments about massive lepton spin vectors.

The new feature for massive fermions is the spin 4-vector  Sp defined by
the particle's 4-momentum vector.

For a lepton or anti-lepton with mass m > 0, energy E,   3-momentum pvec,
3-momentum magnitude |pvec|, and 3-momentum unit vector
phat = pvec / |pvec|,

Sp = [ |pvec|/m, (E/m) phat ],  or in more detail

Sp = [ |pvec|/m, (E/m) phat_x,  (E/m) phat_y,  (E/m) phat_z ]

and this can be defined via our package using the notation
with comp_def (...), and specializing to our specific case
as a concrete example:

comp_def ( p1(E,0,0,p), Sp1 (p/m),0,0,E/m),
                      p2 (E,0,0, -p),  Sp2 (p/m, 0,0,-E/m),
                      p3 (E, k*sin(th), 0, k*cos(th)), Sp3 (k/M, E*sin(th)/M, 0, E*cos(th)/M),
                      p4 (E, -k*sin(th), 0, - k*cos(th)), Sp4 (k/M, - E*sin(th)/M, 0, - E*cos(th)/M) )

In general, for MASSIVE  fermions with mass m > 0 we then have
in explicit matrix notation

u(p,σ) .  ubar(p,σ) = P(σ, Sp) . (sL(p) + m*I4) = (sL(p) + m*I4) . P(σ, Sp)

and for MASSIVE anti-fermions with mass m > 0 we have using 
explicit matrix notation

v(p,σ) .  vbar(p,σ) = P(σ, Sp) . (sL(p) - m*I4)  = (sL(p) - m*I4) . P(σ, Sp)

 Example: RR --> RR case; we first repeat our definition of the relevant Dirac spinors 

up1:UU(E,p,0,0,1)$
up3b:sbar(UU(E,k,th,0,1))$
vp2b:sbar(VV(E,p,%pi,0,1))$
vp4:VV(E,k,%pi-th,%pi,1)$

First consider the square matrix  up1 . up1b , where up1 = u(p,R) = u(p,σ1 = 1).
We need to first define up1b  ( ubar(p,R)  ):

up1b : sbar(up1)$

uub : up1 . up1b;

replace the products sqrt(E-p)*sqrt(E+p)

uub : Ep_to_m(uub,p,E,m);

now define the "right hand side"  of the asserted identity in terms of the appropriate massive
fermion projection operator corresponding to the right-handed electron.

Pp1m :   ratsimp ( P(1, Sp1) . (sL(p1) + m*I4) );

replace E^2 in favor of p^2

Pp1m : Esq_to_psq(Pp1m,p,E,m);

is(equal (uub, Pp1m));

Next consider the square matrix vp4 . vp4b  where vp4 = v(k,R) = v(k, σ4 = 1).
We have already defined vp4 above corresponding to
the right handed anti-muon with four momentum p4
energy E,  3-momentum magnitude k, and mass M.

To confirm the asserted identity involving the projection
operator form will take some massaging of the raw matrix
we initially are presented with.

vp4b : sbar(vp4)$

vvb : vp4 . vp4b;

again we replace the products sqrt(E - k) * sqrt(E + k) in this matrix expression

vvb : Ep_to_m(vvb,k,E,M);

we next convert from f(th/2) to g(th) using fr_ao2( expr, th)

vvb : fr_ao2(vvb,th);

vvb : ratsimp(vvb);

Pp4M :   ratsimp (expand( P(1, Sp4) . (sL(p4) - M*I4) ));

Pp4M : trigsimp(Pp4M);

replace E^2 in favor of k^2 and M^2

Pp4M : Esq_to_psq(Pp4M,k,E,M);

put matrix elements over common denominators:

Pp4M : ratsimp(Pp4M);

is(equal (vvb,Pp4M));

which affirms the asserted identity.

When we evaluate traces symbolically using nc_tr(...),  and also when we
use the symbolic notation for the explicit matrix method m_tr(...), we use the symbolic
projection operator S(σ, Sp) as an argument, which corresponds to the
explicit matrix P(σ, Sp), where Sp is the spin 4-vector corresponding to
the particle 4-momentum p.