Curvature quantities on Riemannian manifolds

In this notebook, we discuss the definition and computation of different curvature quantities of Riemannian manifolds and their relations, see this notebook for a brief introduction to Riemannian manifolds. We focus on the Riemann curvature tensor, curvature operator, Ricci curvature tensor, scalar curvature, and (in 3D) the Einstein tensor.

We use two- and three-dimensional numerical examples to verify the identities. We consider the Poincare disc in 2D and the three-sphere in 3D. For the numerical computations, we interpolate the exact metric into the Regge finite elements.

[1]:

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *
import ngsdiffgeo as dg
from ngsolve.fem import Einsum, LeviCivitaSymbol

mesh2d = Mesh(
    OCCGeometry(Circle((0, 0), 0.7).Face(), dim=2).GenerateMesh(maxh=0.05)
).Curve(3)
pcd = dg.PoincareDisk()
Draw(Norm(pcd.metric), mesh2d)

eps = 2e-1
mesh3d = Mesh(
    OCCGeometry(Box((eps, eps, eps), (pi - eps, pi - eps, 2 * pi - eps))).GenerateMesh(
        maxh=0.3
    )
)
sphere3 = dg.Sphere3()
Draw(Norm(sphere3.metric), mesh3d, order=3)


# Euclidean Levi-Civita permuting symbols in 2D and 3D
Eps2d = LeviCivitaSymbol(2)
Eps3d = LeviCivitaSymbol(3)

# Interpolate exact metric to Regge space
g2d = GridFunction(HCurlCurl(mesh2d, order=5))
g3d = GridFunction(HCurlCurl(mesh3d, order=4))
with TaskManager():
    g2d.Set(pcd.metric, dual=True)
    g3d.Set(sphere3.metric, dual=True)

# Define Christoffel symbols of first and second kind
chr1_2d = g2d.Operator("christoffel")
chr2_2d = g2d.Operator("christoffel2")
chr1_3d = g3d.Operator("christoffel")
chr2_3d = g3d.Operator("christoffel2")

Riemann curvature tensor

We start with the fourth-order Riemann curvature tensor \(\mathfrak{R}\). Its coordinate expression is \begin{align*} \mathfrak{R}_{ijkl} = \partial_i \Gamma_{jkl}-\partial_j \Gamma_{ikl} + \Gamma_{ik}^p\Gamma_{jlp}- \Gamma_{jk}^p\Gamma_{ilp}. \end{align*}

First, we note that the first two terms can be rewritten using the incompatibility operator \(\mathrm{inc}(g)= \varepsilon^{ij}\varepsilon^{kl}\partial_i\partial_k g_{jl}\) in 2D and \(\mathrm{inc}(g)^{ij}= \varepsilon^{ikl}\varepsilon^{jmn}\partial_k\partial_m g_{ln}\) in 3D. There holds \begin{align*} \partial_i \Gamma_{jkl}-\partial_j \Gamma_{ikl} = \begin{cases}\frac{1}{2}\varepsilon_{ij}\varepsilon_{kl}\mathrm{inc}(g), & \text{in 2D},\\ \frac{1}{2}\varepsilon_{ija}\varepsilon_{klb}\mathrm{inc}(g)^{ab}, & \text{in 3D}. \end{cases} \end{align*}

We compare the Riemann differential operator from the Regge finite element with its coordinate expression from the metric \(g\) and the exact Riemann curvature tensor. Note that the numerical computations fit up to numerical rounding error precision. Increasing the resolution of the mesh reduces the relative error between the numerical and exact curvature. We will learn later how to approximate the curvature more accurately from the metric tensor.

[2]:

# NGSolve Riemann curvature tensor
Riem1_2d = g2d.Operator("Riemann")
# manually compute Riemann curvature tensor
Riem2_2d = (
    0.5 * g2d.Operator("inc") * Einsum("ij,kl->ijkl", Eps2d, Eps2d)
    + Einsum("ikp,jlp->ijkl", chr2_2d, chr1_2d)
    - Einsum("jkp,ilp->ijkl", chr2_2d, chr1_2d)
)
# exact Riemann curvature tensor
Riem_ex_2d = pcd.Riemann

# The same in 3D
Riem1_3d = g3d.Operator("Riemann")
Riem2_3d = (
    0.5 * Einsum("ija,klb,ab->ijkl", Eps3d, Eps3d, g3d.Operator("inc"))
    + Einsum("ikp,jlp->ijkl", chr2_3d, chr1_3d)
    - Einsum("jkp,ilp->ijkl", chr2_3d, chr1_3d)
)
Riem_ex_3d = sphere3.Riemann

# check that the expressions are equal
with TaskManager():
    print(
        f"err Riemann curvature tensor 2d = {sqrt(Integrate(InnerProduct(Riem1_2d - Riem2_2d,Riem1_2d - Riem2_2d),mesh2d)):.8f}"
    )
    print(
        f"rel err Riemann curvature tensor 2d exact = {sqrt(Integrate(InnerProduct(Riem1_2d - Riem_ex_2d,Riem1_2d - Riem_ex_2d),mesh2d))/sqrt(Integrate(InnerProduct(Riem_ex_2d, Riem_ex_2d),mesh2d)):.8f}"
    )
    print(
        f"err Riemann curvature tensor 3d = {sqrt(Integrate(InnerProduct(Riem1_3d - Riem2_3d,Riem1_3d - Riem2_3d),mesh3d)):.8f}"
    )

    print(
        f"rel err Riemann curvature tensor 3d exact = {sqrt(Integrate(InnerProduct(Riem1_3d - Riem_ex_3d,Riem1_3d - Riem_ex_3d),mesh3d))/sqrt(Integrate(InnerProduct(Riem_ex_3d, Riem_ex_3d),mesh3d)):.8f}"
    )

err Riemann curvature tensor 2d = 0.00000000
rel err Riemann curvature tensor 2d exact = 0.00009975
err Riemann curvature tensor 3d = 0.00000000
rel err Riemann curvature tensor 3d exact = 0.00162283

Curvature operator

Next, we look at the curvature operator \(\mathcal{Q}\), which incorporates the skew-symmetry property of the Riemann curvature tensor. It is defined via the identity \begin{align*} \mathcal{Q}(X^\flat\wedge Y^\flat, W^\flat \wedge Z^\flat) = \mathfrak{R}(X,Y,Z,W),\qquad \forall X,Y,Z,W\in\mathfrak{X}(M). \end{align*} Note the convention that the \(Z\) and \(W\) arguments are swapped and that e.g. \(X^\flat\wedge Y^\flat \in \Lambda^2(M)\) is a bi-vector. In two dimensions the curvature operator can be identified with a scalar function as \(\mathcal{Q}\,\omega(X,Y)\,\omega(W,Z) \simeq \mathcal{Q}(X^\flat\wedge Y^\flat, W^\flat\wedge Z^\flat)\) because the bi-vectors and the volume form are linearly dependent. Thus, in coordinates, there holds \begin{align*} \mathcal{Q} = -\frac{1}{4}\hat{\varepsilon}^{kl}\hat{\varepsilon}^{mn}\mathfrak{R}_{klmn} = -\frac{1}{4\,\det g}\varepsilon^{kl}\varepsilon^{mn}\mathfrak{R}_{klmn}, \end{align*} where \(\hat{\varepsilon}^{kl}= \frac{1}{\sqrt{\det g}}\varepsilon^{kl}\) denotes the Levi-Civita symbol with respect to the metric. In three dimensions, \(X^\flat\wedge Y^\flat\) can be identified with the cross product \((X\times Y)^\flat\). The curvature operator becomes a symmetric operator acting on two 1-forms and reads in coordinates \begin{align*} \mathcal{Q}^{ij} = -\frac{1}{4}\hat{\varepsilon}^{ikl}\hat{\varepsilon}^{jmn}\mathfrak{R}_{klmn} = -\frac{1}{4\,\det g}\varepsilon^{ikl}\varepsilon^{jmn}\mathfrak{R}_{klmn}. \end{align*} Using the coordinate expression of the Riemann curvature tensor and using the identity \(\hat{\varepsilon}^{ija}\hat{\varepsilon}_{kla}={\varepsilon}^{ija}{\varepsilon}_{kla}= \delta^i_k\delta^j_l-\delta^i_l\delta^j_k\) we have after a short calculation in 2D and 3D \begin{align*} \mathcal{Q}=-\frac{1}{2}\hat{\varepsilon}^{kl}\hat{\varepsilon}^{mn}\left(\partial_k\partial_m g_{ln}+\Gamma^q_{km}\Gamma_{lnq}\right),\qquad \mathcal{Q}^{ij}=-\frac{1}{2}\hat{\varepsilon}^{ikl}\hat{\varepsilon}^{jmn}\left(\partial_k\partial_m g_{ln}+\Gamma^q_{km}\Gamma_{lnq}\right). \end{align*}

On the opposite, we can express the Riemann curvature tensor in terms of the curvature operator \begin{align*} \mathfrak{R}_{ijkl}=\begin{cases} \hat{\varepsilon}_{ij}\hat{\varepsilon}_{kl}\mathcal{Q}, &\text{in 2D}\\ \hat{\varepsilon}_{ija}\hat{\varepsilon}_{klb}\mathcal{Q}^{ab}&\text{in 3D}.\end{cases} \end{align*} We recall that \(\hat{\varepsilon}_{ij} = \sqrt{\det g}\varepsilon_{ij}\).

WARNING: The curvature differential operator .Operator(“curvature”) of a Regge object returns \(\det g\,\mathcal{Q}\) instead of \(\mathcal{Q}\)! Thus, we need to divide by the determinant.

[3]:

Q1_2d = 1 / Det(g2d) * g2d.Operator("curvature")
Q2_2d = (
    -1 / (4 * Det(g2d)) * Einsum("ij,kl,ijkl", Eps2d, Eps2d, g2d.Operator("Riemann"))
)
Q3_2d = (
    -0.5
    / Det(g2d)
    * (g2d.Operator("inc") + Einsum("kl,mn,kmq,lnq", Eps2d, Eps2d, chr2_2d, chr1_2d))
)
Q_ex_2d = pcd.curvature

Q1_3d = 1 / Det(g3d) * g3d.Operator("curvature")
Q2_3d = (
    -1
    / (4 * Det(g3d))
    * Einsum("ikl,jmn,klmn->ij", Eps3d, Eps3d, g3d.Operator("Riemann"))
)
Q3_3d = (
    -0.5
    / Det(g3d)
    * (
        g3d.Operator("inc")
        + Einsum("ikl,jmn,kmq,lnq->ij", Eps3d, Eps3d, chr2_3d, chr1_3d)
    )
)
Q_ex_3d = sphere3.curvature

with TaskManager():
    print(
        f"err curvature operator 2d 1-2 = {sqrt(Integrate((Q1_2d - Q2_2d)**2,mesh2d)):.8f}"
    )
    print(
        f"err curvature operator 2d 1-3 = {sqrt(Integrate((Q1_2d - Q3_2d)**2,mesh2d)):.8f}"
    )
    print(
        f"rel err curvature operator 2d exact = {sqrt(Integrate((Q1_2d - Q_ex_2d)**2,mesh2d))/sqrt(Integrate(Q_ex_2d**2,mesh2d)):.8f}"
    )

    print(
        f"err curvature operator 3d 1-2 = {sqrt(Integrate(InnerProduct(Q1_3d - Q2_3d,Q1_3d - Q2_3d),mesh3d)):.8f}"
    )
    print(
        f"err curvature operator 3d 1-3 = {sqrt(Integrate(InnerProduct(Q1_3d - Q3_3d,Q1_3d - Q3_3d),mesh3d)):.8f}"
    )
    print(
        f"rel err curvature operator 3d exact = {sqrt(Integrate(InnerProduct(Q1_3d - Q_ex_3d,Q1_3d - Q_ex_3d),mesh3d))/sqrt(Integrate(InnerProduct(Q_ex_3d, Q_ex_3d),mesh3d)):.8f}"
    )

err curvature operator 2d 1-2 = 0.00000000
err curvature operator 2d 1-3 = 0.00000000
rel err curvature operator 2d exact = 0.00005460
err curvature operator 3d 1-2 = 0.00000000
err curvature operator 3d 1-3 = 0.00000000
rel err curvature operator 3d exact = 0.01268933

Ricci curvature tensor

The Ricci curvature tensor is defined by contracting two of the indices of the Riemann curvature tensor (i.e., taking a trace) \begin{align*} \mathrm{Ric}_{ij} = \mathfrak{R}_{kijl}g^{kl}, \end{align*} In terms of the curvature operator, we have \begin{align*} \mathrm{Ric} = \begin{cases} \mathcal{Q}\,g, &\text{in 2D}\\ g^{-1}\times \mathcal{Q}, &\text{in 3D}\end{cases}, \end{align*} where \((A \times B)_{ij}= \hat{\varepsilon}_{ikl}\hat{\varepsilon}_{jmn}A^{km}B^{ln}\) is the tensor cross product.

[4]:

Ric1_2d = g2d.Operator("Ricci")
Ric2_2d = 1 / Det(g2d) * g2d.Operator("curvature") * g2d
Ric_ex_2d = pcd.Ricci

Ric1_3d = g3d.Operator("Ricci")
Ric2_3d = Det(g3d) * Einsum(
    "ikl,jmn,km,ln->ij",
    Eps3d,
    Eps3d,
    Inv(g3d),
    1 / Det(g3d) * g3d.Operator("curvature"),
)
Ric_ex_3d = sphere3.Ricci

with TaskManager():
    print(
        f"err Ricci curvature tensor 2d = {sqrt(Integrate(InnerProduct(Ric1_2d - Ric2_2d,Ric1_2d - Ric2_2d),mesh2d)):.8f}"
    )
    print(
        f"rel err Ricci curvature tensor 2d exact = {sqrt(Integrate(InnerProduct(Ric1_2d - Ric_ex_2d,Ric1_2d - Ric_ex_2d),mesh2d))/sqrt(Integrate(InnerProduct(Ric_ex_2d, Ric_ex_2d),mesh2d)):.8f}"
    )
    print(
        f"err Ricci curvature tensor 3d = {sqrt(Integrate(InnerProduct(Ric1_3d - Ric2_3d,Ric1_3d - Ric2_3d),mesh3d)):.8f}"
    )
    print(
        f"rel err Ricci curvature tensor 3d exact = {sqrt(Integrate(InnerProduct(Ric1_3d - Ric_ex_3d,Ric1_3d - Ric_ex_3d),mesh3d))/sqrt(Integrate(InnerProduct(Ric_ex_3d, Ric_ex_3d),mesh3d)):.8f}"
    )

err Ricci curvature tensor 2d = 0.00000000
rel err Ricci curvature tensor 2d exact = 0.00008250
err Ricci curvature tensor 3d = 0.00000000
rel err Ricci curvature tensor 3d exact = 0.00883227

Scalar curvature

The scalar curvature is given by contracting the Riemann curvature tensor twice (or equivalently taking the trace of the Ricci curvature tensor) \begin{align*} S = \mathrm{Ric}_{ij} g^{ij} = \mathfrak{R}_{kijl}g^{kl}g^{ij}. \end{align*}

In two dimensions, the Gauss curvature is given for linear independent \(X,Y\) by \begin{align*} K = \frac{\mathcal{R}(X,Y,Y,X)}{g(X,X)g(Y,Y)-g(X,Y)^2}, \end{align*} which is independent of the choice of \(X,Y\). Further, it is related to the scalar curvature by \(S=2K\).

In two dimensions, the Gauss curvature \(K\) coincides with the curvature operator \(\mathcal{Q}\), \(K=\mathcal{Q}\). In 3D, the scalar curvature can be computed from the curvature operator by taking its trace \begin{align*} S = 2 \mathrm{tr}_g(\mathcal{Q}) = 2 g_{ij}\mathcal{Q}^{ij}. \end{align*}

[5]:

S1_2d = g2d.Operator("scalar")
S2_2d = Einsum("il,jk,ijkl", Inv(g2d), Inv(g2d), g2d.Operator("Riemann"))
S3_2d = 2 / Det(g2d) * g2d.Operator("curvature")
S_ex_2d = pcd.scalar

S1_3d = g3d.Operator("scalar")
S2_3d = Einsum("il,jk,ijkl", Inv(g3d), Inv(g3d), g3d.Operator("Riemann"))
S3_3d = 2 / Det(g3d) * InnerProduct(g3d, g3d.Operator("curvature"))
S_ex_3d = sphere3.scalar

with TaskManager():
    print(
        f"err scalar curvature 2d 1-2 = {sqrt(Integrate((S1_2d - S2_2d)**2,mesh2d)):.8f}"
    )
    print(
        f"err scalar curvature 2d 1-3 = {sqrt(Integrate((S1_2d - S3_2d)**2,mesh2d)):.8f}"
    )
    print(
        f"rel err scalar curvature 2d exact = {sqrt(Integrate((S1_2d - S_ex_2d)**2,mesh2d))/sqrt(Integrate(S_ex_2d**2,mesh2d)):.8f}"
    )

    print(
        f"err scalar curvature 3d 1-2 = {sqrt(Integrate((S1_3d - S2_3d)**2,mesh3d)):.8f}"
    )
    print(
        f"err scalar curvature 3d 1-3 = {sqrt(Integrate((S1_3d - S3_3d)**2,mesh3d)):.8f}"
    )
    print(
        f"rel err scalar curvature 3d exact = {sqrt(Integrate((S1_3d - S_ex_3d)**2,mesh3d))/sqrt(Integrate(S_ex_3d**2,mesh3d)):.8f}"
    )

err scalar curvature 2d 1-2 = 0.00000000
err scalar curvature 2d 1-3 = 0.00000000
rel err scalar curvature 2d exact = 0.00005460
err scalar curvature 3d 1-2 = 0.00000000
err scalar curvature 3d 1-3 = 0.00000000
rel err scalar curvature 3d exact = 0.05574200

Einstein tensor

The Einstein tensor is defined by \begin{align*} G_{ij}=\mathrm{Ric}_{ij}-\frac{1}{2}S\,g_{ij},\qquad G = J\mathrm{Ric}, \end{align*} where \(J\sigma = \sigma-\frac{1}{2}\mathrm{tr}(\sigma)\,g\). In three dimensions, it is the negative of the double lowering of the curvature operator \begin{align*} G = -\mathcal{Q}^{\flat\flat},\qquad G_{ij} = -g_{ik}\mathcal{Q}^{kl}g_{lj}. \end{align*}

[6]:

G1_3d = g3d.Operator("Einstein")
G2_3d = g3d.Operator("Ricci") - 0.5 * g3d.Operator("scalar") * g3d
G3_3d = -g3d * (1 / Det(g3d) * g3d.Operator("curvature")) * g3d
G_ex_3d = sphere3.Einstein

with TaskManager():
    print(
        f"err Einstein tensor 3d 1-2 = {sqrt(Integrate(InnerProduct(G1_3d - G2_3d,G1_3d - G2_3d),mesh3d)):.8f}"
    )
    print(
        f"err Einstein tensor 3d 1-3 = {sqrt(Integrate(InnerProduct(G1_3d - G3_3d,G1_3d - G3_3d),mesh3d)):.8f}"
    )
    print(
        f"rel err Einstein tensor 3d exact = {sqrt(Integrate(InnerProduct(G1_3d - G_ex_3d,G1_3d - G_ex_3d),mesh3d))/sqrt(Integrate(InnerProduct(G_ex_3d, G_ex_3d),mesh3d)):.8f}"
    )

err Einstein tensor 3d 1-2 = 0.00000000
err Einstein tensor 3d 1-3 = 0.00000000
rel err Einstein tensor 3d exact = 0.13515517

[ ]: