Curvature computations

Curvature computations are probably the most important type of computation in GR, and we need to do it in a simple and fast enough way. We present next a number of curvature computations for different metrics corresponding to known exact solutions in GR. The fundamental command in this case is MetricCompute.

<<xAct`xCoba`

?MetricCompute

To have info messages from ComponentValue, use this:

$CVVerbose=True;

To see messages from MetricCompute use this:

SetOptions[MetricCompute,Verbose->True]

Exit

Combination of TensorValues and CTensor frameworks

Load xCoba:

Exit

<<xAct`xCoba`

<<xAct`ShowTime1`

$PrePrint=ScreenDollarIndices;

$DefInfoQ=False;

$CVVerbose=False;

We define a manifold:

DefManifold[M,4,IndexRange[a,p]]

We define a metric and introduce a coordinate chart. We give values to the metric compo:

DefMetric[-1,metricg[-a,-b],CD,PrintAs->"g"]

DefConstantSymbol[Mass,PrintAs->"M"]

We define a metric and introduce a coordinate chart. We give values to the metric components at definition time:

DefChart[SCHW,M,{0,1,2,3},{t[],r[],θ[],ϕ[]},MetricInBasis->{metricg,-SCHW,DiagonalMatrix[{-1+2Mass/r[],1/(1-2Mass/r[]),r[]^2,r[]^2 Sin[θ[]]^2}]},ChartColor->Blue]

TensorValues@metricg

Change the default formatting of the tensor slots when using base indices:

?CIndexForm

SCHW/:CIndexForm[0,SCHW]:="t";
SCHW/:CIndexForm[1,SCHW]:="r";
SCHW/:CIndexForm[2,SCHW]:="θ";
SCHW/:CIndexForm[3,SCHW]:="φ";

{metricg[{0,-SCHW},{0,-SCHW}],metricg[{0,SCHW},{1,SCHW}]}

Use MetricCompute to compute the components of a number of curvature tensors:

?MetricCompute

MetricCompute[metricg,SCHW,All,Parallelize->True,Verbose->True]

Computations are cached:

PrintxCobaCache[]

We check the index configuration of some of the values computed:

TensorValIDs@RiemannCD

Last@TensorValues@RiemannCD

RicciCD[a,b]/.MakeRule[{RicciCD[-a,-b],RicciCD[-a,-b]},MetricOn->All,ContractMetrics->False]

RicciCD[a,b]RicciCD[-a,-b]

SeparateMetric[]@%

%65;

ToCTensor[RicciCD,{-SCHW,-SCHW}]

%//ToValues

%/.RicciCD->Ricci[cd]

cd

?SeparateMetric

RicciCD[a,b]==metricg[a,d]metricg[b,f]RicciCD[-d,-f]

ToBasis[SCHW]/@%

TraceBasisDummy/@%

ToValues@ComponentArray@%

TensorValues@RicciCD

TensorValues@metricg

TensorValues@ChristoffelCDPDSCHW

TensorValues@KretschmannCD

Let us repeat the same computation but now within the CTensor framework. First we transform the metric into a CTensor:

ToCTensor[metricg,{-SCHW,-SCHW}]

ToValues@%

newmetric=%

We can use SetCMetric in the usual way, but now we must bear in mind that the CTensor metric becomes the “second metric”:

SetCMetric[newmetric,SCHW,SignatureOfMetric->{3,1,0}]

$Metrics

MetricsOfVBundle@TangentM

This means that index raising and lowering is carried out with the first metric, whose components are stored according to the TensorValues framework:

newmetric[a,-b]

%//ToValues//Simplify

We can check that in any case SetCMetric carries out the computations as expected:

PrintxCobaCache[]

We have a look at some of the quantities computed by SetCMetric:

cd=CovDOfMetric@newmetric

Christoffel[cd,PDSCHW][a,-b,-c]

Again, remember that indices are raised and lowered with the first metric:

Christoffel[cd,PDSCHW][a,b,-c]

%//ToValues

Usage of MetricCompute with the CTensor metric:

MetricCompute[newmetric,SCHW,All,Parallelize->True,Verbose->True]

PrintxCobaCache[]

We have now components stored in two different ways:

Ricci[cd]

TensorValues@RicciCD

Check that the computations carried out with TensorValues framework and CTensor framework agree:

RiemannDown[cd][-a,-b,-c,-d]//ToBasisExpand

ToValues@TraceBasisDummy@SeparateBasis[SCHW]@RiemannCD[-a,-b,-c,-d]

%%-%

It is possible to use the CTensor metric to raise and lower indices, given that both are in this case the same, but it has to be done under the user responsibility. Change the order of the metrics defined in the tangent vector bundle for that:

TangentM/:MetricsOfVBundle@TangentM=Reverse@MetricsOfVBundle@TangentM;

MetricsOfVBundle@TangentM

newmetric[a,-b]

Simplify@RiemannDown[cd][-a,-b,-c,d]===Simplify@Riemann[cd][-a,-b,-c,d]

Computation of Cartan scalars for the Schwarzschild metric in a null tetrad

Exit

<<xAct`xCoba`

<<xAct`ShowTime1`

$DefInfoQ=False;

$PrePrint=ScreenDollarIndices;

$CVVerbose=False;

DefManifold[M,4,IndexRange[a,p]]

We introduce the Schwarzschild metric:

DefConstantSymbol[Mass,PrintAs->"M"]

DefChart[SCHW,M,{0,1,2,3},{t[],r[],θ[],ϕ[]},ChartColor->Blue]

metricg=CTensor[DiagonalMatrix[{1-2Mass/r[],-1/(1-2Mass/r[]),-r[]^2,-r[]^2 Sin[θ[]]^2}],{-SCHW,-SCHW}];

Delta=r[]^2-2Mass r[];

Introduce the null tetrad:

