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

    e(-,p1,σ1) + e(+,p2,σ2) --> μ(-,p3,σ3) + μ(+,p4,σ4)
   [wxMaxima: title   end   ] */


/* [wxMaxima: comment start ]
ee-mumu-AE.wxm,   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.
   [wxMaxima: comment end   ] */


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


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


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


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

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

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


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


/* [wxMaxima: comment start ]
This first 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.
   [wxMaxima: comment end   ] */


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

These definitions are independent of coordinate frame.
 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
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
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)$
/* [wxMaxima: input   end   ] */


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


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

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


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


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


/* [wxMaxima: comment start ]
We have absorbed e^4 into a 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
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
Then the unpolarized differential scattering cross section is
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
We can recover the high energy limit (m = 0, M = 0, k = p effectively)
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
dsigdo, m = 0, M = 0, k = p;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
which agrees with result in ee-mumu-HE.wxmx.
   [wxMaxima: comment end   ] */


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


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


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


/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
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)) )$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
comp_def creates a separate (zero based)  array for each 4-vector.
Use listarray to see the components
   [wxMaxima: comment end   ] */


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


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


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


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

Use VP(pa,pb) to form the dot product of the pair of 4-vectors. (VP = "vector product").
   [wxMaxima: comment end   ] */


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


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


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


/* [wxMaxima: comment start ]
We use the package trigsimp alternative ts (expr, angle) to replace sin(th) in terms of cos(th):
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
This is dimensionless, since E, k, p each have the dimension of mass in our units.

Look at the high energy limit:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
MfiSQ_th, k = E, p = E, m = 0, M = 0;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
which agrees with the high energy limit.

Returning to the arbitrary energy case:
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: comment start ]
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:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
dsigdo, k = E, p = E, m = 0, M = 0;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
which agrees with the high energy limit.
   [wxMaxima: comment end   ] */


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

   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
Instead of using sub_stu, we use noncov in this sub-section.
   [wxMaxima: comment end   ] */


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


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


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


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


/* [wxMaxima: input   start ] */
Mn_meth2 : factor (Con (tr(p2 - m, mu, p1 + m, nu) * tr(p3 + M, mu, p4 - M, nu), mu,nu) );
/* [wxMaxima: input   end   ] */


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


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


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


/* [wxMaxima: input   start ] */
Mn_meth2_nc : pullfac(%,32);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
to compare this with method 1, express k and p in terms of E, M, and m.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
psq_to_Esq(MfiSQ_th_meth2, k,E,M), expand;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
MfiSQ_th_meth2_E : psq_to_Esq(%,p,E,m), expand;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Remind ourselves what we got for MfiSQ_th in the Mandelstam variable route:
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
psq_to_Esq(MfiSQ_th,k,E,M), expand;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
MfiSQ_th_E : psq_to_Esq (%,p,E,m), expand;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
is (equal (MfiSQ_th_E, MfiSQ_th_meth2_E));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
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:
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
dsigdo_meth2 : (A/4) * factor( MfiSQ_th_meth2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
check the high energy limit
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
dsigdo_meth2, k = E, p = E, M = 0, m = 0, factor;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
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).
   [wxMaxima: comment end   ] */


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


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


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

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


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


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


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


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


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


/* [wxMaxima: input   start ] */
Mn_meth3 : pullfac (%, 32);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
compare with method 2
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
is (equal (MfiSQ_th_meth3, MfiSQ_th_meth2));
/* [wxMaxima: input   end   ] */


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


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

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


/* [wxMaxima: comment start ]
 m_tr(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.
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: comment start ]
replace sin(th)^2  using the Dirac3 package function ts (expr, th)
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
Mn_meth4 : pullfac (%, 32);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
compare with method 3
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
is (equal (MfiSQ_th_meth4, MfiSQ_th_meth3));
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: comment start ]
compare differential scattering cross section with method 2
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
is (equal (dsigdo_meth4, dsigdo_meth2));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
what is the high energy limit?
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
dsigdo_meth4, k = E, p = E, m = 0, M = 0;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
which agrees with the high energy limit.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: input   start ] */
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);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
Mn_meth5 : pullfac (%, 32);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
is (equal (Mn_meth5, Mn_meth4));
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: comment start ]
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.
    
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
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).
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
We start with RR --> RR case,  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.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
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);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
check normalization of vp4; we should have vp4b . vp4 = - 2 M
   [wxMaxima: comment end   ] */


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


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


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


/* [wxMaxima: comment start ]
We first simplify this "by hand."
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
subst(sqrt(E^2 - k^2) = M, %);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
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
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Ep_to_m(vnorm, k, E, M);
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
Ep_to_m(%, k, E, M);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
Ep_to_m(%, p, E, m);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
convert this expression to a function of th instead of the current function of th/2,
using fr_ao2(expr, th)

   [wxMaxima: comment end   ] */


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


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


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


/* [wxMaxima: comment start ]
note that Mfi is real, so |Mfi|^2 = Mfi*Mfi = Mfi^2
   [wxMaxima: comment end   ] */


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


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


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


/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: input   start ] */
up1:UU(E,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);
/* [wxMaxima: input   end   ] */


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


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


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


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


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


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


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


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


/* [wxMaxima: comment start ]
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
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
General Helicity Cases
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
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))$
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: comment start ]
This agrees with Mfi worked out above for the RR --> RR case.

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


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


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


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


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

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

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

Using the function 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.
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: comment start ]
We see that none of the amplitudes is zero.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
save this as mssq_DS for future use.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
define mssq_ts as an expression involving all cosines
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
Let's compare this to MfiSQ_th_E from Sect. 1.3.1
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
is (equal (mssq_ts, MfiSQ_th_E));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
so this method produces an equivalent result.
   [wxMaxima: comment end   ] */


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


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

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

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

In both the nc_tr and m_tr methods, one must use a spin projection operator
 symbolized by S(σ, Sp ) (for the massive lepton or massive anti-lepton case).
 
[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. ]


   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
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.

   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
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) )$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
If we use the symbolic nc_tr(...) method,  we should make use of an invarR list of
4-momentum product replacements.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
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)$
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: input   start ] */
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);
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: input   start ] */
MSQn :  psq_to_Esq(MSQn,p,E,m), expand;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
MSQn :  psq_to_Esq(MSQn ,k,E,M), expand;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
this agrees with  the RR --> RR case
computed in Sec. 2.1 using explicit Dirac spinors.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
general amplitude squared for specified helicities using nc_tr method,
constructed using our experience simplifying the RR --> RR case above.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
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)$
/* [wxMaxima: input   end   ] */


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


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


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


/* [wxMaxima: comment start ]
The m_tr method is more than eight times faster at this job: see below.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
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
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
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);
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: input   start ] */
psq_to_Esq (%,p,E,m),  expand;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
psq_to_Esq (%,k,E,M),  expand;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
this agrees with  the reduced amplitude for the RR --> RR case
computed in Sec. 2.1 using explicit Dirac spinors.
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: input   start ] */
mesq_DS : AE_me(1,1,1,-1)^2;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
ratsimp(mesq_mtr);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
general amplitude squared for specified helicities using m_tr method,
constructed using our experience simplifying the RR --> RR case above.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
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)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
check its operation for RR --> RR case
   [wxMaxima: comment end   ] */


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


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


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


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


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


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


/* [wxMaxima: comment start ]
print out all reduced squared amplitudes and comput the sum using m_tr method
   [wxMaxima: comment end   ] */


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


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


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


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


/* [wxMaxima: input   start ] */
expand(mssq);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
compare this with mssq_DS computed above using explicit Dirac spinors
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
is (equal (%, mssq_DS));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
which shows that the m_tr method agrees with the explcit Dirac spinor method.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: subsubsect start ]
Using mat_trace to get the square of the RR --> RR amplitude
   [wxMaxima: subsubsect end   ] */


/* [wxMaxima: input   start ] */

/*  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);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
psq_to_Esq (%,p,E,m),  expand;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
psq_to_Esq (%,k,E,M),  expand;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
The following section is purely for pedagogical reasons and may be skipped.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
P(s, Sp);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
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) )


   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
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)


   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
 Example: RR --> RR case; we first repeat our definition of the relevant Dirac spinors 
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
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)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
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)  ):
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
up1b : sbar(up1)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
uub : up1 . up1b;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
replace the products sqrt(E-p)*sqrt(E+p)
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
uub : Ep_to_m(uub,p,E,m);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Pp1m :   ratsimp ( P(1, Sp1) . (sL(p1) + m*I4) );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
replace E^2 in favor of p^2
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Pp1m : Esq_to_psq(Pp1m,p,E,m);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
is(equal (uub, Pp1m));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: comment start ]
again we replace the products sqrt(E - k) * sqrt(E + k) in this matrix expression
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
vvb : Ep_to_m(vvb,k,E,M);
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: input   start ] */
vvb : ratsimp(vvb);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
Pp4M :   ratsimp (expand( P(1, Sp4) . (sL(p4) - M*I4) ));
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
replace E^2 in favor of k^2 and M^2
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Pp4M : Esq_to_psq(Pp4M,k,E,M);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
put matrix elements over common denominators:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Pp4M : ratsimp(Pp4M);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
is(equal (vvb,Pp4M));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */



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