# em-em-scatt-AE.wxmx 

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

em-em-scatt-AE.wxmx,  Nov. 2, 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-em-scatt-AE.wxm works out the unpolarized differential scattering
  cross section, 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, psq_to_Esq, 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" is
 
Mr = (u3bar γ^μ u1) ( u4bar γ_μ u2)/ t   -   (u4bar γ^μ u1) ( u3bar γ_μ u2)/ u  = M1 - M2


summed over the Lorentz index μ , where  t = (p1- p3)^2 = (p2 - p4)^2, 
 and u = (p1- p4)^2 = (p2 - p3)^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

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

invarR;

In lowest order perturbation theory, two diagrams involving the exchange of a photon
between the electrons contribute.

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) ( u4bar γ_μ u2)/ t   -   (u4bar γ^μ u1) ( u3bar γ_μ u2)/ u  
      = M1 - M2.

In each term a sum over the repeated Lorentz index μ is assumed, 
and  t = (p1- p3)^2 = (p2 - p4)^2,   u = (p1- p4)^2 = (p2 - p3)^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/u^2 -  M12n/(t*u) - M21n/(t*u) 

   M1n = t^2 * |M1|^2, M2n = u^2 * |M2|^2, 
   M12n = (t*u) * M1*conj(M2),
   M21n = (t*u) * 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.

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 (p3+m,mu,p1+m,nu)*tr (p4+m,mu,p2+m,nu), mu,nu));

M2n : factor ( Con ( tr (p4+m,mu,p1+m,nu)*tr (p3+m,mu,p2+m,nu), 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 (p3+m,mu,p1+m,nu,p4+m,mu,p2+m,nu));

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

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

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

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

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

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

dsigdo_unpol_CM : (A/4)*MfiSQ;

### 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 each of the initial electrons
and each of the final electrons.

The energy of each electron is E = sqrt (m^2 + p^2)

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

assume (E>0, p > 0, th >= 0, th <= %pi )$
  
  comp_def ( p1( E,0,0,p),             
             p2( E,0,0,-p),             
             p3 (E,p*sin(th),0,p*cos(th)),             
             p4 (E,-p*sin(th),0,-p*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
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,
which we call s_th, t_th, and u_th here:

  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.

MfiSQ_th :  trigsimp ( sub_stu (MfiSQ));

MfiSQ_th : trigsimp (psq_to_Esq (MfiSQ_th, p, E, m)) ;

We now define a quantity called MfiSQ_cmp, inspired by the Moller Scattering result in the
text: Berestetskii, Lifshitz, and Pitaevskii (BLP),  Quantum Electrodynamics, (Landau and Lifshitz Course of Theoretical Physics, Vol. 4),
2nd edition, Pergamon Press, 1982, page 323, Eq. (81.10) (which uses absolute Gaussian E & M  units instead of Heaviside units)

MfiSQ_cmp : (16*(E^2+p^2)^2/p^4)*
               (4/sin(th)^4-3/sin(th)^2+(p^2/(E^2+p^2))^2*(1+4/sin(th)^2));

MfiSQ_cmp_th : trigsimp (psq_to_Esq (MfiSQ_cmp, p, E, m));

We want to show the equivalence of MfiSQ_cmp_th with MfiSQ_th

trigsimp (expand (MfiSQ_th  - MfiSQ_cmp_th));

which confirms the equivalence, taking into account the relation between p, E, and m.

Using this equivalent form, and dividing A by 4, we get the unpolarized differential scattering
cross section in the CMS in the form:

dsigdo_CM : (A/4)*MfiSQ_cmp;

grind(%);

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

Ignore the overall factor alpha^2*(E^2+p^2)^2/(p^4*s), and look at the other factor:

taylor (p^4*(4/sin(th)^2+1)/(E^2+p^2)^2-3/sin(th)^2+4/sin(th)^4,th,0,2);

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

Of course that behavior is obvious from the above form of dsigdo.

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

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/u^2 -  M12n/(t*u) - M21n/(t*u) 

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

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

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

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

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

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

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 :  trigsimp ( M1n_nc/t_th^2 + M2n_nc/u_th^2 - M12n_nc/(t_th*u_th) - M21n_nc/(t_th*u_th));

MfiSQ_th_meth2 : trigsimp (psq_to_Esq (MfiSQ_th_meth2, p, E, m));

we want to show the equivalence of this method 2 with our previous method.

 trigsimp ( expand(MfiSQ_th_meth2 -  MfiSQ_th));

which confirms the equivalence between method 1 and method 2

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.

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 :  mcon ( nc_tr (p3+m,mu,p1+m,nu) * nc_tr (p4+m,mu,p2+m,nu), mu,nu);

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

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

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

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

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

MfiSQ_th_meth3 : trigsimp ( psq_to_Esq (MfiSQ_th_meth3, p, E, m));

is (equal (MfiSQ_th_meth3, MfiSQ_th_meth2));

which shows the equality of methods 2 and 3.

### 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);

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] = E, p3[1] = p*sin(th), p3[2] = 0, p3[3] = p*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 :  mcon ( m_tr (p3+m,mu,p1+m,nu) * m_tr (p4+m,mu,p2+m,nu), mu,nu);

time(%);

M2n :  mcon ( m_tr (p4+m,mu,p1+m,nu) * m_tr (p3+m,mu,p2+m,nu), 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 :   mcon (  m_tr (p3+m,mu,p1+m,nu,p4+m,mu,p2+m,nu), mu,nu);

time(%);

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

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

MfiSQ_th_meth4 : trigsimp( psq_to_Esq (MfiSQ_th_meth4, p, E, m));

is (equal (MfiSQ_th_meth4, MfiSQ_th_meth3));

so we get agreement.

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.
               
     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 electrons  have equal and opposite
    3-momentum vectors with magnitude p, and same for the outgoing electrons.
    
    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,  E = sqrt(m^2 + p^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);

needed Dirac spinors for the RR --> RR case.
th is the scattering angle, p is magnitude of 3-momentum of each electron
in CMS.

up1 : UU(E,p,0,0, 1);
up3b : sbar (UU(E,p,th,0, 1));
up2 : UU (E,p,%pi,0, 1);
up4b : sbar( UU (E,p,%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;

our package function Ep_to_m is designed to convert the product of square
roots here to the mass of the electron m

Ep_to_m(%,p,E,m);

which confirms out normalization choice for the Dirac spinors.

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

M(σ1, σ2, σ3, σ4) = -e^2 Mr 
 where the "reduced amplitude" is
 
Mr = (u3bar γ^μ u1) ( u4bar γ_μ u2)/ t   -   (u4bar γ^μ u1) ( u3bar γ_μ u2)/ u  = M1 - M2
      =    M1n/t  -  M2n/u

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

We call Mr the "reduced" amplitude.

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

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

M1n : fr_ao2(M1n,th);

M1n : psq_to_Esq (M1n, p, E, m);

M1 : M1n/t_th;

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

M2n : fr_ao2(M2n,th);

M2n : psq_to_Esq (M2n, p, E, m);

M2 : M2n/u_th;

M_RR_RR : M1 - M2, ratsimp;

can we get a simpler looking answer?

M_RR_RR : ts (M_RR_RR, th);

pullfac ( expand (num(%)), 2)/ denom(%);

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

in this (complicated) definition,  E and p are an unbound parameters, and  t_th and u_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,up2,up3b,up4b, M1, M2],
    up1 : UU(E,p,0,0,s1v),
    up3b : sbar (UU(E,p,th,0,s3v)),
    up2 : UU (E,p,%pi,0,s2v),
    up4b : sbar( UU (E,p,%pi - th,%pi,s4v)),
    trigsimp(   mcon ( (up3b . Gam[mu] . up1)*(up4b . Gam[mu] . up2), mu)),
     fr_ao2 (%%, th) ,
    psq_to_Esq (%%, p, E, m),
    M1 : %% / t_th,
    trigsimp (   mcon ( (up4b . Gam[mu] . up1)*(up3b . Gam[mu] . up2), mu)),
    fr_ao2 (%%, th),
    psq_to_Esq (%%, p, E, m),
    M2 :  %% / u_th,
    ratsimp( M1 - M2 ),
    ts (%%, th))$

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

A_cmp : Amplitude(1,1,1,1);

is (equal (A_cmp, M_RR_RR));

which shows agreement.

## 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 square of the (reduced) matrix elements using our Avsq function.
mssq is global parameter, given 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 ) ),
    mssq : trigsimp (mssq),
    print ("   "))$

time(%);

we have the following (raw) result for mssq:

mssq;

mssq : trigsimp (psq_to_Esq (mssq, p, E, m));

trigsimp (mssq - MfiSQ_th);

which confirms the equivalence of using the explicit Dirac spinor amplitude method
compared 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.

  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 ]   

Use comp_def to define the four particle spin 4-vectors.

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

We can confirm the normalization of the spin 4-vector for each particle,
and its orthogonality to the corresponding particle 4-momentum vector.

VP(Sp1,Sp1);

psq_to_Esq(%, p,E,m);

VP(Sp1, p1);

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

What is the current value of the list invarR?

invarR;

Let's reset the list invarR to include the particle spin 4-vectors

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

express t_th and  u_th in terms of E instead of p

t_th : psq_to_Esq (t_th,p,E,m);

u_th : psq_to_Esq( u_th, p,E,m);

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

We write |Mr|^2 = M1SQn/t^2 + M2SQn/u^2 - M12SQn/(t*u) - M21SQn/(t*u), specialized to
the case RR --> RR, and the fact that all particles are leptons (not antileptons).

