Working with a metric defined in a chart

In this tutorial we will learn how to use charts to solve a basic problems when the metric is defined using inner product.

We start by loading the necessary libraries.

using Manifolds, RecursiveArrayTools, OrdinaryDiffEq, DiffEqCallbacks, CairoMakie

We will work with the exterior of a Schwarzschild black hole with Schwarzschild radius $r_s$. It is a 4-dimensional manifold represented using vectors of length 4.

struct BlackHoleOutside <: AbstractManifold{ℝ}
    rₛ::Float64
end

struct SchwarzschildAtlas <: AbstractAtlas{ℝ} end

Manifolds.manifold_dimension(::BlackHoleOutside) = 4
Manifolds.representation_size(::BlackHoleOutside) = (4,)

To work with this manifold, we need to define a metric. However, we first need to introduce an atlas. For this example let’s use an atlas with one chart defined by Schwarzschild coordinates. We also indicate that we will use the Levi-Civita affine connection.

struct SchwarzschildAtlas <: AbstractAtlas{ℝ} end

function Manifolds.get_parameters!(M::BlackHoleOutside, x, ::SchwarzschildAtlas, i, p)
    x[1] = p[1] # t
    r = norm(p[2:4])
    x[2] = r
    x[3] = acos(p[4] / r) # θ
    x[4] = atan(p[3], p[2]) # ϕ
    return x
end

function Manifolds.get_point!(M::BlackHoleOutside, p, ::SchwarzschildAtlas, i, x)
    p[1] = x[1]
    p[2] = x[2] * sin(x[3]) * cos(x[4])
    p[3] = x[2] * sin(x[3]) * sin(x[4])
    p[4] = x[2] * cos(x[3])
    return p
end

function Manifolds.affine_connection!(M::BlackHoleOutside, Zc, A::SchwarzschildAtlas, i, a, Xc, Yc)
    return Manifolds.levi_civita_affine_connection!(M, Zc, A, i, a, Xc, Yc)
end

Since we use only one chart, we can switch off chart switching functionality:

Manifolds.check_chart_switch(::BlackHoleOutside, A::SchwarzschildAtlas, i, a) = false
Manifolds.get_chart_index(::BlackHoleOutside, ::SchwarzschildAtlas, p) = nothing
Manifolds.get_chart_index(::BlackHoleOutside, ::SchwarzschildAtlas, i, a) = nothing

The metric implemented below is the Schwarzschild metric (signature +, −, −, −, and geometric units $c = G = 1$). In Schwarzschild coordinates at $a=(t, r, \theta, \varphi)$ the inner product of two tangent vectors at point with parameters $a$ with coordinates $X_c$, $Y_c$ is

\[g_a(X_c, Y_c) = X_c^{\top}\operatorname{diag}\!\left(1-\frac{r_s}{r},\; -\frac{1}{1-r_s/r},\; -r^2,\; -r^2\sin^2\theta\right) Y_c.\]

where $r_s$ is the Schwarzschild radius and the coordinates of two tangent vectors are given in the basis ∂/∂t, ∂/∂r, ∂/∂θ, ∂/∂φ. The method inner in this tutorial computes the inner product of two tangent vectors expressed in the chart coordinate basis.

function Manifolds.inner(M::BlackHoleOutside, ::SchwarzschildAtlas, i, a, Xc, Yc)
    t, r, θ, ϕ = a
    r_block = (1 - M.rₛ / r)
    return Xc[1] * r_block * Yc[1] - Xc[2] * Yc[2] / r_block - r^2 * (Xc[3] * Yc[3] + (sin(θ)^2) * Xc[4] * Yc[4])
end

And this is enough to compute geodesics:

M = BlackHoleOutside(1.0)
p0 = [0.0, 10.0, 0.0, 0.0]

A = SchwarzschildAtlas()
i = nothing

X0 =  [1.0, 0.0, 0.18, 0.0]
a_p0 = get_parameters(M, A, i, p0)
B = induced_basis(M, A, a_p0)
c_X0 = get_coordinates(M, p0, X0, B)
final_time = 5000.0
sol = Manifolds.solve_chart_exp_ode(M, a_p0, c_X0, A, i; final_time=final_time, solver=Tsit5())

sampled_solution = sol(range(0.0, final_time; length=20000));

We can also calculate the Kretschmann scalar and display the trajectory. Value of the Kretschmann scalar corresponds to background color.

x_min = -12.0
x_max = 12.0
y_min = -12.0
y_max = 12.0

# sampling Kretschmann scalar
ks_samples = 300
ks_x = range(x_min, x_max; length=ks_samples)
ks_y = range(y_min, y_max; length=ks_samples)
function get_ks(a_x, a_y)
    if a_x^2 + a_y^2 > M.rₛ^2
        return Manifolds.kretschmann_scalar(M, A, i, get_parameters(M, A, i, [0.0, a_x, a_y, 0.0]))
    else
        # the curvature gets very high inside, so we skip computing it
        return 0.0
    end
end

ks_vals = [log.(get_ks(a_x, a_y)) for a_x in ks_x, a_y in ks_y]

# plotting
x_values = [s[1][2] for s in sampled_solution]
y_values = [s[1][3] for s in sampled_solution]

fig = Figure(; size=(800, 800))
ax = Axis(fig[1, 1]; title="2D Plot of Sampled Solution", xlabel="x", ylabel="y", aspect = AxisAspect(1))

xlims!(ax, x_min, x_max)
ylims!(ax, y_min, y_max)
hm = heatmap!(ax, ks_x, ks_y, ks_vals; colormap=:viridis)

Colorbar(fig[:, end+1], hm, label="log(Kretschmann scalar)")

begin
    # filled polygon representing the black hole
    θ = range(0, 2π, length=400)
    xs = cos.(θ) .* M.rₛ
    ys = sin.(θ) .* M.rₛ
    poly!(ax, xs, ys, color = :black)
end

lines!(ax, x_values, y_values, color=:white, label="trajectory")

axislegend(ax)
fig

Below is a pre-recorded animation: