# e-mu-scatt.wxmx 

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

e-mu-scatt.wxmx,   Oct. 17, 2018
 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, e^2 = 4πα
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.

This first section of e-mu-scatt.wxm works out the unpolarized differential scattering
  cross section for arbitrary energy, 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, 
 comp_def, VP, sub_stu, Esq_to_psq, pullfac

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 the "reduced amplitude" Mr = (u3bar γ^μ u1) ( u4bar γ_μ u2)/ t,
summed over the Lorentz index μ , where  t = (p1-p3)^2 = (p2-p4)^2, 
We call Mr the "reduced" amplitude.

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

POPULATE THE LIST invarR OF symbolic D(pa,pb) 4-VEC DOT PRODUCT VALUES,
  using the set_invarR  function. We have convervation of 4-momenta:  p1 + p2 = p3 + p4, 
and the Mandelstam variables are defined by:
s = (p1+p2)^2 = (p3+p4)^2,
 t = (p1-p3)^2 = (p2-p4)^2, 
u = (p1-p4)^2 = (p2-p3)^2 
These definitions are independent of coordinate frame.

The electron (p1 --> p3) has mass m,
the muon (p2 --> p4) has mass M.

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 - m^2 - M^2)/2,
         D(p1,p3) = (2*m^2 - t)/2,
         D(p1,p4) = (m^2 + M^2 - u)/2,
         D(p2,p2) = M^2,
         D(p2,p3) = (M^2 + m^2 - u)/2,
         D(p2,p4) = (2*M^2 - t)/2,
         D(p3,p3) = m^2,
         D(p3,p4) = (s - m^2 - M^2)/2,
         D(p4,p4) = M^2)$

In lowest order perturbation theory, only one
diagram contributes (exchange of a photon
between electron and muon).

omit the factor -e^2 from Mfi for now, leaving Mfi = Mr.

The sum over all helicities of the absolute value squared
of Mr is:
<|Mfi|^2> =  Σ(σ1, σ2, σ3, σ4)  | Mr(σ1, σ2, σ3, σ4) |^2
which we write as
 <|Mfi|^2> = Mn/t^2 , and Mn is the contraction of the product
of two traces.

The arg p3+m inside our  tr(...)   function stands for the square matrix: 

sL(p3) + m*I4  =  Feynman-slash(p3) + m*I4.

The arg mu inside tr function stands for the square Dirac matrix γ^μ.

We let MfiSQ stand for  <|Mfi|^2>

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

MfiSQ : Mn/t^2;

### Ignore the electron mass for now

We consider this result in the limit that we can neglect the small
electron mass m and concentrate on the numerator.

Mn : Mn, m = 0;

### CONVERSION TO AN  EXPLICIT FUNCTION OF THE 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.

Let k be  the magnitude of the 3-momentum of the initial electron and
initial muon and kp be the magnitude of the final electron and final muon in CMS.
Conservation of energy implies
   sqrt(m^2 + k^2) + sqrt(M^2 + k^2) = sqrt(m^2 + kp^2) + sqrt(M^2 + kp^2)
This implies kp = k, so pf = pi and A = α^2/( 4 s).

If we neglect the mass of the electron, then Ee --> k

Let the energy of the muon be E,
 Then  E^2 = k^2 + M^2, and

Then    E(CM) =  E + k, and  s = E(CM)^2  =  (E + k)^2

### Use comp_def(...) to define CMS 4-momenta components

assume ( k > 0, th >= 0, th <= %pi )$

