# Discrete Gravity with Local Lorentz Invariance

Eugene Kur  Lawrence Livermore National Laboratory, Livermore, CA 94550
Alexander S. Glasser  Princeton Plasma Physics Laboratory, Princeton
University, Princeton, New Jersey 08543
Department of Astrophysical Sciences, Princeton University, Princeton, New
Jersey 08544

###### Abstract

A novel structure-preserving algorithm for general relativity in vacuum is
derived from a lattice gauge theoretic discretization of the tetradic Palatini
action. The resulting model of discrete gravity is demonstrated to preserve
local Lorentz invariance and symplectic structure.

##  I Introduction

Since at least the 1990s, structure-preserving algorithms [1] have flourished
in computational physics, having found wide adoption in subfields as diverse
as orbital mechanics [2, 3, 4, 5], geophysics [6, 7] and plasma physics [8, 9,
10, 11, 12, 13, 14, 15, 16]. Such algorithms are generally derived from a
Lagrangian or Hamiltonian formalism and use discretizations that preserve the
symplectic structure, topology, gauge symmetry, and conservation laws of their
underlying physical systems. This preservation of mathematical structure can
substantially improve the accuracy and fidelity of numerical simulations.

Structure-preserving discretizations of general relativity (GR) arguably have
an even longer history. The most widely explored such approach was introduced
in 1961: Regge calculus [17] is a discrete variational approximation of GR
that encodes spacetime data on a simplicial mesh. In four spacetime
dimensions, Regge calculus elegantly approximates the Einstein-Hilbert action
by a sum over areas $A_{h}$ and _deficit angles_ $\delta_{h}$, such that

$\displaystyle\begin{split}S_{\text{Regge}}=\sum\limits_{h}A_{h}\delta_{h}~{}\xrightarrow{A_{h}\rightarrow
0}~{}S_{\text{EH}}=\frac{1}{2}\int\mathrm{d}^{4}x\sqrt{-g}R\end{split}$ (1)

in the continuum limit. Here, $h$ labels each 2-simplex (i.e. triangle) of the
simplicial complex, and $\delta_{h}$ describes the failure of the 4-simplices
adjoining $h$—i.e. ${\\{\sigma^{4}~{}|~{}\sigma^{4}\supset h\\}}$—to tesselate
their embedding in flat $\mathbb{R}^{4}$ spacetime [18].

Since the 1970s, Regge calculus has not only been actively employed as the
basis of many studies in quantum gravity (e.g. [19, 20, 21, 22, 23, 24, 25,
26]), but also as an algorithmic approach to classical numerical relativity
(e.g. [27, 28, 29, 30, 31, 32, 33, 34, 35]). Despite its success as a
numerical tool, however, most studies in numerical relativity continue to
depend upon standard finite difference methods. Two reasons cited for this
include the need to develop (i) a description of matter in Regge calculus, as
well as (ii) a better understanding of its relationship to standard methods in
numerical relativity [36, 37].

In particular, because the degrees of freedom of Regge calculus are quite
distinct from those of continuum GR, it can be challenging to initialize a
Regge calculus simulation with known GR initial conditions, or to test whether
a particular simulation using Regge calculus recovers a known GR solution.
Although various physical solutions have indeed been thoroughly and
successfully benchmarked with Regge calculus [38, 39], it would seem that any
given simulation generally requires a bespoke understanding of the map between
discrete and continuum degrees of freedom.

It is also worth emphasizing that, despite Regge calculus being a variational
method, it nonetheless forfeits—in its complete, nonperturbative
formulation—the local gauge symmetry of GR [23]. While local gauge symmetry is
maintained in a Regge calculus description of flat spacetime—and even in a
linearized Regge calculus of curved spacetimes [40]—this structural feature of
GR is at best only partially preserved overall.

In this paper, an alternative variational approach to simulating general
relativity is developed that ameliorates some of these limitations. Our effort
employs familiar tools of lattice gauge theory [41] to construct a structure-
preserving discretization of the tetradic Palatini action [42]. Using a
Poincaré group-valued connection derived from Cartan geometry, we describe a
novel variational algorithm for numerical relativity that exactly preserves
Lorentz gauge symmetry. We further show the algorithm is (multi)symplectic,
with a symplectic structure analogous to continuum GR.

The approach we take is closely related to Poincaré gauge theoretic studies of
lattice quantum gravity by Menotti, et al. [43, 44]. To our knowledge,
however, the classical physics of these methods, including their equations of
motion, for example, have not previously been explored, nor have they been
extended to define an algorithm for numerical relativity. Moreover, our
construction is general to simplicial and cubical discretizations of
spacetime, and we develop a streamlined construction of the aforementioned
Poincaré connection.

The remainder of this paper is organized as follows: Section II briefly
reviews the tetradic Palatini action and its origins in Cartan geometry;
Section III derives a discretization of this action in a manner that preserves
Lorentz gauge invariance; Section IV derives the discrete, classical equations
of motion that comprise the algorithm; and Section V describes its symplectic
structure. Finally, Section VI summarizes and concludes.

##  II The Tetradic Palatini Action in Continuous Spacetime

Let us first review the tetradic Palatini action in the continuum. We consider
a four-dimensional Lorentzian spacetime with connection, denoted
${(M,g,\Gamma)}$, and employ the following index conventions:

  1. (i)

spacetime coordinate indices $\\{\mu,\nu,\dots\\}$ are raised and lowered by
${g_{\mu\nu}}$, the metric on $M$;

  2. (ii)

internal Lorentz indices $\\{A,B,\dots\\}$ are raised and lowered by
${\eta_{AB}}$, the Minkowski metric; and

  3. (iii)

any other indices $\\{a,b,\dots\\}$ will be specified as needed.

In a coordinate basis ${\\{\partial_{\mu}\\}}$, the affine connection $\Gamma$
has components
${\Gamma^{\sigma}_{\mu\nu}=\mathrm{d}x^{\sigma}(\nabla_{\partial_{\mu}}\partial_{\nu})}$.

Up to local Lorentz gauge, the metric $g$ uniquely determines a tetrad field
$e$ on $M$, a vector-valued 1-form with components
${e^{A}=e^{A}_{\mu}\mathrm{d}x^{\mu}\in\Gamma(T^{*}M)}$ defined to satisfy

$\displaystyle\begin{split}g_{\mu\nu}=e^{A}_{\mu}\eta_{AB}e^{B}_{\nu}.\end{split}$
(2)

Since $g_{\mu\nu}$ is non-degenerate, ${e^{A}_{\mu}(p)}$ defines $\forall$
${p\in M}$ an isomorphism between the tangent space $T_{p}M$ and the ‘internal
Lorentz space’ at $p$. As a result, any vector field ${X\in\Gamma(TM)}$ can be
equally well described in terms of the Lorentz frame
${\\{\partial_{A}=e^{\mu}_{A}\partial_{\mu}\\}}$ such that
${X=X^{A}\partial_{A}=X^{\mu}\partial_{\mu}}$. (Here,
${e^{\mu}_{A}e^{A}_{\nu}=\delta^{\mu}_{\nu}}$ defines a matrix inverse.) In
general, ${\\{\partial_{A}\\}}$ is a non-coordinate basis (since the
commutator ${[\partial_{A},\partial_{B}]}$ need not vanish) and is dual to
${\\{e^{A}\\}}$.

Parallel transport may be defined in the Lorentz frame by the 1-form spin
connection $\omega$, with components ${\omega^{A}_{~{}B}=\omega^{A}_{\mu
B}\mathrm{d}x^{\mu}\in\Gamma(T^{*}M)}$, such that

$\displaystyle\begin{split}\nabla_{\mu}X^{A}&=\partial_{\mu}X^{A}+\omega^{A}_{\mu
B}X^{B}.\end{split}$ (3)

Since the Lorentz frame arises, ultimately, as a change of basis, the spin
connection components have a definite relation to
${\Gamma^{\sigma}_{\mu\nu}}$. In particular, ${\omega^{A}_{\mu
B}=e^{A}(\nabla_{\partial_{\mu}}\partial_{B})}$, which can be more
suggestively expanded as

$\displaystyle\begin{split}\nabla_{\mu}e^{A}_{\nu}=\partial_{\mu}e^{A}_{\nu}+\omega^{A}_{\mu
B}e^{B}_{\nu}-\Gamma^{\sigma}_{\mu\nu}e^{A}_{\sigma}=0.\end{split}$ (4)

The name ‘Lorentz frame’ can be justified by requiring $\eta_{AB}$ to be
invariant under parallel transport—

$\displaystyle\begin{split}0=\nabla_{\mu}\eta_{AB}&=\partial_{\mu}\eta_{AB}-\omega^{C}_{\mu
A}\eta_{CB}-\omega^{C}_{\mu B}\eta_{AC}\\\ &=-(\omega_{\mu BA}+\omega_{\mu
AB}).\end{split}$ (5)

Due to its resulting antisymmetry, $\omega$ is defined by this condition as an
${\mathfrak{so}(3,1)}$-valued 1-form. Studying Eqs. (4) and (5), we also see
that the metric compatibility of ${\Gamma}$ follows immediately from this
$\eta$-compatibility of $\omega$.

In Einstein GR, $\Gamma$ is assumed to be not only metric-compatible but
torsion-free, such that ${\Gamma^{\sigma}_{[\mu\nu]}=0}$. The resulting Levi-
Civita connection $\Gamma_{\text{LC}}$ is then uniquely determined by $g$.
Thus, as a further consequence of vanishing torsion, by Eqs. (2) and (4) the
metric $g$ also uniquely determines—up to Lorentz gauge—the ${\mathbb{R}^{4}}$
and ${\mathfrak{so}(3,1)}$-valued 1-forms $e$ and $\omega$, respectively, on
$M$.

Conversely, the metric $g$ and the connection $\Gamma_{\text{LC}}$ can be
uniquely recovered from the fields $e$ and $\omega$ on a torsion-free
manifold. (More precisely, $g$ can be canonically recovered up to an overall
constant factor.) There is an equivalence, therefore, between a Lorentzian
manifold ${(M,g,\Gamma_{\text{LC}})}$ and its torsion-free _Cartan geometric_
counterpart ${(M,e,\omega)}$ [45]. Let us describe the origin of this
nomenclature.

The 1-forms $e$ and $\omega$ are more economically regarded as components of
the _Cartan connection_ ${A=A_{\mu}\mathrm{d}x^{\mu}}$ on $M$, defined by

$\displaystyle\begin{split}A=\left[\begin{matrix}\omega&e\\\
\mathbf{0}&0\end{matrix}\right]\in\Gamma(\mathfrak{p}\otimes
T^{*}M).\end{split}$ (6)

Here,
${\mathfrak{p}={\mathfrak{so}(3,1)}\ltimes\mathbb{R}^{4}\subset\mathfrak{gl}_{5}(\mathbb{R})}$
denotes the Lie algebra of the Poincaré group, so that $A$ is a
$\mathfrak{p}$-valued 1-form on $M$.111A $\mathfrak{p}$-valued Cartan
connection $A$ is formally defined on an ${SO(3,1)}$-principal bundle $P$ over
$M$ such that ${A:T_{p}P\rightarrow\mathfrak{p}}$ is an isomorphism $\forall$
${p\in P}$. The pair ${(P,A)}$ defines a _Cartan geometry_ [45]. In physical
applications, however, $A$ is conventionally defined by its pullback to $M$
and its overlying bundle is elided. Following the previous discussion, a
torsion-free Lorentzian manifold can be equivalently defined by its metric $g$
or $\mathfrak{p}$-valued Cartan connection $A$, and solving for the dynamics
of $A$ similarly determines dynamics for $g$. In what follows, we therefore
regard the tetradic Palatini action as a dynamical theory of the Cartan
connection $A$ on $M$.

To that end, we first recall the curvature 2-form of the Cartan connection,
defined as

$\displaystyle\begin{split}F=\mathrm{d}A+A\wedge A=\left[\hskip
2.0pt\begin{matrix}[c|c]\mathrm{d}\omega+\omega\wedge\omega&\mathrm{D}e\\\
\hline\cr\mathbf{0}&0\end{matrix}\hskip
2.0pt\right]=\left[\begin{matrix}R&T\\\
\mathbf{0}&0\end{matrix}\right],\end{split}$ (7)

where ${\mathrm{D}e=\mathrm{d}e+\omega\wedge e}$ denotes the exterior
covariant derivative of $e$. In components,

$\displaystyle\begin{split}R^{A}_{~{}B\mu\nu}&=\partial_{\mu}\omega^{A}_{\nu
B}-\partial_{\nu}\omega^{A}_{\mu B}+\omega^{A}_{\mu C}\omega^{C}_{\nu
B}-\omega^{A}_{\nu C}\omega^{C}_{\mu B}\end{split}$ (8)

denotes the _Lorentz curvature_ ${R^{A}_{~{}B}\in\Gamma(\wedge^{2}T^{*}M)}$,
while

$\displaystyle\begin{split}T^{A}_{\mu\nu}=(\mathrm{D}e^{A})_{\mu\nu}=\partial_{\mu}e^{A}_{\nu}-\partial_{\nu}e^{A}_{\mu}+\omega^{A}_{\mu
B}e^{B}_{\nu}-\omega^{A}_{\nu B}e^{B}_{\mu}\end{split}$ (9)

denotes the _torsion_ ${T^{A}\in\Gamma(\wedge^{2}T^{*}M)}$. As previously
noted, torsion is assumed to vanish, ${T^{A}=0}$, a priori in Einstein GR. In
the tetradic Palatini theory, however, torsion does not vanish by assumption,
but rather as a dynamical consequence of the action varied in vacuum, as we
presently demonstrate.

The 4-form Lagrangian ${\mathcal{L}_{\text{Pal}}\in\Gamma(\wedge^{4}T^{*}M)}$
of the tetradic Palatini action
${\mathcal{S}_{\text{Pal}}=\int_{M}\mathcal{L}_{\text{Pal}}}$ is defined in
terms of the tetrad $e$ and spin connection $\omega$ by [42]

$\displaystyle\begin{split}\mathcal{L}_{\text{Pal}}&=\epsilon_{ABCD}\left(e^{A}\wedge
e^{B}\wedge R^{CD}\right).\end{split}$ (10)

It should be noted that Eq. (10) is often called the _Einstein-Cartan-Sciama-
Kibble_ (ECSK) action [46, 47]. However, because ECSK theory prioritizes the
role of torsion in gravity, whereas we will pursue only the torsionless vacuum
equations of Einstein-Cartan gravity, we prefer the nomenclature _tetradic
Palatini action_. We note that the Lorentz invariance of
$\mathcal{L}_{\text{Pal}}$ follows directly from the $SO(3,1)$ invariance of
the Levi-Civita symbol ${\epsilon_{ABCD}}$, that is,

$\displaystyle\begin{split}(\epsilon_{ABCD})^{\prime}=\epsilon_{EFGH}\Lambda^{E}_{~{}A}\Lambda^{F}_{~{}B}\Lambda^{G}_{~{}C}\Lambda^{H}_{~{}D}=\epsilon_{ABCD}\bm{\mathsfit{det}}[\Lambda]\end{split}$
(11)

where ${\bm{\mathsfit{det}}[\Lambda]=1}$.

Unlike the fields ${(g,\Gamma_{\text{LC}})}$ of the Einstein-Hilbert action,
$(e,\omega)$ are taken to be independent in Eq. (10), and varied accordingly.
The variation of each field yields the respective equations of motion [48]

$\displaystyle\begin{split}(\delta e):\hskip
20.0pt&0=\epsilon_{ABCD}e^{B}\wedge R^{CD}\\\ (\delta\omega):\hskip
20.0pt&0=\epsilon_{ABCD}\mathrm{D}(e^{A}\wedge e^{B}).\end{split}$ (12)

Here, we note that ${\mathrm{D}(e^{A}\wedge e^{B})=\mathrm{D}e^{A}\wedge
e^{B}-e^{A}\wedge\mathrm{D}e^{B}}$. Taking
${R^{CD}=\frac{1}{2}R^{CD}_{~{}GH}e^{G}\wedge e^{H}}$ and
${T^{A}=\frac{1}{2}T^{A}_{GH}e^{G}\wedge e^{H}}$, (well-defined expansions for
$e^{A}$ nondegenerate), it is readily established that the former relation
${(\delta e)}$ of Eq. (12) yields Einstein’s vacuum field equations while
${(\delta\omega)}$ yields a zero torsion condition. In particular, since
${\epsilon_{ABCD}e^{A}\wedge e^{G}\wedge e^{H}\wedge
e^{I}=\delta^{GHI}_{BCD}e\bm{\mathsfit{vol}}}$ for a volume form
$\bm{\mathsfit{vol}}$ and ${e=\bm{\mathsfit{det}}[e^{A}_{\mu}]}$, the wedge
product ${(\delta e)\wedge e^{I}}$ gives

$\displaystyle\begin{split}0&=-\frac{1}{2}R^{CD}_{~{}~{}~{}GH}\delta^{GHI}_{ACD}=2R^{HI}_{~{}~{}~{}HA}-R^{GH}_{~{}~{}~{}GH}\delta^{I}_{A},\end{split}$
(13)

which, using Eq. (4), can be demonstrated equivalent to Einstein’s vacuum
equations, ${0=2R_{\mu\nu}-Rg_{\mu\nu}}$. Likewise, ${(\delta\omega)\wedge
e^{I}}$ gives

$\displaystyle\begin{split}0=-\delta^{GHI}_{ACD}T^{A}_{GH}&=2\left(T^{A}_{CA}\delta^{I}_{D}+T^{A}_{AD}\delta^{I}_{C}-T^{I}_{CD}\right).\end{split}$
(14)

Tracing over Eq. (14) with $\delta^{D}_{I}$ in four dimensions leaves
${0=T^{A}_{CA}}$. By Eq. (14), therefore, ${T^{I}_{CD}=0}$ in all components,
and torsion vanishes as desired.

Thus, despite making fewer initial assumptions, the tetradic Palatini action
nevertheless recovers the equations of motion of GR in vacuum; the dynamics of
the Cartan connection indeed recover those of GR.

Before concluding our discussion of continuous spacetime, the following will
be useful for the next section, which discretizes $\mathcal{S}_{\text{Pal}}$.
Evaluated on a 4-tuple of vector
fields—${\mathbf{X}=(X_{1},X_{2},X_{3},X_{4})}$, ${X_{a}\in\Gamma(TM)}$—the
4-form $\mathcal{L}_{\text{Pal}}$ of Eq. (10) yields

