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

    e(-,p1,σ1) + e(+,p2,σ2) --> γ(k1,λ1) + γ(k2,λ2)
   [wxMaxima: title   end   ] */


/* [wxMaxima: comment start ]
ee-gaga-HE.wxm,   Dec. 13, 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

"High Energy"  means center of momentum frame lepton
energy E  is  large compared to the electron mass m, so we
can neglect m (set m = 0).

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

In particular, the diagonal elements of the metric tensor g^{μ, ν} are
diag (g)  = (1, -1, -1, -1).
   [wxMaxima: comment end   ] */


/* [wxMaxima: section start ]
HE 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 an empty list [  ] as 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 section works out the unpolarized differential
  cross section in the high energy (HE) 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
for the process

      e(-,p1,σ1) + e(+,p2,σ2) --> γ(k1,λ1) + γ(k2,λ2).

These symbolic methods use our dirac3 package functions set_invarR, tr, Con, pullfac,
 comp_def, VP, sub_stu, to_ao2, fr_ao2

The amplitude for specific helicities σ1, σ2,  and specific photon polarizations
 λ1, λ2, used in Sec. 2 which employs explicit Dirac spinors and photon polarization
vectors is

M(σ1, σ2, λ1, λ2) = -e^2 Mr
 
and the "reduced amplitude"  Mr = M1 + M2, where
-------------------------------------------------------------------------
M1 =   - (v2bar SL( ε2c ) SL (p1 - k1) SL (ε1c) u1) / t  

M2 =   - (v2bar SL( ε1c ) SL (p1 - k2) SL (ε2c) u1) / u  
----------------------------------------------------------------------------
 where 

 t = (p1 - k1)^2,  
u = (p1 - k2)^2, 
s = (p1 + p2)^2
and  convervation of 4-momentum implies

p1^α + p2^α = k1^α + k2^α. 

In M1 and M2,

u1  = u (p1, σ1),  

v2bar = sbar ( v(p2, σ2)),

ε1c  =  the complex conjugate of the polarization 4-vector ε (k1, λ1),

ε2c  =  the complex conjugate of the polarization 4-vector ε (k2, λ2),

and SL(A) = "Feynman slash of 4-vector A" =  γ_α A^α (summed over the
dummy Lorentz index α for α = 0, 1, 2, 3.),  and SL is a linear operator:

 SL (A + B)  = SL (A)  +  SL (B)

and "sbar" is a Dirac3 package matrix operator with the effect
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
sbar (matrix([a1 + %i*a2],[b1 + %i*b2],[c1 + %i*c2],[d1 + %i*d2]));
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
In terms of ordinary Dirac matrix and spinor notation,

sbar (v2) = v2^{+} γ^0,  where v2^{+} = Hermitian conjugate of the spinor v2
    and  γ^0 is one of the four Dirac gamma matrices γ^α,  α = 0, 1, 2, 3.

 sbar (v2) is a "row vector," while v2 is a "column vector."

   [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  convervation of 4-momentum
    p1^α + p2^α = k1^α + k2^α
and with A^2 (with A being defined as a 4-vector)  standing for A_α A ^α summed over
the dummy Lorentz index  α for α = 0, 1, 2, 3.

    s = (p1 + p2)^2 = (k1 + k2)^2,
     t = (p1 - k1)^2 = (p2 - k2)^2, and  
     u = (p1 - k2)^2 = (p2 - k1)^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.

In this high-energy limit, E >>  m so effectively   m--> 0 and
 D(p1,p1) and D(p2,p2) are  approx 0  instead of m^2, etc. The physical (external)
photons have zero mass.

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

Both D and Gm are symbols which are declared symmetric in the file dirac3.mac, 
with the line:

declare ( [D,Gm], symmetric )$, 

and during symbolic evaluations involving these symbols,
reversing a pair of arguments is assumed to refer to the
same quantity.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
is (equal (Gm(mu,nu),Gm(nu,mu)) );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
is (equal (D(p,q),D(q,p)) );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
set_invarR (D(p1,p1) = 0,
         D(p1, p2) = s/2,
         D(p1, k1) =  - t/2,
         D(p1, k2) =  - u/2,
         D(p2, p2) = 0,
         D(p2, k1) =  - u/2,
         D(p2, k2) =  - t/2,
         D(k1, k1) = 0,
         D(k1, k2) = s/2,
         D(k2, k2) = 0)$
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
Let's check recovering these definitions with an ev(expr, conditions & rules) in disguise  here:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
D(p1,p2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
D(p1,p2), invarR;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
D(p2,p1), invarR;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
so we see that the symmetric property declared for the symbol D allows us to
only include the value of D(p1,p2), but still get the same value for D(p2,p1).
Since D(p1,p2) is a symbol standing for the scalar 4-vector dot product of
p1 and p2,  we need D(p1,p2) = D(p2,p1).

   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]

   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
The sum over all helicities σ and photon polarizations λ  of the absolute value squared
of Mfi is:

<|Mfi|^2> =  Σ(σ1, σ2, λ1, λ2)  | Mfi(σ1, σ2, λ1, λ2) |^2

                 = e^4  Σ(σ1, σ2, λ1, λ2)  | Mr(σ1, σ2, λ1, λ2) |^2 
                 =  e^4 MfiSQ, in the notation used below.

The sum over helicities and photon polarizations of the absolute
value square of the  "reduced" amplitude can be written in the
three term form:

  MfiSQ  =  < |Mr|^2>  =   M11n / t^2  +  M22n / u^2 + 2*M12n / (t*u),

and we concentrate on getting these three quantities M11n, M22n, and M12n  first.

Note that contractions on the repeated Lorentz indices mu and nu inside
a single trace performed with our symbolic trace function tr are automatically
carried out, and we don't have to carry out an explicit contraction using Con.

   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
M11n :  tr (mu,p1 - k1,nu,p1,nu,p1 - k1, mu,p2) ;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M22n :  tr (nu,p1 - k2,mu,p1,mu,p1 - k2, nu,p2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M12n :  tr (mu, p1 - k1, nu, p1, mu, p1 - k2, nu, p2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Now simplify based on the identity (see Sec. 2.1 in pdf) 
s + t + u = the sum of the mass-squared of each particle = 0 here
since the photon is massless and we are  neglecting the mass
of the electron and positron.

   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
M12n : M12n, s = -t -u, expand;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
so ignoring the term ~ M12n, which is zero,
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
MfiSQ : M11n / t^2 + M22n / u^2;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
pullfac (expr, fac) is a Dirac3 package function to pull out fac from each 
term, once expr is expanded out.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
MfiSQ : pullfac (MfiSQ, 8);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Since the Mandelstam variables are Lorentz scalar invariants, we have
in this expression for MfiSQ an expression which can be used in any
chosen frame of reference.  In the following, we use a prefactor A which
applies specifically to the center of momentum frame.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
We have absorbed e^4 into a factor we call A, with e^2 = 4πα.

In  2 --> 2 scattering in CMS,  A = (α^2/(4 s)) (pf/pi) (see p. 90-91 in
Halzen and Martin)  but in the high energy limit we neglect electron and positron masses
 so  |p-vec|   approx = particle energy E .

E(cm) = 2 E, shared equally by electron and positron in CMS,
and then shared equally by the two photons with equal and
opposite 3-momenta in CMS, 

so pf = pi, and A =  α^2/(4 s)
  
 Averaging over initial helicities of the electron and positron means we need to divide A by 4
  to get the unpolarized differential cross section in the CMS and in the  high energy limit.
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: comment start ]
which agrees with Peter Renton, Electroweak Interactions, Sec. 4.4, p. 169, Eq. (4.91), 
and R. Gastmans and Tai Tsun Wu, The Ubiquitous Photon (the Helicity Method for
QED and QCD), Sec. 2.4, p. 19, Eq. (2.39).
   [wxMaxima: comment end   ] */


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


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