comp_def ( p1( k,0,0,k),          
             p2( E,0,0,-k),             
             p3 (k,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 VP(...) to define s_th, t_th, and u_th

now define explicit frame dependent functions of th = scattering angle

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

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

Mn1 : sub_stu (Mn);

We want to verify the expression found in Peskin/Schroeder, p 154, eqn (5.64).
This requires some serious massaging of this function of th.
  
dsigdo is proportional to B defined below,  and
this should be equivalent to the factor Mnfac defined just above
the definition of B.

combine the three terms of Mn1 which are NOT proportional to k^2 or k^3 or k^4 
(the first, second, and fifth terms) in an expression
we call parts, and then show that expression is equal to 16*k^4.

parts : 0;

for j in [1,2,5] do parts : parts +part(Mn1,j);

parts;

parts1 : Esq_to_psq(parts,k,E,M);

now let Mn2 be Mn1 with those 3 terms replaced by parts1:

Mn2 is equivalent to Mn1, but in a different form.

Mn2 : Mn1 -parts + parts1;

now  define Mnfac as what remains of Mn2 if we pull out the factor 32*k^2.

Mnfac :  pullfac(Mn2,32*k^2)  / (32*k^2);

In Peskin/Schroeder, p 154, eqn (5.64),  dsigdo is proportional to what we
call  B here,  and this should be equivalent to the factor Mnfac defined just above.

B : (E+k)^2 + (E+k*cos(th))^2 -M^2*(1-cos(th));

let's try to show B is equivalent to Mnfac if we use the information available.
The following should be zero.

factor( trigsimp(B - Mnfac)) ;

well, let's include the relation between the muon energy E, its
magnitude of 3-momentum k, and its mass M.

Esq_to_psq(%,k,E,M);

which shows the equivalence to Peskin/Schroeder p. 154 m --> 0 limit result.

Now for the differential scattering cross section dsigdo in the CMS.

We have absorbed e^4 into the constant A, with e^2 = 4πα.  

In general 2 --> 2 scattering, A = α^2 * (pf/pi) /(4*s), but we have seen above
that conservation of energy implies that pf = pi = k.
 

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

We then go for the unpolaized electron-muon scattering differential cross section
in the CM frame and in the limit m --> 0

Note we return the factor 32 k^2, and use the form B instead of Mnfac.

Averaging over initial spins means we need to divide A by 4
  to get the unpolarized differential scattering cross section.
(we average over the initial spins σ1 and σ2 by summing over
 σ1  and dividing by 2, and like wise for σ2.

dsigdo : (A/4)*32*k^2*B/t_th^2;

dsigdo, th = theta;

 where s = E(cm)^2 = (E+k)^2.

Next  look at high energy limit, k >> M, so M --> 0, E --> k, E(cm) --> 2*k

dsig_HE : dsigdo, M = 0, E = k, s = 4*k^2;

factor(%);

which gives the  "high energy limit"  dsigdo.

Behavior as a function of θ as θ --> 0?

taylor (%,th,0,2);

hence singular behavior:   dsigdo ~  1/θ^4  near θ = 0.

### UNPOLARIZED CROSS SECTION  WITHOUT Using Mandelstam Variables

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.

As above, <|Mfi|^2> = Mn/t^2 , and Mn is the contraction of the product
of two traces.

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

note we get D(pa,pb) since tr(...) works here without knowledge of the invariants involving s, t, u.

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

Mn_meth2 : noncov (Mn_meth2);

### Ignore the electron mass

Mn_meth2 : Mn_meth2, m = 0;

### Compare result with method 1

compare to the method 1 result Mn1 above

Mn_meth1 : Mn1;

use our package trigsimp replacement function ts(expr,th) to convert sines to cosines.

Mn_meth1 :   ts(Mn_meth1,th), expand;

To compare, consider the difference, after replacing E^2 with M^2 + k^2 using our package
function Esq_to_psq (expr, p,E,m)

diff12 : Esq_to_psq (Mn_meth2 - Mn_meth1, k,E,M);

so Mn_meth2 is an equivalent expression but has a different form.

Factor out 32*k^2 from each term using our package function pullfac(expr, stuff)

Mn_meth2 : pullfac(Mn_meth2, 32*k^2);

This form will yield a different form for the differential scattering cross section:

A;

dsigdo_meth2 : (A/4)*Mn_meth2/t_th^2;

look at the high energy limit of this expression

dsigdo_meth2_HE : dsigdo_meth2, M = 0, E = k, s = 4*k^2;

factor(%);

which agrees with the HE limit of method 1.

To summarize, this non-Mandelstam symbolic method needs either an empty list invarR,
or 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, so we can use noncov(...),  and also,  the evaluation of  t^2  which requires 
 first computing  t_th : VP (p1-p3,p1-p3).
Both noncov and VP use the components of the 4-vectors.

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.

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

invarR;

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. The explicit matrix method m_tr(...) is generally
faster than using the symbolic method nc_tr(...).

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

Mn_meth3 : ts (Mn_meth3, th), expand;

neglect electron mass m as before,

Mn_meth3 : Mn_meth3, m = 0;

compare with method 2

expand (Mn_meth3 - Mn_meth2);

so they agree.

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

m_tr uses the same syntax as  tr or nc_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(p3 + m, mu, p1 + m, nu) is converted into the explicit matrix expression:

mat_trace ( (sL(p3) + 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, p3[0] = k, p3[1] = k*sin(th), p3[2] = 0, p3[3] = k*cos(th),  and
p3 is a "zero based array" holding the components of the
4-momentum vector p3. The components were defined
above using the comp_def function.

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

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

listarray (p3);

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

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

Mn_meth4 : ts(Mn_meth4,th),expand;

neglect electron mass m

Mn_meth4 : Mn_meth4, m = 0;

compare with method 3 above

expand (Mn_meth4 - Mn_meth3);

which confirms the equality.

We could also use "method 5", direct use of mat_trace(....) with explicit matrices as arguments.
See ee-mumu-AE.wxmx for an example in that context.

## POLARIZED AMPLITUDES USING EXPLICIT DIRAC SPINORS

 Using explicit Dirac matrices, in this section we use the package
    functions UU, sbar, Gam, mcon, sub_stu, Avsq, ts, to_ao2,  fr_ao2, Esq_to_psq, Ep_to_m

 case  RR --> RR, ie, (+1,+1) --> (+1,+1) polarized amplitude.
    electron 4-momenta: p1 --> p3, 
     muon 4-momenta:   p2 --> p4
      
     First define the needed Dirac spinors and barred spinors for leptons.

   (See Sec. 2 of ee-mumu-AE.wxmx for corresponding comments about
   Dirac spinors for anti-leptons, which use VV(Ea,pa,θ,φ,σ).)

    If upa is a Dirac spinor (column vector) corresponding to 4-momentum pa, then
    sbar(upa) is the barred Dirac spinor (row vector) upa^{+} γ^0, 
    where (something)^{+} indicates Hermitian conjugate.
    
    Our normalization is  upa_bar . upa = 2m.   
    
    UU(Ea,pa,θ,φ,σ) is a generic Dirac spinor (4-element column vector) corresponding to a lepton with 
    relativistic energy  Ea = sqrt(pa^2 + ma^2), 3-momentum magnitude pa, 
    with direction of 3-momentum defined by the spherical polar angles
     θ and φ, and helicity σ (+1 for R, -1 for L).    
    
    In the center of momentum frame, the incident electron and muon have equal and opposite
    3-momentum vectors with magnitude k, and same for the outgoing electron and muon.
    
    These parameters should agree with those specified in the comp_def statement above.

We neglect the mass of the electron, as above, so  Ee --> k.

Let the energy of the muon be E, so E^2 = k^2 + M^2, and

E(CM) =  E + k,  s = (E + k)^2

Define  th/2 forms

t_th2 : to_ao2 (t_th,th);
  
  u_th2 : to_ao2 (u_th,th);

needed Dirac spinors for the RR --> RR case.
th is the scattering angle of the outgoing electron.

up1 : UU(k,k,0,0,1);
up3b : sbar (UU(k,k,th,0,1));
up2 : UU (E,k,%pi,0,1);
up4b : sbar( UU (E,k,%pi - th,%pi,1)) ;

You can use the grind function to see the display2d=false version of  anything.
The version of wxmaxima we are using omits commas from the row vectors above.

grind ( up4b);

Check our normalization for up1:

up1;

up1b : sbar (up1);

up1b . up1;

which is due to our neglect of the mass of the electron.

Check norm of up2 corresponding to the incoming muon.

up2;

up2b : sbar (up2);

up2b . up2;

Ep_to_m(%,k,E,M);

Now calculate the (reduced) amplitude Mr =  M_RR_RR = Mn/t,  for the RR --> RR case.

The reduced amplitude is dimensionless, since the numerator has the dimensions of (ubar.u)^2 = m^2
and the denominator is  t with dimension m^2.

For contraction, we need to use mcon(expr, indices)

Mn : mcon ( (up3b . Gam[mu] . up1)*(up4b . Gam[mu] . up2), mu);

Mn : trigsimp (Mn);

Mr : Mn/t_th2;

Mr : ratsimp (Mr);

Convert to a function of th (rather than a function of th/2) using fr_ao2(expr, th)

Mr_th  : fr_ao2 (Mr, th);

ratsimp(%);

Mr_th : scanmap(factor,%);

High energy limit,  k >> M,  E --> k

Mr_th, E = k;

%, ratsimp;

Now automate calculation of the reduced amplitudes for given helicity
assignments, defining the function emu_me.

in this definition,  k,  E, and M  are unbound parameters, and  t_th2  was defined above.

This produces the reduced amplitude as a function of θ/2 for given
helicity assignments:  sv = +/-1 corresponds to R/L .

emu_me (s1v,s2v,s3v,s4v) := 
  block ([up1,up2,up3b,up4b],
    up1 : UU(k,k,0,0,s1v),
    up3b : sbar (UU(k,k,th,0,s3v)),
    up2 : UU (E,k,%pi,0,s2v),
    up4b : sbar( UU (E,k,%pi - th,%pi,s4v)),
     mcon ( (up3b . Gam[mu] . up1)*(up4b . Gam[mu] . up2), mu),
     trigsimp(%%)/t_th2,
    ratsimp(%%))$

emu_me(1,1,1,1);

is(equal(%, Mr)) ;

## 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 )
  
After printing out the th/2 form of the amplitude, we
use our package function fr_ao2 to convert the
global quantity mssq to a function of th (instead of th/2).

mssq is the sum over spins of the square of the (reduced) matrix elements using our Avsq function.
mssq is global parameter, given zero value at start of block statement.
For each set of helicity quantum numbers (σ1,σ2,σ3,σ4), we print out
their values and the reduced amplitude Mr(σ1,σ2,σ3,σ4).

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

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

We see that when we use the approximation that the electron mass m is zero,
the amplitudes which involve an electron helicity flip are zero.

As a function of th/2, we have the following result for mssq:

mssq;

convert to a function of th:

mssq_th : fr_ao2(mssq, th);

Simplify this and let Ma be the numerator:

scanmap (factor, ratsimp(mssq_th)) ;

Ma : num(%);

mssq, as a function of th, should be compared to the ratio Mn2/t_th^2 computed above
(in Sect. 1.1.3)  with symbolic methods.  Simplify this ratio and let Mb be the numerator:

ratsimp (Mn2/t_th^2);

scanmap (factor, %);

Since the denominators are obviously equal, we need to show the numerators are equal.

Mb : num(%);

expand(Ma - Mb);

To make the comparison easier, replace E^2 by k^2 + M^2, using our Esq_to_psq function:

Esq_to_psq(%,k,E,M);

which confirms the equivalence of the results of the two methods, which produce
different looking answers.

## 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 the
symbolic method using nc_tr(...) plus mcon or the explicit matrix method
 using m_tr(...)  plus mcon. (Do not use tr(...)  and Con for this type of calculation.)

The matrix method using the m_tr syntax is designed to have
  arguments which look the same as the purely symbolic case which
 uses  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.

  For the massive case, Sp is the particle spin 4-vector determined by the 
 particle's mass and 4-momentum components.

In general, 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|, the associated spin 4-vector has components:

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 ]

The muon, with mass M > 0,  has 4-momentum components that
allow the definition of Sp2 and Sp4.

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

The list invarR is only used by nc_tr, not by m_tr.

What is the current value of the list invarR?

invarR;

so we need to reset this list.

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

invarR;

### Case RR --> RR using nc_tr(...) plus mcon

We write |Mr|^2 = MSQn/t^2, and compute MSQn, the numerator for the case RR --> RR

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

time(%);

We are going to simplify this using trigsimp, ts, expand, and Esq_to_psq
to replace E^2 in favor of k^2 and M^2, producing what we call MSQn_c

MSQn_a : trigsimp(MSQn);

MSQn_b :   ts(MSQn_a,th), expand;

MSQn_c :   Esq_to_psq(MSQn_b,k,E,M);

let's compare this numerator with what we got for the RR --> RR case using explicit Dirac spinors.

We repeat our definition of the needed spinors for this case.

up1 : UU(k,k,0,0,1)$
up3b : sbar (UU(k,k,th,0,1))$
up2 : UU (E,k,%pi,0,1)$
up4b : sbar( UU (E,k,%pi - th,%pi,1))$

Mn : mcon ( (up3b . Gam[mu] . up1)*(up4b . Gam[mu] . up2), mu), trigsimp;

Now we simplify the Dirac spinor (reduced) amplitude numerator by using ts(expr, angle), 
 fr_ao2(expr,th) [ to convert  f(th/2)  to  g(th) ], then squaring the (real) result and then using 
  Esq_to_psq to replace E^2 in favor of k^2 and M^2.

Mn_a : ts (Mn,th/2);

Mn_b : fr_ao2(Mn_a,th);

Mn_b2 : expand(Mn_b^2);

Mn_c2 : Esq_to_psq(Mn_b2,k,E,M);

they look the same, but let's check:

is(equal( MSQn_c, Mn_c2));

which confirms equality for this RR --> RR case.

### Case RR --> RR  using m_tr(...) and mcon

The explicit matrix method m_tr(...)  has the same syntax as the symbolic method
nc_tr(...), but is about five times faster in this example.

MSQmn : mcon (m_tr (mu,S(1),p1,nu,S(1),p3) * m_tr (mu,S(1,Sp2),p2+M,nu,S(1,Sp4),p4+M), mu,nu);

time(%);

The use of the symbolic method nc_tr(...) took about five times longer than using
the explicit matrix method m_tr(....) in this example.

is(equal( MSQmn, MSQn));

So we get the same raw answer as using nc_tr(...)

### Case RR --> LR using nc_tr(...) or m_tr(...)  plus mcon

In the m = 0 approximation, amplitudes which involve a helicity
flip by the electron are zero, according to our Dirac spinor
calculation. So we should get zero using the trace + contraction approach.

We write |Mr|^2 = MSQn/t^2, and compute MSQn, the numerator for the case RL --> RR

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

time(%);

The same numerator using m_tr(...):

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

time(%);