We use, for example:
 
u1 . u1b = P(σ1, Sp1) .  (sL(p1)  + m*I4)   in explicit matrix notation.

Calling the symbolic tr function (by using nc_tr) means the computation will take some time.

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

time(%);

M1SQn :  trigsimp ( psq_to_Esq (M1SQn, p, E, m));

M1SQ : M1SQn / t_th^2, ratsimp;

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

M2SQn :  trigsimp ( psq_to_Esq (M2SQn, p, E, m));

M2SQ : M2SQn / u_th^2, ratsimp;

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 :  trigsimp ( nc_tr(mu,S(1,Sp1),p1+m,nu,S(1,Sp3),p3+m,mu,S(1,Sp2),p2+m,nu,S(1,Sp4),p4+m));

time(%);

M21SQn :  trigsimp ( psq_to_Esq (M21SQn, p,E,m));

M21SQ : M21SQn / (t_th * u_th), ratsimp;

M12SQn :  trigsimp ( nc_tr(mu,S(1,Sp1),p1+m,nu,S(1,Sp4),p4+m,mu,S(1,Sp2),p2+m,nu,S(1,Sp3),p3+m));

time(%);

M12SQn :  trigsimp ( psq_to_Esq (M12SQn, p, E, m));

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

MSQ :  trigsimp ( 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.

MSQ_DS :  trigsimp (Amplitude(1,1,1,1)^2);

MSQ_DS :   trigsimp(MSQ_DS);

trigsimp (expand (MSQ - MSQ_DS));

and we have agreement between the explicit Dirac spinor method and the trace contraction method
for the case RR --> RR, using nc_tr plus mcon.

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

The explicit matrix method m_tr(...)  has the same syntax as the symbolic method nc_tr(...).
The m_tr method does not call the symbolic tr(...) function, nor does it use the list invarR,
and as a result, the m_tr method is much faster than the nc_tr method for a trace and contraction
involving many terms.

For this first expression M1SQn, the m_tr method is about 10 times faster.

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

time(%);

M1SQn  :  trigsimp ( psq_to_Esq (M1SQn, p, E, m));

M1SQ : M1SQn / t_th^2, ratsimp;

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

M2SQn :  trigsimp ( psq_to_Esq (M2SQn, p, E, m));

M2SQ : M2SQn / u_th^2, ratsimp;

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 in this term.

using m_tr plus mcon for this term is about 5 times faster than the nc_tr method.

M21SQn :   trigsimp ( mcon (  m_tr(mu,S(1,Sp1),p1+m,nu,S(1,Sp3),p3+m,mu,S(1,Sp2),p2+m,nu,S(1,Sp4),p4+m), mu, nu));

time(%);

M21SQn :  trigsimp ( psq_to_Esq (M21SQn, p, E, m));

M21SQ : M21SQn / (t_th * u_th), ratsimp;

M12SQn :   trigsimp ( mcon(  m_tr(mu,S(1,Sp1),p1+m,nu,S(1,Sp4),p4+m,mu,S(1,Sp2),p2+m,nu,S(1,Sp3),p3+m), mu,nu));

time(%);

M12SQn :  trigsimp ( psq_to_Esq (M12SQn, p, E, m));

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

MSQ_mtr :   trigsimp ( M1SQ + M2SQ - M12SQ - M21SQ);

Above, we called the corresponding quantity (using nc_tr method) MSQ, so let's
check against the m_tr evaluation method:

expand (MSQ_mtr - MSQ);

which confirms that the m_tr method produces the same value for MSQ as when using nc_tr.

### Case RR --> LR

For the case RR --> LR, 
we have σ1 = 1, σ2 = 1, σ3 = -1, σ4 = 1

We will check the Dirac spinor answer using the
m_tr method.  We need the args  ..... S(-1,Sp3), p3+m, ....
inside m_tr(...) 

MSQ_DS :  trigsimp ( Amplitude(1,1,-1,1)^2);

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

M1SQn  :  trigsimp ( psq_to_Esq (M1SQn, p, E, m));

M1SQ : M1SQn / t_th^2, ratsimp;

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

M2SQn :  trigsimp ( psq_to_Esq (M2SQn, p, E, m));

M2SQ : M2SQn / u_th^2, ratsimp;

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

M21SQn :  trigsimp ( psq_to_Esq (M21SQn, p, E, m));

M21SQ : M21SQn / (t_th * u_th), ratsimp;

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

M12SQn :  trigsimp ( psq_to_Esq (M12SQn, p, E, m));

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

MSQ_mtr :   trigsimp ( M1SQ + M2SQ - M12SQ - M21SQ);

trigsimp ( expand (MSQ_mtr - MSQ_DS));

which confirms the equivalence of the answers returned by the Dirac spinor method
and the m_tr plus mcon method for the square of the real amplitude for the
case RR --> LR.