(in center of momentum frame); th stands for the emergent angle of
the k1 outgoing photon relative to the direction of the incident electron p1.

We replace s, t, and u with explicit functions s_th, t_th, u_th of  the
angle th.

Let p be  the magnitude of the 3-momentum of the initial electron
and positron.  Since we are neglecting the mass of the electrons, the value of p is
the same as the energy E of the electron and the positron, and of the energy
of  each of the final photons.  Then the magnitude of the  3-momentum of each outgoing
photon is also equal to E in the CMS.

   [wxMaxima: comment end   ] */


/* [wxMaxima: subsubsect start ]
Use comp_def(...) to define the 4-momenta components
   [wxMaxima: subsubsect end   ] */


/* [wxMaxima: comment start ]
These are the CMS 4-momentum components.
All the particle 3-momentum vectors lie in the z-x plane.

comp_def ( p1 (E1, px, py, pz), ...) for example assigns to the particle
with 4-momentum p1 the energy E1, and the 3-momentum components the
values px, py, and pz.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
assume ( E > 0, th >= 0, th <= %pi )$

comp_def ( p1( E,0,0,E),
            p2( E,0,0,-E),
            k1 (E,E*sin(th),0,E*cos(th)),
            k2 (E,-E*sin(th),0,-E*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(k2);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: subsubsect start ]
Use VP(...) to define s_th, t_th, and u_th
   [wxMaxima: subsubsect end   ] */


