Quit

xAct - The Basics

Jeremy Sakstein

This is a basic introduction to the Mathematica package xAct, which can be downloaded here: http://www.xact.es/index.html. There are also installation instructions and explanations for some of the more advanced features.

Please note that the contents of this notebook are formatted in large fonts so that they are easily viewed during a live demonstration.

This latest version of this tutorial can be found at www.jeremysakstein.com.

1) The Basics

In this section I will go through the basics of using the package for tensor manipulations. xAct actually has lots of sub packages and the one that handles tensor manipulations is xTensor. The first thing we need to do is to load it:

<<xAct`xTensor`

The first thing we need to do is define a manifold, which we will call M4:

DefManifold[M4,4,{a,b,c,d,e,f,g,h,i,j,k,l}]

This defines a manifold with four dimensions whose tensors carry indices a-l. If you are doing very complicated computations you may want to add more indices. Next, we need to define a metric, which we define using this command:

DefMetric[-1,metric[-a,-b],CD,SymbolOfCovD->{"|","∇"},PrintAs->"g"]

Note the use of the - signs; in xAct, normal letters denote upper indices and negative signs before them denote lower indices. The specification CD tells xAct what command to use for the covariant derivative and we have told it to print the metric using the standard g notation.  As you can see, xAct has also defined other quantities related to the metric such as the Riemann tensor, the connection coefficients and the Ricci scalar. Let’s have a look at what xAct actually does with these things. For example, the Riemann Tensor:

RiemannCD[a,-b,-c,-d]//RiemannToChristoffel

Oh dear! This looks like a right mess. The reason for this is that xTensor uses internal Mathematica indices to do the manipulations. To make this look nice, we use the following command:

$Pre=ScreenDollarIndices;

You also need to tell xTensor to allow upper indices on covariant derivatives

SetOptions[ContractMetric,AllowUpperDerivatives->True];
$CovDFormat= "Prefix";

Let’s try it now

RiemannCD[a,-b,-c,-d]//RiemannToChristoffel

Ah, much better. Let’s go ahead and play around a bit. Our metric is used like this:

metric[-a,-b]

Or, we can use the inverse

metric[a,b]

These should contract to give us a delta-symbol

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

We can also contract other things, like the Ricci tensor

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

Of course, there is no need for this, xTensor already knows what the Ricci scalar is

RicciScalarCD[]

Note you can always convert between these compact expressions and the connection coefficients using the RiemanntoChristoffel command

RicciScalarCD[]//RiemannToChristoffel

The red brackets tell is that xTensor knows that this object is a scalar. Although this is largely aesthetic, there are sometimes circumstances where we removing these can halp with xTensor’s simplifications. The brackets can be removed as follows:

RicciScalarCD[]//RiemannToChristoffel//NoScalar

They can be put back using the PutScalar Command, which can be useful to identify scalar quantities in long expressions

RicciScalarCD[]//RiemannToChristoffel//NoScalar//PutScalar

There is also EinsteinToRicci command that is similar to the RiemannToChristoffel function:

EinsteinCD[-a,-b]

%//EinsteinToRicci

Now let’s do something that is actually useful. Let’s take the Jordan-frame action ∫(d^4) x (Sqrt[-(OverTilde[g])])[(1/2) (Ω^2)(ϕ)(OverTilde[R])-(1/2)h(ϕ)(Subscript[(OverTilde[∇]), μ])ϕ((OverTilde[∇])^μ)ϕ-(Ω^4)(ϕ)V(ϕ)] with h (ϕ) = (Ω^2) (ϕ) (1 - (((3/2)[((d ((lnΩ)^2))/(d ϕ))])^2)) to the Einstein frame. This is done using the Weyl-rescaling (Subscript[(OverTilde[g]), ab=])(Ω^(-2))(Subscript[g, ab]). Following Wald, appendix C, we have (Subscript[((OverTilde[R])^a), bcd]) = (Subscript[(R^a), bcd]) + 2 (Subscript[∇, [c])(Subscript[(K^a), d]b]) +(Subscript[(K^a), e[c])(Subscript[(K^e), d]b]), with (Subscript[(K^a), bc]) =(Subscript[((OverTilde[Γ])^a), bc])-(Subscript[(Γ^a), bc]). Note that I am following 0905.2966 but I have swapped the tildes for each frame. Let’s start by defining the Jordan-frame metric:

DefMetric[-1,metricJ[-a,-b],CDJ,SymbolOfCovD->{"|","D"},PrintAs->"(OverTilde[g])"]

Mathematica warns you that there are other metrics on the manifold but this isn’t a problem computationally. Next, we need to define the scalar ϕ, which we do using the DefTensor command

DefTensor[Phi[],M4,PrintAs->"ϕ"]

This defines a tensor phi and the [] tells xTensor that it is a rank-zero tensor. We can define scalar functions as follows:

DefScalarFunction[Ω]

We will also need to define K

DefTensor[K[a,-b,-c],M4,Symmetric[{-b,-c}]]

Note that K is symmetric and we can tell xTensor this in order to aid in the computations. We can set the indices of one tensor in terms of another using the IndexSet function

IndexSet[metricJ[-a_,-b_],Ω[Phi[]]^(-2)metric[-a,-b]]

Note the use of the _’s. Without this, the index setting would only hold if you use a AND b. We want it to hold for any letters. Because we are using a second metric, we need to define the inverse too:

IndexSet[InvmetricJ[a_,b_],Ω[Phi[]]^2metric[a,b]]

These still give us the delta function:

metricJ[-a,-b]InvmetricJ[b,c]

Now let’s set K:

IndexSet[K[a_,-b_,-c_],ChristoffelCDJ[a,-b,-c]-ChristoffelCD[a,-b,-c]]

We can find an explicit expression for this. First, we change all covariant derivatives to partial derivatives

K[a,-b,-c]//ChristoffelToMetric

Next, we change all partials (PD) to Einstein-frame covariant derivatives using the ChangeCovD command

ChangeCovD[%,PD,CD]//ChristoffelToMetric

Urgh! What a mess. You need to tell xTensor to simplify things using the ToCanonical command. This tell sit to use any symmetry it knows or you tell it to simplify things as much as possible

%//ToCanonical

And there you have it! As someone who’s worked on chameleons, I’ve derived this by hand many times and it always took longer than this! The final thing we need to do is take covariant derivatives. This is done as so:

CD[a][Phi[]]

Now let’s find the Einstein-frame Riemann tensor. It isn’t a good idea to IndexSet xTensors predefined tensors so let’s use any old name for now. xTensor has its own Antisymmetrize function:

RiemannJF=ChangeCovD[RiemannCD[a,-b,-c,-d]+2*Antisymmetrize[CD[-c][K[a,-d,-b]],{-c,-b}]+2*Antisymmetrize[K[a,-e,-c]K[e,-d,-b],{-c,-d}]//ChristoffelToMetric,PD,CD]//ChristoffelToMetric//ToCanonical

Note that I have taken the same steps as before to get it into a nice form. Next we find the Ricci tensor:

RicciTJF=delta[-a,c]RiemannJF

Note that some terms are of the same form but the metric isn’t contracted. With longer expressions, it is best to tell xTensor to contract the metrics explicitly

RicciTJF=delta[-a,c]RiemannJF//ContractMetric

The Ricci scalar is then

RicciJF=InvmetricJ[b,d]RicciTJF//ContractMetric

As you can see, there are several terms that are identical because xTensor han’t used the fact that contracted dummy indices can be replaced by any letter you like. This is taken care of by the ToCanonical command

RicciJF=RicciJF//ToCanonical

but there is also a separate command called ReplaceDummies that is sometimes useful for more complicated expressions

RicciJF=RicciJF//ReplaceDummies

There,  nice and simple. Now let’s define h and then the action:

DefScalarFunction[hh,PrintAs->"h"]

hh[a_]:=Ω[a]^2(1-3/2*(D[Log[Ω[a]^2],a])^2)

hh[Phi[]]

xTensor can’t do the metric determinant so this needs to be done by hand

metricdetJF=Ω[Phi[]]^(-4)*Sqrt[-Detmetric[]]

Don’t worry about the tildes - this is just how xTensor prints determinants:

DetmetricJ[]

Detmetric[]

We need to define the potential

DefScalarFunction[V]

Finally, the Lagrangian is

metricdetJF*(Ω[Phi[]]^2RicciJF/2-InvmetricJ[a,b]CDJ[-a][Phi[]]CDJ[-b][Phi[]]*hh[Phi[]]/2-Ω[Phi[]]^(4)V[Phi[]])

Note that we need to tell xTensor to contract the Jordan-frame covariant derivatives using the Jordan-frame metric. We still need to tidy this up. First we get every covariant derivative into the Einstein-frame

LEF=Collect[ChangeCovD[%,CDJ,CD]//ToCanonical//ContractMetric,Sqrt[-Detmetric[]]]

The last three terms are actually a total derivative:

TD=CD[-a][CD[a][Log[Ω[Phi[]]]]]//ToCanonical

The final action is then

LEF=LEF-3Sqrt[-Detmetric[]]*TD//ToCanonical

After all that, this is simply the action for a canonically normalised scalar field coupled to gravity.

2) Equations of Motion

xAct knows how to do variational derivatives but we need to use a different package: xPert

<<xAct`xPert`

The colouring in this is atrocious so we use this to tidy it up

Unprotect[IndexForm];
IndexForm[LI[x_]]:=ColorString[ToString[x],RGBColor[0,0,1]];
Protect[IndexForm];

We first need to define a metric perturbation, which we’ll call h. This command defines a metric perturbation metpert and to label to order of the perturbation by a perturbation parameter epsilon.

DefMetricPerturbation[metric,metpert,ϵ];
PrintAs[metpert]^="h";

Perturbed quantities can be written as follows:

metpert[LI[4],-a,-b]

This is the fourth-order metric perturbation. Other examples are

Perturbation[metric[-a,-b],1]
Perturbation[metric[a,b],3]

The first one went through okay but the second one caused problems. This is xPert telling us that it doesn’t know how to perturb this so we need to use the ExpandPerturbation command

ExpandPerturbation@Perturbation[metric[a,b],3]

This is the third order perturbation of the inverse metric. We can collect perturbations order by order in the expansion parameter using the Perturbed command:

Perturbed[metric[a,b],4]

ExpandPerturbation@%

Now let’s use this to find the equations of motion by varying the action. We start by defining the perturbation to the scalar field:

DefTensorPerturbation[PertPhi[LI[order]],Phi[],M4,PrintAs->"δφ"]

I’m also going to define the Planck mass as a constant

DefConstantSymbol[mp,PrintAs->"Mp"]

Let’s do the simplest thing possible, general relativity

LGR=Sqrt[-Detmetric[]]mp^2/2*RicciScalarCD[]

Now we perturb this to first order in the metric perturbations

LpertGR=ToCanonical@ContractMetric@ExpandPerturbation@Perturbation@LGR

Now we take a variational derivative with respect to h1:

0==-VarD[metpert[LI[1],a,b],CD][LpertGR]/Sqrt[-Detmetric[]]/.delta[-LI[1],LI[1]]->1//SeparateMetric[metric]//RicciToEinstein//Expand//ContractMetric//ToCanonical

And there we have it, Einstein’s equation. Now let’s do something harder. Let’s define some scalar functions

DefScalarFunction[A]
DefScalarFunction[B]

and some scalar-tensor theories:

LSTmin=LGR+Sqrt[-Detmetric[]]*(-1/2CD[a][Phi[]]CD[-a][Phi[]]-V[Phi[]])

LSTnmin=A[Phi[]]^(-2)LGR+Sqrt[-Detmetric[]]*(-B[Phi[]]CD[a][Phi[]]CD[-a][Phi[]]-V[Phi[]])

The first is a minimally coupled quintessence model and the second a scalar-tensor theory with a non-canonical kinetic term. Let’s find the scalar EOM for both theories

LpertSTmin=ToCanonical@ContractMetric@ExpandPerturbation@Perturbation@LSTmin
LpertSTnmin=ToCanonical@ContractMetric@ExpandPerturbation@Perturbation@LSTnmin

KG=0==VarD[PertPhi[LI[1]],CD][LpertSTmin]/Sqrt[-Detmetric[]]/.delta[-LI[1],LI[1]]->1//ToCanonical//Simplify

The minimally coupled model simply gives us the Klein-Gordon equation.

SEOMnmin=0==VarD[PertPhi[LI[1]],CD][LpertSTnmin]/Sqrt[-Detmetric[]]/.delta[-LI[1],LI[1]]->1//ToCanonical//Simplify

And here is what the non-minimal theory gives us. Now usually we would want to eliminate the Ricci scalar, so let’s find the EOM for the tensor:

0==-VarD[metpert[LI[1],a,b],CD][LpertSTmin]/Sqrt[-Detmetric[]]/.delta[-LI[1],LI[1]]->1//SeparateMetric[metric]//RicciToEinstein//Expand//ContractMetric//ToCanonical

This is just Einstein’s equations sourced by the energy-momentum tensor for a scalar field.

TensEOMnmin={0==-VarD[metpert[LI[1],a,b],CD][LpertSTnmin]/Sqrt[-Detmetric[]]/.delta[-LI[1],LI[1]]->1//SeparateMetric[metric]//RicciToEinstein//Expand//ContractMetric//ToCanonical}

The non-minimal coupling here has given us a horrible mess. Let’s take the trace of this equation:

Tracenmin=metric[a,b]TensEOMnmin//ContractMetric//EinsteinToRicci//ToCanonical

We can now solve this for the Ricci scalar:

RS=Solve[Tracenmin,RicciScalarCD[]][[1]]

We can now eliminate this in the scalar EOM:

SEOMnmin/.RS//ToCanonical

And there we have it, the full scalar EOM with the Ricci scalar eliminated.

xAct is also very useful for dealing with higher-derivative theories such as Galileons or Horndeski Lagrangians. Let’s do an example that encompasses both, the cubic and quartic covariant galileon:

DefConstantSymbol[L3,PrintAs->"(Subscript[Λ, 3])"]
DefConstantSymbol[L4,PrintAs->"(Subscript[Λ, 4])"]

It is often easier to define X=-1/2(Subscript[∇, μ])ϕ(∇^μ)ϕ to make things simpler

DefTensor[X[],M4,PrintAs->"X"]

IndexSet[X[],-1/2CD[a][Phi[]]CD[-a][Phi[]]]

LGal=Sqrt[-Detmetric[]]*(X[]+1/(2*L3^3)*X[]*CD[-b][CD[b][Phi[]]]-(2X[]/L4^6)*(X[]*RicciScalarCD[]+2*(CD[-a][CD[a][Phi[]]]CD[-b][CD[b][Phi[]]]-CD[-a][CD[-b][Phi[]]]CD[a][CD[b][Phi[]]]))+mp^2/2*RicciScalarCD[])

Note that this is a Horndeski Lagrangian with G2=X, G3=-X and G4=1-2(X^2). We can now find the scalar and tensor EOM easily:

LGalpert=ToCanonical@ContractMetric@ExpandPerturbation@Perturbation@LGal

tenseomgal=0==-VarD[metpert[LI[1],a,b],CD][LGalpert]/Sqrt[-Detmetric[]]/.delta[-LI[1],LI[1]]->1//SeparateMetric[metric]//RicciToEinstein//Expand//ContractMetric//ToCanonical

SEOMgal=VarD[PertPhi[LI[1]],CD][LGalpert]/Sqrt[-Detmetric[]]/.delta[-LI[1],LI[1]]->1//ToCanonical

Hang on, what’s going on here, I thought the Horndeski Lagrangian gives second order field equations? It does. The problem is that different higher-order terms have their indices on the covariant derivatives in different orders and we need to commute them past each other, picking up powers of the Riemann tensor as we go, in order to cancel them explicitly and make the second-order nature manifest. We can do this using the SortCovDs command. This tells xAct to commute covariant derivatives, introducing geometric tensors if necessary, until all of the indices are in the same order:

SEOMgal2=ContractMetric[SortCovDs[SEOMgal,CD],AllowUpperDerivatives->True]//ReplaceDummies//ToCanonical

Notice that there are now no scalar derivatives that are more than second-order. What about the metric? There are still terms that look like (Subscript[∇, b])R, which contain third order derivatives of the metric so why are these here? Actually they aren’t, they cancel against each other once the Bianchi identities are applied but Mathematica doesn’t know this so we need to impose the identities by hand. This can be done using the following rule:

BianchiID={CD[-a_]@Phi[]*CD[-b_]@RicciCD[b_,a_]:>1/2 CD[-a]@Phi[] metric[b,a] CD[-b]@RicciScalarCD[],CD[a_][Phi[]]*CD[-b_][RicciCD[-a_,b_]]:>1/2 CD[a]@Phi[]*CD[-a]@RicciScalarCD[]}

SEOM3=SEOMgal2/.BianchiID

There we go, the full equation of motion for the quartic galileon in a curved space-time. Now you could go further and remove the factors of the Ricci tensor and scalar and the Riemann tensor using Einstein’s equation, but I’ll leave that as an exercise.

3) Coordinate Systems

Next, we will have a look at how to implement coordinate systems. For this, we need to use another package called xCobra:

<<xAct`xCoba`

If we want to use coordinates we need a chart, which we’ll call BA

DefChart[BA,M4,{0,1,2,3},{t[],r[],θ[],ϕ[]}]

Let’s work with a Schwarzschild-like metric for now. This will let me show you the power of xAct and we can specify to flat-space later.

MatrixForm[metricarray={{-F[r[]],0,0,0},{0,1/F[r[]],0,0},{0,0,r[]^2,0},{0,0,0,r[]^2 Sin[θ[]]^2}}]

Next, we need to tell xAct to put the metric into this basis

MetricInBasis[metric,-BA,metricarray]

And finally, we need to compute things like the Ricci scalar:

MetricCompute[metric,BA,All,CVSimplify->Simplify]

Now we can calculate nice quantities e.g.

RicciScalarCD[]//ToBasis[BA]//ToValues

We can calculate this for flat space

F[x_]:=1
RicciScalarCD[]//ToBasis[BA]//ToValues

Or a black hole

DefConstantSymbol[rs,PrintAs->"(Subscript[r, s])"]

F[x_]:=1-rs/x

RicciScalarCD[]//ToBasis[BA]//ToValues

Sometimes when you do this Mathematica won’t cancel everything so we need to tell it to do it

%//Simplify

And of course the answer is zero. Another less trivial example is the Gauss-Bonnet term:

GB=RiemannCD[a,b,c,d]RiemannCD[-a,-b,-c,-d]-4RicciCD[-a,-b]RicciCD[a,b]+RicciScalarCD[]^2

In flat space we have

F[x_]:=1
GB//ToBasis[BA]//ToValues

What’s happened? If there are lots of indices involved you need to tell xAct to trace over them:

GB//ToBasis[BA]//TraceBasisDummy//ToValues

And for Schwarzschild

F[x_]:=1-rs/x
GB//ToBasis[BA]//TraceBasisDummy//ToValues//Simplify

Note that this time it still didn’t work even after we traced over the dummy indices. The truth is that xCobra is very temperamental and you sometimes have to play around to get it to work. xCobra defines quantities in terms of other quantities and sometimes these aren’t traced properly so it’s best to separate the metric before tracing:

GB//SeparateMetric[metric]//ToBasis[BA]//TraceBasisDummy//ToValues//Simplify

This takes a while, but it’s done it. Now let’s look at an example involving a scalar. We’ll start with the Klein-Gordon equation in flat space when the field is radially-varying only

Clear[F]
KGEOM=CD[-a][CD[a][Phi[]]]-V[Phi[]]

We need to define a scalar function to represent ϕ when it depends on the coordinates

DefScalarFunction[ϕ0]

ComponentValue[Phi[],ϕ0[r[]]]

KGEOM//SeparateMetric[metric]//ToBasis[BA]//TraceBasisDummy//ToValues

And there you have it. It’s straightforward to apply this to more complicated equations, you just copy and paste.

4) Relativistic Stars

Now that we have all the basics we can look at a solid example: relativistic stars. I’ll begin by quitting the kernel to clear our old coordinate system:

Quit

Next, we set up the problem, including loading the packages we need:

4a) Set up the problem

This time I will just call all the relevant packages at the beginning:

<<xAct`xTensor`
<<xAct`xPert`
<<xAct`xCoba`

We need to redefine the manifold etc. This time we will call the metric g and the manifold M.

DefManifold[M,4,{a,b,c,d,e,f,h,i,j,k,l}];
DefMetric[-1,g[-a,-b],CD,SymbolOfCovD->{"|","∇"},PrintAs->"g"];
$Pre=ScreenDollarIndices;
SetOptions[ContractMetric,AllowUpperDerivatives->True];
Unprotect[IndexForm];
IndexForm[LI[x_]]:=ColorString[ToString[x],RGBColor[0,0,1]];
Protect[IndexForm];
DefMetricPerturbation[g,metpert,ϵ];
PrintAs[metpert]^="h";

Now we set up our actions and find the equations of motion. We will do this for two theories: general relativity, with action

DefConstantSymbol[Mp,PrintAs->"(Subscript[M, pl])"]

ActionGR=Sqrt[-Detg[]]*RicciScalarCD[]*Mp^2/2

And the John term studied by 1603.04876: S = ∫ (d^4)x (Sqrt[-g])((((Subscript[M, pl])^2)/2) R + (η/2) (G^(μν)) (∂_μ)ϕ (∂_ν)ϕ)