$\displaystyle\begin{split}\mathcal{L}_{\text{Pal}}(\mathbf{X})&=\frac{1}{2}\epsilon_{ABCD}(e^{A}_{\mu}e^{B}_{\nu}R^{CD}_{\sigma\tau})\epsilon^{abcd}X_{a}^{\mu}X_{b}^{\nu}X_{c}^{\sigma}X_{d}^{\tau}\\\
&=\frac{1}{2}\epsilon_{ABCD}(e^{A}_{\mu}e^{B}_{\nu}R^{CD}_{\sigma\tau})\epsilon^{\mu\nu\sigma\tau}\bm{\mathsfit{det}}[\mathbf{X}],\end{split}$
(15)

where the function $\bm{\mathsfit{det}}[\mathbf{X}]$ is the matrix determinant
of the 4-tuple, expressed in the coordinate basis induced by ${\\{x^{\mu}\\}}$
and evaluated pointwise over $M$.

##  III The Discrete Action

We now discretize the tetradic Palatini action of the previous section by
methodically mapping its continuum degrees of freedom to their discrete
counterparts on a lattice. Our formalism will be general to orientable
simplicial and cubical discretizations and we take care to preserve the
theory’s Lorentz invariance.

To proceed, we first choose a coordinate chart on the continuum spacetime
manifold $M$, and construct a lattice (simplicial or cubical) on its
coordinate space in $\mathbb{R}^{4}$. As such, the lattice inherits the
Euclidean geometry of the coordinate space (such as straight edges and flat
faces), but this ‘lattice geometry’ will play no role in our description of
spacetime. Topological features of $M$ must be retained in the construction of
the lattice, including via possible identifications of its edges or faces. In
such a case, the lattice should be regarded as only locally embedded in
coordinate space while being globally homeomorphic to the target spacetime
manifold. Such a construction is standard in the triangulation of manifolds
(see e.g. [49]).