/* [wxMaxima: comment start ]
now define explicit frame dependent functions of th, the angle of 
the k1 emerging photon relative to the p1 incident electron direction.

VP(...) is our "4-vector product" function.  For example
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
VP(pa,pa);
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: comment start ]
t_th goes to zero as th --> 0 (the "forward direction") and
u_th goes to zero as the --> π (the "backward direction").

Due to our neglect of the electron mass m, the differential
cross section diverges in both the forward and backward
direction.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
t_th, th = 0;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
u_th, th = %pi;
/* [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 in expr
   [wxMaxima: subsubsect end   ] */


/* [wxMaxima: comment start ]
above, we had MfiSQ in terms of Mandelstam variables.
Let MfiSQ_th be the same quantity, but expressed in terms of
the angle th (of the k1 photon in the CMS) relative to the
direction of travel of the incident electron.

We let MfiSQ_th  be the label for a suitably massaged expression...
this takes some experimentation  with Maxima simplification functions
(as well as our package functions ts and pullfac ):
   [wxMaxima: comment end   ] */


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


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


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


/* [wxMaxima: input   start ] */
 MfiSQ_th :   pullfac (%, 16);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
which agrees with Peter Renton (see ref. above).

Look at the forward and backward direction limits:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
limit (dsigdo_HE, th, 0, plus);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
limit (dsigdo_HE, th, %pi, minus );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Use taylor for a power series expansion thru terms of order th^2
about the point th = 0:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
taylor (dsigdo_HE, th, 0, 2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
so the differential cross section diverges as ~ 1/th^2 near th = 0.
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
Unpolarized Cross Section  WITHOUT  Mandelstam Variables, using symbolic tr and Con
   [wxMaxima: subsect 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) = 0,         
         D(p2, p2) = 0,         
         D(k1, k1) = 0,         
         D(k2, k2) = 0)$
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
The sum over all helicities σ and photon polarizations λ  of the absolute value squared
of Mfi is:

<|Mfi|^2> =  Σ(σ1, σ2, λ1, λ2)  | Mfi(σ1, σ2, λ1, λ2) |^2

                 = e^4  Σ(σ1, σ2, λ1, λ2)  | Mr(σ1, σ2, λ1, λ2) |^2 
                 =  e^4 MfiSQ, in the notation used below.

The sum over helicities and photon polarizations of the absolute
value square of the  "reduced" amplitude can be written in the
three term form:

  MfiSQ  =  < |Mr|^2>  =   M11n / t^2  +  M22n / u^2 + 2*M12n / (t*u),

and we concentrate on getting the three quantities M11n, M22n, and M12n  first.

Note that contractions on the repeated Lorentz indices mu and nu inside
a single trace performed with our symbolic trace function tr are automatically
carried out, and we don't have to carry out an explicit contraction using Con.

We later form MfiSQ using our known values of   s_th, t_th, and u_th.

   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
