Poisson brackets for spinning post-Newtonian dynamics (in xTensor)

Author: Leo C. Stein (leo.stein@gmail.com)
Date: 2019 Apr. 29

Preamble

Needs["xAct`xTras`"];
$PrePrint=ScreenDollarIndices;
DefConstantSymbol/@{α,β,γ};
DefManifold[M,3,    Complement[IndexRange[a,z],{n,r,p}]  ];
DefMetric[1,δ[-a,-b],PD,FlatMetric->True]
PrintAs[epsilonδ]^="ϵ";
DefTensor[#[a],M]&/@{R,P ,L,J,n};

DefLabeledSpatialVector[v_,M_?ManifoldQ]:=Module[{nums,b=DummyIn[TangentBundleOfManifold[M]],d=DummyIn[TangentBundleOfManifold[M]]},
xTensorQ[v[___]]^:=True;
Dagger[v[l___]]^:=v[l];
DefInfo[v[___]]^:={"tensor",""};
DependenciesOfTensor[v[___]]^:={M};
HostsOf[v[___]]^:={M,TangentBundleOfManifold[M]};
PrintAs[v[l__]]^:=SubscriptBox[ToString[v],StringJoin@@ToString/@{l}];
PrintAs[v[]]^:=ToString[v];
SlotsOfTensor[v[___]]^:={TangentBundleOfManifold[M]};
SymmetryGroupOfTensor[v[___]]^:=StrongGenSet[{},GenSet[]];
TensorID[v[___]]^:={};
];
DefLabeledSpatialVector[#,M]&/@{S};

DefConstantSymbol[G];
DefConstantSymbol[ϵ];
DefConstantSymbol[m1,PrintAs->"(Subscript[m, 1])"];
DefConstantSymbol[m2,PrintAs->"(Subscript[m, 2])"];

Poisson bracket structure

PB[a_,b_Plus]:=PB[a,#]&/@b ;
PB[a_Plus,b_]:=PB[#,b]&/@a ;
PB[a_,b_ c_]:=PB[a,b]c+b PB[a,c];
PB[a_ b_, c_]:=PB[a,c]b+a PB[b,c];
PB[a_,func_?ScalarFunctionQ[args___]]:=Plus@@MapThread[(PB[a,#1] (Derivative@@#2)@func@@{args})&,{{args},IdentityMatrix@Length@{args}}];
PB[func_?ScalarFunctionQ[args___],a_]:=Plus@@MapThread[(PB[#1,a] (Derivative@@#2)@func@@{args})&,{{args},IdentityMatrix@Length@{args}}];
PB[a_,Scalar[b_]]:=PB[a,ReplaceDummies[b]];                                                                           
PB[Scalar[a_],b_]:=PB[ReplaceDummies[a],b];                                           
PB[_,_?ConstantQ]=0;                                                                                          
PB[_?ConstantQ,_]=0;
PB[_,_?PSIndepQ]=0;
PB[_?PSIndepQ,_]=0;
PSIndepQ[δ[__]]:=True;
PSIndepQ[Verbatim[delta][__]]:=True;
PSIndepQ[epsilonδ[__]]:=True;
PSVariables={R,P,S[1],S[2]};
PSVarQ[x_]:=MemberQ[PSVariables,x];
PSVarQ[x_[inds___]]:=MemberQ[PSVariables,x];
PB[a_?PSVarQ,b_?PSVarQ]:=PBHelper[a,b];
PBHelper[R[i_],P[j_]]:=δ[i,j];
PBHelper[P[j_],R[i_]]:=-δ[i,j];
PBHelper[S[A_][i_],S[B_][j_]]:=Module[{k},KroneckerDelta[A,B]epsilonδ[i,j,-k]S[A][k]];
PBHelper[_,_]=0;

More definitions and setting up of the Hamiltonian

L[i_]:=ReplaceDummies[epsilonδ[i,-j,-k]R[j]P[k]];
MM=m1+m2;             μ=m1 m2/MM;             ν=μ/MM;
IndexSet[p[i_],((P[i])/μ)];     IndexSet[r[i_],((R[i])/(G MM))];    Rm=Sqrt@Scalar[R[a]R[-a]];    rm=Sqrt@Scalar[r[a]r[-a]];          IndexSet[n[i_], PowerExpand[((r[i])/rm)]];
IndexSet[SS[i_],S[1][i]+S[2][i]];   
IndexSet[J[i_],SS[i]+L[i]];
 Q1=(2+(3/2) (m2/m1));    Q2=(2+(3/2)  (m1/m2)); 
HN=(((Scalar[p[i]p[-i]])/2)-(1/rm))μ;
H1PN=(ϵ((-1)/8)(1-3 ν)(((Scalar[p[i]p[-i]]))^2)-(ϵ/(2 ))(1/rm)((3+ν)(Scalar[p[i]p[-i]])+ ν (((Scalar[n[i]p[-i]]))^2))+ (ϵ/(2(rm^2)))) μ ;         
H2PN=(((ϵ^2)/16)(1-5 ν+5 (ν^2))(((Scalar[p[i]p[-i]]))^3)+((ϵ^2)/8)(1/rm)((5-20 ν-3 (ν^2))  (((Scalar[p[i]p[-i]]))^2)-2 (ν^2)(Scalar[p[i]p[-i]])(((Scalar[n[i]p[-i]]))^2)-3 (ν^2)(((Scalar[n[i]p[-i]]))^4)) + ((ϵ^2)/2)(1/(rm^2))((5+8 ν) (Scalar[p[i]p[-i]]) + 3 ν (((Scalar[n[i]p[-i]]))^2) ) -((ϵ^2)/4)(1/(rm^3))(1+3 ν)    ) μ     ;  
Seff[i_]:=Q1 S[1][i]+Q2 S[2][i] ;
HSO=((ϵ G)/(Rm^3))   ( Seff[i] L[-i]    //PutScalar     )     //ContractMetric //Simplification             ;                
HS1S1= ((ϵ G)/(2 (Rm^3)))(  3Scalar[(S[1][i]n[-i])]Scalar[(S[1][i]n[-i])] -Scalar[S[1][i]S[1][-i]]    )(m2/m1) ;
HS2S2=((ϵ G)/(2 (Rm^3)))    (  3Scalar[(S[2][i]n[-i])]Scalar[(S[2][i]n[-i])] -Scalar[S[2][i]S[2][-i]] )(m1/m2) ;             
HS1S2= ((ϵ G)/(Rm^3)) (3Scalar[(S[1][i]n[-i])]Scalar[(S[2][i]n[-i])] -  Scalar[S[1][i]S[2][-i]]  );            
HSS= HS1S1+HS2S2+HS1S2  ;
HO=HN+H1PN+H2PN;
HS=HSO+HSS;
H=HO+HS;

Computations begin

PBf[vvar1_,vvar2_]:=Module[{var1,var2,var3},
var1=  vvar1                        //ReplaceDummies        //ContractMetric //ToCanonical    ;
var2=vvar2                               //ReplaceDummies                      ;
var3= PB[var1,var2]          //ContractMetric  //PowerExpand//Simplification  //ReplaceDummies
];

var1 = S[1][b] L[-b]   ;
var2 = L[a] ;
var3= R[a];
var4=  P[i];

PBf[ var1, var2]                                  //ContractMetric//PutScalar  //Simplification  
PBf[ var1, var3]                                 //ContractMetric//PutScalar  //Simplification  
PBf[ var3, var4]                                 //ContractMetric//PutScalar  //Simplification