# em-ep-scatt-HE.wxmx     High Energy Limit of

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

em-ep-scatt-HE.wxmx,  Oct. 30, 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 em-ep-scatt-HE.wxm works out the unpolarized differential scattering
  cross section in the high energy 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, 
 comp_def, VP, sub_stu, noncov,  Esq_to_psq, pullfac

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

M(σ1, σ2, σ3, σ4) = -e^2 Mr 
 where the "reduced amplitude" is
 
Mr = (u3bar γ^μ u1) ( v2bar γ_μ v4)/ t   -   (v2bar γ^μ u1) ( u3bar γ_μ v4)/ s  = M1 - M2

summed over the Lorentz index μ , where  t = (p1- p3)^2 = (p2 - p4)^2, 
 and s = (p1 + p2)^2 = (p3 + p4)^2.

We call Mr the "reduced" amplitude.

### Unpolarized Cross Section USING Mandelstam Variables s, t, and u, and using tr and Con

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.

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

In the high energy limit, we neglect the mass m of the electron and positron.

This case of the command for set_invarR can be pasted into your
work from the package file mandelstam.txt.

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

The sum over all helicities of the absolute value squared
of Mfi is:
<|Mfi|^2> =  Σ(σ1, σ2, σ3, σ4)  | Mfi(σ1, σ2, σ3, σ4) |^2
                 = e^4  Σ(σ1, σ2, σ3, σ4)  | Mr(σ1, σ2, σ3, σ4) |^2

The reduced amplitude Mr is defined by

M(σ1, σ2, σ3, σ4) = -e^2 Mr, and
 
Mr = (u3bar γ^μ u1) ( v2bar γ_μ v4)/ t   -   (v2bar γ^μ u1) ( u3bar γ_μ v4)/ s 
      = M1 - M2.

In each term a sum over the repeated Lorentz index μ is assumed, 
and  t = (p1- p3)^2 = (p2 - p4)^2,   s = (p1 + p2)^2 = (p3 + p4)^2. 
 
When we take the absolute value squared of Mr for given helicities,
  we get four terms, and we write the sum over helicities of the
  absolute square of Mr as:

  < |Mr|^2>  =   M1n/t^2  +  M2n/s^2 -  M12n/(t*s) - M21n/(t*s) 

   M1n = t^2 * < |M1|^2 >,
   M2n = s^2 * < |M2|^2 >, 
   M12n = (t*s) * < M1*conj(M2) >,
   M21n = (t*s) * < M2*conj(M1) >

The Dirac package function tr(...) carries out the trace of a product of
Dirac gamma matrices symbolically.  We use this method to allow
the incorporation of the Mandelstam variables in the calculation.

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

sL(p3) + m*I4  =  Feynman-slash(p3) + m*I4.
Since we are neglection the mass of the electron, we replace
an arg p3 + m by just p3, for example.

The arg mu inside tr function stands for the square Dirac matrix γ^μ
which is Gam[mu] if we were using explicit Dirac matrices in the calculation.

We let MfiSQ stand for  <|Mr|^2>.

We use Con (expr, mu,nu) for the case of contraction of an expression
expr containing two lorentz indices mu and nu.

M1n : factor (Con ( tr (mu,p1,nu,p3) * tr(mu,p4,nu,p2), mu,nu));

M2n : factor (Con ( tr(mu,p1,nu,p2)*tr(mu,p4,nu,p3), mu,nu));

In the next two steps, notice how automatic contraction is done for
repeated Lorentz indices within the same tr(...) expression.

M12n : factor ( tr (mu,p1,nu,p2,mu,p4,nu,p3));

M21n : factor ( tr (mu,p1,nu,p3,mu,p4,nu,p2));

MfiSQ :  M1n/t^2 + M2n/s^2 - M12n/(t*s) - M21n/(t*s);

MfiSQ : pullfac (MfiSQ,8);

