/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 17.10.0 ] */
/* [wxMaxima: title   start ]
differential scattering cross section
   for e(-) pi(+) --> e(-) pi(+)
in the High Energy Limit (HE)
   [wxMaxima: title   end   ] */


/* [wxMaxima: comment start ]
ePi-HE.wxm
Maxima by Example, Ch. 12, Dirac Algebra and Quantum Electrodynamics,
Dirac3 package, Ted Woollett, Feb. 12, 2019, 
http://web.csulb.edu/~woollett/
woollett@charter.net
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
First load our package:
   [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: section start ]
UNPOLARIZED CROSS SECTION USING tr and Con
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
e(-,p1,m) Π(+,p2,M) --> e(-,p3,m) Π(+,p4,M) cross section,

averaging over spin σ1 of the initial electron, and summing over the spin σ3 of the
final electron, 

dσ/dΩ (cms) = <|M(f,i)|^2> / (64 π^2 s), 
s = (p1+p2)^2, 
π = 3.141592653589793, where

<|M(f,i)|^2> = (1/2) Σ_(σ1) Σ_(σ3) |M(σ1, σ3)|^2 ,

and the invariant amplitude for given spin assignments is

M(σ1,σ3) =  (e^2/t) ubar(p3,σ3)γ^μ u(p1,σ1) (p2_μ + p4_μ), summed over μ,

where t = (p1 - p3)^2, and  σ = +1 for right-handed  helicity (R),  σ = -1 for left-handed  helicity (L)

   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
As shown in our references, (see, for example, Sec 6.3 in Aitchison and Hey, Gauge Theories
in Particle Physics,  2nd ed, 1989),

<|M(f,i)|^2> = (4 π α/ t )^2 L^{μ, ν} T_{μ, ν}, in terms of t = (p1 - p3)^2,

T^{μ, ν} = pion tensor = (p2 + p4)^μ (p2 + p4)^ν, and
L^{μ, ν} = electron tensor = (1/2) Trace ( (sL(p3) + m) γ^μ (sL(p1) + m) γ^ν ),

where sL(p) = Feynman slash matrix = p_α γ^α, summed over α (summation convention
for repeated indices in the same term).

In the high energy limit, which we consider here, we neglect the electron mass in
L^{μ, ν},  and we call this tensor L in our Maxima work.

In the high energy limit, which we consider here, we neglect the pion mass in
T^{μ, ν},  and we call this tensor T in our Maxima work.

We use UI(p,μ) to symbolically stand for p^μ (Upper Index component).

Our symbolic trace function is tr(...)

D(pa, pb) is the symbolic expression for the 4-vector product of pa and pb, defined
to be a symmetric function of its arguments.

Gm(mu,nu) is the symbolic metric tensor g^{μ, ν}.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
L :  expand ( tr (p3,mu,p1,nu) / 2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
T : ( UI(p2,mu) + UI(p4,mu) ) * ( UI(p2,nu) + UI(p4,nu) );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Con (expr, indices) carries out contraction on the indices supplied.
Con (expr)  carries out contraction on all repeated indices in expr.

For this example, the symbolic contraction function scon(expr, indices) is
automatically called.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
LT : Con (L * T, mu, nu);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Define components of 4-momenta in the center of momentum frame.

High energy limit, (HE), E approx equal to p = magnitude of 3-momentum
 p >> M >> m,  where M = pion mass, m = electron mass.

In CMS equal and opposite 3-momenta plus neglect of masses
means both electron and pion have same energy E both before
and after.

initial electron, 4-momentum p1, and 3-momentum approx. equal to E,
 moves along positive z axis,
initial pion, 4-momentum p2, and 3-momentum approx. equal to E, 
moves in direction of negative z axis,
final electron, 4-momentum p3, and 3-momentum approx. equal to E, 
moves in z-x plane, angle θ relative to z axis
final pion, 4-momentum p4, and 3-momentum approx. equal to E,
 moves in direction opposite to final electron

comp_def (pa (Ea, pax, pay, paz)) assigns frame dependent
energy Ea, and 3-momentum cartesian components pax, pay,
and paz to the 4-momentum pa, and creates an array pa[j]
such that pa[0] --> Ea, pa[1] --> pax, etc.
   [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),
             p3 (E,E*sin(th),0,E*cos(th)),
             p4 (E,-E*sin(th),0,-E*cos(th)) )$
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: comment start ]
Convert the general symbolic expression LT to a frame dependent expression
L_nc by using noncov(expr).
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
LT_nc : noncov (LT), factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
define frame dependent value t_th for the Mandelstam variable t, using
V(pa,pb) --> 4-vector dot product making use of frame dependent
components defined using comp_def above; i.e., uses the defined
arrays pa and pb with our metric tensor conventions.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
t_th : VP (p1-p3, p1-p3), factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Putting the pieces together, to get the approx. CMS differential scattering
cross section in the high energy limit:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
dsigdo_HE : (4*%pi*alpha)^2*LT_nc / t_th^2 / (64*%pi^2*s);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
As above we can calculate the frame dependent value s_th, and re-express
the differential scattering cross section in terms of E and θ.
We need to use the symbol theta in place of th to get the Greek
letter to be printed.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
s_th : VP(p1 + p2,p1+p2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
dsigdo_HE_th : dsigdo_HE, s = s_th;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: section start ]
UNPOLARIZED CROSS SECTION USING nc_tr and mcon
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
Using nc_tr for the electron tensor L,  instead of tr.

The function nc_tr first calls the symbolic trace function tr, and
then immediately calls noncov on each term to get a frame
dependent expression.  

In L_nc, we have an expression, for arbitrary mu and nu, involving
the frame dependent components and the array gmet[mu,nu] which
holds the numerical values of our metric tensor.

   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
L_nc :  expand ( nc_tr (p3,mu,p1,nu) / 2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
gmet[0,0];
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
To evaluate the needed contraction between the electron and pion tensors
we need to write the pion tensor in terms of the frame dependent components.
We call this T_nc.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
T_nc : ( p2[mu] + p4[mu] ) * ( p2[nu] + p4[nu] );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
To carry out the contraction of an expression written in terms of frame
dependent components, we use the syntax:

    mcon (expr, indices)

We must always supply the names of the requested contraction indices
when using mcon.

Call the contraction result LT_nc2 to distinguish the answer from
LT_nc found above.

   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
LT_nc2 : mcon (L_nc * T_nc, mu, nu);
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: input   start ] */
is (equal (LT_nc2, LT_nc));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: section start ]
Unpolarized Cross Section Using Mandelstam Variables Step
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
POPULATE THE LIST invarR OF symbolic D(pa,pb) 4-VEC DOT PRODUCT VALUES
expressed in terms of the Mandelstam variables  s, t, u using the set_invarR  function. 
We have convervation of 4-momenta:  p1 + p2 = p3 + p4, 
and the Mandelstam variables are defined by:

s = (p1+p2)^2 = (p3+p4)^2,
 t = (p1-p3)^2 = (p2-p4)^2, 
u = (p1-p4)^2 = (p2-p3)^2 

These definitions are independent of coordinate frame.

D(pa,pb)  is a symbolic representation of the 4-vector dot product (scalar product)
of 4-vectors pa and pb and is declared to be a symmetric function of its arguments.

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

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

This case of the command for set_invarR can be pasted into your
work from the package file mandelstam.txt.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
This is the resulting list which is used by the symbolic trace function tr(expr)
and ev_Ds(expr).
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
Using our expression for LT in step (%i4) above
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
we can replace all D's in terms of Mandelstam variables by using our
package function ev_Ds(expr).
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
LT_mv : ev_Ds(LT);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
We can then write down the differential scattering cross section in terms
of Mandelstam variables, but with the prefactors chosen to apply to the
center of momentum frame, as
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
dsigdo_HE_mv : (4*%pi*alpha)^2*LT_mv / t^2 / (64*%pi^2*s);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
To convert this expression into an explicit function of the electron scattering angle th and
the particle energy E,  we can use our previous expressions for s_th and t_th, and  also here
define u_th.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
u_th : VP(p1-p4,p1-p4), factor;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
We now use sub_stu (expr) to replace s by s_th, t by t_th, and u by u_th.
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: comment start ]
Without trying to factor the numerator and denominator of this expression,
we show that, after expanding the difference, we get zero, when comparing this
expression to our original explicit cross section above.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
expand (dsigdo_HE_mv_th - dsigdo_HE_th);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
which shows the equivalence.
Now lets try for a factored expression:
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
which reproduces the factored expression found above.
   [wxMaxima: comment end   ] */


