creamrgb.81, .88, 1

0

#  First-order hyperbolic formulation of the
teleparallel gravity theory

Ilya Peshkov∗,1 Héctor Olivares2 & Evgeniy Romenski3

## Abstract

Motivated by numerically solving the Einstein field equations, we derive a
first-order reduction of the second-order $f(T)$-teleparallel gravity field
equations in the pure-tetrad formulation (no spin connection). We then
restrict our attention to the teleparallel equivalent of general relativity
(TEGR) and propose a 3+1 decomposition of the governing equations that can be
used in a computational code. We demonstrate that for the matter-free space-
time the obtained system of first-order equations is equivalent to the tetrad
reformulation of general relativity by Estabrook, Robinson, Wahlquist, and
Buchman and Bardeen and therefore also admits a symmetric hyperbolic
formulation. The structure of the 3+1 equations resembles a lot of
similarities with the equations of relativistic electrodynamics and the
recently proposed dGREM tetrad-reformulation of general relativity.
00footnotetext: * Corresponding author<EMAIL_ADDRESS>00footnotetext:
1 Department of Civil, Environmental and Mechanical Engineering, University of
Trento, Via Mesiano 77, 38123, Trento, Italy 00footnotetext: 2 Department of
Astrophysics/IMAPP, Radboud University Nijmegen, P.O. Box 9010, NL-6500 GL
Nijmegen, Netherlands 00footnotetext: 3 Sobolev Institute of Mathematics,
Novosibirsk, Russia

##  1 Introduction

The teleparallel gravity theory is one of the alternative reformulations of
Einstein’s general relativity (GR) [29, 2, 14]. While for GR gravitational
interaction is a manifestation of the curvature of a torsion-free spacetime,
for teleparallel gravity theory it is realized as a curvature-free linear
connection with non-zero torsion. Although different variables can be taken as
the main dynamical fields in GR (tetrad fields, soldering forms, etc.), the
most extended choice is the metric tensor accompanied by a Levi-Civita
connection. In contrast, for teleparallel theory the metric is trivial and the
main dynamical field is usually the tetrad (or frame) field.

Despite the teleparallel gravity is considered as an alternative physical
theory111Yet, it is guaranteed that TEGR produces all the classical results of
general relativity [2, 5]. to Einstein’s gravity with several promising
features missing in GR, e.g. see the discussion in [2, Sec.18] and [14], in
this paper, we are interested in the teleparallel gravity only from a pure
computational viewpoint and our goal is to use its mathematical structure to
develop an efficient computational framework for numerical relativity. Thus,
the main goal of this paper is to derive a $3+1$-split of the so-called
_teleparallel equivalent of general relativity_ (TEGR) [2, 35] which is know
to pass all standard tests of GR.

From the numerical view point, a proposed $3+1$ formulation has to have well-
posed initial value problem (i.e. a solution exists, the solution is unique
and changes continuously with changes in the initial data) in order to compute
stable evolution of the numerical solution. In other words, the system of
governing equations has to be hyperbolic222Note that in the numerical
relativity the term “strong hyperbolicity” is used emphasizing that not only
eigenvalues must be real but that the full basis of eigenvectors must exist.
Therefore, the second goal of the paper is to test if the proposed $3+1$
formulation of TEGR is hyperbolic. As is the case with other first-order
reductions of the Einstein equations [6], the question of hyperbolicity of
$3+1$ TEGR equations considered here is not trivial and depends on the
delicate use of multiple involution constraints (stationary identities). In
particular, for a vacuum space-time, we have found that the proposed $3+1$
TEGR equations admits an equivalent reformulation which is in turn equivalent
to the symmetric hyperbolic tetrad formulation of GR by Estabrook-Robinson-
Wahlquist [21] and Buchman-Bardeen [11].

We note that we do not consider the most general TEGR formulation as for
example presented in [2]. The linear connection of TEGR is the sum of the two
parts: the Weitzenböck connection (which is the historical connection of TEGR)
and the spin connection (parametrized by Lorentz matrices) representing the
inertial content of the tetrad. The spin connection is necessary to separate
the inertial effect of a chosen frame from its gravity content as well as for
establishing the full covariance of the theory [2, 26, 35], i.e. with respect
to both the diffeomorphisms of the spacetime and the Lorentz transformations
of tangent spaces. However, being important from the theoretical viewpoint and
for extensions of teleparallel gravity [26], the spin connection of TEGR
introduces extra degrees of freedom which do not have evolution equations and
therefore can be treated as parameters (not state variables) of the theory.
Since we are interested in developing a computational framework for GR, the
consideration of this paper is, therefore, restricted to the frames for which
the spin connection is set to zero globally, that is we deliberately work in
the class of inertial frames. We thus consider TEGR in its historical, or
_pure tetrad_ , formulation [26].

Interestingly that another motivation to study the mathematical structure of
the teleparallel gravity is coming from the continuum fluid and solid
mechanics. In particular, the role of the torsion to describe defects in
solids has been known for decades now, i.e. see [32, 30, 50, 37, 8]. Moreover,
the material tetrad field (called also the distortion field in our papers) is
the key field for the unified hyperbolic model of fluid and solid mechanics
[42, 19]. In such a theory, the concept of torsion can be connected with the
inertial effect of small-scale eddies in turbulent flows and with dispersion
effects in heterogeneous solids (e.g. acoustic metamaterials) as discussed in
[43]. Interestingly, that the $3+1$ equations we obtain in this paper resemble
very closely the structure of the equations for continuum fluid and solid
mechanics with torsion [43]. Furthermore, the unified theory of fluids and
solids has been also extended in the general relativistic settings [46] and
therefore, as being a tetrad theory by its nature, it can be straightforwardly
coupled with the $3+1$ TEGR equations discussed in this paper.

Last but not least, as we shall see in Sec. 8, the derived $3+1$ equations of
TEGR are equivalent to the recently proposed dGREM reformulation of GR [39]
which is also a tetrad-based formulation but derived based on the exterior
calculus formalism.

##  2 Definitions

###  2.1 Anholonomic frame field

Throughout this paper, we use the following index convention. Greek letters
$\alpha,\beta\,\gamma$, $\ldots,\lambda,\mu,\nu,...=0,1,2,3$ are used to index
quantities related to the spacetime manifold, the Latin letters from the
beginning of the alphabet $a,b,c,...=\hat{0},\hat{1},\hat{2},\hat{3}$ are used
to index quantities related to the tangent Minkowski space. In the $3+1$
split, the letters $i,j,k,\ldots=1,2,3$ from the middle of the Latin alphabet
are used to denote spatial components of the spacetime tensors, and the upper
case Latin letters $A,B,C,\ldots=\hat{{1}},\hat{{2}},\hat{{3}}$ index spatial
components of the tensors written in a chosen frame of the tangent space.

We shall consider a spacetime manifold $M$ equipped with a general Riemannian
metric $g_{\mu\nu}$ and a coordinate system $x^{\mu}$. At each point of the
spacetime, there is the tangent space $T_{x}M$ spanned by the frame (or
tetrad) ${\bm{\partial}_{\mu}}$ which is the standard coordinate basis. There
is also the cotangent space $T_{x}^{*}M$ spanned by the coframe field
${\bm{\mathrm{d}x}^{\mu}}$. Recall that the frames ${\bm{\partial}_{\mu}}$ and
${\bm{\mathrm{d}x}^{\mu}}$ are _holonomic_ [2].

In addition to $T_{x}M$, we assume that at each point of $M$, there is a
soldered tangent space which is a Minkowski space spanned by an orthonormal
frame (tetrad) $\bm{\eta}_{a}$ and equipped with the Minkowski metric

$\kappa_{ab}=\text{diag}(-1,1,1,1).$ (1)

Similarly, there is the corresponding cotangent space spanned by the co-frames
$\bm{h}^{a}$. It is assumed that the frames $\bm{h}^{a}$ and $\bm{\eta}_{a}$
are independent of the coordinates $x^{\mu}$ and, therefore, are non-holonomic
in general.

The components of the non-holonomic frames $\bm{h}^{a}$ and $\bm{\eta}_{a}$ in
the holonomic frames ${\bm{\mathrm{d}x}^{\mu}}$ and ${\bm{\partial}_{\mu}}$
are denoted by $h^{a}_{\phantom{a}\mu}$ and $\eta^{\mu}_{\phantom{\mu}a}$,
i.e.

$\bm{h}^{a}=h^{a}_{\phantom{a}\mu}{\bm{\mathrm{d}x}^{\mu}},\qquad\text{or}\qquad\bm{\eta}_{a}=\eta^{\mu}_{\phantom{\mu}a}{\bm{\partial}_{\mu}}$
(2)

with $h^{a}_{\phantom{a}\mu}$ being the inverse of
$\eta^{\mu}_{\phantom{\mu}a}$, i.e.

$h^{a}_{\phantom{a}\mu}\eta^{\mu}_{\phantom{\mu}b}=\delta^{a}_{\
b},\qquad\eta^{\nu}_{\phantom{\nu}a}h^{a}_{\phantom{a}\mu}=\delta^{\nu}_{\
\mu},$ (3)

where $\delta^{a}_{\ b}$ and $\delta^{\nu}_{\ \mu}$ are the Kronecker deltas.

The soldering of the spacetime and the tangent Minkowski space means that the
metrics $g_{\mu\nu}$ and $\kappa_{ab}$ are related by

$g_{\mu\nu}=\kappa_{ab}h^{a}_{\phantom{a}\mu}h^{b}_{\phantom{b}\nu}.$ (4)

Note that

$h:=\det(h^{a}_{\phantom{a}\mu})=\sqrt{-g},$ (5)

if $g=\det(g_{\mu\nu})$.

###  2.2 Observer’s 4-velocity

Observer’s 4-velocity is associated with the $0$-th vector of the tetrad basis

$u^{\mu}:=\eta^{\mu}_{\phantom{\mu}\hat{{0}}},\qquad
u_{\mu}=g_{\mu\nu}u^{\nu}.$ (6)

Also, due to

$u_{\mu}=g_{\mu\nu}u^{\nu}=\kappa_{ab}h^{a}_{\phantom{a}\mu}h^{b}_{\phantom{b}\nu}\eta^{\nu}_{\phantom{\nu}\hat{{0}}}=\kappa_{a\hat{{0}}}h^{a}_{\phantom{a}\mu}=-h^{\hat{{0}}}_{\phantom{\hat{{0}}}\mu},$
(7)

the covariant components of the 4-velocity equal to entries of the $0$-th
vector of the co-basis with the opposite sign

$u_{\mu}=-h^{\hat{{0}}}_{\phantom{\hat{{0}}}\mu}.$ (8)

###  2.3 Connection and torsion

Because we work in the framework of the pure-tetrad formulation of TEGR, the
linear connection is set to the pure Weitzenböck connection333Note that we use
a different convention on the positioning of the lower indices of the
Weitzenböck connection
$W^{a}_{\phantom{a}\mu\nu}=\partial_{\mu}h^{a}_{\phantom{a}\nu}$ than in [2].
Namely, the derivative index goes first. [2, 34]:

$W^{a}_{\phantom{a}\mu\nu}:=\partial_{\mu}h^{a}_{\phantom{a}\nu},\qquad\text{or}\qquad
W^{\lambda}_{\phantom{\lambda}\mu\nu}:=\eta^{\lambda}_{\phantom{\lambda}a}\partial_{\mu}h^{a}_{\phantom{a}\nu}.$
(9)

The torsion is then defined as

$T^{a}_{\phantom{a}\mu\nu}:=\partial_{\mu}h^{a}_{\phantom{a}\nu}-\partial_{\nu}h^{a}_{\phantom{a}\mu}=W^{a}_{\phantom{a}\mu\nu}-W^{a}_{\phantom{a}\nu\mu}.$
(10)

Note that while the spacetime derivatives commute

$\displaystyle\partial_{\mu}(\partial_{\nu}V^{\lambda})-\partial_{\nu}(\partial_{\mu}V^{\lambda})$
$\displaystyle=0,$ (11)
$\displaystyle\partial_{\mu}(\partial_{\nu}V^{a})-\partial_{\nu}(\partial_{\mu}V^{a})$
$\displaystyle=0,$ (12)

their tangent space counterparts
$\partial_{a}=\eta^{\mu}_{\phantom{\mu}a}\partial_{\mu}$ do not (for non-
vanishing torsion)

$\partial_{b}(\partial_{c}V^{a})-\partial_{c}(\partial_{b}V^{a})=-T^{d}_{\phantom{d}bc}\partial_{d}V^{a},$
(13)

where
$T^{d}_{\phantom{d}bc}=T^{d}_{\phantom{d}\mu\nu}\eta^{\mu}_{\phantom{\mu}b}\eta^{\nu}_{\phantom{\nu}c}$.

###  2.4 Levi-Civita symbol (tensor density)

We shall also need the Levi-Civita symbol (tensor-density444We use the sign
convention for the tensor density weight according to [49, 28], i.e. under a
general coordinate change $x^{\mu}\to x^{\mu^{\prime}}$ the determinant
$\det(h^{a}_{\phantom{a}\mu})=h$ transforms as
$h^{\prime}=\det\left(\frac{\partial x^{\mu}}{\partial
x^{\mu^{\prime}}}\right)^{W}\cdot h$ with $W=+1$. Therefore, the tetrad’s
determinant $h$ and the square root of the metric determinant $h=\sqrt{-g}$
have weights $+1$, as well as the Lagrangian density in the action integral.
of weight $+1$)

