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ₛ::Float64endstruct SchwarzschildAtlas <: AbstractAtlas{ℝ} endManifolds.manifold_dimension(::BlackHoleOutside) = 4Manifolds.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{ℝ} endfunction 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 xendfunction 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 pendfunction 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) = falseManifolds.get_chart_index(::BlackHoleOutside, ::SchwarzschildAtlas, p) = nothingManifolds.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 = nothingX0 =  [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.0sol = 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.0x_max = 12.0y_min = -12.0y_max = 12.0# sampling Kretschmann scalarks_samples = 300ks_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    endendks_vals = [log.(get_ks(a_x, a_y)) for a_x in ks_x, a_y in ks_y]# plottingx_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)endlines!(ax, x_values, y_values, color=:white, label="trajectory")axislegend(ax)fig

Below is a pre-recorded animation: