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

 π(+, p1) + π(+, p2) --> π(+, p3) + π(+, p4)

Pion Pion Elastic QED Scattering
   [wxMaxima: title   end   ] */


/* [wxMaxima: comment start ]
pi-pi-scatt.wxm,  Dec. 22, 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: comment start ]
We ignore the finite size and quark content (and potential strong 
interaction effects in small impact parameter collisions) in this 
example of scalar QED, and only take into account the first order
electromagnetic interactions of these spinless bosons.
   [wxMaxima: comment 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 ]
Dirac3 package functions used here are set_invarR, ev_Ds, pullfac,  comp_def,
VP, sub_stu, and div_tb.

   [wxMaxima: comment end   ] */


/* [wxMaxima: section start ]
Center of Momentum Cross Section Without Using Mandelstam Variables
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
The fastest approach is to express the reduced amplitude in terms of D(pa,pb)'s, 
define 4-momentum components using comp_def, and then use noncov_ratio on the result.

D(pa, pb) is a symbolic representation of the 4-vector dot product of
4-momenta pa and pb.  D(pb, pa) = D(pa, pb) since we have declared 
D symmetric in our package.

We define Mr as the "reduced amplitude" for the scattering process,
  where Mfi is the amplitude:
            Mfi =  -e^2 Mr

In lowest order there are two diagrams contributing. 
  Mr = M1 + M2.

   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
M1 : D(p1+p3,p2+p4) / D(p1 - p3, p1 - p3);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M2 :  D(p1+p4,p2+p3) / D(p1 - p4, p1 - p4);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
We next define 4-momentum components in the CM frame, using our function
 comp_def  ("components-define") to define arrays which hold the components 
of the 4-momentum vectors in the CM frame,.
   
comp_def ( pa(Ea, pax, pay, paz), ....) is the syntax used here.

   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
assume ( p > 0, th >= 0, th <= %pi )$
comp_def ( p1( E,0,0,p),
             p2( E,0,0,-p),
             p3 (E,p*sin(th),0,p*cos(th)),
             p4 (E,-p*sin(th),0,-p*cos(th)) )$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Now we use noncov_ratio on each of M1 and M2 and add the
result to get Mr in terms of CMS parameters.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Mr_th : noncov_ratio (M1) + noncov_ratio (M2), ratsimp;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Use our Dirac3 package function div_tb(expr, fac) to divide both
the top and the bottom by p^2.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
fundef (div_tb);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
can we simplify this?
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
trigsimp (cos(th)^2-1);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
ts (expr, th) is a Dirac3 package function which is an alternative to trigsimp.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
ts (cos(th)^2-1, th);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
Mr_th : pullfac (num(Mr_th),2) / ts (denom (Mr_th), th ) ;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
The above is the "reduced amplitude" Mr for arbitrary energies.

For high energy collisions in CM frame, neglect masses, and the particle
energy E is the same as the 3-momentum magnitude p, and Mr becomes:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
Mr_th_HE :  Mr_th, p = E;
/* [wxMaxima: input   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) in cms.

Let p be  the magnitude of the 3-momentum of the initial pion and
initial Kaon and pf be the magnitude of the pion and final Kaon in CMS.
Conservation of energy implies pf = pi = p, so  A = α^2/( 4 s).

   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
s is the Mandelstam variable (p1 + p2)^2  = E(cms)^2
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
For arbitrary energies (note Mr_th is real):
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
dsigdo : A * Mr_th^2;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
dsigdo_HE : A * Mr_th_HE^2;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: subsect start ]
Use set_invarR to construct list invarR; Use ev_Ds( expr ) to replace
4-momentum dot products with Mandelstam variables s, t, and u.
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
Use our  set_invarR function to express all invariant dot products of 4-momenta 
in terms of the Mandelstam variables s, t,  and u, as well as the pion mass m.
The use of the function set_invarR in this way  populates a list called invarR, 
which is later used by other package functions to transform expressions.

Recall conservation of 4-momentum: p1^μ + p2^μ = p3^μ + p4^μ
 and our definitions of the Mandelstam variables for this process:

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

  s, t, and u are not all independent since the sum s + t + u = sum of mass-square,
  or here:
   s + t + u = 4 m^2 

The following dot product replacement rules come from the definitions
of s, t, and u, as well as our chosen symbols for the pion mass.

   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: comment start ]
This example shows the use of ev_Ds (expr) function to replace all D(pa,pb) symbols with expressions
in terms of the Mandelstam variables s, t, and u, once the list invarR has
been populated using the function set_invarR. Also the use of pullfac(expr,fac) to
pull a factor (fac) out of expression expr.

We define Mr as the "reduced amplitude" for the scattering process,
  where Mfi is the amplitude:
            Mfi =  -e^2 Mr

In lowest order there are two diagrams contributing. 
  Mr = M1 + M2.

   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
M1 : D(p1+p3,p2+p4)/t;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
M2 :  D(p1+p4,p2+p3)/u;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
M1 : pullfac (M1, 1/t);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M2 : pullfac (M2, 1/u);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: subsect start ]
Use sub_stu (expr) to replace s --> s_th, t --> t_th, u --> u_th
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
We now want to evaluate the (reduced) amplitude Mr in terms of CM frame parameters.
Our package has a function sub_stu(expr) which replaces, in expr, the
  symbol s by s_th, replaces the symbol t by t_th, and replaces the symbol
 u by u_th. We can use fundef to see the function definition:

   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
Note: trigsimp(%%) in this context means evaluate trigsimp of the previous
 expression in this function definition.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
If you have previously bound values to the three symbols s_th, t_th, 
and u_th, the expression produced by sub_stu will reflect those bindings.

We have already used comp_def to define the 4-momentum components.

   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
Use VP (pa, pb) to evaluate 4-vector dot product in terms of defined components.
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
In order to use sub_stu(expr) to substitute values of s_th, t_th, and
 u_th  (Mandelstam variables in terms of frame dependent parameters) for s,
 t, and u, it is easiest to use our function VP (pa,pb)  (Vector Product)
 for evaluation of the vector products of 4-momenta.

   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: comment start ]
The use of VP(pa,pb) requires that the components of the 4-momentum vectors
have been defined, using comp_def(...), as was done above.   

The function comp_def creates a separate array for each 4-momentum vector, and
  you can check the components using listarray(pa), and/or pn[m], (m = 0,1,2,3)
 noting that pn[0] returns the first element of the array (i.e., ``zero based array'').  

The list momRS holds the names of the arrays defined so far by comp_def.

   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: input   start ] */
for n:0 thru 3 do print ('p4[n]," = ", p4[n])$
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: comment start ]
Returning to the task of giving values to the symbols s_th, t_th, and u_th:
From the definitions of the Mandelstam variables:
   [wxMaxima: comment end   ] */


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


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


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


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


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


/* [wxMaxima: input   start ] */
Mr_th :  trigsimp (  sub_stu(Mr));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
We also use here the package function div_tb(expr, fac),  (divide top and bottom
by the same thing) which returns the result of dividing both the numerator 
and denominator (top and bottom) of expr by fac. 

   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
Mr_th :   pullfac (num(Mr_th), -2) / denom(Mr_th);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
The above is the "reduced amplitude" Mr for arbitrary energies.

For high energy collisions in CM frame, neglect masses, and the particle
energy E is the same as the 3-momentum magnitude p, and Mr becomes:
   [wxMaxima: comment end   ] */


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


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


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


/* [wxMaxima: comment start ]
We can replace s by s_th,   E(cms)^2: 
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
For arbitrary energies (note Mr_th is real):
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
dsigdo : A * Mr_th^2;
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
dsigdo_HE : A * Mr_th_HE^2;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]

   [wxMaxima: comment end   ] */



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