$\bm{\in}^{\lambda\mu\nu\rho}=\left\\{\begin{array}[]{ll}+1,&\text{if \
}\lambda\mu\nu\rho\text{ is an even permutation of }0123,\\\\[5.69054pt]
-1,&\text{if \ }\lambda\mu\nu\rho\text{ is an odd \ permutation of
}0123,\\\\[5.69054pt] \phantom{-}0,&\text{otehrwise}.\end{array}\right.$ (14)

Its covariant components $\bm{\in}_{\lambda\mu\nu\rho}$ define a tensor
density of weight $-1$ with the reference value $\bm{\in}_{0123}=-1$. One
could define an absolute Levi-Civita contravariant
$\varepsilon^{\lambda\mu\nu\rho}=h^{-1}\bm{\in}^{\lambda\mu\nu\rho}$ and
covariant $\varepsilon_{\lambda\mu\nu\rho}=h\bm{\in}_{\lambda\mu\nu\rho}$
ordinary tensors but for our further considerations (see Sec. 5), it is
important that the derivatives $\partial_{\sigma}\bm{\in}^{\lambda\mu\nu\rho}$
vanish:

$\partial_{\sigma}\bm{\in}^{\lambda\mu\nu\rho}=0,$ (15)

whereas for $\varepsilon^{\lambda\mu\nu\rho}$ one has

$\partial_{\sigma}\varepsilon^{\lambda\mu\nu\rho}=\partial_{\sigma}(h^{-1}\bm{\in}^{\lambda\mu\nu\rho})=\bm{\in}^{\lambda\mu\nu\rho}\partial_{\sigma}h^{-1}=\\\\[5.69054pt]
-\bm{\in}^{\lambda\mu\nu\rho}h^{-1}\eta^{\eta}_{\phantom{\eta}a}\partial_{\sigma}h^{a}_{\phantom{a}\eta}=-\varepsilon^{\lambda\mu\nu\rho}W^{\eta}_{\phantom{\eta}\sigma\eta},$
(16)

which is not zero in general.

##  3 Variational formulation

We consider a general Lagrangian (scalar-density)
$\Lambda(h^{a}_{\phantom{a}\mu},\partial_{\lambda}h^{a}_{\phantom{a}\nu})$ of
the teleparallel gravity which is a function of the frame field
$h^{a}_{\phantom{a}\mu}$ and its first gradients
$W^{a}_{\phantom{a}\lambda\nu}=\partial_{\lambda}h^{a}_{\phantom{a}\nu}$. In
what follows, we shall not explicitly split $\Lambda$ into the gravity (g) and
matter (m) parts (unless it is explicitly mentioned otherwise), i.e.

$\Lambda=\Lambda^{\text{(m)}}(h^{a}_{\phantom{a}\mu},\partial_{\lambda}h^{a}_{\phantom{a}\mu},\ldots)+\Lambda^{\text{(g)}}(h^{a}_{\phantom{a}\mu},\partial_{\lambda}h^{a}_{\phantom{a}\nu})$
(17)

but the derivation will be performed for the total unspecified Lagrangian
$\Lambda$, so that in principle the derivation can be adopted for the
extensions of the teleparallel gravity such as $f(\mathcal{T})$-teleparallel
gravity, where $f(\mathcal{T})$ is some function of the torsion scalar
$\mathcal{T}$, see Section 12. We shall utilize the explicit form of the TEGR
Lagrangian only in the last part of the paper.

Varying the action of the teleparallel gravity
$\int\Lambda(h^{a}_{\phantom{a}\mu},\partial_{\lambda}h^{a}_{\phantom{a}\mu})\bm{{\bm{\mathrm{d}x}}}$
with respect to the tetrad, one obtains the Euler-Lagrange equations

$\frac{\delta\Lambda}{\delta
h^{a}_{\phantom{a}\mu}}=\partial_{\lambda}(\Lambda_{\partial_{\lambda}h^{a}_{\phantom{a}\mu}})-\Lambda_{h^{a}_{\phantom{a}\mu}}=0,$
(18)

where
$\Lambda_{\partial_{\lambda}h^{a}_{\phantom{a}\mu}}=\frac{\partial\Lambda}{\partial(\partial_{\lambda}h^{a}_{\phantom{a}\mu})}$
and $\Lambda_{h^{a}_{\phantom{a}\mu}}=\frac{\partial\Lambda}{\partial
h^{a}_{\phantom{a}\mu}}$. Equations (18) form a system of 16 _second-order_
partial differential equations for 16 unknowns $h^{a}_{\phantom{a}\mu}$. Our
goal is to replace this second-order system by an equivalent but larger system
of only first-order partial differential equations.

##  4 Equivalence to GR

Before writing system (18) as a system of first-order equations let us first
make a comment on the equivalence of the GR and TEGR formulations.

The Lagrangian density of TEGR (and its extensions) is formed from the torsion
scalar $\mathcal{T}$, see (106). As it is known, e.g. see [2, Eq.(9.30)], the
torsion scalar can be written as

$h\mathcal{T}=-\sqrt{-g}R-\partial_{\mu}(2hT^{\nu}_{\phantom{\nu}\lambda\nu}g^{\lambda\mu}),$
(19)

where $R$ is the Ricci scalar and
$T^{\nu}_{\phantom{\nu}\lambda\mu}=\eta^{\nu}_{\phantom{\nu}a}T^{a}_{\phantom{a}\lambda\mu}$.
In other words, the Lagrangians of TEGR and GR differ by the four-divergence
term (surface term). The latter does not affect the Euler-Lagrange equations
if there are no boundaries which is implied in this paper. Therefore, Euler-
Lagrange equations of TEGR (18) in vacuum is nothing else but the Euler-
Lagrange equations of GR written in terms of the tetrads and hence, their
physical solutions must be equivalent because the information about the
physical interaction is contained in the Euler-Lagrange equations of a theory.
What is different in GR and TEGR is the way one interprets the tetrads and
their first derivatives, i.e. the way one defines the linear connection of the
spacetime from the gradients of tetrads, e.g. torsion-free Levi-Civita
connection of GR and curvature-free Weitzenböck connection of TEGR. These
different geometrical interpretations then define extra evolution equations
(_compatibility constraints/identities_ , e.g. see (23)) that are merely
consequences of the geometrical definitions but do not define the physics of
the gravitational interaction. The critical point for the numerical
relativity, though, is that these extra evolution equations must be solved
simultaneously with the Euler-Lagrange equations and may affect the
mathematical regularity (well-posedness of the Cauchy problem) of the
resulting system.

##  5 First-order extension

Our first goal is to replace second-order system (18) by a larger but first-
order system. This is achieved in this section.

From now on, we shall treat the frame field $h^{a}_{\phantom{a}\mu}$ and its
gradients (the Weitzenböck connection)
$\partial_{\lambda}h^{a}_{\phantom{a}\mu}$ formally as independent variables
and in what follows, we shall rewrite system of second-order PDEs (18) as a
larger system of first-order PDEs for the extended set of unknowns
$\\{h^{a}_{\phantom{a}\mu},\partial_{\lambda}h^{a}_{\phantom{a}\mu}\\}$, or
actually, for their equivalents.

In the setting of the teleparallel gravity, $\Lambda$ is not a function of a
general combination of the gradients
$\partial_{\lambda}h^{a}_{\phantom{a}\nu}$ but of their special combination,
that is torsion. Yet, we shall employ not the torsion directly but its Hodge
dual, i.e. we assume that

$\Lambda(h^{a}_{\phantom{a}\mu},\partial_{\lambda}h^{a}_{\phantom{a}\nu})=L(h^{a}_{\phantom{a}\mu},\accentset{\star}{T}^{a\mu\nu}),$
(20)

where $\accentset{\star}{T}^{a\mu\nu}$ is the Hodge dual to the torsion, i.e.

$\displaystyle\accentset{\star}{T}^{a\mu\nu}$
$\displaystyle:=\frac{1}{2}\bm{\in}^{\mu\nu\rho\sigma}T^{a}_{\phantom{a}\rho\sigma}=\bm{\in}^{\mu\nu\rho\sigma}\partial_{\rho}h^{a}_{\phantom{a}\sigma},$
(21a) $\displaystyle T^{a}_{\phantom{a}\mu\nu}$
$\displaystyle=-\frac{1}{2}\bm{\in}_{\mu\nu\rho\sigma}\accentset{\star}{T}^{a\rho\sigma}.$
(21b)

It is important to emphasize that we deliberately chose to define the Hodge
dual using the Levi-Civita symbol $\bm{\in}^{\lambda\mu\nu\rho}$ and _not_ the
Levi-Civita tensor
$\varepsilon^{\lambda\mu\nu\rho}=h^{-1}\bm{\in}^{\lambda\mu\nu\rho}$ that will
be important later for the so-called integrability condition (23). Remark that
according to definition (21a), $\accentset{\star}{T}^{a\mu\nu}$ is a _tensor-
density_ of weight $+1$.

In terms of the Lagrangian density
$L(h^{a}_{\phantom{a}\mu},\accentset{\star}{T}^{a\mu\nu})$, using notations
(20) and definitions (21a), we can instead rewrite Euler-Lagrange equations
(18) as

$\partial_{\nu}(\bm{\in}^{\mu\nu\lambda\rho}L_{\accentset{\star}{T}^{a\lambda\rho}})=-L_{h^{a}_{\phantom{a}\mu}}.$
(22)

The latter has to be supplemented by the integrability condition

$\partial_{\nu}\accentset{\star}{T}^{a\mu\nu}=0,$ (23)

which is a trivial consequence of the definition of the Hodge dual (21a), i.e.
of the commutativity property of the standard spacetime derivative
$\partial_{\mu}$, and the identity (15). We note that if the Hodge dual was
defined using the Levi-Civita tensor $\varepsilon^{\mu\nu\rho\sigma}$ instead
of the Levi-Civita symbol, then one would have that
$\partial_{\mu}\accentset{\star}{T}^{a\mu\nu}\neq 0$.

Another consequence of the commutative property of $\partial_{\mu}$ and the
definition of the Hodge dual (based on the Levi-Civita symbol) is that the
Noether energy-momentum current density

$J^{\mu}_{\phantom{\mu}a}:=L_{h^{a}_{\phantom{a}\mu}}$ (24)

is conserved in the ordinary sense:

$\partial_{\mu}J^{\mu}_{\phantom{\mu}a}=0.$ (25)

If equations (22), (23) and (25) are accompanied with the torsion definition

$\displaystyle\partial_{\mu}h^{a}_{\phantom{a}\nu}-\partial_{\nu}h^{a}_{\phantom{a}\mu}$
$\displaystyle=T^{a}_{\phantom{a}\mu\nu},\qquad$ (26a)

they form the following system of _first-order_ partial differential equations
(only first-order derivatives are involved)

$\displaystyle\partial_{\nu}(\bm{\in}^{\mu\nu\lambda\rho}L_{\accentset{\star}{T}^{a\lambda\rho}})$
$\displaystyle=-L_{h^{a}_{\phantom{a}\mu}},$ (27a)
$\displaystyle\partial_{\nu}\accentset{\star}{T}^{a\mu\nu}$ $\displaystyle=0,$
(27b) $\displaystyle\partial_{\mu}L_{h^{a}_{\phantom{a}\mu}}$
$\displaystyle=0,$ (27c)
$\displaystyle\partial_{\mu}h^{a}_{\phantom{a}\nu}-\partial_{\nu}h^{a}_{\phantom{a}\mu}$
$\displaystyle=T^{a}_{\phantom{a}\mu\nu},$ (27d)

for the unknowns
$\\{h^{a}_{\phantom{a}\mu},\accentset{\star}{T}^{a\mu\nu}\\}$.

This system forms a base on which we shall build our $3+1$-split of TEGR in
Sections 7 and 8.

##  6 Energy-momentum balance laws

Any conservation law written as a 4-_ordinary_ divergence is a true
conservation law, meaning that it yields a time-conserved “charge” [2]. Hence,
the Noether current $J^{\mu}_{\phantom{\mu}a}=L_{h^{a}_{\phantom{a}\mu}}$ is a
conserved charge in the ordinary sense, see (25), (27c). It expresses the
conservation of the total555The “total” here means the gravitational +
matter/electromagnetic energy-momentum current, i.e.
$J^{a}_{\phantom{a}\mu}=L^{\text{(g)}}_{h^{a}_{\phantom{a}\mu}}+L^{\text{(m)}}_{h^{a}_{\phantom{a}\mu}}$.
energy-momentum current density.

However, its spacetime counterpart

$\sigma^{\mu}_{\ \,\nu}:=h^{a}_{\phantom{a}\nu}L_{h^{a}_{\phantom{a}\mu}},$
(28)

which can be called the total energy-momentum tensor density, is conserved
neither ordinary nor covariantly.

Indeed, after contracting with $h^{a}_{\phantom{a}\nu}$ and adding to it
$0\equiv
L_{h^{a}_{\phantom{a}\mu}}\partial_{\mu}h^{a}_{\phantom{a}\nu}-L_{h^{a}_{\phantom{a}\mu}}\partial_{\mu}h^{a}_{\phantom{a}\nu}=L_{h^{a}_{\phantom{a}\mu}}\partial_{\mu}h^{a}_{\phantom{a}\nu}-L_{h^{a}_{\phantom{a}\mu}}h^{a}_{\phantom{a}\lambda}W^{\lambda}_{\phantom{\lambda}\mu\nu}$,
Noether current conservation law (25) can be rewritten in a pure spacetime
form:

$\partial_{\mu}\sigma^{\mu}_{\ \,\nu}=\sigma^{\mu}_{\
\,\lambda}W^{\lambda}_{\phantom{\lambda}\mu\nu},$ (29)

which has a production term on the right hand-side and, therefore,
$\sigma^{\mu}_{\ \,\nu}$ is not a conserved quantity in the ordinary sense.
Exactly this equation and not the conservation law (27c) will be used in the
$3+1$-split of the TEGR equations. This is our personal preference
conditioned, in particular, by the wish to put the $3+1$ TEGR equations in the
class of _symmetric hyperbolic_ equations and to study their Hamiltonian
structure [41, 40] in a followup paper.

On the other hand, if

$\displaystyle\mathcal{D}_{\lambda}V^{\mu}$
$\displaystyle=\partial_{\lambda}V^{\mu}+V^{\rho}W^{\mu}_{\phantom{\mu}\lambda\rho},$
(30a) $\displaystyle\mathcal{D}_{\lambda}V_{\mu}$
$\displaystyle=\partial_{\lambda}V_{\mu}-V_{\rho}W^{\rho}_{\phantom{\rho}\lambda\mu}$
(30b)

is the covariant derivative of the Weitzenböck connection, and keeping in mind
that $h^{a}_{\phantom{a}\nu}L_{h^{a}_{\phantom{a}\mu}}$ is a tensor density of
weight $+1$, balance law (29) can be rewritten as a covariant divergence with
a production term

$\mathcal{D}_{\mu}\sigma^{\mu}_{\ \,\nu}=-\sigma^{\mu}_{\
\,\nu}T^{\rho}_{\phantom{\rho}\mu\rho},$

and hence, $\sigma^{\mu}_{\ \,\nu}$ does not conserved also in the covariant
sense.

For later needs, the following expression of the energy momentum
$\sigma^{\mu}_{\ \,\nu}=h^{a}_{\phantom{a}\nu}L_{h^{a}_{\phantom{a}\mu}}$ is
required

$\sigma^{\mu}_{\
\,\nu}=2\accentset{\star}{T}^{a\lambda\mu}L_{\accentset{\star}{T}^{a\lambda\nu}}-(\accentset{\star}{T}^{a\lambda\rho}L_{\accentset{\star}{T}^{a\lambda\rho}}-L)\delta^{\mu}_{\
\nu}.$ (31)

which is valid for the TEGR Lagrangian discussed in Sec. 12. This formula will
be used later in the $3+1$-split and is analogous to [2, Eq.(10.13)].

##  7 Preliminaries for $3+1$ split

###  7.1 Transformation of the torsion equations

Before performing a $3+1$-split [1] of system (27), we need to do some
preliminary transformations of every equation in (27).

Similar to electromagnetism, we define the following quantities

$E^{a}_{\phantom{a}\mu}:=T^{a}_{\phantom{a}\mu\nu}u^{\nu},\qquad
B^{a\mu}:=\accentset{\star}{T}^{a\mu\nu}u_{\nu}$ (32)

Note that $E^{a}_{\phantom{a}\mu}$ is a rank-1 spacetime tensor, while
$B^{a\mu}$ is a tensor-density.

It is known that for any skew-symmetric tensor, its Hodge dual, and a time-
like vector $u^{\mu}$ the following decompositions hold

$\displaystyle\accentset{\star}{T}^{a\mu\nu}$
$\displaystyle=u^{\mu}B^{a\nu}-u^{\nu}B^{a\mu}+\bm{\in}^{\mu\nu\lambda\rho}u_{\lambda}E^{a}_{\phantom{a}\rho},$
(33a) $\displaystyle T^{a}_{\phantom{a}\mu\nu}$
$\displaystyle=u_{\mu}E^{a}_{\phantom{a}\nu}-u_{\nu}E^{a}_{\phantom{a}\mu}-\bm{\in}_{\mu\nu\lambda\rho}u^{\lambda}B^{a\rho}.$
(33b)

Furthermore, we assume that the Lagrangian density can be equivalently
expressed in different sets of variables, i.e.

$L(h^{a}_{\phantom{a}\mu},\accentset{\star}{T}^{a\mu\nu})=\mathfrak{L}(h^{a}_{\phantom{a}\mu},T^{a}_{\phantom{a}\mu\nu})=\mathcal{L}(h^{a}_{\phantom{a}\mu},B^{a\mu},E^{a}_{\phantom{a}\nu}).$
(34)

It then can be shown that the derivatives of these different representations
of the Lagrangian are related as

$\displaystyle L_{\accentset{\star}{T}^{a\mu\nu}}u^{\nu}$
$\displaystyle=-\frac{1}{2}\left(\mathcal{L}_{B^{a\mu}}+u_{\mu}\mathcal{L}_{B^{a\lambda}}u^{\lambda}\right),$
(35a) $\displaystyle\mathfrak{L}_{T^{a}_{\phantom{a}\mu\nu}}u_{\nu}$
$\displaystyle=-\frac{1}{2}\left(\mathcal{L}_{E^{a}_{\phantom{a}\mu}}+u^{\mu}\mathcal{L}_{E^{a}_{\phantom{a}\lambda}}u_{\lambda}\right),$
(35b)

and hence, (27a) and (27b) can be written as (see Appendix (C))

$\displaystyle\partial_{\nu}(u^{\mu}\mathcal{L}_{E^{a}_{\phantom{a}\nu}}-u^{\nu}\mathcal{L}_{E^{a}_{\phantom{a}\mu}}+\bm{\in}^{\mu\nu\lambda\rho}u_{\lambda}\mathcal{L}_{B^{a\rho}})$
$\displaystyle=J^{\mu}_{\phantom{\mu}a}$ (36a)
$\displaystyle\partial_{\nu}(u^{\mu}B^{a\nu}-u^{\nu}B^{a\mu}+\bm{\in}^{\mu\nu\lambda\rho}u_{\lambda}E^{a}_{\phantom{a}\rho})$
$\displaystyle=0,$ (36b)

where the source $J^{\mu}_{\phantom{\mu}a}=L_{h^{a}_{\phantom{a}\mu}}$ has yet
to be developed.

Let us now introduce a new potential
$U(h^{a}_{\phantom{a}\mu},B^{a\mu},D_{\phantom{\mu}a}^{\mu})$ as a partial
Legendre transform of the Lagrangian $\mathcal{L}$

$U(h^{a}_{\phantom{a}\mu},B^{a\mu},D_{\phantom{\mu}a}^{\mu}):=E^{a}_{\phantom{a}\lambda}\mathcal{L}_{E^{a}_{\phantom{a}\lambda}}-\mathcal{L}.$
(37)

By abusing a little bit notations for $B^{a\mu}$ (we shall use the same letter
for $B^{a\mu}$ and $-B^{a\mu}$, this is an intermediate change of variables
and will not appear in the final formulation), we introduce the new state
variables

$D_{\phantom{\mu}a}^{\mu}:=\mathcal{L}_{E^{a}_{\phantom{a}\mu}},\qquad
B^{a\mu}:=-B^{a\mu},\qquad h^{a}_{\phantom{a}\mu}:=h^{a}_{\phantom{a}\mu}.$
(38)

Note that both $D_{\phantom{\mu}a}^{\mu}$ and $B^{a\mu}$ are _tensor-
densities_. For derivatives of the new potential, we have the following
relations

$U_{D_{\phantom{\mu}a}^{\mu}}=E^{a}_{\phantom{a}\mu},\quad
U_{B^{a\mu}}=\mathcal{L}_{B^{a\mu}},\quad
U_{h^{a}_{\phantom{a}\mu}}=-\mathcal{L}_{h^{a}_{\phantom{a}\mu}}.$ (39)

This allows us to rewrite equations (36) in the form similar to the non-linear
electrodynamics of moving media [38, 20, 31]

$\displaystyle\partial_{\nu}(u^{\mu}D_{\phantom{\nu}a}^{\nu}-u^{\nu}D_{\phantom{\mu}a}^{\mu}+\bm{\in}^{\mu\nu\lambda\rho}u_{\lambda}U_{B^{a\rho}})$
$\displaystyle=J^{\mu}_{\phantom{\mu}a},$ (40a)
$\displaystyle\partial_{\nu}(u^{\mu}B^{a\nu}-u^{\nu}B^{a\mu}-\bm{\in}^{\mu\nu\lambda\rho}u_{\lambda}U_{D_{\phantom{\rho}a}^{\rho}})$
$\displaystyle=0,$ (40b)

with $B^{a\mu}$ and $D_{\phantom{\mu}a}^{\mu}$ being the analog of the
magnetic and electric fields, accordingly.

Finally, we need to express also the Noether current
$J^{\mu}_{\phantom{\mu}a}=L_{h^{a}_{\phantom{a}\mu}}$ in terms of the new
potential $U$ and the fields $D_{\phantom{\mu}a}^{\mu}$ and $B^{a\mu}$. One
has (see details in Appendix A)

$J^{\mu}_{\phantom{\mu}\hat{{0}}}=-U_{h^{\hat{{0}}}_{\phantom{\hat{{0}}}\mu}}+u^{\lambda}B^{b\mu}U_{B^{b\lambda}}-u^{\mu}B^{b\lambda}U_{B^{b\lambda}}\\\
-u^{\mu}D_{\phantom{\lambda}b}^{\lambda}U_{D_{\phantom{\lambda}b}^{\lambda}}+\bm{\in}^{\mu\lambda\rho\sigma}u_{\rho}U_{B^{b\lambda}}U_{D_{\phantom{\sigma}b}^{\sigma}},$
(41a) $J^{\mu}_{\phantom{\mu}A}=-U_{h^{A}_{\phantom{A}\mu}}\\\
-\eta^{\nu}_{\phantom{\nu}A}u^{\mu}\left(u_{\lambda}D_{\phantom{\lambda}b}^{\lambda}U_{D_{\phantom{\nu}b}^{\nu}}-\bm{\in}_{\nu\lambda\rho\sigma}u^{\rho}B^{b\sigma}D_{\phantom{\lambda}b}^{\lambda}\right).$
(41b)

###  7.2 Transformation of the tetrad equations

Contracting (27d) with the 4-velocity $u^{\mu}$, and then after change of
variables (38) and (39), the resulting equation reads as

$(\partial_{\mu}h^{a}_{\phantom{a}\nu}-\partial_{\nu}h^{a}_{\phantom{a}\mu})u^{\nu}=U_{D_{\phantom{\mu}a}^{\mu}}.$
(42)

Furthermore, using the identity
$\eta^{\mu}_{\phantom{\mu}b}\partial_{\nu}h^{a}_{\phantom{a}\nu}=-h^{a}_{\phantom{a}\nu}\partial_{\nu}\eta^{\mu}_{\phantom{\mu}b}$
and the definition $u^{\mu}=\eta^{\mu}_{\phantom{\mu}\hat{{0}}}$, the latter
equation can be rewritten as

$u^{\nu}\partial_{\nu}h^{a}_{\phantom{a}\mu}+h^{a}_{\phantom{a}\nu}\partial_{\mu}u^{\nu}=-U_{D_{\phantom{\mu}a}^{\mu}},$
(43)

that later will be used in the $3+1$-split.

###  7.3 Transformation of the energy-momentum

Finally, we express the gravitational part of the energy-momentum tensor
${\sigma^{\mu}_{\ \,\nu}}$ (31) in terms of new variables (38) and the
potential $U(h^{a}_{\phantom{a}\mu},B^{a\mu},D_{\phantom{\mu}a}^{\mu})$, while
we keep energy-momentum equation (29) unchanged. It reads

$\displaystyle{\sigma^{\mu}_{\ \,\nu}}=$
$\displaystyle-B^{a\mu}U_{B^{a\nu}}-D_{\phantom{\nu}a}^{\nu}U_{D_{\phantom{\mu}a}^{\mu}}$
$\displaystyle-u^{\lambda}u_{\nu}B^{a\mu}U_{B^{a\lambda}}-u^{\mu}u_{\lambda}D_{\phantom{\lambda}a}^{\lambda}U_{D_{\phantom{\nu}a}^{\nu}}$
$\displaystyle+u^{\mu}u_{\nu}B^{a\lambda}U_{B^{a\lambda}}+u^{\mu}u_{\nu}D_{\phantom{\lambda}a}^{\lambda}U_{D_{\phantom{\lambda}a}^{\lambda}}$
$\displaystyle+\bm{\in}_{\nu\sigma\lambda\rho}u^{\mu}u^{\sigma}B^{a\lambda}D_{\phantom{\rho}a}^{\rho}+\bm{\in}^{\mu\sigma\lambda\rho}u_{\nu}u_{\sigma}U_{B^{a\lambda}}U_{D_{\phantom{\rho}a}^{\rho}}$
$\displaystyle+(B^{a\lambda}U_{B^{a\lambda}}+D_{\phantom{\lambda}a}^{\lambda}U_{D_{\phantom{\lambda}a}^{\lambda}}-U)\delta^{\mu}_{\
\nu}.$ (44)

We shall need this expression for $\sigma^{\mu}_{\ \,\nu}$ in the last part of
the derivation of the $3+1$ equations.

Therefore, the TEGR system in its intermediate form for the unknowns
$\\{h^{a}_{\phantom{a}\mu},D_{\phantom{\mu}a}^{\mu},B^{a\mu}\\}$ reads

$\displaystyle\partial_{\nu}(u^{\mu}D_{\phantom{\nu}a}^{\nu}-u^{\nu}D_{\phantom{\mu}a}^{\mu}+\bm{\in}^{\mu\nu\lambda\rho}u_{\lambda}U_{B^{a\rho}})$
$\displaystyle=J^{\mu}_{\phantom{\mu}a},$ (45a)
$\displaystyle\partial_{\nu}(u^{\mu}B^{a\nu}-u^{\nu}B^{a\mu}-\bm{\in}^{\mu\nu\lambda\rho}u_{\lambda}U_{D_{\phantom{\rho}a}^{\rho}})$
$\displaystyle=0,$ (45b) $\displaystyle\partial_{\mu}J^{\mu}_{\phantom{\mu}a}$
$\displaystyle=0,$ (45c) $\displaystyle
u^{\nu}\partial_{\nu}h^{a}_{\phantom{a}\mu}+h^{a}_{\phantom{a}\nu}\partial_{\mu}u^{\nu}$
$\displaystyle=-U_{D_{\phantom{\mu}a}^{\mu}},$ (45d)

with $J^{\mu}_{\phantom{\mu}a}$ given by (41). After we introduce a particular
choice of the observer’s 4-velocity $u^{\mu}$ at the beginning of the next
section, we shall finalize transformation of system (45) to present the final
$3+1$ equations of TEGR.

##  8 $3+1$ split of the TEGR equations

In this section, we derive a $3+1$ version of system (45) that can be used in
a computational code for numerical relativity.

We first recall that Latin indices from the middle of the alphabet
$i,j,k,\ldots=1,2,3$ are used to denote the spatial components of the
spacetime tensors, and Latin indices
$A,B,C,\ldots=\hat{{1}},\hat{{2}},\hat{{3}}$ to denote the spatial directions
in the tangent Minkowski space. Additionally, we use the hat on top of a
number, e.g. $\hat{{0}}$, for the indices marking the time and space direction
in the tangent space in order to distinguish them from the indices of the
spacetime tensors. Also recall that observer’s 4-velocity $u^{\mu}$ is
associated with the $\hat{0}$-th column of the inverse tetrad
$\eta^{\mu}_{\phantom{\mu}a}$, while the covariant components $u_{\mu}$ of the
4-velocity with the $\hat{0}$-th row of the frame field. For $u^{\mu}$ and
$u_{\mu}$ we standardly assume [1, 45]:

$\displaystyle u^{\mu}$
$\displaystyle=\eta^{\mu}_{\phantom{\mu}\hat{{0}}}=\alpha^{-1}(1,-\beta^{i}),$
(46a) $\displaystyle u_{\mu}$
$\displaystyle=-h^{\hat{{0}}}_{\phantom{\hat{{0}}}\mu}=(-\alpha,0,0,0),$ (46b)

with $\alpha$ being the _lapse function_ , and $\beta^{i}$ being the _shift
vector_. One can write down $h^{a}_{\phantom{a}\mu}$ and
$\eta^{\mu}_{\phantom{\mu}a}$ explicitly :

$\bm{h}=\left(\begin{array}[]{cccc}\alpha&0&0&0\\\\[2.84526pt]
\beta^{\hat{{1}}}&&&\\\ \beta^{\hat{{2}}}&&h^{A}_{\phantom{A}i}&\\\
\beta^{\hat{{3}}}&&&\end{array}\right),$ (47a)
$\bm{h}^{-1}=\bm{\eta}=\left(\begin{array}[]{rlcl}1/\alpha&0&0&0\\\\[2.84526pt]
-\beta^{1}/\alpha&&&\\\
-\beta^{2}/\alpha&&\left(h^{A}_{\phantom{A}i}\right)^{-1}&\\\
-\beta^{3}/\alpha&&&\end{array}\right),$ (47b)

where $\beta^{A}=h^{A}_{\phantom{A}i}\beta^{i}$. The metric tensor and its
inverse are (e.g. see [27])

$\displaystyle g_{\mu\nu}$
$\displaystyle=\left(\begin{array}[]{cc}-\alpha^{2}+\beta_{i}\beta^{i}&\beta_{i}\\\\[2.84526pt]
\beta_{i}&\gamma_{ij}\\\ \end{array}\right),$ (48c) $\displaystyle g^{\mu\nu}$
$\displaystyle=\left(\begin{array}[]{cc}-1/\alpha^{2}&\beta^{i}/\alpha^{2}\\\\[2.84526pt]
\beta^{i}/\alpha^{2}&\gamma^{ij}-\beta^{i}\beta^{j}/\alpha^{2}\\\
\end{array}\right),$ (48f)

where $\gamma_{ij}=\kappa_{AB}h^{A}_{\phantom{A}i}h^{B}_{\phantom{B}j}$,
$\gamma^{ij}=\left(\gamma_{ij}\right)^{-1}$, and
$\beta_{i}=\gamma_{ij}\beta^{j}$.

In the rest of the paper, $h_{3}$ stands for $\det(h^{A}_{\phantom{A}i})$ so
that

$h:=\det(h^{a}_{\phantom{a}\mu})=\alpha h_{3},$ (49)

and we use the following convention for the three-dimensional Levi-Civita
symbol

$\bm{\in}^{0ijk}=\bm{\in}^{ijk},\qquad\bm{\in}_{0ijk}=-\bm{\in}_{ijk}.$ (50)

###  8.1 $3+1$ split of the torsion PDEs

##### Case $\mu=i=1,2,3$.

For our choice of observer’s velocity (46), equations (45a) and (45b) read

$\partial_{t}(\alpha^{-1}D_{\phantom{i}a}^{i})+\partial_{k}(\beta^{i}\alpha^{-1}D_{\phantom{k}a}^{k}-\beta^{k}\alpha^{-1}D_{\phantom{i}a}^{i}\\\
-\bm{\in}^{ikj}\alpha\,U_{B^{aj}})=-J^{i}_{\phantom{i}a},$ (51a)
$-\partial_{t}(\alpha^{-1}B^{ai})+\partial_{k}(\beta^{k}\alpha^{-1}B^{ai}-\beta^{i}\alpha^{-1}B^{ak}\\\
+\bm{\in}^{ikj}\alpha\,U_{D_{\phantom{j}a}^{j}})=0.$ (51b)

Because of the factor $\alpha^{-1}$ everywhere in these equations, it is
convenient to rescale the variables:

$\mathcal{D}_{\phantom{\mu}a}^{\mu}:=\alpha^{-1}D_{\phantom{\mu}a}^{\mu},\qquad\mathcal{B}^{a\mu}:=-\alpha^{-1}B^{a\mu}$
(52)

so that the derivatives of the potential
$\mathcal{U}(h^{a}_{\phantom{a}\mu},\mathcal{D}_{\phantom{\mu}a}^{\mu},\mathcal{B}^{a\mu}):=U(h^{a}_{\phantom{a}\mu},D_{\phantom{\mu}a}^{\mu},B^{a\mu})$
transform as

$\mathcal{U}_{\mathcal{D}_{\phantom{\mu}a}^{\mu}}=\alpha\,U_{D_{\phantom{\mu}a}^{\mu}},\qquad\mathcal{U}_{\mathcal{B}^{a\mu}}=-\alpha\,U_{B^{a\mu}}.$
(53)

Hence, (51) can be rewritten as

$\displaystyle\partial_{t}\mathcal{D}_{\phantom{i}a}^{i}+\partial_{k}(\beta^{i}\mathcal{D}_{\phantom{k}a}^{k}-\beta^{k}\mathcal{D}_{\phantom{i}a}^{i}-\bm{\in}^{ikj}\,\mathcal{U}_{\mathcal{B}^{aj}})$
$\displaystyle=-J^{i}_{\phantom{i}a},$ (54a)
$\displaystyle\partial_{t}\mathcal{B}^{ai}+\partial_{k}(\,\beta^{i}\mathcal{B}^{ak}-\beta^{k}\mathcal{B}^{ai}+\bm{\in}^{ikj}\mathcal{U}_{\mathcal{D}_{\phantom{j}a}^{j}})$
$\displaystyle=0.$ (54b)

Finally, extending formally the shift vector as $\bm{\beta}=(-1,\beta^{k})$
and introducing the change of variables

$\mathtt{D}_{\phantom{\mu}a}^{\mu}:=\mathcal{D}_{\phantom{\mu}a}^{\mu}+\beta^{\mu}\mathcal{D}_{\phantom{0}a}^{0},\qquad\mathtt{B}^{a\mu}:=\mathcal{B}^{a\mu}$
(55)

which are essentially $3$-by-$3$ matrices ($\mathtt{D}_{\phantom{0}a}^{0}=0$,
$\mathtt{B}^{a0}=\mathtt{B}^{\hat{{0}}\mu}=0$, and
$\mathtt{D}_{\phantom{k}\hat{{0}}}^{k}$ are given by (85a)), we arrive at the
final form of the $3+1$ equations for the fields
$\mathtt{D}_{\phantom{i}A}^{i}$ and $\mathtt{B}^{Ai}$

$\displaystyle\partial_{t}\mathtt{D}_{\phantom{i}A}^{i}+\partial_{k}(\beta^{i}\mathtt{D}_{\phantom{k}A}^{k}-\beta^{k}\mathtt{D}_{\phantom{i}A}^{i}-\bm{\in}^{ikj}\,\mathtt{H}_{jA})$
$\displaystyle=-J^{i}_{\phantom{i}A},$ (56a)
$\displaystyle\partial_{t}\mathtt{B}^{Ai}+\partial_{k}(\,\beta^{i}\mathtt{B}^{Ak}-\beta^{k}\mathtt{B}^{Ai}+\bm{\in}^{ikj}\mathtt{E}^{A}_{\phantom{A}j})$
$\displaystyle=0.$ (56b)

where

$\displaystyle\mathtt{E}^{A}_{\phantom{A}j}=\mathcal{U}_{\mathcal{D}_{\phantom{j}A}^{j}}=\alpha
U_{D_{\phantom{j}A}^{j}}=\alpha E^{A}_{\phantom{A}j},$ (57a)
$\displaystyle\mathtt{H}_{jA}=\mathcal{U}_{\mathcal{B}^{Aj}}=-\alpha
U_{B^{Aj}}.$ (57b)

In terms of $\\{\mathtt{B}^{Ak},\mathtt{D}_{\phantom{k}A}^{k}\\}$ the
potential
$\mathcal{U}(h^{a}_{\phantom{a}\mu},\mathcal{D}_{\phantom{\mu}a}^{\mu},\mathcal{B}^{a\mu})$
will be denoted as
$\mathtt{U}(h^{a}_{\phantom{a}\mu},\mathtt{B}^{a\mu},\mathtt{D}_{\phantom{\mu}a}^{\mu})$
and it reads (see (12))

$\displaystyle\mathtt{U}(h^{A}_{\phantom{A}k},\mathtt{B}^{Ak},\mathtt{D}_{\phantom{k}A}^{k})$
$\displaystyle=$ $\displaystyle-\frac{\alpha}{2h_{3}}\bigg{(}$
$\displaystyle\varkappa\left(\mathtt{D}_{\phantom{A}B}^{A}\mathtt{D}_{\phantom{B}A}^{B}-\frac{1}{2}(\mathtt{D}_{\phantom{A}A}^{A})^{2}\right)$
$\displaystyle+$
$\displaystyle\frac{1}{\varkappa}\left(\mathtt{B}^{A}_{\phantom{A}B}\mathtt{B}^{B}_{\phantom{B}A}-\frac{1}{2}(\mathtt{B}^{A}_{\phantom{A}A})^{2}\right)\bigg{)},$
(58)

Alternatively, it can be computed as

$\mathtt{U}=\frac{1}{2}(\mathtt{D}_{\phantom{k}a}^{k}\mathtt{E}^{a}_{\phantom{a}k}+\mathtt{B}^{ak}\mathtt{H}_{ka}).$
(59)

##### Case $\mu=0$.

The $0$-th equations (45a) and (45b) are actually not time-evolution equations
but pure spatial (stationary) constraints

$\partial_{k}\mathtt{D}_{\phantom{k}a}^{k}=J^{0}_{\phantom{0}a},\qquad\partial_{k}\mathtt{B}^{ak}=0.$
(60)

Their possible violation during the numerical integration of the time-
evolution equations (56) is a well-known problem in the numerical analysis of
hyperbolic equations with involution constraints. Different strategies to
preserve such stationary constraints are known, e.g. constraint-cleaning
approach [36, 15, 16] or constraint-preserving integrators [39, 10].

###  8.2 Tetrad PDE

Because of our choice for the tetrad (47), we are only interested in the
evolution equations for the components $h^{A}_{\phantom{A}k}$. Thus, using the
definition of the 4-velocity (46), equation (45d) can be written as

$\partial_{t}h^{A}_{\phantom{A}k}-\beta^{i}\partial_{i}h^{A}_{\phantom{A}k}-h^{A}_{\phantom{A}i}\partial_{k}\beta^{i}=-\mathtt{E}^{A}_{\phantom{A}k}.$
(61)

This equation can be also written in a slightly different form. After adding
$0\equiv-\beta^{i}\partial_{k}h^{A}_{\phantom{A}i}+\beta^{i}\partial_{k}h^{A}_{\phantom{A}i}$
to the left hand-side of (61), one has

$\partial_{t}h^{A}_{\phantom{A}k}-\partial_{k}(\beta^{i}h^{A}_{\phantom{A}i})-\beta^{i}(\partial_{i}h^{A}_{\phantom{A}k}-\partial_{k}h^{A}_{\phantom{A}i})=-\mathtt{E}^{A}_{\phantom{A}k}.$
(62)

Finally, using the definition of $B^{A\mu}$ and $u_{\mu}$, we have that

$B^{A\mu}=\accentset{\star}{T}^{A\mu\nu}u_{\nu}=u_{\nu}\bm{\in}^{\mu\nu\alpha\beta}\partial_{\alpha}h^{A}_{\phantom{A}\beta}=\\\
-u_{\nu}\bm{\in}^{\nu\mu\alpha\beta}\partial_{\alpha}h^{A}_{\phantom{A}\beta}=\alpha\bm{\in}^{0\mu\alpha\beta}\partial_{\alpha}h^{A}_{\phantom{A}\beta}$
(63)

and hence (use that $\bm{\in}^{0ijk}=\bm{\in}^{ijk}$)

$\mathtt{B}^{Ak}=\bm{\in}^{kij}\partial_{i}h^{A}_{\phantom{A}j},$ (64)

or

$-\beta^{i}(\partial_{i}h^{A}_{\phantom{A}k}-\partial_{k}h^{A}_{\phantom{A}i})=\bm{\in}_{klj}\beta^{l}\mathtt{B}^{Aj}.$
(65)

Therefore, (62) can be also rewritten as

$\partial_{t}h^{A}_{\phantom{A}k}-\partial_{k}(\beta^{i}h^{A}_{\phantom{A}i})=-\mathtt{E}^{A}_{\phantom{A}k}-\bm{\in}_{klj}\beta^{l}\mathtt{B}^{Aj}.$
(66)

It is not clear yet which one of these equivalent forms (61), (62), or (66) at
the continuous level is more suited for the numerical discretization. However,
for a structure-preserving integrator, e.g. [10], which is able to preserve
(64) up to the machine precision, all these forms are equivalent.

###  8.3 Energy-momentum PDE

We now turn to the energy-momentum evolution equation

$\partial_{\mu}\sigma^{\mu}_{\ \,\nu}=\sigma^{\mu}_{\
\,\lambda}W^{\lambda}_{\phantom{\lambda}\mu\nu}$ (67)

and provide its $3+1$ version. Because the Lagrangian density in (20)
represents the sum of the gravitational and matter Lagrangian, the energy-
momentum tensor is also assumed to be the sum of the gravity and matter parts:
$\sigma^{\mu}_{\ \,\nu}=\prescript{\text{(m)}}{}{\sigma^{\mu}_{\
\,\nu}}+\prescript{\text{(g)}}{}{\sigma^{\mu}_{\ \,\nu}}$. However, in what
follows, we shall omit the matter part keeping in mind that the energy and
momentum equations discussed below are equations for the total
(matter+gravity) energy-momentum.

We first, explicitly split (67) into the three momentum and one energy
equation:

$\displaystyle\partial_{t}\sigma^{0}_{\ \,i}+\partial_{k}\sigma^{k}_{\ \,i}$
$\displaystyle=\sigma^{\mu}_{\ \,\lambda}W^{\lambda}_{\phantom{\lambda}\mu
i},$ (68a) $\displaystyle\partial_{t}\sigma^{0}_{\
\,0}+\partial_{k}\sigma^{k}_{\ \,0}$ $\displaystyle=\sigma^{\mu}_{\
\,\lambda}W^{\lambda}_{\phantom{\lambda}\mu 0}.$ (68b)

Applying the change of variables (52), (53), and (55), expression (7.3) for
the energy-momentum reduces to the following expressions

  * •

for the total momentum density, $\sigma^{0}_{\ i}=\rho_{i}$ (the matter part
is left unspecified and omitted in this paper):

$\rho_{i}:=-\bm{\in}_{ilj}\mathtt{B}^{al}\mathtt{D}_{\phantom{j}a}^{j},$ (69)

  * •

for the momentum flux, $\sigma^{k}_{\ \,i}$:

$\sigma^{k}_{\
\,i}=-\mathtt{B}^{ak}\mathtt{H}_{ia}-\mathtt{D}_{\phantom{k}a}^{k}\mathtt{E}^{a}_{\phantom{a}i}+\bm{\in}_{ilj}\beta^{k}\mathtt{B}^{al}\mathtt{D}_{\phantom{j}a}^{j}\\\
+(\mathtt{B}^{aj}\mathtt{H}_{ja}+\mathtt{D}_{\phantom{j}a}^{j}\mathtt{E}^{a}_{\phantom{a}j}-\mathtt{U})\delta^{k}_{\
i},$ (70)

  * •

for the energy density:

$\rho_{0}:=\sigma^{0}_{\
\,0}=-\bm{\in}_{ilj}\beta^{i}\mathtt{B}^{al}\mathtt{D}_{\phantom{j}a}^{j}-\mathtt{U}=\beta^{i}\rho_{i}-\mathtt{U},$
(71)

  * •

and for the energy flux, $\sigma^{k}_{\ \,0}$:

${\sigma^{k}_{\
\,0}}=-\beta^{j}(\mathtt{B}^{ak}\mathtt{H}_{ja}+\mathtt{D}_{\phantom{k}a}^{k}\mathtt{E}^{a}_{\phantom{a}j})\\\
+\beta^{k}(\mathtt{B}^{aj}\mathtt{H}_{ja}+\mathtt{D}_{\phantom{j}a}^{j}\mathtt{E}^{a}_{\phantom{a}j})\\\
+\bm{\in}_{ilj}\beta^{k}\beta^{i}\mathtt{B}^{al}\mathtt{D}_{\phantom{j}a}^{j}+\bm{\in}^{klj}\mathtt{H}_{la}\mathtt{E}^{a}_{\phantom{a}j}.$
(72)

It is convenient to split the momentum and energy fluxes $\sigma^{k}_{\ \,i}$
and $\sigma^{k}_{\ \,0}$ in (68) into advective and constitutive parts, so
that the finale form of the $3+1$ equations for the total energy-momentum (68)
reads

$\displaystyle\partial_{t}\rho_{i}$
$\displaystyle+\partial_{k}(-\beta^{k}\rho_{i}+s_{\ i}^{k})=f_{i},$ (73a)
$\displaystyle\partial_{t}\rho_{0}$
$\displaystyle+\partial_{k}(-\beta^{k}\rho_{0}+\beta^{i}s_{\
i}^{k}+\bm{\in}^{klj}\mathtt{H}_{la}\mathtt{E}^{a}_{\phantom{a}j})=f_{0},$
(73b) where the gravity part of $s_{\ i}^{k}$ is given by $\displaystyle{s_{\
i}^{k}}=-\mathtt{B}^{ak}\mathtt{H}_{ia}-\mathtt{D}_{\phantom{k}a}^{k}\mathtt{E}^{a}_{\phantom{a}i}+(\mathtt{B}^{aj}\mathtt{H}_{ja}+\mathtt{D}_{\phantom{j}a}^{j}\mathtt{E}^{a}_{\phantom{a}j}-\mathtt{U})\delta^{k}_{\
i}.$ (73c) The source terms in (73a) and (73b) are given by $\displaystyle
f_{i}=$
$\displaystyle-\rho_{k}\eta^{k}_{\phantom{k}A}\mathtt{E}^{A}_{\phantom{A}i}+\rho_{j}\partial_{i}\beta^{j}+s_{\
j}^{k}W^{j}_{\phantom{j}ki},$ (73d) $\displaystyle f_{0}=$
$\displaystyle-\left(\partial_{t}\ln(\alpha)-\beta^{k}\partial_{k}\ln(\alpha)\right)\mathtt{U}$
$\displaystyle+\left(\partial_{t}\beta^{A}-\beta^{k}\partial_{k}\beta^{A}\right)\rho_{j}\eta^{j}_{\phantom{j}A}$
$\displaystyle+\bm{\in}^{klj}\mathtt{H}_{la}\mathtt{E}^{a}_{\phantom{a}j}\partial_{k}\ln(\alpha)+s_{\
j}^{k}\eta^{j}_{\phantom{j}A}\partial_{k}\beta^{A}.$ (73e)

####  8.3.1 Noether current

Because the energy-momentum $\sigma^{k}_{\ \,i}$ was defined as $\sigma^{k}_{\
\,i}=h^{a}_{\phantom{a}i}J^{k}_{\phantom{k}a}=h^{A}_{\phantom{A}i}J^{k}_{\phantom{k}A}$,
we can use the expression for $\sigma^{k}_{\ \,i}$ to have an explicit formula
for the Noether current. Thus, we have

$J^{k}_{\phantom{k}A}=\eta^{i}_{\phantom{i}A}\sigma^{k}_{\
\,i}=\eta^{i}_{\phantom{i}A}(-\beta^{k}\rho_{i}+s_{\
i}^{k})=-\beta^{k}\rho_{A}+s_{\ A}^{k}.$ (74)

In the same way we can find $J^{0}_{\phantom{0}a}$ necessary in (60):
$J^{0}_{\phantom{0}a}=\eta^{\lambda}_{\phantom{\lambda}a}\sigma^{0}_{\
\,\lambda}=\eta^{0}_{\phantom{0}a}\sigma^{0}_{\
\,0}+\eta^{i}_{\phantom{i}a}\sigma^{0}_{\ \,i}$, and hence

$J^{0}_{\phantom{0}a}=\left\\{\begin{array}[]{rl}-\alpha^{-1}\mathtt{U}=\rho_{\hat{{0}}},&a=\hat{{0}},\\\\[8.53581pt]
\eta^{i}_{\phantom{i}A}\rho_{i}=\rho_{A},&A=\hat{{1}},\hat{{2}},\hat{{3}}.\\\
\end{array}\right.$ (75)

####  8.3.2 Alternative form of the energy-momentum equations

In (67), we deliberately use the energy-momentum tensor with both spacetime
indices because we would like to use the SHTC and Hamiltonian structure [41]
of these equations for designing numerical schemes in the future, e.g. [13,
33]. However, the resulting PDEs (73) do not have a fully conservative form
preferable for example when dealing with the shock waves in the matter fields.
Therefore, one may want to use the true conservation law (25) for the total
(matter+gravity) Noether current $J^{\mu}_{\phantom{\mu}a}$ instead of (73).
Thus, in the new notations introduced above, four conservation laws
$\partial_{\mu}J^{\mu}_{\phantom{\mu}a}=0$ now read

$\displaystyle\partial_{t}\rho_{A}+$
$\displaystyle\partial_{k}(-\beta^{k}\rho_{A}+s_{\ A}^{k})=0,$ (76a)
$\displaystyle\partial_{t}\rho_{\hat{{0}}}+$
$\displaystyle\partial_{k}(-\beta^{k}\rho_{\hat{{0}}}+\alpha^{-1}\bm{\in}^{klj}\mathtt{H}_{la}\mathtt{E}^{a}_{\phantom{a}j})=0,$
(76b)

with $\rho_{\hat{{0}}}=J^{0}_{\phantom{0}\hat{{0}}}=-\alpha^{-1}\mathtt{U}$.

####  8.3.3 Evolution of the space volume

As in the computational Newtonian mechanics [19, 10], the evolution of the
tetrad field at the discrete level has to be performed consistently with the
volume/mass conservation law. In TEGR, the equivalent to the volume
conservation in the Newtonian mechanics is

$\partial_{\mu}(hu^{\mu})=-hE^{\mu}_{\phantom{\mu}\mu}$ (77)

which can be obtained after contracting (45d) with $\partial h/\partial
h^{a}_{\phantom{a}\mu}=h\eta^{\mu}_{\phantom{\mu}a}$, and where
$E^{\mu}_{\phantom{\mu}\mu}=\eta^{\mu}_{\phantom{\mu}a}E^{a}_{\phantom{a}\mu}$.

After using (46), (47), and (80b), this balance law becomes

$\partial_{t}h_{3}-\partial_{k}(h_{3}\beta^{k})=-\frac{\alpha\varkappa}{2}\mathtt{D}_{\phantom{i}i}^{i},$
(78)

where $h_{3}=\det(h^{A}_{\phantom{A}k})=\sqrt{\det(\gamma_{ij})}$.

##  9 Summary of the $3+1$ equations

Here, we summarize the $3+1$ TEGR equations and give explicit expressions for
the constitutive fluxes $\mathtt{E}^{A}_{\phantom{A}k}$ and $\mathtt{H}_{kA}$
which then close the specification of the entire system.

###  9.1 Evolution equations

The system of $3+1$ TEGR governing equations reads

$\displaystyle\partial_{t}\mathtt{D}_{\phantom{i}A}^{i}+\partial_{k}(\beta^{i}\mathtt{D}_{\phantom{k}A}^{k}-\beta^{k}\mathtt{D}_{\phantom{i}A}^{i}-\bm{\in}^{ikj}\,\mathtt{H}_{jA})$
$\displaystyle=\beta^{i}\rho_{A}-s_{\ A}^{i},$ (79a)
$\displaystyle\partial_{t}\mathtt{B}^{Ai}+\partial_{k}(\,\beta^{i}\mathtt{B}^{Ak}-\beta^{k}\mathtt{B}^{Ai}+\bm{\in}^{ikj}\mathtt{E}^{A}_{\phantom{A}j})$
$\displaystyle=0,$ (79b)
$\displaystyle\partial_{t}\rho_{A}+\partial_{k}(-\beta^{k}\rho_{A}+s_{\
A}^{k})$ $\displaystyle=0,$ (79c)
$\displaystyle\partial_{t}\rho_{\hat{{0}}}+\partial_{k}(-\beta^{k}\rho_{\hat{{0}}}+\alpha^{-1}\bm{\in}^{klj}\mathtt{H}_{la}\mathtt{E}^{a}_{\phantom{a}j})$
$\displaystyle=0,$ (79d)
$\displaystyle\partial_{t}h^{A}_{\phantom{A}k}-\beta^{i}\partial_{i}h^{A}_{\phantom{A}k}-h^{A}_{\phantom{A}i}\partial_{k}\beta^{i}$
$\displaystyle=-\mathtt{E}^{A}_{\phantom{A}k},$ (79e)

with the total (matter+gravity) momentum
$\rho_{i}:=\rho^{\text{(m)}}_{i}+\rho^{\text{(g)}}_{i}$ given by

$\rho_{i}:=\rho^{\text{(m)}}_{i}-\bm{\in}_{ilj}\mathtt{B}^{Al}\mathtt{D}_{\phantom{j}A}^{j},$

and the total energy density

$\rho_{0}:=\beta^{i}\rho_{i}-\mathtt{U},$

The gravitational part of the momentum flux reads

$\displaystyle s_{\
i}^{k}:=-\mathtt{B}^{ak}\mathtt{H}_{ia}-\mathtt{D}_{\phantom{k}a}^{k}\mathtt{E}^{a}_{\phantom{a}i}+(\mathtt{B}^{aj}\mathtt{H}_{ja}+\mathtt{D}_{\phantom{j}a}^{j}\mathtt{E}^{a}_{\phantom{a}j}-\mathtt{U})\delta^{k}_{\
i},$

In particular, the structure of this system resembles very much the structure
of the nonlinear electrodynamics of moving media already solved numerically in
[20] for example, as well as the structure of the continuum mechanics
equations with torsion [43]. Moreover, despite deep conceptual differences, it
is identical to the new dGREM formulation of the Einstein equations recently
pushed forward in [39]. However, it is important to mention that the
Lagrangian approach adopted here in principle permits to generalize the
formulation to other $f(\mathcal{T})$ theories. A detailed comparison of these
two formulations will be a subject of a future paper.

It is also important to note that the structure of system (79) remains the
same independently of the Lagrangian in use as it can be seen from system (45)
where all the constitutive parts are defined through the derivatives of the
potential (37). If the Lagrangian changes, then only the constitutive
functions $\mathtt{E}^{a}_{\phantom{a}\mu}$ and $\mathtt{H}_{\mu a}$ have to
be recomputed.

###  9.2 Constitutive relations

For the TEGR Lagrangian (106), $\mathtt{E}^{A}_{\phantom{A}k}$ and
$\mathtt{H}_{kA}$ read

$\mathtt{H}_{kA}:=-\frac{\alpha}{\varkappa
h_{3}}\kappa_{AB}\kappa_{CD}\left(h^{D}_{\phantom{D}k}h^{B}_{\phantom{B}j}-\frac{1}{2}h^{D}_{\phantom{D}j}h^{B}_{\phantom{B}k}\right)\mathtt{B}^{Cj}\\\
-\frac{h_{3}}{\varkappa}\bm{\in}_{kjl}\eta^{j}_{\phantom{j}A}A^{l},$ (80a)
$\mathtt{E}^{A}_{\phantom{A}k}:=-\frac{\alpha\varkappa}{h_{3}}\left(h^{A}_{\phantom{A}j}h^{B}_{\phantom{B}k}-\frac{1}{2}h^{A}_{\phantom{A}k}h^{B}_{\phantom{B}j}\right)\mathtt{D}_{\phantom{j}B}^{j}\\\
-\varkappa h_{3}\bm{\in}_{kjl}\eta^{j}_{\phantom{j}B}\kappa^{AB}\Omega^{l}.$
(80b)

where $A_{k}=\mathtt{E}^{\hat{{0}}}_{\phantom{\hat{{0}}}k}$ and
$\Omega_{k}=\mathtt{H}_{k\hat{{0}}}$.

###  9.3 Stationary differential constraints

System (79) is supplemented by several differential constraints. Thus, as
already was mentioned, the $0$-th equations ($\mu=0$) of (45a) are not
actually time-evolution equations but reduce to the so-called Hamiltonian and
momentum stationary divergence-type constraints that must hold on the solution
to (79) at every time instant:

$\partial_{k}\mathtt{D}_{\phantom{k}\hat{{0}}}^{k}=\rho_{\hat{{0}}},\qquad\partial_{k}\mathtt{D}_{\phantom{k}A}^{k}=\rho_{A},$
(81)

with $\rho_{\hat{{0}}}=-\alpha^{-1}\mathtt{U}=J^{0}_{\phantom{0}\hat{{0}}}$.

Accordingly, the $0$-th equation of (45b) gives the divergence constraint on
the $\mathtt{B}^{Ai}$ field

$\partial_{k}\mathtt{B}^{Ak}=0,$ (82)

Finally, from the definition of $\mathtt{B}^{Ai}$, we also have a curl-type
constraint on the spatial components of the tetrad field:

$\bm{\in}^{kij}\partial_{i}h^{A}_{\phantom{A}j}=\mathtt{B}^{Ak}.$ (83)

###  9.4 Algebraic constraints

As a consequence of our choice of observer’s reference frame (46), and the
fact that $E^{a}_{\phantom{a}\mu}u^{\mu}=0$ and $B^{a\mu}u_{\mu}=0$, we have
the following algebraic constraints

$\mathtt{B}^{\hat{{0}}\mu}=0,\qquad\mathtt{B}^{a0}=0,$ (84a)
$\mathtt{E}^{\hat{{0}}}_{\phantom{\hat{{0}}}0}=\beta^{k}\mathtt{E}^{\hat{{0}}}_{\phantom{\hat{{0}}}k},\qquad\mathtt{E}^{\hat{{0}}}_{\phantom{\hat{{0}}}k}=\partial_{k}\alpha,\qquad\mathtt{E}^{A}_{\phantom{A}0}=\beta^{j}\mathtt{E}^{A}_{\phantom{A}j},$
(84b)

$\displaystyle\mathtt{D}_{\phantom{0}a}^{0}=0,\qquad\mathtt{D}_{\phantom{k}\hat{{0}}}^{k}=-\frac{1}{\varkappa
h_{3}}\bm{\in}^{kli}\mathtt{B}_{\phantom{}li},\qquad\mathtt{D}_{\phantom{}ik}=\mathtt{D}_{\phantom{}ki},$
(85a) $\displaystyle\mathtt{H}_{0A}=\beta^{k}\mathtt{H}_{kA},$ (85b)

where
$\mathtt{D}_{\phantom{}ik}=\gamma_{ij}\mathtt{D}_{\phantom{j}A}^{j}h^{A}_{\phantom{A}k}$.
As already noted in [2], the gravitational part of the energy-momentum
$\sigma^{\mu}_{\ \,\nu}$ is trace-free:

${\sigma^{\mu}_{\ \,\mu}}=0.$ (86)

##  10 Hyperbolicity of the vacuum equations

Despite the structure of system (79) shares a lot of similarities with
equations of nonlinear electrodynamics of moving media which belong to the
class of symmetric hyperbolic thermodynamically compatible (SHTC) equations
[41, 25, 47, 48, 24], we cannot immediately apply the machinery of SHTC
systems to symmetrize the TEGR equations because some elements do not fit well
to the formalism of SHTC systems yet. Instead, we demonstrate the symmetric
hyperbolicity of the TEGR vacuum equations (without momentum (79c) and energy
(79d) equations) via their surprising equivalence to the tetrad reformulation
of GR by Estabrook-Robinson-Wahlquist [21] and Buchman-Bardeen [11] which is
known to be symmetric hyperbolic, and will be referenced to as ERWBB
formulation. We remind that solely hyperbolicity is not enough for well-
posedness of the initial value problem for a general quasi-linear system of
first-order equations. At least, to the best of our knowledge, there is no an
existence and uniqueness theorem for such a class of equations in multiple
dimensions666 Despite this, for a system describing causal propagation of
signals, the hyperbolicity is considered as a minimum requirement for
numerical discretization.. In contrary, such a theorem exists for symmetric
hyperbolic systems, e.g. see [7].

Let us first introduce the main elements of the ERWBB formulation which relies
on the Ricci rotation coefficients

$\displaystyle\mathcal{R}_{abc}=$
$\displaystyle\bm{\eta}_{a}\cdot\nabla_{c}\bm{\eta}_{b}$ $\displaystyle=$
$\displaystyle\eta^{\mu}_{\phantom{\mu}a}g_{\mu\nu}\eta^{\lambda}_{\phantom{\lambda}c}\partial_{\lambda}\eta^{\nu}_{\phantom{\nu}b}+\eta^{\mu}_{\phantom{\mu}a}g_{\mu\nu}\eta^{\lambda}_{\phantom{\lambda}c}\Gamma^{\nu}_{\
\lambda\rho}\eta^{\rho}_{\phantom{\rho}b}$ (87)

where $\nabla_{c}=\eta^{\lambda}_{\phantom{\lambda}c}\nabla_{\lambda}$ and
$\nabla_{\lambda}$ is the standard covariant derivative of GR associated with
the symmetric Levi-Civita connection, and $\Gamma^{\mu}_{\ \nu\lambda}$ are
the Christoffel symbols of the Levi-Civita connection.

If
$T^{a}_{\phantom{a}bc}=T^{a}_{\phantom{a}\mu\nu}\eta^{\mu}_{\phantom{\mu}b}\eta^{\nu}_{\phantom{\nu}c}$
and $T_{\phantom{}abc}=\kappa_{ad}T^{d}_{\phantom{d}bc}$ is the torsion
written in the tetrad basis, then the relation between the Ricci rotation
coefficient and torsion can be expressed by the formula

$\mathcal{R}_{abc}=\frac{1}{2}\left(T_{\phantom{}abc}+T_{\phantom{}bca}-T_{\phantom{}cab}\right),$
(88)

which also shows that the Ricci rotation coefficients equal exactly to the so-
called contortion tensor with the opposite sign $\mathcal{R}_{abc}=-K_{abc}$,
e.g. see [2, Eq.(1.63)].

Then, 24 independent entries of $\mathcal{R}_{abc}$ are organized into the
following state variables

$\mathtt{K}_{{C}{A}}:=\mathcal{R}_{A\hat{{0}}C},\qquad\mathtt{N}_{B}^{\
\,A}:=\frac{1}{2}\varepsilon^{ACD}\mathcal{R}_{CDB}$ (89)

and

$a_{A}=\mathcal{R}_{A\hat{{0}}\hat{{0}}},\qquad\omega^{A}=\frac{1}{2}\varepsilon^{ABC}\mathcal{R}_{CB\hat{{0}}}$
(90)

Here,
$\varepsilon_{ABC}=h_{3}\bm{\in}_{ijk}\eta^{i}_{\phantom{i}A}\eta^{j}_{\phantom{j}B}\eta^{k}_{\phantom{k}C}$.

It can be shown that the following relations between the ERWBB and TEGR state
variables hold

$\displaystyle\mathtt{K}_{{A}{B}}=$
$\displaystyle\frac{\varkappa}{h_{3}}\left(\mathtt{D}_{\phantom{}AB}-\frac{1}{2}\mathtt{D}_{\phantom{C}C}^{C}\kappa_{AB}\right),$
(91a) $\displaystyle-\mathtt{N}_{A}^{\ \,B}=$
$\displaystyle\frac{1}{h_{3}}\left(\mathtt{B}_{A}^{\phantom{A}B}-\frac{1}{2}\mathtt{B}^{C}_{\phantom{C}C}\delta^{B}_{\
A}\right).$ (91b)

We remark that due to our choice of the $3+1$ split (47), i.e. that observer’s
time vector, $\bm{\eta}_{\hat{{0}}}$, is align with the normal vector to the
spatial hypersurfaces, matrices $\mathtt{K}_{{A}{B}}$ and
$\mathtt{D}_{\phantom{}AB}$ are symmetric.

Additionally, we have the relations for $a_{A}$, $\omega^{A}$, and
$A_{k}=\mathtt{E}^{\hat{{0}}}_{\phantom{\hat{{0}}}k}=\partial_{k}\alpha$ and
$\Omega_{k}=\mathtt{H}_{k\hat{{0}}}$

$a_{A}=\alpha^{-1}A_{A},\qquad\omega^{A}=\varkappa\alpha^{-1}\kappa^{AB}\Omega_{B}.$
(92)

Finally, expression of the constitutive fluxes $\mathtt{E}^{A}_{\phantom{A}k}$
and $\mathtt{H}_{kA}$ in terms of
$\\{\mathtt{K}_{{A}{B}},\mathtt{N}_{AB},a_{A},\omega_{A}\\}$ read

$\displaystyle\mathtt{E}_{\phantom{}AB}$
$\displaystyle=-\alpha\,\mathtt{K}_{{A}{B}}+\alpha\varepsilon_{ABC}\omega^{C},$
(93a) $\displaystyle\mathtt{H}_{AB}$
$\displaystyle=\phantom{-}\frac{\alpha}{\varkappa}\mathtt{N}_{AB}-\frac{\alpha}{\varkappa}\varepsilon_{ABC}a^{C}.$
(93b)

Unfortunately, system (79) is not immediately hyperbolic even for the
Minkowski spacetime (it is actually only weakly hyperbolic because some
eigenvectors are missing), and therefore it can not be equivalent to the
symmetric hyperbolic ERWBB formulation of GR right away. In this section, we
discuss what modifications of equations (79) are necessary in order to make it
equivalent to the ERWBB formulation. These modifications rely on the use of
Hamiltonian and momentum constraints (81) and symmetry of
$\mathtt{D}_{\phantom{}AB}$.

Because we are considering a vacuum space-time, the momentum and energy
equations (79c), (79d) can be omitted. Moreover, modifications concern only
torsion equations (79a), (79b), and therefore for the sake of simplicity, we
also omit the tetrad equations (79e)

If one considers 18 equations (79a), (79b) written in a quasilinear form777Due
to the rotational invariance, it is actually sufficient to consider the system
in any preferred coordinate direction, say in $x^{1}$, without lose of
generality.

$\partial_{t}\bm{Q}+\bm{M}^{k}(\bm{Q})\partial_{k}\bm{Q}=\bm{S}$ (94)

for 18 unknowns $\bm{Q}=\\{\mathtt{D}_{\phantom{i}A}^{i},\mathtt{B}^{Ai}\\}$
then the matrix $\bm{M}=\xi_{k}\bm{M}^{k}$ has all real eigenvalues but only
10 linearly independent eigenvectors already for the flat Minkowski spacetime,
that is the system is only weakly hyperbolic. Here,
$\bm{M}^{k}=\partial\bm{F}^{k}/\partial\bm{Q}$ is the Jacobian of the $k$-th
flux vector.

Evolution equations (79a) and (79b) for the torsion fields have a conservative
form (i.e. their differential terms are the 4-divergence) which is in
particular convenient for developing of high-order methods for hyperbolic
equations. However, this flux-conservative form of equations has a flaw in the
eigenvalues already for the flat spacetime with non-zero shift vector
$\beta^{k}$. Namely, some characteristic velocities are always $0$ despite the
non-zero shift vector. The same issue is well known for the
magnetohydrodynamics equations [44]. From our experience with the non-
relativistic equations of nonlinear electrodynamics of moving media [20] and
continuum mechanics with torsion [43] that have the same structure as (79).
Instead, the following non-conservative form does not have such a pathology

$\displaystyle\partial_{t}\mathtt{D}_{\phantom{i}A}^{i}+\partial_{k}(\beta^{i}\mathtt{D}_{\phantom{k}A}^{k}-\beta^{k}\mathtt{D}_{\phantom{i}A}^{i}-\bm{\in}^{ikj}\,\mathtt{H}_{jA})-\beta^{i}\partial_{k}\mathtt{D}_{\phantom{k}A}^{k}$
$\displaystyle=-s_{\ A}^{i},$ (95a)
$\displaystyle\partial_{t}\mathtt{B}^{Ai}+\partial_{k}(\,\beta^{i}\mathtt{B}^{Ak}-\beta^{k}\mathtt{B}^{Ai}+\bm{\in}^{ikj}\mathtt{E}^{A}_{\phantom{A}j})-\beta^{i}\partial_{k}\mathtt{B}^{Ak}$
$\displaystyle=0,$ (95b)

and which is obtained from (79a) and (79b) by turning the momentum $\rho_{A}$
on the right hand-side of (79a) into a differential term according to the
momentum constraint $\rho_{A}=\partial_{k}\mathtt{D}_{\phantom{k}A}^{k}$ and
adding formally a zero, $0=\partial_{k}\mathtt{B}^{Ak}$, to (79b). These
modifications are legitimate for smooth solutions thanks to involution
constraints (81) and (75). After this, the characteristic velocities have the
desired values equal to $\sim\beta^{k}$ or $\sim\beta^{k}\pm c$, if $c$ is the
light speed. However, this only allows to fix the issue with the eigenvalues,
while some eigenvectors are still missing.

In order to recover all missing eigenvectors, we shall use the Hamiltonian
constraint (81)1 and the symmetry of $\mathtt{D}_{\phantom{}ij}$. Thus, we
shall turn the potential $\mathtt{U}=-\alpha\rho_{\hat{{0}}}$ in the
expression of the stress tensor $s_{\ A}^{i}$ into the differential term
$\mathtt{U}=-\alpha\partial_{k}\mathtt{D}_{\phantom{k}\hat{{0}}}^{k}$
according to the Hamiltonian constraint (81)1. Additionally, we use expression
(85a)2 for $\mathtt{D}_{\phantom{k}\hat{{0}}}^{k}$, i.e.
$\mathtt{D}_{\phantom{k}\hat{{0}}}^{k}=-\varkappa^{-1}h_{3}^{-1}\bm{\in}^{klj}\mathtt{B}_{\phantom{}lj}$.
Similarly, we are adding formally a zero
$0=\partial_{k}(h_{3}^{-1}\bm{\in}^{klj}\mathtt{D}_{\phantom{}lj})=-\mathtt{B}^{\hat{{0}}k}$
to the equation for the $\mathtt{B}^{Ai}$ field. This expression is zero due
to the fact that the contraction of antisymmetric Levi-Civita symbol with the
symmetric tensor $\mathtt{D}_{\phantom{}ij}$ is zero. The resulting equations
read

$\partial_{t}\mathtt{D}_{\phantom{i}A}^{i}+\partial_{k}(\beta^{i}\mathtt{D}_{\phantom{k}A}^{k}-\beta^{k}\mathtt{D}_{\phantom{i}A}^{i}-\bm{\in}^{ikj}\,\mathtt{H}_{jA})\\\
-\beta^{i}\partial_{k}\mathtt{D}_{\phantom{k}A}^{k}+\alpha\eta^{i}_{\phantom{i}A}\partial_{k}\mathtt{H}^{k}_{\phantom{k}}=-s_{\
A}^{i}-\mathtt{U}\eta^{i}_{\phantom{i}A},$ (96a)
$\partial_{t}\mathtt{B}^{Ai}+\partial_{k}(\,\beta^{i}\mathtt{B}^{Ak}-\beta^{k}\mathtt{B}^{Ai}+\bm{\in}^{ikj}\mathtt{E}^{A}_{\phantom{A}j})\\\
-\beta^{i}\partial_{k}\mathtt{B}^{Ak}-\alpha\eta^{i}_{\phantom{i}B}\partial_{k}\mathtt{E}^{k}_{\phantom{k}}\kappa^{AB}=0,$
(96b) where we have defined new vectors
$\displaystyle\mathtt{E}^{k}_{\phantom{k}}:=-\varkappa
h_{3}^{-1}\bm{\in}^{klj}\mathtt{D}_{\phantom{}lj},\qquad\mathtt{H}^{k}_{\phantom{k}}:=-\varkappa^{-1}h_{3}^{-1}\bm{\in}^{klj}\mathtt{B}_{\phantom{}lj}.$
(96c)

Now, using relations (91)–(93) and after a lengthy but rather straightforward
sequence of transformations, equations (96) can be written in the ERWBB form

$\displaystyle\partial_{\hat{{0}}}\mathtt{K}_{{A}{B}}$
$\displaystyle-\kappa_{AC}\varepsilon^{CDE}\partial_{D}\mathtt{N}_{EB}-\partial_{A}a_{B}=\text{l.o.t.}$
(97a) $\displaystyle\partial_{\hat{{0}}}\mathtt{N}_{AB}$
$\displaystyle+\kappa_{BC}\varepsilon^{CDE}\partial_{D}\mathtt{K}_{{E}{A}}+\partial_{A}\omega_{B}=\text{l.o.t.}$
(97b)

where ’l.o.t.’ stands for ’low-order terms’ (i.e. algebraic terms that do not
contain space and time derivatives),
$\partial_{\hat{{0}}}=\eta^{\mu}_{\phantom{\mu}\hat{{0}}}\partial_{\mu}=u^{\mu}\partial_{\mu}=\alpha^{-1}\partial_{t}-\alpha^{-1}\beta^{k}\partial_{k}$,
and $\partial_{A}=\eta^{k}_{\phantom{k}A}\partial_{k}$.

Similar to $A_{k}$ and $\Omega_{k}$ in TEGR, the vectors $a_{B}$ and
$\omega_{A}$ are not provided with particular evolution equations and
therefore, one could try to choose their evolution equations in such a way as
to guaranty the well-posedness of the enlarged system for the unknowns
$\\{\mathtt{K}_{{A}{B}},\mathtt{N}_{AB},a_{B},\omega_{B}\\}$. Thus, as shown
in [21, 11], if coupled with equations whose differential part is

$\displaystyle\partial_{\hat{{0}}}a_{A}-\partial_{B}\mathtt{K}^{B}_{\ A}$
$\displaystyle=0,$ (98)
$\displaystyle\partial_{\hat{{0}}}\omega_{A}+\partial_{B}\mathtt{N}^{B}_{\ A}$
$\displaystyle=0,$ (99)

system (97) is symmetric hyperbolic, and thus (96) can also be put into a
symmetric hyperbolic form. However, we can not immediately conclude the same
for the original $3+1$ TEGR equations (79) because the differential operators
of (79) and (96) are different even though the overall equations are
equivalent. We recall that the hyperbolicity (well-posedness of the initial
value problem) is defined by the leading-order terms of the differential
operator (all differential terms in the case of first-order systems).
Nevertheless, this is an important result demonstrating that inside of the
class of equations equivalent to TEGR equations (79) there is a possibility,
at least one given by (96), to rewrite the TEGR vacuum equations as a
symmetric hyperbolic system.

##  11 Relation between the torsion and extrinsic curvature

It is useful to relate the state variables of TEGR to the conventional
quantities used in numerical general relativity [4, 6, 27]. In what follows,
we relate the spatial extrinsic curvature of GR to the
$\mathtt{D}_{\phantom{k}A}^{k}$ field. We remark that the two fields are
conceptually different due to the different geometry interpretations in GR and
TEGR. Therefore, the obtained relation is possible only due to the equivalence
of TEGR and GR.

In GR, the evolution equation of the spatial metric $\gamma_{ij}$ is (see [45,
Eq.(7.64)])

$\partial_{t}\gamma_{ij}-\beta^{l}\partial_{l}\gamma_{ij}-\gamma_{il}\partial_{j}\beta^{l}-\gamma_{jl}\partial_{i}\beta^{l}=-2\alpha
K_{ij},$ (100)

where $K_{ij}$ is the spatial extrinsic curvature. On the other hand, in TEGR,
contracting (61) with $\kappa_{AB}h^{B}_{\phantom{B}j}$, one obtains the
following evolution equation for the spatial metric:

$\partial_{t}\gamma_{ij}-\beta^{l}\partial_{l}\gamma_{ij}-\gamma_{il}\partial_{j}\beta^{l}-\gamma_{jl}\partial_{i}\beta^{l}=\\\
-\kappa_{AB}(h^{A}_{\phantom{A}i}\mathtt{E}^{B}_{\phantom{B}j}+h^{A}_{\phantom{A}j}\mathtt{E}^{B}_{\phantom{B}i}).$
(101)

Therefore, one can deduce an expression for the extrinsic curvature in terms
of $\mathtt{E}^{A}_{\phantom{A}i}$:

$K_{ij}=\frac{1}{2\alpha}\kappa_{AB}(h^{A}_{\phantom{A}i}\mathtt{E}^{B}_{\phantom{B}j}+h^{A}_{\phantom{A}j}\mathtt{E}^{B}_{\phantom{B}i}).$
(102)

To obtain another expression for $K_{ij}$ in terms of the primary state
variable $\mathtt{D}_{\phantom{i}A}^{i}$, one needs to use the the
constitutive relationship (80b) to deduce

$K_{ij}=-\frac{\varkappa}{h_{3}}\left(\mathtt{D}_{\phantom{}ij}-\frac{1}{2}\mathtt{D}_{\phantom{k}k}^{k}\gamma_{ij}\right),$
(103)

which, if contracted, gives the relationship for the traces $K^{i}_{\
i}=\gamma^{ij}K_{ji}$ and
$\mathtt{D}_{\phantom{i}i}^{i}=\mathtt{D}_{\phantom{i}A}^{i}h^{A}_{\phantom{A}i}$

$K^{i}_{\ i}=\frac{\varkappa}{2h_{3}}\mathtt{D}_{\phantom{i}i}^{i}.$ (104)

Remark that if written in the tetrad frame
$K_{AB}=\eta^{i}_{\phantom{i}A}\eta^{j}_{\phantom{j}A}K_{ij}$ then the
extrinsic curvature is exactly $\mathtt{K}_{{A}{B}}$ introduced in (91) apart
from the opposite sign

$K_{AB}=-\mathtt{K}_{{A}{B}}.$ (105)

##  12 Torsion scalar

In TEGR and its $f(\mathcal{T})$-extensions, the Lagrangian density is a
function of the torsion scalar $\mathcal{T}$, e.g. in TEGR, the Lagrangian
density is

$\mathfrak{L}(h^{a}_{\phantom{a}\mu},T^{a}_{\phantom{a}\mu\nu})=\frac{h}{2\,\varkappa}\mathcal{T},$
(106a)
$\displaystyle\mathcal{T}(h^{a}_{\phantom{a}\mu},T^{a}_{\phantom{a}\mu\nu})$
$\displaystyle:=$ (106b)
$\displaystyle\frac{1}{4}g^{\beta\lambda}g^{\mu\gamma}g_{\alpha\eta}T^{\alpha}_{\phantom{\alpha}\lambda\gamma}T^{\eta}_{\phantom{\eta}\beta\mu}$
(106c)
$\displaystyle+\frac{1}{2}g^{\mu\gamma}T^{\lambda}_{\phantom{\lambda}\mu\beta}T^{\beta}_{\phantom{\beta}\gamma\lambda}$
(106d)
$\displaystyle-g^{\mu\lambda}T^{\rho}_{\phantom{\rho}\mu\rho}T^{\gamma}_{\phantom{\gamma}\lambda\gamma},$
(106e)

where $\varkappa=8\pi Gc^{-4}$ is the Einstein gravitational constant, $G$ is
the gravitational constant, $c$ is the speed of light in vacuum.

However, to close system (79), we need not the Lagrangian
$\mathfrak{L}(h^{a}_{\phantom{a}\mu},T^{a}_{\phantom{a}\mu\nu})$ directly but
we need to perform a sequence of variable and potential changes:
$\mathfrak{L}(h^{a}_{\phantom{a}\mu},T^{a}_{\phantom{a}\mu\nu})=\mathcal{L}(h^{a}_{\phantom{a}\mu},E^{a}_{\phantom{a}\mu},B^{a\mu})$
$\longrightarrow$
$E^{a}_{\phantom{a}\mu}D_{\phantom{\mu}a}^{\mu}-\mathcal{L}=U(h^{a}_{\phantom{a}\mu},D_{\phantom{\mu}a}^{\mu},B^{a\mu})=\mathcal{U}(h^{a}_{\phantom{a}\mu},\mathcal{D}_{\phantom{\mu}a}^{\mu},\mathcal{B}^{a\mu})$.
Thus, we have

$\displaystyle\mathcal{L}(h^{a}_{\phantom{a}\mu},E^{a}_{\phantom{a}\mu},B^{a\mu})=$
$\displaystyle\frac{h}{2\,\varkappa}\bigg{(}-\frac{1}{2}(E^{\hat{{0}}\lambda}E^{\hat{{0}}}_{\
\lambda}-2E^{\lambda}_{\ \,\lambda}E^{\beta}_{\ \,\beta}+E^{\
\,\lambda}_{\beta}E^{\beta}_{\ \,\lambda}+E^{\lambda}_{\ \,\beta}E^{\beta}_{\
\,\lambda})$
$\displaystyle+\bm{\in}_{\lambda\gamma\eta\rho}u^{\eta}(E^{\lambda\gamma}B^{\hat{{0}}\rho}+2E^{\hat{{0}}\lambda}_{\phantom{\hat{{0}}\lambda}}B^{\gamma\rho})$
$\displaystyle-\frac{1}{2}h^{-2}(B^{\hat{{0}}\lambda}B^{\hat{{0}}}_{\
\lambda}+B^{\lambda}_{\ \,\lambda}B^{\beta}_{\ \,\beta}-2B^{\lambda}_{\
\,\beta}B^{\beta}_{\ \,\lambda})\bigg{)}.$ (107)

In turn, if we perform the Legendre transform
$U(h^{A}_{\phantom{A}k},D_{\phantom{k}A}^{k},B^{Ak})=E^{a}_{\phantom{a}\mu}D_{\phantom{\mu}a}^{\mu}-\mathcal{L}$
then the new potential $U$ reads

$\displaystyle U(h^{A}_{\phantom{A}k},D_{\phantom{k}A}^{k},B^{Ak})$
$\displaystyle=$ $\displaystyle\frac{1}{4\,h}\bigg{(}$
$\displaystyle\varkappa\left(D_{\phantom{i}i}^{i}D_{\phantom{k}k}^{k}-2D_{\phantom{k}i}^{k}D_{\phantom{i}k}^{i}\right)$
$\displaystyle+$ $\displaystyle\frac{1}{\varkappa}\left(B^{i}_{\ \,i}B^{k}_{\
\,k}-2B^{i}_{\ \,k}B^{k}_{\ \,i}\right)$
$\displaystyle+\varkappa(\beta^{j}D_{\phantom{0}j}^{0}+$ $\displaystyle
2\beta^{j}D_{\phantom{0}j}^{0}\mathtt{D}_{\phantom{k}k}^{k}-4\beta^{j}D_{\phantom{0}k}^{0}\mathtt{D}_{\phantom{k}j}^{k})\bigg{)},$
(108)

where $D_{\phantom{i}k}^{i}=D_{\phantom{i}A}^{i}h^{A}_{\phantom{A}k}$ and
$B^{k}_{\ i}=B^{Aj}\eta^{k}_{\phantom{k}A}\gamma_{ji}$, and in the last two
terms, one should pay attention that the new field
$\mathtt{D}_{\phantom{k}A}^{k}$ introduced in (55) appears there.

Apparently, apart from the last terms in (12) depending on
$D_{\phantom{0}k}^{0}$, the potential $U$ is more symmetric in the variables
$D_{\phantom{\mu}a}^{\mu}$ and $B^{a\mu}$ rather than the Lagrangian
$\mathcal{L}$ in $E^{a}_{\phantom{a}\mu}$ and $B^{a\mu}$. This in particular,
justifies the introduction of the new and final variables (55)

$\mathtt{D}_{\phantom{k}A}^{k}=D_{\phantom{k}A}^{k}+\beta^{k}D_{\phantom{0}A}^{0},\qquad\mathtt{B}^{Ak}=B^{Ak},$
(109)

so that the potential
$U(h^{A}_{\phantom{A}k},D_{\phantom{k}A}^{k},B^{Ak})=\mathtt{U}(h^{A}_{\phantom{A}k},\mathtt{D}_{\phantom{k}A}^{k},\mathtt{B}^{Ak})$
becomes just

$\displaystyle\mathtt{U}(h^{A}_{\phantom{A}k},\mathtt{D}_{\phantom{k}A}^{k},\mathtt{B}^{Ak})$
$\displaystyle=$ $\displaystyle-\frac{\alpha}{2h_{3}}\bigg{(}$
$\displaystyle\varkappa\left(\mathtt{D}_{\phantom{A}B}^{A}\mathtt{D}_{\phantom{B}A}^{B}-\frac{1}{2}(\mathtt{D}_{\phantom{A}A}^{A})^{2}\right)$
$\displaystyle+$
$\displaystyle\frac{1}{\varkappa}\left(\mathtt{B}^{A}_{\phantom{A}B}\mathtt{B}^{B}_{\phantom{B}A}-\frac{1}{2}(\mathtt{B}^{A}_{\phantom{A}A})^{2}\right)\bigg{)}.$
(110)

##  13 Conclusion and discussion

We have presented a $3+1$-split of the TEGR equations in their historical pure
tetrad version, i.e. with the spin connection set to zero. To the best of our
knowledge there were no attempt to obtain a $3+1$-split of the TEGR equations
and therefore no attempts to solve numerically the full TEGR system of
equations for general spacetimes. This paper, therefore, may help to cover
this gap.

The derivation was done for a general Lagrangian and therefore extensions of
TEGR such as $f(\mathcal{T})$-teleparallel theories also can be covered in the
future.

The $3+1$ governing partial differential equations have appeared to have the
same structure as equations of nonlinear electrodynamics [20] and equations
for continuum fluid and solid mechanics with torsion [43]. Moreover, it has
appeared that the derived equations are equivalent to the recently proposed
dGREM tetrad formulation of GR [39].

The derived $3+1$ TEGR equations are not immediately hyperbolic as usually the
case for many first-order reductions of the Einstein equations. We
demonstrated that the differential operator of the equations can be
transformed into a different but equivalent form which, in turn and
surprisingly is equivalent to the symmetric hyperbolic tetrad reformulation of
GR by Estabrook-Robinson-Wahlquist [21] and Buchman-Bardeen [11]. The main
idea of these transformations relies on the use of multiple stationary
identities of the theory, e.g. Hamiltonian and momentum constraints. It,
therefore, should be very critical to fulfill these constraints on the
discrete level.

Despite it is argued that TEGR is fully equivalent to Einstein’s general
relativity, the proposed $3+1$ TEGR equations have yet to be carefully tested
in a computational code and have yet to be proved to pass all the standard
benchmark tests of GR. Therefore, further research will concern implementation
of the TEGR equations in a high-order discontinuous Galerkin code [17, 12],
with a possibility of constraint cleaning [16] and well-balancing [23]. This
in particular would allow a direct comparison of the TEGR with other $3+1$
equations of GR, such as Z4 formulations [3] forwarded by Bona et al in [9],
and FO-CCZ4 by Dumbser et al [18], a strongly hyperbolic formulations of GR,
within the same computational code. Another numerical strategy would be to use
staggered grids and to develop a structure-preserving discretization [10, 39,
22] that should allow to keep errors of div and curl-type involution
constraints of TEGR at the machine precision.

##### Acknowledgments

I.P. is grateful to A. Golovnev for sharing his thoughts about various aspects
of teleparallel gravity that were important for the appearance of this work.
Also, I.P. would like to thank M. Dumbser, O. Zanotti, E.R. Most for their
support and insightful discussions. I.P. is a member of the Gruppo Nazionale
per il Calcolo Scientifico of the Istituto Nazionale di Alta Matematica (INdAM
GNCS) and acknowledges the financial support received from the Italian
Ministry of Education, University and Research (MIUR) in the frame of the
Departments of Excellence Initiative 2018–2022 attributed to the Department of
Civil, Environmental and Mechanical Engineering (DICAM) of the University of
Trento (Grant No. L.232/2016) and in the frame of the Progetti di Rilevante
Interesse Nazionale (PRIN) 2017, Project No. 2017KKJP4X, “Innovative numerical
methods for evolutionary partial differential equations and applications”.
Funding for H.O. comes from Radboud University Nijmegen through a Virtual
Institute of Accretion (VIA) postdoctoral fellowship from the Netherlands
Research School for Astronomy (NOVA). The work of E.R. is supported by
Mathematical Center in Akademgorodok, the agreement with Ministry of Science
and High Education of the Russian Federation number 075-15-2022-281.

##  Appendix A Transformation of the Noether current

Here, we express Noether’s current
$J^{\mu}_{\phantom{\mu}a}=L_{h^{a}_{\phantom{a}\mu}}$ for the gravitational
part of the Lagrangian (i.e. the matter part is ignored in this section) in
terms of the potential $U$ and new variables $D_{\phantom{\mu}a}^{\mu}$ and
$B^{a\mu}$.

Thus, for the parametrization
$L(h^{a}_{\phantom{a}\mu},\accentset{\star}{T}^{a\mu\nu})=\mathcal{L}(h^{a}_{\phantom{a}\mu},B^{a\mu},E^{a}_{\phantom{a}\mu})$,
one has

$L_{h^{a}_{\phantom{a}\mu}}=\mathcal{L}_{h^{a}_{\phantom{a}\mu}}+\mathcal{L}_{B^{b\lambda}}\frac{\partial
B^{b\lambda}}{\partial
h^{a}_{\phantom{a}\mu}}+\mathcal{L}_{E^{b}_{\phantom{b}\lambda}}\frac{\partial
E^{b}_{\phantom{b}\lambda}}{\partial h^{a}_{\phantom{a}\mu}}.$ (111)

Then, using the definitions of the frame 4-velocity
$u^{\nu}=\eta^{\nu}_{\phantom{\nu}\hat{{0}}}$ and
$u_{\nu}=-h^{\hat{{0}}}_{\phantom{\hat{{0}}}\nu}$ and the torsion fields
$B^{b\lambda}=\accentset{\star}{T}^{b\lambda\nu}u_{\nu}=-\accentset{\star}{T}^{b\lambda\nu}h^{\hat{{0}}}_{\phantom{\hat{{0}}}\nu}$
and
$E^{b}_{\phantom{b}\lambda}=T^{b}_{\phantom{b}\lambda\nu}u^{\nu}=-\frac{1}{2}\bm{\in}_{\lambda\nu\alpha\beta}\accentset{\star}{T}^{b\alpha\beta}\eta^{\nu}_{\phantom{\nu}\hat{{0}}}$,
and the fact that $\partial\eta^{\lambda}_{\phantom{\lambda}b}/\partial
h^{a}_{\phantom{a}\mu}=-\eta^{\lambda}_{\phantom{\lambda}a}\eta^{\mu}_{\phantom{\mu}b}$,
we can rewrite (111) as

$L_{h^{a}_{\phantom{a}\mu}}=\mathcal{L}_{h^{a}_{\phantom{a}\mu}}+\mathcal{L}_{B^{b\lambda}}\frac{\partial
B^{b\lambda}}{\partial
h^{a}_{\phantom{a}\mu}}+\mathcal{L}_{E^{b}_{\phantom{b}\lambda}}\frac{\partial
E^{b}_{\phantom{b}\lambda}}{\partial h^{a}_{\phantom{a}\mu}}=\\\
\mathcal{L}_{h^{a}_{\phantom{a}\mu}}-\mathcal{L}_{B^{b\lambda}}\delta^{\hat{{0}}}_{\
a}(u^{\lambda}B^{b\mu}-u^{\mu}B^{b\lambda}+\bm{\in}^{\lambda\mu\alpha\beta}u_{\alpha}E^{b}_{\phantom{b}\beta})\\\
-\mathcal{L}_{E^{b}_{\phantom{b}\lambda}}(u_{\lambda}E^{b}_{\phantom{b}\nu}-u_{\nu}E^{b}_{\phantom{b}\lambda}-\bm{\in}_{\lambda\nu\alpha\rho}u^{\alpha}B^{b\rho})\eta^{\nu}_{\phantom{\nu}a}u^{\mu}.$
(112)

Using the definitions (38) and (39), the latter can be rewritten as

$\displaystyle L_{h^{a}_{\phantom{a}\mu}}=-U_{h^{a}_{\phantom{a}\mu}}$
$\displaystyle-U_{B^{b\lambda}}\delta^{\hat{{0}}}_{\
a}(-u^{\lambda}B^{b\mu}+u^{\mu}B^{b\lambda}+\bm{\in}^{\lambda\mu\alpha\beta}u_{\alpha}U_{D_{\phantom{\beta}b}^{\beta}})$
$\displaystyle-
D_{\phantom{\lambda}b}^{\lambda}\eta^{\nu}_{\phantom{\nu}a}u^{\mu}(u_{\lambda}U_{D_{\phantom{\nu}b}^{\nu}}-u_{\nu}U_{D_{\phantom{\lambda}b}^{\lambda}}+\bm{\in}_{\lambda\nu\alpha\rho}u^{\alpha}B^{b\rho}),$
(113)

which, after some term rearrangements, reads

$\displaystyle L_{h^{a}_{\phantom{a}\mu}}=-U_{h^{a}_{\phantom{a}\mu}}$
$\displaystyle+\delta^{\hat{{0}}}_{\
a}\bigg{(}u^{\lambda}B^{b\mu}U_{B^{b\lambda}}-u^{\mu}B^{b\lambda}U_{B^{b\lambda}}$
$\displaystyle-u^{\mu}D_{\phantom{\lambda}b}^{\lambda}U_{D_{\phantom{\lambda}b}^{\lambda}}+\bm{\in}^{\mu\lambda\rho\sigma}u_{\rho}U_{B^{b\lambda}}U_{D_{\phantom{\sigma}b}^{\sigma}}\bigg{)}$
$\displaystyle-\eta^{\nu}_{\phantom{\nu}a}u^{\mu}\left(u_{\lambda}D_{\phantom{\lambda}b}^{\lambda}U_{D_{\phantom{\nu}b}^{\nu}}-\bm{\in}_{\nu\lambda\rho\sigma}u^{\rho}B^{b\sigma}D_{\phantom{\lambda}b}^{\lambda}\right),$
(114)

and exactly is (41).

This formula, in particular, can be used to get the following expression for
the energy-momentum $\sigma^{\mu}_{\
\,\nu}=h^{a}_{\phantom{a}\nu}J^{\mu}_{\phantom{\mu}a}$:

$\displaystyle\sigma^{\mu}_{\ \,\nu}=$
$\displaystyle-h^{a}_{\phantom{a}\nu}U_{h^{a}_{\phantom{a}\mu}}$
$\displaystyle-u^{\lambda}u_{\nu}B^{a\mu}U_{B^{a\lambda}}-u^{\mu}u_{\lambda}D_{\phantom{\lambda}a}^{\lambda}U_{D_{\phantom{\nu}a}^{\nu}}$
$\displaystyle+u^{\mu}u_{\nu}B^{a\lambda}U_{B^{a\lambda}}+u^{\mu}u_{\nu}D_{\phantom{\lambda}a}^{\lambda}U_{D_{\phantom{\lambda}a}^{\lambda}}$
$\displaystyle+\bm{\in}_{\nu\sigma\lambda\rho}u^{\mu}u^{\sigma}B^{a\lambda}D_{\phantom{\rho}a}^{\rho}$
$\displaystyle+\bm{\in}^{\mu\sigma\lambda\rho}u_{\nu}u_{\sigma}U_{B^{a\lambda}}U_{D_{\phantom{\rho}a}^{\rho}}.$
(115)

##  Appendix B Expression for the energy-momentum

In this section, we derive expression (7.3) for the gravitational part fof the
energy-momentum

$\sigma^{\mu}_{\
\,\nu}=2\accentset{\star}{T}^{a\lambda\mu}L_{\accentset{\star}{T}^{a\lambda\nu}}-(\accentset{\star}{T}^{a\lambda\rho}L_{\accentset{\star}{T}^{a\lambda\rho}}-L)\delta^{\mu}_{\
\nu}$ (116)

Because $\accentset{\star}{T}^{a\mu\nu}$ is antisymmetric tensor, to compute
the derivative $L_{\accentset{\star}{T}^{a\lambda\nu}}$ one needs to use its
parametrization via the Weitzenböck connection, which is not symmetric, i.e.
$\accentset{\star}{T}^{a\mu\nu}=\bm{\in}^{\mu\nu\rho\sigma}W^{a}_{\phantom{a}\rho\sigma}$.
Thus, for Lagrangians
$\Lambda(h^{a}_{\phantom{a}\mu},W^{a}_{\phantom{a}\mu\nu})=L(h^{a}_{\phantom{a}\mu},\accentset{\star}{T}^{a\mu\nu})$
one can write

$\Lambda_{W^{a}_{\phantom{a}\lambda\mu}}=\bm{\in}^{\lambda\mu\rho\sigma}L_{\accentset{\star}{T}^{a\rho\sigma}},$
(117)

or using the identity
$\bm{\in}_{\lambda\mu\alpha\beta}\bm{\in}^{\lambda\mu\rho\sigma}=-2(\delta^{\rho}_{\
\alpha}\delta^{\sigma}_{\ \beta}-\delta^{\rho}_{\ \beta}\delta^{\sigma}_{\
\alpha})$,

$\bm{\in}_{\alpha\beta\rho\sigma}\Lambda_{W^{a}_{\phantom{a}\lambda\mu}}=-4L_{\accentset{\star}{T}^{a\rho\sigma}}.$
(118)

On the other hand, using the definitions
$E^{a}_{\phantom{a}\mu}=(W^{a}_{\phantom{a}\mu\nu}-W^{a}_{\phantom{a}\nu\mu})u^{\nu}$,
$B^{a\mu}=\bm{\in}^{\mu\nu\rho\sigma}W^{a}_{\phantom{a}\rho\sigma}u_{\nu}$,
and the parametrization
$\Lambda(h^{a}_{\phantom{a}\mu},W^{a}_{\phantom{a}\mu\nu})=\mathcal{L}(h^{a}_{\phantom{a}\mu},B^{a\mu},E^{a}_{\phantom{a}\mu})$,
one can write

$\displaystyle\Lambda_{W^{b}_{\phantom{b}\lambda\gamma}}$
$\displaystyle=\mathcal{L}_{E^{a}_{\phantom{a}\mu}}\frac{\partial
E^{a}_{\phantom{a}\mu}}{\partial
W^{b}_{\phantom{b}\lambda\gamma}}+\mathcal{L}_{B^{a\mu}}\frac{\partial
B^{a\mu}}{\partial W^{b}_{\phantom{b}\lambda\gamma}}$
$\displaystyle=u^{\gamma}\mathcal{L}_{E^{b}_{\phantom{b}\lambda}}-u^{\lambda}\mathcal{L}_{E^{b}_{\phantom{b}\gamma}}-\bm{\in}^{\lambda\gamma\nu\mu}u_{\nu}\mathcal{L}_{B^{a\mu}},$
(119)

and hence, from (118) and (B), one can deduce

$\displaystyle L_{\accentset{\star}{T}^{a\lambda\nu}}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left(u_{\lambda}\mathcal{L}_{B^{a\nu}}-u_{\nu}\mathcal{L}_{B^{a\lambda}}-\bm{\in}_{\lambda\nu\rho\sigma}u^{\rho}\mathcal{L}_{E^{a}_{\phantom{a}\sigma}}\right).$
(120)

Then, after contracting the later equation with
$\accentset{\star}{T}^{a\lambda\mu}=u^{\lambda}B^{a\mu}-u^{\mu}B^{a\lambda}+\bm{\in}^{\lambda\mu\alpha\beta}u_{\alpha}E^{a}_{\phantom{a}\beta}$,
one obtains

$\displaystyle 2\accentset{\star}{T}^{a\lambda\mu}$ $\displaystyle
L_{\accentset{\star}{T}^{a\lambda\nu}}=$ $\displaystyle
B^{a\mu}\mathcal{L}_{B^{a\nu}}-E^{a}_{\phantom{a}\nu}\mathcal{L}_{E^{a}_{\phantom{a}\mu}}$
$\displaystyle+u^{\lambda}u_{\nu}B^{a\mu}\mathcal{L}_{B^{a\lambda}}-u^{\mu}u_{\nu}B^{a\lambda}\mathcal{L}_{B^{a\lambda}}$
$\displaystyle-\bm{\in}_{\nu\lambda\rho\sigma}u^{\mu}u^{\sigma}B^{a\lambda}\mathcal{L}_{E^{a}_{\phantom{a}\rho}}$
$\displaystyle-\bm{\in}^{\mu\lambda\rho\sigma}u_{\nu}u_{\rho}E^{a}_{\phantom{a}\sigma}\mathcal{L}_{B^{a\lambda}}$
$\displaystyle+(u^{\mu}u_{\nu}+\delta^{\mu}_{\
\nu})E^{a}_{\phantom{a}\lambda}\mathcal{L}_{E^{a}_{\phantom{a}\lambda}}-u^{\mu}u_{\lambda}E^{a}_{\phantom{a}\nu}\mathcal{L}_{E^{a}_{\phantom{a}\lambda}}.$
(121)

This can be used to demonstrate that the full contraction
$\accentset{\star}{T}^{a\lambda\rho}L_{\accentset{\star}{T}^{a\lambda\rho}}$
results in

$\accentset{\star}{T}^{a\lambda\rho}L_{\accentset{\star}{T}^{a\lambda\rho}}=B^{a\lambda}\mathcal{L}_{B^{a\lambda}}+E^{a}_{\phantom{a}\lambda}\mathcal{L}_{E^{a}_{\phantom{a}\lambda}}.$
(122)

Collecting together (B) and (122) and using the change of variables and
potential (37)–(39), we arrive at

$\displaystyle 2\accentset{\star}{T}^{a\lambda\mu}$ $\displaystyle
L_{\accentset{\star}{T}^{a\lambda\nu}}-(\accentset{\star}{T}^{a\lambda\rho}L_{\accentset{\star}{T}^{a\lambda\rho}}-L)\delta^{\mu}_{\
\nu}=$
$\displaystyle-B^{a\mu}U_{B^{a\nu}}-D_{\phantom{\nu}a}^{\nu}U_{D_{\phantom{\mu}a}^{\mu}}$
$\displaystyle-u^{\lambda}u_{\nu}B^{a\mu}U_{B^{a\lambda}}-u^{\mu}u_{\lambda}D_{\phantom{\lambda}a}^{\lambda}U_{D_{\phantom{\nu}a}^{\nu}}$
$\displaystyle+u^{\mu}u_{\nu}B^{a\lambda}U_{B^{a\lambda}}+u^{\mu}u_{\nu}D_{\phantom{\lambda}a}^{\lambda}U_{D_{\phantom{\lambda}a}^{\lambda}}$
$\displaystyle+\bm{\in}_{\nu\sigma\lambda\rho}u^{\mu}u^{\sigma}B^{a\lambda}D_{\phantom{\rho}a}^{\rho}$
$\displaystyle+\bm{\in}^{\mu\sigma\lambda\rho}u_{\nu}u_{\sigma}U_{B^{a\lambda}}U_{D_{\phantom{\rho}a}^{\rho}}$
$\displaystyle+(B^{a\lambda}U_{B^{a\lambda}}+D_{\phantom{\lambda}a}^{\lambda}U_{D_{\phantom{\lambda}a}^{\lambda}}-U)\delta^{\mu}_{\
\nu}.$ (123)

##  Appendix C Transformation of the torsion PDE

In this appendix, we demonstrate how the Euler-Lagrange equation (35a)

$\partial_{\nu}(\bm{\in}^{\mu\nu\lambda\rho}L_{\accentset{\star}{T}^{a\lambda\rho}})=-L_{h^{a}_{\phantom{a}\mu}}$
(124)

can be transformed to the form (36a).

Based on the different parametrization of the Lagrangian (we omit for the
moment dependence of the Lagrangian on the tetrad field)
$\Lambda(W^{a}_{\phantom{a}\mu})=L(\accentset{\star}{T}^{a\mu\nu})=\mathfrak{L}(T^{a}_{\phantom{a}\mu\nu})$,
one can obtain

$\Lambda_{W^{a}_{\phantom{a}\lambda\mu}}=\bm{\in}^{\lambda\mu\rho\sigma}L_{\accentset{\star}{T}^{a\rho\sigma}},\quad\Lambda_{W^{a}_{\phantom{a}\lambda\mu}}=\mathfrak{L}_{T^{a}_{\phantom{a}\lambda\mu}}-\mathfrak{L}_{T^{a}_{\phantom{a}\mu\lambda}},$
(125)

from which it follows that the objects

$\accentset{\star}{\mathbb{T}}_{a\mu\nu}:=L_{\accentset{\star}{T}^{a\mu\nu}},\qquad\mathbb{T}_{a}^{\phantom{a}\mu\nu}:=\frac{1}{2}\left(\mathfrak{L}_{T^{a}_{\phantom{a}\mu\nu}}-\mathfrak{L}_{T^{a}_{\phantom{a}\nu\mu}}\right)$
(126)

are Hodge duals of each others:

$\displaystyle\accentset{\star}{\mathbb{T}}_{a\mu\nu}$
$\displaystyle:=-\frac{1}{2}\bm{\in}_{\mu\nu\rho\sigma}\mathbb{T}_{a}^{\phantom{a}\rho\sigma},$
(127a) $\displaystyle\mathbb{T}_{a}^{\phantom{a}\mu\nu}$
$\displaystyle:=\phantom{+}\frac{1}{2}\bm{\in}^{\mu\nu\rho\sigma}\accentset{\star}{\mathbb{T}}_{a\rho\sigma}.$
(127b)

Therefore, one can write the following identities

$\displaystyle\accentset{\star}{\mathbb{T}}_{a\mu\nu}$
$\displaystyle=u_{\mu}\mathbb{H}_{a\nu}-u_{\nu}\mathbb{H}_{a\mu}-\bm{\in}_{\mu\nu\rho\sigma}u^{\rho}\mathbb{D}_{a}^{\phantom{a}\sigma},$
(128a) $\displaystyle\mathbb{T}_{a}^{\phantom{a}\mu\nu}$
$\displaystyle=u^{\mu}\mathbb{D}_{a}^{\phantom{a}\nu}-u^{\nu}\mathbb{D}_{a}^{\phantom{a}\mu}+\bm{\in}^{\mu\nu\rho\sigma}u_{\rho}\mathbb{H}_{a\sigma},$
(128b)

where

$\mathbb{H}_{a\mu}:=\accentset{\star}{\mathbb{T}}_{a\mu\nu}u^{\nu},\qquad\mathbb{D}_{a}^{\phantom{a}\mu}:=\mathbb{T}_{a}^{\phantom{a}\mu\nu}u_{\nu},$
(129)

Hence, the Euler-Lagrange equation (124) now reads

$\partial_{\nu}\mathbb{T}_{a}^{\phantom{a}\mu\nu}=-\frac{1}{2}L_{h^{a}_{\phantom{a}\mu}},$
(130)

or, according to (128b), it can be written as

$\partial_{\nu}(u^{\mu}\mathbb{D}_{a}^{\phantom{a}\nu}-u^{\nu}\mathbb{D}_{a}^{\phantom{a}\mu}+\bm{\in}^{\mu\nu\alpha\beta}u_{\alpha}\mathbb{H}_{a\beta})=-\frac{1}{2}L_{h^{a}_{\phantom{a}\mu}}.$
(131)

It remains to express $\mathbb{D}_{a}^{\phantom{a}\mu}$ and
$\mathbb{H}_{a\mu}$ in terms of
$\mathcal{L}(B^{a\mu},E^{a}_{\phantom{a}\mu})$. Thus, using the fact that
$\Lambda_{W^{a}_{\phantom{a}\mu\nu}}=2\mathbb{T}_{a}^{\phantom{a}\mu\nu}$ and
$\Lambda_{W^{a}_{\phantom{a}\mu\nu}}=\bm{\in}_{\mu\nu\rho\sigma}\accentset{\star}{\mathbb{T}}_{a\rho\sigma}$,
and the expression (B), one can derive that

$\displaystyle\mathbb{H}_{a\mu}$
$\displaystyle=-\frac{1}{2}\left(\mathcal{L}_{B^{a\mu}}+u^{\lambda}\mathcal{L}_{B^{a\lambda}}u_{\mu}\right),$
(132a) $\displaystyle\mathbb{D}_{a}^{\phantom{a}\mu}$
$\displaystyle=-\frac{1}{2}\left(\mathcal{L}_{E^{a}_{\phantom{a}\mu}}+u_{\lambda}\mathcal{L}_{E^{a}_{\phantom{a}\lambda}}u^{\mu}\right).$
(132b)

Finally, plugging these expressions in (131), one obtains the desired result

$\partial_{\nu}(u^{\mu}\mathcal{L}_{E^{a}_{\phantom{a}\nu}}-u^{\nu}\mathcal{L}_{E^{a}_{\phantom{a}\mu}}+\bm{\in}^{\mu\nu\rho\sigma}u_{\rho}\mathcal{L}_{B^{a\sigma}})=L_{h^{a}_{\phantom{a}\mu}}.$
(133)

##  Appendix D Some expressions of the torsion scalar

Denoting the scalars in the right-hand side of (106b) as
($\mathcal{T}=\mathcal{T}_{1}+\mathcal{T}_{2}+\mathcal{T}_{3}$)

$\displaystyle\mathcal{T}_{1}=$
$\displaystyle\frac{1}{4}\kappa_{ab}g^{\beta\lambda}g^{\mu\gamma}T^{a}_{\phantom{a}\lambda\gamma}T^{b}_{\phantom{b}\beta\mu},$
(134) $\displaystyle\mathcal{T}_{2}=$
$\displaystyle\frac{1}{2}g^{\mu\gamma}\eta^{\lambda}_{\phantom{\lambda}a}\eta^{\beta}_{\phantom{\beta}b}T^{a}_{\phantom{a}\mu\beta}T^{b}_{\phantom{b}\gamma\lambda},$
(135) $\displaystyle\quad\mathcal{T}_{3}=$
$\displaystyle-g^{\mu\lambda}\eta^{\rho}_{\phantom{\rho}a}\eta^{\gamma}_{\phantom{\gamma}b}T^{a}_{\phantom{a}\mu\rho}T^{b}_{\phantom{b}\lambda\gamma},$
(136)

we can also write

$\displaystyle\mathcal{T}_{1}$ $\displaystyle=Q_{3}+\frac{1}{2}C_{4},$ (137)
$\displaystyle\mathcal{T}_{2}$
$\displaystyle=-Q_{1}+Q_{4}+L_{1}-C_{1}+C_{2}+\frac{1}{2}C_{4},$ (138)
$\displaystyle\mathcal{T}_{3}$
$\displaystyle=2Q_{1}+Q_{2}+L_{2}+2C_{1}+C_{3}-C_{4},$ (139)

where the scalars $Q$, $L$, and $C$ are scalars made of
$E^{a}_{\phantom{a}\mu}$ and $B^{a\mu}$.

Quadratic in $E^{a}_{\phantom{a}\mu}$:

$\displaystyle Q_{1}$
$\displaystyle=-\frac{1}{2}E^{\hat{{0}}}_{\phantom{\hat{{0}}}\alpha}g^{\alpha\beta}E^{\hat{{0}}}_{\phantom{\hat{{0}}}\beta},$
(140) $\displaystyle Q_{2}$
$\displaystyle=\eta^{\alpha}_{\phantom{\alpha}a}E^{a}_{\phantom{a}\alpha}\eta^{\beta}_{\phantom{\beta}b}E^{b}_{\phantom{b}\beta}=E^{\alpha}_{\phantom{\alpha}\alpha}E^{\beta}_{\phantom{\beta}\beta},$
(141) $\displaystyle Q_{3}$
$\displaystyle=-\frac{1}{2}\kappa_{ab}E^{a}_{\phantom{a}\alpha}g^{\alpha\mu}E^{b}_{\phantom{b}\mu},$
(142) $\displaystyle Q_{4}$
$\displaystyle=-\frac{1}{2}\eta^{\lambda}_{\phantom{\lambda}a}E^{a}_{\phantom{a}\beta}\eta^{\beta}_{\phantom{\beta}b}E^{b}_{\phantom{b}\lambda}=E^{\lambda}_{\phantom{\lambda}\beta}E^{\beta}_{\phantom{\beta}\lambda}.$
(143)

Mixed scalars (linear in $E^{a}_{\phantom{a}\mu}$):

$\displaystyle L_{1}$
$\displaystyle=\bm{\in}_{\lambda\tau\varphi\gamma}u^{\varphi}g^{\tau\beta}\eta^{\lambda}_{\phantom{\lambda}a}E^{a}_{\phantom{a}\beta}B^{\hat{{0}}\gamma},$
(144) $\displaystyle L_{2}$
$\displaystyle=2\bm{\in}_{\lambda\tau\varphi\gamma}u^{\varphi}g^{\lambda\beta}\eta^{\tau}_{\phantom{\tau}a}E^{\hat{{0}}}_{\phantom{\hat{{0}}}\beta}B^{a\gamma}.$
(145)

Quadratic in $B^{a\mu}$ (constant in $E^{a}_{\phantom{a}\mu}$)

$\displaystyle C_{1}$
$\displaystyle=-\frac{1}{2}h^{-2}g_{\lambda\rho}B^{\hat{{0}}\rho}B^{\hat{{0}}\lambda},$
(146) $\displaystyle C_{2}$
$\displaystyle=-\frac{1}{2}h^{-2}g_{\rho\varphi}\eta^{\varphi}_{\phantom{\varphi}a}B^{a\rho}\eta^{\beta}_{\phantom{\beta}b}B^{b\lambda}g_{\lambda\beta},$
(147) $\displaystyle C_{3}$
$\displaystyle=h^{-2}\eta^{\varphi}_{\phantom{\varphi}a}B^{a\sigma}g_{\varphi\beta}\eta^{\lambda}_{\phantom{\lambda}b}B^{b\beta}g_{\lambda\sigma},$
(148) $\displaystyle C_{4}$
$\displaystyle=h^{-2}\eta^{\varphi}_{\phantom{\varphi}a}B^{a\sigma}g_{\varphi\lambda}\eta^{\lambda}_{\phantom{\lambda}b}B^{b\beta}g_{\beta\sigma}.$
(149)

In terms of the Hodge dual $\accentset{\star}{T}^{a\mu\nu}$, the torsion
scalar $\mathcal{T}$ can be rewritten as

$\mathcal{T}=\mathcal{H}_{1}+\mathcal{H}_{2}+\mathcal{H}_{3}+\mathcal{H}_{4},$
(150)

where

$\displaystyle\mathcal{H}_{1}$
$\displaystyle=\frac{1}{2h^{2}}\kappa_{ac}\kappa_{bd}h^{c}_{\phantom{c}\lambda}h^{d}_{\phantom{d}\sigma}g_{\tau\rho}\accentset{\star}{T}^{a\rho\sigma}\accentset{\star}{T}^{b\lambda\tau},$
(151a) $\displaystyle\mathcal{H}_{2}$
$\displaystyle=\frac{1}{2h^{2}}\kappa_{ac}\kappa_{bd}h^{c}_{\phantom{c}\tau}h^{d}_{\phantom{d}\rho}g_{\lambda\sigma}\accentset{\star}{T}^{a\rho\sigma}\accentset{\star}{T}^{b\lambda\tau},$
(151b) $\displaystyle\mathcal{H}_{3}$
$\displaystyle=-\frac{1}{4h^{2}}\kappa_{ac}\kappa_{bd}h^{c}_{\phantom{c}\sigma}h^{d}_{\phantom{d}\lambda}g_{\tau\rho}\accentset{\star}{T}^{a\rho\sigma}\accentset{\star}{T}^{b\lambda\tau},$
(151c) $\displaystyle\mathcal{H}_{4}$
$\displaystyle=-\frac{1}{4h^{2}}\kappa_{ac}\kappa_{bd}h^{c}_{\phantom{c}\rho}h^{d}_{\phantom{d}\tau}g_{\lambda\sigma}\accentset{\star}{T}^{a\rho\sigma}\accentset{\star}{T}^{b\lambda\tau}.$
(151d)

Finally, in terms of the Weitzenböck connection, the torsion scalar reads:

$\mathcal{T}=\sum_{n=1}^{8}\mathcal{W}_{n}$ (152a)
$\displaystyle\mathcal{W}_{1}$
$\displaystyle=\frac{1}{2}\kappa_{ab}g^{\beta\lambda}g^{\mu\gamma}W^{a}_{\phantom{a}\lambda\gamma}W^{b}_{\phantom{b}\beta\mu},$
(152b) $\displaystyle\mathcal{W}_{2}$
$\displaystyle=-\frac{1}{2}\kappa_{ab}g^{\beta\lambda}g^{\mu\gamma}W^{a}_{\phantom{a}\gamma\lambda}W^{b}_{\phantom{b}\beta\mu},$
(152c) $\displaystyle\mathcal{W}_{3}$
$\displaystyle=\frac{1}{2}g^{\mu\gamma}\eta^{\beta}_{\phantom{\beta}a}W^{a}_{\phantom{a}\lambda\gamma}\eta^{\lambda}_{\phantom{\lambda}b}W^{b}_{\phantom{b}\beta\mu},$
(152d) $\displaystyle\mathcal{W}_{4}$
$\displaystyle=\frac{1}{2}g^{\mu\gamma}\eta^{\beta}_{\phantom{\beta}a}W^{a}_{\phantom{a}\gamma\lambda}\eta^{\lambda}_{\phantom{\lambda}b}W^{b}_{\phantom{b}\mu\beta},$
(152e)
$\mathcal{W}_{5}=-g^{\mu\gamma}\eta^{\beta}_{\phantom{\beta}a}W^{a}_{\phantom{a}\gamma\lambda}\eta^{\lambda}_{\phantom{\lambda}b}W^{b}_{\phantom{b}\beta\mu},$
(152f) $\displaystyle\mathcal{W}_{6}$
$\displaystyle=-g^{\mu\lambda}w^{1}_{\lambda}w^{1}_{\mu},$ (152g)
$\displaystyle\mathcal{W}_{7}$
$\displaystyle=-g^{\mu\lambda}w^{2}_{\lambda}w^{2}_{\mu},$ (152h)
$\displaystyle\mathcal{W}_{8}$
$\displaystyle=2g^{\mu\lambda}w^{1}_{\mu}w^{2}_{\lambda},$ (152i) where
$w^{1}_{\mu}=\eta^{\lambda}_{\phantom{\lambda}a}W^{a}_{\phantom{a}\mu\lambda}$,
$w^{2}_{\mu}=\eta^{\lambda}_{\phantom{\lambda}a}W^{a}_{\phantom{a}\lambda\mu}$.

## References

  * [1] Miguel Alcubierre “Introduction to 3+1 Numerical Relativity” Oxford: Oxford University Press, 2008 DOI: 10.1093/acprof:oso/9780199205677.001.0001
  * [2] Ruben Aldrovandi and José Geraldo Pereira “Teleparallel Gravity” Dordrecht: Springer Netherlands, 2013 DOI: 10.1007/978-94-007-5143-9
  * [3] Daniela Alic, Carles Bona-Casas, Carles Bona, Luciano Rezzolla and Carlos Palenzuela “Conformal and covariant formulation of the Z4 system with constraint-violation damping” In _Physical Review D - Particles, Fields, Gravitation and Cosmology_ 85.6, 2012, pp. 1–12 DOI: 10.1103/PhysRevD.85.064040
  * [4] Richard Arnowitt, Stanley Deser and Charles W. Misner “Republication of: The dynamics of general relativity” In _General Relativity and Gravitation_ 40.9, 2008, pp. 1997–2027 DOI: 10.1007/s10714-008-0661-1
  * [5] Sebastian Bahamonde and Jackson Levi Said “Teleparallel Gravity: Foundations and Observational Constraints—Editorial” In _Universe_ 7.8, 2021, pp. 269 DOI: 10.3390/universe7080269
  * [6] T Baumgarte “Numerical relativity and compact binaries” In _Physics Reports_ 376.2, 2003, pp. 41–131 DOI: 10.1016/S0370-1573(02)00537-9
  * [7] Sylvie Benzoni-Gavage and Denis Serre “Multi-dimensional hyperbolic partial differential equations” 325, Grundlehren der mathematischen Wissenschaften Berlin, Heidelberg: Oxford University Press, 2006 DOI: 10.1093/acprof:oso/9780199211234.001.0001
  * [8] Christian G. Böhmer and Yongjo Lee “Compatibility conditions of continua using Riemann–Cartan geometry” In _Mathematics and Mechanics of Solids_ 26.4, 2021, pp. 513–529 DOI: 10.1177/1081286520961453
  * [9] C. Bona, T. Ledvinka, C. Palenzuela and M. Žáček “General-covariant evolution formalism for numerical relativity” In _Physical Review D_ 67.10, 2003, pp. 104005 DOI: 10.1103/PhysRevD.67.104005
  * [10] W. Boscheri, M. Dumbser, M. Ioriatti, I. Peshkov and E. Romenski “A structure-preserving staggered semi-implicit finite volume scheme for continuum mechanics” In _Journal of Computational Physics_ 424, 2021, pp. 109866 DOI: 10.1016/j.jcp.2020.109866
  * [11] L.. Buchman and J.. Bardeen “Hyperbolic tetrad formulation of the Einstein equations for numerical relativity” In _Physical Review D_ 67.8, 2003, pp. 084017 DOI: 10.1103/PhysRevD.67.084017
  * [12] Saray Busto, Simone Chiocchetti, Michael Dumbser, Elena Gaburro and Ilya Peshkov “High Order ADER Schemes for Continuum Mechanics” In _Frontiers in Physics_ 8.32, 2020 DOI: 10.3389/fphy.2020.00032
  * [13] Saray Busto, Michael Dumbser, Ilya Peshkov and Evgeniy Romenski “On Thermodynamically Compatible Finite Volume Schemes for Continuum Mechanics” In _SIAM Journal on Scientific Computing_ 44.3, 2022, pp. A1723–A1751 DOI: 10.1137/21M1417508
  * [14] Yi-Fu Cai, Salvatore Capozziello, Mariafelicia De Laurentis and Emmanuel N. Saridakis “f ( T ) teleparallel gravity and cosmology” In _Reports on Progress in Physics_ 79.10 IOP Publishing, 2016, pp. 106901 DOI: 10.1088/0034-4885/79/10/106901
  * [15] A Dedner, F Kemm, D Kröner, C.-D. Munz, T Schnitzer and M Wesenberg “Hyperbolic Divergence Cleaning for the MHD Equations” In _Journal of Computational Physics_ 175.2, 2002, pp. 645–673 DOI: 10.1006/jcph.2001.6961
  * [16] Michael Dumbser, Francesco Fambri, Elena Gaburro and Anne Reinarz “On GLM curl cleaning for a first order reduction of the CCZ4 formulation of the Einstein field equations” In _Journal of Computational Physics_ , 2019, pp. 109088 DOI: 10.1016/j.jcp.2019.109088
  * [17] Michael Dumbser, Francesco Fambri, Maurizio Tavelli, Michael Bader and Tobias Weinzierl “Efficient Implementation of ADER Discontinuous Galerkin Schemes for a Scalable Hyperbolic PDE Engine” In _Axioms_ 7.3 Multidisciplinary Digital Publishing Institute, 2018, pp. 63 DOI: 10.3390/axioms7030063
  * [18] Michael Dumbser, Federico Guercilena, Sven Köppel, Luciano Rezzolla and Olindo Zanotti “Conformal and covariant Z4 formulation of the Einstein equations: Strongly hyperbolic first-order reduction and solution with discontinuous Galerkin schemes” In _Physical Review D_ 97.8, 2018, pp. 084053 DOI: 10.1103/PhysRevD.97.084053
  * [19] Michael Dumbser, Ilya Peshkov, Evgeniy Romenski and Olindo Zanotti “High order ADER schemes for a unified first order hyperbolic formulation of continuum mechanics: Viscous heat-conducting fluids and elastic solids” In _Journal of Computational Physics_ 314, 2016, pp. 824–862 DOI: 10.1016/j.jcp.2016.02.015
  * [20] Michael Dumbser, Ilya Peshkov, Evgeniy Romenski and Olindo Zanotti “High order ADER schemes for a unified first order hyperbolic formulation of Newtonian continuum mechanics coupled with electro-dynamics” In _Journal of Computational Physics_ 348 Elsevier Inc., 2017, pp. 298–342 DOI: 10.1016/j.jcp.2017.07.020
  * [21] Frank B. Estabrook, R. Robinson and Hugo D. Wahlquist “Hyperbolic equations for vacuum gravity using special orthonormal frames” In _Classical and Quantum Gravity_ 14.5, 1997, pp. 1237–1247 DOI: 10.1088/0264-9381/14/5/025
  * [22] Francesco Fambri “A novel structure preserving semi‐implicit finite volume method for viscous and resistive magnetohydrodynamics” In _International Journal for Numerical Methods in Fluids_ 93.12, 2021, pp. 3447–3489 DOI: 10.1002/fld.5041
  * [23] Elena Gaburro, Manuel J. Castro and Michael Dumbser “A Well Balanced Finite Volume Scheme for General Relativity” In _SIAM Journal on Scientific Computing_ 43.6, 2021, pp. B1226–B1251 DOI: 10.1137/21M1399154
  * [24] S K Godunov “An interesting class of quasilinear systems” In _Dokl. Akad. Nauk SSSR_ 139(3), 1961, pp. 521–523
  * [25] S.. Godunov, T.. Mikhaîlova and E.. Romenskiî “Systems of thermodynamically coordinated laws of conservation invariant under rotations” In _Siberian Mathematical Journal_ 37.4 Springer, 1996, pp. 690–705 DOI: 10.1007/BF02104662
  * [26] Alexey Golovnev, Tomi Koivisto and Marit Sandstad “On the covariance of teleparallel gravity theories” In _Classical and Quantum Gravity_ 34.14 IOP Publishing, 2017, pp. 145013 DOI: 10.1088/1361-6382/aa7830
  * [27] Eric Gourgoulhon “3+1 Formalism in General Relativity” In _Lecture Notes in Physics_ 846, Lecture Notes in Physics Berlin Heidelberg: Springer-Verlag, 2012 DOI: 10.1007/978-3-642-24525-1
  * [28] Pavel Grinfeld “Determinants and the Levi-Civita Symbol” In _Introduction to Tensor Analysis and the Calculus of Moving Surfaces_ New York, NY: Springer New York, 2013, pp. 133–157 DOI: 10.1007/978-1-4614-7867-6˙9
  * [29] Friedrich W. Hehl, Paul Heyde, G. Kerlick and James M. Nester “General relativity with spin and torsion: Foundations and prospects” In _Reviews of Modern Physics_ 48.3, 1976, pp. 393–416 DOI: 10.1103/RevModPhys.48.393
  * [30] Friedrich W. Hehl and Yuri N. Obukhov “Elie Cartan’s torsion in geometry and in field theory, an essay” In _Annales de la Fondation Louis de Broglie_ 32.2-3, 2007, pp. 157–194 arXiv: http://arxiv.org/abs/0711.1535
  * [31] Manuel Hohmann, Laur Järv, Martin Krššák and Christian Pfeifer “Teleparallel theories of gravity as analogue of nonlinear electrodynamics” In _Physical Review D_ 97.10, 2018, pp. 104042 DOI: 10.1103/PhysRevD.97.104042
  * [32] M.O Katanaev and I.V Volovich “Theory of defects in solids and three-dimensional gravity” In _Annals of Physics_ 216.1, 1992, pp. 1–28 DOI: 10.1016/0003-4916(52)90040-7
  * [33] Ondřej Kincl, Ilya Peshkov, Michal Pavelka and Václav Klika “Unified description of fluids and solids in Smoothed Particle Hydrodynamics”, 2022 arXiv: http://arxiv.org/abs/2207.00063
  * [34] Hagen Kleinert “Multivalued Fields” WORLD SCIENTIFIC, 2008 DOI: 10.1142/6742
  * [35] M. Krššák, R.. Hoogen, J.. Pereira, C.. Böhmer and A.. Coley “Teleparallel theories of gravity: illuminating a fully invariant approach” In _Classical and Quantum Gravity_ 36.18, 2019, pp. 183001 DOI: 10.1088/1361-6382/ab2e1f
  * [36] C.-D. Munz, P. Omnes, R. Schneider, E. Sonnendrücker and U. Voß “Divergence Correction Techniques for Maxwell Solvers Based on a Hyperbolic Model” In _Journal of Computational Physics_ 161.2, 2000, pp. 484–511 DOI: 10.1006/jcph.2000.6507
  * [37] Van Hoi Nguyen, Guy Casale and Loïc Le Marrec “On tangent geometry and generalised continuum with defects” In _Mathematics and Mechanics of Solids_ 27.7, 2022, pp. 1255–1283 DOI: 10.1177/10812865211059222
  * [38] Y.N. Obukhov “Electromagnetic energy and momentum in moving media” In _Annalen der Physik_ 17.9-10, 2008, pp. 830–851 DOI: 10.1002/andp.200810313
  * [39] H. Olivares, I.. Peshkov, E.. Most, F.. Guercilena and L.. Papenfort “A new first-order formulation of the Einstein equations exploiting analogies with electrodynamics” In _Physical Review D_ 105.124038, 2022 DOI: 10.1103/PhysRevD.105.124038
  * [40] Michal Pavelka, Václav Klika and Miroslav Grmela “Multiscale Thermo-Dynamics” Berlin, Boston: De Gruyter, 2018 DOI: 10.1515/9783110350951
  * [41] Ilya Peshkov, Michal Pavelka, Evgeniy Romenski and Miroslav Grmela “Continuum mechanics and thermodynamics in the Hamilton and the Godunov-type formulations” In _Continuum Mechanics and Thermodynamics_ 30.6, 2018, pp. 1343–1378 DOI: 10.1007/s00161-018-0621-2
  * [42] Ilya Peshkov and Evgeniy Romenski “A hyperbolic model for viscous Newtonian flows” In _Continuum Mechanics and Thermodynamics_ 28.1-2, 2016, pp. 85–104 DOI: 10.1007/s00161-014-0401-6
  * [43] Ilya Peshkov, Evgeniy Romenski and Michael Dumbser “Continuum mechanics with torsion” In _Continuum Mechanics and Thermodynamics_ 31.5 Springer Berlin Heidelberg, 2019, pp. 1517–1541 DOI: 10.1007/s00161-019-00770-6
  * [44] Kenneth G Powell, Philip L Roe, Timur J Linde, Tamas I Gombosi, Darren L De Zeeuw and W M Keck “A Solution-Adaptive Upwind Scheme for Ideal Magnetohydrodynamics” In _Journal of Computational Physics_ 154, 1999, pp. 284–309 URL: http://www.idealibrary.com
  * [45] Luciano Rezzolla and Olindo Zanotti “Relativistic Hydrodynamics” Oxford: Oxford University Press, 2013 DOI: 10.1093/acprof:oso/9780198528906.001.0001
  * [46] Evgeniy Romenski, Ilya Peshkov, Michael Dumbser and Francesco Fambri “A new continuum model for general relativistic viscous heat-conducting media” In _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 378.2170, 2020, pp. 20190175 DOI: 10.1098/rsta.2019.0175
  * [47] E I Romensky “Hyperbolic systems of thermodynamically compatible conservation laws in continuum mechanics” In _Mathematical and computer modelling_ 28.10, 1998, pp. 115–130 DOI: 10.1016/S0895-7177(98)00159-9
  * [48] Evgeniy I Romensky “Thermodynamics and hyperbolic systems of balance laws in continuum mechanics” In _Godunov Methods: Theory and Applications_ New York: Springer US, 2001, pp. 745–761 DOI: 10.1007/978-1-4615-0663-8
  * [49] Lewis Ryder “Introduction to General Relativity” Cambridge: Cambridge University Press, 2009 DOI: 10.1017/CBO9780511809033
  * [50] Arash Yavari and Alain Goriely “Riemann–Cartan Geometry of Nonlinear Dislocation Mechanics” In _Archive for Rational Mechanics and Analysis_ 205.1, 2012, pp. 59–118 DOI: 10.1007/s00205-012-0500-0