/* [wxMaxima: section start ]
Electron Helicity Amplitudes for this problem
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
General comments about Dirac spinors and barred spinors for leptons.

   (See Sec. 2 of ee-mumu-AE.wxmx for corresponding comments about
   Dirac spinors for anti-leptons, which use VV(Ea,pa,θ,φ,σ).)

    If upa is a Dirac spinor (column vector) corresponding to 4-momentum pa, then
    sbar(upa) is the barred Dirac spinor (row vector) upa^{+} γ^0, 
    where (something)^{+} indicates Hermitian conjugate.

Both upa and sbar(upa) are  Maxima matrices.
    
    Our normalization is  upa_bar . upa = 2m.   (In this worksheet, we get zero, since we neglect
        the mass of the electron.)
    
    UU(Ea,pa,θ,φ,σ) is a generic Dirac spinor (4-element column vector) corresponding to a lepton with 
    relativistic energy  Ea = sqrt(pa^2 + ma^2), 3-momentum magnitude pa, 
    with direction of 3-momentum defined by the spherical polar angles
     θ and φ, and helicity σ (+1 for R, -1 for L).    

In the HE limit, E is approximately equal to p, ie., p >> m or M (m = electron mass, M = pion mass)
for both the electron and pion (center of momentum frame).    

   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
The invariant amplitude for given helicity assignments is

M(σ1,σ3) =  (e^2/t) ubar(p3,σ3)γ^μ u(p1,σ1) (p2_μ + p4_μ), summed over μ,

where t = (p1 - p3)^2, and  σ = +1 for right-handed  helicity (R),  σ = -1 for left-handed  helicity (L).

In the HE limit, E is approximately equal to p, ie., p >> m or M (m = electron mass, M = pion mass).

For contractions involving matrices, we need to use mcon (expr, indices) syntax.
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: subsect start ]
Case R --> R
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
Start with the case  R --> R:   e_R +  pi --> e_R + pi ; right-handed electron to right-handed elecron.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
up1 : UU(E,E,0,0,1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
up3 : UU(E,E,th,0,1);
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: comment start ]
We need to explore how to get to the simplest expression.
p2 and p4 will introduce cos(th) and sin(th).  
up1 and up3 will introduce cos(th/2) and sin(th/2).
First try, use t_th, which is a function of th.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
M_RR : (4*%pi*alpha / t_th) * mcon ( (sbar(up3) . Gam[mu] . up1) *  (p2[mu] + p4[mu]), mu);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Second try:  use t_th2, converting t_th by using to_ao2 (expr, a).
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
t_th2 : to_ao2( t_th, th);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M_RR : (4*%pi*alpha / t_th2) * mcon ( (sbar(up3) . Gam[mu] . up1) *  (p2[mu] + p4[mu]), mu);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
To cancel E^2  between the numerator and the denominator, we could either expand the expression followed by ratsimp
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
or use div_tb (expr, fac)  ( divide top and bottom by the same factor).
This is one of our package functions.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
div_tb(M_RR,E^2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Next try to convert f(th) to g(th/2) in this expression
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
So, this is a relatively simple expression for the amplitude M_RR.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Let's define a function Amplitude(s1v,s3v) to return the amplitude
corresponding to the values of s1v and s2v, respectively standing
for the values of the helicity σ1 = +/- 1 and σ2 = +/- 1.

+1 helicity corresponds to R (right-handed lepton)
-1 helicity corresponds to L (left-handed lepton)
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Amplitude (s1v,s3v) := 
  block ([up1,up3b],
    up1 : UU(E,E,0,0,s1v),
    up3b : sbar (UU(E,E,th,0,s3v)),
    (4*%pi*alpha / t_th2) * mcon ( (up3b . Gam[mu] . up1) *  (p2[mu] + p4[mu]), mu),
    ratsimp (expand (%%)),
     to_ao2 (%%, th))$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
We then find that M_LL = M_RR, and M_RL = 0, and M_LR = 0.
   [wxMaxima: comment end   ] */


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


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


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


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


/* [wxMaxima: comment start ]
dσ/dΩ (cms) = <|M(f,i)|^2> / (64 π^2 s), 

s = (p1+p2)^2, 

    and <|M(f,i)|^2> = (1/2) Σ_(σ1) Σ_(σ3) |M(σ1, σ3)|^2  = M_RR^2
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
dsigdo : Mfi_sqr / (64*%pi^2*s);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
%, s = s_th;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Let's compare the above value of dsigdo in terms of th/2 
with the expression dsigdo_HE_th in terms of th derived
above:
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: comment start ]
which shows they are equivalent forms.
   [wxMaxima: comment end   ] */



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