Since  pf = pi, our differential scattering cross section prefactor
 is just A = α^2/( 4 s).

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

Averaging over initial spins means we need to divide A by 4
 to get the unpolarized differential cross section (CM, HE)

dsigdo_unpol_CM_HE : (A/4)*MfiSQ;

which agrees with Renton's function of s,t, and u in Sec. 4.3.2,
page 159,  (4.57),  (Peter Renton, Electroweak Interactions, Cambridge Univ. Press, 1990)

### 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 p be  the magnitude of the 3-momentum of the initial electron
and positron, and of the final electron and positron.

Since we are neglecting the mass of the electrons, the value of p is
the same as the energy E of each electron

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

These are the CMS 4-momentum components.

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

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, the scattering angle of the electron.

VP(...) is our "4-vector product" function.  For example

VP(pa,pa);

Use VP to compute CMS values of the Mandelstam variables s, t, and u.
We call them  s_th, t_th, and u_th here, notation which is later used by our
package function sub_stu.

  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

above, we had MfiSQ in terms of Mandelstam variables.
Let MfiSQ_th be the same quantity, but expressed in terms of
the scattering angle th (of the electron in the CMS).

MfiSQ_th : sub_stu(MfiSQ), factor;

dsigdo_HE : (A/4) * MfiSQ_th;

which agrees with Peter Renton's HE result on
 p. 160, Eq. (4.58b).

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

taylor (%,th,0,2);

hence singular behavior:   dsigdo ~  1/θ^4  near θ = 0 (i.e., near cos(th) = 1).

### Unpolarized Cross Section  WITHOUT  Mandelstam Variables, using tr and Con

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) = 0,         
         D(p2,p2) = 0,         
         D(p3,p3) = 0,         
         D(p4,p4) = 0)$

invarR;

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

When we take the absolute value squared of Mr for given helicities,
  we get four terms, and we write the sum over helicities of the
  absolute square of Mr as:

  < |Mr|^2>  =   M1n/t^2  +  M2n/s^2 -  M12n/(t*s) - M21n/(t*s) 

   M1n = t^2 * < |M1|^2 >,
   M2n = s^2 * < |M2|^2 >, 
   M12n = (t*s) * < M1*conj(M2) >,
   M21n = (t*s) * < M2*conj(M1) >

M1n : factor (Con ( tr (mu,p1,nu,p3) * tr(mu,p4,nu,p2), mu,nu));

M2n : factor (Con ( tr(mu,p1,nu,p2)*tr(mu,p4,nu,p3), mu,nu));

note automatic trace for repeated indices inside a tr(....) function in the next two pieces.

M12n : factor ( tr (mu,p1,nu,p2,mu,p4,nu,p3));

M21n : factor ( tr (mu,p1,nu,p3,mu,p4,nu,p2));

note we get D(pa,pb)'s  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.

M1n_nc : noncov(M1n);

M2n_nc : noncov(M2n);

M12n_nc : noncov(M12n);

M21n_nc : noncov(M21n);

MfiSQ_th_meth2 :  M1n_nc/t_th^2 + M2n_nc/s_th^2 - M12n_nc/(t_th*s_th) - M21n_nc/(t_th*s_th), factor;

this is the same result as MfiSQ_th which we obtained in the first method
(Con (tr(...))  via Mandelstam variable route.

is(equal (MfiSQ_th_meth2, MfiSQ_th));

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_th and u_th
using VP.

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(...) and mcon

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

M1n : factor (mcon ( nc_tr (mu,p1,nu,p3) * nc_tr(mu,p4,nu,p2), mu,nu));

M2n : factor (mcon ( nc_tr(mu,p1,nu,p2)*nc_tr(mu,p4,nu,p3), mu,nu));

note automatic trace for repeated indices inside a nc_tr(....) function in the next two pieces.

M12n : factor ( nc_tr (mu,p1,nu,p2,mu,p4,nu,p3));

M21n : factor ( nc_tr (mu,p1,nu,p3,mu,p4,nu,p2));

MfiSQ_th_meth3 :  M1n/t_th^2 + M2n/s_th^2 - M12n/(t_th*s_th) - M21n/(t_th*s_th);

 MfiSQ_th_meth3 :   trigsimp(MfiSQ_th_meth3), factor;

is(equal (MfiSQ_th_meth3, MfiSQ_th));

which shows that method 3 using mcon ( nc_tr(...) ) gives the same result as method 1 using Con (tr (...)  )
with Mandelstam variables.

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

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.

listarray (p3);

grind(%)$

m_tr ( mu, p1,  nu, p3 ,=) is converted into the explicit matrix expression:

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

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] = E, p3[1] = E*sin(th), p3[2] = 0, p3[3] = E*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.

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

M1n : factor (mcon ( m_tr (mu,p1,nu,p3) * m_tr(mu,p4,nu,p2), mu,nu));

M2n : factor (mcon ( m_tr(mu,p1,nu,p2) * m_tr(mu,p4,nu,p3), mu,nu));

When using m_tr method, there is no automatic contraction of repeated indices inside m_tr(...),
so we must explicitly call mcon(expr, indices)

M12n : factor ( mcon (  m_tr (mu,p1,nu,p2,mu,p4,nu,p3), mu, nu));

For example, here is what you get from m_tr in this step, if you
do NOT call mcon for contraction of repeated indices:

m_tr (mu,p1,nu,p2,mu,p4,nu,p3);

M21n : factor ( mcon (  m_tr (mu,p1,nu,p3,mu,p4,nu,p2), mu, nu));

MfiSQ_th_meth4 :  M1n/t_th^2 + M2n/s_th^2 - M12n/(t_th*s_th) - M21n/(t_th*s_th);

MfiSQ_th_meth4 : trigsimp (MfiSQ_th_meth4), factor;

is(equal (MfiSQ_th_meth4, MfiSQ_th));

