GRQUICK source: https://library.wolfram.com/infocenter/MathSource/8329/
This version can be executed on Mathematica 13 and higher.

Basics

This notebook gives examples of the functions in GRQUICK.   First we load the package. 

To install the package go to:
1.) go to File->Install
2.) in the dialogue box that pops up select:
3.) Type of Item to Install: Package
3.) Source->From File
4.) find and select the file GRQUICK-1.wl

To load the package use the “<<” command below.

<<GRQUICK-1.wl

Here we load a list of functions and obtain information on the Ricci function.

?helpGRQUICK

?Ricci

The starting point of GRQUICK is to input a metric using “Metin”.  We first use the auto-loader and then we input our own spherically symmetric metric.

?Metin

The Schwarzschild metric can be loaded as

Metin["Schwar"]

A list of defualt metrics can be loaded using ?DefaultMet

?DefaultMet

Or we can use our own metric matrices

Metin[DiagonalMatrix[{a[t,x],b[t,x]}],{t,x}]

Note, dimension does not matter.  We shall focus on 4-d spherically symmetric space-times for now.

Metin[DiagonalMatrix[{m[t,r],n[t,r],(r^2),(r^2)((Sin[θ])^2)}],{t,r,θ,ϕ}]

Now that we have used “Metin” we can look at the metric tensor using GMET.

MatrixForm[GMET]

Similarly, “GINV” can be used to display all of the inverse metric.

GINV

After a metric has been input we can calculate relevant geometric objects.  For example we calculate (Subscript[(Γ^0), 00])

Christoffel[{{0},{0,0}}]

, and or (Subscript[(R^0), 101])  (notice the index structure of the tensor!!!).

Riemann[{{0},{1,0,1}}]

We can also obtain specific components of the metric and inverse metric using “Metric” and “Imetric”.

Metric[0,0]

IMetric[0,0]

We can also calculate the Ricci and Einstein tensors.

Simplify[Ricci[0,0]]

Einstein[0,0]

Rscalar gives the Ricci scalar.

Rscalar

If we wish we can display all of the Einstein tensor as a NxN matrix using GARRAY.  Similar arrays exist for the other tensors as well but they can only be called either if 
1) said tensor has been computed or 
2) the Einstein tensor has been computed (see ?helpGRQUICK).

Simplify[GARRAY]

GRQUICK can also perform basic operations on contravariant vectors. For example, the norm of the position vector:

Metnorm[{t,r,θ,ϕ}]

The scalar product of two contravariant vectors is:

Metprod[{a,b,c,d},{t,r,θ,ϕ}]

?Metprod

We can also normalize vectors:

Normalizevec[{t,r,θ,ϕ}]

If entering the covariant vector is easier then one can use RaiseVector to switch, for example

RaiseVector[{a,b,c,d}]

Or we can lower the vector

LowerVector[{a,b,c,d}]

Geodesics

Now we look at what GRQUICK can do with geodesics. We focus on the Schwarzschild metric for simplicity.

Metin[DiagonalMatrix[{-(1-((2 M)/r)),(((1-((2 M)/r)))^(-1)),(r^2),(r^2)((Sin[θ])^2)}],{t,r,θ,ϕ}]

A faster way of loading default metrics is by inputing a string

Metin["Schwar"]

We have just loaded the Schwarzschild metric instantly.  For a list of default metrics see ?DefaultMet.  (Future releases will allow saving of metrics).

We populate the arrays GARRAY, CHARRAY, RsARRAY, and RARRAY by using Einstein[0,0].

Einstein[0,0];

Now we can find the geodesic equation, namely (((d^2)(x^μ))/((dτ)^2))+(Subscript[(Γ^μ), αβ])(((d (x^α))/(dτ))((d (x^β))/(dτ))) by using the “Geodesics” function.

MatrixForm[Geodesics[τ]]

Similarly we can use a different parameter for the geodesic equations (such as an affine parameter).

MatrixForm[Geodesics[λ]]

If the metric does not have any non-numerical variables other than the coordinates, we can plot the geodesic solutions.  Let us set M=1 and the speed of light to 1:

Setlight[1];M=1

