Quit[]

Some general info and notation convention

Author: H. Witek
This notebook has been developed to perform the $3+1$ decomposition of Einstein’s equations, in the ``ADM-York’’ formulation. To do this calculation we use the xtensor package (J.M. Martin Garcia et al ‘02-‘21) for mathematica. The package can be downloaded here: http://www.xact.es/ This website also offers good documentation and example notebooks.
- xtensor uses the following notation for tensors: 
$X_{ab...}$: X[-a,-b,...]
$X^{ab...}$: X[a,b,...]
- For the metrics I will use the following notation:
g: spacetime metric (with signature (-+++) )
n: unit vector normal to the spatial hypersurfaces
γ: induced 3-metric 
- I use geometric units c=1, G=1

Loading xtensor

Let’s start by loading the xtensor package into mathematica:

<<xAct`xTensor`
<<xAct`ExpressionManipulation`

With the following specifications
(i)   the indices are output in a nice form,
(ii)  the sign convention for the extrinsic curvature is $K_{ab} = - 1/2 ℒ_{n} γ_{ab}$ 
(recall: $n$ denotes the time-like normal vector)
(iii) we enforce the sign convention for the acceleration that is usually used in NumRel

$PrePrint=ScreenDollarIndices;
$ExtrinsicKSign=-1;
$AccelerationSign=-1;

Definitions of manifold, 4- and 3-metric, tensors

Here we will define the 4-dimensional manifold, metric and tensors.

Definition of manifold and 4-metric

Here, we define the 4-dimensional spacetime {M, gg_{ab} } with manifold $M$ and metric $gg_{ab} $. 
The signature is (-,+,+,...,+) The covariant derivative with respect to $gg_{ab}$ is denoted as “CD” and (by default) printed as $\nabla$.
We also specify a list of characters used as indices throughout our computations.

Index=Union[Table[ToExpression["q"<>ToString[i]],{i,1,100}],IndexRange[a,m]];

DefManifold[M,4,Index](*Definition of manifold M, in D=4 *)
DefMetric[-1,gg[-a,-b],CD,PrintAs->"g"];(*definition of spacetime metric g*)

Define symbols of 4-tensors for the action and EoMs

Lagrangian and action

DefTensor[ℒEH4D[],M]
DefTensor[ℒTmn4D[],M]
DefTensor[𝒮EH4D[],M]

Note, that we will write the EoMs in the form Eq=0, which just makes the following treatment simpler.
First we define a set of tensor representing the EoMs in this form.
The gravitational coupling constant is kGR = 16 π G/c^4, and we use geometric units (G=1, c=1)

DefTensor[Tmn[-a,-b],M,Symmetric[{-a,-b}]](*energy-momentum tensor*)
DefTensor[EoM4DTen[-a,-b],M,Symmetric[{-a,-b}]](*tensor EoMs*)
DefTensor[EoM4DTenTR[-a,-b],M,Symmetric[{-a,-b}]](*trace-reversed form of tensor EoMs*)
DefTensor[EoM4DTenVac[-a,-b],M,Symmetric[{-a,-b}]](*tensor EoMs*)
DefTensor[EoM4DTenTRVac[-a,-b],M,Symmetric[{-a,-b}]](*trace-reversed form of tensor EoMs*)
DefConstantSymbol[Λ](*cosmological constant*)
DefConstantSymbol[kGR](*GR coupling to the energy-momentum tensor, give kGR=16π as rule*)
rulekGR=kGR:>16 π;
ruleΛ=Λ:>0;(*choice for asymptotically flat spacetimes*)

Define normal vector, 3-metric and lapse

We define the normal vector $n^a $ and specify the normalization $n_{a} n^{a} = -1$.

DefTensor[n[a],M]; (*timelike normal-vector*)
n/:n[a_] n[-a_]:=-1 (*normalization*)

We define the induced metric $γ_{ab}$ on the spatial slice $Σ_t$  which is induced from $gg_{ab}$ using the vector $n^a$ normal to the hypersurface. Because $γ_{ab}$ is a spatial metric the signature now is (+,+,...,+) . We denote the covariant derivative with respect to $γ_{ab}$ as “cd”, while throughout the notebook it will be shown as “D”.
The projection operator is $γ_{a}^{b} = δ_{a}^{b} + n_{a} n^{b} $. We define a rule which allows us to replace the projector symbol with the actual projection operator.