M11n :  tr (mu,p1 - k1,nu,p1,nu,p1 - k1, mu,p2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M22n :  tr (nu,p1 - k2,mu,p1,mu,p1 - k2, nu,p2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M12n :  tr (mu, p1 - k1, nu, p1, mu, p1 - k2, nu, p2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
note we get D(pa,pb)'s  since tr(...) works here without knowledge of the invariants involving s, t, u.
   [wxMaxima: comment end   ] */


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


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


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


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


/* [wxMaxima: comment start ]
so, ignoring the zero contribution of M12n_nc, and using known t_th and u_th values,
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
MfiSQ_th_meth2 : M11n_nc / t_th^2 + M22n_nc / u_th^2;
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: input   start ] */
MfiSQ_th_meth2 : pullfac (%, 16);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
this is the same result as MfiSQ_th which we obtained in the first method
using tr(...) method via the Mandelstam variable route.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
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.
   [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(...) and mcon
   [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: input   start ] */
M11n :  nc_tr (mu,p1 - k1,nu,p1,nu,p1 - k1, mu,p2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M22n :  nc_tr (nu,p1 - k2,mu,p1,mu,p1 - k2, nu,p2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M12n :  nc_tr (mu, p1 - k1, nu, p1, mu, p1 - k2, nu, p2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
again, ignoring M12n term,
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
MfiSQ_th_meth3 : M11n / t_th^2 + M22n / u_th^2;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
MfiSQ_th_meth3 : pullfac (%, 16);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
which confirms the correctness of this method.
   [wxMaxima: comment end   ] */


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


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


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


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


/* [wxMaxima: comment start ]
As an example of the m_tr method,

m_tr ( mu, p1,  nu, k1 )  is automatically converted into the explicit matrix expression:

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

where mat_trace is defined in the linearalgebra package and is
automatically loaded when we first try to use mat_trace.

In this example, k1[0] = E, k1[1] = E*sin(th), k1[2] = 0, k1[3] = E*cos(th),  and
k1 is a "zero based array" holding the components of the
4-momentum vector k1. The components were defined
above using the comp_def function.

sL(k1) is the Feynman slash  matrix corresponding to 
the 4-momentum k1.

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 ] */
M11n :   mcon ( m_tr (mu,p1 - k1,nu,p1,nu,p1 - k1, mu,p2), mu, nu);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Since m_tr is a purely matrix method (in disguise), no use is made of the symbolic tr function,
so the contents of the list invarR do not matter, and no automatic contraction is carried out
on repeated indices.  Instead you should call mcon (expr, indices).

Here is an example of what happens with m_tr if you forget to call mcon:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
 m_tr (mu,p1 - k1,nu,p1,nu,p1 - k1, mu,p2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M22n :   mcon ( m_tr (nu,p1 - k2,mu,p1,mu,p1 - k2, nu,p2), mu, nu );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M12n :   mcon ( m_tr (mu, p1 - k1, nu, p1, mu, p1 - k2, nu, p2), mu, nu);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
MfiSQ_th_meth4 : M11n / t_th^2 + M22n / u_th^2;
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: input   start ] */
MfiSQ_th_meth4 :  pullfac( expand(%), 16);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
which confirms the correctness of the mcon ( m_tr,  indices) method.

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


/* [wxMaxima: section start ]
HELICITY AMPLITUDES USING EXPLICIT DIRAC SPINORS AND POLARIZATION VECTORS
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
Using explicit Dirac matrices, in this section we use the package
    functions UU, VV, sbar,  Avsq, ts.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
We use explicit Dirac spinors, photon polarization vectors and Dirac
matrices to compute the amplitudes for the process

      e(-,p1,σ1) + e(+,p2,σ2) --> γ(k1,λ1) + γ(k2,λ2).


The amplitude for specific helicities σ1, σ2,  and specific photon polarizations
 λ1, λ2, employs explicit Dirac spinors and photon polarization
vectors is

M(σ1, σ2, λ1, λ2) = -e^2 Mr
 
and the "reduced amplitude"  Mr = M1 + M2, where
-------------------------------------------------------------------------
M1 =   - (v2bar SL( ε2c ) SL (p1 - k1) SL (ε1c) u1) / t  

M2 =   - (v2bar SL( ε1c ) SL (p1 - k2) SL (ε2c) u1) / u  
----------------------------------------------------------------------------
 where 

 t = (p1 - k1)^2,  
u = (p1 - k2)^2, 
s = (p1 + p2)^2
and  convervation of 4-momentum implies

p1^α + p2^α = k1^α + k2^α. 

In M1 and M2,

u1  = u (p1, σ1),  

v2bar = sbar ( v(p2, σ2)),

ε1c  =  the complex conjugate of the polarization 4-vector ε (k1, λ1),

ε2c  =  the complex conjugate of the polarization 4-vector ε (k2, λ2),

and SL(A) = "Feynman slash of 4-vector A" =  γ_α A^α (summed over the
dummy Lorentz index α for α = 0, 1, 2, 3.),  and SL is a linear operator:

 SL (A + B)  = SL (A)  +  SL (B)

and "sbar" is a Dirac3 package matrix operator with the effect

   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
To calculate M1 and M2, we need to define circular polarization 4-vectors
(actually we need the complex conjugates, indicated here by "c")  associated
with the two outgoing photons (4-momenta k1 and k2), which, in the center of momentum frame
have equal and opposite 3-momentum vectors k1_vec = - k2_vec.

Since we have two outgoing photons, calculating the amplitude requires
the complex conjugate of the polarization 4-vectors, called above ε1c and ε2c.

So ek1Rc is the complex conjugate of ek1R = right-handed circular polarization
4-vector associated with the photon with 4-momentum k1,  so ek1Rc corresponds
to outgoing photon k1 having helicity = +1.

Likewise ek1Lc corresponds to the outgoing photon k1 having helicity = -1,
(left-handed circular polarization).

The explicit circular polarization 4-vectors (used here) are derived in detail in the
worksheet photon2.wxm.

   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
comp_def (ek1Rc (0, - cos (th)/sqrt(2), %i/sqrt(2), sin(th)/sqrt(2) ),
                  ek1Lc (0, cos (th)/sqrt(2), %i/sqrt(2), - sin(th)/sqrt(2)),
                  ek2Rc (0,  cos (th)/sqrt(2),  %i/sqrt(2), - sin(th)/sqrt(2) ),
                  ek2Lc (0, - cos (th)/sqrt(2),  %i/sqrt(2),  sin(th)/sqrt(2)) )$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Define two functions Ek1c( lam ) and Ek2c (lam), where lam is short for "lambda",
which takes on values 1 and -1.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Ek1c(lam) :=  if lam = 1 then ek1Rc 
                else if lam = -1 then ek1Lc
              else (print ("error"), done) $
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
Ek2c(lam) :=  if lam = 1 then ek2Rc 
                else if lam = -1 then ek2Lc
              else (print ("error"), done) $
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: comment start ]
Then we can get the same correct Feynman slash matrix either of two ways:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
sL1 : sL (ek1Rc);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
sL2 : sL ( Ek1c (1) );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
is (equal (sL1, sL2));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
We expect to get 0 for the helicity amplitude RR --> RR

for e(-, p1, R) + e (+,p2,R) --> γ(k1, R) + γ(k2, R)
We have σ1 = 1, σ2 = 1,  λ1 = 1, λ2 = 1

define needed spinors  up1 and vp2b (recall definitions of
UU and VV above):
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
up1 : UU (E,E,0, 0, 1);
vp2b : sbar ( VV (E,E,%pi, 0, 1));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
 set M1 = M1n/t  and M2 = M2n/u,  where we can use t_th and u_th from our definitions above.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
M1n : - vp2b . sL(Ek2c(1)) . sL(p1 - k1) . sL(Ek1c(1)) . up1;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M2n : - vp2b . sL(Ek1c(1)) . sL(p1 - k2) . sL(Ek2c(1)) . up1;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
so the helicity amplitude RR --> RR is 0, as we expected.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Next consider the case RL --> RL,  up1 is same but need to redefine vp2b

for e(-, p1, R) + e (+,p2,L) --> γ(k1, R) + γ(k2, L)

We have σ1 = 1, σ2 = -1,  λ1 = 1, λ2 = -1
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
vp2b : sbar ( VV (E,E,%pi, 0, -1));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M1n : -  trigsimp ( vp2b . sL(Ek2c(-1)) . sL(p1 - k1) . sL(Ek1c(1)) . up1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M1 : M1n / t_th, factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M2n : - trigsimp (  vp2b . sL(Ek1c(1)) . sL(p1 - k2) . sL(Ek2c(-1)) . up1 );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M2 : M2n / u_th, factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
so the (reduced) helicity amplitude for R L --> R L is 
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Mr : trigsimp (M1 + M2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Now automate calculation of the reduced amplitudes for given helicity
assignments, defining the function Amplitude(s1v,s2v,lam1,lam2), which
stands for  Mr (σ1, σ2, λ1, λ2).

In this definition,  E is an unbound parameter, and  t_th and u_th have been  defined above.

This produces the reduced helicity amplitude as a function of θ for given
helicity assignments:  sv = +/-1 corresponds to R/L , likewise lam = +/- 1 
corresponds to R/L.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Amplitude (s1v,s2v,lam1,lam2) := 
  block ([up1,vp2b, M1, M2],
    up1 : UU (E, E, 0, 0, s1v),    
    vp2b : sbar ( VV (E, E, %pi, 0, s2v)),
    -  trigsimp ( vp2b . sL(Ek2c(lam2)) . sL(p1 - k1) . sL(Ek1c(lam1)) . up1),
   M1 : factor ( %% / t_th ),
   - trigsimp (  vp2b . sL(Ek1c(lam1)) . sL(p1 - k2) . sL(Ek2c(lam2)) . up1 ),
  M2 : factor (%% / u_th),
  trigsimp (M1 + M2))$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Test this function with the two cases we did "by hand."
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: section start ]
Helicity amplitudes --> unpolarized cross section
   [wxMaxima: section 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 pdf file) 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 )  
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
block([s1, s2, l1, l2, temp],    
    print ("  "),
    print (" s1 s2  l1  l2     amplitude     "), 
    for s1 in sL do
        for s2 in sL do
            for l1 in sL do
                for l2 in sL do (
                     temp :  Amplitude (s1, s2, l1, l2),                     
                     mssq : mssq + Avsq(temp),
                     print ("  "),
                      print (s1, s2, l1, l2," ",temp ) ),
    mssq : trigsimp (mssq),
    print ("   "))$
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
We will compare this with MfiSQ_th found in method 1 (First method using
Mandelstam variable route to unpolarized case).
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
trigsimp (mssq - MfiSQ_th);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
which confirms the equality, taking into account trig identities.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Let's factor out the "16" from the numerator of mssq.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
pullfac( num(mssq), 16) / denom (mssq);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: section start ]
Squared Helicity Amplitudes via mat_trace methods
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
We used comp_def above to define the 4-vector arrays ek1Rc, ek1Lc, ek2Rc, and ek2Lc,
corresponding the the complex conjugates, as well as the function Ek1c(lam)
and Ek2c(lam).

To allow computations with both real and complex polarization 4-vectors, we need
to add the definitions (pasted from photon2.wxm):
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
comp_def ( ek1R (0, - cos (th)/sqrt(2), - %i/sqrt(2), sin(th)/sqrt(2) ),
                  ek1L (0, cos (th)/sqrt(2), - %i/sqrt(2), - sin(th)/sqrt(2)),
                   ek2R (0,  cos (th)/sqrt(2), - %i/sqrt(2), - sin(th)/sqrt(2) ),
                  ek2L (0, - cos (th)/sqrt(2), - %i/sqrt(2),  sin(th)/sqrt(2)) )$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
Ek1(lam) :=  if lam = 1 then ek1R 
                else if lam = -1 then ek1L
              else (print ("error"), done) $
Ek2(lam) :=  if lam = 1 then ek2R 
                else if lam = -1 then ek2L
              else (print ("error"), done) $
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
For given helicity assignments, we write the square of the reduced amplitude 
  Mr = M1 + M2     as  (using  complex conjugate of a matrix expression is the
    same as the hermitian conjugate of the matrix exprression):

MfiSQ  =    |Mr|^2 = |M1|^2 + |M2|^2 + M2 * M1^{+} + M1 * M2^{+} 

              =  M11n/t^2  + M22n/u^2 +  M21n/(t*u)  + M12n/(t*u), where ^{+} indicates hermitian conjugate.

   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
case RR --> RR,  or more specifically

 e(-, p1, R) + e(+,p2,P) --> γ(k1,R)  +  γ(k2,R)

so we get factors such as P(1) . sL(p1)  from u1 . u1bar,   
P(-1) . sL(p2)  from v2 . v2bar,    sL (Ek1(1)), and  sL (Ek2(1))  for example.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
M11n : mat_trace (sL(Ek2c(1)) . sL(p1 - k1) . sL(Ek1c(1)) . P(1) . sL(p1) . sL(Ek1(1)) . 
        sL(p1 - k1) . sL(Ek2(1)) . P(-1) . sL(p2) )  ;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M22n : mat_trace (sL(Ek1c(1)) . sL(p1 - k2) . sL(Ek2c(1)) . P(1) . sL(p1) . sL(Ek2(1)) . 
        sL(p1 - k2) . sL(Ek1(1)) . P(-1) . sL(p2) );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M21n : mat_trace (sL(Ek1c(1)) . sL(p1 - k2) . sL(Ek2c(1)) . P(1) . sL(p1) . sL(Ek1(1)) . 
        sL(p1 - k1)  . sL(Ek2(1)) . P(-1) . sL(p2) );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M12n : mat_trace (sL(Ek2c(1)) . sL(p1 - k1) . sL(Ek1c(1)) . P(1) . sL(p1) . sL(Ek2(1)) . 
        sL(p1 - k2) . sL(Ek1(1)) . P(-1) . sL(p2) );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
So the square helicity amplitude for the case RR ---> RR is zero.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Case  RL --> RL,  more specifically

 e(-, p1, R) + e(+,p2,L) --> γ(k1,R)  +  γ(k2,L)

so we get factors such as P(1) . sL(p1)  from u1 . u1bar,
 P(1) . sL(p2)   from v2 . v2bar,   sL (Ek1(1)), and  sL (Ek2(-1))  for example.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
M11n :  trigsimp (   mat_trace (sL(Ek2c(-1)) . sL(p1 - k1) . sL(Ek1c(1)) . P(1) . sL(p1) . sL(Ek1(1)) . 
        sL(p1 - k1) . sL(Ek2(-1)) . P(1) . sL(p2) ) ) ;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M22n :  trigsimp (  mat_trace (sL(Ek1c(1)) . sL(p1 - k2) . sL(Ek2c(-1)) . P(1) . sL(p1) . sL(Ek2(-1)) . 
        sL(p1 - k2) . sL(Ek1(1)) . P(1) . sL(p2) )) ;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M21n :  trigsimp (  mat_trace (sL(Ek1c(1)) . sL(p1 - k2) . sL(Ek2c(-1)) . P(1) . sL(p1) . sL(Ek1(1)) . 
        sL(p1 - k1) . sL(Ek2(-1)) . P(1) . sL(p2) )) ;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M12n : trigsimp (  mat_trace (sL(Ek2c(-1)) . sL(p1 - k1) . sL(Ek1c(1)) . P(1) . sL(p1) . sL(Ek2(-1)) . 
        sL(p1 - k2) . sL(Ek1(1)) . P(1) . sL(p2) )) ;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
MfiSQ : trigsimp ( M11n / t_th^2 + M22n / u_th^2 + M21n / (t_th*u_th)  + M12n / (t_th*u_th) ) ;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Compare this result which used the mat_trace method with the square of the
helicity amplitude found using explicit Dirac spinors.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
trigsimp (MfiSQ - MfiSQ_DS);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
which confirms the equality with the square of the helicity amplitude computed
using the explicit Dirac spinor methods.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
For a short tutorial on the explicit spin projection matrix P(σ)  when used
in the high energy case (neglecting the mass of the particles), see section (3.4)
of the worksheet ee-mumu-HE.wxm.

   [wxMaxima: comment end   ] */


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



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