/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 17.10.0 ] */
/* [wxMaxima: title   start ]
e-mu-scatt.wxm 

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


/* [wxMaxima: comment start ]
e-mu-scatt.wxm,  Oct. 17, 2018
 Edwin (Ted) Woollett,   Maxima by Example Ch. 12,
Dirac Algebra and Quantum Electrodynamics, ver. 3, "Dirac3"
http://web.csulb.edu/~woollett/

 We use Heaviside-Lorentz electromagnetic units, e^2 = 4πα
In general, we follow the conventions in Peskin & Schroeder: Quantum Field Theory.
   [wxMaxima: comment end   ] */


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


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


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


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

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

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


/* [wxMaxima: comment start ]
This first section of e-mu-scatt.wxm works out the unpolarized differential scattering
  cross section for arbitrary energy, first in an arbitrary frame in terms
  of the Mandelstam variables s, t, and u, and then  in the CMS frame, using symbolic methods.

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

The amplitude for specific helicities σ1, σ2, σ3, σ4, used in Sec. 2 which employs
explicit Dirac spinors, is

M(σ1, σ2, σ3, σ4) = -e^2 Mr 
 where the "reduced amplitude" Mr = (u3bar γ^μ u1) ( u4bar γ_μ u2)/ t,
summed over the Lorentz index μ , where  t = (p1-p3)^2 = (p2-p4)^2, 
We call Mr the "reduced" amplitude.
   [wxMaxima: comment end   ] */


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


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

The electron (p1 --> p3) has mass m,
the muon (p2 --> p4) has mass M.

D(pa,pb)  is a symbolic representation of the 4-vector dot product (scalar product)
of 4-vectors pa and pb.

The dimensions of s, t and u are all the same
[s]  = [t]  = [u]  = E^2
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
In lowest order perturbation theory, only one
diagram contributes (exchange of a photon
between electron and muon).

omit the factor -e^2 from Mfi for now, leaving Mfi = Mr.

The sum over all helicities of the absolute value squared
of Mr is:
<|Mfi|^2> =  Σ(σ1, σ2, σ3, σ4)  | Mr(σ1, σ2, σ3, σ4) |^2
which we write as
 <|Mfi|^2> = Mn/t^2 , and Mn is the contraction of the product
of two traces.

The arg p3+m inside our  tr(...)   function stands for the square matrix: 

sL(p3) + m*I4  =  Feynman-slash(p3) + m*I4.

The arg mu inside tr function stands for the square Dirac matrix γ^μ.

We let MfiSQ stand for  <|Mfi|^2>
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: subsubsect start ]
Ignore the electron mass for now
   [wxMaxima: subsubsect end   ] */


/* [wxMaxima: comment start ]
We consider this result in the limit that we can neglect the small
electron mass m and concentrate on the numerator.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Mn : Mn, m = 0;
/* [wxMaxima: input   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 scattering angle.
Replace s, t, and u with explicit functions s_th, t_th, u_th of th.

Let k be  the magnitude of the 3-momentum of the initial electron and
initial muon and kp be the magnitude of the final electron and final muon in CMS.
Conservation of energy implies
   sqrt(m^2 + k^2) + sqrt(M^2 + k^2) = sqrt(m^2 + kp^2) + sqrt(M^2 + kp^2)
This implies kp = k, so pf = pi and A = α^2/( 4 s).

If we neglect the mass of the electron, then Ee --> k

Let the energy of the muon be E,
 Then  E^2 = k^2 + M^2, and

Then    E(CM) =  E + k, and  s = E(CM)^2  =  (E + k)^2
   [wxMaxima: comment end   ] */


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


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

comp_def ( p1( k,0,0,k),          
             p2( E,0,0,-k),             
             p3 (k,k*sin(th),0,k*cos(th)),             
             p4 (E,-k*sin(th),0,-k*cos(th)))$
/* [wxMaxima: input   end   ] */


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


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


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


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


/* [wxMaxima: comment start ]
now define explicit frame dependent functions of th = scattering angle
   [wxMaxima: comment end   ] */


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


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


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