DefTensor[Phi[],M,PrintAs->"ϕ"]

DefConstantSymbol[et,PrintAs->"η"]

ActionJohn=Sqrt[-Detg[]]*(Mp^2/2*RicciScalarCD[]+ et/2*EinsteinCD[b,c]CD[-b][Phi[]]CD[-c][Phi[]])

We can find the equations of motion using the techniques in section 2. First we define the tensor perturbation:

DefTensorPerturbation[PertPhi[LI[order]],Phi[],M,PrintAs->"δϕ"]

Then define the perturbed actions

pertGR=ToCanonical@ContractMetric@ExpandPerturbation@Perturbation@ActionGR
pertJohn=ToCanonical@ContractMetric@ExpandPerturbation@Perturbation@ActionJohn

And find the equations of motion using a variational derivative. We have a tensor equation for GR:

GRtenseqn=-VarD[metpert[LI[1],b,c],CD][pertGR]/Sqrt[-Detg[]]/.delta[-LI[1],LI[1]]->1//SeparateMetric[g]//RicciToEinstein//Expand//ContractMetric//ToCanonical

And both a scalar and tensor equation for the John theory:

Johntenseqn=-VarD[metpert[LI[1],b,c],CD][pertJohn]/Sqrt[-Detg[]]/.delta[-LI[1],LI[1]]->1//SeparateMetric[g]//RicciToEinstein//Expand//ContractMetric//SortCovDs//ToCanonical

Johnscalareqn=-VarD[PertPhi[LI[1]],CD][pertJohn]/Sqrt[-Detg[]]/.delta[-LI[1],LI[1]]->1//SeparateMetric[g]//RicciToEinstein//Expand//ContractMetric//ToCanonical

Now let’s use xCoba to set up the coordinate system we want. Let’s use a Schwarzchild-like system. We first define the chart:

DefChart[sch,M,{0,1,2,3},{t[],r[],θ[],ϕ[]}]

Then two scalar functions:

DefScalarFunction[ν]
DefScalarFunction[λ]

Then the metric

MatrixForm[metricarray={{-Exp[ν[r[]]],0,0,0},{0,Exp[λ[r[]]],0,0},{0,0,r[]^2,0},{0,0,0,r[]^2 Sin[θ[]]^2}}]

Then compure everything we need:

MetricInBasis[g,-sch,metricarray]

MetricCompute[g,sch,All,CVSimplify->Simplify]

Before moving on to put our equations of motion into this coordinate system we first define the energy-momentum tensor. We begin by defining both (T^(μν)) and (u^μ):

DefTensor[u[b],M,PrintAs->"u"]
DefTensor[T[b,c],M,PrintAs->"T"]

In this coordinate system, we have (u^μ) = ((e^((-ν(r))/2)),0,0,0). We can define vectors in this basis using the CTensor command:

Uup=CTensor[{Exp[-ν[r[]]/2],0,0,0},{sch}]

This tells xCoba that Uup is a tensor in the sch basis with indices up given by the vector I specified. We can also define the down tensor (Subscript[u, μ]):

Udown=CTensor[{Exp[ν[r[]]/2],0,0,0},{-sch}]

Now we need to tell xCoba that (u^μ)  and (Subscript[u, μ])  are equal to these in this coordinate system. Note that we need to do this for both the upper and lower index and we need to use an underscore the tell xCoba that b is a dummy index:

ComponentValue[u[b_],Uup[b]]
ComponentValue[u[-b_],Udown[-b]]

We can check this by raising the lower index or lowering the raised index:

g[b,c]u[-c]//ToBasis[sch]//ToBasis[sch]//TraceBasisDummy//ComponentArray//ToValues//ToValues

g[-b,-c]u[c]//ToBasis[sch]//ToBasis[sch]//TraceBasisDummy//ComponentArray//ToValues//ToValues

And we can check the normalisation

g[-b,-c]u[b]u[c]//ToBasis[sch]//ToBasis[sch]//TraceBasisDummy//ToValues//ToValues

u[-c]u[c]//ToBasis[sch]//ToBasis[sch]//TraceBasisDummy//ToValues//ToValues

As we can see, (u^μ) is a time-like vector, as we expect. Let’s define the energy-momentum tensor now. First, we define the pressure and density:

DefScalarFunction[den,PrintAs->"ρ"]
DefScalarFunction[pres,PrintAs->"P"]

and use the usual definition:

IndexSet[T[b_,c_],(den[r[]]+pres[r[]])u[b]u[c]+pres[r[]]g[b,c]]

Let’s check that this works

MatrixForm[EMup=T[b,c]//ToBasis[sch]//ToBasis[sch]//TraceBasisDummy//ComponentArray//ToValues//ToValues]

MatrixForm[Emmixed=T[b,-c]//ToBasis[sch]//ToBasis[sch]//TraceBasisDummy//ComponentArray//ToValues//ToValues]

MatrixForm[EMdown=T[-b,-c]//ToBasis[sch]//ToBasis[sch]//TraceBasisDummy//ComponentArray//ToValues//ToValues]

This is what we expect. Next, we can get some theory-independent equations from (Subscript[∇, μ]) (T^(μν))=0:

emcons=(CD[-b][T[b,c]]//ToBasis[sch]//ToBasis[sch]//TraceBasisDummy//ComponentArray//ToValues//ToValues)

We can solve this for P’

emcons[[2]]

emsol=Solve[emcons[[2]]==0,pres'[r[]]][[1]]

This gives us P’(ν’).

This is as far as we can go without a theory so let’s start by looking at GR. We’ll do John later.

4b) GR

We begin by putting the Einstein equation into our coordinates (remember the indices are down):

GRtenseqnsch=GRtenseqn//ToBasis[sch]//ToBasis[sch]//TraceBasisDummy//ComponentArray//ToValues//ToValues//Simplify

This gives us a matrix. We can extract the relevant components using Mathematica’s usual functions e.g. the ϕϕ-component is

GRtenseqnsch[[4]][[4]]

Let’s start with the 00-equation, (Subscript[G, 00]) = 8(Subscript[πGT, 00]) (note that the requisite factor of G is already contanined in the factor of (Subscript[M, pl])). At this point, I will define a constant G and a rule for turning (Subscript[M, pl]) into G. Note that we need to account for the definition of (Subscript[T, μν])= -(2/(Sqrt[-g])) ((δS)/((δg)^(μν))) in the equation.

DefConstantSymbol[G]

Grule={Mp^2->1/(8*Pi*G)}

GReqn00=GRtenseqnsch[[1]][[1]]==EMdown[[1]][[1]]/2/.Grule

This isn’t as nice as it could be because Mathematica doesn’t know anything about the metric potentials. Let’s change that:

$Assumptions=ν[x_]>0&&λ[x_]>0&&r[]>0&&G>0

GReqn001=GReqn00//Simplify

That’s better but Mathematica doesn’t know what form we want this in so it is best to manipulate it by hand:

GReqn002=-Exp[-λ[r[]]]+GReqn001*Exp[-λ[r[]]]//Expand

We can now define the mass:

DefScalarFunction[mass,PrintAs->"M"]

We can integrate the 00-equation once:

intGReqn00=2G*mass[r[]]==r[]-r[]*Exp[-λ[r[]]]

and check that this indeed solves the equation

testeqnGR=D[intGReqn00,r[]]-GReqn002

The only thing we are missing is the fact that M'=4((πr)^2)ρ:

mrule={mass'[r[]]->4*Pi*r[]^2*den[r[]]}

testeqnGR/.mrule

We can then solve the integrated equation for λ (remove the rule at the end to see why it is necessary):

λsolGR=Solve[intGReqn00,λ[r[]]][[1]]/.C[1]->0

Next, we have the rr-equation

GReqnrr=GRtenseqnsch[[2]][[2]]==EMdown[[2]][[2]]/2/.Grule

We can use our solution for λ to simplify this

GReqnrr2=GReqnrr/.λsolGR//Expand

Finally, we can solve this for ν’:

νsolGR=Solve[GReqnrr2,ν'[r[]]][[1]]

Recall that we also have P’(ν’):

emsol

So we can just subsistitue this in:

emsol/.νsolGR//Simplify

This is the Tolman-Oppenheimer-Volkov equation.

4c) John

Now let’s do something more fun, the John equation of motion. I will start with the scalar equation of motion so we need to tell xCoba that ϕ depends on r only. We first define a new scalar function and a new constant q

DefScalarFunction[φ,PrintAs->"ϕ"]
DefConstantSymbol[q]

Next, we need to set the scalar to be qt + φ(r)

ComponentValue[Phi[],q*t[]+φ[r[]]]

Johnscalareqnsch=Johnscalareqn//SeparateMetric[g]//ToBasis[sch]//ToBasis[sch]//TraceBasisDummy//ToValues//ToValues

Now we know we can integrate this once but this can be a pain so let me show you a trick for shift symmetric scalar-tensor theories. Note that the equation of motion is

(Subscript[∇, μ]) ((∂L)/(∂((Subscript[∇, μ])ϕ)))- (Subscript[∇, μ])((Subscript[∇, ν]) ((∂L)/(∂((Subscript[∇, μ])(Subscript[∇, ν])ϕ)))) = (Subscript[∇, μ]) (J^μ) = 0  so that (J^μ) = ((∂L)/(∂((Subscript[∇, μ])ϕ)))-(Subscript[∇, ν])((∂L)/(∂((Subscript[∇, μ])(Subscript[∇, ν])ϕ))) but this is nothing more than the equation of motion for a vector (Subscript[Z, μ]) = (Subscript[∇, μ])ϕ. To see how this works, let’s take the John action and replace  (Subscript[∇, μ])ϕ with (Subscript[Z, μ]):

DefTensor[Z[-b],M]

DefTensorPerturbation[PertZ[LI[order],-b],Z[-b],M,PrintAs->"δZ"]

ActionJohn2=ActionJohn/.CD[-c_][Phi[]]->Z[-c]/.CD[c_][Phi[]]->Z[c]

Now we find the EOM:

pertJohn2=ToCanonical@ContractMetric@ExpandPerturbation@Perturbation@ActionJohn2

Johnscalareqn2=VarD[PertZ[LI[1],-b],CD][pertJohn2]/Sqrt[-Detg[]]/.delta[-LI[1],LI[1]]->1//SeparateMetric[g]//RicciToEinstein//Expand//ContractMetric//ToCanonical

This looks a lot nicer than the other equation. In this case it was obvious from the action but it is not always so simple for more complicated theories. Let’s put this equation our coordinate system by first putting the factors of (Subscript[∇, μ])ϕ back:

J=Johnscalareqn2/.Z[b_]->CD[b][Phi[]]//SeparateMetric[g]//ToBasis[sch]//ToBasis[sch]//TraceBasisDummy//ComponentArray//ToValues//ToValues

If this really is (J^μ)then the 0-component does not contribute to the EOM because it is time-independent and so we have (J^r) = 0:

Johnscalareqnsch2=J[[2]]==0

This has two possible soltions: the trivial one, where ϕ’=0 and the non-trivial one where ϕ’ is finite but not zero. Let’s take the second on:

SS1=Solve[Johnscalareqnsch2,ν'[r[]]][[1]]

This is our first stellar structure equation. Next we need the tensor equations.

The tensor equation in this coordinate system is (note this takes a few minutes so now is a good time to ask questions):

Johntenseqnsch=Johntenseqn//SeparateMetric[g]//ToBasis[sch]//ToBasis[sch]//TraceBasisDummy//ComponentArray//ToValues//ToValues//Simplify

The 00- and rr- equations are

Johneqn00=Johntenseqnsch[[1]][[1]]==EMdown[[1]][[1]]/2/.SS1//Simplify

Johneqnrr=Johntenseqnsch[[2]][[2]]==EMdown[[2]][[2]]/2/.SS1//Simplify

Note that I have applied the equation for ν’ to these to simplify them. We already have an equation for ν’ and so we need one for ϕ’ and one for λ’:

SS2=Solve[{Johneqn00&&Johneqnrr}/.φ'[r[]]^2->U,{U,λ'[r[]]}][[1]]/.U->φ'[r[]]^2//Simplify

Finally, we need to solve for ϕ’’, which we do as follows

SS3=Map[1/2D[#,r[]]&,{SS2[[1]]},{2}]//Expand//Simplify

SS4=Join[SS1,SS2,SS3]

The full set of stellar structure equations are then:

MatrixForm[SSeqns={ν'[r[]]==(ν'[r[]]/.SS4),λ'[r[]]==(λ'[r[]]/.SS4/.SS4),φ'[r[]]^2==(φ'[r[]]^2/.SS4)}]//Simplify

In terms of the variables used in 1603.04876 we have A = Exp[ν] and B = Exp[-λ]:

DefScalarFunction[A]
DefScalarFunction[B]

ABrules={ν->Function[r,Log[A[r]]],λ->Function[r,-Log[B[r]]]}

SSeqnsAB=SSeqns/.ABrules

This can be tidied up by solving for A’ and B’

MatrixForm[Solve[SSeqnsAB/.φ'[r[]]^2->U,{A'[r[]],B'[r[]],U}]/.U->φ'[r[]]^2//FullSimplify//Transpose]

One can check by eye that these match. Finally, the full TOV system is composed of this set of equations and the one relating the pressure gradient to ν’ (or A’) that comes from energy-momentum conservation:

emsol

emsol/.ABrules//Simplify

5) Perturbation Theory

Now for the fun part - perturbation theory about an arbitrary space time. This requires you to load a new package called xPand, which you have to download separately. This contain s all of the packages we loaded previously so it is always best to quit the Kernel before you load it.

Quit

<<xAct`xPand`

Next, we need to define a background manifold, which is done in the same was as defining a regular manifold. We’ll define the metric here while we are at it.

DefManifold[M,4,{b,c,d,e,f,i,j}];
DefMetric[-1,g[-b,-c],CD,{";","\!\(\*OverscriptBox[\(∇\), \(_\)]\)"},PrintAs->"g"];

Next, we need to do a 3+1 decomposition. We do this using the SetSlicing Command. This tells xPand to decomose the metric (g) into a normal vector (n) and a metric on 3-dimensional hypersurfaces (h). It also defines a covariant derivative on these hypersurfaces cd. At this point we need to choose which background metric we want to perturb around. The possible options are: Anisotropic, BianchiB, BianchiA, BianchiI, FLCurved, FLFlat and Minkowski. Let’s stick to what we know best: Flat FRW.

SetSlicing[g,n,h,cd,{"|","\!\(\*OverscriptBox[\(D\), \(_\)]\)"},"FLFlat"];

Next, we need to define the matter perturbations and tell xPand which symbols to use for metric and matter perturbations:

DefMetricFields[g,dg,h];
DefMatterFields[u,du,h];

Now we are almost ready to calculate. Before that however, there are two useful pieces of book keeping. The first is to redifine the colours to make them easier to read. The second has to do with ease of use. The fundamental function is ToExpand but this requires 7 arguments, four of which are always the same and so it is useful to define a MyToExpand function which makes things quicker. The ord function collects things order by order in the perturbation parameter.

MakeBoxes[ϵ,StandardForm]:=StyleBox[ToString[$PerturbationParameter],FontColor->RGBColor[0.0,0.0,1]]
MyToxPand[expr_,gauge_,order_]:=ToxPand[expr,dg,u,du,h,gauge,order]
org[expr_]:=NoScalar@Collect[ContractMetric[expr],$PerturbationParameter,ToCanonical]

For ease, I will only look at scalar perturbations in this tutorial so I will turn off vectors and tensors now. I will also use coordinate time for now.

$FirstOrderVectorPerturbations=False;
$FirstOrderTensorPerturbations=False;
$ConformalTime=False;

Now we can begin. The MyToExpand[exp,gauge,order] function takes an expression exp and expands it in the gauge to a specified order. As an example, let’s expand the metric to 1st order:

metricpert=MyToxPand[g[-b,-c],"AnyGauge",1]//org

Okay, this is a bit messy so what’s happened. Well you see there are two parts: the zeroth order part (not containing ϵ’s) contains the background metric and the part proportional to ϵ is the first order perturbation. xPand perturbs along the normal n and in the orthogonal hypersurfaces described by h[-a,-b]. What about the letters E, B, ϕ and ψ? These are the metric perturbations. Usually we call ϕ Φ and we can rename it as follows:

PrintAs[ϕh]^="Φ";

xPand has a useful function for visualising tensors called VisualizeTensor:

VisualizeTensor[metricpert,h]

The 0-0 component is the n-n direction, the 0i the n-h and the ij the h-h. This is starting to look more familiar. Of course we know that things are less messy once one fixes a gauge but this is useful if one needs a gauge invariant formalism. The possible choices of gauge are

?$ListOfGauges

Let’s work in the Conformal Newtonian Gauge and switch back to conformal time.

$ConformalTime=True;

The numbers above the metric perturbations indicates their order. For example, 1 is the first order perturbation and 2 is the second order. In xPand, these are written in the form ϕh[order,time-derivative], for example

ϕh[LI[1],LI[2]]

is the second time-derivative of the first-order perturbation of ϕ. As an example, we can try

MyToxPand[RicciScalarCD[],"NewtonGauge",2]

which is the Ricci scalar expanded to second order in perturbation theory. Okay, now that we understand this, let’s try to derive something we know. Let’s look at the evolution of matter perturbations. Let’s start with the left hand side of the Einstein equations:

einsteinpert=MyToxPand[EinsteinCD[-b,-c],"NewtonGauge",1]

This gives us the background equations and the perturbations. Let’s begin by extracting the Friedmann equations. We can do this using the ExtractOrder command, which extracts a specific order in pertrubation theory and the ExtractComponents command, which pulls out components of tensor equations.

FriedmannLHS1=ExtractComponents[ExtractOrder[einsteinpert,0],h,{"Time","Time"}]

FriedmannLHS2=h[b,c]ExtractComponents[ExtractOrder[einsteinpert,0],h,{"Space","Space"}]//ContractMetric

Note that in the second Friedmann equation I have taken the space-space trace of the ij equation. Next, we need to set up the energy-momentum tensor. We do this as follows (note I will use only dust in what follows but I have left the option of putting pressure pertrubations in).

DefTensor[T[b,c],M]
$Dust=True;
IndexSet[T[b_,c_],((ρ[u][]+If[$Dust,0,P[u][]])u[b]u[c]+(If[$Dust,0,P[u][]])g[b,c])]

We need to perturb this to first order:

empert=MyToxPand[T[-b,-c],"NewtonGauge",1]
VisualizeTensor[empert,h]

We also need the conservation equation (Subscript[∇, μ])(Subscript[(T^μ), ν])=0.

conspert=MyToxPand[CD[-b][T[b,-c]],"NewtonGauge",1]

DefConstantSymbol[G]

Friedmann1=FriedmannLHS1==8*Pi*G*ExtractComponents[ExtractOrder[empert,0],h,{"Time","Time"}]

Friedmann2=FriedmannLHS2==h[b,c]8*Pi*G*ExtractComponents[ExtractOrder[empert,0],h,{"Space","Space"}]

We also have the contunuity equation

Cont=ExtractComponents[ExtractOrder[conspert,0],h,{"Time"}]==0

We can solve these for quantities that may simplify the equations for the perturbations:

backgroundsol1=Last[Solve[Friedmann1,ρu[]]]
backgroundsol2=Last[Solve[Friedmann2,Hh[LI[0],LI[1]]]]
contsol=Last[Solve[Cont,ρu[LI[0],LI[1]]]]

These should be familiar. Now for the perturbations. Our ultimate goal is to reproduce the equation for the evolution of CDM perturbations inside the horizon during the matter era: (Overscript[δ, ..])+H(OverDot[δ])-4πGa^2ρ=0 . This means we need to find the perturbed Einstein and continuity equations, simplify them using the background equations of motion and then eliminate all of the variables in terms of the density perturbation. We will also need to fourier transform the equations which can be done using the following rule:

DefConstantSymbol[k,PrintAs->"k"]

FourierRules={IndexRule[cd[-b_][cd[b_][x_]],-(k^2) x],IndexRule[cd[b_][cd[-b_][x_]],-(k^2) x]}

Let’s begin with the perturbed continuity equation. This has both a 0 and an i component:

pertcont0=ExtractComponents[ExtractOrder[conspert,1],h,{"Time"}]==0

pertconti=ExtractComponents[ExtractOrder[conspert,1],h,{"Space"}]==0

These have spatial derivatives that need to be Fourier transformed:

pertcont0new=pertcont0/.FourierRules

pertcontinew=cd[c][pertconti[[1]]]==0/.FourierRules/.contsol

One can use the first equation to eliminate the velocity perturbation:

velsol=Last[Solve[pertcont0new,Vshu[LI[1],LI[0]]]]

We need to take a derivative of this in order to eliminate it in the other equation. This is done as follows:

velsol2=Map[LieD[n[b]],velsol,{2}]/.contsol//Simplify

Now we can put this into the other equation:

rhoev=pertcontinew2=pertcontinew/.velsol/.velsol2//Simplify

This is the evolution equation for the density perturbation but it contains second derivatives of ψ, which must be eliminated. This can be done by looking at the Einstein equations. We start with the 00 component.

g00=ExtractComponents[ExtractOrder[einsteinpert,1],h,{"Time","Time"}]==8*Pi*G*ExtractComponents[ExtractOrder[empert,1],h,{"Time","Time"}]/.FourierRules//Simplify

And in a similar way for ij

g0i=ExtractComponents[ExtractOrder[einsteinpert,1],h,{"Time","Space"}]==8*Pi*G*ExtractComponents[ExtractOrder[empert,1],h,{"Time","Space"}]/.FourierRules//Simplify

This has derivatives so let’s take one more derivative and fourier transform. We can also eliminate the velocity:

g0i2=cd[c][g0i[[1]]]==0/.FourierRules/.velsol//Simplify

People tend to use the symmetric trace free part of the Einstein equations rather than the ij components directly and xPand has a method for finding this:

gij=ExtractComponents[ExtractOrder[einsteinpert,1],h,{"Space","Space"}]==8*Pi*G*ExtractComponents[ExtractOrder[empert,1],h,{"Space","Space"}]/.FourierRules//Simplify

stf=STFPart[gij,h]//ContractMetric//ToCanonical

We know that this is solved by Φ=ψ:

stf/.ψh[LI[1],LI[0]]->ϕh[LI[1],LI[0]]

and so we can make the rule

grrule={ψh[LI[1],LI[0]]->ϕh[LI[1],LI[0]],ψh[LI[1],LI[1]]->ϕh[LI[1],LI[1]],ψh[LI[1],LI[2]]->ϕh[LI[1],LI[2]]}

Now we can use this in the evolution equation for ρ along with the background solutions:

rhoev2=rhoev/.grrule/.backgroundsol2//Simplify

We are left with a set of equations for ϕ and ρ. In general, these are coupled but we are interested in sub-horizon perturbations. Looking at the 00 equation, we have:

g002=g00/.grrule/.backgroundsol1//Simplify

One can see that the sub-horizon limit can be achieved by sending H->0:

g003=g002/.Hh[]->0

Solving we find

phisol=Last[Solve[g003,ϕh[LI[1],LI[0]]]]

Looking at the evolution equation, (ϕ^('' ))can be neglected, as can ϕ'since it is supressed relative to the k^2 term:

rhoev3=rhoev2/.ϕh[LI[1],LI[2]]->0/.ϕh[LI[1],LI[1]]->0

Substituting in for ϕ we find

rhoev4=rhoev3/.phisol

The final thing to do is to switch from ρ to δ=δρ/ρ. We do this by first finding relations between the derivatives:

Δ=ρu[LI[1],LI[0]]/ρu[LI[0],LI[0]]

Δ'=Map[LieD[n[b]],Δ,{0}]/.contsol//Simplify

Δ''=Map[LieD[n[b]],Δ',{0}]/.contsol/.backgroundsol2//Simplify

Next, we invert these relations:

δ1=Last[Solve[Δ==δ,ρu[LI[1],LI[0]]]]
δ2=Last[Solve[Δ'==δ',ρu[LI[1],LI[1]]]]/.δ1
δ3=Last[Solve[Δ''==δ'',ρu[LI[1],LI[2]]]]/.δ1/.δ2/.δ1

Finally, we subsititute these in to find

rhoev4/ρu[LI[0],LI[0]]/.δ3/.δ2/.δ1//Simplify

And there you have it, as promised. Of course, this was a simple thing to do but you can do more complicated things like derive the evolution of CMB photons and the evolution of matter perturbatons outside the horizon, this just means you need to subsititute more equations into each other until they close.

As a final example, and to show you how scalars work, let’s perturb the covariant quartic galileon. We quit the kernel earlier so let’s just re-run the commands to generate the equations of motion. Note that xPand comes with its own tensor φ. You need to use this since it already knows how to slice covariant derivatives etc. If you want more scalars you can look up the DefProjectedTensor command. Now for the equations:

DefConstantSymbol[L3,PrintAs->"(Subscript[Λ, 3])"]
DefConstantSymbol[L4,PrintAs->"(Subscript[Λ, 4])"]
DefConstantSymbol[Mp,PrintAs->"(Subscript[M, p])"]
DefTensor[X[],M,PrintAs->"X"];
BianchiID={CD[-d_]@φ[]*CD[-b_]@RicciCD[b_,d_]:>1/2 CD[-d]@φ[] g[b,d] CD[-b]@RicciScalarCD[],CD[d_][φ[]]*CD[-b_][RicciCD[-d_,b_]]:>1/2 CD[d]@φ[]*CD[-d]@RicciScalarCD[]};
IndexSet[X[],-1/2CD[-b][φ[]]CD[b][φ[]]];
LGal=Sqrt[-Detg[]]*(X[]+1/(2*L3^3)*X[]*CD[-b][CD[b][φ[]]]-(2X[]/L4^6)*(X[]*RicciScalarCD[]+2*(CD[-d][CD[d][φ[]]]CD[-b][CD[b][φ[]]]-CD[-d][CD[-b][φ[]]]CD[d][CD[b][φ[]]]))+Mp^2/2*RicciScalarCD[])//ToCanonical;
LGalpert=ToCanonical@ContractMetric@ExpandPerturbation@Perturbation@LGal;
tensoreom=2(-VarD[dg[LI[1],b,c],CD][LGalpert]/Sqrt[-Detg[]]/.delta[-LI[x_],LI[y_]]->KroneckerDelta[x,y]/.BianchiID//SortCovDs//SeparateMetric[g]//RicciToEinstein//Expand//ContractMetric//ToCanonical)//ToCanonical//SortCovDs//ContractMetric//ToCanonical
scalareom=((VarD[φ[LI[1]],CD][LGalpert]/Sqrt[-Detg[]]/.delta[-LI[x_],LI[y_]]->KroneckerDelta[x,y]//SortCovDs//ToCanonical)/.BianchiID//ContractMetric//ToCanonical)//ContractMetric//ToCanonical

These are both absolutely horrific but xPand can deal with it. Let’s switch back to coordinate time.

$ConformalTime=False;

And perturb:

tenspert=MyToxPand[tensoreom,"NewtonGauge",1]

115 seconds is about 2 minutes - certainly faster than I could do this. From this we get the Friedmann equation:

FriedLHS=ExtractComponents[ExtractOrder[tenspert,0],h,{"Time","Time"}]

And the Einstein equations, for example, the STF part:

STFPart[ExtractComponents[ExtractOrder[tenspert,1],h,{"Space","Space"}],h]//ContractMetric//ToCanonical

Similarly, we can expand the scalar eom:

scalarpert=MyToxPand[scalareom,"NewtonGauge",1]

43 seconds, not bad. Here we get the scalar eom in the background

ExtractOrder[scalarpert,0]==0

And the perturbed EOM

ExtractOrder[scalarpert,1]==0

From here, one can repeat the process above and find the new equation for the evolution of density perturbations and more. The possibilities are limitless.

Epilogue: What’s Next?

We’ve done all of our computations, so what’s next? Well, we write up of course. It turns out you don’t need to copy and paste and put in the latex yourself, xAct can do it for you! You need to load the TexAct package:

<<xAct`TexAct`

Let’s tex up the Klein-Gordon equation

DefScalarFunction[V]
DefTensor[Phi[],M,PrintAs->"ϕ"]

KG=CD[-b][CD[b][Phi[]]]+V'[Phi[]]==0

Tex[KG]

AddEquationMarks[%,"equation"]

You can even add the equation environment!