DefBasis[NP,TangentM,{1,2,3,4},BasisColor->RGBColor[1,0,0]]

We specify the behaviour of the elements of this basis under complex conjugation. This is achieved with the command SetDaggerMatrix.

SetDaggerMatrix[NP,{{1,0,0,0},{0,1,0,0},{0,0,0,1},{0,0,1,0}}]

Relation between the null tetrad and the Schwarzschild holonomic frame:

{Basis[{1,-NP},a]==1/Delta r[]^2 Basis[{0,-SCHW},a]+Basis[{1,-SCHW},a],Basis[{2,-NP},a]==1/2Basis[{0,-SCHW},a]-Delta/(2r[]^2)Basis[{1,-SCHW},a],Basis[{3,-NP},a]==1/(r[]Sqrt[2])(Basis[{2,-SCHW},a]+I Csc[θ[]] Basis[{3,-SCHW},a]),Dagger/@(Basis[{3,-NP},a]==1/(r[]Sqrt[2])(Basis[{2,-SCHW},a]+I Csc[θ[]] Basis[{3,-SCHW},a]))}

Map[Dagger,%,{2}]

Map[ToBasis[SCHW],%,{2}]

Flatten@ComponentArray@%

%/.Equal->Rule

Components of the change basis matrix:

changebasis=ComponentArray@Basis[-{a,NP},{b,SCHW}]/.%

Set the relation between the two bases, in order that we can transform CTensor objects from one basis to another :

SetBasisChange[CTensor[changebasis,{-NP,SCHW}],SCHW]

PrintxCobaCache[]

We check that the frame transformation has the expected properties:

metricNP=ToCTensor[metricg,{-NP,-NP}]

ToCTensor[metricNP,{-SCHW,-SCHW}]

Simplify@%===Simplify@metricg

Choose the metric in the null tetrad as the metric for doing computations:

?SetCMetric

SetCMetric[metricNP,SCHW,SignatureOfMetric->{1,3,0}]

MetricCompute[metricNP,SCHW,All,Parallelize->True,Verbose->True]

Ricci rotation coefficients in the null tetrad:

cd=CovDOfMetric[metricNP]

cd2=ToCCovD[cd,PDSCHW,{SCHW,-SCHW,-SCHW}];

Change of basis to compute the Christoffel symbols in the Schwarzschild coordinates:

ToCCovD[cd,PDSCHW,{SCHW,-SCHW,-SCHW},{-SCHW,-SCHW}]

Check of the vacuum condition:

Ricci[cd]

Riemann[cd]

Computation of higher order covariant derivatives of the Cartan scalars:

$xCobaCacheVerbose=True;

RiemannDown[cd][-a,-b,-c,-d]

TensorDerivative[RiemannDown[cd],cd]

When computer higher order covariant derivatives of the Riemann tensor the components of the change basis matrix do not get automatically replaced. We enforce the automatic replacement:

RuleToSet@Basis;

TensorDerivative[RiemannDown[cd],cd,cd]//Simplify

%[-a,-b,-c,-d,-e,-f]

TensorDerivative[RiemannDown[cd],cd,cd,cd];

TensorDerivative[RiemannDown[cd],cd,cd,cd,cd];

TensorDerivative[RiemannDown[cd],cd,cd,cd,cd,cd];

TensorDerivative[RiemannDown[cd],cd,cd,cd,cd,cd]//ByteCount

Have a look at one of the Cartan scalars:

TensorDerivative[RiemannDown[cd],cd,cd,cd,cd,cd][{1,-NP},{2,-NP},{1,-NP},{2,-NP},{2,-NP},{2,-NP},{2,-NP},{2,-NP},{2,-NP}]

Extract cache information:

DownValues@xCobaCacheTable//Length

ByteCount/@DownValues@xCobaCacheTable

pp-Wave: computation of Cartan scalars

Needs["xAct`xCoba`"]

<<xAct`ShowTime1`

$DefInfoQ=$UndefInfoQ=False;
$CVVerbose=False;
$PrePrint=ScreenDollarIndices;

Define a 4d manifold:

DefManifold[M,4,{a,b,c,d,e,f}]

Coordinates & wave amplitude:

DefChart[PP,M,{0,1,2,3},{u[],v[],x[],y[]},FormatBasis->{"Partials","Differentials"}]

DefScalarFunction[H]

Formatting of coordinate tensor slots:

PP/:CIndexForm[0,PP]="u";
PP/:CIndexForm[1,PP]="v";
PP/:CIndexForm[2,PP]="x";
PP/:CIndexForm[3,PP]="y";

The metric tensor:

Basis[-a,{0,PP}]Basis[-b,{0,PP}]H[x[],y[]]+Basis[-a,{0,PP}]Basis[-b,{1,PP}]+Basis[-a,{1,PP}]Basis[-b,{0,PP}]-Basis[-a,{2,PP}]Basis[-b,{2,PP}]-Basis[-a,{3,PP}]Basis[-b,{3,PP}]

ComponentArray@ToBasis[PP]@%

metric=CTensor[%,{-PP,-PP}]

SetCMetric[metric,PP,SignatureOfMetric->{1,3,0}]

PrintxCobaCache[]

MetricCompute[metric,PP,All,Parallelize->True,Verbose->True]

cd=CovDOfMetric[metric];

Einstein[cd][-a,-b]

Now it is easy to compute the components of the higher order covariant derivatives of the Riemann tensor:

IndicesOfVBundle@TangentM

$xCobaCacheVerbose=True;

TensorDerivative[RiemannDown[cd],cd][-a,-b,-c,-d,-e];

TensorDerivative[RiemannDown[cd],cd,cd][-a,-b,-c,-d,-e,-f];

TensorDerivative[RiemannDown[cd],cd,cd,cd][-a,-b,-c,-d,-e,-f,-f1];

TensorDerivative[RiemannDown[cd],cd,cd,cd,cd];