which confirms the equality of results of method 3 using mcon ( nc_tr (...),  indices  )
with method 4 using mcon ( m_tr(...),  indices )

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, VV, sbar, Gam, mcon, sub_stu, Avsq, ts, to_ao2,  fr_ao2, Esq_to_psq, Ep_to_m

 
   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.
    
 Our normalization is  upa_bar . upa = 2m and  vpa_bar . vpa = - 2m (for electrons and
positrons respectively  (in this high energy case 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.

 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 electron and positron.
    
    These parameters should agree with those specified in the comp_def statement above.
E(CM) =  E + E = 2 E,  and   s = E(cm)^2 = 4 E^2

Define  th/2 forms  for t and u  using to_ao2 (expr, th)

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

The amplitude for specific helicities σ1, σ2, σ3, σ4  is (Griffiths, p235)

M(σ1, σ2, σ3, σ4) = -e^2 Mr 
 where the "reduced amplitude" is
 
Mr = (u3bar γ^μ u1) ( v2bar γ_μ v4)/ t   -   (v2bar γ^μ u1) ( u3bar γ_μ v4)/ s  = M1 - M2

summed over the Lorentz index μ , where  t = (p1- p3)^2 = (p2 - p4)^2, 
 and s = (p1 + p2)^2 = (p3 + p4)^2.

We call Mr the "reduced" amplitude.

case  RR --> RR, ie, (+1,+1) --> (+1,+1) polarized amplitude.
               
     First define the needed Dirac spinors and barred spinors for leptons.

th is the electron scattering angle; here p = E since
we neglect the mass m of the electron.

up1 : UU (E,E,0, 0, 1);
up3b : sbar (UU (E,E,th, 0, 1));
vp4 : VV (E,E,%pi-th, %pi, 1);
vp2b : sbar ( VV (E,E,%pi, 0, 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 ( up3b)$

Check our normalization for up1:

up1;

up1b : sbar (up1);

up1b . up1;

which is correct in the limit m --> 0 (the high energy limit).

Now calculate the (reduced) amplitude Mr =  M_RR_RR
 
The amplitude for specific helicities σ1, σ2, σ3, σ4  is (Griffiths, p235)

M(σ1, σ2, σ3, σ4) = -e^2 Mr 
 where the "reduced amplitude" is
 
Mr = (u3bar γ^μ u1) ( v2bar γ_μ v4)/ t   -   (v2bar γ^μ u1) ( u3bar γ_μ v4)/ s  = M1 - M2

summed over the Lorentz index μ , where  t = (p1- p3)^2 = (p2 - p4)^2, 
 and s = (p1 + p2)^2 = (p3 + p4)^2.

We call Mr the "reduced" amplitude.

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

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

M1n : trigsimp (M1n);

M1 : M1n/t_th;

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

M2 : M2n/s_th;

M_RR_RR : M1 - M2;

to_ao2(%,th);

Now automate calculation of the reduced amplitudes for given helicity
assignments, defining the function Amplitude(s1v,s2v,s3v,s4v)

in this definition,  E is an unbound parameter, and  t_th and s_th have been  defined above.

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

Amplitude (s1v,s2v,s3v,s4v) := 
  block ([up1,vp4,up3b,vp2b, M1, M2],
    up1 : UU (E, E, 0, 0, s1v),
    up3b : sbar (UU (E, E, th, 0, s3v)),
    vp2b : sbar ( VV (E, E, %pi, 0, s2v)),
    vp4 : VV (E, E, %pi - th, %pi, s4v),
    trigsimp(   mcon ( (up3b . Gam[mu] . up1) * (vp2b . Gam[mu] . vp4), mu)),
    M1 : fr_ao2 (%%, th) / t_th,
    trigsimp (   mcon ( (vp2b . Gam[mu] . up1) * (up3b . Gam[mu] . vp4), mu)),
    M2 :  fr_ao2 (%%, th) / s_th,
    ratsimp( M1 - M2 ))$

Try this function out on the case RR --> RR

Amplitude(1,1,1,1);

is(equal(%, M_RR_RR)) ;

## 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 is the sum over spins of the absolute value squared of the (reduced) matrix elements using our Avsq function.
mssq is global parameter, given a zero value at start.
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 :  Amplitude (s1,s2,s3,s4),                     
                     mssq : mssq + Avsq(temp),
                     print ("  "),
                      print (s1,s2,s3,s4," ",temp ) ),
    print ("   "))$

we have the following (raw) result for mssq:

mssq;

mssq :  trigsimp (mssq);

mssq : scanmap ('factor, mssq);

is (equal (mssq, MfiSQ_th));

dsigdo_DS : (A/4) * mssq;

is (equal (dsigdo_DS, dsigdo_HE));

which confirms the equivalence of using the explicit Dirac spinor amplitude method
to the more usual  trace, contraction method of obtaining the differential scattering
cross section.

## 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
  (which is actually the case for this worksheet, since we are ignoring the
   mass m of the electron here).

[See the worksheet em-ep-scatt_AE.wxm for discussion of the lepton
  spin projection operator S(σ, Sp ) needed for massive particles.]

[  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 ]   ]

For the massless case (HE) we are considering,  no definition of
particle spin 4-vector is available or needed.

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

What is the current value of the list invarR?

invarR;

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

We write |Mr|^2 = M1SQn/t^2 + M2SQn/s^2 - M12SQn/(t *s) - M21SQn/(t *s), 

  Since all particles are considered massless,
we have, for example,  u1 . u1b = P(σ1) . sL(p1) and v2 . v2b = P(-σ2) . sL(p2) in explicit matrix notation.

we consider then the case RR --> RR first:

 M1SQn :   mcon ( nc_tr (mu,S(1),p1,nu,S(1),p3) *nc_tr(mu,S(-1),p4,nu,S(-1),p2), mu,nu);

time(%);

M1SQn : trigsimp(M1SQn);

M1SQ : M1SQn / t_th^2;

 M2SQn :  mcon ( nc_tr (mu,S(1),p1,nu,S(-1),p2) *nc_tr(mu,S(-1),p4,nu,S(1),p3), mu,nu);

M2SQ : M2SQn / s_th^2;

Sinc nc_tr uses the symbolic tr(..) function (followed by noncov on each term),
and the symbolic tr(...) function automatically contracts on repeated indices inside
a given trace, we don't need to call mcon for contraction for this term.

M21SQn :   nc_tr (mu,S(1),p1,nu,S(1),p3,mu,S(-1),p4,nu,S(-1),p2);

M21SQ : M21SQn / (t_th * s_th);

M12SQn : nc_tr (mu,S(1),p1,nu,S(-1),p2,mu,S(-1),p4,nu,S(1),p3);

M12SQ : M12SQn/(t_th * s_th);

MfiSQ : M1SQ + M2SQ - M12SQ - M21SQ;

we compare this with the square (the amplitudes here are real numbers)
of the Dirac spinor amplitude given by the function Amplitude defined
above.

Amplitude(1,1,1,1)^2;

and we have agreement.

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

The explicit matrix method m_tr(...)  has the same syntax as the symbolic method nc_tr(...)

M1SQn :   mcon ( m_tr (mu,S(1),p1,nu,S(1),p3) *m_tr(mu,S(-1),p4,nu,S(-1),p2), mu,nu);

time(%);

We see that the m_tr(...) method for M1SQn is about 17 times faster than the nc_tr(...) method for this example.

M1SQn : trigsimp(M1SQn);

M1SQ : M1SQn / t_th^2;

M2SQn :  mcon ( m_tr (mu,S(1),p1,nu,S(-1),p2) *m_tr(mu,S(-1),p4,nu,S(1),p3), mu,nu);

M2SQ : M2SQn / s_th^2;

Because m_tr(...) doesn't call the symbolic trace function tr(...), no automatic contractions
are carried out in repeated Lorentz indices inside m_tr(...), and we must explicitly use
mcon (expr, indices) for contraction.

M21SQn :   mcon (   m_tr (mu,S(1),p1,nu,S(1),p3,mu,S(-1),p4,nu,S(-1),p2), mu,nu);

M21SQ : M21SQn / (t_th * s_th);

M12SQn :  mcon (  m_tr (mu,S(1),p1,nu,S(-1),p2,mu,S(-1),p4,nu,S(1),p3), mu, nu );

M12SQ : M12SQn / (t_th * s_th), ratsimp;

MfiSQ :  M1SQ + M2SQ - M12SQ - M21SQ;

which produces the same value for MfiSQ as when using nc_tr.

### Case RR --> LR

The Dirac spinor amplitude for this case is 0.

We have σ1 = 1, σ2 = 1, σ3 = -1, σ4 = 1

Amplitude(1,1,-1,1)^2;

M1SQn :   mcon ( m_tr (mu,S(1),p1,nu,S(-1),p3) *m_tr(mu,S(-1),p4,nu,S(-1),p2), mu,nu);

M2SQn :  mcon ( m_tr (mu,S(1),p1,nu,S(-1),p2) *m_tr(mu,S(-1),p4,nu,S(-1),p3), mu,nu);

M21SQn :   mcon (   m_tr (mu,S(1),p1,nu,S(-1),p3,mu,S(-1),p4,nu,S(-1),p2), mu,nu);

M12SQn :  mcon (  m_tr (mu,S(1),p1,nu,S(-1),p2,mu,S(-1),p4,nu,S(-1),p3), mu, nu );

so we  get a zero value for the square of the RR -->  LR amplitude, in agreement with
the Dirac spinor method.