DefMetric[1,γ[-a,-b],cd,SymbolOfCovD->{"|","D"},InducedFrom->{gg,n},PrintAs->"γ"]
ruleProj=Projectorγ->ProjectWith[γ];

We next introduce the lapse function α and give a prescription how to substitute the acceleration vector with the lapse
a_{i} = \frac{1}{α} D_{i} α

DefTensor[α[],M]
Accelerationn/:Accelerationn[a_]:=cd[a][α[]]/α[]

Example for decomposition that will be used in the following sections

Any tensor can be decomposed into its spatial, normal and mixed components, as can be seen in this example:

DefTensor[TMP[-a,-b],M]

InducedDecomposition[TMP[-a,-b],{γ,n}](*/.ruleProj*)//NoScalar

Using the replacement rule defined above we can express the projection operator in terms of the induced metric and normal vector

InducedDecomposition[TMP[-a,-b],{γ,n}]/.ruleProj//NoScalar

In the following we will decompose the Einstein equations and energy momentum tensor in a similar manner.

Define parameters for explicit time evolution (lapse and shift)

Let’s perform a last change of variables to bring the time evolution equations in an explicit form. Therefore, remember that the time parameter $t$ is related to the lapse function α and shift vector β^i :
t^a = α u^a + β^a
Then the Lie derivative along u^a is given by 
ℒ_u = 1/α (\partial_t - ℒ_β )

DefParameter[t]
DefTensor[β[a],M,OrthogonalTo->{n[-a]},ProjectedWith->{γ[-a,b]}]

ruleLien=LieD[n[c_]][X_]:>1/α[] (LieD[t][X]-LieD[β[c]][X])

Action & eoms (4D)

Substitution and simplification rules

Auxiliary tensors

DefTensor[𝒮4Dstep1[],M]
DefTensor[EoM4Dstep1[],M]

Recover Ricci tensor

Define some rules to recover R_{ab} from derivative of Christoffels

ruleGam2Ric={PD[-e_][ChristoffelCD[e_,-a_,-b_]]:>Module[{c,d},RicciCD[-a,-b]+PD[-a][ChristoffelCD[c,-b,-c]]+ChristoffelCD[c,-a,-d]ChristoffelCD[d,-b,-c]-ChristoffelCD[c,-c,-d]ChristoffelCD[d,-a,-b]],
CD[-e_][ChristoffelCD[e_,-a_,-b_]]:>Module[{c,d},RicciCD[-a,-b]+CD[-a][ChristoffelCD[c,-b,-c]]-ChristoffelCD[c,-a,-d]ChristoffelCD[d,-b,-c]+ChristoffelCD[c,-c,-d]ChristoffelCD[d,-a,-b]]};

ruleGam2Ric2={PD[-e_][ChristoffelCD[e_,a_,b_]]:>Module[{c,d},RicciCD[a,b]+PD[a][ChristoffelCD[c,b,-c]]+ChristoffelCD[c,a,-d]ChristoffelCD[d,b,-c]-ChristoffelCD[c,-c,-d]ChristoffelCD[d,a,b]],
CD[-e_][ChristoffelCD[e_,a_,b_]]:>Module[{c,d},RicciCD[a,b]+CD[a][ChristoffelCD[c,b,-c]]-ChristoffelCD[c,a,-d]ChristoffelCD[d,b,-c]+ChristoffelCD[c,-c,-d]ChristoffelCD[d,a,b]]};

ruleGamSym={CD[-a_][ChristoffelCD[d_,-b_,-d_]]-CD[-b_][ChristoffelCD[d_,-a_,-d_]]:>0,
PD[-a_][ChristoffelCD[d_,-b_,-d_]]-PD[-b_][ChristoffelCD[d_,-a_,-d_]]:>0};

Define Lagrangians and action

Einstein-Hilbert Lagrangian including the cosmological constant

Here we define the Einstein-Hilbert Lagrangian

ℒEH4D[]=RicciScalarCD[]-2Λ

Action (remember notation: w/o integral)

NOTE: in xtensor the action is represented w/o the integral (this is just notation, so remember that it is actually there!)
Here, we assume a vacuum spacetime. You can, of course, add matter in the usual manner

(*𝒮EH4D[]=Sqrt[-Detgg[]](1/kGR ℒEH4D[]+ℒTmn4D[])*)
𝒮EH4D[]=Sqrt[-Detgg[]] ℒEH4D[]

Derive tensor  EoMs

We will make use of Einstein’s eqs in standard and trace-reversed form:
(i)  EoMT_{ab}:=R_{ab} - 1/2 g_{ab} R + Λ g_{ab} - kGR T_{ab} = 0
(ii) EoMTtr_{ab}:=R_{ab} -Λ g_{ab} - kGR (T_{ab} - 1/(D-2) g_{ab} T ) = 0 
By varying the action above (recall that the integral is define implicitly) we can derive the eoms:

Since we do not have an explicit form for the matter Lagrangian, we set it to 0 here, and add the general energy-momentum tensor afterwards.

(*𝒮4Dstep1[]=kGR 𝒮EH4D[]/.ℒTmn4D[]:>0//Expand//RiemannToChristoffel//NoScalar*)
𝒮4Dstep1[]=𝒮EH4D[]//RiemannToChristoffel//NoScalar

EoM4Dstep1[-a_,-b_]=1/Sqrt[-Detgg[]]VarD[gg[a,b]][𝒮4Dstep1[]]//Expand

EoM4DTenVac[-a_,-b_]=(EoM4Dstep1[-a,-b]/.ruleGam2Ric//Simplification)/.ruleGamSym//Simplification//ToCanonical

Now we add the energy momentum tensor again. Recall, that we write the eoms in the form
EoMT_{ab} =G_{ab} + Λ g_{ab} - kGR T_{ab} = 0

EoM4DTen[-a_,-b_]=EoM4DTenVac[-a,-b]-kGR/2 Tmn[-a,-b]/.rulekGR

Trace-reversed form (used for derivation of dt K_{ij})

EoM4DTenTR[-a_,-b_]=Module[{c,d},Collect[EoM4DTen[-a,-b]-1/2gg[-a,-b]gg[c,d]EoM4DTen[-c,-d]//Simplification//ContractMetric//ToCanonical, π]]

Projections of the energy-momentum tensor

Here, we decompose the energy-momentum tensor and define its the various projections
(a) energy density:             ρ = T_{ab} u^{a} u^{b}
(b) energy current:            j_i = -T_{ab} u^a h^{b}_{i}
(c) spatial components:  S_{ij} = h^a_i h^b_j T_{ab}

Define tensors

We define the decomposed energy momentum tensor $DTmn$, the energy density  $ρ$ and flux $j_{i}$ and spatial components $S_{ij}$ of the energy momentum tensor.

DefTensor[DTmn[-a,-b],M,Symmetric[{-a,-b}]]
DefTensor[DTmnstep1[-a,-b],M,Symmetric[{-a,-b}]]
DefTensor[ρ[],M]
DefTensor[ji[-a],M,OrthogonalTo->{n[a]},ProjectedWith->{γ[a,-c]} ]
DefTensor[Sij[-a,-b],M,Symmetric[{-a,-b}],OrthogonalTo->{n[a],n[b]},ProjectedWith->{γ[a,-c],γ[b,-d]} ]

Rules: Tab to ρ,  ji, Sij

We define some substitution rules to rewrite projections & contractions of the energy-momentum tensor in terms of the energy density, current, etc

rule2ρ=Tmn[-a_,-b_]n[a_]n[b_]:>ρ[];
rule2ji={γ[-a_,c_]n[b_]Tmn[-b_,-c_]:>-ji[-a],γ[-a_,c_]n[b_]Tmn[-c_,-b_]:>-ji[-a]};
rule2Sij=γ[-a_,c_]γ[-b_,d_]Tmn[-c_,-d_]:>Sij[-a,-b];
ruleTrTmn=Tmn[a_,-a_]:>Module[{b},Sij[b,-b]-ρ[]];

We can convince ourselves that the last relation holds

gg[a,b](InducedDecomposition[Tmn[-a,-b],{γ,n}]/.ruleProj//NoScalar)/.rule2ρ/.rule2ji/.rule2Sij//MetricToProjector//Simplification

Decomposition

Decomposition of the energy-momentum tensor (as in general example above)

DTmnstep1[-a_,-b_]=InducedDecomposition[Tmn[-a,-b],{γ,n}]/.ruleProj//NoScalar

Let’s replace the projections & contractions of the energy-momentum tensor with the energy density etc

DTmn[-a_,-b_]=DTmnstep1[-a,-b]/.rule2ρ/.rule2ji/.rule2Sij

Decomposition of tensor EoMs

Now, we perform the 3+1 split of the tensor EoMs. The Hamiltonian and momentum constraints as well as the time evolution eq. for the extrinsic curvature are obtained from the appropriate contractions along the normal vector and projections onto the hypersurface.
Conveniently, xtensor already knows the Gauss-Codazzi relations for the Riemann and Ricci tensor, which we can easily apply.

Define tensors and subsitution rules

First, we introduce two additional tensor fields which will denote the decomposed equations of motion, again in ``standard’’ and trace-reversed form

DefTensor[DEoM4DTen[-a,-b],M,Symmetric[{-a,-b}]]
DefTensor[DEoM4DTenstep1[-a,-b],M,Symmetric[{-a,-b}]]
DefTensor[DEoM4DTenTR[-a,-b],M,Symmetric[{-a,-b}]]

Here we define the tensors for the normal contraction (-> Hamiltonian), mixed projection (-> momentum constraint) and spatial projection (in the notation $𝒫_{ab}=0$). We also define the trace and trace-free part of the spatial projection

DefTensor[ℋ[],M]
DefTensor[ℋstep1[],M]
DefTensor[ℋstep2[],M]
DefTensor[ℳ[-a],M,OrthogonalTo->{n[a]},ProjectedWith->{γ[a,-b]} ]
DefTensor[𝒫[-a,-b],M,Symmetric[{-a,-b}],OrthogonalTo->{n[a],n[b]},ProjectedWith->{γ[a,-c],γ[b,-d]} ]
DefTensor[𝒫tr[],M]
DefTensor[𝒫tf[-a,-b],M,Symmetric[{-a,-b}],OrthogonalTo->{n[a],n[b]},ProjectedWith->{γ[a,-c],γ[b,-d]} ]

definitions for the time evol eqs of the physical 3-metric $h_{ab}$ and extrinsic curvature $K_{ab}$

DefTensor[tmpdtγ[-a,-b],M,Symmetric[{-a,-b}]]
DefTensor[dtγ[-a,-b],M,Symmetric[{-a,-b}]]
DefTensor[dtKij[-a,-b],M,Symmetric[{-a,-b}]]

Simplification and substitution rules

These are some simplification rules that we will need

ruleSimADM1={cd[-a_][cd[a_][α[]]]:>Module[{b,c},γ[b,c]cd[-b][cd[-c][α[]]]],cd[a_][cd[-a_][α[]]]:>Module[{b,c},γ[b,c]cd[-b][cd[-c][α[]]]]};
ruleSimADM2=γ[a_,b_]CD[c_][ExtrinsicKγ[-a_,-b_]]:>Module[{d,e},CD[c][ExtrinsicKγ[d,-d]]-ExtrinsicKγ[-d,-e]CD[c][γ[d,e]]];
ruleSimADM3=γ[-a_,d_]γ[-b_,e_]CD[c_][ExtrinsicKγ[-d_,-e_]]:>Module[{f,g},CD[c][ExtrinsicKγ[-a,-b]]-ExtrinsicKγ[-f,-g]CD[c][γ[-a,f]γ[-b,g]]];
ruleSimADM4=γ[a_,b_]LieD[n[c_]][ExtrinsicKγ[-a_,-b_]]:>Module[{d,e},LieD[n[c]][ExtrinsicKγ[d,-d]]-ExtrinsicKγ[-d,-e]LieD[n[c]][γ[d,e]]];

Geometric relation -> evol eq for metric

tmpdtγ[-a_,-b_]=Module[{c},LieD[n[c]][γ[-a,-b]]-1/α[](LieD[t][γ[-a,-b]]-LieD[β[c]][γ[-a,-b]])//Simplification]

ruledtγ=Solve[tmpdtγ[-a,-b]==0,LieD[t][γ[-a,-b]]][[1,1]];
dtγ[-a_,-b_]=LieD[t][γ[-a,-b]]/.ruledtγ

Decomposition of EoMs

Decompose tensor for eoms and re-express contractions/projections of energy momentum tensor

DEoM4DTenstep1[-a_,-b_]=(InducedDecomposition[EoM4DTen[-a,-b],{γ,n}]/.ruleProj//NoScalar//MetricToProjector)/.Tmn[a_,b_]:>0(*/.rule2ρ/.rule2ji/.rule2Sij*)//ToCanonical

Apply the Gauss-Codazzi relations, and some further simplifications

DEoM4DTen[-a_,-b_]=(GaussCodazzi[DEoM4DTenstep1[-a,-b]]/.ruleProj//NoScalar//MetricToProjector//GradNormalToExtrinsicK//Simplification//ToCanonical)/.ruleSimADM1/.ruleSimADM2//Simplification//ToCanonical

We also decompose the trace-reversed version of the eoms. Those will be useful to obtain $\p_{t} K_{ij}$

DEoM4DTenTR[-a_,-b_]=(GaussCodazzi[(InducedDecomposition[EoM4DTenTR[-a,-b],{γ,n}]/.ruleProj//NoScalar//MetricToProjector)/.Tmn[a_,b_]:>0//Simplification//ToCanonical]/.ruleProj//NoScalar//MetricToProjector//GradNormalToExtrinsicK//Simplification//ToCanonical)/.ruleSimADM1/.ruleSimADM2//Simplification//ToCanonical;

Normal projection -> Hamiltonian constraint

Recall: ℋ = 2 E_{ab} n^{a} n^{b}

ℋstep1[]=Module[{a,b},2n [a]n[b]DEoM4DTen[-a,-b]//Simplification//ToCanonical]

ℋstep2[]=ToInducedDerivative[ℋstep1[],CD,cd]//Simplification

ℋ[]=ℋstep2[]//GradNormalToExtrinsicK//Simplification//ToCanonical

Mixed projection -> momentum constraint

Recall: ℳ_{a} = - γ^{b}{}_{a} n^{c} E_{bc}

ℳ[-a_]=Module[{b,c},-γ[-a,b]n[c]DEoM4DTen[-b,-c]//Simplification//ToCanonical]

Spatial projection -> dt Kij

Here we take the fully spatial projection of the trace-reversed eoms, i.e., γ^{c}_{a} γ^{d}_{b} Etr_{cd}

𝒫[-a_,-b_]=Module[{c,d},ToInducedDerivative[(γ[-a,c]γ[-b,d]DEoM4DTenTR[-c,-d]//Simplification)/.ruleSimADM3//Simplification,CD,cd]//GradNormalToExtrinsicK//Simplification//ContractMetric//ToCanonical]

ruledtKij=Solve[(𝒫[-a,-b]/.ruleLien)==0,LieD[t][ExtrinsicKγ[-a_,-b_]]][[1,1]];
dtKij[-a_,-b_]=Collect[LieD[t][ExtrinsicKγ[-a,-b]]/.ruledtKij//Simplification//ContractMetric//ToCanonical,π α[]]

Take trace and tracefree part

𝒫tr[]=Module[{c,d},(γ[c,d]𝒫[-c,-d]//Simplification//ToCanonical)/.ruleSimADM4//Simplification//ContractMetric//ToCanonical]

𝒫tf[-a_,-b_]=Module[{c,d},(𝒫[-a,-b]-1/3γ[-a,-b]γ[c,d]𝒫[-c,-d]//Simplification)/.ruleSimADM4//Simplification//ContractMetric//ToCanonical]

Output as tex files

Now we will use the package TexAct (part of the xAct bundle) to export our results as tex-files

<<xAct`TexAct`

NOTE: adjust the folder names to your local paths

constraints:

Export["/Users/helviwitek/Teaching/ETK_SummerSchool_UIUC_July2021/GR_Split/HamADM.tex",TexPrint[ℋ[]],"Text"];
Export["/Users/helviwitek/Teaching/ETK_SummerSchool_UIUC_July2021/GR_Split/MomADM.tex",TexPrint[ℳ[-i]],"Text"];

evolution equations for metric and extrinsic curvature

Export["/Users/helviwitek/Teaching/ETK_SummerSchool_UIUC_July2021/GR_Split/dtgammaADM.tex",TexPrint[dtγ[-i,-j]],"Text"];
Export["/Users/helviwitek/Teaching/ETK_SummerSchool_UIUC_July2021/GR_Split/dtKijADM.tex",TexPrint[dtKij[-i,-j]],"Text"];