TensorDerivative[RiemannDown[cd],cd,cd,cd,cd,cd];

TensorDerivative[RiemannDown[cd],cd,cd,cd,cd,cd,cd];

$xCobaCacheVerbose=False;

Check that higher order scalar curvature invariants are zero:

TensorDerivative[RiemannDown[cd],cd][-a,-b,-c,-d,-e]TensorDerivative[RiemannDown[cd],cd][a,b,c,d,e]

TensorDerivative[RiemannDown[cd],cd,cd][-a,-b,-c,-d,-e,-f]TensorDerivative[RiemannDown[cd],cd,cd][a,b,c,d,e,f]

TensorDerivative[RiemannDown[cd],cd,cd,cd][-a,-b,-c,-d,-e,-f,-f1]TensorDerivative[RiemannDown[cd],cd,cd,cd][a,b,c,d,e,f,f1]

Repeat the computations in a null tetrad:

DefBasis[NP,TangentM,{1,2,3,4},BasisColor->RGBColor[0,0,1]]

We specify the behaviour of the elements of this basis under complex conjugation. This is achieved with the command SetDaggerMatrix.

SetDaggerMatrix[NP,{{1,0,0,0},{0,1,0,0},{0,0,0,1},{0,0,1,0}}]

Computation of the null tetrad:

H[x[],y[]]du^2 +2Sqrt@H[x[],y[]]du dv/Sqrt@H[x[],y[]]+dv^2/H[x[],y[]]

(Sqrt@H[x[],y[]]du+dv/Sqrt@H[x[],y[]])^2-dv^2/H[x[],y[]]

{e0->Sqrt@H[x[],y[]]du+dv/Sqrt@H[x[],y[]],e1->dv/Sqrt[H[x[],y[]]]}

{l==1/Sqrt[2](e0+e1),n==1/Sqrt[2](e0-e1),m==1/Sqrt[2](e2+I  e3),mb==1/Sqrt[2](e2-I  e3)}/.%

%/.{l->Basis[-a,{1,NP}],du->Basis[-a,{0,PP}],dv->Basis[-a,{1,PP}],n->Basis[-a,{2,NP}],m->Basis[-a,{3,NP}],mb->Basis[-a,{4,NP}],e2->Basis[-a,{2,PP}],e3->Basis[-a,{3,PP}]}

Map[ToBasis[PP],%,{2}]

Flatten@ComponentArray@%

%/.Equal->ComponentValue

basischange=ToValues@ComponentArray@Basis[-{a,PP},{b,NP}]

H/:WeightOf@H[x[],y[]]=0;

SetBasisChange[CTensor[change,{-PP,NP}],PP]

TensorValues@Basis

metricNP=ToCTensor[metric,{-NP,-NP}]

metricNP[-a,-b]

%//ToBasisExpand

Compute curvature tensors for the metric in the non-holonomic frame:

UnsetCMetric[metric]

$Metrics

SetCMetric[metricNP,PP,SignatureOfMetric->{1,3,0}]

MetricCompute[metricNP,PP,All,Parallelize->True,Verbose->True]

$LargeComponentSize=10000;

Extract information:

cd=CovDOfMetric[metricNP];

Ricci rotation coefficients in the non-holonomic frame:

Christoffel[cd,PDNP][a,-b,-c]

%//ToBasisExpand

Consistency checks:

RicciScalar[cd][]

Einstein[cd]

Kretschmann[cd][]

$xCobaCacheVerbose=True;

$CVSimplify=Together;

RuleToSet@Basis;

TensorDerivative[RiemannDown[cd],cd];

TensorDerivative[RiemannDown[cd],cd,cd];

TensorDerivative[RiemannDown[cd],cd,cd,cd];

TensorDerivative[RiemannDown[cd],cd,cd,cd,cd];

TensorDerivative[RiemannDown[cd],cd,cd,cd,cd,cd];

TensorDerivative[RiemannDown[cd],cd,cd,cd,cd,cd,cd];

ByteCount/@DownValues@xCobaCacheTable

Kerr-Newman in an orthogonal frame

Needs["xAct`xCoba`"]

<<xAct`ShowTime1`

$DefInfoQ=$UndefInfoQ=False;
$CVVerbose=False;
$PrePrint=ScreenDollarIndices;

Define a 4d manifold:

DefManifold[M,4,{a,b,c,d,e,f}]

Boyer-Lindquist coordinates & Kerr-Newman parameters:

DefChart[BL,M,{0,1,2,3},{t[],r[],θ[],ϕ[]},FormatBasis->{"Partials","Differentials"}]

DefConstantSymbol[Mass,PrintAs->"M"]
DefConstantSymbol[Rotation,PrintAs->"a"]
DefConstantSymbol[Charge,PrintAs->"Q"]

Change of basis between the non-holonomic and the holonomic frame:

ρ2=r[]^2+Rotation^2Cos[θ[]]^2;

change={{((r[]^2+Rotation^2)/(ρ2)),0,0,(Rotation/(ρ2))},{0,1,0,0},{0,0,1,0},{Rotation Sin[θ[]]^2,0,0,1}}

Now define a frame in terms of the BL coordinates. Give the relation to the BL coordinates at definition time:

DefBasis[frame,TangentM,{0,1,2,3},BasisColor->RGBColor[0,0,1],BasisChange->CTensor[change,{-frame,BL}]]

This explains the change we have provided:

Basis[-{a,frame},b]==TraceBasisDummy@SeparateBasis[BL]@Basis[-{a,frame},b]

ToValues@ComponentArray@%

The metric tensor is diagonal in the non-holonomic frame:

factor=((r[]^2)/(ρ2))(1-((2Mass)/(r[]))+((Rotation^2+Charge^2)/(r[]^2)));