/* [wxMaxima: comment start ]
We want to verify the expression found in Peskin/Schroeder, p 154, eqn (5.64).
This requires some serious massaging of this function of th.
  
dsigdo is proportional to B defined below,  and
this should be equivalent to the factor Mnfac defined just above
the definition of B.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
combine the three terms of Mn1 which are NOT proportional to k^2 or k^3 or k^4 
(the first, second, and fifth terms) in an expression
we call parts, and then show that expression is equal to 16*k^4.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
for j in [1,2,5] do parts : parts +part(Mn1,j);
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: comment start ]
now let Mn2 be Mn1 with those 3 terms replaced by parts1:

Mn2 is equivalent to Mn1, but in a different form.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Mn2 : Mn1 -parts + parts1;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
now  define Mnfac as what remains of Mn2 if we pull out the factor 32*k^2.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Mnfac :  pullfac(Mn2,32*k^2)  / (32*k^2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
In Peskin/Schroeder, p 154, eqn (5.64),  dsigdo is proportional to what we
call  B here,  and this should be equivalent to the factor Mnfac defined just above.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
B : (E+k)^2 + (E+k*cos(th))^2 -M^2*(1-cos(th));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
let's try to show B is equivalent to Mnfac if we use the information available.
The following should be zero.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
factor( trigsimp(B - Mnfac)) ;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
well, let's include the relation between the muon energy E, its
magnitude of 3-momentum k, and its mass M.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
which shows the equivalence to Peskin/Schroeder p. 154 m --> 0 limit result.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Now for the differential scattering cross section dsigdo in the CMS.

We have absorbed e^4 into the constant A, with e^2 = 4πα.  

In general 2 --> 2 scattering, A = α^2 * (pf/pi) /(4*s), but we have seen above
that conservation of energy implies that pf = pi = k.
 
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
We then go for the unpolaized electron-muon scattering differential cross section
in the CM frame and in the limit m --> 0

Note we return the factor 32 k^2, and use the form B instead of Mnfac.

Averaging over initial spins means we need to divide A by 4
  to get the unpolarized differential scattering cross section.
(we average over the initial spins σ1 and σ2 by summing over
 σ1  and dividing by 2, and like wise for σ2.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
dsigdo : (A/4)*32*k^2*B/t_th^2;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
 where s = E(cm)^2 = (E+k)^2.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Next  look at high energy limit, k >> M, so M --> 0, E --> k, E(cm) --> 2*k
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
dsig_HE : dsigdo, M = 0, E = k, s = 4*k^2;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
which gives the  "high energy limit"  dsigdo.

Behavior as a function of θ as θ --> 0?
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
hence singular behavior:   dsigdo ~  1/θ^4  near θ = 0.
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
UNPOLARIZED CROSS SECTION  WITHOUT Using Mandelstam Variables
   [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) = m^2,
         D(p2,p2) = M^2,
         D(p3,p3) = m^2,
         D(p4,p4) = M^2)$
/* [wxMaxima: input   end   ] */


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


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

As above, <|Mfi|^2> = Mn/t^2 , and Mn is the contraction of the product
of two traces.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
note we get D(pa,pb) 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: comment start ]
We now use noncov (expr) to evaluate this expression in the center of momentum
frame. This works because we have already used comp_def (...) to define the components
of the relevant 4-momenta above.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: subsubsect start ]
Ignore the electron mass
   [wxMaxima: subsubsect end   ] */


/* [wxMaxima: input   start ] */
Mn_meth2 : Mn_meth2, m = 0;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
compare to the method 1 result Mn1 above
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Mn_meth1 : Mn1;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
use our package trigsimp replacement function ts(expr,th) to convert sines to cosines.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
To compare, consider the difference, after replacing E^2 with M^2 + k^2 using our package
function Esq_to_psq (expr, p,E,m)
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
diff12 : Esq_to_psq (Mn_meth2 - Mn_meth1, k,E,M);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
so Mn_meth2 is an equivalent expression but has a different form.

Factor out 32*k^2 from each term using our package function pullfac(expr, stuff)
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Mn_meth2 : pullfac(Mn_meth2, 32*k^2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
This form will yield a different form for the differential scattering cross section:
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
dsigdo_meth2 : (A/4)*Mn_meth2/t_th^2;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
look at the high energy limit of this expression
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
dsigdo_meth2_HE : dsigdo_meth2, M = 0, E = k, s = 4*k^2;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
which agrees with the HE limit of method 1.
   [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  t^2  which requires 
 first computing  t_th : VP (p1-p3,p1-p3).
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(...)
   [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: subsubsect start ]
For Contraction, Use  mcon (expr, index1,index2,...) with traces produced using nc_tr
   [wxMaxima: subsubsect end   ] */


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


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


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


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


/* [wxMaxima: comment start ]
neglect electron mass m as before,
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Mn_meth3 : Mn_meth3, m = 0;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
expand (Mn_meth3 - Mn_meth2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
so they agree.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
m_tr uses the same syntax as  tr 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: comment start ]
m_tr(p3 + m, mu, p1 + m, nu) is converted into the explicit matrix expression:

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

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

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

sL(p3) is the Feynman slash  matrix corresponding to 
the 4-momentum p3, and I4 is the explicit 4 x 4 unit matrix.

Gam[mu] is the explicit matrix representation of γ^μ, for example.
   [wxMaxima: comment end   ] */


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


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


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


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


/* [wxMaxima: comment start ]
neglect electron mass m
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Mn_meth4 : Mn_meth4, m = 0;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
expand (Mn_meth4 - Mn_meth3);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
which confirms the equality.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
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 ]
POLARIZED AMPLITUDES USING EXPLICIT DIRAC SPINORS
   [wxMaxima: section end   ] */


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


/* [wxMaxima: comment start ]
 case  RR --> RR, ie, (+1,+1) --> (+1,+1) polarized amplitude.
    electron 4-momenta: p1 --> p3, 
     muon 4-momenta:   p2 --> p4
      
     First define the needed Dirac spinors and barred spinors for leptons.

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

    If upa is a Dirac spinor (column vector) corresponding to 4-momentum pa, then
    sbar(upa) is the barred Dirac spinor (row vector) upa^{+} γ^0, 
    where (something)^{+} indicates Hermitian conjugate.
    
    Our normalization is  upa_bar . upa = 2m.   
    
    UU(Ea,pa,θ,φ,σ) is a generic Dirac spinor (4-element column vector) corresponding to a lepton with 
    relativistic energy  Ea = sqrt(pa^2 + ma^2), 3-momentum magnitude pa, 
    with direction of 3-momentum defined by the spherical polar angles
     θ and φ, and helicity σ (+1 for R, -1 for L).    
    
    In the center of momentum frame, the incident electron and muon have equal and opposite
    3-momentum vectors with magnitude k, and same for the outgoing electron and muon.
    
    These parameters should agree with those specified in the comp_def statement above.

We neglect the mass of the electron, as above, so  Ee --> k.

Let the energy of the muon be E, so E^2 = k^2 + M^2, and

E(CM) =  E + k,  s = (E + k)^2
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Define  th/2 forms
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
needed Dirac spinors for the RR --> RR case.
th is the scattering angle of the outgoing electron.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
You can use the grind function to see the display2d=false version of  anything.
The version of wxmaxima we are using omits commas from the row vectors above.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
Check our normalization for up1:
   [wxMaxima: comment end   ] */


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


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


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


/* [wxMaxima: comment start ]
which is due to our neglect of the mass of the electron.

Check norm of up2 corresponding to the incoming muon.
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: input   start ] */
up2b . up2;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
Now calculate the (reduced) amplitude Mr =  M_RR_RR = Mn/t,  for the RR --> RR case.

The reduced amplitude is dimensionless, since the numerator has the dimensions of (ubar.u)^2 = m^2
and the denominator is  t with dimension m^2.

For contraction, we need to use mcon(expr, indices)
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: input   start ] */
Mr : Mn/t_th2;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
Convert to a function of th (rather than a function of th/2) using fr_ao2(expr, th)
   [wxMaxima: comment end   ] */


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


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


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


/* [wxMaxima: comment start ]
High energy limit,  k >> M,  E --> k
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Mr_th, E = k;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
%, ratsimp;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Now automate calculation of the reduced amplitudes for given helicity
assignments, defining the function emu_me.

in this definition,  k,  E, and M  are unbound parameters, and  t_th2  was defined above.

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


/* [wxMaxima: input   start ] */
emu_me (s1v,s2v,s3v,s4v) := 
  block ([up1,up2,up3b,up4b],
    up1 : UU(k,k,0,0,s1v),
    up3b : sbar (UU(k,k,th,0,s3v)),
    up2 : UU (E,k,%pi,0,s2v),
    up4b : sbar( UU (E,k,%pi - th,%pi,s4v)),
     mcon ( (up3b . Gam[mu] . up1)*(up4b . Gam[mu] . up2), mu),
     trigsimp(%%)/t_th2,
    ratsimp(%%))$
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: section start ]
Polarized 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) 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 )
  
After printing out the th/2 form of the amplitude, we
use our package function fr_ao2 to convert the
global quantity mssq to a function of th (instead of th/2).
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
mssq is the sum over spins of the square of the (reduced) matrix elements using our Avsq function.
mssq is global parameter, given zero value at start of block statement.
For each set of helicity quantum numbers (σ1,σ2,σ3,σ4), we print out
their values and the reduced amplitude Mr(σ1,σ2,σ3,σ4).
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: comment start ]
We see that when we use the approximation that the electron mass m is zero,
the amplitudes which involve an electron helicity flip are zero.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
As a function of th/2, we have the following result for mssq:
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
convert to a function of th:
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
Simplify this and let Ma be the numerator:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
scanmap (factor, ratsimp(mssq_th)) ;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
mssq, as a function of th, should be compared to the ratio Mn2/t_th^2 computed above
(in Sect. 1.1.3)  with symbolic methods.  Simplify this ratio and let Mb be the numerator:

   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
ratsimp (Mn2/t_th^2);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
Since the denominators are obviously equal, we need to show the numerators are equal.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
expand(Ma - Mb);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
To make the comparison easier, replace E^2 by k^2 + M^2, using our Esq_to_psq function:
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
which confirms the equivalence of the results of the two methods, which produce
different looking answers.
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
We can check the (square) of the polarized amplitude calculated above 
  (using Dirac spinors for given helicity assignments) by using either the
