Computations associated to Article “Heisenberg-invariant self-dual Einstein manifolds”, by V. Cortés and Á. Murcia.

We divide all the computations into different sections. Please access the corresponding section to execute the code and check how it has been done

(Timelike) quaternionic (para)Kähler Heisenberg four-manifolds

(*I clear all previously defined variables and I make the necessary definitions to call the package "diffgeo", which we will extensively use*)
Clear["Global`*"];
coord={t,x,y,z};
X=coord;
metric=DiagonalMatrix[{ϵ,ϵ,1,1}];
<<diffgeo.m

n=4; (*dimensions*)

eucmetric=DiagonalMatrix[{1,1,1,1}];
VB[i_]:=Table[eucmetric[[i,k]],{k,1,n}];

(*This L will contain the Lie brackets of the time-evolving orthonormal frame*)
L=Normal@SymmetrizedArray[{i_,j_,k_}->l[i,j,k],{n,n,n},{Antisymmetric[{1,2}]}];
Do[Do[Do[l[i,j,k]=0,{i,1,n}],{j,1,n}],{k,1,n}];

(*I define the matrix Ut, which relates the frame at t=0 with the frame at t*)

Ut=Transpose[{{1,0,0,0},{0,a[t],p[t],g[t]},{0,v[t],b[t],h[t]},{0,w[t],f[t],c[t]}}];

(*As justified in the paper, we can set:*)

v[t_]:=0;w[t_]:=0;p[t_]:=0;

(*I write the time derivatives of the time-evolving frame*)

det=D[Ut,t].Inverse[Ut]//FullSimplify;

(*Then:*)

Do[Do[l[1,i,k]=det[[i,k]],{k,1,n}],{i,2,n}];

(*In order to write the Lie brackets between the elements of the frame, I define the Lie brackets at t=0*)

L0=Normal@SymmetrizedArray[{i_,j_,k_}->l0[i,j,k],{n,n,n},{Antisymmetric[{1,2}]}];
Do[Do[Do[l0[i,j,k]=0,{i,1,n}],{j,1,n}],{k,1,n}];
l0[3,4,2]=-2 q0;

(*At any time:*)

