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

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


/* [wxMaxima: comment start ]
compton1.wxm,   Dec. 6, 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).

 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 ]
 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: 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 arbitrary energy (AE) case, first  in terms
  of the Mandelstam variables s, t, and u, and then  in the lab frame 
  of the initial electron, 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. 3 which employs explicit Dirac spinors and photon polarization
vectors is

M(σ1, σ2, λ1, λ2) = -e^2 Mr
 
and the "reduced amplitude"  Mr = M1 + M2, where
-------------------------------------------------------------------------
M1 =   - (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: input   start ] */
set_invarR (
         D(p1,p1) = m^2,
         D(p1, p2) =  m^2 - t/2,
         D(p1, k1) =   s/2 - m^2/2,
         D(p1, k2) =  m^2/2  - u/2,
         D(p2, p2) = m^2,
         D(p2, k1) =  m^2/2  - u/2,
         D(p2, k2) =  s/2 -  m^2/2,
         D(k1, k1) = 0,
         D(k1, k2) = - t/2,
         D(k2, k2) = 0)$
/* [wxMaxima: input   end   ] */


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


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


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


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


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


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

For later use in solving for the magnitude of the 3-momentum (kp) of the
outgoing photon (k2), in terms of the magnitude of the 3-momentum (k)
of the incident photon (k1) in the rest frame of the initial electron, we 
note that in a general frame, D(p1,k1) is equal to D(p2,k2).
   [wxMaxima: comment end   ] */


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


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


/* [wxMaxima: comment start ]

   [wxMaxima: comment end   ] */


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

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

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

The sum over helicities and photon polarizations of the absolute
value square of the  "reduced" amplitude is written as

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

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


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


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


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


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

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


/* [wxMaxima: input   start ] */
MfiSQ : M11n / (u - m^2)^2  + M22n / (s - m^2)^2 + 
                    M21n / ((u - m^2)*(s - m^2)) + M12n / ((u - m^2)*(s - m^2)), 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.
   [wxMaxima: comment end   ] */


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


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


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

We assume here  the initial electron is at rest.
 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.
Let kp  ("k-prime")  be the energy and also the 3-momentum magnitude of the final photon k2.

   [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 initial electron rest frame.

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

comp_def ( p1 (E1, px, py, pz), ...) for example assigns to the particle
with 4-momentum p1 the energy E1, and the 3-momentum components the
values px, py, and pz.

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

We write the components of the emerging electron p2 using
conservation of 4-momentum: p1 + k1 = p2 + k2.
   [wxMaxima: comment end   ] */


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

comp_def (
            p1( m, 0, 0, 0),
            k1 (k, 0, 0, k),
            k2 (kp, kp*sin(th), 0, kp*cos(th)),
            p2( m + k - kp, - kp*sin(th), 0, k - kp*cos(th)) )$
/* [wxMaxima: input   end   ] */


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


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


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


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


/* [wxMaxima: comment start ]
now define explicit frame dependent functions of th, the angle of 
the 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: comment start ]
Use VP to compute explicit values of the Mandelstam variables s, t, and u
in the rest frame of the initial electron ("lab frame").

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: comment start ]
We can solve for the magnitude of the 3-momentum ( kp ) of the outgoing 
photon k2, in terms of the magnitude of the 3-momentum (k) of the incident
photon k1, by noting that both D(p1,k1) and D(p2,k2) are equal to (s - m^2)/2.

   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
eqn : VP(p1,k1) = VP(p2,k2);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
solns : solve (eqn,kp);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
kp_soln : rhs (solns[1]);
/* [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));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
trigsimp ( ratsubst (kp_soln, kp, %));
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
A simpler way of writing this result is found many places.
We have consulted Peskin & Schroeder, Quantum Field
Theory, pp. 162 - 163.  In making a comparison, note that there
is no φ dependence, and that we can therefore write
  d Ω =  2 π  d ( cos (θ) ) under an integral sign.

This comparison also depends on the expression A_lab which
multiplies < |Mr|^2 > for this process; see below.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
MfiSQ_cmp : 8*(kp/k + k/kp - sin(th)^2 );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
ratsubst (kp_soln,kp,%);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
which shows that the simple form of MfiSQ_cmp is equivalent to our expression
MfiSQ_th, when we make use of the expression kp_soln which connects kp
to k and th.
   [wxMaxima: comment end   ] */


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

We have worked through problem 6.10, parts a and b,  in David Griffiths "Introduction to
Elementary Particles",  pp. 211 - 212. The simplification in part b applies
to our problem. We have identified particle 2, initially at rest, as our initial
electron with mass m2 = m, and particle 4 as the recoiling electron with
mass m4 = m. The incident photon is particle 1 and the scattered photon
is particle 3.  k = magnitude of the incident photon's 3-momentum, and
kp  (k') is the magnitude of the scattered photon's 3-momentum.

Also note M = -e^2 Mr in our work, and e^2 = 4 π α in our units.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
A_lab :  (alpha^2 / (4*m^2) ) * (kp/k)^2;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
Averaging over the two spin states of the incident photon and
the two spin states of the initial electron at rest means that we
must divide the prefactor A_lab by 4.

We then get the "unpolarized differential scattering cross section"
in the rest frame of the initial electron as:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
dsigdo_lab : (A_lab/4) * MfiSQ_cmp;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
In section 2 we obtain this form more straightforwardly by starting with
unpolarized electrons but polarized photons, and starting with the use
of a minimum set of D(pa,pb)'s (and later using D_sub(...) ) rather
than explicit Mandelstam variables. The use of D_sub allows more
control over the final form, and is especially useful when we are using
the rest frame of one of the particles.
   [wxMaxima: comment end   ] */


/* [wxMaxima: subsect start ]
Unpolarized Cross Section  WITHOUT  Mandelstam Variables, using symbolic tr and Con
   [wxMaxima: subsect end   ] */


/* [wxMaxima: comment start ]
reset the list invarR of invariant replacement rules to avoid the Mandelstam variables s, t, and u,
using set_invarR.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
set_invarR (
         D(p1,p1) = m^2,         
         D(p2, p2) = m^2,         
         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 - m^2)^2  +  M22n / (s - m^2)^2  + 
                                                   M21n / ( (u - m^2)*(s - m^2) )  + M12n / ( (u - m^2)*(s - m^2) )

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


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


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


/* [wxMaxima: input   start ] */
M12n :  tr (mu, p1 - k2 + m, nu, p1 + m, mu, p1 + k1 + m, nu, p2 + m);
/* [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, 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 - m^2)^2  + M22n_nc / (s_th - m^2)^2 + 
                    M21n_nc / ((u_th - m^2)*(s_th - m^2)) + M12n_nc / ((u_th - m^2)*(s_th - m^2)));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
As we did in method 1, we replace kp by kp_soln everywhere:
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
MfiSQ_th_meth2 :   trigsimp ( ratsubst (kp_soln, kp, MfiSQ_th_meth2));
/* [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 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 + m, nu, p1 + m, nu, p1 - k2 +m, mu, p2 + m);
/* [wxMaxima: input   end   ] */


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


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


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


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


/* [wxMaxima: input   start ] */
MfiSQ_th_meth3 :   trigsimp ( ratsubst (kp_soln, kp, MfiSQ_th_meth3));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
trigsimp (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] = kp, k1[1] = kp*sin(th), k1[2] = 0, k1[3] = kp*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 + m, nu, p1 + m, nu, p1 - k2 +m, mu, p2 + m), mu, nu));
/* [wxMaxima: input   end   ] */


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


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


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


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


/* [wxMaxima: input   start ] */
MfiSQ_th_meth4 :   trigsimp ( ratsubst (kp_soln, kp, MfiSQ_th_meth4));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
trigsimp (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: comment start ]

   [wxMaxima: comment end   ] */


/* [wxMaxima: section start ]
LINEARLY POLARIZED PHOTONS but UNPOLARIZED ELECTRONS Using D_sub
   [wxMaxima: section end   ] */


/* [wxMaxima: comment start ]
In this section we use REAL linear polarization 4-vectors for the photons, which we
call e1 and e2,  so we start with

M(σ1, σ2, λ1, λ2) = -e^2 Mr
 
and the "reduced amplitude"  Mr = M1 + M2, where
-------------------------------------------------------------------------
M1 =   - (u2bar SL( e1 ) (SL (p1 - k2) + m*I4) SL (e2) u1) / [ (p1 - k2)^2 - m^2]
      =    (u2bar SL( e1 ) (SL (p1 - k2) + m*I4) SL (e2) u1) / 2*D(p1,k2)

M2 =   - (u2bar SL( e2 ) (SL (p1 + k1) + m*I4) SL (e1) u1) / [ (p1 + k1)^2 - m^2]
      =    - (u2bar SL( e2 ) (SL (p1 + k1) + m*I4) SL (e1) u1) / 2*D(p1,k1)
----------------------------------------------------------------------------

If we sum over the helicities σ1 and σ2 of the initial and final electrons, but
leave the spins  λ1 and λ2 of the initial and final photons arbitrary, we get
the square of the reduced amplitude in a form (we absorb the minus sign in
the cross terms into the definition of M21n and M12n)

MfiSQ_pol = M11n / (4*D(p1,k2)^2) + M22n  / (4*D(p1,k1)^2) + M21n / (4*D(p1,k2)*D(p1,k1))  + M12n / (4*D(p1,k2)*D(p1,k1))

In the above form, the photon linear polarization 4-vectors remain within Feynman slashes.

We will use our symbolic tr function which will produce expressions which are functions
of D(pa,pb).
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
The symbolic tr function will make use of the list invarR, so we can provide simplifications
for all of the D(pa,pb)'s except for two in the list invarR.

Just as we can express all products of 4-momenta D(pa,pb) in terms of any
pair of the Mandelstam variables, say u and s, we can instead  express 
all such products in terms of D(k2,p1) and D(k1,p1), with the goal in view
to find the differential scattering cross section in the rest frame of the
initial electron p1, in which the components of p1 are [m,0,0,0].

Thus since s = (p1 + k1)^2 = (p2 + k2)^2 by conservation of 4-momentum, we
have s = m^2 + 2*D(p1,k1) = m^2 + 2*D(p2,k2) so we can replace
   D(p2,k2) >> D(p1,k1). 

Likewise since u = (p1 - k2)^2 = (p2 - k1)^2, we have u = m^2 -2*D(p1,k2) = 
  m^2 -2*D(p2,k1), we can replace D(k1,p2)   >> D(k2,p1).

Also t = (p1 - p2)^2 = (k1 - k2)^2 implies:
    D(k1,k2) = D(p1,p2) - m^2,  come back to this and first evaluate:

  D(p1,p2) = D (p2, p1 + k1 - k2) = m^2 + D(p1,k1) - D(p1,k2),
  and we finally get

   D(k1,k2) = D(p1,k1) -  D(p1,k2)

   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
set_invarR (
         D(p1,p1)=m^2,
         D(p2,p2)=m^2,
         D(k1,k1) = 0, 
         D(k2,k2) = 0,
         D(p1,p2) = m^2 + D(k1,p1) - D(k2,p1),
         D(k1,p2) = D(k2,p1),
         D(k2,p2) = D(k1,p1),
         D(k1,k2) = D(k1,p1) - D(k2,p1) )$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
length(invarR);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
Next we need to ADD to the list invarR replacements for the products of
the (real) linear polarization 4-vectors e1 and e2 with themselves and with
the particle 4-momenta p1, p2, k1, k2.

In the following, the first four reflect basic properties of the two photon polarization
4-vectors (see photon2.wxm), the second two reflect the facts that we are ANTICIPATING
that  p1 = [m,0,0,0],  as well as the fact that we are  using transverse gauge real linear 
polarization 4-vectors, so that e1[0] = 0 and e2[0] = 0, so D(p1,e1) = 0 = D(p1,e2).

Finally, the last two replacements come from using conservation of 4-momentum in the form

D(e1,p2) = D(e1,p1 + k1 - k2), and D(e2,p2) = D(e2,p1 + k1 - k2),  plus the properties
of the polarization 4-vectors.

We are letting the real 4-vector e1 represent ek1 and e2 represent ek2.

By using the syntax invar (....) we can ADD to the list invarR, as we will check.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
invar (D(e1,e1) = -1,
        D(k1,e1) = 0,
        D(e2,e2) = -1,
        D(k2,e2) = 0,
        D(p1,e1) = 0,
        D(p1,e2) = 0,
        D(e1,p2) = -D(e1,k2),
        D(e2,p2) = D(e2,k1))$ 
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
length (invarR);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: input   start ] */
M11n : tr (e1,p1 - k2 + m, e2, p1 + m, e2, p1 - k2 + m, e1, p2 + m);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
M22n : tr (e2,p1 + k1 + m, e1, p1 + m, e1, p1 + k1 + m, e2, p2 + m);
/* [wxMaxima: input   end   ] */


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


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


/* [wxMaxima: input   start ] */
expand (M21n - M12n);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
The  quantities M21n and M12n are equal,  so we can write
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
MfiSQ_pol : expand ( M11n / (4*D(p1,k2)^2)  + M22n / (4*D(p1,k1)^2) + 
                              2*M21n / (4*D(p1,k2)*D(p1,k1)) );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
we have alredy used comp_def above to define the components of the momentum
4-vectors, for example:
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: comment start ]
The dirac3 package function D_sub (expr, D or list of D's) 
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
D_sub (MfiSQ_pol, D(k1,p1) );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
MfiSQ_pol : D_sub (MfiSQ_pol, [D(k1,p1),D(k2,p1)]);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
MfiSQ_pol : pullfac (MfiSQ_pol, 2);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
For fixed polarization states of the initial and final photons, and
averaging over spin states of initial electron, and summing over
final spin states of final electron, we need to divide A_lab by 2.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
dsigdo_lab_pol : (A_lab/2)*MfiSQ_pol;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
which agrees with Weinberg (8.7.39) [The Quantum Theory
of Fields, Vol. I]

We can call the above dsigdo[s, r] where s = (1,2) specifies the
linear polarization state of the incident photon, and r = (1,2)
specifies the linear polarization state of the outgoing photon.

If we now average over the unknown polarization state s of the
incident photon, we get dsigdo[r] for a given polarization state
of the outgoing photon r.

dsigdo[r] = (1/2) sum( dsigdo[s,r],s,1, 2)  == < dsigdo[s,r] >_{s}

Then use < B* (C + f(s,r) ) >_{s} = B * (C  + < f(s,r)>_{s} ) where
B and C do not depend on s or r.

Using lists to represent vectors, we show in photon1.wxmx that
< 4 * ( ek[s] . ekp[r] )^2 >_{s} = 2 - 2 * ( khat . ekp-vec[r] )^2

We then recover formula (8.7.40) in Weinberg. 

If we also sum over the final polarization state of the outgoing photon,
we can use Theorem 1 in photon1.wxmx, which looks like

sum (sum ( (ek[s] . ekp[r])^2,s,1,2),r,1,2) =  1 + cos(th)^2

and can then recover Weinberg (8.7.41).
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
dsigdo_e2 : subst (4*D(e1,e2)^2 = 2 - 2*dot3(ku,e2_vec)^2, dsigdo_lab_pol);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
which agrees with formula (8.7.40) in Weinberg, The Quantum Theory of Fields, Vol. I.
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
To recover the unpolarized result we considered in Sec. 1 above, we use Theorem 1
proved in photon1.wxm :

 sum ( sum (  (e1[s] . e2[r] )^2, s,1,2), r,1,2)  = 1 + cos(th)^2

We start with 

dsigdo_unpol =    (1/4)  Σ_{σ1, s, σ2, r}  dsigdo (σ1, s; σ2, r) = (A_lab/4) Σ_{σ1, s, σ2, r}  |Mr (σ1, s; σ2, r)|^2

Our result for MfiSQ_pol above was the result of summing just over σ1, and σ2.

MfiSQ_pol =  Σ_{σ1, σ2 }  |Mr (σ1, s; σ2, r)|^2
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
Msq_const : Mtemp - part(Mtemp,3);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
Msq_pol : subst(D(e1,e2)^2 = 1 + cos(th)^2, part(Mtemp,3));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
MfiSQ_unpol : expand ( 4*Msq_const + Msq_pol);
/* [wxMaxima: input   end   ] */


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


/* [wxMaxima: comment start ]
Since we are averaging over both the two polarization states of the initial electron and
the two pol. states of the  initial photon, we must divide A_lab by 2*2 = 4.
   [wxMaxima: comment end   ] */


/* [wxMaxima: input   start ] */
dsigdo_unpol : (A_lab/4) * MfiSQ_unpol;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: comment start ]
which agrees with Weinberg (8.7.41)
   [wxMaxima: comment end   ] */


/* [wxMaxima: comment start ]
Our previous expression had sin(th)^2 in it
   [wxMaxima: comment end   ] */


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


/* [wxMaxima: input   start ] */
trigsimp (dsigdo_unpol - dsigdo_lab);
/* [wxMaxima: input   end   ] */


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

In compton2.wxm we consider calculating helicity amplitudes, in the
center of momentum frame, using explicit Dirac spinors and photon
polarization vectors. 
   [wxMaxima: comment end   ] */



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