symbolic method using nc_tr(...) plus mcon or the explicit matrix method
 using m_tr(...)  plus mcon. (Do not use tr(...)  and Con for this type of calculation.)

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

In both the nc_tr and m_tr methods, one must use a spin projection operator
 symbolized by S(σ, Sp ) (for the massive lepton or massive anti-lepton case) 
 or S(σ) for the massless lepton case and  S(-σ)  for the massless anti-lepton case. 

We are ignoring the mass m of the electron.

  For the massive case, Sp is the particle spin 4-vector determined by the 
 particle's mass and 4-momentum components.

In general, for a lepton or anti-lepton with mass m > 0, energy E,   3-momentum pvec,
3-momentum magnitude |pvec|, and 3-momentum unit vector
phat = pvec / |pvec|, the associated spin 4-vector has components:

Sp = [ |pvec|/m, (E/m) phat ],  or in more detail

Sp = [ |pvec|/m, (E/m) phat_x,  (E/m) phat_y,  (E/m) phat_z ]

The muon, with mass M > 0,  has 4-momentum components that
allow the definition of Sp2 and Sp4.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
comp_def (   Sp2 (k/M, 0,0, -E/M),
                     Sp4 (k/M, - E*sin(th)/M, 0, - E*cos(th)/M) )$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
The list invarR is only used by nc_tr, not by m_tr.

What is the current value of the list invarR?
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
so we need to reset this list.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
set_invarR(D(p1,p1) = 0,
             D(p2,p2) = M^2,
             D(p3,p3) = 0,
             D(p4,p4) = M^2,
             D(p2,Sp2) = 0,
             D(Sp2,Sp2) = -1,
             D(p4,Sp4) = 0,
             D(Sp4,Sp4) = -1)$
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: subsect start ]
Case RR --> RR using nc_tr(...) plus mcon
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
We write |Mr|^2 = MSQn/t^2, and compute MSQn, the numerator for the case RR --> RR
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
MSQn : mcon (nc_tr (mu,S(1),p1,nu,S(1),p3) * nc_tr (mu,S(1,Sp2),p2+M,nu,S(1,Sp4),p4+M), mu,nu);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
We are going to simplify this using trigsimp, ts, expand, and Esq_to_psq
to replace E^2 in favor of k^2 and M^2, producing what we call MSQn_c
   [wxMaxima: comment end   ] */


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


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


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


/* [wxMaxima: comment start ]
let's compare this numerator with what we got for the RR --> RR case using explicit Dirac spinors.

We repeat our definition of the needed spinors for this case.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
up1 : UU(k,k,0,0,1)$
up3b : sbar (UU(k,k,th,0,1))$
up2 : UU (E,k,%pi,0,1)$
up4b : sbar( UU (E,k,%pi - th,%pi,1))$
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
Now we simplify the Dirac spinor (reduced) amplitude numerator by using ts(expr, angle), 
 fr_ao2(expr,th) [ to convert  f(th/2)  to  g(th) ], then squaring the (real) result and then using 
  Esq_to_psq to replace E^2 in favor of k^2 and M^2.
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: input   start ] */
Mn_b2 : expand(Mn_b^2);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
they look the same, but let's check:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
is(equal( MSQn_c, Mn_c2));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
which confirms equality for this RR --> RR case.
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
Case RR --> RR  using m_tr(...) and mcon
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
The explicit matrix method m_tr(...)  has the same syntax as the symbolic method
nc_tr(...), but is about five times faster in this example.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
MSQmn : mcon (m_tr (mu,S(1),p1,nu,S(1),p3) * m_tr (mu,S(1,Sp2),p2+M,nu,S(1,Sp4),p4+M), mu,nu);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
The use of the symbolic method nc_tr(...) took about five times longer than using
the explicit matrix method m_tr(....) in this example.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
is(equal( MSQmn, MSQn));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
So we get the same raw answer as using nc_tr(...)
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
Case RR --> LR using nc_tr(...) or m_tr(...)  plus mcon
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
In the m = 0 approximation, amplitudes which involve a helicity
flip by the electron are zero, according to our Dirac spinor
calculation. So we should get zero using the trace + contraction approach.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
We write |Mr|^2 = MSQn/t^2, and compute MSQn, the numerator for the case RL --> RR
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: comment start ]
The same numerator using m_tr(...):
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
MSQmn : mcon (m_tr (mu,S(1),p1,nu,S(-1),p3) * m_tr (mu,S(1,Sp2),p2+M,nu,S(1,Sp4),p4+M), mu,nu);
/* [wxMaxima: input   end   ] */


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



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