metric=CTensor[DiagonalMatrix[{-factor,1/factor,ρ2,ρ2 Sin[θ[]]^2}],{-frame,-frame},0]

The coordinated metric tensor is non-diagonal:

metricBL=ToCTensor[metric,{-BL,-BL}]

%[-a,-b]

BasisExpand[%,BL]

Compute curvature tensors for the metric in the non-holonomic frame:

MetricCompute[metric,BL,All,Parallelize->True,Verbose->True]

That computed and simplified everything, including the Kretschmann scalar. All those results are now cached. If we run it again the computation is much faster:

MetricCompute[metric,BL,All,Verbose->True]

Ricci rotation coefficients in the non-holonomic frame:

cd=CovDOfMetric[metric];

Christoffel[cd,PDframe][a,-b,-c]

%//ToBasisExpand

Christoffel symbols in BL coordinates:

cd2=ToCCovD[cd,PDBL,{BL,-BL,-BL},{-BL,-BL}]//Simplify;

Christoffel[cd2,PDBL][a,-b,-c]

%//ToBasisExpand

Consistency check: we compute the Ricci scalar using the Ricci rotation coefficients of the orthogonal frame (first computation) and the Christoffel symbols in BL coordinates (second computation). The second computation takes longer because we did not use MetricCompute on the metric written in the holonomic frame:

RicciScalar[cd][]

$xCobaCacheVerbose=True;

RicciScalar[cd2][]

The Ricci tensor is again diagonal in the orthonormal frame:

Ricci[cd]

%[-a,-b]

Kretschmann[cd][]

Check the Bianchi identity in the orthogonal frame:

cdriemdown=TensorDerivative[RiemannDown[cd],cd];

cdriemdown[-a,-b,-c,-d,-e]+cdriemdown[-b,-c,-a,-d,-e]+cdriemdown[-c,-a,-b,-d,-e]//Simplify//Timing

Clean up:

Remove[metric,cd,factor,ρ2,change,cdriemdown]

UndefConstantSymbol[{Mass,Rotation,Charge}];

UndefTensor[JacobianBLframe]

UndefBasis[frame]

UndefChart[BL]

UndefManifold[M]

Clear the cache:

ClearxCobaCache[All]

Exit

Kerr solution in Boyer Lindqvist coordinates and a null tetrad

Exit

Needs["xAct`xTerior`"]

<<xAct`ShowTime1`

$DefInfoQ=$UndefInfoQ=False;
$CVVerbose=False;
$PrePrint=ScreenDollarIndices;

Define a 4d manifold:

DefManifold[M,4,{a,b,c,d,e,f}]

Boyer-Lindquist coordinates & Kerr-Newman parameters:

DefChart[BL,M,{0,1,2,3},{t[],r[],θ[],ϕ[]},FormatBasis->{"Partials","Differentials"}]

DefConstantSymbol[Mass,PrintAs->"M"]
DefConstantSymbol[Rotation,PrintAs->"a"]

Introducing the null tetrad

Introduce the null tetrad:

DefBasis[NP,TangentM,{1,2,3,4},BasisColor->RGBColor[0,0,1]]

We specify the behaviour of the elements of this basis under complex conjugation. This is achieved with the command SetDaggerMatrix.

SetDaggerMatrix[NP,{{1,0,0,0},{0,1,0,0},{0,0,0,1},{0,0,1,0}}]

eta=(r[]-Mass)/Sqrt[(Mass^2-Rotation^2)];

Delta=(Mass^2-Rotation^2)(eta^2-1)

rho2=r[]^2+Rotation^2 Cos[θ[]]^2

rhobar=r[]+I Rotation Cos[θ[]]

Relation between the null tetrad and BL holonomic frame:

{Basis[{1,-NP},a]==1/Delta((r[]^2+Rotation^2)Basis[{0,-BL},a]+Delta Basis[{1,-BL},a]+Basis[{3,-BL},a]Rotation),Basis[{2,-NP},a]==1/(2 rho2)((r[]^2+Rotation^2)Basis[{0,-BL},a]-Delta Basis[{1,-BL},a]+Basis[{3,-BL},a]Rotation),Basis[{3,-NP},a]==1/(Sqrt[2]rhobar)(I Basis[{0,-BL},a]Rotation Sin[θ[]]+Basis[{2,-BL},a]+I Basis[{3,-BL},a]Csc[θ[]]),Dagger/@(Basis[{3,-NP},a]==1/(Sqrt[2]rhobar)(I Basis[{0,-BL},a]Rotation Sin[θ[]]+Basis[{2,-BL},a]+I Basis[{3,-BL},a]Csc[θ[]]))}

Map[ToBasis[BL],%,{2}]

Flatten@ComponentArray@%

Map[Simplify,%,{2}]

%/.Equal->ComponentValue;

TensorValues@Basis

TensorValIDs@Basis

matrixchange=ToValues@ToCTensor[Basis,{-NP,BL}]

$LargeComponentSize=8000;

SetBasisChange[matrixchange,BL]

PrintxCobaCache[]

metricNP=CTensor[
 {{0,1,0,0},{1,0,0,0},{0,0,0,-1},{0,0,-1,0}}, 
 {-NP, -NP}, 0]

Metric in the null tetrad:

metricNP[-a,-b]

Metric in the BL coordinates:

ToCTensor[metricNP,{-BL,-BL}]//Simplify

%//FullSimplify

metricBL=%;

Computation of connection & curvature using the null tetrad:

SetCMetric[metricNP,BL,SignatureOfMetric->{1,3,0}]

MetricCompute[metricNP,BL,All,Parallelize->True,Verbose->True]

cd=CovDOfMetric[metricNP];

Ricci[cd]

Weyl[cd][-a,-b,-c,-d]//FullSimplify

Weyl scalar (Subscript[Ψ,  2]):

Weyl[cd][{1,-NP},{3,-NP},{2,-NP},{4,-NP}]

Computation of the surface gravity

We work now in the holonomic frame of BL coordinates:

UnsetCMetric@metricNP

SetCMetric[metricBL,BL,SignatureOfMetric->{1,3,0}]

Axial & stationary Killing vectors:

{vt=CTensor[{1,0,0,0},{BL}],vφ=CTensor[{0,0,0,1},{BL}]}

Killing vector matrix

matrixK={{metricBL[-a,-b]vt[a]vt[b],metricBL[-a,-b]vt[a]vφ[b]},{metricBL[-a,-b]vφ[a]vt[b],metricBL[-a,-b]vφ[a]vφ[b]}}

Computation of the event horizon=Killing horizon:

Simplify@Det@%

%[[2]]

Solve[%==0,r[]]//Flatten

EventHorizon=Last@%

Killing matrix at the event horizon:

matrixH=matrixK/.EventHorizon

Simplify//@%

The Killing matrix is degenerate at the Killing horizon. This means that at least one of its eigenvalues is zero:

Eigenvectors@%

The first eigenvector corresponds to the zero eigenvalue:

matrixH.(zerov=%[[1]])

Simplify@%

The eigenvector does not depend on the angle θ. It yields the null generator of the horizon:

kv=CTensor[{zerov[[1]],0,0,zerov[[-1]]},{BL}]

%//TrigExpand//Simplify

%//FullSimplify

Normalization:

Rotation/(2Mass(Mass+Sqrt[Mass^2-Rotation^2]))%

kv=%;

We read off the angular rotation of the horizon:

Ω=kv[{3,BL}]

BL coordinates are not well-behaved at the event horizon. To avoid that problem we use the NP tetrad:

UnsetCMetric@metricBL

SetCMetric[metricNP,BL,SignatureOfMetric->{1,3,0}]

cd=CovDOfMetric@metricNP;

Null generator in the NP tetrad:

kvNP=ToCTensor[kv,{NP}]

kvNP[a]

Automatic replacement of the components of the Jacobian matrix relating the NP and the BL frames:

RuleToSet@Basis

We apply the formula that gives the surface gravity: if (ξ^a) is the Killing vector corresponding to the null generator of the horizon then the surface gravity κ is given by  (ξ^a)(Subscript[∇, a])(ξ^b)=κ (ξ^b):

kvNP[a]cd[-a]@kvNP[b]

%/.EventHorizon

Simplify//@%

FullSimplify//@%

Extract the non-vanishing component:

k1=(Head@%)[{2,NP}]

kvNP[b]/.EventHorizon;

Simplify//@%

Extract the non-vanishing component:

k2=(Head@%)[{2,NP}]

Surface gravity:

k1/k2

FullSimplify@%

Computations with xTerior

Relation between the co-bases:

Coframe[M][a]==SeparateBasis[BL]@dx[M][a]

TraceBasisDummy/@%

ToBasis[NP]/@%

ToValues@ComponentArray@%

%/.Equal->ComponentValue;

Solve[%%,Union@Cases[%%,_Diff,Infinity]]//Flatten//Simplify

diffToNP=%;

Cartan structure equations:

Diff/@ComponentArray@Coframe[M][{a,NP}]

Equal[#,UseCartan[#,cd]]&/@%;

Map[Simplification,%,{2}]

Map[Collect[#,_Wedge,Simplify]&,%,{2}]

CartanStructureEqs=%;

Complete Integrability of the Cartan system

Map[Diff,CartanStructureEqs,{2}]

%/.diffToNP;

Map[ToCanonical,%,{2}]

Map[UseCartan[#,cd]&,%,{2}];

Simplification//@%

Connection 1-forms:

ChristoffelForm[cd][a,-b]

Curvature 2-forms (not simplified by default):

RiemannForm[cd][a,-b]

%/.CTensor[expr_,slots__,weight_]:>CTensor[Simplification@expr,slots,weight]

Hodge Laplacian computed with the metric in the null tetrad

DefScalarFunction[Φ]

UseDimensionStart[]

Codiff[metricNP]@Diff@Φ[t[],r[],θ[],ϕ[]]+Diff@Codiff[metricNP]@Φ[t[],r[],θ[],ϕ[]]

Expand the co-differential first:

CodiffToDiff@%

%/.diffToNP

Work out the Hodge duals. Note that we need to use ExpandHodgeDual several times:

ExpandHodgeDual[%,Coframe[M],metricNP]

ExpandHodgeDual[%,Coframe[M],metricNP]

%/.diffToNP

UseCartan[%,cd]

ExpandHodgeDual[%,Coframe[M],metricNP]

Simplify@%

FullSimplify@%

hp=%

Hodge Laplacian computed with the metric in the BL coordinates (simpler):

SignatureOfMetric@metricBL^={1,3,0};

Determinant@metricBL

SignDetOfMetric@metricBL^=-1;

Codiff[metricBL]@Diff@Φ[t[],r[],θ[],ϕ[]]+Diff@Codiff[metricBL]@Φ[t[],r[],θ[],ϕ[]]

Expand the co-differential first:

CodiffToDiff@%

ExpandHodgeDual[%,dx[M],metricBL]

ExpandHodgeDual[%,dx[M],metricBL]

PowerExpand@%;

Simplify@%

Check that the result agrees with the result of the previous subsection:

%-hp

Simplify@%

General 2d metric connection with torsion

Needs["xAct`xCoba`"]

$DefInfoQ=$UndefInfoQ=False;
$CVVerbose=False;

Define a 2d manifold:

DefManifold[M,2,{a,b,c,d,e,f}]

with a general chart:

DefChart[gen,M,{1,2},{x[],y[]},FormatBasis->{"Partials","Differentials"}]

Let us now work with a fully general 2-dimensional metric:

DefScalarFunction/@{gxx,gxy,gyy};

PrintAs@gxx^="(Subscript[g, xx])";PrintAs@gxy^="(Subscript[g, xy])";PrintAs@gyy^="(Subscript[g, yy])";

(metric=CTensor[{{gxx[x[],y[]],gxy[x[],y[]]},{gxy[x[],y[]],gyy[x[],y[]]}},{-gen,-gen}])[-a,-b]

SymmetryGroupOfTensor[metric]

In this case we will also have torsion, a general torsion tensor:

DefScalarFunction/@{T112,T212};

PrintAs@T112^="(Subscript[(T^x), xy])";PrintAs@T212^="(Subscript[(T^y), xy])";

(torsion=CTensor[{{{0,T112[x[],y[]]},{-T112[x[],y[]],0}},{{0,T212[x[],y[]]},{-T212[x[],y[]],0}}},{gen,-gen,-gen}])[a,-b,-c]

SymmetryGroupOfTensor[torsion]

Construct a covariant derivative that is associated to the metric and has the given torsion tensor:

cd=CovDOfMetric[metric,torsion]

It has the torsion we provided:

Torsion[cd]===torsion

Therefore this covariant derivative is not the Levi-Civita connection of the metric. We can obtain the Levi-Civita connection specifying that torsion is Zero (which is actually the default):

lc=CovDOfMetric[metric,Zero]

Torsion[lc]

There is a shortcut to get the Levi-Civita connection:

LC[metric]===lc

Both covariant derivatives are associated to the metric, in the sense that:

cd[-a][metric[-b,-c]]

lc[-a][metric[-b,-c]]

Same computation done with TensorDerivative:

TensorDerivative[metric,cd]

TensorDerivative[metric,lc]

The Christoffel tensor relating cd to the partial derivative of the chart is not symmetric:

christoffel=Christoffel[cd,PDgen]

SymmetryGroupOfTensor@%

Twice its antisymmetric part is precisely given by the torsion of cd:

christoffel[a,-b,-c]-christoffel[a,-c,-b]==torsion[a,-b,-c]//Simplify

Christoffel tensor relating the covariant derivative with torsion and the Levi-Civita connection:

chris=Christoffel[cd,lc]//Simplify

Its antisymmetric part is the original torsion tensor (torsion of cd):

chris[a,-b,-c]-chris[a,-c,-b]==torsion[a,-b,-c]//Simplify

Its symmetric part is a complicated object:

chris[a,-b,-c]+chris[a,-c,-b]//Simplify

Head[%]

SymmetryGroupOfTensor@%

Let us now look at the curvature tensors:

riemann=Riemann[cd]

To analyze the symmetries of the result we need the Riemann with four indices down:

riemanndown=RiemannDown[cd]

In this 2-dimensional case we get the usual symmetries of the Riemann tensor:

SymmetryGroupOfTensor[riemann]

SymmetryGroupOfTensor[riemanndown]

Now compute everything for the covariant derivative cd. Note that we need to include the torsion as one of the arguments of MetricCompute. Most of the time is spent in the simplification of the results:

MetricCompute[metric,torsion,gen,All,Parallelize->True,Verbose->True]//AbsoluteTiming

Ricci-scalar curvature:

RicciScalar[cd][]

Einstein tensor is zero, even having torsion:

Einstein[cd][-a,-b]

The Levi-Civita connection has different curvature (note that curvature objects not computed for the Levi-Civita connection must be computed first):

$xCobaCacheVerbose=True;

RicciScalar[lc]

RicciScalar[cd]-RicciScalar[lc]

%//Simplify

Full cache table:

PrintxCobaCache[]

Clean up:

Remove[metric,torsion,cd,christoffel,riemann,riemanndown,lc,chris]

UndefScalarFunction[$ScalarFunctions];

UndefChart[gen]

UndefManifold[M]

ClearxCobaCache[All]

PrintxCobaCache[]

Exit

A hard example: The Pomeransky-Senkov metric

The Pomeransky-Senkov metric is a vacuum solution of the Einstein field equations in 5 dimensions. It represents a doubly rotating black ring.

Needs["xAct`xCoba`"]

<<xAct`ShowTime1`

$DefInfoQ=$UndefInfoQ=False;

DefManifold[M,5,{a,b,c,d,e,f}]

DefChart[PS,M,{0,1,2,3,4},{t[],x[],y[],ϕ[],ψ[]}]

Formulas taken from [arXiv:0911.0802]:

DefConstantSymbol[{λ,ν,k}]

ds2=((2 H[x,y] (k^2) )/((((1-ν))^2) (((x-y))^2)))(((dx^2)/(G[x]))-((dy^2)/(G[y])))-2((J[x,y])/(H[y,x]))dϕ dψ-((H[y,x])/(H[x,y]))(((dt+Ω))^2)-((F[x,y])/(H[y,x])) ((dψ)^2)+((F[y,x])/(H[y,x])) ((dϕ)^2);

Ω:=Mf[x,y] dψ+P[x,y] dϕ;

Mf[x_,y_]:=((2 k λ (Sqrt[(((ν+1))^2)-(λ^2)]) (y+1) (-λ+ν-2 ν x+ν x ((λ+ν-1) x+2) y-1))/((1-λ+ν) H[y,x]));
P[x_,y_]:=((2 k λ (Sqrt[ν]) (Sqrt[(((ν+1))^2)-(λ^2)]) ((x^2)-1) y)/(H[y,x]));

H[x_,y_]:=(λ^2)+2 ν (1-(x^2)) y λ+2 x (1-(ν^2) (y^2)) λ-(ν^2)+ν (-(λ^2)-(ν^2)+1) (x^2) (y^2)+1;
F[x_,y_]:=((2 (k^2))/((((x-y))^2) (((1-ν))^2))) ((1-(y^2)) (((((1-ν))^2)-(λ^2)) (ν+1)+y λ (-(λ^2)-3 (ν^2)+2 ν+1)) G[x]+(-(1-ν) ν ((λ^2)+(ν^2)-1) (x^4)+λ (2 (ν^3)-3 (ν^2)-(λ^2)+1) (x^3)+((((1-ν))^2)-(λ^2)) (ν+1) (x^2)+λ ((λ^2)+(((1-ν))^2)) x+2 (λ^2)) G[y]);
J[x_,y_]:=((2 (k^2) (1-(x^2)) (1-(y^2)) λ (Sqrt[ν]) ((λ^2)+2 (x+y) ν λ-(ν^2)-x y ν (-(λ^2)-(ν^2)+1)+1))/((x-y) (((1-ν))^2)));
G[x_]:=(1-(x^2)) (ν (x^2)+λ x+1);

Computation of the metric tensor:

(matrix=CoefficientArrays[ds2,{dt,dx,dy,dϕ,dψ},"Symmetric"->True][[3]]//Normal//Simplify;)

A check:

ds2-{dt,dx,dy,dϕ,dψ}.matrix.{dt,dx,dy,dϕ,dψ}//Simplify//AbsoluteTiming

met=CTensor[matrix/.{x->x[],y->y[]},{-PS,-PS}];

met[-a,-b]

Simplify is way too slow for this metric. In this case it is more useful to simplify with Together:

$CVSimplify=Together;

Inverse has a hard time with this matrix. Specify explicitly the signature:

SetCMetric[met,PS,SignatureOfMetric->{4,1,0},Verbose->True]

Now we compute all the way down to the Ricci tensor.

MetricCompute[met,PS,"Ricci"[-1,-1],Verbose->True,Parallelize->True]

It is a Ricci-flat metric:

cd=CovDOfMetric[met];

Ricci[cd]

Let us compute the location of the event horizon. In a stationary BH the event horizon is a Killing horizon which is a set where the matrix formed by the scalar product of the Killing vectors is degenerate:

KillingMatrix={{met[{0,-PS},{0,-PS}],met[{0,-PS},{3,-PS}],met[{0,-PS},{4,-PS}]},{met[{3,-PS},{0,-PS}],met[{3,-PS},{3,-PS}],met[{3,-PS},{4,-PS}]},{met[{4,-PS},{0,-PS}],met[{4,-PS},{3,-PS}],met[{4,-PS},{4,-PS}]}};

Det@KillingMatrix;

Simplify@%

FullSimplify@%

The actual event horizon is among the zeros of this polynomial:

Last@%==0

Solve[%,y[]]//Simplify

This is the event horizon (the other root corresponds to a Cauchy horizon):

Last@%

eventhorizon=First@%

KillingMatrix/.eventhorizon;

KillingMatrixEH=(Simplify//@%);

The angular momenta parameters are solutions of the following equation:

MapThread[Equal,{KillingMatrixEH.{1,Jϕ,Jψ},{0,0,0}}];

%[[Range[1,2]]];

Solve[%,{Jϕ,Jψ}]//Simplify//Flatten;

angularmomenta=%;

The angular momenta are parameters that do not depend on the coordinate x. The only way of showing it is by computing the derivative and checking that it is zero:

angularmomenta[[All,-1]];

D[#,x[]]&/@angularmomenta[[All,-1]];

Simplify@%

Hence, the actual values of the angular momenta parameters are:

angularmomenta/.x[]->0

Simplify//@%

FullSimplify//@%

angularmomentavalues=%;

Remove[met,cd,dt,dx,dy,dϕ,dψ,ds2,F,G,H,J,P,Ω,matrix]

UndefConstantSymbol[{λ,ν,k}]

UndefChart[PS]

UndefManifold[M]

ClearxCobaCache[All]

Vector analysis in 3 dimensions

<<xAct`xCoba`

<<xAct`ShowTime1`

$DefInfoQ=False;

$CVVerbose=False;

The default simplification function is Simplify, but it may be a bit too slow with lots of radicals. We need some intermediate simplification because otherwise intermediate steps grow too fast. But we don’t need to lose too much time simplifying intermediate steps. Using Together is usually a good compromise:

$CVSimplify=Together;

Increase the size of objects displayed as CTensor components:

$LargeComponentSize=10000;

Define a 3-dimensional manifold with two charts (Cartesian and Spherical coordinates):

DefManifold[R3,3,{a,b,c,d,e,f,i,j,k}]

Spherical and Cartesian charts:

DefChart[sph,R3,{1,2,3},{r[],θ[],ϕ[]},ChartColor->Red]
DefChart[cart,R3,{1,2,3},{x[],y[],z[]},ChartColor->Green]

Change the default formatting:

sph/:CIndexForm[1,sph]:="r";
sph/:CIndexForm[2,sph]:="θ";
sph/:CIndexForm[3,sph]:="ϕ";
cart/:CIndexForm[1,cart]:="x";
cart/:CIndexForm[2,cart]:="y";
cart/:CIndexForm[3,cart]:="z";

SymbolOfCovD@PDcart^=List["\!\(\*StyleBox[\"¸\",FontColor->RGBColor[0., 1., 0.]]\)","\!\(\*StyleBox[\"\!\(\*StyleBox[\(\\[PartialD]\)]\)\",FontColor->RGBColor[0., 1., 0.]]\)"];

SymbolOfCovD@PDsph^=List["\!\(\*StyleBox[\"¸\",FontColor->RGBColor[1., 0., 0.]]\)","\!\(\*StyleBox[\"\!\(\*StyleBox[\(\\[PartialD]\)]\)\",FontColor->RGBColor[1., 0., 0.]]\)"];

The change of coordinates:

cart2sph={Sqrt[x[]^2+y[]^2+z[]^2],ArcCos[z[]/Sqrt[x[]^2+y[]^2+z[]^2]],ArcTan[x[],y[]]};
sph2cart=r[]{Sin[θ[]]Cos[ϕ[]],Sin[θ[]]Sin[ϕ[]],Cos[θ[]]};

cart2sphRules=Thread[{x[],y[],z[]}->sph2cart]

sph2cartRules=Thread[{r[],θ[],ϕ[]}->cart2sph]

Properties of the coordinates:

$Assumptions=And[
(x[]|y[]|z[]|r[]|θ[]|ϕ[])∈Reals,
r[]>0,
0<θ[]<Pi,
0≤ϕ[]<2Pi
];

We check that these rules are indeed correct

cart2sphRules/.sph2cartRules//Simplify

sph2cartRules/.cart2sphRules//FullSimplify

Compute Jacobian matrices:

jc2s=Simplify@Outer[D,cart2sph,{x[],y[],z[]}];
js2c=Simplify@Outer[D,sph2cart,{r[],θ[],ϕ[]}];

Transformation of the Jacobian matrices into a CTensor :(Important: beware of the order of indices e.g. {sph,-cart} and not {-cart,sph})

cart2sphJacobian=CTensor[jc2s,{sph,-cart}];
sph2cartJacobian=CTensor[js2c,{cart,-sph}];

{cart2sphJacobian[a,-b],sph2cartJacobian[a,-b]}

Check of the elementary properties:

(cart2sphJacobian[b,-a]sph2cartJacobian[c,-b]/.cart2sphRules)//Simplify

(sph2cartJacobian[b,-a]cart2sphJacobian[c,-b]/.sph2cartRules)//Simplify

Use SetBasisChange to set a change of basis. In the one-argument form it only pre-computes algebraic expressions, and not derivatives of the components of the Jacobian matrices or Christoffels:

SetBasisChange[cart2sphJacobian]

PrintxCobaCache[]

Consistency checks:

ToCTensor[Basis,{sph,-cart}]===cart2sphJacobian

ToCTensor[Basis,{cart,-sph}]===(sph2cartJacobian/.sph2cartRules//Simplify)

Determinants of the Jacobian matrices. Note that they are given in Cartesian coordinates because SetBasisChange used a CTensor written in Cartesian coordinates:

Jacobian[cart,sph][]//ToValues

Jacobian[sph,cart][]//ToValues

Use SetBasisChange with two arguments. The last argument selects the default coordinate chart. We now also compute derivatives, with respect to the coordinates of that given chart:

SetBasisChange[cart2sphJacobian,cart]

PrintxCobaCache[]

Orthonormal basis in spherical coordinates

We define an orthonormal basis related to spherical coordinates.

DefBasis[Osph,TangentR3,{1,2,3},BasisColor->Blue]

SetBasisChange[CTensor[{{1,0,0},{0,r[],0},{0,0,r[]Sin[θ[]]}},{-sph,Osph}],sph]

Basis[-{a,Osph},b]==TraceBasisDummy@SeparateBasis[sph]@Basis[-{a,Osph},b]

ToValues@ComponentArray@%

Metric

Introduce the Euclidean metric in the orthonormal basis :

gOsph=CTensor[DiagonalMatrix[{1,1,1}],{-Osph,-Osph},0]

gOsph[-i,-j]

In the spherical coordinates the metric is:

g=ToCTensor[gOsph,{-sph,-sph}]

We set this metric as our default metric:

SetCMetric[g,sph,SignatureOfMetric->{3,0,0}]

gcart=ToCTensor[g,{-cart,-cart}]/.cart2sphRules//Simplify

Connection

We compute the Levi-Civita connection of the metric:

CD=CovDOfMetric[g]

Curvature and Torsion

This computes all components of curvature tensors

MetricCompute[g,sph,All]//AbsoluteTiming

(riem=Simplify@Riemann[CD])[-a,-b,-c,d]

We check that there is torsion for the derivative associated with the orthonormal basis, but not for the derivative associated with the spherical coordinates

ToCTensor[TorsionPDOsph,{Osph,-Osph,-Osph}]
ToCTensor[TorsionPDsph,{sph,-sph,-sph}]

Write the Levi-Civita connection with respect to PDcart instead of PDsph.

CDcart=CovDOfMetric[ToCTensor[g,{-cart,-cart}]/.sph2cartRules//Simplify]

Some formulae from vector analysis

We would like to reproduce the usual formulas of vector calculus in spherical coordinates. We define a vector field and a scalar field:

DefTensor[S[],R3]
DefTensor[V[i],R3]

Associated CTensor vector, in different bases:

Vsph=ToCTensor[V,{sph},0]
Vcart=ToCTensor[V,{cart},0]

VOsph1=ToCTensor[V,{Osph},0]

Sometimes, even working in spherical coordinates we usually express results in terms of vector components in the orthonormal basis. Hence we transform the components inside VOsph1:

VOsph=ToCTensor[VOsph1,{sph}]

The gradient

With Cartesian coordinates:

CDcart[-i][S[]]

With spherical coordinates:

CD[-i][S[]]

Write components with respect to the orthonormal basis:

ToCTensor[Head@%,{-Osph}]

Compare with the classical formula:

Grad[S[r,θ,ϕ],{r,θ,ϕ},"Spherical"]

Divergence

Cartesian coordinates:

CDcart[-i]@Vcart[i]

Spherical coordinates:

CD[-i]@Vsph[i]

Use the components with respect to the orthonormal basis:

CD[-i]@VOsph[i]//Expand

Compare with the classical formula:

Div[{v1[r,θ,ϕ],v2[r,θ,ϕ],v3[r,θ,ϕ]},{r,θ,ϕ},"Spherical"]//Expand

Laplacian of scalars

Spherical coordinates:

CD[-i][S[]]

CD[-j][%]

Inv[g][i,j]%

%//Expand

Compare with the classical formula:

Laplacian[S[r,θ,ϕ],{r,θ,ϕ},"Spherical"]//Expand