Do[Do[Do[l[i,j,k]=Sum[Sum[Sum[Ut[[i,r]]Ut[[j,ll]] L0[[r,ll,mm]]Inverse[Ut][[mm,k]]//Simplify,{r,2,n}],{ll,2,n}],{mm,2,n}],{i,2,n}],{j,2,n}],{k,2,n}]

(*I define the Lie bracket between two vectors*)

Lieb[v_,w_]:=Table[Sum[Sum[v[[i]]w[[j]] L[[i,j,k]],{j,1,n}],{i,1,n}]+v[[1]]D[w[[k]],t]-w[[1]] D[v[[k]],t],{k,1,n}];

(*Now we move into the computation of the covariant derivatives of the frame*)

imet=Inverse[metric];
Covdevs=Normal@SymmetrizedArray[{i_,j_,k_}->cd[i,j,k],{n,n,n}];
Do[Do[Do[cd[i,j,q]=1/2 Sum[Sum[L[[i,j,s]] metric[[s,k]] -L[[i,k,s]] metric[[j,s]]-L[[j,k,s]] metric[[i,s]],{s,1,n}]imet[[k,q]],{k,1,n}],{q,1,n}],{j,1,n}],{i,1,n}]

(*The covariant derivative operation is:*)

CD[v_,w_]:=Table[Sum[Sum[v[[i]]w[[j]] Covdevs[[i,j,k]],{j,1,n}],{i,1,n}]+v[[1]] D[w[[k]],t],{k,1,n}];

(*Now we can compute the Riemann tensor as follows*)

Do[Do[Do[Do[RR[i,j,k,q]=CD[VB[i],CD[VB[j],VB[k]]][[q]]-CD[VB[j],CD[VB[i],VB[k]]][[q]]-CD[Lieb[VB[i],VB[j]],VB[k]][[q]],{q,1,n}],{k,1,n}],{j,1,n}],{i,1,n}]//Simplify;

(*The Riemann tensor with components g(R(e_i,e_j)e_k,e_q)σ_q, where σ_q is the signature of the spacetime character of e_q, is*)

Riemcomp=Normal@SymmetrizedArray[{i_,j_,k_,q_}->RR[i,j,k,q],{n,n,n,n}];

(*With all indices down:*)

Do[Do[Do[Do[RRdown[i,j,k,q]=Sum[RR[i,j,k,ll] metric[[ll,q]],{ll,1,n}],{q,1,n}],{k,1,n}],{j,1,n}],{i,1,n}]//Simplify;
Riemcompdown=Normal@SymmetrizedArray[{i_,j_,k_,q_}->RRdown[i,j,k,q],{n,n,n,n}];

(*The Ricci tensor is therefore*)

Ricci=Table[Sum[Riemcomp[[k,i,j,k]],{k,1,n}],{i,1,n},{j,1,n}]//Simplify;

(*We observe*)

Ricci[[1,4]]

(*Then, as justified in the paper*)

f[t_]:=0;

(*Now...*)

Ricci[[1,3]]

(*So...*)

g[t_]:=0;

Ricci[[4,4]]//FullSimplify

(*Now I derive the Weyl tensor by decomposing the Riemann tensor in the common way*)

Ricci=Ricci/.{ϵ^2->1}//Simplify;
Scal=contract[Ricci]//Simplify;

Z=Ricci-1/4 Scal metric//Simplify;
S=Scal/6 antisymmetrize[swapIndices[metric**metric,{2,4}],{1,2}]//Simplify;

EE=2*antisymmetrize[swapIndices[Z**metric,{2,4}],{1,2},{3,4}]//Simplify;

(*The Weyl tensor is*)

W=((Riemcompdown-EE-S))/.{ϵ^2->1}//Simplify;
(*Consistency check*)

contract[W,{1,4}]/.{ϵ^2->1}//Simplify

(*We compute the dual Weyl tensor*)

eps=LeviCivita/.{ϵ^2->1,Sign[ϵ]^2->1};
starW=1/2 contract[eps**W,{3,5},{4,6}]//Simplify;
(*The self-duality tensor is*)

selfduality=(starW- W)/.{ϵ^2->1}//Simplify;

(*Now we proceed to complete the classification*)

hpp=Solve[Ricci[[3,4]]==0,h''[t]]//Simplify

selfduality=selfduality/.hpp[[1]]//Simplify;
selfduality[[2,3,2,4]]

(*Assume q0 bc+a'=0. Then:*)

Ricci[[2,2]]/.{a'[t]->-q0 b[t]c[t],a''[t]->-q0(b'[t]c[t]+b[t]c'[t])}//Simplify

(*The Einstein constant should vanish in that case, so...*)

h[t_]:=0;

(*Now:*)

Ricci=Ricci//Simplify;
solRicci=Solve[Ricci==Λ metric,{a'[t],a''[t],b''[t],c''[t]}]//Simplify

selfduality=selfduality/.solRicci[[1]]//Simplify;
Wsim=W/.solRicci[[1]]//Simplify;
selfduality[[2,3,2,3]]//Simplify

lambdasol=Solve[selfduality[[2,3,2,3]]==0,Λ]//Simplify

selfduality=selfduality/.lambdasol[[1]]//Simplify;
Wsim=Wsim/.lambdasol[[1]]//Simplify;
selfduality[[1,2,1,2]]//Simplify

(*If the first or second terms vanish...*)

(Λ/.lambdasol[[1]])/.{a[t]->q0 b[t]b[t] c[t]/b'[t]}//Simplify

(Λ/.lambdasol[[1]])/.{a[t]->q0 c[t]c[t] b[t]/c'[t]}//Simplify

(*So necessarily*)

c[t_]:=b[t];

(*Y el tensor de Weyl es autodual:*)

(selfduality)/.{ϵ^2->1}//Simplify

lambdaexpfinal=(Λ/.lambdasol[[1]])//Simplify

(*The expression of a' is...*)

a'[t]/.solRicci[[1]]//Simplify

Weyl tensor and its eigenvectors

(*The Weyl tensor is:*)
((Wsim-1/2 contract[eps**Wsim,{3,5},{4,6}])//Simplify/.{ϵ^2->1})/.{ϵ^2->1}
(*I want to diagonalize the Weyl tensor. I define the following basis of self-dual 2-forms*)

bb1=((FormToTensor[Wedge[dd[t],dd[x]]+HodgeStar[Wedge[dd[t],dd[x]]]])//Simplify/.{ϵ^2->1,(Sign[ϵ])^2->1})/.{ϵ^2->1,(Sign[ϵ])^2->1,1/ϵ^2->1}

bb2=((FormToTensor[Wedge[dd[t],dd[y]]+HodgeStar[Wedge[dd[t],dd[y]]]])//Simplify/.{ϵ^2->1,(Sign[ϵ])^2->1})/.{ϵ^2->1,(Sign[ϵ])^2->1,1/ϵ^2->1}

bb3=((FormToTensor[Wedge[dd[t],dd[z]]+HodgeStar[Wedge[dd[t],dd[z]]]])//Simplify/.{ϵ^2->1,(Sign[ϵ])^2->1})/.{ϵ^2->1,(Sign[ϵ])^2->1,1/ϵ^2->1}

bb={bb1,bb2,bb3};

met2form=((Table[contract[bb[[i]]**bb[[j]],{1,3},{2,4}],{i,1,3},{j,1,3}])/.{1/ϵ^2->1,1/ϵ^3->1/ϵ}//Simplify)/.{1/ϵ->ϵ}

(*The (local) matrix elements of the Weyl tensor on the bundle of self-dual 2-forms is*)

W2=(contract[Wsim**Wsim,{1,5},{2,6},{3,7},{4,8}])/.{ϵ^2->1,ϵ^4->1}//Simplify

Wmat=Table[(met2form[[j,j]])^(-1)(contract[contract[Wsim**bb[[i]],{3,5},{4,6}]**bb[[j]],{1,3},{2,4}])/.{ϵ^2->1}//Simplify,{i,1,3},{j,1,3}]

(*This is already the diagonal basis!*)
(*By direct inspection, we observe that Wmat already has two equal eigenvalues*)

Wmat/.{ϵ^4->1}

(*This implies our space is (locally) conformally Kähler. In fact:*)

ωk1=bb1*conf1[t];
ωk2=bb2*conf2[t];
ωk3=bb3*conf2[t];

W2

extder[ω_]:=Table[D[ω[[j,k]],X[[i]]]-D[ω[[i,k]],X[[j]]]+D[ω[[i,j]],X[[k]]]+Sum[-ω[[s,k]]Lieb[VB[i],VB[j]][[s]]+ω[[s,j]] Lieb[VB[i],VB[k]][[s]]-ω[[s,i]] Lieb[VB[j],VB[k]][[s]],{s,1,n}],{i,1,n},{j,1,n},{k,1,n}]

dωk1=extder[ωk1]//Simplify
dωk2=extder[ωk2]//Simplify
dωk3=extder[ωk3]//Simplify

conf1[t_]:=(((1536 (q0^6) ((b[t])^14))/((ϵ^2) ((a[t])^4) (((-3 q0 ((b[t])^3)+a[t] b'[t]))^2))))^(1/3);

((dωk1/.solRicci[[1]])/.lambdasol[[1]])//Simplify

(contract[ωk1**ωk1,{1,3},{2,4}]*(W2)^(-2/3))/.{ϵ^2->1,1/ϵ^2->1}//Simplify

(*Then, ωk1 is Kähler with respect to the rescaled metric \bar{g}=W2^(1/3)g. Interestingly enough: *)

conf2[t_]:=a[t] b[t];

dωk2//Simplify

dωk3//Simplify

contract[ωk2**ωk2,{1,3},{2,4}]*1/(a[t]b[t])^2/.{ϵ^2->1,1/ϵ^3->1/ϵ}//Simplify
contract[ωk3**ωk3,{1,3},{2,4}]*1/(a[t]b[t])^2/.{ϵ^2->1,1/ϵ^3->1/ϵ}//Simplify

(*We check that ω1 is a Kähler, while ω2 and ω3 are almost Kähler*)

J1=(raise[bb1,{1}]//Simplify)/.{1/ϵ^2->1}
J2=(raise[bb2,{1}]//Simplify)/.{1/ϵ^2->1}
J3=(raise[bb3,{1}]//Simplify)/.{1/ϵ^2->1}

J1.J1
J2.J2
J3.J3

NJ1=(Table[-(J1.J1).Lieb[VB[i],VB[j]]+J1.(Lieb[J1.VB[i],VB[j]]+Lieb[VB[i],J1.VB[j]])-Lieb[J1.VB[i],J1.VB[j]],{i,1,n},{j,1,n}]//Simplify)/.{ϵ^2->1}

NJ2=(((Table[-(J2.J2).(Lieb[VB[i],VB[j]])+J2.(Lieb[J2.VB[i],VB[j]]+Lieb[VB[i],J2.VB[j]])-Lieb[J2.VB[i],J2.VB[j]],{i,1,n},{j,1,n}]//Simplify)/.solRicci[[1]])/.lambdasol[[1]])//Simplify/.{ϵ^2->1}

NJ3=(Table[-(J3.J3).(Lieb[VB[i],VB[j]])+J3.(Lieb[J3.VB[i],VB[j]]+Lieb[VB[i],J3.VB[j]])-Lieb[J3.VB[i],J3.VB[j]],{i,1,n},{j,1,n}]//Simplify)/.{ϵ^2->1}

Non stationary solutions

(*We first observe:*)

Bl=1/(3a[t]^3)*(3q0 a[t]^2 b[t]^3-(ϵ Λ a[t]^6 b[t]^2 Exp[-2i π ll/3])/(9 q0 ϵ Λ a[t]^8 b[t]^5+Sqrt[Λ^2 a[t]^(16) b[t]^6*(ϵ Λ a[t]^2+81 q0^2 b[t]^4)])^(1/3)+Exp[2i π ll/3]*(9 q0 ϵ Λ a[t]^8 b[t]^5+Sqrt[Λ^2 a[t]^(16) b[t]^6*(ϵ Λ a[t]^2+81 q0^2 b[t]^4)])^(1/3))//Simplify;

((((3 (((q0 ((b[t])^3)-a[t] b'[t]))^3))/(ϵ ((a[t])^2) ((b[t])^2) (-3 q0 ((b[t])^3)+a[t] b'[t])))-Λ)/.{b'[t]->Bl}//Simplify)/.{ll->0,ϵ^2->1}//Simplify

(((((3 (((q0 ((b[t])^3)-a[t] b'[t]))^3))/(ϵ ((a[t])^2) ((b[t])^2) (-3 q0 ((b[t])^3)+a[t] b'[t])))-Λ)/.{b'[t]->Bl}//Simplify)/.{ll->1,ϵ^2->1}//Simplify)/.{(E^(2 i π))->1}

(((((3 (((q0 ((b[t])^3)-a[t] b'[t]))^3))/(ϵ ((a[t])^2) ((b[t])^2) (-3 q0 ((b[t])^3)+a[t] b'[t])))-Λ)/.{b'[t]->Bl}//Simplify)/.{ll->2,ϵ^2->1}//Simplify)/.{(E^(4 i π))->1}

(*So these are the solutions to the cubic equation for Λ*)

Clear[a,b]

(*Negative (timelike) quaternionic (para)Kähler Heisenberg four-manifolds*)
a[t_]:=s q0 Sqrt[-2 ϵ Λ/3] ρ[t]Sqrt[(ρ[t]+2γ)/(ρ[t]+γ)];
b[t_]:=s Sqrt[-ϵ Λ/3] ρ[t]/Sqrt[ρ[t]+2γ];

rhoder=Sqrt[-2 ϵ Λ/3] ρ[t]*Sqrt[(ρ[t]+γ)/(ρ[t]+2γ)];
rhoder2=D[rhoder,t]/.{ρ'[t]->rhoder}//Simplify

((Ricci/.{ρ'[t]->rhoder,ρ''[t]->rhoder2})/.{s^2->1,s^4->1}//Simplify)/.{s^2->1,s^4->1}//Simplify/.{1/s^2->1}

(Wsim-1/2 contract[eps**Wsim,{3,5},{4,6}]//Simplify)/.{ϵ^2->1}

(*Let us explore the issue of the initial conditions*)

γsol=Solve[(b[t]/.{t->0})==1,γ]//Simplify

conda=((a[t]/.{t->0})/.γsol[[1]]//Simplify)/.{s^2->1}

eqconda=(9-3 ϵ Λ ρ[0])*(conda^2-1)/.{s^2->1}//FullSimplify

rhol=-Exp[-(4-2ll)I Pi/3]/(2q0 (ϵ Λ)^(1/3) (9 q0+Sqrt[81 q0^2+ ϵ Λ])^(1/3))+Exp[(4-2ll)I Pi/3](9 q0+Sqrt[81 q0^2+ ϵ Λ])^(1/3)/(2 q0 (ϵ Λ)^(2/3));

(eqconda/.{ρ[0]->rhol}//Simplify)/.{ll->1}//Simplify

(eqconda/.{ρ[0]->rhol}//Simplify)/.{ll->2}//Simplify

(eqconda/.{ρ[0]->rhol}//Simplify)/.{ll->3}//Simplify

(*By direct inspection... We proceed numerically, because Mathematica has some difficulties with fractional powers. In order to grasp the qualitative behaviour*)

N[rhol/.{q0->1,ϵ->1,s->1,Λ->-95,ll->3}]

N[((((b[t]/.{t->0})/.γsol[[1]])/.{ρ[0]->rhol}//Simplify)/.{ll->1}//Simplify)/.{q0->1,ϵ->-1,s->1,Λ->80}//Simplify]

N[((γ/.γsol[[1]])/.{ρ[0]->rhol})/.{q0->1,ϵ->-1,s->1,Λ->80,ll->1}]

(*So, if s=1, three real solutions for ϵΛ<-81q0^2 and one real solution (ll=1) if 0> ϵΛ>-81 q0^2*)

(*By some numerical results, we have that, for s=-1:*)

N[rhol/.{q0->1,ϵ->1,s->-1,Λ->-97,ll->3}]

Plot[(rhol/.{q0->1,ϵ->1,s->-1,ll->3}),{Λ,-1000,-81}]

N[5/80]

N[(((b[t]/.{t->0})/.γsol[[1]])/.{ρ[0]->rhol})/.{q0->1,ϵ->1,s->-1,Λ->-80,ll->1}//Simplify]

(*All the necessary graphs and studies can be done to verify what the paper claims*)

(*Positive (timelike) quaternionic (para)Kähler Heisenberg four-manifolds*)

a[t_]:= q0 Sqrt[2 ϵ Λ/3] ρ[t]Sqrt[-(ρ[t]+2γ)/(ρ[t]+γ)];
b[t_]:= Sqrt[ϵ Λ/3] ρ[t]/Sqrt[-ρ[t]-2γ];

rhoder=Sqrt[2 ϵ Λ/3] ρ[t]*Sqrt[-(ρ[t]+γ)/(ρ[t]+2γ)];
rhoder2=D[rhoder,t]/.{ρ'[t]->rhoder}//Simplify

(Ricci/.{ρ'[t]->rhoder,ρ''[t]->rhoder2})//Simplify

(Wsim-1/2 contract[eps**Wsim,{3,5},{4,6}]//Simplify)/.{ϵ^2->1}

Some geodesics

e1t0=FormToTensor[dd[z]+q0 x dd[y]-q0  y dd[x]];
e2t0={0,1,0,0};
e3t0={0,0,1,0};
metgeo=ϵ{1,0,0,0}**{1,0,0,0}+ϵ(a[t])^(-2)e1t0**e1t0+(b[t])^(-2)(e2t0**e2t0+e3t0**e3t0)//Simplify

Clear["Global`*"];
coord={ρ,x,y,z};
metric=({{ 3/(2*Sqrt[Λ^2]*ρ^2)*(ρ+2γ)/(ρ+γ)ϵ,0,0,0},{0,(((q0^2) (y^2) ϵ)/((a[t])^2))+(1/((b[t])^2)),-(((q0^2) x y ϵ)/((a[t])^2)),-((q0 y ϵ)/((a[t])^2))},{0,-(((q0^2) x y ϵ)/((a[t])^2)),(((q0^2) (x^2) ϵ)/((a[t])^2))+(1/((b[t])^2)),((q0 x ϵ)/((a[t])^2))},{0,-((q0 y ϵ)/((a[t])^2)),((q0 x ϵ)/((a[t])^2)),(ϵ/((a[t])^2))}})/.{t->ρ};
<<diffgeo.m

c[t_]:=b[t];

((RicciTensor- Λ metric)/.{a'[t]->-(((q0^2) ((b[t])^3) ((c[t])^3)+((a[t])^2) (ϵ Λ b[t] c[t]+b'[t] c'[t]))/(a[t] (c[t] b'[t]+b[t] c'[t]))),a''[t]->(2 (q0^4) ((b[t])^7) ((c[t])^7)+(q0^2) ((a[t])^2) ((b[t])^3) ((c[t])^3) (-3 ((c[t])^2) ((b'[t])^2)-2 b[t] c[t] b'[t] c'[t]+((b[t])^2) (4 ϵ Λ ((c[t])^2)-3 ((c'[t])^2)))+((a[t])^4) (2 (ϵ^2) (Λ^2) ((b[t])^3) ((c[t])^3)-((c[t])^2) ((b'[t])^3) c'[t]-((b[t])^2) b'[t] c'[t] (-4 ϵ Λ ((c[t])^2)+((c'[t])^2))))/(((a[t])^3) b[t] c[t] (((c[t] b'[t]+b[t] c'[t]))^2)),b''[t]->((q0^2) ((b[t])^4) ((c[t])^4) b'[t]+2 ((a[t])^2) ((c[t])^2) ((b'[t])^3)+ϵ Λ ((a[t])^2) ((b[t])^3) c[t] c'[t]+2 (q0^2) ((b[t])^5) ((c[t])^3) c'[t]+2 ((a[t])^2) b[t] c[t] ((b'[t])^2) c'[t]+((a[t])^2) ((b[t])^2) b'[t] ((c'[t])^2))/(((a[t])^2) b[t] c[t] (c[t] b'[t]+b[t] c'[t])),c''[t]->((q0^2) ((b[t])^3) ((c[t])^4) (2 c[t] b'[t]+b[t] c'[t])+((a[t])^2) (((c[t])^2) ((b'[t])^2) c'[t]+2 ((b[t])^2) ((c'[t])^3)+b[t] c[t] b'[t] (ϵ Λ ((c[t])^2)+2 ((c'[t])^2))))/(((a[t])^2) b[t] c[t] (c[t] b'[t]+b[t] c'[t]))}//Simplify)/.{ϵ^2->1};

(*I want to derive the geodesic equation, assuming a curve is affinely parametrized*)

(*Negative (timelike) quaternionic (para)Kähler Heisenberg four-manifolds*)

curv={ρρ[τ],xx[τ],yy[τ],zz[τ]};
Vec=D[curv,τ]//Simplify;
geodeq=(D[curv,{τ,2}]/.{x->xx[τ],y->yy[τ],z->zz[τ],ρ->ρρ[τ]})+contract[((Christoffel)/.{x->xx[τ],y->yy[τ],z->zz[τ],ρ->ρρ[τ]})**lower[Vec]**lower[Vec],{2,4},{3,5}]//Simplify

a[ρ_]:= q0 Sqrt[-2 ϵ Λ/3] ρ Sqrt[(ρ+2γ)/(ρ+γ)];
b[ρ_]:= Sqrt[-ϵ Λ/3] ρ/Sqrt[ρ+2γ];

(*Trials*)

Clear[xx,yy,ρρ,zz,q0,Λ]

zz[τ_]:=0;yy[τ_]:=0;

xx[τ_]:=vx Integrate[ ((( ((ρρ[σ])^2))/(2 γ+ρρ[σ]))),{σ,0,τ}];

ϵ=-1;

geodeqsim=Assuming[Λ>0,geodeq[[1]]//FullSimplify]

rhodos=Solve[geodeqsim==0,ρρ''[τ]]//FullSimplify

Clear[ϵ,q0,Λ]

(*Now I provide three examples of numerical resolutions which provide a strong argument in favour of geodesic incompleteness *)
Clear[ϵ,q0,Λ]
ϵ=-1;
q0=1;
Λ=7;
vx=1;
nsol=NSolve[{((a[ρ])/.{ρ->ρ0})==1,((b[ρ])/.{ρ->ρ0})==1},{ρ0,γ}]

nsoll=NDSolve[{ρρ''[τ]==((ρρ''[τ]/.rhodos[[1]])/.nsol[[1]]),ρρ[0]==(ρ0/.nsol[[1]]),ρρ'[0]==1},ρρ[τ],{τ,0,2}]

Plot[(ρρ[τ]/.nsoll[[1]]),{τ,0,1.5}]

Clear[ϵ,q0,Λ]
ϵ=-1;
q0=1;
Λ=20;
vx=1;
nsol=NSolve[{((a[ρ])/.{ρ->ρ0})==1,((b[ρ])/.{ρ->ρ0})==1},{ρ0,γ}]

nsoll=NDSolve[{ρρ''[τ]==((ρρ''[τ]/.rhodos[[1]])/.nsol[[1]]),ρρ[0]==(ρ0/.nsol[[1]]),ρρ'[0]==1},ρρ[τ],{τ,0,2}]

Plot[(ρρ[τ]/.nsoll[[1]]),{τ,0,1.2}]

Clear[ϵ,q0,Λ]
ϵ=-1;
q0=1;
Λ=12;
vx=-2;
nsol=NSolve[{((a[ρ])/.{ρ->ρ0})==1,((b[ρ])/.{ρ->ρ0})==1},{ρ0,γ}]

nsoll=NDSolve[{ρρ''[τ]==((ρρ''[τ]/.rhodos[[1]])/.nsol[[1]]),ρρ[0]==(ρ0/.nsol[[1]]),ρρ'[0]==1},ρρ[τ],{τ,0,2}]

Plot[(ρρ[τ]/.nsoll[[1]]),{τ,0,1}](*It diverges*)

Stationary solutions

(*This analysys can be easili carried out from the Ricci tensor:*)

Ricstat=Ricci/.{a'[t]->a[t] μ[t], a''[t]->a[t] (μ[t])^2+a[t] μ'[t], b'[t]->b[t]λ[t], b''[t]->b[t](λ[t])^2 +b[t]λ'[t]}//Simplify

(*A true stationary solution requires...*)

Ricstat=(Ricstat/.{a[t]->(b[t])^2,μ[t]->2 λ[t]})/.{μ[t]->μ,λ[t]->λ,μ'[t]->0,λ'[t]->0}//Simplify

(*Then:*)

Solve[Ricstat==Λ metric,{λ,Λ}]//Simplify

(*Consequently, picking the minus sign for convenience:*)
a[t_]:=Exp[-2q0 (t-t0)];b[t_]:=Exp[-q0 (t-t0)];

Ricci//Simplify

(*I will need later the Riemann tensor*)

((Riemcomp)//Simplify/.{ϵ^2->1})/.{ϵ^2->1}

(*The remaining task is to prove that it is (locally if ϵ=-1, globally if ϵ=1) isometric to a left-invariant metric on the simply transitive solvable Iwasawa subgroup of SU(1,2) and SU(2,1), respectively if ϵ=1 or ϵ=-1*)

(*We first check that it is a locally symmetric space by computing the covariant derivative of the Riemann tensor and observing that it vanishes*)

Clear[a,b,nablaR]

nablaRiem=Normal@SymmetrizedArray[{i_,j_,k_,q_,s_}->nablaR[i,j,k,q,s],{n,n,n,n,n}];

Do[Do[Do[Do[Do[nablaR[i,j,k,q,s]=eucmetric[[i,1]] D[RRdown[j,k,q,s],t]-Sum[CD[VB[i],VB[j]][[ll]] RRdown[ll,k,q,s],{ll,1,n}]-
Sum[CD[VB[i],VB[k]][[ll]] RRdown[j,ll,q,s],{ll,1,n}]-
Sum[CD[VB[i],VB[q]][[ll]] RRdown[j,k,ll,s],{ll,1,n}]-
Sum[CD[VB[i],VB[s]][[ll]] RRdown[j,k,q,ll],{ll,1,n}],{i,1,n}],{j,1,n}],{k,1,n}],{q,1,n}],{s,1,n}]

a[t_]:=Exp[-2q0 (t-t0)];b[t_]:=Exp[-q0 (t-t0)];

((nablaRiem)//Simplify/.{ϵ^2->1})/.{ϵ^2->1}

(*We observe it is covariantly constant*)

(*Now we want to identify the local (isometry type)*)

Riemannian case

(*Lie algebra of SU(2,1) embedded in GL(6,R)*)

Clear["Global`*"];
e1={{0,-1,0,0,0,0},{1,0,0,0,0,0},{0,0,0,0,0,0},{0,0,0,0,0,0},{0,0,0,0,0,1},{0,0,0,0,-1,0}};
e2={{0,0,0,0,0,0},{0,0,0,0,0,0},{0,0,0,1,0,0},{0,0,-1,0,0,0},{0,0,0,0,0,-1},{0,0,0,0,1,0}};
e3={{0,0,0,0,0,0},{0,0,0,0,0,0},{0,0,0,0,1,0},{0,0,0,0,0,1},{0,0,-1,0,0,0},{0,0,0,-1,0,0}};
e4={{0,0,0,0,0,0},{0,0,0,0,0,0},{0,0,0,0,0,1},{0,0,0,0,-1,0},{0,0,0,1,0,0},{0,0,-1,0,0,0}};
e5={{0,0,0,0,1,0},{0,0,0,0,0,1},{0,0,0,0,0,0},{0,0,0,0,0,0},{1,0,0,0,0,0},{0,1,0,0,0,0}};
e6={{0,0,0,0,0,1},{0,0,0,0,-1,0},{0,0,0,0,0,0},{0,0,0,0,0,0},{0,-1,0,0,0,0},{1,0,0,0,0,0}};
e7={{0,0,1,0,0,0},{0,0,0,1,0,0},{1,0,0,0,0,0},{0,1,0,0,0,0},{0,0,0,0,0,0},{0,0,0,0,0,0}};
e8={{0,0,0,1,0,0},{0,0,-1,0,0,0},{0,-1,0,0,0,0},{1,0,0,0,0,0},{0,0,0,0,0,0},{0,0,0,0,0,0}};
Ve1=1/Sqrt[2]*(e1+ e2);
Ve2=1/Sqrt[6]*(e1-e2);
Ve=μ/Sqrt[4]{Sqrt[2] Ve1, Sqrt[2] Ve2,e3,e4,e5,e6,e7,e8};

n=8;
KForm[A_,B_]:=1/μ^2 Tr[A.B];
metric=Table[KForm[Ve[[i]],Ve[[j]]],{i,1,n},{j,1,n}]

L=Normal@SymmetrizedArray[{i_,j_,k_}->l[i,j,k],{n,n,n},{Antisymmetric[{1,2}]}];
Do[Do[Do[l[i,j,k]=0,{i,1,n}],{j,1,n}],{k,1,n}]

Do[Do[Do[l[i,j,k]=Inverse[metric][[k,k]] KForm[(Ve[[i]].Ve[[j]]-Ve[[j]].Ve[[i]]),Ve[[k]]],{k,1,n}],{i,1,n}],{j,1,n}]

(*I notice:*)

Table[L[[i,j,k]],{i,1,4},{j,5,n},{k,1,4}]//Simplify

(*So SU(2,1)/S(U(1)xU(2)) is, at the Lie algebra level at least, a reductive space. This means that we can compute its Riemann in the usally easier form:*)

lres=Table[L[[i,j,k]],{i,5,n},{j,5,n},{k,1,4}]//Simplify;

Riemred=-Table[Sum[lres[[i,j,k]] Inverse[metric][[s,s]] KForm[(Ve[[k]].Ve[[ll]]-Ve[[ll]].Ve[[k]]),Ve[[s]]],{k,1,4}],{i,1,4},{j,1,4},{ll,5,8},{s,5,8}]//Simplify;

(*Now want to see if this Riemann tensor, after simple (orthogonal) transformations, can be set to be equal to the components of the Riemann curvature of our (Riemannian) stationary solution*)

(*We copy what we found before, for ϵ=1*)

Riemcomp={{{{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}},{{0,4 (q0^2),0,0},{-4 (q0^2),0,0,0},{0,0,0,2 (q0^2) ϵ},{0,0,-2 (q0^2) ϵ,0}},{{0,0,(q0^2),0},{0,0,0,(q0^2) ϵ},{-((q0^2)/ϵ),0,0,0},{0,-(q0^2),0,0}},{{0,0,0,(q0^2)},{0,0,-(q0^2) ϵ,0},{0,(q0^2),0,0},{-((q0^2)/ϵ),0,0,0}}},{{{0,-4 (q0^2),0,0},{4 (q0^2),0,0,0},{0,0,0,-2 (q0^2) ϵ},{0,0,2 (q0^2) ϵ,0}},{{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}},{{0,0,0,-(q0^2) ϵ},{0,0,(q0^2),0},{0,-((q0^2)/ϵ),0,0},{(q0^2),0,0,0}},{{0,0,(q0^2) ϵ,0},{0,0,0,(q0^2)},{-(q0^2),0,0,0},{0,-((q0^2)/ϵ),0,0}}},{{{0,0,-(q0^2),0},{0,0,0,-(q0^2) ϵ},{((q0^2)/ϵ),0,0,0},{0,(q0^2),0,0}},{{0,0,0,(q0^2) ϵ},{0,0,-(q0^2),0},{0,((q0^2)/ϵ),0,0},{-(q0^2),0,0,0}},{{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}},{{0,2 (q0^2),0,0},{-2 (q0^2),0,0,0},{0,0,0,((4 (q0^2))/ϵ)},{0,0,-((4 (q0^2))/ϵ),0}}},{{{0,0,0,-(q0^2)},{0,0,(q0^2) ϵ,0},{0,-(q0^2),0,0},{((q0^2)/ϵ),0,0,0}},{{0,0,-(q0^2) ϵ,0},{0,0,0,-(q0^2)},{(q0^2),0,0,0},{0,((q0^2)/ϵ),0,0}},{{0,-2 (q0^2),0,0},{2 (q0^2),0,0,0},{0,0,0,-((4 (q0^2))/ϵ)},{0,0,((4 (q0^2))/ϵ),0}},{{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}}}}/.{ϵ->1};

(*We must choose the parameter μ appropiately. By observing*)

Ricred=Table[Sum[Riemred[[i,j,k,i]],{i,1,4}],{j,1,4},{k,1,4}]//Simplify

(*Therefore...*)

μ=2 q0;

(*And now I trivially notice*)

Riemcomp- Riemred//Simplify

(*Now we have to find the Heisenberg algebra within su(2,1))*)
(*By direct inspection*)

(*We observe*)

h1=Ve[[3]]+Ve[[5]];
h2=Ve[[4]]+Ve[[6]];
h3=-(Ve[[1]]+Ve[[8]]);

(h1.h2-h2.h1+a h3)/.{a->2q0}//Simplify

h2.h3-h3.h2

h1.h3-h3.h1//Simplify

(*We have found the Heisenberg group  within the su(2,1) algebra*)
(*Since we are in Riemannian signature, and noticing that the intersection between this embedded Heisenberg group with the algebra of the isotropy group at the origin has empty intersection, the conclusions of our work follow*)

Neutral-signature case

(*Six dimensional representation of su(1,2)*)

Clear["Global`*"];
e1={{0,1,0,0,0,0},{-1,0,0,0,0,0},{0,0,0,-1,0,0},{0,0,1,0,0,0},{0,0,0,0,0,0},{0,0,0,0,0,0}};
e2={{0,-1,0,0,0,0},{1,0,0,0,0,0},{0,0,0,0,0,0},{0,0,0,0,0,0},{0,0,0,0,0,1},{0,0,0,0,-1,0}};
e3={{0,0,0,0,-1,0},{0,0,0,0,0,-1},{0,0,0,0,0,0},{0,0,0,0,0,0},{-1,0,0,0,0,0},{0,-1,0,0,0,0}};
e4={{0,0,0,0,0,1},{0,0,0,0,-1,0},{0,0,0,0,0,0},{0,0,0,0,0,0},{0,-1,0,0,0,0},{1,0,0,0,0,0}};
e5={{0,0,1,0,0,0},{0,0,0,1,0,0},{-1,0,0,0,0,0},{0,-1,0,0,0,0},{0,0,0,0,0,0},{0,0,0,0,0,0}};
e6={{0,0,0,-1,0,0},{0,0,1,0,0,0},{0,-1,0,0,0,0},{1,0,0,0,0,0},{0,0,0,0,0,0},{0,0,0,0,0,0}};
e7={{0,0,0,0,0,0},{0,0,0,0,0,0},{0,0,0,0,1,0},{0,0,0,0,0,1},{0,0,1,0,0,0},{0,0,0,1,0,0}};
e8={{0,0,0,0,0,0},{0,0,0,0,0,0},{0,0,0,0,0,1},{0,0,0,0,-1,0},{0,0,0,-1,0,0},{0,0,1,0,0,0}};
Ve1=1/Sqrt[2]*(e1+ e2);
Ve2=1/Sqrt[6]*(e1-e2);
Ve=μ/Sqrt[4]{Sqrt[2] Ve1, Sqrt[2] Ve2,e3,e4,e5,e6,e7,e8};

n=8;
KForm[A_,B_]:=-1/μ^2 Tr[A.B];
metric=Table[KForm[Ve[[i]],Ve[[j]]],{i,1,n},{j,1,n}]

L=Normal@SymmetrizedArray[{i_,j_,k_}->l[i,j,k],{n,n,n},{Antisymmetric[{1,2}]}];
Do[Do[Do[l[i,j,k]=0,{i,1,n}],{j,1,n}],{k,1,n}]

Do[Do[Do[l[i,j,k]=Inverse[metric][[k,k]] KForm[(Ve[[i]].Ve[[j]]-Ve[[j]].Ve[[i]]),Ve[[k]]],{k,1,n}],{i,1,n}],{j,1,n}]

(*I notice:*)

Table[L[[i,j,k]],{i,1,4},{j,5,n},{k,1,4}]//Simplify

(*So SU(2,1)/S(U(1)xU(2)) is, at the Lie algebra level at least, a reductive space. This means that we can compute its Riemann in the usally easier form:*)

lres=Table[L[[i,j,k]],{i,5,n},{j,5,n},{k,1,4}]//Simplify;

Riemred=-Table[Sum[lres[[i,j,k]] Inverse[metric][[s,s]] KForm[(Ve[[k]].Ve[[ll]]-Ve[[ll]].Ve[[k]]),Ve[[s]]],{k,1,4}],{i,1,4},{j,1,4},{ll,5,8},{s,5,8}]//Simplify;

(*Now want to see if this Riemann tensor, after simple (orthogonal) transformations, can be set to be equal to the components of the Riemann curvature of our (Riemannian) stationary solution*)

(*We copy what we found before, for ϵ=1*)

Riemcomp={{{{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}},{{0,4 (q0^2),0,0},{-4 (q0^2),0,0,0},{0,0,0,2 (q0^2) ϵ},{0,0,-2 (q0^2) ϵ,0}},{{0,0,(q0^2),0},{0,0,0,(q0^2) ϵ},{-((q0^2)/ϵ),0,0,0},{0,-(q0^2),0,0}},{{0,0,0,(q0^2)},{0,0,-(q0^2) ϵ,0},{0,(q0^2),0,0},{-((q0^2)/ϵ),0,0,0}}},{{{0,-4 (q0^2),0,0},{4 (q0^2),0,0,0},{0,0,0,-2 (q0^2) ϵ},{0,0,2 (q0^2) ϵ,0}},{{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}},{{0,0,0,-(q0^2) ϵ},{0,0,(q0^2),0},{0,-((q0^2)/ϵ),0,0},{(q0^2),0,0,0}},{{0,0,(q0^2) ϵ,0},{0,0,0,(q0^2)},{-(q0^2),0,0,0},{0,-((q0^2)/ϵ),0,0}}},{{{0,0,-(q0^2),0},{0,0,0,-(q0^2) ϵ},{((q0^2)/ϵ),0,0,0},{0,(q0^2),0,0}},{{0,0,0,(q0^2) ϵ},{0,0,-(q0^2),0},{0,((q0^2)/ϵ),0,0},{-(q0^2),0,0,0}},{{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}},{{0,2 (q0^2),0,0},{-2 (q0^2),0,0,0},{0,0,0,((4 (q0^2))/ϵ)},{0,0,-((4 (q0^2))/ϵ),0}}},{{{0,0,0,-(q0^2)},{0,0,(q0^2) ϵ,0},{0,-(q0^2),0,0},{((q0^2)/ϵ),0,0,0}},{{0,0,-(q0^2) ϵ,0},{0,0,0,-(q0^2)},{(q0^2),0,0,0},{0,((q0^2)/ϵ),0,0}},{{0,-2 (q0^2),0,0},{2 (q0^2),0,0,0},{0,0,0,-((4 (q0^2))/ϵ)},{0,0,((4 (q0^2))/ϵ),0}},{{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}}}}/.{ϵ->-1};

(*We must choose the parameter μ appropiately. By observing*)

Ricred=Table[Sum[Riemred[[i,j,k,i]],{i,1,4}],{j,1,4},{k,1,4}]//Simplify

(*Therefore...*)

μ=2 q0;

(*And now I trivially notice that after performing the following permutation of the basis elements*)

ord={3,4,2,1}

Riemcomp- Riemred[[ord,ord,ord,ord]]//Simplify

(*Now we have to find the Heisenberg algebra within su(2,1))*)
(*By direct inspection*)

(*We observe*)

h1=Ve[[3]]+Ve[[5]];
h2=Ve[[4]]+Ve[[6]];
h3=(Ve[[1]]+Ve[[8]]);

(h1.h2-h2.h1+a h3)/.{a->2q0}//Simplify

h2.h3-h3.h2

h1.h3-h3.h1//Simplify

(*We just have to study the signature of the metric in the Heisenberg group orbits*)

(*To this aim, I project orthogonally the Heisenberg group into the orthogonal complement to the stabilizator*)

h1ort=Ve[[5]];
h2ort=Ve[[6]];
h3ort=Ve[[8]];

KForm[h1ort,h1ort]//Simplify
KForm[h1ort,h2ort]//Simplify
KForm[h1ort,h3ort]//Simplify
KForm[h2ort,h2ort]//Simplify
KForm[h3ort,h2ort]//Simplify
KForm[h3ort,h3ort]//Simplify

(*It has the appropiate signature*)

Spacelike quaternionic paraKähler Heisenberg four-manifolds

(*I clear all previously defined variables and I make the necessary definitions to call the package "diffgeo", which we will extensively use*)
Clear["Global`*"];
coord={t,x,y,z};
X=coord;
metric=DiagonalMatrix[{ϵ,ϵ,1,1}];
<<diffgeo.m

ϵ=-1;

n=4; (*dimensions*)

eucmetric=DiagonalMatrix[{1,1,1,1}];
VB[i_]:=Table[eucmetric[[i,k]],{k,1,n}];

(*This L will contain the Lie brackets of the time-evolving orthonormal frame*)
L=Normal@SymmetrizedArray[{i_,j_,k_}->l[i,j,k],{n,n,n},{Antisymmetric[{1,2}]}];
Do[Do[Do[l[i,j,k]=0,{i,1,n}],{j,1,n}],{k,1,n}];

(*I define the matrix Ut, which relates the frame at t=0 with the frame at t*)

Ut=Transpose[{{1,0,0,0},{0,c[t],f[t],v[t]},{0,h[t],b[t],w[t]},{0,g[t],p[t],a[t]}}];

(*As justified in the paper, we can set:*)

v[t_]:=0;w[t_]:=0;

(*I write the time derivatives of the time-evolving frame*)

det=D[Ut,t].Inverse[Ut]//FullSimplify;

(*Then:*)

Do[Do[l[1,i,k]=det[[i,k]],{k,1,n}],{i,2,n}];

(*In order to write the Lie brackets between the elements of the frame, I define the Lie brackets at t=0*)

L0=Normal@SymmetrizedArray[{i_,j_,k_}->l0[i,j,k],{n,n,n},{Antisymmetric[{1,2}]}];
Do[Do[Do[l0[i,j,k]=0,{i,1,n}],{j,1,n}],{k,1,n}];
l0[2,3,4]=-2 q0;

(*At any time:*)

Do[Do[Do[l[i,j,k]=Sum[Sum[Sum[Ut[[i,r]]Ut[[j,ll]] L0[[r,ll,mm]]Inverse[Ut][[mm,k]]//Simplify,{r,2,n}],{ll,2,n}],{mm,2,n}],{i,2,n}],{j,2,n}],{k,2,n}]

(*I define the Lie bracket between two vectors*)

Lieb[v_,w_]:=Table[Sum[Sum[v[[i]]w[[j]] L[[i,j,k]],{j,1,n}],{i,1,n}]+v[[1]]D[w[[k]],t]-w[[1]] D[v[[k]],t],{k,1,n}];

(*Now we move into the computation of the covariant derivatives of the frame*)

imet=Inverse[metric];
Covdevs=Normal@SymmetrizedArray[{i_,j_,k_}->cd[i,j,k],{n,n,n}];
Do[Do[Do[cd[i,j,q]=1/2 Sum[Sum[L[[i,j,s]] metric[[s,k]] -L[[i,k,s]] metric[[j,s]]-L[[j,k,s]] metric[[i,s]],{s,1,n}]imet[[k,q]],{k,1,n}],{q,1,n}],{j,1,n}],{i,1,n}]

(*The covariant derivative operation is:*)

CD[v_,w_]:=Table[Sum[Sum[v[[i]]w[[j]] Covdevs[[i,j,k]],{j,1,n}],{i,1,n}]+v[[1]] D[w[[k]],t],{k,1,n}];

(*Now we can compute the Riemann tensor as follows*)

Do[Do[Do[Do[RR[i,j,k,q]=CD[VB[i],CD[VB[j],VB[k]]][[q]]-CD[VB[j],CD[VB[i],VB[k]]][[q]]-CD[Lieb[VB[i],VB[j]],VB[k]][[q]],{q,1,n}],{k,1,n}],{j,1,n}],{i,1,n}]//Simplify;

(*The Riemann tensor with components g(R(e_i,e_j)e_k,e_q)σ_q, where σ_q is the signature of the spacetime character of e_q, is*)

Riemcomp=Normal@SymmetrizedArray[{i_,j_,k_,q_}->RR[i,j,k,q],{n,n,n,n}];

(*With all indices down:*)

Do[Do[Do[Do[RRdown[i,j,k,q]=Sum[RR[i,j,k,ll] metric[[ll,q]],{ll,1,n}],{q,1,n}],{k,1,n}],{j,1,n}],{i,1,n}]//Simplify;
Riemcompdown=Normal@SymmetrizedArray[{i_,j_,k_,q_}->RRdown[i,j,k,q],{n,n,n,n}];

(*The Ricci tensor is therefore*)

Ricci=Table[Sum[Riemcomp[[k,i,j,k]],{k,1,n}],{i,1,n},{j,1,n}]//Simplify;

(*We observe*)

Solve[{Ricci[[1,2]]==0,Ricci[[1,3]]==0},{p'[t],g'[t]}]//Simplify

(*Then, as justified in the paper, we set...*)

p[t_]:=0;g[t_]:=0;f[t_]:=-h[t];

Ricci=Ricci//Simplify;

(*Now we proceed to the classification*)

Scal=contract[Ricci]//Simplify;

Z=Ricci-1/4 Scal metric//Simplify;
S=Scal/6 antisymmetrize[swapIndices[metric**metric,{2,4}],{1,2}]//Simplify;

EE=2*antisymmetrize[swapIndices[Z**metric,{2,4}],{1,2},{3,4}]//Simplify;

(*The Weyl tensor is*)

W=((Riemcompdown-EE-S))//Simplify;
(*Consistency check*)

contract[W,{2,4}]//Simplify

(*We compute the dual Weyl tensor*)

eps=LeviCivita;
starW=1/2 contract[eps**W,{3,5},{4,6}]//Simplify;
(*The self-duality tensor is*)

selfduality=(starW- W)//Simplify;

contract[starW,{1,4}]//Simplify

hppsol=Solve[Ricci[[2,3]]==0,h''[t]]//Simplify

selfduality=selfduality/.hppsol[[1]]//Simplify;

selfduality[[2,1,4,2]]//Simplify

(*If the first factor is zero...*)

Ricci[[4,4]]/.{a'[t]->q0 b[t]c[t]+q0 h[t]^2,a''[t]->q0 (b'[t]c[t]+b[t]c'[t]+2 h[t]h'[t])}//Simplify

(*Then the second factor must be zero:*)

h[t_]:=0;

Ricci=(Ricci/.hppsol[[1]])//Simplify;

riccisol=Solve[Ricci==Λ metric,{a'[t],a''[t],b''[t],c''[t]}]//Simplify

selfduality=selfduality/.riccisol[[1]]//Simplify;
Wsim=W/.riccisol[[1]]//Simplify;

selfduality[[3,4,3,4]]//Simplify

lambdasol=Solve[selfduality[[3,4,3,4]]==0,Λ]//Simplify

selfduality=selfduality/.lambdasol[[1]]//Simplify;
Wsim=Wsim/.lambdasol[[1]]//Simplify;
selfduality[[2,3,1,4]]//Simplify

(*If the first terms vanish:*)

(Λ/.lambdasol[[1]])/.{a[t]->-q0 (b[t])^2 c[t]/b'[t]}//Simplify

(Λ/.lambdasol[[1]])/.{a[t]->-q0 (c[t])^2 b[t]/c'[t]}//Simplify

(*Then, the third term must vanish and...*)

c[t_]:=b[t];

(*Y el tensor de Weyl es autodual:*)

selfduality/.lambdasol[[1]]//Simplify

lambdaexpfinal=(Λ/.lambdasol[[1]])//Simplify

(*The expression of a' is...*)

a'[t]/.riccisol[[1]]//Simplify

(*The Weyl tensor is:*)
(Wsim-1/2 contract[eps**Wsim,{3,5},{4,6}])//Simplify
(*I want to diagonalize the Weyl tensor. I define the following basis of self-dual 2-forms*)

bb1=(FormToTensor[Wedge[dd[t],dd[x]]+HodgeStar[Wedge[dd[t],dd[x]]]])//Simplify

bb2=(FormToTensor[Wedge[dd[t],dd[y]]+HodgeStar[Wedge[dd[t],dd[y]]]])//Simplify

bb3=(FormToTensor[Wedge[dd[t],dd[z]]+HodgeStar[Wedge[dd[t],dd[z]]]])//Simplify

bb={bb1,bb2,bb3};

met2form=((Table[contract[bb[[i]]**bb[[j]],{1,3},{2,4}],{i,1,3},{j,1,3}])/.{1/ϵ^2->1,1/ϵ^3->1/ϵ}//Simplify)

(*The (local) matrix elements of the Weyl tensor on the bundle of self-dual 2-forms is*)

W2=(contract[Wsim**Wsim,{1,5},{2,6},{3,7},{4,8}])/.{ϵ^2->1,ϵ^4->1}//Simplify

Wmat=Table[(met2form[[j,j]])^(-1)(contract[contract[Wsim**bb[[i]],{3,5},{4,6}]**bb[[j]],{1,3},{2,4}])/.{ϵ^2->1}//Simplify,{i,1,3},{j,1,3}]

(*This is already the diagonal basis!*)
(*By direct inspection, we observe that Wmat already has two equal eigenvalues*)

(*This implies our space is (locally) conformally Kähler. In fact:*)

ωk1=bb1*conf1[t];
ωk2=bb2*conf1[t];
ωk3=bb3*conf2[t];

W2

extder[ω_]:=Table[D[ω[[j,k]],X[[i]]]-D[ω[[i,k]],X[[j]]]+D[ω[[i,j]],X[[k]]]+Sum[-ω[[s,k]]Lieb[VB[i],VB[j]][[s]]+ω[[s,j]] Lieb[VB[i],VB[k]][[s]]-ω[[s,i]] Lieb[VB[j],VB[k]][[s]],{s,1,n}],{i,1,n},{j,1,n},{k,1,n}]

dωk1=extder[ωk1]//Simplify
dωk2=extder[ωk2]//Simplify
dωk3=extder[ωk3]//Simplify

conf2[t_]:=(((1536 (q0^6) ((b[t])^14))/(((a[t])^4) (((3 q0 ((b[t])^3)+a[t] b'[t]))^2))))^(1/3);

((dωk3/.riccisol[[1]])/.lambdasol[[1]])//Simplify

(contract[ωk3**ωk3,{1,3},{2,4}]*(W2)^(-2/3))/.{pot->1/3}//Simplify

(*Then, ωk1 is Kähler with respect to the rescaled metric \bar{g}=W2^(1/3)g. Interestingly enough: *)

conf1[t_]:=a[t] b[t];

dωk1//Simplify

dωk2//Simplify

contract[ωk1**ωk1,{1,3},{2,4}]*1/(a[t]b[t])^2/.{ϵ^2->1,1/ϵ^3->1/ϵ}//Simplify
contract[ωk2**ωk2,{1,3},{2,4}]*1/(a[t]b[t])^2/.{ϵ^2->1,1/ϵ^3->1/ϵ}//Simplify

J1=(raise[bb1,{1}]//Simplify)/.{1/ϵ^2->1}
J2=(raise[bb2,{1}]//Simplify)/.{1/ϵ^2->1}
J3=(raise[bb3,{1}]//Simplify)/.{1/ϵ^2->1}

J1.J1
J2.J2
J3.J3

NJ1=(((Table[-(J1.J1).Lieb[VB[i],VB[j]]+J1.(Lieb[J1.VB[i],VB[j]]+Lieb[VB[i],J1.VB[j]])-Lieb[J1.VB[i],J1.VB[j]],{i,1,n},{j,1,n}]//Simplify)/.riccisol[[1]])/.lambdasol[[1]])//Simplify

NJ2=(((Table[-(J2.J2).(Lieb[VB[i],VB[j]])+J2.(Lieb[J2.VB[i],VB[j]]+Lieb[VB[i],J2.VB[j]])-Lieb[J2.VB[i],J2.VB[j]],{i,1,n},{j,1,n}]//Simplify)/.riccisol[[1]])/.lambdasol[[1]])//Simplify

NJ3=(Table[-(J3.J3).(Lieb[VB[i],VB[j]])+J3.(Lieb[J3.VB[i],VB[j]]+Lieb[VB[i],J3.VB[j]])-Lieb[J3.VB[i],J3.VB[j]],{i,1,n},{j,1,n}]//Simplify)

Stationary solutions

Ricstat=Ricci/.{a'[t]->a[t] μ[t], a''[t]->a[t] (μ[t])^2+a[t] μ'[t], b'[t]->b[t]λ[t], b''[t]->b[t](λ[t])^2 +b[t]λ'[t]}//Simplify

(*A true stationary solution requires...*)

Ricstat=(Ricstat/.{a[t]->(b[t])^2,μ[t]->2 λ[t]})/.{μ[t]->μ,λ[t]->λ,μ'[t]->0,λ'[t]->0}//Simplify

(*Then:*)

Solve[Ricstat==Λ metric,{λ,Λ}]//Simplify

(*Consequently, picking the minus sign for convenience:*)
a[t_]:=Exp[-2q0 (t-t0)];b[t_]:=Exp[-q0 (t-t0)];

Ricci//Simplify

(*I will need later the Riemann tensor*)

(Riemcomp)//Simplify

(*The remaining task is to prove that it is (locally if ϵ=-1, globally if ϵ=1) isometric to a left-invariant metric on the simply transitive solvable Iwasawa subgroup of SU(1,2) and SU(2,1), respectively if ϵ=1 or ϵ=-1*)

(*We first check that it is a locally symmetric space by computing the covariant derivative of the Riemann tensor and observing that it vanishes*)

Clear[a,b,nablaR]

nablaRiem=Normal@SymmetrizedArray[{i_,j_,k_,q_,s_}->nablaR[i,j,k,q,s],{n,n,n,n,n}];

Do[Do[Do[Do[Do[nablaR[i,j,k,q,s]=eucmetric[[i,1]] D[RRdown[j,k,q,s],t]-Sum[CD[VB[i],VB[j]][[ll]] RRdown[ll,k,q,s],{ll,1,n}]-
Sum[CD[VB[i],VB[k]][[ll]] RRdown[j,ll,q,s],{ll,1,n}]-
Sum[CD[VB[i],VB[q]][[ll]] RRdown[j,k,ll,s],{ll,1,n}]-
Sum[CD[VB[i],VB[s]][[ll]] RRdown[j,k,q,ll],{ll,1,n}],{i,1,n}],{j,1,n}],{k,1,n}],{q,1,n}],{s,1,n}]

a[t_]:=Exp[-2q0 (t-t0)];b[t_]:=Exp[-q0 (t-t0)];

(nablaRiem)//Simplify

(*Now we want to check that these solutions are locally isometric to SL(3,R)/(S(GL(1)xGL(2))*)

Clear["Global`*"];
e1={{-1,0,0},{0,1,0},{0,0,0}};
e2={{-1,0,0},{0,0,0},{0,0,1}};
e3={{0,0,0},{0,0,0},{0,1,0}};
e4={{0,0,0},{0,0,1},{0,0,0}};
e5={{0,1,0},{0,0,0},{0,0,0}};
e6={{0,0,1},{0,0,0},{0,0,0}};
e7={{0,0,0},{1,0,0},{0,0,0}};
e8={{0,0,0},{0,0,0},{1,0,0}};

Ve1=1/Sqrt[6]*(e1+e2);
Ve2=1/Sqrt[2]*(e1-e2);
Ve3=1/Sqrt[2]*(e3+e4);
Ve4=1/Sqrt[2]*(e3-e4);
Ve5=1/Sqrt[2]*(e5+e7);
Ve6=1/Sqrt[2]*(e5-e7);
Ve7=1/Sqrt[2]*(e6+e8);
Ve8=1/Sqrt[2]*(e6-e8);

Ve=μ{Ve1,Ve2,Ve3,Ve4,Ve5,Ve6,Ve7,Ve8};
n=8;
KForm[A_,B_]:=-1/μ^2 Tr[A.B];
metric=Table[KForm[Ve[[i]],Ve[[j]]],{i,1,n},{j,1,n}]//Simplify

L=Normal@SymmetrizedArray[{i_,j_,k_}->l[i,j,k],{n,n,n},{Antisymmetric[{1,2}]}];
Do[Do[Do[l[i,j,k]=0,{i,1,n}],{j,1,n}],{k,1,n}]

Do[Do[Do[l[i,j,k]=Inverse[metric][[k,k]] KForm[(Ve[[i]].Ve[[j]]-Ve[[j]].Ve[[i]]),Ve[[k]]],{k,1,n}],{i,1,n}],{j,1,n}]

(*I notice:*)

Table[L[[i,j,k]],{i,1,4},{j,5,n},{k,1,4}]//Simplify

(*So SL(3,R)/S(GL(1)xGL(2)) is, at the Lie algebra level at least, a reductive space. This means that we can compute its Riemann in the usally easier form:*)

lres=Table[L[[i,j,k]],{i,5,n},{j,5,n},{k,1,4}]//Simplify;

Riemred=-Table[Sum[lres[[i,j,k]] Inverse[metric][[s,s]] KForm[(Ve[[k]].Ve[[ll]]-Ve[[ll]].Ve[[k]]),Ve[[s]]],{k,1,4}],{i,1,4},{j,1,4},{ll,5,8},{s,5,8}]//Simplify;

Riemcomp={{{{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}},{{0,(q0^2),0,0},{-(q0^2),0,0,0},{0,0,0,-(q0^2)},{0,0,(q0^2),0}},{{0,0,(q0^2),0},{0,0,0,(q0^2)},{(q0^2),0,0,0},{0,(q0^2),0,0}},{{0,0,0,4 (q0^2)},{0,0,2 (q0^2),0},{0,2 (q0^2),0,0},{4 (q0^2),0,0,0}}},{{{0,-(q0^2),0,0},{(q0^2),0,0,0},{0,0,0,(q0^2)},{0,0,-(q0^2),0}},{{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}},{{0,0,0,2 (q0^2)},{0,0,4 (q0^2),0},{0,4 (q0^2),0,0},{2 (q0^2),0,0,0}},{{0,0,(q0^2),0},{0,0,0,(q0^2)},{(q0^2),0,0,0},{0,(q0^2),0,0}}},{{{0,0,-(q0^2),0},{0,0,0,-(q0^2)},{-(q0^2),0,0,0},{0,-(q0^2),0,0}},{{0,0,0,-2 (q0^2)},{0,0,-4 (q0^2),0},{0,-4 (q0^2),0,0},{-2 (q0^2),0,0,0}},{{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}},{{0,(q0^2),0,0},{-(q0^2),0,0,0},{0,0,0,-(q0^2)},{0,0,(q0^2),0}}},{{{0,0,0,-4 (q0^2)},{0,0,-2 (q0^2),0},{0,-2 (q0^2),0,0},{-4 (q0^2),0,0,0}},{{0,0,-(q0^2),0},{0,0,0,-(q0^2)},{-(q0^2),0,0,0},{0,-(q0^2),0,0}},{{0,-(q0^2),0,0},{(q0^2),0,0,0},{0,0,0,(q0^2)},{0,0,-(q0^2),0}},{{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}}}};

Ricred=Table[Sum[Riemred[[i,j,k,i]],{i,1,4}],{j,1,4},{k,1,4}]//Simplify

(*Therefore...*)

μ=Sqrt[2] q0;

(*And now I trivially notice*)

ord={1,3,4,2};

Table[metric[[i,k]],{i,5,8},{k,5,8}][[ord,ord]]

Riemcomp- Riemred[[ord,ord,ord,ord]]//Simplify

(*Now we have to find the Heisenberg algebra within sl(3,R))*)
(*By direct inspection*)

(*We observe*)

h1=e4;
h2=2q0 e5;
h3=e6;

(h1.h2-h2.h1+a h3)/.{a->2q0}//Simplify

h2.h3-h3.h2

h1.h3-h3.h1//Simplify

(*At the origin, the orbit of the Heisenberg group is not 3-dim. However, if we conjugate the isotropy group at the origin to obtain the isotropy group at a different point:*)

(*Group element with which I conjugate:*)

g={{(1/v),0,0},{0,(1/v),(1/v)},{(1/v),0,(v^2)}}/.{v->1};

Det[g]

matK={{a,0,0},{0,c,d},{0,e,-c-a}};

Tr[matK]

(*A general element of the Heisenberg algebra*)

heis=b1 e4+b2 e5+ b3 e6;

heisconj=(g.heis.Inverse[g])//Simplify

Solve[heisconj==matK,{a,b1,b2,b3,a,c,d,e}]

(*We observe the conjugated Heisenberg algebra has null intersection with the isotropy algebra*)

(*The conjugated Heisenberg algebra is spanned by:*)

h1conj=Sum[KForm[Ve[[k]],(heisconj/.{b2->0,b3->0,b1->1})] Ve[[k]] metric[[k,k]],{k,5,8}]

h2conj=Sum[KForm[Ve[[k]],(heisconj/.{b1->0,b3->0,b2->1})] Ve[[k]] metric[[k,k]],{k,5,8}]

h3conj=Sum[KForm[Ve[[k]],(heisconj/.{b1->0,b2->0,b3->1})] Ve[[k]] metric[[k,k]],{k,5,8}]

h={h1conj,h2conj,h3conj};

(*The restricted metric is:*)

metrestringida=μ^2Table[KForm[h[[i]],h[[j]]],{i,1,3},{j,1,3}]//Simplify

Det[metrestringida]

(*It has signature (2,1) or (1,2)*)

(*On observing...*)

a2={0,1,0};
a1={1,1,0};

a={a1,a2};

metpos=Table[a[[i]].metrestringida.a[[j]],{i,1,2},{j,1,2}]

(*It has signature (2,1) consequently and we conclude*)

Lightlike quaternionic paraKähler Heisenberg four-manifolds

(*I clear all previously defined variables and I make the necessary definitions to call the package "diffgeo", which we will extensively use*)
Clear["Global`*"];
coord={t,x,y,z};
X=coord;
metric={{ϵ,0,0,0},{0,0,1,0},{0,1,0,0},{0,0,0,1}};
<<diffgeo.m

ϵ=-1;

n=4; (*dimensions*)

eucmetric=DiagonalMatrix[{1,1,1,1}];
VB[i_]:=Table[eucmetric[[i,k]],{k,1,n}];

(*This L will contain the Lie brackets of the time-evolving orthonormal frame*)
L=Normal@SymmetrizedArray[{i_,j_,k_}->l[i,j,k],{n,n,n},{Antisymmetric[{1,2}]}];
Do[Do[Do[l[i,j,k]=0,{i,1,n}],{j,1,n}],{k,1,n}];

(*I define the matrix Ut, which relates the frame at t=0 with the frame at t*)

Ut=Transpose[{{1,0,0,0},{0,c[t],f[t],v[t]},{0,h[t],b[t],w[t]},{0,g[t],p[t],a[t]}}];

(*As justified in the paper, we can set:*)

v[t_]:=0;w[t_]:=0;

(*I write the time derivatives of the time-evolving frame*)

det=D[Ut,t].Inverse[Ut]//FullSimplify;

(*Then:*)

Do[Do[l[1,i,k]=det[[i,k]],{k,1,n}],{i,2,n}];

(*In order to write the Lie brackets between the elements of the frame, I define the Lie brackets at t=0*)

L0=Normal@SymmetrizedArray[{i_,j_,k_}->l0[i,j,k],{n,n,n},{Antisymmetric[{1,2}]}];
Do[Do[Do[l0[i,j,k]=0,{i,1,n}],{j,1,n}],{k,1,n}];
l0[3,4,2]=-2 q0;

(*At any time:*)

Do[Do[Do[l[i,j,k]=Sum[Sum[Sum[Ut[[i,r]]Ut[[j,ll]] L0[[r,ll,mm]]Inverse[Ut][[mm,k]]//Simplify,{r,2,n}],{ll,2,n}],{mm,2,n}],{i,2,n}],{j,2,n}],{k,2,n}]

(*I define the Lie bracket between two vectors*)

Lieb[v_,w_]:=Table[Sum[Sum[v[[i]]w[[j]] L[[i,j,k]],{j,1,n}],{i,1,n}]+v[[1]]D[w[[k]],t]-w[[1]] D[v[[k]],t],{k,1,n}];

(*Now we move into the computation of the covariant derivatives of the frame*)

imet=Inverse[metric];
Covdevs=Normal@SymmetrizedArray[{i_,j_,k_}->cd[i,j,k],{n,n,n}];
Do[Do[Do[cd[i,j,q]=1/2 Sum[Sum[L[[i,j,s]] metric[[s,k]] -L[[i,k,s]] metric[[j,s]]-L[[j,k,s]] metric[[i,s]],{s,1,n}]imet[[k,q]],{k,1,n}],{q,1,n}],{j,1,n}],{i,1,n}]

(*The covariant derivative operation is:*)

CD[v_,w_]:=Table[Sum[Sum[v[[i]]w[[j]] Covdevs[[i,j,k]],{j,1,n}],{i,1,n}]+v[[1]] D[w[[k]],t],{k,1,n}];

(*Now we can compute the Riemann tensor as follows*)

Do[Do[Do[Do[RR[i,j,k,q]=CD[VB[i],CD[VB[j],VB[k]]][[q]]-CD[VB[j],CD[VB[i],VB[k]]][[q]]-CD[Lieb[VB[i],VB[j]],VB[k]][[q]],{q,1,n}],{k,1,n}],{j,1,n}],{i,1,n}]//Simplify;

(*The Riemann tensor with components g(R(e_i,e_j)e_k,e_q)σ_q, where σ_q is the signature of the spacetime character of e_q, is*)

Riemcomp=Normal@SymmetrizedArray[{i_,j_,k_,q_}->RR[i,j,k,q],{n,n,n,n}];

(*With all indices down:*)

Do[Do[Do[Do[RRdown[i,j,k,q]=Sum[RR[i,j,k,ll] metric[[ll,q]],{ll,1,n}],{q,1,n}],{k,1,n}],{j,1,n}],{i,1,n}]//Simplify;
Riemcompdown=Normal@SymmetrizedArray[{i_,j_,k_,q_}->RRdown[i,j,k,q],{n,n,n,n}];

(*The Ricci tensor is therefore*)

Ricci=Table[Sum[Riemcomp[[k,i,j,k]],{k,1,n}],{i,1,n},{j,1,n}]//Simplify;

(*We observe...*)

ghsol=Solve[{Ricci[[1,3]]==0,Ricci[[1,4]]==0},{g'[t],h'[t]}]//Simplify

ghsol[[1]]/.{h[t]->0,g[t]->0}//Simplify

(*Hence...*)

g[t_]:=0;h[t_]:=0;c[t_]:=1;

Ricci=Ricci//Simplify;

Ricci[[3,4]]

(*We follow the methods of the paper*)

appsol=Solve[Ricci[[4,4]]==Λ metric[[4,4]],a''[t]]//Simplify

Ricci=Ricci/.appsol[[1]]//Simplify;

fpsol=Solve[{Ricci[[3,4]]==0,Ricci[[3,3]]==0},{f''[t],p''[t]}]

Ricci=Ricci/.fpsol[[1]]//Simplify;

bppsol=Solve[Ricci==Λ metric,{b''[t],a'[t]}]//Simplify

Ricci=Ricci/.bppsol[[1]]//Simplify

(*We see that the metric is Einstein*)

(*Now we compute the Weyl self-duality tensor*)

Ricci=Ricci//Simplify;
Scal=contract[Ricci]//Simplify;

Z=Ricci-1/4 Scal metric//Simplify;
S=Scal/6 antisymmetrize[swapIndices[metric**metric,{2,4}],{1,2}]//Simplify;

EE=2*antisymmetrize[swapIndices[Z**metric,{2,4}],{1,2},{3,4}]//Simplify;

(*The Weyl tensor is*)

W=(((Riemcompdown-EE-S))/.appsol[[1]]/.fpsol[[1]]/.bppsol[[1]])//Simplify;
(*Consistency check*)

contract[W,{1,4}]//Simplify

(*We compute the dual Weyl tensor*)

eps=LeviCivita/.{ϵ^2->1,Sign[ϵ]^2->1};
starW=1/2 contract[eps**W,{3,5},{4,6}]//Simplify;
(*The self-duality tensor is*)

selfduality=(starW- W)/.{ϵ^2->1}//Simplify;

selfduality[[2,1,1,3]]//Simplify

lambdasol=Solve[selfduality[[2,1,1,3]]==0,Λ]//Simplify

selfduality=selfduality/.lambdasol[[1]]//Simplify;

selfduality[[1,4,3,4]]//Simplify

(*If the first factor vanishes...*)

(Λ/.lambdasol[[1]])/.{b'[t]->0}//Simplify

(*So the second factor must vanish*)

p[t_]:=0;

(*Similarly*)

selfduality[[3,1,3,1]]

(*So*)

f[t_]:=0;

bppsol/.lambdasol[[1]]//Simplify

(*We observe it is actually conformally-flat*)

W/.lambdasol[[1]]//Simplify

DSolve[{a'[t]==((a[t] b'[t])/(2 b[t])),((3 ((b'[t])^2))/(4 ((b[t])^2)))==Λ},{a[t],b[t]},t]

(*We pick one solution, because the other is obtained by time inversion*)

a[t_]:=Exp[(t-t0)*Sqrt[Λ/3]];
b[t_]:=Exp[2(t-t0)*Sqrt[Λ/3]];

Geodesics

eut0=FormToTensor[dd[z]+q0 x dd[y]-q0  y dd[x]];
evt0={0,1,0,0};
e3t0={0,0,1,0};
metlight=-{1,0,0,0}**{1,0,0,0}+Exp[-2(t-t0)*Sqrt[Λ/3]]eut0**evt0+Exp[-2(t-t0)*Sqrt[Λ/3]]evt0**eut0+Exp[-2(t-t0)*Sqrt[Λ/3]]e3t0**e3t0//Simplify

Clear["Global`*"];
coord={t,x,y,z};
metric={{-1,0,0,0},{0,-2 (E^(-((2 (t-t0) (Sqrt[Λ]))/(Sqrt[3])))) q0 y,(E^(-((2 (t-t0) (Sqrt[Λ]))/(Sqrt[3])))) q0 x,(E^(-((2 (t-t0) (Sqrt[Λ]))/(Sqrt[3]))))},{0,(E^(-((2 (t-t0) (Sqrt[Λ]))/(Sqrt[3])))) q0 x,(E^(-((2 (t-t0) (Sqrt[Λ]))/(Sqrt[3])))),0},{0,(E^(-((2 (t-t0) (Sqrt[Λ]))/(Sqrt[3])))),0,0}};
<<diffgeo.m

RicciTensor- Λ metric//Simplify

(*I want to derive the geodesic equation, assuming a curve is affinely parametrized*)

curv={tt[τ],xx[τ],yy[τ],zz[τ]};
Vec=D[curv,τ]//Simplify;
geodeq=D[curv,{τ,2}]+contract[((Christoffel)/.{x->xx[τ],y->yy[τ],z->zz[τ],t->tt[τ]})**lower[Vec]**lower[Vec],{2,4},{3,5}]//Simplify

(*By direct insprection, looking for particular solutions*)

xx[τ_]:=0;
zz[τ_]:=0;

geodeq//Simplify

DSolve[-((2 (Sqrt[Λ]) tt'[τ] yy'[τ])/(Sqrt[3]))+yy''[τ]==0,yy[τ],τ]//Simplify

yy[τ_]:=Integrate[Exp[((2 (Sqrt[Λ]) tt[σ])/(Sqrt[3]))],{σ,0,τ}]

geodeq//Simplify

t0=0;

DSolve[-(((E^((2 (Sqrt[Λ]) (t0+tt[τ]))/(Sqrt[3]))) (Sqrt[Λ]))/(Sqrt[3]))+tt''[τ]==0,tt[τ],τ]//Simplify

tt[τ_]:=-Sqrt[3]/(2Sqrt[Λ]) Log[4 Λ/3 (Sinh[τ/2+B])^2];
yy[τ_]:=3/(4Λ)Integrate[1/(Sinh[σ/2+B])^2,{σ,0,τ}];
yyd[τ_]:=3/(4Λ)1/(Sinh[τ/2+B])^2; (*I define separately the derivative of yy because Mathematica has some strange problems dealing with derivatives of Integrals*)

-((2 (Sqrt[Λ]) tt'[τ] yyd[τ])/(Sqrt[3]))+yyd'[τ]//Simplify

-(((E^((2 (Sqrt[Λ]) (t0-tt[τ]))/(Sqrt[3]))) (Sqrt[Λ]) ((yyd[τ])^2))/(Sqrt[3]))+tt''[τ]

(*Choosing B=0, this geodesic cannot be extended further away from τ=-B, and we conclude*)

(Timelike) (para)hyperKähler Heisenberg four-manifolds

(*I clear all previously defined variables and I make the necessary definitions to call the package "diffgeo", which we will extensively use*)
Clear["Global`*"];
coord={t,x,y,z};
X=coord;
metric=DiagonalMatrix[{ϵ,ϵ,1,1}];
<<diffgeo.m

n=4; (*dimensions*)

eucmetric=DiagonalMatrix[{1,1,1,1}];
VB[i_]:=Table[eucmetric[[i,k]],{k,1,n}];

(*This L will contain the Lie brackets of the time-evolving orthonormal frame*)
L=Normal@SymmetrizedArray[{i_,j_,k_}->l[i,j,k],{n,n,n},{Antisymmetric[{1,2}]}];
Do[Do[Do[l[i,j,k]=0,{i,1,n}],{j,1,n}],{k,1,n}];

(*I define the matrix Ut, which relates the frame at t=0 with the frame at t*)

Ut=Transpose[{{1,0,0,0},{0,a[t],p[t],g[t]},{0,v[t],b[t],h[t]},{0,w[t],f[t],c[t]}}];

v[t_]:=0;w[t_]:=0;p[t_]:=0;g[t_]:=0;f[t_]:=0;

(*I write the time derivatives of the time-evolving frame*)

det=D[Ut,t].Inverse[Ut]//FullSimplify;

(*Then:*)

Do[Do[l[1,i,k]=det[[i,k]],{k,1,n}],{i,2,n}];

(*In order to write the Lie brackets between the elements of the frame, I define the Lie brackets at t=0*)

L0=Normal@SymmetrizedArray[{i_,j_,k_}->l0[i,j,k],{n,n,n},{Antisymmetric[{1,2}]}];
Do[Do[Do[l0[i,j,k]=0,{i,1,n}],{j,1,n}],{k,1,n}];
l0[3,4,2]=-2 q0;

(*At any time:*)

Do[Do[Do[l[i,j,k]=Sum[Sum[Sum[Ut[[i,r]]Ut[[j,ll]] L0[[r,ll,mm]]Inverse[Ut][[mm,k]]//Simplify,{r,2,n}],{ll,2,n}],{mm,2,n}],{i,2,n}],{j,2,n}],{k,2,n}]

(*I define the Lie bracket between two vectors*)

Lieb[v_,w_]:=Table[Sum[Sum[v[[i]]w[[j]] L[[i,j,k]],{j,1,n}],{i,1,n}]+v[[1]]D[w[[k]],t]-w[[1]] D[v[[k]],t],{k,1,n}];

(*Now we move into the computation of the covariant derivatives of the frame*)

imet=Inverse[metric];
Covdevs=Normal@SymmetrizedArray[{i_,j_,k_}->cd[i,j,k],{n,n,n}];
Do[Do[Do[cd[i,j,q]=1/2 Sum[Sum[L[[i,j,s]] metric[[s,k]] -L[[i,k,s]] metric[[j,s]]-L[[j,k,s]] metric[[i,s]],{s,1,n}]imet[[k,q]],{k,1,n}],{q,1,n}],{j,1,n}],{i,1,n}]

(*The covariant derivative operation is:*)

CD[v_,w_]:=Table[Sum[Sum[v[[i]]w[[j]] Covdevs[[i,j,k]],{j,1,n}],{i,1,n}]+v[[1]] D[w[[k]],t],{k,1,n}];

(*Now we can compute the Riemann tensor as follows*)

Do[Do[Do[Do[RR[i,j,k,q]=CD[VB[i],CD[VB[j],VB[k]]][[q]]-CD[VB[j],CD[VB[i],VB[k]]][[q]]-CD[Lieb[VB[i],VB[j]],VB[k]][[q]],{q,1,n}],{k,1,n}],{j,1,n}],{i,1,n}]//Simplify;

(*The Riemann tensor with components g(R(e_i,e_j)e_k,e_q)σ_q, where σ_q is the signature of the spacetime character of e_q, is*)

Riemcomp=Normal@SymmetrizedArray[{i_,j_,k_,q_}->RR[i,j,k,q],{n,n,n,n}];

(*With all indices down:*)

Do[Do[Do[Do[RRdown[i,j,k,q]=Sum[RR[i,j,k,ll] metric[[ll,q]],{ll,1,n}],{q,1,n}],{k,1,n}],{j,1,n}],{i,1,n}]//Simplify;
Riemcompdown=Normal@SymmetrizedArray[{i_,j_,k_,q_}->RRdown[i,j,k,q],{n,n,n,n}];

(*The Ricci tensor is therefore*)

Ricci=Table[Sum[Riemcomp[[k,i,j,k]],{k,1,n}],{i,1,n},{j,1,n}]//Simplify;

Ricci=Ricci/.{ϵ^2->1}//Simplify;

(*I define a self-dual two-form*)

selfomega={{0,ω12[t,x,y,z],ω13[t,x,y,z],ω14[t,x,y,z]},{-ω12[t,x,y,z],0,ϵ ω14[t,x,y,z],-ϵ ω13[t,x,y,z]},{-ω13[t,x,y,z],-ϵ ω14[t,x,y,z],0,ω12[t,x,y,z]},{-ω14[t,x,y,z],ϵ ω13[t,x,y,z],-ω12[t,x,y,z],0}};
(selfomega-FormToTensor[HodgeStar[TensorToForm[selfomega]]]/.{ϵ^2->1,Sign[ϵ]->1}//Simplify)/.{ϵ^2->1,1/ϵ->ϵ,1/ϵ^2->1}
(*If it represents the Heisenberg-invariant Kähler form:*)
ω12[t_,x_,y_,z_]:=ω12[t];ω13[t_,x_,y_,z_]:=ω13[t];ω14[t_,x_,y_,z_]:=ω14[t];

(*I compute the covariant detivative:*)

nablaw=Table[D[contract[selfomega**VB[j]**VB[k],{1,3},{2,4}],X[[i]]]-contract[selfomega**CD[VB[i],VB[j]]**VB[k],{1,3},{2,4}]-contract[selfomega**VB[j]**CD[VB[i],VB[k]],{1,3},{2,4}],{i,1,4},{j,1,4},{k,1,4}];

(*Parallelicity implies:*)

(*ϵ=1*)
ϵ=1;

sol1=Solve[nablaw[[1]]==0,{ω12'[t],ω13'[t],ω14'[t]}]//Simplify;
sol2=Solve[nablaw[[2]]==0,a'[t]]//Simplify;
sol3=Solve[nablaw[[3]]==0,{b'[t],h'[t]}]//Simplify;
sol4=Solve[nablaw[[4]]==0,{c'[t],h'[t]}]//Simplify;

sol1

sol2

sol3

sol4

(*Then...*)

h[t_]:=0;

(*The remaining ODEs are solved by...*)

DSolve[{-q0 b[t] c[t]+a'[t]==0,b'[t]==-((q0 ((b[t])^2) c[t])/(a[t])),c'[t]==-((q0 b[t] ((c[t])^2))/(a[t]))},{a[t],b[t],c[t]},t]

a[t_]:=(1+3q0(t-t0))^(1/3);
b[t_]:=(1+3q0(t-t0))^(-1/3);
c[t_]:=(1+3q0(t-t0))^(-1/3);

nablaw//Simplify

(*ϵ=1*)
Clear[ϵ,a,b,c,h]
ϵ=-1;

sol1=Solve[nablaw[[1]]==0,{ω12'[t],ω13'[t],ω14'[t]}]//Simplify;
sol2=Solve[nablaw[[2]]==0,a'[t]]//Simplify;
sol3=Solve[nablaw[[3]]==0,{b'[t],h'[t]}]//Simplify;
sol4=Solve[nablaw[[4]]==0,{c'[t],h'[t]}]//Simplify;

sol1

sol2

sol3

sol4

(*Then...*)

h[t_]:=0;

(*The remaining ODEs are solved by...*)

DSolve[{-q0 b[t] c[t]+a'[t]==0,b'[t]==-((q0 ((b[t])^2) c[t])/(a[t])),c'[t]==-((q0 b[t] ((c[t])^2))/(a[t]))},{a[t],b[t],c[t]},t]

a[t_]:=(1+3q0(t-t0))^(1/3);
b[t_]:=(1+3q0(t-t0))^(-1/3);
c[t_]:=(1+3q0(t-t0))^(-1/3);

nablaw//Simplify

Ricci

(*I define the exterior derivative operation*)

Clear[ϵ]

der[u_]:=-Table[Sum[u[[k]]Lieb[VB[i],VB[j]][[k]],{k,1,n}],{i,1,n},{j,1,n}];

extder[ω_]:=Table[D[ω[[j,k]],X[[i]]]-D[ω[[i,k]],X[[j]]]+D[ω[[i,j]],X[[k]]]+Sum[-ω[[s,k]]Lieb[VB[i],VB[j]][[s]]+ω[[s,j]] Lieb[VB[i],VB[k]][[s]]-ω[[s,i]] Lieb[VB[j],VB[k]][[s]],{s,1,n}],{i,1,n},{j,1,n},{k,1,n}]

(*We prove that we get a (para)hyperKähler structure*)

ω1=(Wedge[dd[t],dd[x]]+HodgeStar[Wedge[dd[t],dd[x]]])/.{ϵ^2->1,Sign[ϵ]->1,1/ϵ^2->1}//Simplify
ω2=(Wedge[dd[t],dd[y]]+HodgeStar[Wedge[dd[t],dd[y]]])/.{ϵ^2->1,Sign[ϵ]->1,1/ϵ^2->1}//Simplify
ω3=(Wedge[dd[t],dd[z]]+HodgeStar[Wedge[dd[t],dd[z]]])/.{ϵ^2->1,Sign[ϵ]->1,1/ϵ^2->1}//Simplify

ω1t=FormToTensor[ω1]//Simplify;
ω2t=FormToTensor[ω2]//Simplify;
ω3t=FormToTensor[ω3]//Simplify;

(*I compute directly the exterior derivatives*)

Clear[a,b,c]

dω1t=extder[ω1t]//Simplify;

dω2t=extder[ω2t]//Simplify;

dω3t=extder[ω3t]//Simplify;

a[t_]:=(1+3q0(t-t0))^(1/3);
b[t_]:=(1+3q0(t-t0))^(-1/3);
c[t_]:=(1+3q0(t-t0))^(-1/3);

dω1t//Simplify

dω2t//Simplify

dω3t//Simplify

(*I check antiself-duality*)

Ricci=Ricci/.{ϵ^2->1}//Simplify;
Scal=contract[Ricci]//Simplify;

Z=Ricci-1/4 Scal metric//Simplify;
S=Scal/6 antisymmetrize[swapIndices[metric**metric,{2,4}],{1,2}]//Simplify;

EE=2*antisymmetrize[swapIndices[Z**metric,{2,4}],{1,2},{3,4}]//Simplify;

(*The Weyl tensor is*)

W=((Riemcompdown-EE-S))/.{ϵ^2->1}//Simplify;
(*Consistency check*)

contract[W,{1,4}]/.{ϵ^2->1}//Simplify

(*We compute the dual Weyl tensor, and we observe it is antiself-dual*)

eps=LeviCivita/.{ϵ^2->1,Sign[ϵ]^2->1};
starW=1/2 contract[eps**W,{3,5},{4,6}]//Simplify;
(*The self-duality tensor is*)

antiselfduality=(starW+ W)/.{ϵ^2->1}//Simplify;

antiselfduality/.{ϵ^2->1}//Simplify

Spacelike parahyperKähler Heisenberg four-manifolds

(*I clear all previously defined variables and I make the necessary definitions to call the package "diffgeo", which we will extensively use*)
Clear["Global`*"];
coord={t,x,y,z};
X=coord;
metric=DiagonalMatrix[{ϵ,ϵ,1,1}];
<<diffgeo.m

ϵ=-1;

n=4; (*dimensions*)

eucmetric=DiagonalMatrix[{1,1,1,1}];
VB[i_]:=Table[eucmetric[[i,k]],{k,1,n}];

(*This L will contain the Lie brackets of the time-evolving orthonormal frame*)
L=Normal@SymmetrizedArray[{i_,j_,k_}->l[i,j,k],{n,n,n},{Antisymmetric[{1,2}]}];
Do[Do[Do[l[i,j,k]=0,{i,1,n}],{j,1,n}],{k,1,n}];

(*I define the matrix Ut, which relates the frame at t=0 with the frame at t*)

Ut=Transpose[{{1,0,0,0},{0,c[t],f[t],w[t]},{0,h[t],b[t],v[t]},{0,g[t],p[t],a[t]}}];

f[t_]:=-h[t];w[t_]:=0;v[t_]:=0;g[t_]:=0;p[t_]:=0;

(*I write the time derivatives of the time-evolving frame*)

det=D[Ut,t].Inverse[Ut]//FullSimplify;

(*Then:*)

Do[Do[l[1,i,k]=det[[i,k]],{k,1,n}],{i,2,n}];

(*In order to write the Lie brackets between the elements of the frame, I define the Lie brackets at t=0*)

L0=Normal@SymmetrizedArray[{i_,j_,k_}->l0[i,j,k],{n,n,n},{Antisymmetric[{1,2}]}];
Do[Do[Do[l0[i,j,k]=0,{i,1,n}],{j,1,n}],{k,1,n}];
l0[2,3,4]=-2 q0;

(*At any time:*)

Do[Do[Do[l[i,j,k]=Sum[Sum[Sum[Ut[[i,r]]Ut[[j,ll]] L0[[r,ll,mm]]Inverse[Ut][[mm,k]]//Simplify,{r,2,n}],{ll,2,n}],{mm,2,n}],{i,2,n}],{j,2,n}],{k,2,n}]

(*I define the Lie bracket between two vectors*)

Lieb[v_,w_]:=Table[Sum[Sum[v[[i]]w[[j]] L[[i,j,k]],{j,1,n}],{i,1,n}]+v[[1]]D[w[[k]],t]-w[[1]] D[v[[k]],t],{k,1,n}];

(*Now we move into the computation of the covariant derivatives of the frame*)

imet=Inverse[metric];
Covdevs=Normal@SymmetrizedArray[{i_,j_,k_}->cd[i,j,k],{n,n,n}];
Do[Do[Do[cd[i,j,q]=1/2 Sum[Sum[L[[i,j,s]] metric[[s,k]] -L[[i,k,s]] metric[[j,s]]-L[[j,k,s]] metric[[i,s]],{s,1,n}]imet[[k,q]],{k,1,n}],{q,1,n}],{j,1,n}],{i,1,n}]

(*The covariant derivative operation is:*)

CD[v_,w_]:=Table[Sum[Sum[v[[i]]w[[j]] Covdevs[[i,j,k]],{j,1,n}],{i,1,n}]+v[[1]] D[w[[k]],t],{k,1,n}];

(*Now we can compute the Riemann tensor as follows*)

Do[Do[Do[Do[RR[i,j,k,q]=CD[VB[i],CD[VB[j],VB[k]]][[q]]-CD[VB[j],CD[VB[i],VB[k]]][[q]]-CD[Lieb[VB[i],VB[j]],VB[k]][[q]],{q,1,n}],{k,1,n}],{j,1,n}],{i,1,n}]//Simplify;

(*The Riemann tensor with components g(R(e_i,e_j)e_k,e_q)σ_q, where σ_q is the signature of the spacetime character of e_q, is*)

Riemcomp=Normal@SymmetrizedArray[{i_,j_,k_,q_}->RR[i,j,k,q],{n,n,n,n}];

(*With all indices down:*)

Do[Do[Do[Do[RRdown[i,j,k,q]=Sum[RR[i,j,k,ll] metric[[ll,q]],{ll,1,n}],{q,1,n}],{k,1,n}],{j,1,n}],{i,1,n}]//Simplify;
Riemcompdown=Normal@SymmetrizedArray[{i_,j_,k_,q_}->RRdown[i,j,k,q],{n,n,n,n}];

(*The Ricci tensor is therefore*)

Ricci=Table[Sum[Riemcomp[[k,i,j,k]],{k,1,n}],{i,1,n},{j,1,n}]//Simplify;

(*I define a self-dual two-form*)

selfomega={{0,ω12[t,x,y,z],ω13[t,x,y,z],ω14[t,x,y,z]},{-ω12[t,x,y,z],0,ϵ ω14[t,x,y,z],-ϵ ω13[t,x,y,z]},{-ω13[t,x,y,z],-ϵ ω14[t,x,y,z],0,ω12[t,x,y,z]},{-ω14[t,x,y,z],ϵ ω13[t,x,y,z],-ω12[t,x,y,z],0}};
(selfomega-FormToTensor[HodgeStar[TensorToForm[selfomega]]])//Simplify
(*If it represents the Heisenberg-invariant Kähler form:*)
ω12[t_,x_,y_,z_]:=ω12[t];ω13[t_,x_,y_,z_]:=ω13[t];ω14[t_,x_,y_,z_]:=ω14[t];

(*I compute the covariant detivative:*)

nablaw=Table[D[contract[selfomega**VB[j]**VB[k],{1,3},{2,4}],X[[i]]]-contract[selfomega**CD[VB[i],VB[j]]**VB[k],{1,3},{2,4}]-contract[selfomega**VB[j]**CD[VB[i],VB[k]],{1,3},{2,4}],{i,1,4},{j,1,4},{k,1,4}];

(*Parallelicity implies:*)

sol1=Solve[nablaw[[1]]==0,{ω12'[t],ω13'[t],ω14'[t]}]//Simplify;
sol2=Solve[nablaw[[2]]==0,{c'[t],h'[t]}]//Simplify;
sol3=Solve[nablaw[[3]]==0,{b'[t],h'[t]}]//Simplify;
sol4=Solve[nablaw[[4]]==0,a'[t]]//Simplify;

sol1

sol2

sol3

sol4

(*I see*)

h[t_]:=0;

(*The remaining EDOs are solved by*)

DSolve[{q0 b[t] c[t]+a'[t]==0,b'[t]==((q0 ((b[t])^2) c[t])/(a[t])),c'[t]==((q0 b[t] ((c[t])^2))/(a[t]))},{a[t],b[t],c[t]},t]

a[t_]:=1/(1-3q0(t-t0))^(-1/3);
b[t_]:=(1-3q0(t-t0))^(-1/3);
c[t_]:=(1-3q0(t-t0))^(-1/3);

nablaw//Simplify

Ricci//Simplify

(*We prove that we get a (para)hyperKähler structure*)

der[u_]:=-Table[Sum[u[[k]]Lieb[VB[i],VB[j]][[k]],{k,1,n}],{i,1,n},{j,1,n}];

extder[ω_]:=Table[D[ω[[j,k]],X[[i]]]-D[ω[[i,k]],X[[j]]]+D[ω[[i,j]],X[[k]]]+Sum[-ω[[s,k]]Lieb[VB[i],VB[j]][[s]]+ω[[s,j]] Lieb[VB[i],VB[k]][[s]]-ω[[s,i]] Lieb[VB[j],VB[k]][[s]],{s,1,n}],{i,1,n},{j,1,n},{k,1,n}]

ω1=(Wedge[dd[t],dd[x]]+HodgeStar[Wedge[dd[t],dd[x]]])//Simplify
ω2=(Wedge[dd[t],dd[y]]+HodgeStar[Wedge[dd[t],dd[y]]])//Simplify
ω3=(Wedge[dd[t],dd[z]]+HodgeStar[Wedge[dd[t],dd[z]]])//Simplify

ω1t=FormToTensor[ω1]//Simplify;
ω2t=FormToTensor[ω2]//Simplify;
ω3t=FormToTensor[ω3]//Simplify;

(*I compute directly the exterior derivatives*)

Clear[a,b,c]

dω1t=extder[ω1t]//Simplify;

dω2t=extder[ω2t]//Simplify;

dω3t=extder[ω3t]//Simplify;

a[t_]:=(1-3q0(t-t0))^(1/3);
b[t_]:=(1-3q0(t-t0))^(-1/3);
c[t_]:=(1-3q0(t-t0))^(-1/3);

dω1t//Simplify

dω2t//Simplify

dω3t//Simplify

(*I check antiself-duality*)

Ricci=Ricci//Simplify;
Scal=contract[Ricci]//Simplify;

Z=Ricci-1/4 Scal metric//Simplify;
S=Scal/6 antisymmetrize[swapIndices[metric**metric,{2,4}],{1,2}]//Simplify;

EE=2*antisymmetrize[swapIndices[Z**metric,{2,4}],{1,2},{3,4}]//Simplify;

(*The Weyl tensor is*)

W=((Riemcompdown-EE-S))//Simplify;
(*Consistency check*)

contract[W,{1,4}]//Simplify

(*We compute the dual Weyl tensor, and we observe it is antiself-dual*)

eps=LeviCivita;
starW=1/2 contract[eps**W,{3,5},{4,6}]//Simplify;
(*The self-duality tensor is*)

antiselfduality=(starW+ W)//Simplify;

antiselfduality//Simplify

Lightlike parahyperKähler Heisenberg four-manifolds

(*I clear all previously defined variables and I make the necessary definitions to call the package "diffgeo", which we will extensively use*)
Clear["Global`*"];
coord={t,x,y,z};
X=coord;
metric={{ϵ,0,0,0},{0,0,1,0},{0,1,0,0},{0,0,0,1}};
<<diffgeo.m

ϵ=-1;

n=4; (*dimensions*)

eucmetric=DiagonalMatrix[{1,1,1,1}];
VB[i_]:=Table[eucmetric[[i,k]],{k,1,n}];

(*This L will contain the Lie brackets of the time-evolving orthonormal frame*)
L=Normal@SymmetrizedArray[{i_,j_,k_}->l[i,j,k],{n,n,n},{Antisymmetric[{1,2}]}];
Do[Do[Do[l[i,j,k]=0,{i,1,n}],{j,1,n}],{k,1,n}];

(*I define the matrix Ut, which relates the frame at t=0 with the frame at t*)

Ut=Transpose[{{1,0,0,0},{0,c[t],f[t],w[t]},{0,h[t],b[t],v[t]},{0,g[t],p[t],a[t]}}];

c[t_]:=1;w[t_]:=0;g[t_]:=0;v[t_]:=0;h[t_]:=0;

(*I write the time derivatives of the time-evolving frame*)

det=D[Ut,t].Inverse[Ut]//FullSimplify;

(*Then:*)

Do[Do[l[1,i,k]=det[[i,k]],{k,1,n}],{i,2,n}];

(*In order to write the Lie brackets between the elements of the frame, I define the Lie brackets at t=0*)

L0=Normal@SymmetrizedArray[{i_,j_,k_}->l0[i,j,k],{n,n,n},{Antisymmetric[{1,2}]}];
Do[Do[Do[l0[i,j,k]=0,{i,1,n}],{j,1,n}],{k,1,n}];
l0[3,4,2]=-2 q0;

(*At any time:*)

Do[Do[Do[l[i,j,k]=Sum[Sum[Sum[Ut[[i,r]]Ut[[j,ll]] L0[[r,ll,mm]]Inverse[Ut][[mm,k]]//Simplify,{r,2,n}],{ll,2,n}],{mm,2,n}],{i,2,n}],{j,2,n}],{k,2,n}]

(*I define the Lie bracket between two vectors*)

Lieb[v_,w_]:=Table[Sum[Sum[v[[i]]w[[j]] L[[i,j,k]],{j,1,n}],{i,1,n}]+v[[1]]D[w[[k]],t]-w[[1]] D[v[[k]],t],{k,1,n}];

(*Now we move into the computation of the covariant derivatives of the frame*)

imet=Inverse[metric];
Covdevs=Normal@SymmetrizedArray[{i_,j_,k_}->cd[i,j,k],{n,n,n}];
Do[Do[Do[cd[i,j,q]=1/2 Sum[Sum[L[[i,j,s]] metric[[s,k]] -L[[i,k,s]] metric[[j,s]]-L[[j,k,s]] metric[[i,s]],{s,1,n}]imet[[k,q]],{k,1,n}],{q,1,n}],{j,1,n}],{i,1,n}]

(*The covariant derivative operation is:*)

CD[v_,w_]:=Table[Sum[Sum[v[[i]]w[[j]] Covdevs[[i,j,k]],{j,1,n}],{i,1,n}]+v[[1]] D[w[[k]],t],{k,1,n}];

(*Now we can compute the Riemann tensor as follows*)

Do[Do[Do[Do[RR[i,j,k,q]=CD[VB[i],CD[VB[j],VB[k]]][[q]]-CD[VB[j],CD[VB[i],VB[k]]][[q]]-CD[Lieb[VB[i],VB[j]],VB[k]][[q]],{q,1,n}],{k,1,n}],{j,1,n}],{i,1,n}]//Simplify;

(*The Riemann tensor with components g(R(e_i,e_j)e_k,e_q)σ_q, where σ_q is the signature of the spacetime character of e_q, is*)

Riemcomp=Normal@SymmetrizedArray[{i_,j_,k_,q_}->RR[i,j,k,q],{n,n,n,n}];

(*With all indices down:*)

Do[Do[Do[Do[RRdown[i,j,k,q]=Sum[RR[i,j,k,ll] metric[[ll,q]],{ll,1,n}],{q,1,n}],{k,1,n}],{j,1,n}],{i,1,n}]//Simplify;
Riemcompdown=Normal@SymmetrizedArray[{i_,j_,k_,q_}->RRdown[i,j,k,q],{n,n,n,n}];

(*The Ricci tensor is therefore*)

Ricci=Table[Sum[Riemcomp[[k,i,j,k]],{k,1,n}],{i,1,n},{j,1,n}]//Simplify;

(*I define a self-dual two-form*)

selfomega={{0,ω12[t,x,y,z],ω13[t,x,y,z],ω14[t,x,y,z]},{-ω12[t,x,y,z],0,- ω14[t,x,y,z],ω12[t,x,y,z]},{-ω13[t,x,y,z], ω14[t,x,y,z],0,-ω13[t,x,y,z]},{-ω14[t,x,y,z],-ω12[t,x,y,z],ω13[t,x,y,z],0}};;
(selfomega-FormToTensor[HodgeStar[TensorToForm[selfomega]]])//Simplify
(*If it represents the Heisenberg-invariant Kähler form:*)
ω12[t_,x_,y_,z_]:=ω12[t];ω13[t_,x_,y_,z_]:=ω13[t];ω14[t_,x_,y_,z_]:=ω14[t];

(*I compute the covariant detivative:*)

nablaw=Table[D[contract[selfomega**VB[j]**VB[k],{1,3},{2,4}],X[[i]]]-contract[selfomega**CD[VB[i],VB[j]]**VB[k],{1,3},{2,4}]-contract[selfomega**VB[j]**CD[VB[i],VB[k]],{1,3},{2,4}],{i,1,4},{j,1,4},{k,1,4}];

(*Parallelicity implies:*)

sol1=Solve[nablaw[[1]]==0,{(ω12)'[t],(ω13)'[t],(ω14)'[t]}]//Simplify;
sol2=Solve[nablaw[[2]]==0,b'[t]]//Simplify;
sol3=Solve[nablaw[[3]]==0,{f'[t],p'[t]}]//Simplify;
sol4=Solve[nablaw[[4]]==0,{a'[t],p'[t]}]//Simplify;

sol1

sol2

sol3

sol4

(*We observe:*)
a[t_]:=1;p[t_]:=0;b[t_]:=1;f[t_]:=-2 q0 (t-t0);

nablaw//Simplify

Ricci//Simplify

(*We prove that we get a (para)hyperKähler structure*)

der[u_]:=-Table[Sum[u[[k]]Lieb[VB[i],VB[j]][[k]],{k,1,n}],{i,1,n},{j,1,n}];

extder[ω_]:=Table[D[ω[[j,k]],X[[i]]]-D[ω[[i,k]],X[[j]]]+D[ω[[i,j]],X[[k]]]+Sum[-ω[[s,k]]Lieb[VB[i],VB[j]][[s]]+ω[[s,j]] Lieb[VB[i],VB[k]][[s]]-ω[[s,i]] Lieb[VB[j],VB[k]][[s]],{s,1,n}],{i,1,n},{j,1,n},{k,1,n}]

(*The Kähler forms*)

ω1=(Wedge[dd[t],dd[x]]+HodgeStar[Wedge[dd[t],dd[x]]])//Simplify
ω2=(Wedge[dd[t],dd[y]]+HodgeStar[Wedge[dd[t],dd[y]]])//Simplify
ω3=(Wedge[dd[t],dd[z]]+HodgeStar[Wedge[dd[t],dd[z]]])//Simplify

ω1t=FormToTensor[ω1]//Simplify;
ω2t=FormToTensor[ω2]//Simplify;
ω3t=FormToTensor[ω3]//Simplify;

(*I compute directly the exterior derivatives*)

dω1t=extder[ω1t]//Simplify;

dω2t=extder[ω2t]//Simplify;

dω3t=extder[ω3t]//Simplify;

dω1t//Simplify

dω2t//Simplify

dω3t//Simplify

(*We observe the space is actually flat*)

Riemcomp//Simplify

Geodesics

eut0=FormToTensor[dd[z]+q0 x dd[y]-q0  y dd[x]];
evt0={0,1,0,0};
e3t0={0,0,1,0};
metlight=-{1,0,0,0}**{1,0,0,0}+eut0**evt0+evt0**eut0+4 q0 (t-t0) evt0**evt0+e3t0**e3t0//Simplify

Clear["Global`*"];
coord={t,x,y,z};
metric={{-1,0,0,0},{0,2 q0 (2 t-2 t0-y),q0 x,1},{0,q0 x,1,0},{0,1,0,0}};
<<diffgeo.m

RicciTensor

(*I want to derive the geodesic equation, assuming a curve is affinely parametrized*)

Vec

curv={tt[τ],xx[τ],yy[τ],zz[τ]};
Vec=D[curv,τ]//Simplify;
geodeq=D[curv,{τ,2}]+contract[((Christoffel)/.{x->xx[τ],y->yy[τ],z->zz[τ],t->tt[τ]})**lower[Vec]**lower[Vec],{2,4},{3,5}]//Simplify

(*By direct insprection*)

xx[τ_]:=a4+a3 τ;

geodeq//Simplify

(*By direct insprection*)

yy[τ_]:=a5+a6 τ-a3^2 q0 τ^2;

geodeq//Simplify

tt[τ_]:=a1+a2 τ-a3^2 q0 τ^2;

geodeq//Simplify

zz[τ_]:=a7+a8 τ+ q0 a3 τ^2(q0 a3 xx[τ]-2 a2+a6);

geodeq//Simplify