"""Curvature tensors of general relativity, as a small symbolic library. Specification given (28 Oct 2026): "Write the Riemann curvature tensor as part of a Python library: classes for the metric, the Christoffel symbols and the Riemann tensor, with documentation, conversion between upper and lower index positions, and Lorentz transformations." Written by Claude Opus 5.5 from that specification. Unlike the first version, the author had already read EinsteinPy's code earlier in the same session. Conventions ----------- Index positions are written as a string, one letter per index: 'u' for an upper (contravariant) index and 'l' for a lower (covariant) one. The Riemann tensor is stored as R^rho_{sigma mu nu} ('ulll') with R^rho_{sigma mu nu} = d_mu Gamma^rho_{nu sigma} - d_nu Gamma^rho_{mu sigma} + Gamma^rho_{mu lam} Gamma^lam_{nu sigma} - Gamma^rho_{nu lam} Gamma^lam_{mu sigma} """ from itertools import product import sympy as sp class Tensor: """A tensor of any rank with symbolic components and named index positions. Parameters ---------- components : sympy.Array or nested list The components, one axis per index. coords : list of sympy.Symbol The coordinates, e.g. [t, r, theta, phi]. config : str Index positions, 'u' (upper) or 'l' (lower) per index. metric : Metric, optional The metric used to raise and lower indices. """ def __init__(self, components, coords, config, metric=None): self.components = sp.MutableDenseNDimArray(components) self.coords = list(coords) self.config = config self.metric = metric if self.components.rank() != len(config): raise ValueError("config must name one position per index") if set(config) - {"u", "l"}: raise ValueError("config may contain only 'u' and 'l'") if any(size != len(self.coords) for size in self.components.shape): raise ValueError("every index must run over all coordinates") @property def dims(self): """Number of space-time dimensions.""" return len(self.coords) @property def rank(self): """Number of indices.""" return len(self.config) def __getitem__(self, index): return self.components[index] def simplify(self): """Simplify every component in place and return the tensor.""" for index in product(range(self.dims), repeat=self.rank): self.components[index] = sp.simplify(self.components[index]) return self def change_config(self, new_config, metric=None): """Return a copy with indices raised or lowered to match new_config. Lowering uses g_{ab}; raising uses the inverse metric g^{ab}. """ metric = metric or self.metric if metric is None: raise ValueError("a metric is needed to raise or lower indices") if len(new_config) != self.rank: raise ValueError("new_config must have one letter per index") result = self.components for position, (old, new) in enumerate(zip(self.config, new_config)): if old == new: continue converter = metric.lower if new == "l" else metric.inverse result = _contract_index(result, converter, position, self.dims) return type(self)._from_parts(result, self.coords, new_config, metric) def lorentz_transform(self, matrix): """Return the components in a frame related by the matrix Lambda. Upper indices transform with Lambda^a'_b, lower indices with the inverse, so contractions such as scalars are unchanged. """ matrix = sp.Matrix(matrix) if matrix.shape != (self.dims, self.dims): raise ValueError("the transformation must be a dims x dims matrix") inverse = matrix.inv() result = self.components for position, kind in enumerate(self.config): # new upper index: Lambda^a'_b T^b ; new lower index: T_b (Lambda^-1)^b_a' converter = matrix if kind == "u" else inverse.T result = _contract_index(result, converter, position, self.dims) new_metric = self.metric.lorentz_transform(matrix) if self.metric else None return type(self)._from_parts(result, self.coords, self.config, new_metric) @classmethod def _from_parts(cls, components, coords, config, metric): obj = Tensor.__new__(cls) Tensor.__init__(obj, components, coords, config, metric) return obj def _contract_index(array, matrix, position, dims): """Contract one index of array with a matrix: A'[..a..] = M[a, b] A[..b..].""" rank = array.rank() result = sp.MutableDenseNDimArray.zeros(*([dims] * rank)) for index in product(range(dims), repeat=rank): total = 0 for b in range(dims): source = index[:position] + (b,) + index[position + 1:] term = array[source] if term != 0: total += matrix[index[position], b] * term result[index] = total return result class Metric(Tensor): """The space-time metric g_{ab}, stored with both indices lower.""" def __init__(self, matrix, coords): matrix = sp.Matrix(matrix) if matrix.shape != (len(coords), len(coords)): raise ValueError("the metric must be square and match the coordinates") if matrix != matrix.T: raise ValueError("the metric must be symmetric") super().__init__(sp.Array(matrix.tolist()), coords, "ll", metric=None) self.metric = self self.lower = matrix self.inverse = matrix.inv() def lorentz_transform(self, matrix): """Return the metric in the transformed frame: g' = L^-T g L^-1.""" inverse = sp.Matrix(matrix).inv() return Metric(sp.simplify(inverse.T * self.lower * inverse), self.coords) class ChristoffelSymbols(Tensor): """Christoffel symbols of the second kind, Gamma^a_{bc} ('ull').""" @classmethod def from_metric(cls, metric): """Compute Gamma^a_{bc} = 1/2 g^{ad} (d_b g_{dc} + d_c g_{db} - d_d g_{bc}).""" n, g, ginv, x = metric.dims, metric.lower, metric.inverse, metric.coords gamma = sp.MutableDenseNDimArray.zeros(n, n, n) for a, b in product(range(n), repeat=2): for c in range(b, n): # symmetric in the lower indices value = sum( ginv[a, d] * (sp.diff(g[d, c], x[b]) + sp.diff(g[d, b], x[c]) - sp.diff(g[b, c], x[d])) for d in range(n) ) / 2 gamma[a, b, c] = gamma[a, c, b] = sp.simplify(value) return cls(gamma, x, "ull", metric) class RiemannTensor(Tensor): """The Riemann curvature tensor R^rho_{sigma mu nu} ('ulll').""" @classmethod def from_metric(cls, metric): """Compute the Riemann tensor of a metric.""" return cls.from_christoffels(ChristoffelSymbols.from_metric(metric)) @classmethod def from_christoffels(cls, gamma): """Compute the Riemann tensor from Christoffel symbols Gamma^a_{bc}.""" n, x = gamma.dims, gamma.coords riemann = sp.MutableDenseNDimArray.zeros(n, n, n, n) for rho, sigma, mu in product(range(n), repeat=3): for nu in range(mu + 1, n): # antisymmetric in the last two indices value = ( sp.diff(gamma[rho, nu, sigma], x[mu]) - sp.diff(gamma[rho, mu, sigma], x[nu]) + sum( gamma[rho, mu, lam] * gamma[lam, nu, sigma] - gamma[rho, nu, lam] * gamma[lam, mu, sigma] for lam in range(n) ) ) value = sp.simplify(value) riemann[rho, sigma, mu, nu] = value riemann[rho, sigma, nu, mu] = -value return cls(riemann, x, "ulll", gamma.metric) def kretschmann(self): """The Kretschmann scalar R_{abcd} R^{abcd}, the same in every frame.""" lower = self.change_config("llll") upper = self.change_config("uuuu") total = sum( lower[index] * upper[index] for index in product(range(self.dims), repeat=4) if lower[index] != 0 ) return sp.simplify(total)