/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 17.10.0 ] */
/* [wxMaxima: title   start ]
compton-CMS-HE.wxm   
High Energy Compton Scattering in the
     Center of Momentum Frame 

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


/* [wxMaxima: comment start ]
compton-CMS-HE.wxm,   Dec. 12, 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

compton1.wxm evaluates Compton scattering in the rest frame
of the initial electron e(-, p1, σ1), whereas THIS WORKSHEET uses the
center of momentum frame, and also evaluates helicity amplitudes
and their squares.

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

Following Griffiths and Halzen/Martin, we interpret the Feynman rules
as giving an expression for - i M.
   [wxMaxima: comment end   ] */


/* [wxMaxima: section start ]
 HIGH ENERGY UNPOLARIZED DIFFERENTIAL SCATTERING CROSS SECTION 
   [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 section works out the unpolarized differential
  cross section in the high energy (HE) case,  in terms
  of the Mandelstam variables s, t, and u,  using symbolic methods
for the process

      e(-, p1, σ1) +  γ(k1, λ1)    -->    e(-, p2, σ2) + γ(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 (before we
take the HE limit and set m = 0)
-------------------------------------------------------------------------
M1 =   - (u2bar SL( ε1 ) (SL (p1 - k2) + m*I4) SL (ε2c) u1) / ( u - m^2)

M2 =   - (u2bar SL( ε2c ) (SL (p1 + k1) + m*I4) SL (ε1) u1) / (s - m^2)
----------------------------------------------------------------------------
 where 

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

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

In M1 and M2,

u1  = u (p1, σ1),  

u2bar = sbar ( u(p2, σ2)),

ε1  =   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 (u2) = u2^{+} γ^0,  where u2^{+} = Hermitian conjugate of the spinor u2
    and  γ^0 is one of the four Dirac gamma matrices γ^α,  α = 0, 1, 2, 3.

 sbar (u2) is a "row vector," while u2 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^α + k1^α = p2^α + 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 + k1)^2 = (p2 + k2)^2,
     t = (p1 - p2)^2 = (k1 - 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.

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: comment start ]
At this point we set m = 0 for the high energy limit E >> m
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
this is our list of replacement rules for D(pa,pb)
   [wxMaxima: comment 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 ]
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 is written as

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

and we concentrate on getting the four numerator quantities M11n, M22n, M21n, 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 - k2 , nu, p1 , nu, p1 - k2 , mu, p2 );
/* [wxMaxima: input   end   ] */


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


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


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


/* [wxMaxima: comment start ]
We then add these terms up, each divided by the appropriate denominator,

MfiSQ  =  < |Mr|^2>  =   M11n / u^2  +  M22n / s^2  + 
                                                   M21n / ( u*s )  + M12n / ( u*s )
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
MfiSQ : M11n / u^2  + M22n / s^2 + 
                    M21n / (u*s) + M12n / (u*s), ratsimp;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
We can use the identity relating the Mandelstam variables
  s + t + u = sum of mass-squared for each particle involved,
so here s + t + u = 2 m^2,  or s = 2 m^2 - t - u --> - t - u   in the high 
energy limit.
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: subsubsect start ]
CONVERSION TO A CENTER OF MOMENTUM FRAME FUNCTION OF THE SCATTERING ANGLE θ
   [wxMaxima: subsubsect end   ] */


/* [wxMaxima: comment start ]
  CONVERSION TO A CENTER OF MOMENTUM FRAME FUNCTION OF THE PHOTON SCATTERING ANGLE θ 

 th stands for the emergent angle of the k2 outgoing 
photon relative to the direction of the incident photon k1
(that is, relative to the positive z-axis).

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

Let k be  the energy and also the 3-momentum magnitude of the initial photon k1
with 4-momentum components [k,0,0,k].

Conservation of 4-momentum p1 + k1 = p2 + k2  breaks down into

1. conservation of 3-momentum:  In the center of momentum frame (CMS), the total
3-momentum is zero, so the 3-momentum of the initial electron is equal in magnitude
and opposite in direction to that of the initial photon. Since we here neglect the mass m
of the electron, the initial electron has 4-momentum components [k,0,0,-k].

2. The total initial energy Ei = 2 k = Ef  = kf + pf.  But in CMS pf = kf, so kf = k = pf.
In this frame, the magnitude of 3-momentum is unchanged, only the directions of
the 3-momenta are rotated by the scattering angle of the photon th.

   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
These define the 4-momenta in the iCMS.

All the particle 3-momentum vectors lie in the z-x plane.

In general, 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.

th, standing for θ, is the angle between the emergent photon k2's
3-momentum vector and the positive z-axis.

   [wxMaxima: comment end   ] */


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

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


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


/* [wxMaxima: input   start ] */
map ('listarray, [k1, p1, k2, p2] );
/* [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 k2 emerging photon relative to the k1 incident photon direction.

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


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


/* [wxMaxima: input   start ] */
grind(%)$
/* [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.

 s = (p1 + k1)^2 = (p2 + k2)^2,
     t = (p1 - p2)^2 = (k1 - k2)^2, and  
     u = (p1 - k2)^2 = (p2 - k1)^2 

   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
  s_th : factor(VP (p1 + k1, p1 + k1));
  
  t_th : factor (VP (p1 - p2, p1 - p2));
  
  u_th : factor (VP (p1 - k2, p1 - k2));
/* [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 have MfiSQ in terms of Mandelstam variables.
Let MfiSQ_th be the same quantity, but expressed in terms of
the angle th (of the k2 photon in the initial electron rest frame) relative to the
direction of travel of the incident photon k1.

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 ] */
MfiSQ_th : trigsimp (  sub_stu (MfiSQ)), factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
The form of the prefactor A appropriate to this process, which will
mulitply |Mr|^2 to determine the differential scattering cross section for
in the CMS, specific given spin states, is derived in various references. 

See Halzen/Martin, Sec. 4.3, or Mark Thomson, Sec. 3.5.1.

dsigdo(cms) = (1/ 64 π^2 s) (pf/pi) |Mfi|^2, where s = E(cms)^2 = (2 k)^2 = 4 k^2,
and here pf = pi = k.  We want to write dsigdo(cms) = A |Mr|^2 for specified
helicities, where Mr = "reduced amplitude" = -e^2 Mfi.  We also have e^2 = 4 π α.

So we get dsigdo(cms) = (α^2/ (4 s) |Mr|^2  = A |Mr|^2 for specified helicities.

Then our prefactor A = (α^2/ (4 s).

   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
To get the UNPOLARIZED differential cross section
we average over the two polarization states of the incident photon and
the two spin states of the initial electron, and sum over the polarization state
of the final photon and  the helicity of the final electron. This means that we must
divide the prefactor by (2) (2) = 4.

We then get the "unpolarized differential scattering cross section"
in the CMS  as
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
dsigdo : (A/4) * MfiSQ_th;
/* [wxMaxima: input   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 is written as

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

and we concentrate on getting the four numerator quantities M11n, M22n, M21n, 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 - k2 , nu, p1 , nu, p1 - k2 , mu, p2 );
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: input   start ] */
M12n :  tr (mu, p1 - k2 , nu, p1 , mu, p1 + k1 , 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 ] */
M21n_nc : noncov (M21n);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
so we can ignore the cross terms, which are zero.
 Including the denominators with the u --> u_th and s --> s_th, and using trigsimp,
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
MfiSQ_th_meth2 :  trigsimp (  M11n_nc / u_th^2  + M22n_nc / s_th^2 ), factor
;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Compare with method 1:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
expand (MfiSQ_th_meth2 - MfiSQ_th);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
which confirms that the method 2 answer is correct.
   [wxMaxima: comment 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 s_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.

In this example, the contractions on mu and nu are done automatically
when nc_tr calls the symbolic trace function tr, so we don't need to call
mcon for contraction.

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 - k2 , nu, p1 , nu, p1 - k2 , mu, p2 );
/* [wxMaxima: input   end   ] */


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


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


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


/* [wxMaxima: input   start ] */
MfiSQ_th_meth3 :  trigsimp (  M11n / u_th^2  + M22n / s_th^2), factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
expand (MfiSQ_th_meth3 - MfiSQ_th_meth2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
which confirms the correctness of  method 3 using nc_tr.
   [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 (k2);
/* [wxMaxima: input   end   ] */


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


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


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

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

mat_trace (  Gam[mu] .  ( sL(p1) + m*I4)   . Gam[nu]  . sL(k2)   )


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

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

sL(k2) is the Feynman slash  matrix corresponding to 
the 4-momentum k1, and I4 is the unit matrix in four dimensions.

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: comment start ]
Because m_tr does NOT call the symbolic trace function tr, we must explicitly
call mcon (expr, indices) to get the sum over the Lorentz indices (contraction).
   [wxMaxima: comment end   ] */


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


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


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


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


/* [wxMaxima: input   start ] */
MfiSQ_th_meth4 :  trigsimp (  M11n / u_th^2  + M22n / s_th^2), factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
expand (MfiSQ_th_meth4 - MfiSQ_th_meth3);
/* [wxMaxima: input   end   ] */


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

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 ]
HIGH ENERGY HELICITY AMPLITUDES USING EXPLICIT DIRAC SPINORS AND PHOTON POLARIZATION VECTORS
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
Using explicit Dirac matrices, in this section we use the package
    functions UU,  sbar,  Avsq, ts, and fr_ao2(expr, th)
   [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).
    
    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. 
    
    The parameters E, p, θ, φ for each electron should agree with those specified in the comp_def statement above.

   [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) +  γ(k1, λ1)    -->    e(-, p2, σ2) + γ(k2, λ2).

The amplitude for specific helicities σ1, σ2,  and specific photon polarizations
 λ1, λ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
( before we set m = 0 in the HE limit)
-------------------------------------------------------------------------
M1 =   - (u2bar SL( ε1 ) (SL (p1 - k2) + m*I4) SL (ε2c) u1) / ( u - m^2)

M2 =   - (u2bar SL( ε2c ) (SL (p1 + k1) + m*I4) SL (ε1) u1) / (s - m^2)
----------------------------------------------------------------------------
 where 

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

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

In M1 and M2,

u1  = u (p1, σ1),  

u2bar = sbar ( u(p2, σ2)),

ε1  =   the polarization 4-vector ε (k1, λ1),

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

SL indicates symbolically here the Feynman slash operator, which
is implemented in matrix notation by sL(a).
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Let ek1R and ek1L represent the polarization 4-vectors ε (k1, λ1 = 1), ε (k1, λ1 = -1), respectively.
 and then let ek1Rc and ek1Lc represent the complex conjugates.
 
Likewise, let ek2R and ek2L represent the polarization 4-vectors ε (k2, λ2 = 1), ε (k2, λ2 = -1), respectively.
 and then let ek2Rc and ek2Lc represent the complex conjugates.

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

In this worksheet, photon k1 is moving along the positive z-axis. In photon2.wxm,
Sec. 1, we derive circular polarization 4-vectors for a photon with k-hat = [0,0,1].
(In Sec. 1 of photon2.wxm, we used notation  epsR, epsL, which we alter
here to call ek1R and ek1L.)

In Sec. 3 of photon2.wxm, we derive transverse circular polarization 4-vectors
for a photon with k=hat = [sin(th), 0, cos(th)], appropriate for photon k2 here.
( in Sec. 3 of photon2.wxm, we called the photon's 4-momentum k1; we alter
  our definitions here to reflect that we are calling our outgoing photon k2.)

   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
comp_def ( ek1R (0, -1/sqrt(2),  -%i/sqrt(2), 0 ),
                    ek1L (0,  1/sqrt(2),  -%i/sqrt(2),  0),
                    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: input   start ] */
map ('listarray, [ek1R,ek1L,ek2Rc,ek2Lc] );
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
Ek1(lam) :=  if lam = 1 then ek1R 
                else if lam = -1 then ek1L
              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: comment start ]
Then we can get the same correct Feynman slash matrix in either of two ways:
   [wxMaxima: comment end   ] */


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


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


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


/* [wxMaxima: comment start ]
edited to here


   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
Amplitude for Case RR --> RR Using Explicit Dirac Spinors
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
As a reminder, we repeat that

M(σ1, σ2, λ1, λ2) = -e^2 Mr
 
and the "reduced amplitude"  Mr = M1 + M2, where, in the high energy
limit (m = 0) 
-------------------------------------------------------------------------
M1 =   - (u2bar SL( ε1 ) SL (p1 - k2)  SL (ε2c) u1) / u

M2 =   - (u2bar SL( ε2c ) SL (p1 + k1)  SL (ε1) u1) / s
----------------------------------------------------------------------------
 where 

 t = (p1 - p2)^2,  
u = (p1 - k2)^2, 
s = (p1 + k1)^2

and

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.

   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
We also remind ourselves of the components declared
for the 4-momentum vectors of the particles above.

comp_def (            
            k1 (k, 0, 0, k),
            p1( k, 0, 0, - k),
            k2 (k, k*sin(th), 0, k*cos(th)),
            p2 (k, - k*sin(th), 0, - k*cos(th)))$

The spherical polar angles appropriate to p1 are th1 = %pi , and phi1 = 0
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
[sin(th1)*cos(phi1), sin(th1)*sin(phi1), cos(th1)], [th1 = %pi, phi1 = 0];
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
The spherical polar angles appropriate to p2 are th2 = %pi - th, and phi2 = %pi
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
[sin(th2)*cos(phi2), sin(th2)*sin(phi2), cos(th2)], [th2 = %pi - th, phi2 = %pi];
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
We first consider the helicity amplitude RR --> RR
which means the process:

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

Then we have σ1 = 1, σ2 = 1,  λ1 = 1, λ2 = 1

Define the needed spinors  up1 and up2b (recall  the definition of
UU above)  where E^2 = m^2 + p^2,  and we work in the center of
momentum frame (CMS).


   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
up1 : UU (k, k, %pi, 0, 1);
up2b : sbar ( UU (k, k, %pi - th, %pi, 1));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
 set M1 = M1n/ u  and M2 = M2n/ s,  where we can use  u_th and s_th from our definitions above,
and the numerators are:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
M1n : -  trigsimp ( up2b . sL(Ek1(1)) . sL(p1 - k2) . sL(Ek2c(1)) . up1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M2n : -  trigsimp ( up2b . sL(Ek2c(1)) . sL(p1 + k1) . sL(Ek1(1)) . up1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Now get the helicity amplitude M_RR_RR, including the angle
dependent denominators:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
M_RR_RR : trigsimp ( M1n/ u_th + M2n / s_th );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M_RR_RR : trigsimp (fr_ao2 (expand (%), th)),factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: subsect start ]
|Amplitude|^2 for Case RR --> RR Using mat_trace
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
To use mat_trace to compute |Mr|^2 for specified helicities, we need to
add definitions of the photon polarization 4-vectors ek1Rc, ek1Lc (the
complex conjugates of ek1R and ek1L) and ek2R and ek2L,  and then
define the functions Ek1c(λ) and  Ek2(λ).
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
comp_def ( ek1Rc (0, -1/sqrt(2),  %i/sqrt(2), 0 ),
                    ek1Lc (0,  1/sqrt(2),  %i/sqrt(2),  0),
                    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 ] */
map ('listarray, [ek1Rc, ek1Lc, ek2R, ek2L] );
/* [wxMaxima: input   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 ] */
Ek2(lam) :=  if lam = 1 then ek2R 
                else if lam = -1 then ek2L
              else (print ("error"), done) $
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
We then calculate, for specified helicities, MfiSQ = |Mr|^2 = |M1 + M2\^2 
 =  M1*conjugate(M1) + M2*conjugate(M2) + M2*conjugate(M1) + M1*conjugate(M2)  = 
M11n / u^2 + M22n/ s^2 + M21n/ (u*s) + M12n / (u*s).

For m = 0 we can use the matrix identities:  u1 . u1b = P(σ1) . sL(p1)  and u2 . u2b = P(σ2) . sL(p2)

P(sv) is the lepton spin projection matrix for m = 0, where sv = +/- 1
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: input   start ] */
P(-1);
/* [wxMaxima: input   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.wxmx.

For the case RR --> RR, we have σ1 = 1, σ2 = 1 for the electrons and
λ1 = 1,  λ2 = 1  for the photons.
   [wxMaxima: comment end   ] */


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


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


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


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


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


/* [wxMaxima: subsect start ]
Comparison of Square of Dirac spinor amplitude with mat_trace method for RR --> RR
   [wxMaxima: subsect end   ] */


/* [wxMaxima: subsubsect start ]
Compare M11n with M11n_DS  = M1n^2
   [wxMaxima: subsubsect end   ] */


/* [wxMaxima: comment start ]
compare M11n_DS  = M1n^2 computed using Dirac spinors above with M11n computed using mat_trace.
Note that M1n is real.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
M11n_DS : M1n^2;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
M11n, factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
which shows the equivalence of M11n_DS with M11n.
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsubsect start ]
Compare M22n with M22n_DS  = M2n^2
   [wxMaxima: subsubsect end   ] */


/* [wxMaxima: comment start ]
compare  M22n_DS = M2n^2 computed above using Dirac spinors with
M22n computed with mat_trace. Note that M2n is real.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
M22n_DS : trigsimp (M2n^2);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
M22n, factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
which shows the equivalence of M22n with M2n^2
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsubsect start ]
Compare M21n with M21n_DS  = M2n * M1n
   [wxMaxima: subsubsect end   ] */


/* [wxMaxima: comment start ]
Compare M21n computed with mat_trace with M21n_DS = M2n*M1n where M2n and M1n
were computed above using explicit Dirac spinors, and M1n and M2n are both real.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
M2n;
M1n;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M21n_DS : expand ( M2n * M1n );
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: comment start ]
which shows the equivalence of M21n and M2n * M1n
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsubsect start ]
Can We Massage M1n and M2n into a More Convenient Form?
   [wxMaxima: subsubsect end   ] */


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


/* [wxMaxima: input   start ] */
M1n_alt : trigsimp ( fr_ao2 (M1n, th)), factor;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
M2n_alt : trigsimp ( fr_ao2 (M2n, th)), factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M1_alt : trigsimp (M1n_alt / u_th), factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M2_alt : trigsimp (M2n_alt / s_th), factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: subsect start ]
Construct  |M_RR_RR|^2 from Dirac Spinor Amplitude for RR --> RR
   [wxMaxima: subsect end   ] */


/* [wxMaxima: input   start ] */
M11n_DS : M1n*conjugate(M1n);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
M11_DS : trigsimp (M11n_DS / u_th^2), factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M22n_DS : M2n*conjugate(M2n);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M22n_DS : trigsimp( fr_ao2(expand(%),th) ), factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M22_DS : trigsimp (M22n_DS / s_th^2) , factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M21n_DS : M2n*conjugate(M1n);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M21n_DS : trigsimp ( fr_ao2(expand(%),th ) ) , factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M21_DS : trigsimp (M21n_DS / (u_th * s_th) ) , factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M12n_DS : M1n*conjugate(M2n);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M12n_DS : trigsimp ( fr_ao2 (expand(%), th) ) ,factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M12_DS : trigsimp (M12n_DS / (u_th * s_th)) , factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
MSQ_DS : trigsimp (M11_DS + M22_DS + M21_DS + M12_DS) , factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
which agrees with our mat_trace result above, MfiSQ. for case RR --> RR
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: subsect start ]
AmplitudeSQ Function using Dirac Spinor Method
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
construct Dirac spinor amplitude squared |M_DS|^2  for given helicities

We can use the ideas that worked in the last subsubsection.

AmplitudeSQ (σ1, λ1, σ2, λ2) for p1 + k1 --> p2 + k2

s1v and s2v are respectively the values of σ1 and σ1 for the electrons p1 and p2.
lam1 and lam2 are respectively the values of λ1 and λ2 for the photons k1 and k2.

In this function, u_th and s_th are globally known expressions.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
AmplitudeSQ (s1v, lam1, s2v, lam2) := 
  block ([up1,up2b, M1n, M2n, M11n, M22n,M21n,M12n,
                    M11,M22,M21,M12],
    up1 : UU (k, k, %pi, 0, s1v),
    up2b : sbar ( UU (k, k, %pi - th, %pi, s2v)),
    
    M1n : -  trigsimp ( up2b . sL(Ek1(lam1)) . sL(p1 - k2) . sL(Ek2c(lam2)) . up1),
    M11n : M1n*conjugate(M1n),
    M11n :  factor (trigsimp ( fr_ao2( expand (%%),th))),
    M11  :  factor ( trigsimp (M11n / u_th^2)),
        
    M2n : -  trigsimp ( up2b . sL(Ek2c(lam2)) . sL(p1 + k1) . sL(Ek1(lam1)) . up1),
    M22n : M2n*conjugate(M2n),
    M22n :  factor (  trigsimp( fr_ao2(expand(%%),th) )),
    M22 :   factor (  trigsimp (M22n / s_th^2)) ,
    
    M21n : M2n*conjugate(M1n),
    M21n :  factor ( trigsimp ( fr_ao2(expand(%%),th ) )) ,
    M21 :  factor ( trigsimp (M21n / (u_th * s_th) )),

   M12n : M1n*conjugate(M2n),
    M12n :  factor ( trigsimp ( fr_ao2(expand(%%),th ) )) ,
    M12 :  factor ( trigsimp (M12n / (u_th * s_th) )),    
     factor ( trigsimp (M11 + M22 + M21 + M12)))$

/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Look at "non-flip" helicity cases.

Let's try AmplitudeSQ for the case RR --> RR
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
 LR --> LR
   [wxMaxima: comment end   ] */


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


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


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


/* [wxMaxima: comment start ]
LL  --> LL
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
 An example of a "photon  helicity flip" squared amplitude:
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
An example of a "electron  helicity flip" squared amplitude:
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: subsubsect start ]
Polarized Amplitudes to Unpolarized Cross Section
   [wxMaxima: subsubsect end   ] */


/* [wxMaxima: comment start ]
We now compute  MfiSQ_th =   Σ(σ1, λ1, σ2, λ2)  | Mr(σ1, λ1, σ2, λ2) |^2  
                                                     =    Σ(σ1, λ1, σ2, λ2)  AmplitudeSQ (σ1, λ1, σ2, λ2)
and call the result mssq, to avoid overlap with our previous definition of MfiSQ_th above
(using the symbolic tr method with Mandelstam variables).

We start by setting the global variable mssq to 0, and defining a list sL.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
Now we acculumate the sum the the reduced Dirac spinor amplitudes squared
(temp) for specified helicities, and the sum of the squares is mssq.
   [wxMaxima: comment end   ] */


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


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


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


/* [wxMaxima: comment start ]
After some experimentation, we get the same form as MfiSQ_th:
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
We found above, using four different methods,  the expression MfiSQ_th
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
so we have agreement.
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
AmplitudeSQ Function using the mat_trace method
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
For comparison purposes, we construct a mat_trace version of the square
of amplitudes for specified helicities.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
AmplitudeSQ_MT (s1v, lam1, s2v, lam2) := 
   block ( [M11n, M22n, M21n, M12n  ],
       M11n :  trigsimp (   mat_trace (sL(Ek1(lam1)) . sL(p1 - k2) . sL(Ek2c(lam2)) . P(s1v) . sL(p1) . 
            sL(Ek2(lam2)) . sL(p1 - k2) . sL(Ek1c(lam1)) . P(s2v) . sL(p2) ) ),
       M22n :  trigsimp (  mat_trace (sL(Ek2c(lam2)) . sL(p1 + k1) . sL(Ek1(lam1)) . P(s1v) . sL(p1) . 
            sL(Ek1c(lam1)) . sL(p1 + k1) . sL(Ek2(lam2)) . P(s2v) . sL(p2) )),
       M21n :  trigsimp (  mat_trace (sL(Ek2c(lam2)) . sL(p1 + k1) . sL(Ek1(lam1)) . P(s1v) . sL(p1) . 
            sL(Ek2(lam2)) .  sL(p1 - k2) . sL(Ek1c(lam1)) . P(s2v) . sL(p2) )),
       M12n : trigsimp (  mat_trace (sL(Ek1(lam1)) . sL(p1 - k2) . sL(Ek2c(lam2)) . P(s1v) . sL(p1) . 
            sL(Ek1c(lam1)) . sL(p1 + k1) . sL(Ek2(lam2)) . P(s2v) . sL(p2) )),
     factor (trigsimp ( M11n / u_th^2 + M22n / s_th^2 + M21n / (s_th*u_th)  + M12n / (s_th*u_th) )))$
    
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Try out this function for non-zero amplitude cases.
First case RR --> RR
   [wxMaxima: comment end   ] */


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


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


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


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


/* [wxMaxima: comment start ]
and for a helicity flip case:
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: subsubsect start ]
Polarized Amplitudes to Unpolarized Cross Section
   [wxMaxima: subsubsect end   ] */


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


/* [wxMaxima: input   start ] */
block([s1, l1, s2, l2, temp],    
    print (" mat_trace method "),
    print ("    "),
    print (" s1  l1  s2  l2     |amplitude|^2     "), 
    for s1 in sL do
        for l1 in sL do
            for s2 in sL do
                for l2 in sL do ( 
                     temp :  AmplitudeSQ_MT (s1, l1, s2, l2),
                     mssq : mssq + temp,
                     print ("  "),
                      print (s1, l1, s2, l2," ",temp ) ),
    mssq :  factor ( trigsimp (mssq)),
    print ("   "))$
/* [wxMaxima: input   end   ] */


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


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


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


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


/* [wxMaxima: comment start ]
so we have agreement.

The mat_trace method takes a little longer than the explicit Dirac spinor method,
but takes less coding time, with less simplification steps needed than the
Dirac spinor method.

   [wxMaxima: comment end   ] */



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