# Hodge-Dirac wave systems and structure-preserving discretizations of the linearized Einstein equations

Marien-Lorenzo Hanot, Kaibo Hu

► To cite this version:

Marien-Lorenzo Hanot, Kaibo Hu. Hodge-Dirac wave systems and structure-preserving discretizations of the linearized Einstein equations. 2025. hal-05363159

HAL Id: hal-05363159

https://hal.science/hal-05363159v1

Preprint submitted on 13 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Hodge-Dirac wave systems and structure-preserving discretizations of the linearized Einstein equations

Marien-Lorenzo Hanot,\* Kaibo Hu†

November 12, 2025

Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Reformulation of the ADM formulation</td><td>6</td></tr><tr><td>3</td><td>Continuous well-posedness</td><td>7</td></tr><tr><td>4</td><td>Variational formulation and discretization</td><td>9</td></tr><tr><td>4.1</td><td>Spatial discretization</td><td>10</td></tr><tr><td>4.2</td><td>Time discretization</td><td>10</td></tr><tr><td>4.3</td><td>Proof of the error estimate</td><td>12</td></tr><tr><td>4.4</td><td>Examples: a specific scheme</td><td>14</td></tr><tr><td>5</td><td>Numerical results</td><td>15</td></tr><tr><td>5.1</td><td>Discrete complex</td><td>15</td></tr><tr><td>5.2</td><td>Wave-like solutions</td><td>16</td></tr><tr><td>5.3</td><td>Results</td><td>16</td></tr><tr><td>6</td><td>Conclusions</td><td>17</td></tr><tr><td>A</td><td>Vector calculus identities</td><td>17</td></tr></table>

Abstract

We derive a reformulation of the linearized Arnowitt-Deser-Misner (ADM) equations as a Hodge-Dirac wave system with the divdiv complex, addressing challenges in numerical relativity such as gauge fixing, constraint propagation, and tensor symmetries. The differential and algebraic structures of the divdiv complex ensure the well-posedness of the formulation and facilitate structure-preserving discretization via finite element exterior calculus. We establish the well-posedness of this Hodge-Dirac wave equation and develop a discretization scheme applicable to both conforming and non-conforming discrete complexes, deriving error estimates under minimal assumptions.

**Key words.** Discrete divdiv complex, Hodge-Dirac wave, ADM equations, finite element exterior calculus

MSC2020. 65M12, 65N30, 65M60, 83C27

<footnote>\*Université de Lille, UMR 8524 - Laboratoire Paul Painlevé, CNRS, Inria, France. Email: marien-lorenzo.hanot@univ-lille.fr</footnote>

<footnote>†Mathematical Institute, University of Oxford. Radcliffe Observatory, Andrew Wiles Building, Oxford OX2 6GG, United Kingdom. Email: kaibo.hu@maths.ox.ac.uk</footnote>
# 1 Introduction

Numerical relativity is crucial in gravitational wave detection by providing essential templates. A central challenge in this field is solving the Einstein equations numerically. Significant advancements have been made over recent decades, notably the 2005 breakthrough in achieving accurate, long-term evolutions of black hole systems [32]. This success highlighted the critical role of the mathematical properties of the Einstein equations, particularly their hyperbolicity, in enabling robust numerical methods.

Despite these achievements, the increasing sophistication of gravitational wave detectors demands greater precision and long-term stability in numerical computations. Moreover, fundamental algorithmic challenges, such as a rigorous understanding of convergence properties, remain unresolved. The absence of comprehensive numerical analytic studies of the Einstein equations also hinders progress in exploring modified gravity models [26, 30].

The challenges of numerically solving the Einstein equations stem from at least the following issues:

1. Gauge fixing and hyperbolicity. As a geometric PDE, the Einstein equations permit reformulations through coordinate choices, known as gauge freedom. The mathematical properties, and thus the numerical performance, of these equations depend critically on the chosen gauge.

2. Constraint propagation and preservation. In 3+1 decompositions, geometric quantities must satisfy constraint equations at each time step, which are inherently preserved by the evolution equations. Maintaining these constraints numerically is both challenging and essential for ensuring numerical stability.

3. Nonlinearity. As characteristic of geometric PDEs, the Einstein equations exhibit strong nonlinearity.

Gauge fixing and constraint preservation necessitate a deep understanding of the differential and algebraic structures of both continuous and discrete equations. Over recent decades, Finite Element Exterior Calculus (FEEC) [2, 5, 6] has made significant strides in structure-preserving discretization and efficient solvers, with applications in electromagnetism and continuum mechanics. Differential and cohomological structures have proven critical for accurate and efficient numerical solutions. In the context of the Einstein equations, FEEC-based approaches were pioneered in [33] using the Einstein-Bianchi (first-order) formulation and in [24] with Regge finite elements, inspired by Regge calculus. The Einstein-Bianchi formulation in [33] employs traceless-transverse (TT, i.e., symmetric, traceless, and divergence-free) matrix fields, which naturally arise in conformal complexes [3, 10]. However, enforcing both symmetry and tracelessness constraints simultaneously poses significant challenges, although some conforming finite element conformal complexes with supersmoothness were recently constructed [16, 20, 22]. To encode the tensor symmetries, [33] utilized the Hessian (and div/div) complex with weaker algebraic constraints, imposing them via Lagrange multipliers, following approaches in elasticity [4]. Similarly, [17] explored the Einstein-Bianchi system with weakly imposed symmetries. The work in [24] extended Regge finite elements and analyzed their properties. We also mention the recent numerical scheme based on a differential form formulation [28]. In this paper, we build on the paradigm of structure-preserving discretizations for the Einstein equations. Specifically, we propose a Hilbert complex-based reformulation and its discretization, which naturally incorporates all algebraic and differential constraints *strongly* and ensures well-posedness.

We begin with the Arnowitt-Deser-Misner (ADM) equations, which govern the evolution of a metric and its curvature, satisfying the Einstein field equations. These equations arise from a (3 + 1)-decomposition of the field equations, separating temporal and spatial derivatives. A detailed exposition of the (3 + 1) formalism and the ADM equations can be found, for example, in [1]. In a (3 + 1)-decomposition, one assumes a foliation of spacetime with spacial slices. The distance between the
spacial slices is described by a function $\alpha$, referred to the *lapse*; and the *shift* between slices is described by a vector quantity $\beta$.

In this work, we focus on the York version of the ADM equations [35], particularly their linearization, given by

$$ \gamma_{tt} + S \mathbf{inc} \gamma - 2 \mathbf{hess} \alpha - 2 \mathbf{def} \beta_t = 0, \quad (1a) $$

$$ \mathbf{div} S(\gamma_t - 2 \mathbf{def} \beta) = 0, \quad (1b) $$

$$ \mathbf{div} \mathbf{div} S\gamma = 0, \quad (1c) $$

where $S: M_{3 \times 3}(\mathbb{R}) \to M_{3 \times 3}(\mathbb{R})$ is an operator defined for any $A \in M_{3 \times 3}(\mathbb{R})$ as $S(A) := A^\top - \text{tr}(A)I_3$, $\gamma: \mathbb{R}^4 \to \mathbb{S}$ represents the perturbation of the spatial metric, $1 + \alpha: \mathbb{R}^4 \to \mathbb{R}$ is the lapse function, and $\beta: \mathbb{R}^4 \to \mathbb{R}^3$ is the shift vector. The linearization is performed around $I + \gamma \approx I$ (Euclidean spacial metric), $\beta \approx 0$ (zero shift between slices), and $1 + \alpha \approx 1$ (uniform distance between slices). In adapted coordinates, the full spacetime metric $g$ is recovered as:

$$ g := \begin{pmatrix} -(1+\alpha)^2 + \beta \cdot \beta & \beta^\top \\ \beta & I+\gamma \end{pmatrix}. $$

The lapse function $\alpha$ and shift vector $\beta$ represent coordinate choices rather than physical quantities, serving as gauge functions. For the operators, **inc** is a row-wise curl composed with a column-wise curl; **def** := $\text{sym } \text{grad}$ is the symmetric gradient; **hess** is the Hessian operator of a scalar function; $S\sigma := \sigma^T - \frac{1}{2} \text{tr}(\sigma)I$ is a bijective algebraic operator (see [3]).

The ADM formulation (1) comprises two sets of equations: (1a) represents the evolutionary equation, while (1b)–(1c) are constraint equations that must be satisfied throughout the evolution. The system exhibits constraint propagation, meaning that if (1b)–(1c) hold for the initial data, they remain satisfied during the evolution. The standard ADM formulation (1) and its nonlinear counterpart lack hyperbolicity. Moreover, straightforward discretization of the evolutionary equations in (1) introduces numerical errors in the constraints, which accumulate over time and lead to instability.

**Hodge-Dirac reformulation.** The first contribution of this paper is to address the loss of constraint propagation, well-posedness and tensor symmetries in numerical discretization by reformulating (1) as a mixed formulation based on a differential complex. Then in the framework of finite element exterior calculus [2, 5, 6], we can discretize such formulations by discretizing the corresponding complexes and ensure well-posedness and structure-preserving properties.

To encode the symmetries of tensors in (1), we will reformulate the ADM formulation (1) into a canonical form based on the $\mathbf{div} \mathbf{div}$ complex [3, 10]:

$$ 0 \longrightarrow H(\text{dev grad}; \mathbb{V}) \xrightarrow{\text{dev grad}} H(\text{sym curl}; \mathbb{T}) \xrightarrow{\text{sym curl}} H(\text{div div}; \mathbb{S}) \xrightarrow{\text{div div}} L^2(\mathbb{R}) \longrightarrow 0. \quad (2) $$

More specifically, denoting by $d$ the differential operators in (2), i.e.,

$$ d := \begin{pmatrix} 0 & 0 & 0 & 0 \\ \text{dev grad} & 0 & 0 & 0 \\ 0 & \text{sym curl} & 0 & 0 \\ 0 & 0 & \text{div div} & 0 \end{pmatrix} $$

with its (formal) adjoint

$$ d^* := \begin{pmatrix} 0 & -\mathbf{div} & 0 & 0 \\ 0 & 0 & \mathbf{curl} & 0 \\ 0 & 0 & 0 & \mathbf{hess} \\ 0 & 0 & 0 & 0 \end{pmatrix} $$
and defining $J := \text{diag}(1, -1, 1, -1)$, we will show that (1) can be reformulated as $\partial_t U = J(d + d^*)U$ with properly chosen variables $U$ with components from all the spaces in (2) (see (8)). This formulation closely resembles a time-dependent, skew-symmetric Hodge-Dirac problem for the de Rham complex. See [23] for a study for stationary problems. Analogous to how the Hodge-Dirac problem can be viewed as a “square root” of a Hodge-Laplace problem, the proposed problem can be regarded as a “square root” of the Hodge(-Laplace) wave equation presented in [2, Section 8.5] and [34] with the divdiv complex (2). By analogy, we refer to the problem (8) as the *Hodge-Dirac wave*. The primary distinction between the two wave problems is that the Hodge-Laplace equation involves only three consecutive spaces in the complex, whereas the Hodge-Dirac wave considers the entire complex. For example, in a complex comprising four spaces, the mixed formulations of the Hodge-Laplace (HL) wave equation and the Hodge-Dirac (HD) wave equation take the forms:

$$ (\text{HL}): \frac{d}{dt} U = \begin{pmatrix} 0 & \delta & 0 & 0 \\ -d & 0 & -\delta & 0 \\ 0 & d & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} U, \quad (\text{HD}): \frac{d}{dt} U = \begin{pmatrix} 0 & \delta & 0 & 0 \\ -d & 0 & -\delta & 0 \\ 0 & d & 0 & \delta \\ 0 & 0 & -d & 0 \end{pmatrix} U, $$

respectively. For time-independent equations, stability requires considering the entire complex [23]. The Hodge-Dirac wave equation, and thus the reformulation of the ADM formulation, is well-posed. Each of the four components in the Hodge-Dirac system satisfies a Hodge-Laplacian wave equation; while the fourth component in the Hodge-Laplacian wave equation is constant in time, which may exhibit numerical instability under certain perturbations in the nonlinear case. Therefore, one may expect that the Hodge-Dirac system has better stability than the Hodge-Laplacian wave.

Compared to previous works that used differential complexes to address systems in numerical relativity, such as [16, 17, 34], to the best of our knowledge, this is the first work that encode all the algebraic symmetries of the tensors and the constraint equations in the ADM-Einstein equations with the divdiv complex. These constraints are thus preserved in numerical discretization.

**Discretization of Hodge-Dirac systems.** As another main contribution of this paper, we propose a scheme for solving the Hodge-Dirac wave, establish its well-posedness, and derive error estimates for the time-dependent problem under minimal assumptions. Although presented in three dimensions, the approach is straightforward to extend to any number of dimensions. More importantly, the construction and results do not require conforming discretizations. Any discrete complex may be used, provided it possesses (not necessarily bounded) commuting interpolators and (discrete) Poincaré inequalities. The error estimates then depend on the consistency properties of the chosen spaces. The ability to use non-conforming spaces is particularly valuable for complexes more intricate than the de Rham complex, where conforming discretizations are often challenging to construct and come with other limitations, such as high minimal polynomial degrees or the need for special meshes.

The ADM-type formulation discussed in this paper involves the divdiv complex (2). The divdiv complex incorporates a second-order differential operator and tensor-valued elements with additional symmetry constraints (e.g., traceless or symmetric elements), posing significant challenges for discretization. We briefly review existing approaches.

* Some approaches relax strong symmetry constraints and enforce them weakly, as done for the linear elasticity [4] and the linearized Einstein-Bianchi system [34]. We choose, however, to preserve strong symmetries within the spaces.

* Conforming discretizations of the spaces in (2) were developed in [8, 12, 17, 18]. Under certain assumptions, conforming spaces and complexes are viable for our scheme. However, most constructions on simplicial meshes are involved due to supersmoothness constraints, which also exclude the use of low-order polynomials. In the numerical results presented in this paper, we
use a div**div** discrete complex based on tensor-product splines [8]. This conforming complex is restricted to Cartesian meshes but is significantly simpler to implement than other alternatives considered. It also yields a highly structured matrix system for numerical schemes. Although the overall degree of local polynomials in the complex can be high (reaching degree 7 in some spaces), the degree in any Cartesian direction remains below 3. We discuss this complex in greater detail in Section 5.

* To reduce the overall complexity of the method, non-conforming discretizations are also considered. A “fully discrete” div**div** complex was developed in [31], using collections of local polynomials associated with various mesh entities to represent discrete data. This approach supports arbitrary polyhedral elements, not just simplicial meshes, and allows the use of lower-degree polynomials compared to conforming methods. Nevertheless, the spaces remain complex.

* Another approach involves distributional elements. This concept, explored in [9, 13, 15, 27, 29], incorporates Dirac deltas into the design of numerical schemes. Distributional spaces can be viewed either through duality with a conforming space or as functions attached to lower-dimensional entities. Such schemes rely on intrinsic finite elements with weaker regularity and their complexes, which have gained attention since Braess and Schöberl's work on a posteriori estimators [9] and Christiansen's finite element reinterpretation of Regge calculus [13], as well as in recent works [7, 14, 19, 21, 25]. In particular, a distributional div**div** complex was introduced in [21]. While distributional finite elements offer a promising approach, verifying the conditions arising from our analysis remains an open task for future work.

The rest of the paper is organized as follows. In Section 2, we reformulate the linearized ADM equations on the div**div** complex. In Section 3, we abstract the resulting formulation to a more general setting of “Hodge-Dirac wave” and show the well-posedness of the continuous problem. In Section 4, we study the discretization of this problem on discrete, not necessary conforming, complexes. Obtaining both the well-posedness, and error estimates under fairly general assumptions. Lastly, in Section 5, we introduce some classes of exact solutions, and provide and discuss numerical results obtained with the scheme.

**Notation.** In this paper, we define operators acting on matrix fields column-wise. This follows the convention in [3]. We introduce some notation following [3]. Moreover, mskw : $\mathbb{R}^3 \mapsto \mathbb{K}$ defined by

<table><caption>Table 1: Notations</caption><tr><td>V</td><td>R<sup>n</sup></td></tr><tr><td>M</td><td>the space of all n × n-matrices</td></tr><tr><td>S</td><td>symmetric matrices</td></tr><tr><td>K</td><td>skew symmetric matrices</td></tr><tr><td>T</td><td>trace-free matrices</td></tr><tr><td>skw : M → K</td><td>skew symmetric part of a matrix</td></tr><tr><td>sym : M → S</td><td>symmetric part of a matrix</td></tr><tr><td>tr : M → R</td><td>matrix trace</td></tr><tr><td>ι : R → M</td><td>the map u := uI identifying a scalar with a scalar matrix</td></tr><tr><td>dev : M → T</td><td>deviator (trace-free part of a matrix) given by dev u := u − 1/n tr(u)I</td></tr><tr><td>S : M → M</td><td>the map given by Su = uT − tr(u)I</td></tr></table>

$(mskw V)_{ij} := -\epsilon_{ijk}V_k$ maps an axial vector to its matrix representation; and $(vskw M)_i := -\frac{1}{2}\epsilon_{ijk}M_{jk}$ takes the skew-symmetric part of a matrix and maps it to the axial vector.
## 2 Reformulation of the ADM formulation

As mentioned in the introduction, we want to find a reformulation of the linearized ADM equations (1) based on the divdiv complex (2). Incorporating the constraint equations of (1), we propose the following strong formulation:

$$ \mathrm{div} A = 0 \qquad (3a) $$

$$ A_t + \mathrm{curl} S\gamma = 0 \qquad (3b) $$

$$ S\gamma_t - \mathrm{sym} \mathrm{curl} A = 0 \qquad (3c) $$

$$ \mathrm{divdiv} S\gamma = 0. \qquad (3d) $$

The choice of variables is summarized in the following complex:

<figure><img src="image_2.png" /></figure>

Let us assume for now that (3) is well-posed and verify that (3) is equivalent to the ADM formulation.

**Theorem 1.** The solution $\gamma$ of (3) with initial conditions ($A_0, \gamma_0$) satisfies the linearized York version of the ADM formulation (1) for the gauge $\alpha := \frac{1}{2} \mathrm{tr} \gamma$ and $\beta := \frac{1}{2} \int_0^t \mathrm{div} S\gamma + \mathrm{vskw} A_0$.

*Proof.* The proof is divided into Lemma 4 and Lemma 5. $\square$

*Remark 2.* The main difficulty is to derive the correct gauge for this formulation. The skew symmetric part of $A$ introduces a non-zero value of $\mathrm{div} S\gamma_t$ which must be corrected using the shift vector $\beta$ (see the proof of Lemma 5). While the lapse function $\alpha$ is used to absorb the contribution of $\mathrm{tr} \gamma$ (see the proof of Lemma 4).

*Remark 3.* The main advantage of this system over an inc based formulation is that, without the time derivatives, it reduces to a simple Hodge-Dirac problem. The latter is well studied in the context of FEEC (see [23]), and can be discretized in a stable way using any divdiv-complex. Moreover, it strongly enforces both constraints. The $\mathrm{divdiv} S\gamma = 0$ constraint is directly enforced, while the $\mathrm{div} S\gamma_t = 2 \mathrm{div} S \mathrm{def} \beta$ can be taken as a definition for the shift vector $\beta$ (see Lemma 4).

**Lemma 4.** If $(A, \gamma)$ is a solution of (3) then

$$ \gamma_{tt} + S \mathrm{inc} \gamma - 2 \mathrm{def} \beta_t - 2 \mathrm{hess} \alpha = 0, \qquad (4) $$

where $\alpha = \frac{1}{2} \mathrm{tr} \gamma$ and $\beta = \frac{1}{2} \int_0^t \mathrm{div} S\gamma + \mathrm{vskw} A_0$.

*Proof.* Taking the time derivative of (3c), we have

$$ \begin{align*} \gamma_{tt} - S^{-1} \mathrm{sym} \mathrm{curl} A_t &= 0 \\ \gamma_{tt} + S^{-1} \mathrm{sym} \mathrm{curl curl} S\gamma &= 0 \tag{3b} \\ \gamma_{tt} + S^{-1} (\mathrm{inc} \gamma - S \mathrm{def} \mathrm{div} S\gamma - S \mathrm{hess} \mathrm{tr} \gamma) &= 0 \tag{29} \\ \gamma_{tt} + S^{-1} \mathrm{inc} \gamma - 2 \mathrm{def} \beta_t - 2 \mathrm{hess} \alpha &= 0 \\ \gamma_{tt} + S \mathrm{inc} \gamma + \frac{1}{2} \mathrm{tr} \mathrm{inc} \gamma - 2 \mathrm{def} \beta_t - 2 \mathrm{hess} \alpha &= 0 \\ \gamma_{tt} + S \mathrm{inc} \gamma - \frac{1}{2} \mathrm{div} \mathrm{div} S\gamma - 2 \mathrm{def} \beta_t - 2 \mathrm{hess} \alpha &= 0 \tag{3d} \end{align*} $$

where we used the fact that $S^{-1} = S + \frac{1}{2} \mathrm{tr}$ on the fifth line, and that $\mathrm{tr} \mathrm{inc} = - \mathrm{divdiv} S$ on the last. $\square$
**Lemma 5.** If $(A, \gamma)$ is solution of (3) then

$$ \mathrm{div} S\gamma_t = 2 \mathrm{div} S \mathrm{def} \beta, \qquad (5) $$

$$ \text{where } \beta = \frac{1}{2} \int_{0}^{t} \mathrm{div} S\gamma + \mathrm{vskw} A_{0}. $$

*Proof.* Taking the divergence of (3c), we have

$$ \mathrm{div} S\gamma_t = \mathrm{div} \, \mathrm{sym} \, \mathrm{curl} A $$

$$ \stackrel{(3a),(32)}{=} - \mathrm{curl} \, \mathrm{curl} \, \mathrm{vskw} A $$

$$ = - \mathrm{curl} \, \mathrm{curl} \, \mathrm{vskw} \left( A_0 + \int_0^t A_t \right) $$

$$ \stackrel{(3b)}{=} - \mathrm{curl} \, \mathrm{curl} \left( \mathrm{vskw} A_0 - \int_0^t \mathrm{vskw} \, \mathrm{curl} S\gamma \right) $$

$$ = - \mathrm{curl} \, \mathrm{curl} \left( \mathrm{vskw} A_0 + \frac{1}{2} \int_0^t \mathrm{div} S^2 \gamma \right) $$

$$ = - \mathrm{curl} \, \mathrm{curl} \left( \mathrm{vskw} A_0 + \frac{1}{2} \int_0^t \mathrm{div} S\gamma - \frac{1}{2} \int_0^t \mathrm{grad} \, \mathrm{tr} S\gamma \right) $$

$$ = - \mathrm{curl} \, \mathrm{curl} \beta + \frac{1}{2} \int_0^t \mathrm{curl} \, \mathrm{curl} \, \mathrm{grad} \, \mathrm{tr} S\gamma $$

$$ \stackrel{(34)}{=} 2 \mathrm{div} S \mathrm{def} \beta, $$

where $A_0$ is the value of $A$ at $t = 0$. We used the identity $2 \mathrm{vskw} \, \mathrm{curl} = - \mathrm{div} S$ on the fifth line, and the fact the $S\gamma \in \mathbb{S}$ on the sixth. $\square$

*Remark 6* (Recovering the shift vector). The shift vector $\beta$ may be computed as the solution of

$$ \mathrm{curl} \, \mathrm{curl} \, \beta = -\frac{1}{2} \mathrm{div} S\gamma_t $$

$$ \mathrm{div} \, \beta = \mathrm{div} \, \mathrm{vskw} A_0. $$

### 3 Continuous well-posedness

To simplify the notations and emphasize the role played by the divdiv complex (2), we will denote the differential operators by either $d$ or $d^*$. In the following, the domain of $d^*$ is determined from the domain of $d$ to make $(d^*, D(d^*))$ the adjoint of $(d, D(d))$. Using this notation, and writing $g := S\gamma$, the system (3) becomes:

$$ d^* A = 0 $$

$$ A_t + d^* g = 0 $$

$$ g_t - dA = 0 $$

$$ dg = 0. $$

(6)

In order to show the well-posedness of (6), we will first complete it into the following system:

$$ \partial_t \lambda_0 = d^* A $$

$$ \partial_t A = -d^* g - d\lambda_0 $$

$$ \partial_t g = dA + d^* \lambda_3 $$

$$ \partial_t \lambda_3 = -dg. $$

(7)
We can readily verify that we retrieve (6) from (7) if $\lambda_0 = 0$ and $\lambda_3 = 0$.

Let

$$J := \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & -1 \end{pmatrix}, \quad \mathcal{L} := J(d + d^*) := \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & -1 \end{pmatrix} \begin{pmatrix} 0 & d^* & 0 & 0 \\ d & 0 & d^* & 0 \\ 0 & d & 0 & d^* \\ 0 & 0 & d & 0 \end{pmatrix}, \quad U := \begin{pmatrix} \lambda_0 \\ A \\ g \\ \lambda_3 \end{pmatrix}.$$

The system (7) can be written in a more condensed form as

$$\partial_t U = \mathcal{L} U. \qquad (8)$$

We define the domain of $\mathcal{L}$ as

$$D(\mathcal{L}) := D(d) \cap D(d^*).$$

**Lemma 7 (Skew adjointness).** The operator $\mathcal{L}$ is skew-adjoint as an unbounded operator.

*Proof.* Noticing that, for all $0 \le k < 3$, $v \in D(d)$ such that only the $k$-th component of $v$ is non-zero, $dJv = (-1)^{k+1}dv$, while $Jdv = (-1)^k dv$, we have $dJ = -Jd$. Therefore, since $J$ is self-adjoint, we have formally

$$(Jd)^* = d^*J^* = d^*J = -Jd^*.$$

Using the same identity for $Jd^*$, we find:

$$\mathcal{L}^* = (Jd + Jd^*)^* = -(Jd^* + Jd) = -\mathcal{L}.$$

Viewing $\mathcal{L}$ as an unbounded operator defined on $D(d) \cap D(d^*)$, we find the domain of $\mathcal{L}^*$ to be $D(d^*) \cap D(d) = D(\mathcal{L})$. Therefore $(\mathcal{L}, D(\mathcal{L}))$ is skew-adjoint as an unbounded operator. $\square$

The main tool to study the system (8) is the following result from [11, Proposition 4.1.6, Corollary 2.4.9]:

**Theorem 8 (Hille-Yosida).** Given a Hilbert space $X$, a skew-adjoint operator $A$ with domain $D(A)$ dense in $X$, $U_0 \in D(A)$ and $f \in C([0, T], X)$ such that $f \in W^{1,1}([0, T], X)$ or $f \in L^1((0, T), D(A))$. There exists a unique solution

$$U \in C^0([0, T], D(A)) \cap C^1([0, T], X)$$

to the system $\partial_t U = AU + f$ with $U(0) = U_0$. Moreover, if $f \in L^1((0, T), X)$ then by [11, Lemma 4.1.5],

$$\|U\|_{C([0,T],X)} \le \|U_0\|_X + \|f\|_{L^1((0,T),X)}.$$

*Remark 9.* From [11, Proposition 6.1.1], we also have the energy estimate for skew-symmetric operators

$$\frac{1}{2} \frac{d}{dt} \|U\|^2 = \langle f, U \rangle.$$

$$\text{Let } X = L^2(\Omega, \mathbb{V}) \oplus L^2(\Omega, \mathbb{T}) \oplus L^2(\Omega, \mathbb{S}) \oplus L^2(\Omega, \mathbb{R}).$$

**Theorem 10 (Well-posedness).** If $D(\mathcal{L})$ is dense in $X$, and the initial condition $(\lambda_0^0, A_0, S\gamma_0, \lambda_3^0) \in D(\mathcal{L})$, then there exists a unique solution $(\lambda_0, A, S\gamma, \lambda_3)$ to the system (7). Moreover, its norm is bounded by

$$\|(\lambda_0, A, S\gamma, \lambda_3)\|_{C([0,T],X)} \le \|(\lambda_0^0, A_0, S\gamma_0, \lambda_3^0)\|_X.$$

If the initial conditions are compatible with (3) (i.e. $\lambda_0^0 = 0$, $\lambda_3^0 = 0$, $\mathrm{div} A_0 = 0$, and $\mathrm{div} \mathrm{div} S\gamma_0 = 0$), then $(A, S\gamma)$ is the unique solution of the system (3).
Proof. If the initial condition $(\lambda_0^0, A_0, S\gamma_0, \lambda_3^0) \in D(\mathcal{L})$, then Lemma 7 allows us to apply Theorem 8 with $f = 0$ to show the existence and uniqueness of a solution $U := (\lambda_0, A, S\gamma, \lambda_3)$ of (8), together with the bound on the norm.

Let us now assume that $\lambda_0^0 = 0$, $\lambda_3^0 = 0$, $\mathrm{div} A_0 = 0$, and $\mathrm{div} \, \mathrm{div} S\gamma_0 = 0$. It remains to prove that $\lambda_0 = 0$ and $\lambda_3 = 0$ to retrieve (6). To this end, we will first show that $W := (0, A - P_{\mathrm{Im} \, \mathrm{dev} \, \mathrm{grad}} A, S\gamma - P_{\mathrm{Im} \, \mathrm{hess}} S\gamma, 0)$ is another solution and concludes with the uniqueness of the solution. Using the Hodge decomposition

$$L^2 \otimes \mathbb{T} = \mathrm{Im} \, \mathrm{dev} \, \mathrm{grad} \oplus \mathrm{Im} \, \mathrm{curl} \oplus \mathfrak{S}^1, \quad L^2 \otimes \mathbb{S} = \mathrm{Im} \, \mathrm{hess} \oplus \mathrm{Im} \, \mathrm{sym} \, \mathrm{curl} \oplus \mathfrak{S}^2,$$

where $\mathfrak{S}^1$ and $\mathfrak{S}^2$ are the (possibly empty) spaces of harmonic forms, we deduce that $\mathrm{div}(A - P_{\mathrm{Im} \, \mathrm{dev} \, \mathrm{grad}} A) = \mathrm{div}(P_{\mathrm{Im} \, \mathrm{curl}} \oplus \mathfrak{S}^1 A) = 0$, and $\mathrm{div} \, \mathrm{div}(S\gamma - P_{\mathrm{Im} \, \mathrm{hess}} S\gamma) = \mathrm{div} \, \mathrm{div} \, P_{\mathrm{Im} \, \mathrm{sym} \, \mathrm{curl}} \oplus \mathfrak{S}^2 S\gamma = 0$. Moreover, we have $\mathrm{sym} \, \mathrm{curl}(A - P_{\mathrm{Im} \, \mathrm{dev} \, \mathrm{grad}} A) = \mathrm{sym} \, \mathrm{curl} A$, and $\mathrm{curl}(S\gamma - P_{\mathrm{Im} \, \mathrm{hess}} S\gamma) = \mathrm{curl} S\gamma$. Therefore, at all time $t$, $W \in D(\mathcal{L})$. Since the orthogonal projections are continuous and commute with the time derivative, we infer that $W \in C^0([0, T], D(\mathcal{L})) \cap C^1([0, T], X)$. Since $\mathrm{div} A_0 = 0$, we have $(A - P_{\mathrm{Im} \, \mathrm{dev} \, \mathrm{grad}} A)_0 = A_0$, and, likewise $\mathrm{div} \, \mathrm{div} S\gamma_0 = 0$ implies $(S\gamma - P_{\mathrm{Im} \, \mathrm{hess}} S\gamma)_0 = S\gamma_0$. It only remains to check that $\partial_t W = \mathcal{L}W$, that we infer from the Hodge decomposition giving $\partial_t A = \partial_t P_{\mathrm{Im} \, \mathrm{dev} \, \mathrm{grad}} A - \mathrm{curl} S\gamma$, and $\partial_t S\gamma = \partial_t P_{\mathrm{Im} \, \mathrm{hess}} S\gamma + \mathrm{sym} \, \mathrm{curl} A$.

Therefore $(0, A - P_{\mathrm{Im} \, \mathrm{dev} \, \mathrm{grad}} A, S\gamma - P_{\mathrm{Im} \, \mathrm{hess}} S\gamma, 0)$ is a solution of (8), and by uniqueness we must have $U = W$, hence $\lambda_0 = 0$ and $\lambda_3 = 0$. Thus $(A, \gamma)$ is a solution of (3). We readily verify that any solution of (3) gives a solution of (8), showing the uniqueness of $(A, \gamma)$. $\square$

## 4 Variational formulation and discretization

Consider the problem: Given $U_0 \in D(d) \cap D(d^*)$, find $U \in C^1([0, T], X) \cap C^0([0, T], D(d))$, such that $U(0) = U_0$, and for all $t \in [0, T]$, all $V \in D(d)$,

$$(\partial_t U, V) = (JdU, V) - (U, JdV). \qquad (9)$$

**Lemma 11.** A function $U$ is a solution of (9) if and only if $U$ is a solution of (8).

The well-posedness of (9) follows from Theorem 10.

Let us now consider the discretization of this problem. For the discretization in space, we consider a discrete complex:

$$L^2(\Omega, \mathbb{V}) \xrightarrow{d} L^2(\Omega, \mathbb{T}) \xrightarrow{d} L^2(\Omega, \mathbb{S}) \xrightarrow{d} L^2(\Omega, \mathbb{R}) \\ \xrightarrow{\substack{\downarrow \\ \underline{\mathbf{X}}_{r,h}^0}} \xrightarrow{d_h} \underline{\mathbf{X}}_{r,h}^1 \xrightarrow{d_h} \underline{\mathbf{X}}_{r,h}^2 \xrightarrow{d_h} \underline{\mathbf{X}}_{r,h}^3$$

We denote by $\underline{\mathbf{X}}_{r,h} := \times_{i=0}^3 \underline{\mathbf{X}}_{r,h}^i$ the global discrete space, and introduce two norms on the discrete space: the $L^2$-like norm $\|V\|_h := \sqrt{\sum_{i=0}^3 (v_i, v_i)_h}$, and the graph norm $\|V\|_{1,h} := \|V\|_h + \|d_h V\|_h$.

**Assumption 12.** The discrete complex must satisfy the following properties:

(A1) The discrete complex admits uniform Poincaré inequalities: there is $c_p$ independent of the mesh size $h$, such that for all $0 \le i < 3$,

$$\forall \underline{\tau}_h \in \underline{\mathbf{X}}_{r,h}^i \cap (\mathrm{Ker} \, d_h)^{\perp h}, \quad \|\underline{\tau}_h\|_h \le c_p \|d_h \underline{\tau}_h\|_h,$$

where $(\mathrm{Ker} \, d_h)^{\perp h}$ denotes the orthogonal complement of $\mathrm{Ker} \, d_h$ with respect to the $\|\cdot\|_h$ inner product.
(A2) The interpolator is a cochain map, i.e. $\forall 0 \le i < 3$, $d_h \underline{I}_{r,h}^i = \underline{I}_{r,h}^{i+1} d$.

Notice that the right-hand side of (9) is, in general, degenerate. If $W$ is such that $dW = 0$ and $d^*W = 0$ (if $W$ is an harmonic form), then (9) becomes $\partial_t W = 0$. Since the equation is linear, this simply means that the harmonic component of a solution is preserved by the evolution, and it does not require a particular treatment in the implementation. We introduce the notation $P_S$ for the $L^2$ orthonormal projector into the subspace of harmonic forms. We use the same notation for the projector acting on the continuous and the discrete spaces; its meaning can be deduced from its argument. With this convention of notation, we have on both continuous and discrete levels

$$ \partial_t P_S U = 0. \qquad (10) $$

We denote by $\mathfrak{S}_h \subset \underline{X}_{r,h}$ the space of discrete harmonic forms.

## 4.1 Spatial discretization

We first consider the discretization of the spatial derivative appearing in (9). We define the bilinear operator $\mathcal{L}_h : \underline{X}_{r,h} \times \underline{X}_{r,h} \to \mathbb{R}$ for all $V_h, W_h \in \underline{X}_{r,h}$ by

$$ \mathcal{L}_h (V_h, W_h) := (J d_h V_h, W_h)_h - (V_h, J d_h W_h)_h. \qquad (11) $$

Notice that $\mathcal{L}_h$ is continuous for the graph norm.

**Lemma 13 (Partial Inf-Sup stability of $\mathcal{L}_h$).** The bilinear form $\mathcal{L}_h$ is inf-sup stable on $\mathfrak{S}^{\perp h}$: there is $C_L > 0$ depending only on $c_p$ such that, for all $V_h \in \mathfrak{S}^{\perp h}$,

$$ \sup_{W_h \in \mathfrak{S}^{\perp h}, W_h \neq 0} \frac{\mathcal{L}_h(V_h, W_h)}{\|W_h\|_{1,h}} \geq C_L \|V_h\|_{1,h}. $$

*Proof.* Let $V_h \in \mathfrak{S}^{\perp h}$. We infer from the discrete Poincaré inequality Assumption (A1) the existence of $\rho_h \in \underline{X}_{r,h} \cap (\text{Ker } d_h)^{\perp h}$ such that $d_h \rho_h = -JP_{\text{Im } d_h} V_h$ and $\|\rho_h\|_h \le c_p \|V_h\|_h$. Since $V_h \in \mathfrak{S}^{\perp h}$, we have $P_{\text{Im } d_h} V_h = P_{\text{Ker } d_h} V_h$. Setting $W_h = J d_h V_h + \rho_h$, we have

$$ \begin{aligned} \mathcal{L}_h (V_h, W_h) &:= (J d_h V_h, J d_h V_h)_h + (J d_h V_h, \rho_h)_h - (V_h, J d_h \rho_h)_h \\ &= \|d_h V_h\|_h^2 + 0 + \|P_{\text{Ker } d_h} V_h\|_h^2 \\ &\ge \frac{1}{2} \|d_h V_h\|_h^2 + \frac{1}{2c_p^2} \|(\text{Id} - P_{\text{Ker } d_h}) V_h\|_h^2 + \|P_{\text{Ker } d_h} V_h\|_h^2 \quad (12) \\ &\ge \frac{1}{2 \max\{1, c_p^2\}} \|V_h\|_{1,h}^2. \end{aligned} $$

Moreover, from the construction of $W_h$, we have

$$ \|W_h\|_{1,h} = \|d_h V_h\|_h + \|\rho_h\|_h + \|d_h \rho_h\|_h \le (1+c_p)\|U\|_h + \|d_h V_h\|_h \le (1+c_p)\|V_h\|_{1,h}. \quad (13) $$

We infer the result dividing (12) by (13). $\square$

## 4.2 Time discretization

For the discretization in time, we consider a splitting of the time interval $t_0 = 0 < t_1 < \dots < t_N = T$ and, for each space of the complex (2), a linear operator approximating the time derivative at $t_n$, $\mathcal{T}^n := \mathcal{T}_B^n + \mathcal{T}_L^n$, together with its discrete counterpart $\mathcal{T}_h^n := \mathcal{T}_{h,B}^n + \mathcal{T}_{h,L}^n$. The splitting is done to separate the implicit part $\mathcal{T}_B^n$ from the explicit part $\mathcal{T}_L^n$ of the operator. Explicitly, the component $\mathcal{T}_{h,L}^n (V_h^m)_{0 \le m \le n} = \mathcal{T}_{h,L}^n (V_h^m)_{0 \le m < n}$ is a vector of $\underline{X}_{r,h}^i$ that does not depend on $V_h^n$, while the component $\mathcal{T}_{h,B}^n (V_h^m)_{0 \le m \le n} = \mathcal{T}_{h,B}^n V_h^n$ is linear in $V_h^n$.
**Example 14.** For the Euler method, the time discretization operator reads:

$$ \mathcal{T}_{h,B}^{n}((V_{h}^{m})_{m \le n}) := \frac{V^{n}}{t_{n} - t_{n-1}}, \quad \mathcal{T}_{h,L}^{n}((V_{h}^{m})_{m \le n}) := \frac{-V^{n-1}}{t_{n} - t_{n-1}}. $$

Moreover, the time discretization must satisfy the following properties:

**Assumption 15.** For all $0 \le n \le N$, and all $0 \le i \le 3$,

(A3) The operators commute with the interpolator, i.e. $\underline{I}_{r,h}^{i} \mathcal{T}_{B}^{n} = \mathcal{T}_{h,B}^{n} \underline{I}_{r,h}^{i}$ and $\underline{I}_{r,h}^{i} \mathcal{T}_{L}^{n} = \mathcal{T}_{h,L}^{n} \underline{I}_{r,h}^{i}$.

(A4) The operator commutes with the spatial derivatives: $d\mathcal{T}^n = \mathcal{T}^n d$ and $d^*\mathcal{T}^n = \mathcal{T}^n d^*$.

(A5) There exists $\theta_n > 0$, such that for all $V_h \in \underline{X}_{r,h}^i$, $(\mathcal{T}_{h,B}^n V_h, V_h)_h \ge \theta_n \|V_h\|_h^2$. Moreover, we must have $\sum_{n=1}^N \frac{1}{\theta_n} \le CT$ for some constant $C$ independent of $N$.

(A6) The discrete operator commutes with $P_{\mathfrak{S}}: P_{\mathfrak{S}}\mathcal{T}_{h}^{n} = \mathcal{T}_{h}^{n}P_{\mathfrak{S}}$, and the continuous operator is zeroth order consistent: if $W$ is constant in time then $\mathcal{T}^{n}W = 0$.

The time discretization operators are defined on each space $\underline{X}_{r,h}^i$ and extend to $\underline{X}_{r,h}$ diagonally. In order to quantify the error on the discrete solution, we introduce the following notations:

$$ \epsilon_{\Delta t, n}(V) := \mathcal{T}^n V - \partial_t V(t_n), \quad (14) $$

$$ \epsilon_h(V, W_h) := (\underline{I}_{r,h}V, d_hW_h)_h - (\underline{I}_{r,h}d^*V, W_h)_h. \quad (15) $$

The discrete problem is: Given $U_h^0 \in \underline{X}_{r,h}^i$, find $(U_h^n)_{1 \le n \le N} \in (\underline{X}_{r,h}^i)^N$ such that for all $1 \le n \le N$, and all $V_h \in \underline{X}_{r,h}^i$,

$$ a_h^n(U_h^n, V_h) = l_h^n(V_h), \quad (16) $$

where,

$$ a_h^n(U_h^n, V_h) := (\mathcal{T}_{h,B}^n U_h^n, V_h) - \mathcal{L}_h(U_h^n, V_h), \quad l_h^n(V_h) := -(\mathcal{T}_{h,L}^n(U_h^n), V_h)_h. $$

**Theorem 16 (Error estimate).** There exists a unique solution $(U_h^n)$ to the problem (16). Moreover, denoting by $U$ the solution of the continuous problem (9), if for all $1 \le n \le N$, $(V_h^n)_{n \le N} \in (\underline{X}_{r,h}^i)^N$,

$$ \begin{align*} & \|\mathcal{T}_{h,L}^n(V_h^m)_{m<n}\|_h \le \theta_n \|V_h^{n-1}\|_h, \\ & \|\underline{I}_{r,h} \epsilon_{\Delta t, N}(U)\|_h + \|\underline{I}_{r,h} \mathcal{T}^n(\epsilon_{\Delta t, m}(U))_{m \le n}\|_h \le (\Delta t)^l E(|U|), \\ & \|\epsilon_h(U, \cdot)\|_{1,h} + \|\epsilon_h(\mathcal{T}^n U, \cdot)\|_{1,h} \le h^{r+1} E(|U|), \end{align*} $$

with $r, l \in \mathbb{N}$ the order of convergence in space and time, $h$ the characteristic size of the discrete spaces, and $E(|U|)$ proportional to some semi-norm of $U$, then,

$$ \|U_h^N - \underline{I}_{r,h} U(T)\|_h \le ((\Delta t)^l + h^{r+1})CE(|U|) + \|(P_{\mathfrak{S}}\underline{I}_{r,h} - \underline{I}_{r,h}P_{\mathfrak{S}})U\|_h, $$

for some constant $C$ independent of $\Delta t$ and $h$.

*Proof.* The proof is detailed in Section 4.3 □

*Remark 17.* The last term can be seen as the harmonic gap. In many cases, it can be estimated from other argument. For instance, if the domain is contractible, the only harmonic forms are in the space of 0-forms (or 3-forms if enforcing Dirichlet boundary conditions). Since the 0 and 3-forms components of $U$ are zero, we have $\|(P_{\mathfrak{S}}\underline{I}_{r,h} - \underline{I}_{r,h}P_{\mathfrak{S}})U\|_h = 0$.

*Remark 18.* In order to simplify the proof, we restrict ourselves to first order time discretizations. However, only a simple adaptation is necessary to consider more general time discretizations.
## 4.3 Proof of the error estimate

**Lemma 19 (Coercivity for the L²-norm).** The bilinear form $a_h^n$ is coercive for the norm $|| \cdot ||_h$:

$$ \forall V_h \in \underline{X}_{r,h}, \quad a_h^n(V_h, V_h) \ge \theta_n ||V_h||_h^2. $$

*Proof.* The result stems from the skew-symmetry of the operator $\mathcal{L}$, and of the coercivity (A5) of $\mathcal{T}_{h,B}^n$. Evaluating $a_h^n(V_h, V_h)$, we have:

$$ \begin{aligned} a_h^n(V_h, V_h) &= (\mathcal{T}_{h,B}^n V_h, V_h)_h - (Jd_h V_h, V_h)_h + (V_h, Jd_h V_h)_h \\ &\ge \theta_n ||V_h||_h^2. \end{aligned} $$

□

Let $U$ be the solution of (9). In order to derive the error estimate, we first introduce another problem: For all $0 < n \le N$, find $\hat{U}_h^n \in \underline{X}_{r,h}$, such that for all $V_h \in \underline{X}_{r,h}$,

$$ \mathcal{L}_h((\mathrm{Id} - P_{\mathcal{S}})\hat{U}_h^n, V_h) = ((\mathrm{Id} - P_{\mathcal{S}})\underline{I}_{r,h}\mathcal{T}^n U, V_h)_h, \quad P_{\mathcal{S}}\hat{U}_h^n = \underline{I}_{r,h}P_{\mathcal{S}}U(t_n). \quad (17) $$

*Remark 20.* The difficulty here is that we must use the graph norm to get the error estimates on the spatial discretization, but we also need to sharply control the $L^2$ norm to preserve the convergence in time. Hence the need of this auxiliary problem.

**Lemma 21 (Error estimate from previous steps).** Let $(U_h^n)_n$ be solution of (16), and $(\hat{U}_h^n)_n$ be solution of (17). For any $0 < n \le N$, it holds:

$$ \|\hat{U}_h^n - U_h^n\|_h \le \frac{1}{\theta_n} \|(\mathrm{Id} - P_{\mathcal{S}}) (\mathcal{T}_h^n (\hat{U}_h^m)_{m \le n} - \underline{I}_{r,h} \mathcal{T}^n U)\|_h + \frac{1}{\theta_n} \|\mathcal{T}_{h,L}^n ((\hat{U}_h^m)_{m < n} - (U_h^m)_{m < n})\|_h. $$

*Proof.* Inserting the definition (17) of $\hat{U}_h^n$ into the definition of $a_h^n$, we infer that for all $V_h \in \underline{X}_{r,h}$,

$$ a_h^n(\hat{U}_h^n, V_h) = (\mathcal{T}_{h,B}^n \hat{U}_h^n, V_h)_h - ((\mathrm{Id} - P_{\mathcal{S}}) \underline{I}_{r,h} \mathcal{T}^n U, V_h)_h. $$

Moreover, we infer from Assumption (A6), (10) and (17) that

$$ P_{\mathcal{S}}\mathcal{T}_h^n(\hat{U}_h^m)_{m \le n} = \mathcal{T}_h^n(P_{\mathcal{S}}\hat{U}_h^m)_{m \le n} = \mathcal{T}_h^n(\underline{I}_{r,h}P_{\mathcal{S}}U) = \underline{I}_{r,h}\mathcal{T}^n(P_{\mathcal{S}}U) = 0. $$

Therefore, from the definition (16) of $U_h^n$, we have

$$ \begin{aligned} a_h^n(\hat{U}_h^n - U_h^n, V_h) &= (\mathcal{T}_{h,B}^n \hat{U}_h^n - (\mathrm{Id} - P_{\mathcal{S}}) \underline{I}_{r,h} \mathcal{T}^n U + \mathcal{T}_{h,L}^n (U_h^m)_{m<n}, V_h)_h \\ &= (\mathcal{T}_{h,B}^n \hat{U}_h^n + \mathcal{T}_{h,L}^n (\hat{U}_h^m)_{m<n} - (\mathrm{Id} - P_{\mathcal{S}}) \underline{I}_{r,h} \mathcal{T}^n U, V_h)_h + (\mathcal{T}_{h,L}^n (U_h^m)_{m<n} - \mathcal{T}_{h,L}^n (\hat{U}_h^m)_{m<n}, V_h)_h \\ &= (\mathcal{T}_h^n (\hat{U}_h^m)_{m \le n} - (\mathrm{Id} - P_{\mathcal{S}}) \underline{I}_{r,h} \mathcal{T}^n U, V_h)_h + (\mathcal{T}_{h,L}^n ((\hat{U}_h^m)_{m < n} - (U_h^m)_{m < n}), V_h)_h \\ &= ((\mathrm{Id} - P_{\mathcal{S}}) (\mathcal{T}_h^n (\hat{U}_h^m)_{m \le n} - \underline{I}_{r,h} \mathcal{T}^n U), V_h)_h + (\mathcal{T}_{h,L}^n ((\hat{U}_h^m)_{m < n} - (U_h^m)_{m < n}), V_h)_h. \end{aligned} \quad (18) $$

We conclude from Lemma 19 and (18), writting

$$ \begin{aligned} \|\hat{U}_h^n - U_h^n\|_h^2 &\le \frac{1}{\theta_n} a_h^n(\hat{U}_h^n - U_h^n, \hat{U}_h^n - U_h^n) \\ &\le \frac{1}{\theta_n} \left( \| (\mathrm{Id} - P_{\mathcal{S}}) (\mathcal{T}_h^n (\hat{U}_h^m)_{m \le n} - \underline{I}_{r,h} \mathcal{T}^n U) \|_h + \| \mathcal{T}_{h,L}^n ((\hat{U}_h^m)_{m < n} - (U_h^m)_{m < n}) \|_h \right) \|\hat{U}_h^n - U_h^n\|_h. \end{aligned} $$

□
**Lemma 22** (Error estimate at a single step). Let $(U_h^n)_n$ be solution of (16), and $(\hat{U}_h^n)_n$ be solution of (17). For any $0 < n \le N$, it holds:

$$ \|(\mathrm{Id} - P_S) (\underline{\mathbf{I}}_{r,h} \mathcal{T}^n U - \mathcal{T}_h^n (\hat{U}_h^m)_{m \le n})\|_{1,h} \le \frac{1}{C_L} \left( \|\|\epsilon_h(\mathcal{T}^n U, \cdot)\|\|_{1,h} + \|\underline{\mathbf{I}}_{r,h} \mathcal{T}^n (\epsilon_{\Delta t, m}(U))_{m \le n}\|_h \right), $$

where $C_L$ is the constant given by Lemma 13.

*Proof.* Applying the linearity of the time discretization operator $\mathcal{T}_h^n$ to (17), we have for any $V_h \in \mathfrak{S}^{\perp h}$:

$$ \mathcal{L}_h((\mathrm{Id} - P_S) \mathcal{T}_h^n (\hat{U}_h^m)_{m \le n}, V_h) = (\mathcal{T}_h^n \underline{\mathbf{I}}_{r,h} (\mathcal{T}^m U)_{m \le n}, V_h). \quad (19) $$

Therefore, we have

$$ \begin{align} \mathcal{L}_h((\mathrm{Id} - P_S) (\underline{\mathbf{I}}_{r,h} \mathcal{T}^n U - \mathcal{T}_h^n (\hat{U}_h^m)_{m \le n}), V_h) & \stackrel{(19)}{=} (J \mathrm{d}h \underline{\mathbf{I}}_{r,h} \mathcal{T}^n U, V_h)_h - (\underline{\mathbf{I}}_{r,h} \mathcal{T}^n U, J \mathrm{d}h V_h)_h - (\mathcal{T}_h^n \underline{\mathbf{I}}_{r,h} (\mathcal{T}^m U)_{m \le n}, V_h)_h \nonumber \\ & \stackrel{(15)}{=} (J \mathrm{d}h \underline{\mathbf{I}}_{r,h} \mathcal{T}^n U, V_h)_h - (\underline{\mathbf{I}}_{r,h} \mathrm{d}^* J \mathcal{T}^n U, V_h)_h - \epsilon_h (J \mathcal{T}^n U, V_h) - (\mathcal{T}_h^n \underline{\mathbf{I}}_{r,h} (\mathcal{T}^m U)_{m \le n}, V_h)_h \nonumber \\ & \stackrel{(\mathrm{A2}), (\mathrm{A3})}{=} (\underline{\mathbf{I}}_{r,h} J \mathrm{d} \mathcal{T}^n U, V_h)_h + (\underline{\mathbf{I}}_{r,h} J \mathrm{d}^* \mathcal{T}^n U, V_h)_h - \epsilon_h (J \mathcal{T}^n U, V_h) - (\underline{\mathbf{I}}_{r,h} \mathcal{T}^n (\mathcal{T}^m U)_{m \le n}, V_h)_h \nonumber \\ & \stackrel{(\mathrm{A4}), (14)}{=} (\underline{\mathbf{I}}_{r,h} \mathcal{T}^n (J(\mathrm{d} + \mathrm{d}^*)U - \partial_t U), V_h)_h - \epsilon_h (J \mathcal{T}^n U, V_h) - (\underline{\mathbf{I}}_{r,h} \mathcal{T}^n (\epsilon_{\Delta t, m}(U))_{m \le n}, V_h)_h \nonumber \\ & \le (\|\|\epsilon_h(\mathcal{T}^n U, \cdot)\|\|_{1,h} + \|\underline{\mathbf{I}}_{r,h} \mathcal{T}^n (\epsilon_{\Delta t, m}(U))_{m \le n}\|_h) \|V_h\|_{1,h}, \tag{20} \end{align} $$

where we removed the projector $P_S$ using the fact that $dP_S = 0$ and $d^*P_S = 0$, the cancellation occurs because $U$ is solution of (9), and we used a Cauchy-Schwartz inequality together with $\|V_h\|_h \le \|V_h\|_{1,h}$ on the last line. We conclude evaluating (20) for $V_h$ realizing the inf-sup inequality of Lemma 13. $\square$

**Lemma 23.** Let $\hat{U}_h^N$ be solution of (17), and $C_L$ be the constant given by Lemma 13. It holds

$$ \|(\mathrm{Id} - P_S) (\underline{\mathbf{I}}_{r,h} U(T) - \hat{U}_h^N)\|_{1,h} \le \frac{1}{C_L} \left( \|\|\epsilon_h(U, \cdot)\|\|_{1,h} + \|\underline{\mathbf{I}}_{r,h} \epsilon_{\Delta t, N}(U)\|_h \right). $$

*Proof.* Using the definition of $\hat{U}_h^N$, we have for any $V_h \in \mathfrak{S}^{\perp h}$:

$$ \begin{align} \mathcal{L}_h((\mathrm{Id} - P_S) (\underline{\mathbf{I}}_{r,h} U(T) - \hat{U}_h^N), V_h) &= (J \mathrm{d} \underline{\mathbf{I}}_{r,h} U, V_h)_h - (\underline{\mathbf{I}}_{r,h} U, J \mathrm{d} h V_h)_h - (\underline{\mathbf{I}}_{r,h} \mathcal{T}^N U, V_h)_h \nonumber \\ &\stackrel{(\mathrm{A2})}{=} (\underline{\mathbf{I}}_{r,h} (J(\mathrm{d} + \mathrm{d}^*)U - \partial_t U), V_h)_h - \epsilon_h (J U, V_h) - (\underline{\mathbf{I}}_{r,h} \epsilon_{\Delta t, N} U, V_h)_h \nonumber \\ &\le (\|\|\epsilon_h(U, \cdot)\|\|_{1,h} + \|\underline{\mathbf{I}}_{r,h} \epsilon_{\Delta t, N}(U)\|_{m \le n}\|_h) \|V_h\|_{1,h}, \tag{21} \end{align} $$

where the cancellation occurs because $U$ is solution of (9). We conclude using evaluating (21) for $V_h$ realizing the inf-sup inequality of Lemma 13. $\square$

*Proof of Theorem 16.* The well-posedness of the discrete problem (16) readily follows from the coercivity of $a_h^n$ proven in Lemma 19. Let us now prove the error estimate: Let $U$ be the solution of (9). For any $0 < n \le N$, applying Lemma 21 together with the assumption $\|\mathcal{T}_{h,L}^n ((\hat{U}_h^m)_{m<n} - (U_h^m)_{m<n})\|_h \le \theta_n \|\hat{U}_h^{n-1} - U_h^{n-1}\|_h$ gives:

$$ \|\hat{U}_h^n - U_h^n\|_h \le \frac{1}{\theta_n} \|(\mathrm{Id} - P_S) (\mathcal{T}_h^n (\hat{U}_h^m)_{m \le n} - \underline{\mathbf{I}}_{r,h} \mathcal{T}^n U)\|_h + \|\hat{U}_h^{n-1} - U_h^{n-1}\|_h. $$
Recalling that $U_h^0 = \hat{U}_h^0$, we have

$$ \|\hat{U}_h^N - U_h^N\|_h \le \sum_{n=1}^N \frac{1}{\theta_n} \|(\mathrm{Id} - P_S) (\mathcal{T}_h^n (\hat{U}_h^m)_{m \le n} - \mathbf{I}_{r,h} \mathcal{T}^n U)\|_h. \quad (22) $$

Applying Lemma 22 to (22), we find

$$ \|\hat{U}_h^N - U_h^N\|_h \le \frac{1}{C_L} \sum_{n=1}^N \frac{1}{\theta_n} \left( \|\epsilon_h(\mathcal{T}^n U, \cdot)\|_{1,h} + \|\mathbf{I}_{r,h} \mathcal{T}^n (\epsilon_{\Delta t, m}(U))_{m \le n}\|_h \right) \quad (23) $$

Introducing the assumed bound on $\|\epsilon_h(\mathcal{T}^n U, \cdot)\|_{1,h}$ and $\|\mathbf{I}_{r,h} \mathcal{T}^n (\epsilon_{\Delta t, m}(U))_{m \le n}\|_h$ into (23), we have

$$ \begin{aligned} \|\hat{U}_h^N - U_h^N\|_h &\le \frac{1}{C_L} \sum_{n=1}^N \frac{1}{\theta_n} \left( (\Delta t)^l + h^{r+1} \right) E(|U|) \quad (24) \\ &\le \left( (\Delta t)^l + h^{r+1} \right) \frac{C}{C_L} E(|U|). \end{aligned} $$

We conclude using the triangle inequality:

$$ \|\underline{\mathbf{I}}_{r,h} U(T) - U_h^N\|_h \le \|(\mathrm{Id} - P_S) (\underline{\mathbf{I}}_{r,h} U(T) - \hat{U}_h^N)\|_h + \|P_S (\underline{\mathbf{I}}_{r,h} U(T) - \hat{U}_h^N)\|_h + \|\hat{U}_h^N - U_h^N\|_h, $$

then using Lemma 23 to bound the first term, (17) on the second term, and (24) to bound the last term. □

## 4.4 Examples: a specific scheme

In this section we consider a specific choice of discretization in order to prove more properties of the scheme. We introduce the notation $a \le b$, meaning that there is $C > 0$ depending only on the chosen discrete complex such that $a \le Cb$. We consider a backward Euler time stepping with constant time step $\Delta t$, setting

$$ \mathcal{T}_h^n (V_h^m)_{m \le n} := \frac{V_h^n - V_h^{n-1}}{\Delta t}, \quad \theta_n = \frac{1}{\Delta t}. \quad (25) $$

Let $Y \subset D(d) \cap D(d^*)$ denote a subset of $X$ over which $\underline{\mathbf{I}}_{r,h}$ is continuous for the $L^2$-norm (i.e. $\forall V \in Y$, $\|\underline{\mathbf{I}}_{r,h}V\|_h \lesssim \|V\|_{L^2(\Omega)}$), and let $U$ be the solution of (9).

**Lemma 24 (Backward Euler time stepping).** If $U \in C^3([0, T], Y)$, then for all $1 \le n \le N$, $(V_h^n)_{n \le N} \in (\underline{\mathbf{X}}_{r,h})^N$, it holds

$$ \begin{aligned} \|\mathcal{T}_{h,L}^n (V_h^m)_{m<n}\|_h &\le \theta_n \|V_h^{n-1}\|_h, \\ \|\underline{\mathbf{I}}_{r,h} \epsilon_{\Delta t, N}(U)\|_h + \|\underline{\mathbf{I}}_{r,h} \mathcal{T}^n (\epsilon_{\Delta t, m}(U))_{m \le n}\|_h &\lesssim \Delta t (|U|_{C^2([0,T],Y)} + |U|_{C^3([0,T],Y)}). \end{aligned} $$

*Proof.* The first bound is trivial from the definition (25) of $\mathcal{T}_{h,L}^n$ giving

$$ \|\mathcal{T}_{h,L}^n (V_h^m)_{m<n}\|_h = \frac{\|V_h^{n-1}\|_h}{\Delta t} = \theta_n \|V_h^{n-1}\|_h. $$

Assuming the $C^3$ regularity in time of $U$, we can write its Taylor expansion at $t_n$ for $n \ge 2$:

$$ \begin{aligned} U(t_{n-1}) &= U(t_n) - \Delta t \partial_t U(t_n) + (\Delta t)^2 \frac{1}{2} \partial_t^2 U(t_n) - (\Delta t)^3 \frac{1}{6} \partial_t^3 U(c) \\ \partial_t^2 U(t_{n-1}) &= \partial_t^2 U(t_n) - \Delta t \partial_t^3 U(c), \end{aligned} \quad (26) $$
where $c \in [t_{n-1}, t_n]$. Injecting (26) into the definition of $\epsilon_{\Delta t, n}(U)$ gives:

$$ \epsilon_{\Delta t, n}(U) = \frac{U(t_n) - U(t_{n-1})}{\Delta t} - \partial_t U(t_n) = -\Delta t \frac{1}{2} \partial_t^2 U(t_n) + (\Delta t)^2 \frac{1}{6} \partial_t^3 U(c). $$

Using the same formula for $\epsilon_{\Delta t, n-1}(U)$, we find

$$ \begin{aligned} \mathcal{T}^n(\epsilon_{\Delta t, m}(U))_{m \le n} &= \frac{\epsilon_{\Delta t, n}(U) - \epsilon_{\Delta t, n-1}(U)}{\Delta t} \\ &= -\Delta t \frac{1}{2} \frac{\partial_t^2 U(t_n) - \partial_t^2 U(t_{n-1})}{\Delta t} + \Delta t \frac{1}{6} (\partial_t^3 U(c') - \partial_t^3 U(c'')) \\ &= -\Delta t \frac{1}{2} \partial_t^3 U(c) + \Delta t \frac{1}{6} (\partial_t^3 U(c') - \partial_t^3 U(c'')). \end{aligned} $$

Taking the $L^2$ norm on both side and bounding $\|\partial_t^2 U(c)\|_{L^2(\Omega)}$ by $\|U\|_{C^2([0,T],Y)}$, we have

$$ \|\epsilon_{\Delta t, n}(U)\|_{L^2(\Omega)} \lesssim \Delta t \|U\|_{C^2([0,T],Y)}, \quad \|\mathcal{T}^n(\epsilon_{\Delta t, m}(U))_{m \le n}\|_{L^2(\Omega)} \lesssim \Delta t \|U\|_{C^3([0,T],Y)}. $$

We infer the result from the continuity of $I_{r,h}$ on $Y$. $\square$

**Remark 25.** The result is straightforward to extend to higher-order schemes approximating the time derivative at $t_n$. Other schemes such as the Crank-Nicolson time stepping require some slight modification. For instance, since the Crank-Nicolson time stepping approximates the time derivative at $t_{n-\frac{1}{2}} := \frac{t_n+t_{n-1}}{2}$, to preserve the second order accuracy, we should use $\epsilon_{\Delta t, n-\frac{1}{2}}(U) := \frac{U(t_n)-U(t_{n-1})}{\Delta t} - \partial_t U(t_{n-\frac{1}{2}})$. Then we can show that the error

$$ \mathcal{T}^n(\epsilon_{\Delta t, m-\frac{1}{2}}(U))_{m \le n} = \frac{\epsilon_{\Delta t, n-\frac{1}{2}}(U) - \epsilon_{\Delta t, n-\frac{3}{2}}(U)}{\Delta t} $$

is second order accurate.

# 5 Numerical results

## 5.1 Discrete complex

We use a conforming divdiv complex based on tensor product of splines [8]. The main ingredients in its construction are one dimensional finite elements. Let $\mathcal{E}_h$ and $\mathcal{V}_h$ denote respectively the set of edges and vertices dividing a segment $[0, 1]$. We consider the following spaces:

$$ \bullet S^1 := \{v \in C^1([0, 1]) : \forall E \in \mathcal{E}_h, v_E \in \mathcal{P}_3(E)\} $$

$$ \bullet S^0 := \{v \in C^0([0, 1]) : \forall E \in \mathcal{E}_h, v_E \in \mathcal{P}_2(E)\} $$

$$ \bullet S^{-1} := \{v \in L^2([0, 1]) : \forall E \in \mathcal{E}_h, v_E \in \mathcal{P}_1(E)\} $$

The associated interpolator are:

$$ \bullet I^1 \text{ such that } \forall v \in C^1([0, 1]), \forall V \in \mathcal{V}_h, I^1(v)(V) = v(V), (I^1(v))'(V) = v'(V). $$

$$ \bullet I^0 \text{ such that } \forall v \in C^0([0, 1]), \forall V \in \mathcal{V}_h, I^0(v)(V) = v(V), \forall E \in \mathcal{E}_h, \int_E I^0(v) = \int_E v. $$

$$ \bullet I^{-1} \text{ such that } \forall v \in L^2([0, 1]), \forall E \in \mathcal{E}_h, \forall p \in \mathcal{P}_1(E), \int_E p I^{-1}(v) = \int_E pv. $$

The basis for the discrete complex are the tensor product of these spaces. We denote by $S^{a,b,c} := S^a \otimes S^b \otimes S^c$. The spaces of the discrete complex are:

$$ \underline{\mathbf{X}}_{r,h}^0 := \begin{pmatrix} S^{1,0,0} \\ S^{0,1,0} \\ S^{0,0,1} \end{pmatrix}, \quad \underline{\mathbf{X}}_{r,h}^1 := \begin{pmatrix} S^{0,0,0} & S^{1,-1,0} & S^{1,0,-1} \\ S^{-1,1,0} & S^{0,0,0} & S^{0,1,-1} \\ S^{-1,0,1} & S^{0,-1,1} & S^{0,0,0} \end{pmatrix}, $$

$$ \underline{\mathbf{X}}_{r,h}^2 := \begin{pmatrix} S^{1,-1,-1} & S^{0,0,-1} & S^{0,-1,0} \\ S^{0,0,-1} & S^{-1,1,-1} & S^{-1,0,0} \\ S^{0,-1,0} & S^{-1,0,0} & S^{-1,-1,1} \end{pmatrix}, \quad \underline{\mathbf{X}}_{r,h}^3 := S^{-1,-1,-1}. $$
All the spaces are conforming, hence we can take the restriction of the continuous differential as the discrete differential [8].

## 5.2 Wave-like solutions

In order to numerically validate our scheme, we consider two classes of wave-like solutions. The first depends on three parameters $\lambda_1, \lambda_2, c \in \mathbb{R}$. We define

$$ k := \begin{pmatrix} 0 \\ 0 \\ c \end{pmatrix}, \quad A := \begin{pmatrix} 0 & \lambda_1 & 0 \\ -\lambda_2 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} e^{i(k \cdot x - ct)}, \quad \gamma := \begin{pmatrix} \frac{\lambda_1 - \lambda_2}{2} & 0 & 0 \\ 0 & \frac{\lambda_2 - \lambda_1}{2} & 0 \\ 0 & 0 & -\frac{\lambda_1 + \lambda_2}{2} \end{pmatrix} e^{i(k \cdot x - ct)}. \quad (27) $$

The second class depends on two parameters $\lambda, c \in \mathbb{R}$. We define

$$ k := \begin{pmatrix} c \\ c \\ 0 \end{pmatrix}, \quad A := \begin{pmatrix} 0 & 0 & -\lambda \\ 0 & 0 & -\lambda \\ 0 & 0 & 0 \end{pmatrix} e^{i(k \cdot x - ct)}, \quad \gamma := \begin{pmatrix} \lambda & 0 & 0 \\ 0 & -\lambda & 0 \\ 0 & 0 & 0 \end{pmatrix} e^{i(k \cdot x - ct)}. \quad (28) $$

Notice that in the second case $k \cdot k = 2c^2$.

We readily verify that both (27) and (28) satisfy (8).

## 5.3 Results

We have implemented two test cases: The first one is given by (27) with $c := \pi$, $\lambda_1 := 2$, $\lambda_2 := 1$, and the second one is given by (28) with $c := \pi$ and $\lambda := 2$. In either case, we initialized the solution with the reference one at $t = 0$, and let the system evolve until $t = 2\pi$. The domain consists of a unit cube, and the mesh is a Cartesian grid. We enforced a Dirichlet boundary condition on the whole boundary, deriving the value from the analytical solution. With the setting, the only harmonic forms are the 3-forms spanned by the linear polynomials. Since the component in the space of 3-forms of our analytical solutions corresponds to $\lambda_3$, and is taken to be zero, the $\|(P_S I_{r,h} - I_{r,h} P_S)U\|_h$ term vanishes in Theorem 16. We explored various time steps and spatial subdivisions, ranging between $\Delta t = 10^{-2}$ and $\Delta t = 10^{-5}$ for the time steps, and between $N = 2$ to $N = 10$ subdivision (hence between 8 and 1000 cells). In order to save on resources, we did not run the simulation for all possible combinations.

The error computed for the various fields with respect to the cell size $h$ using the backward Euler time stepping is given in Figure 1. The value showed for the error is $L^1([0, 2\pi]) \times L^2([0, 1]^3)$ norm of the local error,

$$ E := \int_{t=0}^{2\pi} \left( \int_{x \in [0,1]^3} \|u_h - I_{r,h}(u)\|^2 \right)^{\frac{1}{2}} . $$

We notice that the error attributed to the spatial discretization converges quickly below the contribution of the temporal discretization. The error on $\lambda_0$ is initially of the order of machine precision and grows due to accumulation errors (between $10^{-13}$ and $10^{-7}$). The evolution of the error during the simulation is given in Figure 3 for a case dominated by the time discretization, one dominated by the spatial discretization, and one intermediary.

Since we could not neglect the error due to the time discretization, we modeled our error as $E = \gamma ((\Delta t)^\alpha + \delta h^\beta)$, and computed the convergence rates $\alpha$ and $\beta$ as a best fit for the data obtained in our simulation. The results obtained are shown in Table 2. The results are consistent with a first order time discretization. Since there exists bounded cochain projections to the spline complex [8], we expect to see a spatial convergence of order $k+1$ where $k$ is the polynomial degree. The discrete spaces for the space of 1, 2, and 3-forms contain all polynomial of degree 1, but not higher in a certain direction. Hence, the standard convergence theory gives a convergence rate of 2. The computed convergence rate is compatible with the result, even significantly higher. Two possible explanations for this higher convergence rate are
<table><caption>Table 2: Computed convergence rates using the Euler time stepping.</caption><thead><tr><th rowspan="2"></th><th colspan="2">First case</th><th colspan="2">Second case</th></tr><tr><th>α</th><th>β</th><th>α</th><th>β</th></tr></thead><tbody><tr><th>g</th><td>0.986</td><td>3.31</td><td>0.987</td><td>3.20</td></tr><tr><th>A</th><td>0.990</td><td>3.13</td><td>0.996</td><td>3.05</td></tr><tr><th>λ<sub>3</sub></th><td>1.01</td><td>4.40</td><td>1.06</td><td>4.11</td></tr></tbody></table>

<table><caption>Table 3: Computed convergence rates (in space) using the Crank-Nicolson time stepping.</caption><thead><tr><th></th><th>First case</th><th>Second case</th></tr></thead><tbody><tr><th>g</th><td>3.22</td><td>3.12</td></tr><tr><th>A</th><td>3.09</td><td>3.03</td></tr><tr><th>λ<sub>3</sub></th><td>4.16</td><td>4.07</td></tr></tbody></table>

the alignment between our analytical solutions, and the directions containing higher order polynomials in the discrete spaces, or a superconvergence due to the use of a smooth solution, and the use of the discrete norm to compute the error.

To confirm the convergence rates, we also implemented the Crank-Nicolson time stepping that is second order accurate in time. We compute the error using the same formula on the same test cases. The results are given in Figure 2. Since the contributions to the error from the time discretization are neglectable, we computed the convergence rate in space only using a timestep of $\Delta t = 10^{-4}$. The results are shown in Table 3. They are coherent with the value obtained using the backward Euler time stepping.

# 6 Conclusions

We aim to extend this stability to the full, nonlinear equations. The linearized version provides a promising first step, since the mathematical properties of nonlinear hyperbolic systems mainly rely on their principal parts. However, deriving the associated nonlinear mixed formulation is beyond the scope of this paper and is left as a future direction.

## A Vector calculus identities

**Lemma 26.** If $\gamma \in \mathbb{S}$ then

$$ \mathrm{sym} \, \mathrm{curl} \, \mathrm{curl} \, S\gamma = \mathrm{inc} \, \gamma - S \, \mathrm{def} \, \mathrm{div} \, S\gamma - S \, \mathrm{hess} \, \mathrm{tr} \, \gamma. \qquad (29) $$

*Proof.* The proof relies on the identities

$$ \begin{gathered} \mathrm{curl} \, \iota = \mathrm{mskw} \, \mathrm{grad} \\ 2 \, \mathrm{skw} \, \mathrm{curl} = \mathrm{mskw} \, 2 \, \mathrm{vskw} \, \mathrm{curl} = - \mathrm{mskw} \, \mathrm{div} \, S, \end{gathered} \qquad (30) $$

and

$$ \begin{gathered} \mathrm{curl} \, \mathrm{mskw} = S \, \mathrm{grad} \\ \mathrm{sym} \, \mathrm{curl} \, \mathrm{mskw} = S \, \mathrm{def}. \end{gathered} \qquad (31) $$

If $\gamma \in \mathbb{S}$, then $\mathrm{curl} \, \gamma \in \mathbb{T}$ and $S^{-1} \mathrm{curl} \, \gamma = T \mathrm{curl} \, \gamma = \mathrm{curl} \, \gamma - 2 \, \mathrm{skw} \, \mathrm{curl} \, \gamma$. Expanding the
<figure><img src="image_3.png" /><figcaption>(a) Error on A in the first case.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) Error on A in the second case.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) Error on g in the first case.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(d) Error on g in the second case.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(e) Norm of $\lambda_3$ in the first case.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(f) Norm of $\lambda_3$ in the second case.</figcaption></figure>

<caption>Figure 1: Error with respect to the mesh size for various time steps using the Euler time stepping.</caption>
<figure><img src="image_9.png" /><figcaption>(a) Error on A in the first case.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Error on A in the second case.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c) Error on g in the first case.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(d) Error on g in the second case.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(e) Norm of $\lambda_3$ in the first case.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(f) Norm of $\lambda_3$ in the second case.</figcaption></figure>

Figure 2: Error with respect to the mesh size for various time steps using the Crank-Nicolson time stepping.
<figure><img src="image_15.png" /><figcaption>(a) Error dominated by the time approximation.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) Error with comparable contribution from time and space.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(c) Error dominated by the spatial approximation.</figcaption></figure>

<caption>Figure 3: Evolution of the error during the simulated time for the first case.</caption>
definition of S, we have

$$
\begin{align*}
\operatorname{sym} \operatorname{curl} \operatorname{curl} S\gamma &= \operatorname{sym} \operatorname{curl} (\operatorname{curl} \gamma - \operatorname{curl} \iota \operatorname{tr} \gamma) \\
&= \operatorname{sym} \operatorname{curl} (S^{-1} \operatorname{curl} \gamma + 2 \operatorname{skw} \operatorname{curl} \gamma - \operatorname{curl} \iota \operatorname{tr} \gamma) \\
&\stackrel{(30)}{=} \operatorname{sym} \operatorname{curl} (S^{-1} \operatorname{curl} \gamma - \operatorname{mskw} \operatorname{div} S\gamma - \operatorname{mskw} \operatorname{grad} \operatorname{tr} \gamma) \\
&\stackrel{(31)}{=} \operatorname{sym} \operatorname{inc} \gamma - S \operatorname{def} \operatorname{div} S\gamma - \operatorname{sym} S \operatorname{hess} \operatorname{tr} \gamma \\
&= \operatorname{inc} \gamma - S \operatorname{def} \operatorname{div} S\gamma - S \operatorname{hess} \operatorname{tr} \gamma.
\end{align*}
$$

□

**Lemma 27.** If div A = 0 then

$$
\operatorname{div} \operatorname{sym} \operatorname{curl} A = - \operatorname{curl} \operatorname{curl} \operatorname{vskw} A. \tag{32}
$$

*Proof.* The proof relies on the identity

$$
\operatorname{div} \operatorname{mskw} = \operatorname{curl} . \tag{33}
$$

Expanding the definition of sym, we have

$$
\begin{align*}
\operatorname{div} \operatorname{sym} \operatorname{curl} A &= \frac{1}{2} \operatorname{div} \operatorname{curl} A + \frac{1}{2} \operatorname{div} T \operatorname{curl} A \\
&\stackrel{(35)}{=} \frac{1}{2} \operatorname{curl} \operatorname{div} A^\top \\
&= \frac{1}{2} \operatorname{curl} \operatorname{div} A - \operatorname{curl} \operatorname{div} \operatorname{skw} A \\
&= - \operatorname{curl} \operatorname{div} \operatorname{mskw} \operatorname{vskw} A \\
&\stackrel{(33)}{=} - \operatorname{curl} \operatorname{curl} \operatorname{vskw} A,
\end{align*}
$$

where we used the decomposition $A^\top = A - 2 \operatorname{skw} A$ on the third line.

□

**Lemma 28.** In general, it holds

$$
\operatorname{curl} \operatorname{curl} = -2 \operatorname{div} S \operatorname{def}, \quad (34)
$$

$$
\operatorname{div} T \operatorname{curl} = \operatorname{curl} \operatorname{div} T, \tag{35}
$$

where *T* is the transpose operator.

*Proof.* The result is a straightforward computation.

Proof of (34): Using Einstein notations, for any vector field v, we have (curl v)_i = ε_{ijk} ∂_j v_k where ε is the fully skew-symmetric tensor. Therefore,

$$
(\mathrm{curl} \ \mathrm{curl} \ v)_i = \epsilon_{ijk} \partial_j (\mathrm{curl} \ v)_k = \epsilon_{ijk} \epsilon_{klm} \partial_j \partial_l v_m = (\delta_i^l \delta_j^m - \delta_i^m \delta_j^l) \partial_j \partial_l v_m = \partial_j \partial_i v_j - \partial_j \partial_i v_i, \quad (36)
$$

which is the vector Laplacian identity $-\Delta = \text{curl curl} - \text{grad div}$. On the other hand, writing $S \text{ def } v = \text{def } v - \iota \text{ tr } \text{grad } v$, we have

$$
2(S \operatorname{def} v)_{ji} = 2(\operatorname{def} v)_{ji} - 2(\iota \operatorname{tr} \operatorname{grad} v)_{ji} = \partial_j v_i + \partial_i v_j - 2\delta_j^i \partial_k v_k.
$$
We infer (34) equating (36) with

$$-2(\operatorname{div} S \operatorname{def} v)_i = -\partial_j (2S \operatorname{def} v)_{ji} = 2\partial_i \partial_j v_j - \partial_j \partial_i v_j - \partial_j \partial_i v_i = \partial_i \partial_j v_j - \partial_j \partial_i v_i.$$

**Proof of (35):** For any matrix field $M$, we have by convention $(\operatorname{curl} M)_{ij} = \epsilon_{ikl}\partial_k M_{lj}$. Thus

$$(\operatorname{div} T \operatorname{curl} M)_i = \partial_j (T \operatorname{curl} M)_{ji} = \partial_j (\operatorname{curl} M)_{ij} = \epsilon_{ikl} \partial_j \partial_k M_{lj}. \qquad (37)$$

On the other hand, we have $(\operatorname{div} M^T)_i = \partial_j (M^T)_{ji} = \partial_j M_{ij}$. We infer (35) equating (37) with

$$(\operatorname{curl} \operatorname{div} M^T)_i = \epsilon_{ikl} \partial_k (\operatorname{div} M^T)_l = \epsilon_{ikl} \partial_k \partial_j M_{lj}.$$

□

## Acknowledgement

The work was supported by a Royal Society University Research Fellowship (URF\R1\221398, RF\ERE\221047), an ERC Starting Grant (project 101164551, GeoFEM) and a Royal Society International Exchanges Grant (IEC\NSFC\233594). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

## References

[1] M. Alcubierre. *Introduction to 3+1 Numerical Relativity*. Oxford University Press, Apr. 2008. ISBN: 9780199205677. DOI: 10.1093/acprof:oso/9780199205677.001.0001.

[2] D. Arnold. *Finite Element Exterior Calculus*. SIAM, 2018. DOI: 10.1137/1.9781611975543.

[3] D. Arnold and K. Hu. "Complexes from Complexes". In: *Foundations of Computational Mathematics* 6 (2021), pp. 1739–1774. DOI: 10.1007/s10208-021-09498-9.

[4] D. Arnold, R. Falk, and R. Winther. "Mixed finite element methods for linear elasticity with weakly imposed symmetry". In: *Mathematics of Computation* 76.260 (2007), pp. 1699–1723.

[5] D. N. Arnold, R. S. Falk, and R. Winther. "Finite element exterior calculus, homological techniques, and applications". In: *Acta Numer.* 15 (2006), pp. 1–155. DOI: 10.1017/S0962492906210018.

[6] D. N. Arnold, R. S. Falk, and R. Winther. "Finite element exterior calculus: from Hodge theory to numerical stability". In: *Bull. Amer. Math. Soc.* (N.S.) 47.2 (2010), pp. 281–354. DOI: 10.1090/S0273-0979-10-01278-4.

[7] Y. Berchenko-Kogan and E. S. Gawlik. "Finite element spaces of double forms". In: *arXiv preprint arXiv:2505.17243* (2025).

[8] F. Bonizzoni, K. Hu, G. Kanschat, and D. Sap. "Discrete tensor product BGG sequences: Splines and finite elements". In: *Mathematics of Computation* 94.352 (2025), pp. 517–549.

[9] D. Braess and J. Schöberl. "Equilibrated residual error estimator for edge elements". In: *Math. Comp.* 77.262 (2008), pp. 651–672. DOI: 10.1090/S0025-5718-07-02080-7.

[10] A. Čap and K. Hu. "BGG sequences with weak regularity and applications". In: *Foundations of Computational Mathematics* (2023), pp. 1–40.

[11] T. Cazenave and A. Haraux. *An Introduction to Semilinear Evolution Equations*. Oxford University Press, 1998.

[12] L. Chen and X. Huang. *Finite elements for divdiv-conforming symmetric tensors*. 2020. arXiv: 2005.01271v3 [math.NA].
[13] S. H. Christiansen. "On the linearization of Regge calculus". In: *Numerische Mathematik* 119.4 (2011), pp. 613–640.

[14] J. Gopalakrishnan, K. Hu, and J. Schöberl. "A 2-complex containing Sobolev spaces of matrix fields". In: *arXiv preprint arXiv:2507.11869* (2025).

[15] J. Gopalakrishnan, P. L. Lederer, and J. Schöberl. "A mass conserving mixed stress formulation for the Stokes equations". In: *IMA Journal of Numerical Analysis* 40.3 (2020), pp. 1838–1874.

[16] Y. Guo, J. Hu, and T. Lin. "Discretizing linearized Einstein-Bianchi system by symmetric and traceless tensors". In: *arXiv preprint arXiv:2508.04560* (2025).

[17] J. Hu, Y. Liang, and R. Ma. *Conforming finite element DIVDIV complexes and the application for the linearized Einstein-Bianchi system*. 2021. arXiv: 2103.00088 [math.NA].

[18] J. Hu, Y. Liang, R. Ma, and M. Zhang. "A family of conforming finite element divdiv complexes on cuboid meshes". In: *Numerische Mathematik* 156.4 (2024), pp. 1603–1638.

[19] K. Hu and T. Lin. "Finite element form-valued forms: Construction". In: *arXiv preprint arXiv:2503.03243* (2025).

[20] K. Hu, T. Lin, and B. Shi. "Finite elements for symmetric and traceless tensors in three dimensions". In: *arXiv preprint arXiv:2311.16077* (2023).

[21] K. Hu, T. Lin, and Q. Zhang. "Distributional Hessian and divdiv complexes on triangulation and cohomology". In: *SIAM Journal on Applied Algebra and Geometry* 9.1 (2025), pp. 108–153.

[22] X. Huang. "Finite element conformal complexes in three dimensions". In: *arXiv preprint arXiv:2508.01238* (2025).

[23] P. Leopardi and A. Stern. "The abstract Hodge–Dirac operator and its stable discretization". In: *SIAM J. Numer. Anal.* 54.6 (2016), pp. 3258–3279. DOI: 10.1137/15m1047684. URL: https://doi.org/10.1137/15m1047684.

[24] L. Li. "Regge finite elements with applications in solid mechanics and relativity". PhD thesis. University of Minnesota, 2018.

[25] M. W. Licht. "Complexes of discrete distributional differential forms and their homology theory". In: *Found. Comput. Math.* 17.4 (2017), pp. 1085–1122. DOI: 10.1007/s10208-016-9315-y.

[26] C. Llinares. "Simulation techniques for modified gravity". In: *International Journal of Modern Physics D* 27.15 (2018), p. 1848003.

[27] M. Neunteufel and J. Schöberl. "The Hellan–Herrmann–Johnson and TDNNS methods for linear and nonlinear shells". In: *Computers & Structures* 305 (2024), p. 107543.

[28] T. A. Oliynyk and J. J. Qian. "A polytopal discrete de Rham scheme for the exterior calculus Einstein's equations". In: *arXiv preprint arXiv:2505.00286* (2025).

[29] A. Pechstein and J. Schöberl. "Tangential-displacement and normal-normal-stress continuous mixed finite elements for elasticity". In: *Mathematical Models and Methods in Applied Sciences* 21.08 (2011), pp. 1761–1782.

[30] A. Petrov, J. R. Nascimento, and P. Porfirio. *Introduction to modified gravity*. Springer Nature, 2023.

[31] D. A. Di-Pietro and M.-L. Hanot. *A discrete three-dimensional divdiv complex on polyhedral meshes with application to a mixed formulation of the biharmonic problem*. 2023. arXiv: 2305.05729 [math.NA].

[32] F. Pretorius. "Evolution of binary black-hole spacetimes". In: *Physical review letters* 95.12 (2005), p. 121101.
[33] V. Quenneville-Belair. "A new approach to finite element simulations of general relativity". In: (2015).

[34] V. Quenneville-Bélair. "A New Approach to Finite Element Simulations of General Relativity". PhD thesis. University of Minnesota, 2015.

[35] J. W. York. "Kinematics and Dynamics of General Relativity". In: *Sources of Gravitational Radiation*. Ed. by L. L. Smarr. Cambridge University Press, 1979, pp. 83–126.