To establish notation for lattice degrees of freedom, we denote the set of
lattice $k$-cells by ${\Sigma^{k}=\\{\sigma^{k}\\}}$, such that an arbitrary
oriented $k$-cell will be denoted $\sigma^{k}$, or will otherwise be specified
by an ordered label of its vertices. ${\sigma_{ij}\in\Sigma^{1}}$, for
example, denotes an edge oriented from vertex ${\sigma_{i}\in\Sigma^{0}}$ to
vertex ${\sigma_{j}\in\Sigma^{0}}$. We define ${N_{i}(\sigma^{k})=\\{j\neq
i~{}|~{}\sigma_{ij}\subset\sigma^{k}\\}}$ as the set of labels of neighboring
vertices in the cell $\sigma^{k}$ that share an edge with basepoint
$\sigma_{i}$. In both simplicial and cubical discretizations in four
dimensions, for example, ${\\#N_{i}(\sigma^{4})=4}$ if
${\sigma_{i}\subset\sigma^{4}}$ and 0 otherwise. We denote the permutation set
of these neighboring vertex labels as
${\Pi_{i}(\sigma^{4})=S[N_{i}(\sigma^{4})]}$.

As described in Section II, the geometric information of $M$ is encoded in its
Cartan connection—the fields $e^{A}_{\mu}(x)$ and $\omega^{A}_{\mu
B}(x)$—which may be regarded as defined on coordinate space. A natural first
(provisional) approximation of the tetradic Palatini action, as it appears in
Eq. (15), therefore follows by mapping these fields to the lattice, such that

$\displaystyle\begin{split}\mathcal{S}_{\text{Pal}}&=\sum_{\sigma^{4}\in\Sigma^{4}}\int_{\sigma^{4}}\mathcal{L}_{\text{Pal}}(x)\mathrm{d}^{4}x\\\
&\approx\sum_{\begin{subarray}{c}\sigma^{4}\in\Sigma^{4}\\\
\sigma_{i}\in\sigma^{4}\end{subarray}}\frac{(-1)^{|\pi|}V_{f}}{2n_{v}}\epsilon^{\mu\nu\sigma\tau}\epsilon_{ABCD}\Big{(}e^{A}_{\mu}e^{B}_{\nu}R^{CD}_{\sigma\tau}\Big{)}\Big{|}_{\sigma_{i}}\bm{\mathsfit{det}}[\mathbf{V}_{\sigma_{i}}]\end{split}$
(16)

where ${\mathcal{L}_{\text{Pal}}(x)}$ denotes $\mathcal{L}_{\text{Pal}}$ on
coordinate space. The second line approximates the integral over
${\sigma^{4}}$ by averaging the value of its integrand, as expressed in Eq.
(15), at each of its vertices. More specifically:

  * •

$n_{v}$ denotes the number of vertices in $\sigma^{4}$ over which the
integrand is averaged. ${n_{v}=5}$ on a simplicial lattice and ${n_{v}=16}$ on
a cubical lattice.

  * •

${\mathbf{V}_{\sigma_{i}}=(V_{i1},V_{i2},V_{i3},V_{i4})}$ is a 4-tuple of
‘edge vectors’ emanating from vertex $\sigma_{i}$. These point in the
directions of neighboring vertices in $\sigma^{4}$, with magnitudes set by the
edges’ coordinate lengths.

  * •

$(-1)^{|\pi|}$ accounts for the relative orientation between the 4-tuple
$\mathbf{V}_{\sigma_{i}}$ and the cell $\sigma^{4}$, whose orientation is
inherited from $M$. This factor is expressed in terms of a permutation $\pi$
to be defined more concretely below.

  * •

The volume factor $V_{f}$ corrects for the fact that
${\bm{\mathsfit{det}}[\mathbf{V}_{\sigma_{i}}]}$ implicitly evaluates
${\mathcal{L}_{\text{Pal}}}$ not on $\sigma^{4}$, but on a (hyper-)
parallelapiped specified by $\mathbf{V}_{\sigma_{i}}$. On a cubical lattice,
these volumes coincide and ${V_{f}=1}$, but on a simplicial lattice,
${\bm{\mathsfit{det}}[\mathbf{V}_{\sigma_{i}}]}$ overcounts the volume of
$\sigma^{4}$ by the ratio of a normalized hypercube to one of its corners,
such that ${V_{f}}=1/4!$.

This approximation of $\int_{\sigma^{4}}$ by the average of vertex evaluations
is, in effect, a second order accurate multi-dimensional trapezoid rule (see
[50] and references therein).

Eq. (16) instructively approximates the continuum action, but it is
insufficient to determine dynamics for a discrete theory. In particular, Eq.
(16) discretely samples degrees of freedom that are manifestly defined in the
continuum—e.g. ${\partial\omega}$. (If this continuous derivative were avoided
by regarding $R$ itself as a Lie-algebra-valued degree of freedom, rather than
$\omega$, the resulting action would not recover the equations of motion of
Einstein GR.) To reformulate Eq. (16) with bona fide discrete degrees of
freedom, we now proceed hewing more closely to the underlying differential
geometry of the tetradic Palatini action.

In particular, the study of structure-preserving discretizations (such as
discrete exterior calculus (DEC) [51] and finite element exterior calculus
(FEEC) [52, 53]) has demonstrated the importance of preserving the degrees of
discrete differential forms. Therefore, rather than sampling continuum fields
at vertices, as we do in Eq. (16), we will instead map 1-forms to data
associated with edges, and 2-forms to data associated with faces.

However, an additional challenge we must overcome is the gauge-dependent
character of the fields we are modeling, which thwarts conventional approaches
such as DEC and FEEC. The preservation of Lorentz invariance in our theory
will require that we express discrete fields in a definite (if arbitrary)
Lorentz gauge, which is associated in a continuum theory to each point of
spacetime, and in a discrete theory to each vertex. This pointwise gauge
selection is in tension with the desire to characterize differential forms
over edges and faces of finite extent. For example, a scalar-valued 1-form is
conventionally approximated on an edge by its integral over that edge. Here,
such an integral involves a continuum of different gauge choices in spacetime
that prevent the simple summation of fields defined at disparate points.

A resolution to this tension is naturally found in the holonomy of a
connection. The path-ordered integral of a 1-form connection produces the
means to parallel transport between different gauge choices. It is an object
that can be naturally associated with an edge, and which by construction
accounts for a difference in gauge between two vertices. In this sense, the
Cartan connection—which retains the geometric data of a Lorentzian manifold—
provides a natural approach to a structure-preserving discretization of GR.

To map the Cartan connection $A$ on $M$ to holonomies on the discrete lattice
in coordinate space, we associate to each edge $\sigma_{ij}$ the following
path-ordered integral:

$\displaystyle\begin{split}U_{ij}=\mathcal{P}\left\\{\exp\int_{\sigma_{ij}}A\right\\}=\mathcal{P}\left\\{\exp\int_{\sigma_{ij}}\left[\begin{matrix}\omega&e\\\
\mathbf{0}&0\end{matrix}\right]\right\\}=\left[\begin{matrix}\Lambda_{ij}&\ell_{ij}\vspace{2pt}\\\
\mathbf{0}&1\end{matrix}\right].\end{split}$ (17)

This is a standard construction of lattice gauge theory [41]. $U_{ij}$
constitutes the Poincaré group-valued holonomy associated with edge
$\sigma_{ij}$, expressed in the representation
${SO(3,1)\ltimes\mathbb{R}^{4}\subset GL_{5}(\mathbb{R})}$ and characterized
by Lorentz and translation group elements, ${\Lambda_{ij}\in SO(3,1)}$ and
${\ell_{ij}\in\mathbb{R}^{4}}$, respectively.

We denote the ${(A,B)^{\text{th}}}$ component of the Lorentz connection along
edge $\sigma_{ij}$ by $\Lambda^{A}_{ijB}$, and the ${A^{\text{th}}}$ component
of the corresponding translation connection by $\ell^{A}_{ij}$. We also adopt
a notation for a Lorentz holonomy with an arbitrary number of edges. In
particular, for the holonomy comprised of ${(n-1)}$ connections between the
vertices ${\sigma_{i_{1}},\dots,\sigma_{i_{n}}}$, we write

$\displaystyle\begin{split}\Lambda^{A}_{i_{1}\cdots
i_{n}B}&=\big{(}\Lambda_{i_{1}i_{2}}\Lambda_{i_{2}i_{3}}\cdots\Lambda_{i_{n-1}i_{n}}\big{)}^{A}_{~{}B}\\\
&=\Lambda^{A}_{i_{1}i_{2}C}\Lambda^{C}_{i_{2}i_{3}D}\cdots\Lambda^{E}_{i_{n-1}i_{n}B}\end{split}$
(18)

where intermediate Lorentz indices ${\\{C,D,\dots,E\\}}$ are all contracted.
Here, we have implicitly defined the holonomy to act from the right, and note
that the matrix multiplication of holonomies effects the concatenation of
path-ordered integrals, as defined in Eq.(17).

$U_{ij}$ is seen to ‘mediate’ between Lorentz gauges at $\sigma_{i}$ and
$\sigma_{j}$, as desired. In particular, given an arbitrary Lorentz gauge
transformation defined at each vertex, say

$\displaystyle\begin{split}\Big{\\{}g_{i}=g(\sigma_{i})\in
SO(3,1)~{}\forall~{}\sigma_{i}\in\Sigma^{0}\Big{\\}},\end{split}$ (19)

the gauge transformation of $U_{ij}$ readily follows from Eq. (17), such that

$\displaystyle\begin{split}U_{ij}^{\prime}&=\left[\begin{matrix}g_{i}&\mathbf{0}\vspace{2pt}\\\
\mathbf{0}&1\end{matrix}\right]^{-1}U_{ij}\left[\begin{matrix}g_{j}&\mathbf{0}\vspace{2pt}\\\
\mathbf{0}&1\end{matrix}\right]=\left[\begin{matrix}[c|c]g_{i}^{-1}\Lambda_{ij}g_{j}&g_{i}^{-1}\ell_{ij}\\\
\hline\cr\mathbf{0}&1\end{matrix}\right].\end{split}$ (20)

From this calculation, we note that ${\ell_{ij}}$ can be regarded as if ‘based
at’ $\sigma_{i}$. By examining Eq. (17), the holonomy $U_{ji}$ is also readily
calculated to be ${U_{ji}=U_{ij}^{-1}}$. In particular,

$\displaystyle\begin{split}\Lambda_{ji}=(\Lambda_{ij})^{-1}~{}~{}\text{and}~{}~{}\ell_{ji}=-\Lambda_{ji}\ell_{ij}.\end{split}$
(21)

We must consider holonomies on closed paths as well. For example, given a loop
${(\partial\sigma^{2})_{i}=\sigma_{i}\sigma_{j}\cdots\sigma_{k}\sigma_{i}}$
around a single face $\sigma^{2}$ with basepoint $\sigma_{i}$, we define

$\displaystyle\begin{split}\left[\begin{matrix}\Omega_{ijk}&\Theta_{ijk}\\\
0&1\end{matrix}\right]=\mathcal{P}\left\\{\exp\oint_{(\partial\sigma^{2})_{i}}A\right\\}.\end{split}$
(22)

Note, when a holonomy is comprised of the connections along the edges of a
single face ${\sigma^{2}\in\Sigma^{2}}$ (i.e. when it is a “minimal”
nontrivial loop), we use the symbol $\Omega$ for its Lorentz holonomy rather
than $\Lambda$ as in Eq. (18), and we suppress some of its indices. This
notation is general to the simplicial and cubical setting, such that, for
example,

Simplicial: $\displaystyle\Omega^{AB}_{ijk}=$
$\displaystyle\big{(}\Lambda_{ij}\Lambda_{jk}\Lambda_{ki}\big{)}^{A}_{~{}C}\eta^{CB}$
Cubical: $\displaystyle\Omega^{AB}_{ijk}=$
$\displaystyle\big{(}\Lambda_{ij}\Lambda_{ji^{\prime}}\Lambda_{i^{\prime}k}\Lambda_{ki}\big{)}^{A}_{~{}C}\eta^{CB}.$
(23)

Here, $i^{\prime}$ labels the vertex diagonal to $i$ on the appropriate face
of a cubical lattice; in a more typical notation,
${(i,j,i^{\prime},k)=(\mathbf{n},\mathbf{n}+\hat{a},\mathbf{n}+\hat{a}+\hat{b},\mathbf{n}+\hat{b})}$.
$\Omega$ thereby characterizes Lorentz curvature over a face $\sigma^{2}$,
while $\Theta$ characterizes the corresponding torsion.

To see how these holonomies can be substituted for the fields of Eq. (16), let
us examine their continuous limit. We Taylor expand around $\sigma_{i}$ to
find [41]

$\displaystyle\begin{split}\Lambda_{ij}&\approx\mathbb{1}+\omega_{ij}(\sigma_{i})\Delta+\omega_{ij}(\sigma_{i})^{2}\frac{\Delta^{2}}{2}+\mathcal{O}(\Delta^{3})\\\
\ell_{ij}&\approx
e_{ij}(\sigma_{i})\Delta+\omega_{ij}(\sigma_{i})e_{ij}(\sigma_{i})\frac{\Delta^{2}}{2}+\mathcal{O}(\Delta^{3})\\\
\Omega_{ijk}-\Omega_{ikj}&\approx
2A_{f}R_{ijk}(\sigma_{i})\Delta^{2}+\mathcal{O}(\Delta^{3}).\end{split}$ (24)

Here, ${\omega_{ij}(\sigma_{i})}$ denotes the component of the continuum
Lorentz connection along the lattice edge ${\sigma_{ij}}$, evaluated at
$\sigma_{i}$. ${e_{ij}(\sigma_{i})}$ is defined analogously.
$R_{ijk}(\sigma_{i})$ denotes the component of the continuum Lorentz curvature
at $\sigma_{i}$ corresponding to edge vectors ${\sigma_{ij}}$ and
${\sigma_{ik}}$. $\Delta$ denotes the length of $\sigma_{ij}$ and
$\sigma_{ik}$ in coordinate space (in this expansion we assume these to be
equal for simplicity, though they need not be in general), and $\omega$, $e$,
and $R$ are implicitly expressed in the corresponding coordinate basis. The
area factor $A_{f}$ is analogous to $V_{f}$ in Eq. (16)—it corrects for the
implicit overcounting of area in the simplicial setting on the parallelogram
formed by $\sigma_{ij}$ and $\sigma_{ik}$. In particular, ${A_{f}=1}$ on a
cubical lattice and ${A_{f}=1/2}$ on a simplicial lattice. It is further worth
noting that the difference ${\Omega-\Omega^{-1}}$ in Eq. (24) is, in fact,
${\mathfrak{so}(3,1)}$-valued, since
${(\Omega-\Omega^{-1})^{T}\eta=\eta(\Omega^{-1}-\Omega)}$ $\forall$
${\Omega\in SO(3,1)}$.

With the expansions of Eq. (24) in mind, it is now straightforward to
reconstruct a discrete tetradic Palatini action using edge holonomies, such
that Eq. (16) is recovered to least order in the continuum limit. In
particular, we define the following action summed over lattice hypercells
${\\{\sigma^{4}\\}=\Sigma^{4}}$:

$\displaystyle\begin{split}S&=\smashoperator[]{\sum_{\sigma^{4}\in\Sigma^{4}}^{}}L(\sigma^{4})\\\
L(\sigma^{4})&=\smashoperator[]{\sum_{\begin{subarray}{c}\sigma_{i}\subset\sigma^{4}\\\
\pi\in\Pi_{i}(\sigma^{4})\end{subarray}}^{}}\frac{(-1)^{\left|\pi\right|}}{2\rho_{f}n_{v}}\epsilon_{ABCD}\left(\ell^{A}_{i\pi(1)}\ell^{B}_{i\pi(2)}\Omega^{CD}_{i\pi(3)\pi(4)}\right).\end{split}$
(25)

With factors $V_{f}$ and $A_{f}$ as defined above, the quantity
${\rho_{f}=V_{f}/A_{f}}$ satisfies ${\rho_{f}=1}$ (${\rho_{f}=12}$) for
cubical (simplicial) lattices. Note, we need not explicitly antisymmetrize
$\Omega$ and $\Omega^{-1}$ because the Levi-Civita symbol does this for us.
This cancels the factor of $2$ appearing in Eq. (24). The sum over
permutations $\pi$ replaces $\epsilon^{\mu\nu\sigma\tau}$ in Eq. (16) and the
corresponding summation of spacetime indices. In particular, as first
introduced in Eq. (16), ${\pi\in\Pi_{i}(\sigma^{4})}$ is now explicitly
defined as a permutation of vertices neighboring $\sigma_{i}$ in $\sigma^{4}$.
We define the parity $|\pi|$ to correct for any disagreement between the
overall orientation of coordinate space and the orientation of edge vectors in
$\sigma^{4}$, emanating from $\sigma_{i}$ and ordered by $\pi$.

It is worth emphasizing the following important features of this discrete
action:

  * •

The discrete Lagrangian $L(\sigma^{4})$ is locally Lorentz invariant. Under an
arbitrary gauge transformation ${\\{g_{i}\in
SO(3,1)\\}_{\sigma_{i}\in\Sigma^{0}}}$ using Eq. (20) we find

$\displaystyle\begin{split}\Big{(}&\epsilon_{ABCD}\ell^{A}_{i\pi(1)}\ell^{B}_{i\pi(2)}\Omega^{CD}_{i\pi(3)\pi(4)}\Big{)}^{\prime}\\\
&=\epsilon_{ABCD}\big{(}g_{i}^{-1}\ell_{i\pi(1)}\big{)}^{A}\big{(}g_{i}^{-1}\ell_{i\pi(2)}\big{)}^{B}\big{(}g_{i}^{-1}\Omega_{i\pi(3)\pi(4)}g_{i}\big{)}^{CD}\\\
&=\epsilon_{ABCD}\ell^{A}_{i\pi(1)}\ell^{B}_{i\pi(2)}\Omega^{CD}_{i\pi(3)\pi(4)}.\end{split}$
(26)

The last equality above follows from the Lorentz group relation
${(g_{i})^{E}_{~{}F}\eta^{FD}=(g_{i}^{-1})^{D}_{~{}F}\eta^{FE}}$ and the
$SO(3,1)$-invariance of the Levi-Civita symbol.

  * •

The Poincaré holonomies ${U_{ij}=(\Lambda_{ij},\ell_{ij})}$ are not to be
confused with the Poincaré symmetry group of Minkowski spacetime. There is a
gauge symmetry transformation that acts on our Poincaré holonomies, but the
gauge group is Lorentz, not Poincaré. Such ‘internal’ or ‘vertical’ (e.g.
Lorentz) gauge groups are typical in Cartan geometries, despite their
connections’ ‘external’ or ‘horizontal’ (e.g. translation) components [45].
Even as the internal Lorentz gauge symmetry of our theory transforms the
tetrad, it leaves spacetime geometry (i.e. the metric) completely
unaffected—regardless of what (global) symmetries the geometry may or may not
possess. By contrast, the Poincaré symmetry group of Minkowski spacetime is
comprised of global transformations of the spacetime. It is a subgroup of the
full diffeomorphism group that leaves the Minkwoski metric invariant. (I.e.,
the metric is invariant only if it happens to be Minkowski). This distinction
means, in particular, that our use of Poincaré holonomies should not be taken
to imply that we are describing Minkowski spacetime. Indeed, our theory is
capable of describing any (discrete) spacetime.

##  IV The Discrete Equations of Motion

We now compute equations of motion (EOM) by varying the discrete action with
respect to the connection. To compactify notation, when an element of the
permutation $\pi$ appears in an index, it will hereafter be denoted only by a
corresponding underlined number, for example, ${\underline{1}=\pi(1)}$. As
usual in a first-order formalism, we assume $\Lambda_{ij}$ and $\ell_{ij}$ to
be independent. Varying the action with respect to $\ell^{A}_{ij}$, and
applying the expression for $U_{ji}$ from Eq. (21) where appropriate, we find

$\displaystyle\begin{split}0=\frac{\partial
S}{\ell^{A}_{ij}}=&\sum\limits_{\sigma^{4}\supset\sigma_{ij}}\Biggr{[}\sum\limits_{\begin{subarray}{c}\pi\in\Pi_{i}(\sigma^{4})\\\
\pi(1)=j\end{subarray}}\frac{(-1)^{\left|\pi\right|}}{\rho_{f}n_{v}}\epsilon_{ABCD}\ell^{B}_{i\underline{2}}\Omega^{CD}_{i\underline{3}\underline{4}}\\\
&-\sum\limits_{\begin{subarray}{c}\pi\in\Pi_{j}(\sigma^{4})\\\
\pi(1)=i\end{subarray}}\frac{(-1)^{\left|\pi\right|}}{\rho_{f}n_{v}}\epsilon_{EBCD}\Lambda^{E}_{jiA}\ell^{B}_{j\underline{2}}\Omega^{CD}_{j\underline{3}\underline{4}}\Biggr{]}.\end{split}$
(27)

The first sum of Eq. (27) arises from terms with basepoint $i$ and the second
from terms with basepoint ${j}$. We note that although frames are permuted at
distinct basepoints in these two lines, their parities are understood to be
induced by a global orientation and are therefore mutually consistent. Eq.
(27) is counterpart to $(\delta e)$ of Eq. (12), and constitutes a discrete
reformulation of Einstein’s vacuum equations.

In particular, we may examine the continuum limit of Eq. (27) on a cubical
lattice coordinatized by $\\{x^{\mu}\\}$ with regular lattice spacing
$\Delta$. Expanding each degree of freedom near $\sigma_{i}$ as in Eq. (24)
and taking edge vectors along coordinate directions (e.g.
${\sigma_{ij}\parallel\partial_{\mu}}$), we find at leading order
$\mathcal{O}(\Delta^{3})$,

$\displaystyle
0=\epsilon^{\mu\nu\sigma\tau}\epsilon_{ABCD}e^{B}_{\nu}R^{CD}_{~{}~{}\sigma\tau}$

—mirroring ($\delta e$) of Eq. (12).

We now derive the the Lorentz connection EOM, exercising caution to ensure
that the variation of $\Lambda_{ij}$ is constrained to the $SO(3,1)$ manifold.
In particular, ${\Lambda^{T}\eta\Lambda=\eta}$ implies
${(\Lambda^{-1}\delta\Lambda)^{T}\eta+\eta(\Lambda^{-1}\delta\Lambda)=0}$, so
that ${\Lambda^{-1}\delta\Lambda\in\mathfrak{so}(3,1)}$ for a variation
$\delta\Lambda$. We can impose this constraint by taking a variation that
satisfies
${(\Lambda^{-1}\delta\Lambda)^{AB}=(\Lambda^{-1}\delta\Lambda)^{[AB]}}$, but
is otherwise arbitrary.

To that end, we consider as an example the variation of
${\Omega^{AB}_{i\underline{3}\underline{4}}=\big{(}\Lambda_{i\underline{3}}\Lambda_{\underline{3}\underline{4}}\Lambda_{\underline{4}i}\big{)}^{AB}}$
in the simplicial setting with respect to the Lorentz connection on the edge
$\sigma_{i\underline{3}}$:

$\displaystyle\delta\Omega^{AB}_{i\underline{3}\underline{4}}$
$\displaystyle=\left(\Lambda^{A}_{i\underline{3}C}\Lambda^{C}_{\underline{3}iD}\right)\delta\Lambda^{D}_{i\underline{3}E}\Lambda^{E}_{\underline{3}\underline{4}F}\Lambda^{FB}_{\underline{4}i}$
$\displaystyle=\Lambda^{A}_{i\underline{3}C}\Lambda^{B}_{i\underline{4}\underline{3}E}\left(\Lambda^{[C|}_{\underline{3}iD}\delta\Lambda^{D|E]}_{i\underline{3}}\right)$
$\displaystyle=\Lambda^{A}_{i\underline{3}[C|}\Omega^{B}_{i\underline{43}F}\Lambda^{F}_{i\underline{3}|E]}\left(\Lambda^{[C|}_{\underline{3}iD}\delta\Lambda^{D|E]}_{i\underline{3}}\right).$

The term in parentheses on the first line is a conveniently chosen form of
$\delta^{A}_{D}$—the Kronecker delta. The second line follows from the
notation of Eq. (18), the identity ${(\Lambda^{-1})^{AB}=\Lambda^{BA}}$, and
from asserting the antisymmetry of the variation
${(\Lambda^{-1}\delta\Lambda)^{[CE]}}$. The third line follows after inserting
another Kronecker delta to form a closed loop holonomy. (In general, when
$\Omega$ based at $\sigma_{i}$ is varied with respect to its Lorentz
connection along $\sigma_{jk}$, the result can be expressed in terms of
$\Omega$ or $\Omega^{-1}$ along with two antisymmetrized Lorentz
transformations that effect a parallel transport from $\sigma_{i}$ to
$\sigma_{k}$.)

To further facilitate the variation of the action, we introduce a couple
concise notations. For brevity, we denote

$\displaystyle
a^{\left|\pi\right|}=\frac{(-1)^{\left|\pi\right|}}{\rho_{f}n_{v}}$

and also define

$\displaystyle\mathcal{P}^{k\underline{12}}_{CD}=\epsilon_{ABCD}\ell_{k\underline{1}}^{A}\ell_{k\underline{2}}^{B}.$

$\mathcal{P}$ is antisymmetric both in its Lorentz indices and vertex
permutation indices. It can be roughly regarded as a (non-idempotent)
projection that annihilates any
$\ell\in\bm{\mathsfit{span}}\\{\ell_{k\underline{1}},\ell_{k\underline{2}}\\}$.

Continuing in this way, we vary $S$ with respect to
${(\Lambda_{ji}\delta\Lambda_{ij})^{[MN]}}$ to find, in the simplicial case:

$\displaystyle\begin{split}0&=\frac{\partial
S_{\text{simplicial}}}{(\Lambda_{ji}\delta\Lambda_{ij})^{[MN]}}\\\
&=\sum\limits_{\sigma^{4}\supset{\sigma_{ij}}}\Biggr{[}\sum\limits_{\begin{subarray}{c}\pi\in\Pi_{i}(\sigma^{4})\\\
\pi(3)=j\end{subarray}}a^{\left|\pi\right|}\Big{(}\mathcal{P}^{i\underline{12}}_{CD}\Omega^{D}_{i\underline{4}jE}\Big{)}\Lambda^{E}_{ij[M|}\Lambda^{C}_{ij|N]}\\\
+&\sum\limits_{\begin{subarray}{c}k\in\sigma^{4}\\\ k\neq
i,j\end{subarray}}\sum\limits_{\begin{subarray}{c}\pi\in\Pi_{k}(\sigma^{4})\\\
\pi(3)=i\\\
\pi(4)=j\end{subarray}}a^{\left|\pi\right|}\Big{(}\mathcal{P}^{k\underline{12}}_{CD}\Omega^{D}_{kijE}\Big{)}\Lambda^{E}_{kj[M|}\Lambda^{C}_{kj|N]}\\\
+&\hskip 5.0pt\sum\limits_{\begin{subarray}{c}\pi\in\Pi_{j}(\sigma^{4})\\\
\pi(4)=i\end{subarray}}a^{\left|\pi\right|}\Big{(}\mathcal{P}^{j\underline{12}}_{CD}\Omega^{D}_{j\underline{3}iE}\Big{)}\delta^{E}_{[M|}\delta^{C}_{|N]}\Biggr{]}.\end{split}$
(28)

The first line of Eq. (28) arises from terms with basepoint $i$, the middle
line from terms with basepoint ${k\neq i,j}$ in $\sigma^{4}$ and the last from
terms with basepoint ${j}$. The Lorentz EOM for a cubical discretization
follows similarly:

$\displaystyle\begin{split}0&=\frac{\partial
S_{\text{cubic}}}{(\Lambda_{ji}\delta\Lambda_{ij})^{[MN]}}\\\
&=\sum\limits_{\sigma^{4}\supset{\sigma_{ij}}}\Biggr{[}\sum\limits_{\begin{subarray}{c}\pi\in\Pi_{i}(\sigma^{4})\\\
\pi(3)=j\end{subarray}}a^{|\pi|}\Big{(}\mathcal{P}^{i\underline{12}}_{CD}\Omega^{D}_{i\underline{4}jE}\Big{)}\Lambda^{E}_{ij[M|}\Lambda^{C}_{ij|N]}\\\
&+\sum\limits_{\begin{subarray}{c}k\in\sigma^{4}\\\ k\neq
i,j\end{subarray}}\sum\limits_{\begin{subarray}{c}\pi\in\Pi_{k}(\sigma^{4})\\\
\pi(3)=i\\\
k^{\prime}=j\end{subarray}}a^{|\pi|}\Big{(}\mathcal{P}^{k\underline{12}}_{CD}\Omega^{D}_{k\underline{4}iE}\Big{)}\Lambda^{E}_{kij[M|}\Lambda^{C}_{kij|N]}\\\
&+\sum\limits_{\begin{subarray}{c}k\in\sigma^{4}\\\ k\neq
i,j\end{subarray}}\sum\limits_{\begin{subarray}{c}\pi\in\Pi_{k}(\sigma^{4})\\\
\pi(4)=j\\\
k^{\prime}=i\end{subarray}}a^{|\pi|}\Big{(}\mathcal{P}^{k\underline{12}}_{CD}\Omega^{D}_{k\underline{3}jE}\Big{)}\Lambda^{E}_{kj[M|}\Lambda^{C}_{kj|N]}\\\
&+\sum\limits_{\begin{subarray}{c}\pi\in\Pi_{j}(\sigma^{4})\\\
\pi(4)=i\end{subarray}}a^{\left|\pi\right|}\Big{(}\mathcal{P}^{j\underline{12}}_{CD}\Omega^{D}_{j\underline{3}iE}\Big{)}\delta^{E}_{[M|}\delta^{C}_{|N]}\Biggr{]}.\end{split}$
(29)

Eqs. (28-29) enforce discrete zero-torsion conditions analogous to
$(\delta\omega)$ of Eq. (12). Again employing the ordering of Eq. (24) on a
cubical lattice, it is readily computed that the least nontrivial contribution
to Eq. (29) is ${\mathcal{O}(\Delta^{3})}$ and arises from its second and
third lines alone, with ${\Omega^{D}_{~{}E}=\delta^{D}_{E}}$. This leading
order expression is given by

$\displaystyle 0$
$\displaystyle=\epsilon^{\mu\nu\alpha\beta}\epsilon_{ABNM}\Big{[}\partial_{\nu}\left(e^{A}_{\alpha}e^{B}_{\beta}\right)+\omega^{A}_{\nu
I}e^{I}_{\alpha}e^{B}_{\beta}+\omega^{B}_{\nu
J}e^{A}_{\alpha}e^{J}_{\beta}\Big{]},$

mirroring ($\delta\omega$) of Eq. (12).

Eqs. (27-29) define the desired algorithm for vacuum numerical relativity.
However, while these equations suffice to compute simulation steps in the
_bulk_ , the evolution of _boundary connections_ —including connections along
both spacelike and timelike boundaries—still requires some explanation. In
particular, even if initial and boundary connections are known a priori, Eqs.
(27-29) involve data from holonomies that generally extend outside of the
boundary wall, and are therefore underspecified on the boundary.

The strategy we adopt [54] to derive equations of motion for boundary
connections, therefore, is to extend all spacelike and timelike boundary
surfaces outward from the bulk, creating a narrow ‘double wall’ of some
fiducial thickness $\epsilon$ around the simulation domain. This double wall
is then populated with cells of width $\epsilon$, such that connections
between an inner wall vertex $\sigma_{i_{\text{in}}}$ and an outer wall vertex
$\sigma_{i_{\text{out}}}$ will have
${\Lambda_{i_{\text{in}}i_{\text{out}}}\sim\mathbb{1}+\mathcal{O}(\epsilon)}$
and ${\ell_{i_{\text{in}}i_{\text{out}}}\sim\mathcal{O}(\epsilon)}$. The
connections lying along the outer wall itself are chosen to copy the initial
or boundary conditions of the inner wall. Then, equations of motion for the
inner wall connections can be derived as usual from Eqs. (27-29), as they now
behave as connections in the bulk. Finally, we take ${\epsilon\rightarrow 0}$
in the resulting equations of motion for the (inner wall) boundary
connections.

It is worth noting that not all boundary and initial conditions will satisfy
the discrete equations of motion. Just as boundary constraints must be
satisfied in the continuum theory, care must be taken to ensure that Eqs.
(27-29) are satisfied on the initial surfaces of the discrete theory.

##  V Symplectic Structure of the Discrete Action

Variational integrators for field theories have a natural multisymplectic
structure (see e.g. Refs. [55, 56] and references therein), generalizing the
ordinary symplectic structure possessed by variational integrators in particle
mechanics [57]. Here we review the proof that variational integrators are
naturally (multi)symplectic, thereby confirming the multisymplectic structure
of Eqs. (27-29).

In a variational integrator for particle mechanics, the action evaluated on a
temporal cell ${[t_{i},t_{i+1}]}$ provides a generating function for a
canonical (symplectic) transformation across the cell [57]. Specifically, the
discrete action $S_{i}$ for cell $i$ is a generating function for the
canonical transformation $(q_{i},p_{i})\to(q_{i+1},p_{i+1})$, where $q_{i}$,
$q_{i+1}$ are the particle coordinates at the left and right end of the cell
respectively, $p_{i}=-\frac{\partial S_{i}}{\partial q_{i}}$, and
$p_{i+1}=\frac{\partial S_{i}}{\partial q_{i+1}}$. The equations of motion
(e.g. $\frac{\partial S_{i-1}}{\partial q_{i}}+\frac{\partial S_{i}}{\partial
q_{i}}=0$) guarantee that the momentum at a point is identical whether using
the left or right cell to define it. In this way the symplectic
transformations inside the cells are glued consistently across cells to
produce a global symplectic evolution.

In field theory, the situation is slightly different. A field $\phi$ has a
multimomentum $\pi^{\mu}$ (one for each dimension of space-time) [58], which
in the case of a scalar field can be recast as a 3-form
$\pi=\pi^{\mu}\mathrm{d}^{3}x_{\mu}=\frac{1}{3!}\epsilon_{\mu\alpha\beta\gamma}\pi^{\mu}\,\mathrm{d}x^{\alpha}\wedge\mathrm{d}x^{\beta}\wedge\mathrm{d}x^{\gamma}$.
For any spacetime region $R$, we then have the boundary fields
$\phi(\sigma),\pi(\sigma)$ living on $\Sigma=\partial R$, where
$\pi(\sigma)=\pi\big{|}_{\Sigma}$ can be regarded as a pseudo-scalar field.
The action evaluated over $R$, $S(R)$, is a generating function for a
submanifold $\left\\{\left(\phi(\sigma),\pi(\sigma)=\frac{\delta
S(R)}{\delta\phi}\right)\right\\}$ in this “boundary phase space”. One may
regard $S(R)$ (imprecisely) as a generating function for a canonical
transformation between any two parts of the boundary. In the case when the
boundary of $R$ consists of two disconnected pieces corresponding to two
different times, $S(R)$ is the generating function for a canonical
transformation between those times.

In the discrete setting, we take $R=\sigma^{d}$, a hypercell of maximum
dimension in our lattice ($d$ is the spacetime dimension). The boundary phase
space no longer consists of fields, but of pairs
${\left(\phi_{i},\pi_{i}(\sigma^{d})\right)}$ for each vertex
$\sigma_{i}\subset\sigma^{d}$. The discrete action over $\sigma^{d}$,
$L(\sigma^{d})$, is a generating function for a manifold in the boundary phase
space: $\left\\{\left(\phi_{i},\pi_{i}(\sigma^{d})=\frac{\partial
L(\sigma^{d})}{\partial\phi_{i}}\right)\right\\}$, in agreement with the
continuum multisymplectic structure discussed above. Note that in this case,
the momentum at a vertex is not unique, but rather depends on the hypercell
$\sigma^{d}$ used to compute it (the same holds true in the continuum: the
momentum depends on both the location and the boundary used to define it). The
equations of motion (e.g. $\sum_{\sigma^{d}\supset\sigma_{i}}\frac{\partial
L(\sigma^{d})}{\partial\phi_{i}}=0$) do not guarantee a unique momentum at
each vertex, but rather that the sum of momenta defined for each
region/boundary containing that vertex vanishes. This guarantees the
integrator will be symplectic when stepping in time.

To see this, let vertex $\sigma_{i}$ be associated with time $t_{0}$ and
define $\sigma^{d}_{+}=\\{\sigma^{d}\supset\sigma_{i}\text{ with }t>t_{0}\\}$
and $\sigma^{d}_{-}=\\{\sigma^{d}\supset\sigma_{i}\text{ with }t<t_{0}\\}$.
Then if $\pi_{i}^{+}=-\sum_{\sigma^{d}\in\sigma^{d}_{+}}\pi_{i}(\sigma^{d})$
(the minus sign takes care of the orientation for convenience) and
$\pi_{i}^{-}=\sum_{\sigma^{d}\in\sigma^{d}_{-}}\pi_{i}(\sigma^{d})$, we get
the usual gluing of symplectic transformations under time-stepping:
$\pi_{i}^{-}=\pi_{i}^{+}$. Furthermore, this will hold true no matter how we
choose to define our time and associated time-stepping (i.e. if there are
multiple ways to perform time-stepping in our cellular complex, each of them
will be guaranteed to result in symplectic evolution). This argument neglects
subtleties that may arise at boundaries or when the number of vertices changes
between time slices. To resolve these, a more global perspective is necessary,
following for example the presentation in [59, 60, 61].

In the case of gravity in four dimensions, we are using 1-form fields rather
than scalar fields, so the multimomentum is more naturally a 2-form.
Additionally, our two primary fields ($e$ and $\omega$) are conjugate to each
other (in the sense that the multimomentum of $\omega$ is a function of $e$,
while the multimomentum of $e$ vanishes, leaving behind an ${(e,\omega)}$
phase space). All of this is captured by the discrete Eqs. (27-29). The
boundary phase space of a cell $\sigma^{4}$ consists of pairs
$(\ell_{ij}^{A},\Lambda^{A}_{ijB})$ for each edge
$\sigma_{ij}\subset\sigma^{4}$ (rather than for each vertex, as in the case of
a scalar field). The discrete action $L(\sigma^{4})$ is a generating function,
which defines momenta conjugate to $\ell_{ij}^{A}$ and $\Lambda^{A}_{ijB}$ as
the bracketed summands of Eqs. (27) and Eqs. (28-29), respectively. The
momentum conjugate to $\Lambda$ is a function of $\ell$, while the
initialization of the algorithm (see the end of Sec. IV) ensures the vanishing
of the momentum conjugate to $\ell$ on the inner wall of the double wall
boundary (i.e. torsion is made to vanish by construction in the
$\epsilon$-width cells making up the double wall). The symplectic structure of
the time-stepping then ensures the momentum conjugate to $\ell$ vanishes for
the entire complex. In this way the discrete multisymplectic structure of our
gravitational integrator reproduces the continuous multisymplectic structure
of GR.

##  VI Conclusion

We have presented a new numerical scheme for general relativity, detailed in
Eqs. (27) and (28)-(29). This scheme preserves both the (multi)symplectic
structure and local Lorentz invariance of the tetrad formulation of GR.
Furthermore, its discrete variables have a clear relationship with their
continuum counterparts. As such, this scheme holds promise as an integrator
for numerical relativity (its structure preservation maintains exact
conservation laws and bounded errors in simulation) and for studying the
classical limits of certain quantum gravity theories (such as loop quantum
gravity, spin foams, etc.). In these roles, the scheme’s symplectic structure
promises an improvement over non-symplectic finite-difference and spectral
methods, while its natural association with continuum variables makes it a
more viable alternative to other symplectic approaches to discrete gravity
(most notably Regge calculus). In future work, implementations of this
algorithm will be needed to demonstrate its practical utility. Furthermore,
like Regge calculus, further study is required to incorporate matter into our
approach (though the way forward seems clearer).

It may also be of interest to explore the potential union between the
algorithm defined here and other structure-preserving discretizations suitable
for numerical relativity. For example, it may be useful to explore the
relationship between our holonomy-centric approach and the recently developed
technique of group-equivariant interpolation in symmetric spaces [62, 63]. It
may also be useful to compare our effort with finite element cochain complexes
suitable for applications in numerical relativity [64]. In this way, the
algorithm we have introduced can be an advantageous starting point for
explorations into structure-preserving discrete gravity theories.

##  VII Acknowledgments

Thank you to Hong Qin for helpful discussions and encouragement. Thank you to
Robert Littlejohn for discussions and inspiration in the lead-up to this work.
This work was performed under the auspices of the U.S. Department of Energy by
Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344. This
research was further supported by the U.S. Department of Energy Fusion Energy
Sciences Postdoctoral Research Program administered by the Oak Ridge Institute
for Science and Education (ORISE) for the DOE. ORISE is managed by Oak Ridge
Associated Universities (ORAU) under DOE contract number DE-SC0014664. All
opinions expressed in this paper are the authors’ and do not necessarily
reflect the policies and views of DOE, ORAU, or ORISE. A.S.G. further
acknowledges the generous support of the Princeton University Charlotte
Elizabeth Procter Fellowship.

## References

  * Hairer  _et al._ [2006] E. Hairer, C. Lubich, and G. Wanner,  _Geometric Numerical Integration_ , 2nd ed. (Springer, Berlin, 2006).
  * Kinoshita  _et al._ [1991] H. Kinoshita, H. Yoshida, and H. Nakai, Symplectic integrators and their application to dynamical astronomy, CELESTIAL MECHANICS AND DYNAMICAL ASTRONOMY 50, 59 (1991).
  * Gladman  _et al._ [1991] B. Gladman, M. Duncan, and J. Candy, Symplectic integrators for long-term integrations in celestial mechanics, Celestial Mechanics and Dynamical Astronomy 52, 221 (1991).
  * Chambers  _et al._ [2002] J. E. Chambers, E. V. Quintana, M. J. Duncan, and J. J. Lissauer, Symplectic Integrator Algorithms for Modeling Planetary Accretion in Binary Star Systems, The Astronomical Journal 123, 2884 (2002).
  * Bravetti  _et al._ [2020] A. Bravetti, M. Seri, M. Vermeeren, and F. Zadra, Numerical integration in celestial mechanics: a case for contact geometry, Celestial Mechanics and Dynamical Astronomy 132 (2020).
  * Li  _et al._ [2012] X. Li, W. Wang, M. Lu, M. Zhang, and Y. Li, Structure-preserving modelling of elastic waves, Geophysical Journal International 188, 1382 (2012).
  * Liu  _et al._ [2015] S. Liu, X. Li, W. Wang, L. Xu, and B. Li, A modified symplectic scheme for seismic wave modeling, Journal of Applied Geophysics 116, 110 (2015).
  * Squire  _et al._ [2012] J. Squire, H. Qin, and W. M. Tang, Geometric integration of the Vlasov-Maxwell system with a variational particle-in-cell scheme, Physics of Plasmas 19, 084501 (2012).
  * Xiao  _et al._ [2015] J. Xiao, H. Qin, J. Liu, Y. He, R. Zhang, and Y. Sun, Explicit high-order non-canonical symplectic particle-in-cell algorithms for Vlasov-Maxwell systems, Physics of Plasmas 22, 112504 (2015).
  * He  _et al._ [2015] Y. He, H. Qin, Y. Sun, J. Xiao, R. Zhang, and J. Liu, Hamiltonian time integrators for Vlasov-Maxwell equations, Physics of Plasmas 22, 124503 (2015).
  * Crouseilles  _et al._ [2015] N. Crouseilles, L. Einkemmer, and E. Faou, Hamiltonian splitting for the Vlasov–Maxwell equations, Journal of Computational Physics 283, 224 (2015).
  * Qin  _et al._ [2016] H. Qin, J. Liu, J. Xiao, R. Zhang, Y. He, Y. Wang, Y. Sun, J. W. Burby, L. Ellison, and Y. Zhou, Canonical symplectic particle-in-cell method for long-term large-scale simulations of the Vlasov-Maxwell equations, Nuclear Fusion 56, 014001 (2016).
  * Kraus  _et al._ [2017] M. Kraus, K. Kormann, P. J. Morrison, and E. Sonnendrücker, GEMPIC: Geometric ElectroMagnetic Particle-In-Cell Methods, Journal of Plasma Physics 83 (2017).
  * Morrison [2017] P. J. Morrison, Structure and structure-preserving algorithms for plasma physics, Physics of Plasmas 24, 055502 (2017).
  * Glasser and Qin [2020] A. S. Glasser and H. Qin, The geometric theory of charge conservation in particle-in-cell simulations, Journal of Plasma Physics 86, 835860303 (2020).
  * Glasser and Qin [2021] A. S. Glasser and H. Qin, A gauge-compatible Hamiltonian splitting algorithm for particle-in-cell simulations using finite element exterior calculus, arXiv:2110.10346 [physics.plasm-ph] (2021).
  * Regge [1961] T. Regge, General relativity without coordinates, Il Nuovo Cimento 19 (1961).
  * Misner  _et al._ [1973] C. W. Misner, K. S. Thorne, and J. A. Wheeler,  _Gravitation_ (W. H. Freeman, 1973).
  * Hawking [1978] S. Hawking, Spacetime foam, Nuclear Physics B 144, 349 (1978).
  * Caselle  _et al._ [1989] M. Caselle, A. D’Adda, and L. Magnea, Regge calculus as a local theory of the Poincaré group, Physics Letters B 232 (1989).
  * Williams and Tuckey [1992] R. M. Williams and P. A. Tuckey, Regge calculus: a brief review and bibliography, Classical and Quantum Gravity 9, 1409 (1992).
  * Immirzi [1997] G. Immirzi, Quantum Gravity and Regge Calculus, Nuclear Physics B - Proceedings Supplements 57, 65 (1997).
  * Loll [1998] R. Loll, Discrete Approaches to Quantum Gravity in Four Dimensions, Living Reviews in Relativity 1 (1998).
  * Ambjörn  _et al._ [2000] J. Ambjörn, J. Jurkiewicz, and R. Loll, Nonperturbative Lorentzian Path Integral for Gravity, Physical Review Letters 85, 924 (2000).
  * Gionti [2005] G. Gionti, Discrete gravity as a local theory of the Poincaré group in the first-order formalism, Classical and Quantum Gravity 22, 4217 (2005).
  * Dittrich  _et al._ [2021] B. Dittrich, S. Gielen, and S. Schander, Lorentzian quantum cosmology goes simplicial, Classical and Quantum Gravity (2021).
  * Collins and Williams [1972] P. A. Collins and R. M. Williams, Application of Regge Calculus to the Axially Symmetric Initial-Value Problem in General Relativity, Physical Review D 5, 1908 (1972).
  * Collins and Williams [1973] P. A. Collins and R. M. Williams, Dynamics of the Friedmann Universe Using Regge Calculus, Physical Review D 7, 965 (1973).
  * Sorkin [1975] R. Sorkin, Time-evolution problem in regge calculus, Physical Review D 12, 385 (1975).
  * Porter [1987] J. Porter, A new approach to the Regge calculus, Classical and Quantum Gravity 4, 375 (1987).
  * Dubal [1989] M. R. Dubal, Relativistic collapse using Regge calculus. I. Spherical collapse equations, Classical and Quantum Gravity 6, 1925 (1989).
  * Barrett  _et al._ [1997] J. W. Barrett, M. Galassi, W. A. Miller, R. D. Sorkin, P. A. Tuckey, and R. M. Williams, A Parallelizable Implicit Evolution Scheme for Regge Calculus, International Journal of Theoretical Physics 36, 815 (1997).
  * Gentle [2002] A. P. Gentle, Regge Calculus: A Unique Tool for Numerical Relativity, General Relativity and Gravitation 34, 1701 (2002).
  * Khavari [2009] P. Khavari,  _Regge Calculus as a Numerical Approach to General Relativity_ , Ph.D. thesis, University of Toronto, Toronto, ON (2009).
  * Gentle [2013] A. P. Gentle, A cosmological solution of Regge calculus, Classical and Quantum Gravity 30, 085004 (2013).
  * Gentle and Miller [2002] A. P. Gentle and W. A. Miller, A brief review of Regge calculus in classical numerical relativity, arXiv:gr-qc/0101028 , 1467 (2002).
  * Barrett  _et al._ [2019] J. W. Barrett, D. Oriti, and R. M. Williams, Tullio Regge’s legacy: Regge calculus and discrete gravity, arXiv:1812.06193 [gr-qc, physics:hep-th] (2019).
  * Gentle  _et al._ [1999] A. P. Gentle, D. E. Holz, W. A. Miller, and J. A. Wheeler, Apparent horizons in simplicial Brill wave initial data, Classical and Quantum Gravity 16, 1979 (1999).
  * Gentle [1999] A. P. Gentle, Simplicial Brill wave initial data, Classical and Quantum Gravity 16, 1987 (1999).
  * Dittrich and Höhn [2010] B. Dittrich and P. A. Höhn, From covariant to canonical formulations of discrete gravity, Classical and Quantum Gravity 27, 155001 (2010).
  * Kogut [1979] J. B. Kogut, An introduction to lattice gauge theory and spin systems, Reviews of Modern Physics 51, 659 (1979).
  * Palatini [1919] A. Palatini, Deduzione invariantiva delle equazioni gravitazionali dal principio di Hamilton, Rendiconti del Circolo Matematico di Palermo (1884-1940) 43, 203 (1919).
  * Menotti and Pelissetto [1987a] P. Menotti and A. Pelissetto, Poincare, de Sitter, and conformal gravity on the lattice, Physical Review D 35, 1194 (1987a).
  * Menotti and Pelissetto [1987b] P. Menotti and A. Pelissetto, Gauge invariance and functional integration measure in lattice gravity, Nuclear Physics B 288, 813 (1987b).
  * Sharpe [1997] R. Sharpe,  _Differential Geometry \- Cartan’s Generalization of Klein’s Erlangen Program_ (Springer-Verlag, New York, 1997).
  * Hehl  _et al._ [1976] F. W. Hehl, P. von der Heyde, G. D. Kerlick, and J. M. Nester, General relativity with spin and torsion: Foundations and prospects, Reviews of Modern Physics 48, 393 (1976).
  * Hehl [1980] F. W. Hehl, Four Lectures on Poincaré Gauge Field Theory, in  _Cosmology and Gravitation_ , NATO Advanced Study Institutes Series, edited by P. G. Bergmann and V. De Sabbata (Springer US, Boston, MA, 1980) pp. 5–61.
  * Menotti [2017] P. Menotti, Lectures on Gravitation (2017).
  * Nakahara [2018] M. Nakahara,  _Geometry, topology and physics_ (CRC press, 2018).
  * Haber [1970] S. Haber, Numerical evaluation of multiple integrals, SIAM review 12, 481 (1970).
  * Desbrun  _et al._ [2005] M. Desbrun, A. N. Hirani, M. Leok, and J. E. Marsden, Discrete exterior calculus, arXiv:math/0508341 (2005).
  * Arnold  _et al._ [2006] D. N. Arnold, R. S. Falk, and R. Winther, Finite element exterior calculus, homological techniques, and applications, Acta Numerica 15 (2006).
  * Arnold  _et al._ [2010] D. N. Arnold, R. S. Falk, and R. Winther, Finite element exterior calculus: from Hodge theory to numerical stability, Bulletin of the American Mathematical Society 47, 281 (2010).
  * Stern [2009] A. Stern,  _Geometric Discretization of Lagrangian Mechanics and Field Theories_ , Ph.D. thesis, California Institute of Technology, Pasadena, California (2009).
  * Marsden  _et al._ [1998] J. E. Marsden, G. W. Patrick, and S. Shkoller, Multisymplectic geometry, variational integrators, and nonlinear pdes, Communications in Mathematical Physics 199, 351 (1998).
  * Gotay and Marsden [2006] M. Gotay and J. Marsden, Momentum maps and classical relativistic fields. part iii: Gauge symmetries and initial value constraints (2006).
  * Marsden and West [2001] J. E. Marsden and M. West, Discrete mechanics and variational integrators, Acta Numerica 10, 357 (2001).
  * Gotay  _et al._ [1998] M. J. Gotay, J. Isenberg, J. E. Marsden, and R. Montgomery, Momentum maps and classical relativistic fields. part i: Covariant field theory, arXiv preprint physics/9801019 (1998).
  * Dittrich and Höhn [2012] B. Dittrich and P. A. Höhn, Canonical simplicial gravity, Classical and quantum gravity 29, 115009 (2012).
  * Dittrich and Höhn [2013] B. Dittrich and P. A. Höhn, Constraint analysis for variational discrete systems, Journal of Mathematical Physics 54, 093505 (2013).
  * Höhn [2015] P. A. Höhn, Canonical linearized regge calculus: counting lattice gravitons with pachner moves, Physical Review D 91, 124034 (2015).
  * Gawlik and Leok [2018] E. S. Gawlik and M. Leok, Interpolation on Symmetric Spaces Via the Generalized Polar Decomposition, Foundations of Computational Mathematics 18, 757 (2018).
  * Leok [2019] M. Leok, Variational Discretizations of Gauge Field Theories Using Group-Equivariant Interpolation, Foundations of Computational Mathematics 19, 965 (2019).
  * Arnold and Hu [2021] D. N. Arnold and K. Hu, Complexes from complexes, arXiv:2005.12437 [cs, math] (2021).