Program for calculating the Dirac equation on curved spacetimes.
Here the program is adapted to a Taub-NUT spacetime !
See gr-qc/0209096 and gr-qc/0010085
> restart;
> grtw();
Define first the Pauli matrices (denoted with "sigma") algebra and their product ("pr") and direct product ("pd"). Then we need an product operator for the direct product (!) ("&p"). Remark that this product is NOT automatically linear, so it is necessary to define all his properties!
> define(sigma,sigma(0)=1);
> define(pr,pr(1,1)=1,pr(1,sigma(1))=sigma(1),pr(1,sigma(2))=sigma(2),pr(1,sigma(3))=sigma(3),pr(sigma(1),1)=sigma(1),pr(sigma(2),1)=sigma(2),pr(sigma(3),1)=sigma(3),pr(sigma(2),sigma(1))=-pr(sigma(1),sigma(2)),pr(sigma(1),sigma(2))=I*sigma(3),pr(sigma(3),sigma(1))=-pr(sigma(1),sigma(3)),pr(sigma(1),sigma(3))=-I*sigma(2),pr(sigma(3),sigma(2))=-pr(sigma(2),sigma(3)),pr(sigma(2),sigma(3))=I*sigma(1),pr(sigma(1),sigma(1))=1,pr(sigma(2),sigma(2))=1,pr(sigma(3),sigma(3))=1,pr((a::integer)*sigma(b::algebraic),(d::integer)*sigma(c::integer))=a*d*pr(sigma(b),sigma(c)));
> define(pd,pd(0,a::algebraic)=0,pd(a::algebraic,0)=0,pd(1,1)=1,pd(I*a::algebraic,b::algebraic)=I*pd(a,b),pd(-I*a::algebraic,b::algebraic)=-I*pd(a,b),pd(a::algebraic,I*b::algebraic)=I*pd(a,b),pd(a::algebraic,-I*b::algebraic)=-I*pd(a,b),pd(-a::algebraic,b::algebraic)=-pd(a,b),pd(a::algebraic,-b::algebraic)=-pd(a,b));
> define(`&p`,`&p`(-a::algebraic,-b::algebraic)=`&p`(a,b),(-pd(sigma(a::algebraic),sigma(b::algebraic))) &p (-pd(sigma(c::algebraic),sigma(d::algebraic)))=pd(sigma(a),sigma(b)) &p pd(sigma(c),sigma(d)),(-pd(a::algebraic,b::algebraic)) &p (-m::function*pd(c::algebraic,d::algebraic))=m*(pd(a,b) &p pd(c,d)),pd(a::algebraic,b::algebraic) &p (-m::function*pd(c::algebraic,d::algebraic))=-m*(pd(a,b) &p pd(c,d)),(-pd(a::algebraic,b::algebraic)) &p pd(c::algebraic,d::algebraic)=-(pd(a,b) &p pd(c,d)),pd(sigma(a::algebraic),sigma(b::algebraic)) &p pd(sigma(c::algebraic),sigma(d::algebraic))=pd(pr(sigma(a),sigma(c)),pr(sigma(b),sigma(d))),pd(1,sigma(a::algebraic)) &p pd(sigma(b::algebraic),sigma(c::algebraic))=pd(pr(1,sigma(b)),pr(sigma(a),sigma(c))),pd(sigma(a::algebraic),sigma(b::algebraic)) &p pd(1,sigma(c::algebraic))=pd(pr(sigma(a),1),pr(sigma(b),sigma(c))),pd(sigma(a::algebraic),1) &p pd(sigma(b::algebraic),sigma(c::algebraic))=pd(pr(sigma(a),sigma(b)),pr(1,sigma(c))),pd(sigma(a::algebraic),sigma(b::algebraic)) &p pd(sigma(c::algebraic),1)=pd(pr(sigma(a),sigma(c)),pr(sigma(b),1)),pd(1,sigma(a::algebraic)) &p pd(1,sigma(b::algebraic))=pd(pr(1,1),pr(sigma(a),sigma(b))),pd(sigma(a::algebraic),1) &p pd(sigma(b::algebraic),1)=pd(pr(sigma(a),sigma(b)),pr(1,1)),I*(a::algebraic) &p (I*b::algebraic)=-(a &p b),(I*a::algebraic) &p (b::algebraic)=I*(a &p b),(a::algebraic) &p (I*b::algebraic)=I*(a &p b));
>
> definemore(`&p`,`&p`(c::algebraic*pd(sigma(a::algebraic),sigma(b::algebraic)),d::algebraic*pd(sigma(e::algebraic),sigma(f::algebraic)))=c*d*`&p`(pd(sigma(a),sigma(b)),pd(sigma(e),sigma(f))),`&p`(c::algebraic*pd(a::algebraic,sigma(b::algebraic)),d::algebraic*pd(sigma(e::algebraic),sigma(f::algebraic)))=c*d*`&p`(pd(a,sigma(b)),pd(sigma(e),sigma(f))),`&p`(c::algebraic*pd(sigma(a::algebraic),b::algebraic),d::algebraic*pd(sigma(e::algebraic),sigma(f::algebraic)))=c*d*`&p`(pd(sigma(a),b),pd(sigma(e),sigma(f))),`&p`(c::algebraic*pd(sigma(a::algebraic),sigma(b::algebraic)),d::algebraic*pd(e::algebraic,sigma(f::algebraic)))=c*d*`&p`(pd(sigma(a),sigma(b)),pd(e,sigma(f))),`&p`(c::algebraic*pd(sigma(a::algebraic),sigma(b::algebraic)),d::algebraic*pd(sigma(e::algebraic),f::algebraic))=c*d*`&p`(pd(sigma(a),sigma(b)),pd(sigma(e),f)),`&p`(pd(sigma(a::algebraic),1),pd(1,sigma(b::algebraic)))=pd(pr(sigma(a),1),pr(1,sigma(b))),`&p`(c::algebraic*pd(sigma(a::algebraic),e::algebraic),d::algebraic*pd(f::algebraic,sigma(b::algebraic)))=c*d*`&p`(pd(sigma(a),e),pd(f,sigma(b))),`&p`(a::algebraic,0)=0 );
Just some checkings
> (-3*x*pd(sigma(1),sigma(2))) &p (-r*pd(sigma(0),sigma(3)));
Define now the Dirac gamma matrices and their commutator with "&p "
> define(gam,gam(1)=I*pd(sigma(2),sigma(1)),gam(2)=I*pd(sigma(2),sigma(2)),gam(3)=I*pd(sigma(2),sigma(3)),gam(0)=pd(sigma(1),1),gam(5)=-pd(sigma(3),1));
> define(comu,comu(a::algebraic,b::algebraic)=a &p b - b &p a);
Again a small test
> comu(gam(1),gam(2));
> I*comu(gam(2),gam(1))/4;
Download now the metric !
> qload(tnb_1);
> grcalc(metric);
Calculated e(bup,dn) for tnb_1 (.009000 sec.)
Calculated g(dn,dn) for tnb_1 (.012000 sec.)
> gralter(metric,simplify);
Component simplification of a GRTensorII object:
Applying routine simplify to object g(dn,dn)
> grdisplay(metric);
> grcalc(ds);
Calculated ds for tnb_1 (.002000 sec.)
> gralter(ds,trigsin,expand,normal);
Component simplification of a GRTensorII object:
Applying routine `simplify[trigsin]` to object ds
Applying routine expand to object ds
Applying routine normal to object ds
> grdisplay(ds);
Define and compute (one by one) the components of the SL(2C) generators
> grdef(`SS{ ^a ^b }`);
Created definition for SS(up,up)
> grcalc(SS(up,up)):
Enter components for object SS(up,up)
If you wish to quit at any point and leave this object
uninitialized, enter the string exit.
REMEMBER to complete each entry with a semicolon.
grcalc> (I/4)*comu(gam(1),gam(1));
grcalc> (I/4)*comu(gam(2),gam(1));
grcalc> (I/4)*comu(gam(3),gam(1));
grcalc> (I/4)*comu(gam(0),gam(1));
grcalc> (I/4)*comu(gam(1),gam(2));
grcalc> (I/4)*comu(gam(2),gam(2));
grcalc> (I/4)*comu(gam(3),gam(2));
grcalc> (I/4)*comu(gam(0),gam(2));
grcalc> (I/4)*comu(gam(1),gam(3));
grcalc> (I/4)*comu(gam(2),gam(3));
grcalc> (I/4)*comu(gam(3),gam(3));
grcalc> (I/4)*comu(gam(0),gam(3));
grcalc> (I/4)*comu(gam(1),gam(0));
grcalc> (I/4)*comu(gam(2),gam(0));
grcalc> (I/4)*comu(gam(3),gam(0));
grcalc> (I/4)*comu(gam(0),gam(0));
Calculated SS(up,up) for tnb_1 (.579000 sec.)
> grdisplay(SS(up,up));
The Christoffell symbols
> grcalc(Chr(dn,dn,dn));
Calculated g(dn,dn,pdn) for tnb_1 (.005000 sec.)
Calculated Chr(dn,dn,dn) for tnb_1 (.005000 sec.)
> gralter(Chr(dn,dn,dn),trigsin);
Component simplification of a GRTensorII object:
Applying routine `simplify[trigsin]` to object Chr(dn,dn,dn)
> grdisplay(Chr(dn,dn,dn));
> grcalc(Chr(bdn,bdn,bdn));
Created definition for Chr(bdn,bdn,bdn)
Calculated detg for tnb_1 (.003000 sec.)
Calculated g(up,up) for tnb_1 (.024000 sec.)
Calculated e(bdn,up) for tnb_1 (.008000 sec.)
Calculated Chr(dn,dn,up) for tnb_1 (.025000 sec.)
Calculated Chr(bdn,bdn,bdn) for tnb_1 (.234000 sec.)
> gralter(Chr(bdn,bdn,bdn),trigsin,expand,simplify);
Component simplification of a GRTensorII object:
Applying routine `simplify[trigsin]` to object Chr(bdn,bdn,bdn)
Applying routine expand to object Chr(bdn,bdn,bdn)
Applying routine simplify to object Chr(bdn,bdn,bdn)
> grdisplay(Chr(bdn,bdn,bdn));
The special term from eq. (9) of gr-qc/0209096
> grdef(`de{ i }:=(I/2)*SS{ ^a ^b }*Chr{ (a) (i) (b) }`);
Created definition for de(dn)
> grcalc(de(dn));
Calculated de(dn) for tnb_1 (.014000 sec.)
> gralter(de(dn),factor,trigsin,expand);
Component simplification of a GRTensorII object:
Applying routine factor to object de(dn)
Applying routine `simplify[trigsin]` to object de(dn)
Applying routine expand to object de(dn)
> grdisplay(de(dn));
The Dirac gamma matrices and then testing their properties
> grdef(`ga{ ^a }:=[gam(1),gam(2),gam(3),gam(0)]`);
Components assigned for metric: tnb_1
Created definition for ga(up)
> grdisplay(ga(up));
> grdef(`ver{ ^a ^b }`);
Created definition for ver(up,up)
> grcalc(ver(up,up));
Enter components for object ver(up,up)
If you wish to quit at any point and leave this object
uninitialized, enter the string exit.
REMEMBER to complete each entry with a semicolon.
grcalc>
grcomponent(ga(up),[y]) &p grcomponent(ga(up),[y]) +grcomponent(ga(up),[y]) &p grcomponent(ga(up),[y])-2*grcomponent(eta(bup,bup),[y,y]);
grcalc>
grcomponent(ga(up),[theta]) &p grcomponent(ga(up),[y]) +grcomponent(ga(up),[y]) &p grcomponent(ga(up),[theta])-2*grcomponent(eta(bup,bup),[theta,y]);
grcalc>
grcomponent(ga(up),[phi]) &p grcomponent(ga(up),[y]) +grcomponent(ga(up),[y]) &p grcomponent(ga(up),[phi])-2*grcomponent(eta(bup,bup),[phi,y]);
grcalc>
grcomponent(ga(up),[t]) &p grcomponent(ga(up),[y]) +grcomponent(ga(up),[y]) &p grcomponent(ga(up),[t])-2*grcomponent(eta(bup,bup),[t,y]);
grcalc>
grcomponent(ga(up),[y]) &p grcomponent(ga(up),[theta]) +grcomponent(ga(up),[theta]) &p grcomponent(ga(up),[y])-2*grcomponent(eta(bup,bup),[y,theta]);
grcalc>
grcomponent(ga(up),[theta]) &p grcomponent(ga(up),[theta]) +grcomponent(ga(up),[theta]) &p grcomponent(ga(up),[theta])-2*grcomponent(eta(bup,bup),[theta,theta]);
grcalc>
grcomponent(ga(up),[phi]) &p grcomponent(ga(up),[theta]) +grcomponent(ga(up),[theta]) &p grcomponent(ga(up),[phi])-2*grcomponent(eta(bup,bup),[phi,theta]);
grcalc>
grcomponent(ga(up),[t]) &p grcomponent(ga(up),[theta]) +grcomponent(ga(up),[theta]) &p grcomponent(ga(up),[t])-2*grcomponent(eta(bup,bup),[t,theta]);
grcalc>
grcomponent(ga(up),[y]) &p grcomponent(ga(up),[phi]) +grcomponent(ga(up),[phi]) &p grcomponent(ga(up),[y])-2*grcomponent(eta(bup,bup),[y,phi]);
grcalc>
grcomponent(ga(up),[theta]) &p grcomponent(ga(up),[phi]) +grcomponent(ga(up),[phi]) &p grcomponent(ga(up),[theta])-2*grcomponent(eta(bup,bup),[theta,phi]);
grcalc>
grcomponent(ga(up),[phi]) &p grcomponent(ga(up),[phi]) +grcomponent(ga(up),[phi]) &p grcomponent(ga(up),[phi])-2*grcomponent(eta(bup,bup),[phi,phi]);
grcalc>
grcomponent(ga(up),[t]) &p grcomponent(ga(up),[phi]) +grcomponent(ga(up),[phi]) &p grcomponent(ga(up),[t])-2*grcomponent(eta(bup,bup),[t,phi]);
grcalc>
grcomponent(ga(up),[y]) &p grcomponent(ga(up),[t]) +grcomponent(ga(up),[t]) &p grcomponent(ga(up),[y])-2*grcomponent(eta(bup,bup),[y,t]);
grcalc>
grcomponent(ga(up),[theta]) &p grcomponent(ga(up),[t]) +grcomponent(ga(up),[t]) &p grcomponent(ga(up),[theta])-2*grcomponent(eta(bup,bup),[theta,t]);
grcalc>
grcomponent(ga(up),[phi]) &p grcomponent(ga(up),[t]) +grcomponent(ga(up),[t]) &p grcomponent(ga(up),[phi])-2*grcomponent(eta(bup,bup),[phi,t]);
grcalc>
grcomponent(ga(up),[t]) &p grcomponent(ga(up),[t]) +grcomponent(ga(up),[t]) &p grcomponent(ga(up),[t])-2*grcomponent(eta(bup,bup),[t,t]);
Calculated ver(up,up) for tnb_1 (.545000 sec.)
> grdisplay(ver(up,up));
The wave function as a covariant vector
> grdef(`Psid{ a }:=[diff(psi(t),y),diff(psi(t),theta),diff(psi(t),phi),diff(psi(t),t)]`);
Components assigned for metric: tnb_1
Created definition for Psid(dn)
> grcalc(Psid(dn));
> grdisplay(Psid(dn));
> grcalc(Psid(bdn));
Created definition for Psid(bdn)
Calculated Psid(bdn) for tnb_1 (.003000 sec.)
> gralter(Psid(bdn),trigsin,simplify);
Component simplification of a GRTensorII object:
Applying routine `simplify[trigsin]` to object Psid(bdn)
Applying routine simplify to object Psid(bdn)
> grdisplay(Psid(bdn));
A series of tricks to obtain the Dirac operator ("ddd" see below) in monomial form
> a0:=expand(grcomponent(de(dn),[t]));a00:=0;
> u0:=whattype(a0);u0;
> nops(a0);
> if u0=`+` then for i from 1 to nops(a0) do a00:=a00+I*h*grcomponent(ga(up),[t]) &p op(i,a0) od else a00:=I*h*grcomponent(ga(up),[t]) &p a0 fi;
>
> a00;
> a1:=expand(grcomponent(de(dn),[y]));a11:=0;
> u1:=whattype(a1);u1;
> nops(a1);
> if u1=`+` then for i from 1 to nops(a1) do a11:=a11+I*h*grcomponent(ga(up),[y]) &p op(i,a1) od else a11:=I*h*grcomponent(ga(up),[y]) &p a1 fi;
> a11;
> a2:=expand(grcomponent(de(dn),[theta]));a22:=0;
> u2:=whattype(a2);u2;
> nops(a2);
> if u2=`+` then for i from 1 to nops(a2) do a22:=a22+I*h*grcomponent(ga(up),[theta]) &p op(i,a2) od else a22:=I*h*grcomponent(ga(up),[theta]) &p a2 fi;
> a22;
> a3:=expand(grcomponent(de(dn),[phi]));a33:=0;
>
u3:=whattype(a3);u3;
> nops(a3);
> if u3=`+` then for i from 1 to nops(a3) do a33:=a33+I*h*grcomponent(ga(up),[phi]) &p op(i,a3) od else a33:=I*h*grcomponent(ga(up),[phi]) &p a3 fi;
> a33;
> grdef(`ddd`);
Created definition for ddd
> grcalc(ddd);
Enter components for object ddd
If you wish to quit at any point and leave this object
uninitialized, enter the string exit.
REMEMBER to complete each entry with a semicolon.
grcalc> a00+a11+a22+a33;
Calculated ddd for tnb_1 (.022000 sec.)
> grdisplay(ddd);
> gralter(ddd,trigsin,factor,expand);
Component simplification of a GRTensorII object:
Applying routine `simplify[trigsin]` to object ddd
Applying routine factor to object ddd
Applying routine expand to object ddd
> grdisplay(ddd);
> grdisplay(ga(up));
And now finally the Dirac equation
> grdef(`dirac:=I*h*ga{ ^l }*Psid{ (l) }+ddd*psi(t)-m*c*psi(t)`);
Created definition for dirac
> grdisplay(ddd);
> grcalc(dirac);
Calculated dirac for tnb_1 (.021000 sec.)
> gralter(dirac,trigsin,factor,power,expand);
Component simplification of a GRTensorII object:
Applying routine `simplify[trigsin]` to object dirac
Applying routine factor to object dirac
Applying routine `simplify[power]` to object dirac
Applying routine expand to object dirac
> grdisplay(dirac);
New definitions for providing the Dirac equation in a more comprehensible form
> define(`gen`);
> define(`gama`);
> grmap(dirac,subs,pd(sigma(1),1)=gama(0),pd(sigma(2),1)=-I*gama(0)*gama(5),pd(sigma(3),1)=-gama(5),pd(1,sigma(1))=2*gen(2,3),pd(1,sigma(2))=2*gen(3,1),pd(1,sigma(3))=2*gen(1,2),pd(sigma(2),sigma(1))=-I*gama(1),pd(sigma(2),sigma(2))=-I*gama(2),pd(sigma(2),sigma(3))=-I*gama(3),pd(sigma(1),sigma(1))=gama(1)*gama(5),pd(sigma(1),sigma(2))=gama(2)*gama(5),pd(sigma(1),sigma(3))=gama(3)*gama(5),pd(sigma(3),sigma(1))=2*I*gen(0,1),pd(sigma(3),sigma(2))=2*I*gen(0,2),pd(sigma(3),sigma(3))=2*I*gen(0,3),`x`);
Applying routine subs to dirac
> grdisplay(dirac);
> nops(grcomponent(dirac));
> gralter(dirac,normal,simplify);
Component simplification of a GRTensorII object:
Applying routine normal to object dirac
Applying routine simplify to object dirac
> grdisplay(dirac);
>
>