Here the Setlight function is used to set the speed of light to 1.  The default light speed is c=1 but it is necessary to properly define c for the units used for the metric.


If we know the initial conditions for the geodesic, we can plot it.  Let us take the initial position and 4-velocity to be.

init={0,3,0,0};

vinit={(1/(Sqrt[1-(2/3)])),0,0,0};

For the above initial conditions, we plot “t” and “r” with τ∈{0,3}. (NOTE: just like any differential equation, one needs to pay attention to the initial conditions and watch out for any coordinate singularities)

Geoplot[τ,init,vinit,{0,3},{t,r}]

NOTE: Geoplot requires a properly normalized vector to function.  For example:

Geoplot[τ,init,{(4/(Sqrt[1-(2/3)])),0,0,0},{0,3},{t,r}]

For time-like vectors we can input ("Max",((dx^1)/(dτ)),((dx^2)/(dτ)),((dx^3)/(dτ))) and have Geoplot properly normalize ((dx^0)/(dτ)) (NOTE:  the normalization must be accurate to 10^-10 precision). 

As an example we look at an object orbiting in the plane θ=(π/2).

init={0,4,(π/2),0};

vinit={"Max",0,0,1};

Geoplot[τ,init,vinit,{0,10},{t,r,ϕ}]

The normalization has been automatically done for us.  The text “Max” tells us to use the maximum solution, ((dx^0)/(dτ)), of
-(c^2)-(Subscript[g, ij])(((dx^i)/(dτ))((dx^j)/(dτ)) )=(((Subscript[g, 00])(((dx^0)/(dτ))))^2)+2(Subscript[g, 0i])(((dx^i)/(dτ))((dx^0)/(dτ))).
If we want the minimum solution one replaces “Max” with “Min”.

Once Geoplot has been used, an interpolation object called “GEOSOL” is created that contains the solutions to the geodesic equations.

GEOSOL

“GEOSOL” be used to manipulate the numerical geodesic solutions.  For example, for the previous initial conditions, we can plot the particle trajectory in the x-y plane:

ParametricPlot[Evaluate[{r[τ] Cos[ϕ[τ]],r[τ] Sin[ϕ[τ]]}/.GEOSOL],{τ,0,10}]

Clear[M]

It is also useful to have a function which checks whether a specific (x^μ)(u) is a solution.  Issol[{(Subscript[x, 0])(u),(Subscript[x, 1])(u),(Subscript[x, 2])(u),(Subscript[x, 3])(u)},u] returns either True, False, or an expression which must be equal to zero if the four vector is to be a geodesic. Let us see if the trajectory
 (x^μ)(τ)={τ, R, (π/2),0}  is a solution to the geodesic equation

Issol[{τ,R, (π/2),0},τ]

Thus (x^μ) is a solution if R=∞  (Note: R=2M should always be treated with care because the metric is singular in these coordinates. Any expression given by Issol that gives g((x^μ))=∞ should be checked clearly.)

?Issol

DISCLAIMER: Issol is by no means perfect. If Issol gives false it does not mean your geodesic is not a solution it might mean that the program has not done enough simplifications or taken into account algebraic identities.

Finally

?ReplaceCoords

This function is useful when wanting to input a particular value to your coordinates. For example

ReplaceCoords[GMET,{τ, R, (π/2),0}]

ReplaceCoords[Christoffel[{{1},{1,1}}],{τ, R, (π/2),0}]

Christoffel[{{1},{1,1}}]

Covariant Derivatives and Tensor Definitions

The index structure in GRQUICK is {{top indices},{bottom indices}}.  For a purely covariant tensor one can use {bottom indices}, and for a pure contravariant one must use {{top indices},{}} (this will be fixed in the future).

We switch to our original spherical symmetric metric.

Metin[DiagonalMatrix[{m[t,r],n[t,r],(r^2),(r^2)((Sin[θ])^2)}],{t,r,θ,ϕ}]

Einstein[0,0] Re-Populates arrays

Einstein[0,0];

We calculate (Subscript[∇, 0])(Subscript[G, 00])

Coderv["Einstein",{0,0},0]

Coderv also works with more complex tensor such as:

Coderv["Riemann",{{0},{1,0,1}},0]

We now do a quick check of Coderv.

check=(D[Riemann[{{0},{1,0,1}}],t]+
Sum[Christoffel[{{0},{0,ss}}] Riemann[{{ss},{1,0,1}}],{ss,0,3}]
-Sum[Christoffel[{{ss},{0,1}}]Riemann[{{0},{ss,0,1}}],{ss,0,3}]
-Sum[Christoffel[{{ss},{0,0}}]Riemann[{{0},{1,ss,1}}],{ss,0,3}]
-Sum[Christoffel[{{ss},{0,1}}]Riemann[{{0},{1,0,ss}}],{ss,0,3}]);

FullSimplify[Expand[check-Coderv["Riemann",{{0},{1,0,1}},0]]]

Along with the standard tensors we can define new tensors exactly the way one defines functions in Mathematica (I am working on a define tensor function but for now this works).  Here is an example of a tensor definition of (Subscript[(G^μ), ν]):

Gupdn[{{a_},{b_}}]:=Sum[IMetric[a,jj]Einstein[jj,b],{jj,0,3}]

Now we can take covariant derivatives of the new object.

Coderv["Gupdn",{{0},{0}},0]

We now define a rank 2 contravariant tensor and take a covariant derivative (the bad notation for pure contravariant tensors will be fixed later).:

Gupup[{{a_,b_},{}}]:=Sum[IMetric[b,ss]Gupdn[{{a},{ss}}],{ss,0,3}]

Coderv["Gupup",{{0,0},{}},0]

Finally for rank 2 covariant tensor we can take their trace as :

Tenstr["Einstein"]-Sum[Gupdn[{{aa},{aa}}],{aa,0,3}]

Vectors

Vectors have a special place in GRQUICK because one can either put them directly into “Coderv” or functionally define them.

For example let us take (Subscript[∇, 0]) (x^0):

x={t,r,θ,ϕ}

X[a_]:=x[[a+1]]

NOTE: the +1 comes from the fact that we want the index to start at zero but mathematica starts the index at 1.

We can now take a covariant derivative in two ways:

Coderv["X",{{0},{}},0]

Coderv[{t,r,θ,ϕ},{{0},{}},0]

Coderv["X",{{0},{}},0]-Coderv[{t,r,θ,ϕ},{{0},{}},0]

NOTE: a vector in “Coderv” is defined to be covariant or contravariant based solely on the input of the index structure. For example:

Coderv[{t,r,θ,ϕ},{0},0]

takes the input vector to be covariant

Clear[x,X]

Double Derivatives

Double (or higher) covariant derivatives require a bit more work to obtain but it’s not that bad. (Future work will automate higher covariant derivatives).

In order to show the abilities of GRQUICK we now switch to an arbitrary diagonal matrix (we don’t use a fully arbitrary metric because it takes a fair amount of time and memory).

Metin[DiagonalMatrix[{f00[t,x,y,z],f11[t,x,y,z],f22[t,x,y,z],f33[t,x,y,z]}],{t,x,y,z}]

MatrixForm[GMET]

For the new metric, we populate the relevant arrays

Einstein[0,0];

Now we define a vector and the covariant derivative of the vector: (Subscript[(R^s), a b v])

A={(((t+x+y+z))^1),(((t+x+y+z))^2),(((t+x+y+z))^3),(((t+x+y+z))^4)}

Dvec[{v_,p_}]:=Coderv[A,{v},p]

Now we prove (Subscript[∇, σ]) (Subscript[∇, p])(Subscript[A, v]) -(Subscript[∇, p]) (Subscript[∇, σ])(Subscript[A, v])= (Subscript[A, s])(Subscript[(R^s), v p σ])

FullSimplify[Table[(Coderv["Dvec",{v,p},σ]-Coderv["Dvec",{v,σ},p])
-Sum[A[[kk]] Riemann[{{kk-1},{v,p,σ}}],{kk,1,4}],{v,0,3},{p,0,3},{σ,0,3}]]

These are all of the functions GRQUICK has.  Future plans include a define tensor function, which makes defining a tensor easier (i.e. implicit dummy summation) as well as a covariant derivative function which can take multiple co-derivatives.