The unique continuation problem for the wave equation
discretized with a high-order space-time nonconforming
method

Erik Burman* Guillaume Delay† Alexandre Ern‡

March 19, 2025

Abstract

We are interested in solving the unique continuation problem for the wave equation, i.e., we want to reconstruct the solution of the wave equation given its (noised) value in a subset of the computational domain. Homogeneous Dirichlet boundary conditions are imposed, whereas the initial datum is unknown. We discretize this problem using a space-time discontinuous Galerkin method (including hybrid variables in space and in time) and look for the solution corresponding to the saddle-point of a discrete Lagrangian. We establish discrete inf-sup stability and bound the consistency error, leading to a priori estimates on the residual. Our main result proves the convergence of the discrete solution to the exact solution in a shifted energy norm involving weaker Sobolev norms than the standard energy norm for the wave equation. The proof combines the above a priori bound with a conditional stability estimate at the continuous level. Finally, we run numerical simulations to assess the performance of the method in practice. A static condensation procedure is used to eliminate the cell unknowns and reduce the size of the linear system.

**Keywords:** unique continuation, data assimilation, wave equation, hybridized discontinuous Galerkin, error estimate, static condensation

# 1 Introduction

In the present work, we are interested in solving numerically a data assimilation problem subject to the wave equation. Homogeneous Dirichlet boundary conditions are imposed, whereas initial conditions are unknown. In order to compensate for the lack of initial datum, we use the knowledge of the solution in a subdomain. We also investigate the influence of noise on this additional datum. Specifically, we consider a bounded open domain $\Omega \subset \mathbb{R}^d$, $d \in \{1, 2, 3\}$, with a smooth boundary $\partial\Omega$, a subset $\varpi \subset \Omega$, and a time interval $J := (0, T_f)$ with final time $T_f > 0$. The assumption that $\partial\Omega$ is smooth has its

<footnote>*Department of Mathematics, University College London, London, UK-WC1E 6BT, UK.</footnote>

<footnote>†Sorbonne Université, Université Paris Cité, CNRS, Inria, Laboratoire Jacques-Louis Lions (LJLL), Paris, France.</footnote>

<footnote>‡CERMICS, ENPC, Institut Polytechnique de Paris, 6&8 av. Blaise Pascal, 77455 Marne-la-Vallée cedex 2, and Inria Paris, 48, rue Barrault, CS 61534, 75647 Paris Cedex France.</footnote>
limits in practice, owing to the hybrid nature of the discretization method; see Remark 3.1. We denote by $Q := J \times \Omega$ the computational domain. Our goal is to approximate the function $u: Q \to \mathbb{R}$ that satisfies

$$L(u) := \partial_{tt}^2 u - \Delta u = f \quad \text{in } Q, \qquad (1a)$$

$$u = 0 \quad \text{on } J \times \partial\Omega, \qquad (1b)$$

$$u = g \quad \text{in } J \times \bar{\omega}, \qquad (1c)$$

where $f \in L^2(Q)$ is a given source term and $g \in H^1(J; L^2(\bar{\omega})) \cap L^2(J; H^1(\bar{\omega}))$ is the restriction to $J \times \bar{\omega}$ of a solution to the wave equation in $Q$. Notice that $f$ and $g$ have to be chosen in a compatible way such that there exists a solution to the problem (1). In the sequel, we consider perturbed data $g_\delta := g + \delta$ to account for some noise $\delta \in L^2(J \times \bar{\omega})$ in the measurements.

This work is carried out under the so-called assumption of geometric control condition.

**Assumption 1 (Geometric control condition).** Every compressed generalized bicharacter-istic intersects the set $J \times \bar{\omega}$ when projected to $Q$.

Roughly this means that every ray has to cross the domain where measurements are available (accounting for reflexions on the computational domain boundary). We refer to [1] for a technical definition of a compressed generalized bicharacteristic, see also [28] where a time-dependent measurement domain is considered. This condition is necessary and sufficient for the exact controllability of the wave equation by a boundary control, see [19]. Notice that the assumption that $\partial\Omega$ is smooth is used to fit the setting of [28].

We consider the following shifted energy norm:

$$\|v\|_{\text{sft}} := \{\|v\|_{L^{\infty}(J;L^2(\Omega))}^2 + \|\nabla v\|_{H^{-1}(J;L^2(\Omega))}^2 + \|\partial_t v\|_{L^2(J;H^{-1}(\Omega))}^2\}^{\frac{1}{2}}. \quad (2)$$

The terminology is motivated by the use of a weaker Sobolev norm than the standard energy norm for the wave equation, which is $\{\|\nabla v\|_{L^2(J;L^2(\Omega))}^2 + \|\partial_t v\|_{L^2(J;L^2(\Omega))}^2\}^{\frac{1}{2}}$. Under the above geometric control condition, we have the following stability estimate.

**Lemma 2 (Stability estimate).** Let the geometric control condition (Assumption 1) be satisfied. There exists $C_{\text{stb}} > 0$ such that for all $v \in H^1(J; H^{-1}(\Omega)) \cap H^{-1}(J; H_0^1(\Omega)) \cap L^\infty(J; L^2(\Omega))$, we have

$$\|v\|_{\text{sft}} \le C_{\text{stb}} (\|v\|_{L^2(J \times \bar{\omega})} + \|L(v)\|_{H^{-1}(Q)}). \qquad (3)$$

The proof of this lemma is given in Appendix A. It is one of the building blocks to prove the convergence of our approximation method.

Several works already deal with the numerical resolution of data assimilation subject to a transient PDE. For instance, data assimilation for the heat and wave equations in 1d has already been considered in [3] using a quasi-reversibility method. The case of the wave equation in a multi-dimensional configuration is addressed in [21], where a mixed formula- tion is proposed and proven to be well-posed at the continuous level only. A discretization of this mixed problem using a $C^1$ space-time finite element method is proposed as well. The reader is also referred to [2] for the reconstruction of the wave celerity based on Carle- man estimates. Also notice the work [23] where the authors develop a general framework
for least-squares methods to solve ill-posed problems that are conditionally stable. The norm in the least-squares method and the regularization are determined by the components of the conditional stability estimate. The method is applied to data assimilation problems subject to the heat and wave equations.

It is well-known that data assimilation problems and control problems are closely linked. The discretization of the control problem for the wave equation has been studied in [22]. A mixed formulation is proposed and $C^1$-conforming space-time finite elements are used for the discretization. Another work on discretizing the control problem of the wave equation is [13], where a low-order finite element method is considered for space discretization and a low-order finite difference method for time discretization. This work has been extended in [18], where space-time arbitrary high-order finite elements were considered. More examples of control and data assimilation methods for the wave equation can be found in [25].

We are interested in a class of methods that consists in first discretizing the ill-posed continuous problem and then adding stabilization terms at the discrete level. The discrete system to be solved corresponds to the equations characterizing the saddle point of a Lagrangian functional minimizing the solution discrepancy with respect to the data under the constraint of the wave equation. Some stabilization terms are also considered to enforce the well-posedness of the discrete problem. This framework for data assimilation has been developed in [4, 5] for elliptic problems. Several other application problems have been considered since then, such as the Stokes [7] and Helmholtz [11, 15] equations. Unstationary problems have also been considered in this framework; see for instance [8, 10, 17] for the heat equation. The wave equation has been dealt with in [12], where a low-order finite element method for the space discretization and a finite difference for time have been considered. This work has been continued in [16], with an arbitrary high-order space-time finite element method. Notice though that the above framework is computationally challenging since all the space-time unknowns are globally coupled, thereby adding one more dimension compared to a discretization of a well-posed problem by a time-marching scheme. Finally, in the recent work [9], a dG approach in time is combined with a standard FEM approach in space. The wave equation is written in a mixed formulation. Two preconditioning approaches are also proposed and studied numerically.

The main asset of the above framework is that a priori error bounds in a shifted energy norm can be derived. The method to prove these bounds follows the following two steps:

1. stability, consistency, and a priori error bounds in a residual norm are first established;

2. then, using the stability estimate of the continuous problem (see Lemma 2), a priori error bounds in the shifted energy norm are obtained.

In the present work, we follow these two steps. Stability, consistency and a priori error bounds in residual norm are respectively proved in Lemmas 4, 7 and 8. The a priori error bound in the shifted energy norm is proved in Theorems 11 and 12.

The discretization method we propose here is based on a space-time Galerkin discretization. Discontinuous polynomials are considered in time and in space. A distinctive feature to increase computational effectivity is that interfacial elements in space and in
time are added to reduce the size of the linear system through a static condensation pro-
cedure that eliminates the cell unknowns. This is a significant difference with respect
to [17], where the authors have proposed and analyzed a nonconforming discretization
for the data assimilation problem subject to the heat equation without the use of time-
interfacial unknowns, thereby precluding the possibility of performing static condensation
in time. Actually, the numerical tests in [17] emphasized that the method was computa-
tionally expensive. In the present work, we benefit from the static condensation procedure
to reduce the size of the linear problem to solve. Moreover, the present method can use
arbitrary high-order polynomials, and we expect the static condensation procedure to save
more degrees of freedom for higher-order polynomials.

The rest of this work is organized as follows. In Section 2, we present the discretization
method. In Section 3, we perform the error analysis. Finally, in Section 4, we present
some numerical results that corroborate our theoretical error estimates.

2 Problem discretization

In this section, we describe the time and space discretization, and we present the numerical
scheme studied in the present work.

## 2.1 Time-space meshes

We discretize the space-time cylinder $Q = J \times \Omega$ by a tensor-product space-time mesh.

Let $\mathcal{T}_{\tau}^{t}$ be a quasi-uniform mesh of $J$. We denote by $\tau$ the diameter of its largest
cell. We also define the outward time-cell normal $n_{T^t}$ associated with a generic time cell
$T^t \in \mathcal{T}_{\tau}^{t}$ which is conventionally taken to be an open subset of $\mathbb{R}$. The mesh faces in time
(which coincide with the discrete time nodes) are collected in the set $\mathcal{F}_{\tau}^{t}$, and every generic
time face $F^t \in \mathcal{F}_{\tau}^{t}$ is oriented by the fixed unit normal vector $n^t \in \mathbb{R}^1$ positively orienting
the time line.

Let $\mathcal{T}_h^x$ be a shape-regular mesh of $\Omega$, where $h$ refers to the maximal diameter of the
mesh cells. We assume that $\mathcal{T}_h^x$ covers $\Omega$ exactly and that it is fitted to the measurement
subset $\varpi$ (so that, in particular, at least one mesh cell is contained in $\varpi$). In principle,
the mesh cells in $\mathcal{T}_h^x$ can be polyhedra with planar faces in $\mathbb{R}^d$, and hanging nodes are also
possible. However, the analysis below requires the mesh to be such that the underlying
discontinuous polynomial approximation space has a global $H^1$-conforming subspace with
optimal approximation properties. For simplicity, we will therefore restrict the discus-
sion to meshes composed of simplices (one can also readily consider meshes composed of
cuboids). The mesh cells are conventionally taken to be open subsets of $\mathbb{R}^d$, and $n_{T^x} \in \mathbb{R}^d$
denotes the unit outward normal of the generic mesh cell $T^x \in \mathcal{T}_h^x$. The mesh faces in
space are collected in the set $\mathcal{F}_h^x$, and every generic face $F^x \in \mathcal{F}_h^x$ is oriented by a fixed
unit normal vector $n_{F^x} \in \mathbb{R}^d$. To avoid technicalities, we assume henceforth that the mesh
$\mathcal{T}_h^x$ is quasi-uniform. Therefore, we will use $h$ to measure the diameter of any cell in $\mathcal{T}_h^x$
or any face in $\mathcal{F}_h^x$.

The space-time mesh of $Q$ is defined as $Q_{h,\tau} := \{q := T^t \times T^x \mid T^t \in \mathcal{T}_{\tau}^t, T^x \in \mathcal{T}_h^x\}$.
For the sake of simplicity, we assume that there exists a constant $C > 0$ such that $C^{-1}h \le \tau \le Ch$.
We could lift this assumption, but we would have to include the ratio $h/\tau$
in several estimates in the error analysis. With the above assumption, we will use $h$ as the mesh parameter of the space-time mesh instead of $\max(h, \tau)$. Thus, we simply write $Q_h$ instead of $Q_{h,\tau}$ for the space-time mesh. We notice that, for every space-time cell $q := T^t \times T^x \in Q_h$, its boundary $\partial q$ can be decomposed as $\partial q = \partial^t q \cup \partial^x q$ with

$$ \partial^t q := \bigcup_{\substack{F^t \in \mathcal{F}_\tau^t \\ F^t \subset \partial T^t}} \{F^t \times T^x\}, \quad \partial^x q := \bigcup_{\substack{F^x \in \mathcal{F}_h^x \\ F^x \subset \partial T^x}} \{T^t \times F^x\}. \quad (4) $$

We loosely say that $\partial^x q$ collects the space faces of $q$ and that $\partial^t q$ collects the time faces of $q$. We let $n_q$ be the unit outward normal to $q$. We also set $n_{\partial^t q} := n_q|_{\partial^t q} := n^t$ and $n_{\partial^x q} := n_q|_{\partial^x q}$ with $n_{\partial^x q}|_{T^t \times F^x} := n_{F^x}^t$ for all $T^t \times F^x \in \partial^x q$.

The collection of all the faces of $Q_h$ is denoted by $\mathcal{G}_h$, which is partitioned as $\mathcal{G}_h = \mathcal{G}_h^x \cup \mathcal{G}_h^t$, where $\mathcal{G}_h^x$ collects all the space faces and $\mathcal{G}_h^t$ collects all the time faces of the space-time mesh. Thus, $\partial^x q \subset \mathcal{G}_h^x$ and $\partial^t q \subset \mathcal{G}_h^t$ for all $q \in Q_h$. Moreover, the set $\mathcal{G}_h^x$ is further partitioned as $\mathcal{G}_h^x = \mathcal{G}_h^{x,\text{int}} \cup \mathcal{G}_h^{x,\text{ext}}$, where $\mathcal{G}_h^{x,\text{int}}$ collects all the space faces that are subsets of $J \times \Omega$ (called interior space faces), whereas $\mathcal{G}_h^{x,\text{ext}}$ collects all the space faces that are subsets of $J \times \partial\Omega$ (called exterior space faces). Similarly, we consider the partition $\mathcal{G}_h^t = \mathcal{G}_h^{t,\text{int}} \cup \mathcal{G}_h^{t,\text{ext}}$ and observe that the time faces in $\mathcal{G}_h^{t,\text{ext}}$ correspond either to the initial or to the final time.

Finally, we set $\mathcal{G}_h^{\text{int}} := \mathcal{G}_h^{t,\text{int}} \cup \mathcal{G}_h^{x,\text{int}}$. For all $g \in \mathcal{G}_h^{\text{int}}$, we define $n_g := n^t$ if $g \in \mathcal{G}_h^{t,\text{int}}$ and $n_g := n_{F^x}$ if $g := T^t \times F^x \in \mathcal{G}_h^{x,\text{int}}$, and $[v]_g$ denotes the jump of the piecewise smooth, space-time function $v$ across $g$ in the direction of $n_g$. The notation is summarized in Figure 1.

<figure><img src="image_1.png" /><figcaption>Figure 1: Mesh notation</figcaption></figure>

## 2.2 Time-space polynomial spaces

Let $k \ge 1$ be the polynomial degree of the hybridized dG method in space and let $\ell \ge 1$ be the polynomial degree of the hybridized dG method in time. We denote by $\mathbb{P}^k(S)$ the set of polynomials of total degree at most $k$ on the subset $S \subseteq \Omega$. Moreover, for a linear space $U$ composed of functions defined on $S$, we denote by $\mathbb{P}^\ell(I; U)$ the set of $U$-valued polynomials of degree at most $\ell$ on $I \subseteq \bar{J} = [0, T_f]$.
The discrete unknowns are tensor-product space-time polynomials attached to every space-time cell, $q \in Q_h$, and to every space and time face, $g^x \in \mathcal{G}_h^x$ and $g^t \in \mathcal{G}_h^t$, respectively. For all $q := T^t \times T^x \in Q_h$, we set

$$U_q := \mathbb{P}^\ell(T^t; \mathbb{P}^k(T^x)). \qquad (5)$$

For all $g^x := T^t \times F^x \in \mathcal{G}_h^x$ and all $g^t := F^t \times T^x \in \mathcal{G}_h^t$, we set

$$U_{g^x} := \mathbb{P}^\ell(T^t; \mathbb{P}^k(F^x)), \quad U_{g^t} := \mathbb{P}^k(T^x). \qquad (6)$$

(There is no time variation in $U_{g^t}$ since $F^t$ is just a discrete time node). We then set

$$U_{\partial q} := U_{\partial^x q} \times U_{\partial^t q}, \quad U_{\partial^x q} := \prod_{g^x \in \partial^x q} U_{g^x}, \quad U_{\partial^t q} := \prod_{g^t \in \partial^t q} U_{g^t}. \qquad (7)$$

The discrete problem involves cell polynomials $u_Q := (u_q)_{q \in Q_h}$, and face polynomials $u_G := (u_g)_{g \in \mathcal{G}_h} := ((u_{g^t})_{g^t \in \mathcal{G}_h^t}, (u_{g^x})_{g^x \in \mathcal{G}_h^x})$. Notice that we partitioned the face polynomials into those attached to the time faces and those attached to the space faces. We introduce the space-time discrete space

$$\hat{U}_Q := U_Q \times U_G, \quad U_G := U_{G^t} \times U_{G^x}, \qquad (8)$$

with

$$U_Q := \{v_Q \in L^2(Q) \mid v_q := v_Q|_q \in U_q, \forall q \in Q_h\}, \qquad (9a)$$

$$U_{G^t} := \{v_{G^t} \in L^2(\mathcal{G}_h^t) \mid v_{g^t} := v_{G^t}|_{g^t} \in U_{G^t}, \forall g^t \in \mathcal{G}_h^t\}, \qquad (9b)$$

$$U_{G^x} := \{v_{G^x} \in L^2(\mathcal{G}_h^x) \mid v_{g^x} := v_{G^x}|_{g^x} \in U_{G^x}, \forall g^x \in \mathcal{G}_h^{x, \text{int}}, v_{g^x} = 0, \forall g^x \in \mathcal{G}_h^{x, \text{ext}}\}. \qquad (9c)$$

Notice that functions in $U_{G^x}$ vanish identically on $J \times \partial \Omega$. For a generic function $\hat{v}_Q \in \hat{U}_Q$, we write $\hat{v}_Q := (v_Q, v_G)$ with $v_Q \in U_Q$ and $v_G \in U_G$. Finally, we denote by

$$\hat{U}_{Q,0} := U_Q \times U_{G,0}, \quad U_{G,0} := U_{G^x} \times U_{G^t,0}, \qquad (10a)$$

$$U_{G^t,0} := \{v_{G^t} \in U_{G^t} \mid v_{g^t} = 0, \forall g^t \in \mathcal{G}_h^{t, \text{ext}}\}, \qquad (10b)$$

the subspace of $\hat{U}_Q$ composed of pairs whose time-face polynomials are null at the initial and at the final times. We will use the space $\hat{U}_Q$ for the discrete primal unknowns and the space $\hat{U}_{Q,0}$ for the discrete dual unknowns.

For a space-time cell polynomial $v_Q \in U_Q$, we define the broken gradient and Laplacian operators, $\nabla_Q$ and $\Delta_Q$, such that $(\nabla_Q v_Q)|_q := \nabla(v_Q|_q)$ and $(\Delta_Q v_Q)|_q := \Delta(v_Q|_q)$ for all $q \in Q_h$. A similar notation is adopted for the first-order and second-order time-derivatives, $\partial_{t,Q}$ and $\partial_{tt,Q}^2$. Finally, we set

$$L_Q(v_Q) := \partial_{tt,Q}^2 v_Q - \Delta_Q v_Q, \quad \forall v_Q \in U_Q. \qquad (11)$$

**Remark 2.1. (Cell polynomials)** The space-time polynomial space defined in (5) has a space-time tensor-product form. Indeed, any polynomial $v_q \in U_q$ can be written as

$$v_q = \sum_{\substack{\beta \in \{0, \ell\} \\ \alpha := (\alpha_1, \dots, \alpha_d) \in \mathbb{N}^d \\ 0 \le \alpha_1 + \dots + \alpha_d \le k}} \gamma_{\alpha\beta} t^\beta x^\alpha,$$

with $x^\alpha := x_1^{\alpha_1} \dots x_d^{\alpha_d}$.
## 2.3 Bilinear forms

We can now introduce the bilinear forms needed to formulate the discrete problem. Let $\hat{v}_Q$, $\hat{w}_Q$ be generic functions in $\hat{U}_Q$ (primal variables) and let $\hat{\zeta}_Q$, $\hat{\eta}_Q$ be generic functions in $\hat{U}_{Q,0}$ (dual variables). The bilinear form associated with the discretization of the wave equation in a generic space-time cell $q \in Q_h$ is

$$ a_q(\hat{v}_q, \hat{\eta}_q) := (\nabla v_q, \nabla \eta_q)_q - (\nabla v_q \cdot \mathbf{n} \partial^{x_q}, \eta_q - \eta \partial^{x_q}) \partial^{x_q} - (v_q - v \partial^{x_q}, \nabla \eta_q \cdot \mathbf{n} \partial^{x_q}) \partial^{x_q} \\ - (\partial_t v_q, \partial_t \eta_q)_q + (\partial_t v_q \mathbf{n} \partial^{t_q}, \eta_q - \eta \partial^{t_q}) \partial^{t_q} + (v_q - v \partial^{t_q}, \partial_t \eta_q \mathbf{n} \partial^{t_q}) \partial^{t_q}. \quad (12) $$

In order to allow for a more compact notation, we define

$$ B := \begin{pmatrix} -1 & 0_{1,d} \\ 0_{d,1} & \mathbb{I}_{d,d} \end{pmatrix} \in \mathbb{R}^{d+1,d+1}, \quad (13) $$

together with the space-time gradient operator $\tilde{\nabla} v := (\partial_t v, \nabla v)^T$ and the space-time unit outward normal vector $\tilde{n}_q \in \mathbb{R}^{d+1}$ for all $q \in Q_h$ such that $\tilde{n}_q|_{\partial^t q} := (\mathbf{n}_{\partial^t q}, \mathbf{0})^T \in \mathbb{R}^{d+1}$ and $\tilde{n}_q|_{\partial^{x_q}} := (\mathbf{0}, \mathbf{n}_{\partial^{x_q}})^T \in \mathbb{R}^{d+1}$. We also use a similar notation $\tilde{n}_g \in \mathbb{R}^{d+1}$ for all $g \in \mathcal{G}_h^{\text{int}}$ and we consider the broken version of the above space-time gradient operator, $\tilde{\nabla}_Q$. We can then rewrite the local bilinear form $a_q$ as follows:

$$ a_q(\hat{v}_q, \hat{\eta}_q) = (B\tilde{\nabla}v_q, \tilde{\nabla}\eta_q)_q - (B\tilde{\nabla}v_q \cdot \tilde{n}_q, \eta_q - \eta\partial_q)\partial_q - (v_q - v\partial_q, B\tilde{\nabla}\eta_q \cdot \tilde{n}_q)\partial_q. $$

The global bilinear form is assembled by summing over the space-time cells:

$$ a_Q(\hat{v}_Q, \hat{\eta}_Q) := \sum_{q \in Q_h} a_q(\hat{v}_q, \hat{\eta}_q). \quad (14) $$

Moreover, the local stabilization bilinear forms read as follows:

$$ s_q(\hat{v}_q, \hat{w}_q) := h^{-1}(v_q - v_{\partial q}, w_q - w_{\partial q})\partial_q, \quad (15a) $$

$$ \sigma_q(\hat{\zeta}_q, \hat{\eta}_q) := (\tilde{\nabla}\zeta_q, \tilde{\nabla}\eta_q)_q + s_q(\hat{\zeta}_q, \hat{\eta}_q), \quad (15b) $$

leading to the following global stabilization bilinear forms:

$$ s_Q(\hat{v}_Q, \hat{w}_Q) := \sum_{q \in Q_h} s_q(\hat{v}_q, \hat{w}_q), \quad (16a) $$

$$ \sigma_Q(\hat{\zeta}_Q, \hat{\eta}_Q) := \sum_{q \in Q_h} \sigma_q(\hat{\zeta}_q, \hat{\eta}_q) = (\tilde{\nabla}_Q \zeta_Q, \tilde{\nabla}_Q \eta_Q)_Q + s_Q(\hat{\zeta}_Q, \hat{\eta}_Q). \quad (16b) $$

Notice that $s_Q$ consists of least-squares penalties on the difference between the trace of the cell polynomials and the polynomials on the space and time faces of every space-time cell $q \in Q_h$. For the dual variables, we enforce the same penalty together with a least-squares penalty on the broken space-time derivatives.
## 2.4 Lagrangian and discrete problem

We want to find the saddle-point of the Lagrangian defined for all $(\hat{v}_Q, \hat{\zeta}_Q) \in \hat{U}_Q \times \hat{U}_{Q,0}$ by

$$ \mathcal{L}_Q(\hat{v}_Q, \hat{\zeta}_Q) := \frac{1}{2} \|v_Q - g_\delta\|_{J \times \bar{\omega}}^2 + \frac{1}{2} s_Q(\hat{v}_Q, \hat{v}_Q) - \frac{1}{2} \sigma_Q(\hat{\zeta}_Q, \hat{\zeta}_Q) \\ + a_Q(\hat{v}_Q, \hat{\zeta}_Q) - (f, \zeta_Q)_Q, \qquad (17) $$

where $g_\delta := g + \delta$ denotes the available perturbed measurement of $g$. Notice that there is no initial nor final condition on the primal variable, whereas the dual variable is null at initial and final time.

The discrete problem is derived by seeking a critical point of the Lagrangian and reads as follows: Find $(\hat{u}_Q, \hat{\xi}_Q) \in \hat{U}_Q \times \hat{U}_{Q,0}$ such that

$$ (u_Q, w_Q)_{J \times \bar{\omega}} + s_Q(\hat{u}_Q, \hat{w}_Q) + a_Q(\hat{w}_Q, \hat{\xi}_Q) = (g_\delta, w_Q)_{J \times \bar{\omega}}, \qquad (18a) $$

$$ a_Q(\hat{u}_Q, \hat{\eta}_Q) - \sigma_Q(\hat{\xi}_Q, \hat{\eta}_Q) = (f, \eta_Q)_Q, \qquad (18b) $$

where the first equation holds for all $\hat{w}_Q \in \hat{U}_Q$ and the second for all $\hat{\eta}_Q \in \hat{U}_{Q,0}$. For all $(\hat{v}_Q, \hat{\zeta}_Q)$ and $(\hat{w}_Q, \hat{\eta}_Q)$ in $\hat{U}_Q \times \hat{U}_{Q,0}$, we define the bilinear form

$$ A_Q((\hat{v}_Q, \hat{\zeta}_Q), (\hat{w}_Q, \hat{\eta}_Q)) := (v_Q, w_Q)_{J \times \bar{\omega}} + s_Q(\hat{v}_Q, \hat{w}_Q) + a_Q(\hat{w}_Q, \hat{\zeta}_Q) \\ + a_Q(\hat{v}_Q, \hat{\eta}_Q) - \sigma_Q(\hat{\zeta}_Q, \hat{\eta}_Q). \qquad (19) $$

The discrete problem (18) can be rewritten as follows: Find $(\hat{u}_Q, \hat{\xi}_Q) \in \hat{U}_Q \times \hat{U}_{Q,0}$ such that, for all $(\hat{w}_Q, \hat{\eta}_Q) \in \hat{U}_Q \times \hat{U}_{Q,0}$,

$$ A_Q((\hat{u}_Q, \hat{\xi}_Q), (\hat{w}_Q, \hat{\eta}_Q)) = (g_\delta, w_Q)_{J \times \bar{\omega}} + (f, \eta_Q)_Q. \qquad (20) $$

We will see in Corollary 5 that the discrete problem (18) (or equivalently, (20)) is well-posed.

The discrete problem (18) can be solved efficiently by eliminating locally all the cell degrees of freedom using a static condensation procedure. The global problem to be solved then involves only the face degrees of freedom. The stencil of the global problem couples unknowns attached to neighboring time and space faces (in the sense that two neighbors share a common space-time cell). To eliminate the cell degrees of freedom, one needs to solve the following linear system: Find $(u_Q, \xi_Q) \in U_Q \times U_Q$ such that, for all $(w_Q, \eta_Q) \in U_Q \times U_Q$,

$$ A_Q(((u_Q, 0), (\xi_Q, 0)), ((w_Q, 0), (\eta_Q, 0))) = (g_\delta, w_Q)_{J \times \bar{\omega}} + (f, \eta_Q)_Q \\ - A_Q(((0, u_G), (0, \xi_G)), ((w_Q, 0), (\eta_Q, 0))), \qquad (21) $$

where $u_G \in U_G$ and $\xi_G \in U_{G,0}$ are given. We will see in Corollary 5 that the discrete problem (21) is well-posed. At this stage, we only record that the linear system (21) is easy to solve since it amounts to a block linear system on the cell unknowns. Indeed, (21)
is equivalent to solving, for all $q := T^t \times T^x \in Q_h$, all $u_{\partial q} \in U_{\partial q}$, and all $\xi_{\partial q} \in U_{\partial q}$, the following local system: Find $(u_q, \xi_q) \in U_q \times U_q$ such that, for all $(w_q, \eta_q) \in U_q \times U_q$,

$$ (u_q, w_q)_{T^t \times (\varpi \cap T^x)} + s_q((u_q, 0), (w_q, 0)) + a_q((w_q, 0), (\xi_q, 0)) = \Phi_q^1(w_q), \quad (22a) $$

$$ a_q((u_q, 0), (\eta_q, 0)) - \sigma_q((\xi_q, 0), (\eta_q, 0)) = \Phi_q^2(\eta_q), \quad (22b) $$

with $\Phi_q^1(w_q) := (g_\delta, w_q)_{T^t \times (\varpi \cap T^x)} - s_q((0, u_{\partial q}), (w_q, 0)) - a_q((w_q, 0), (0, \xi_{\partial q}))$ and $\Phi_q^2(\eta_q) := (f, \eta_q)_q - a_q((0, u_{\partial q}), (\eta_q, 0)) + \sigma_q((0, \xi_{\partial q}), (\eta_q, 0))$.

# 3 Analysis

In this section, we first study the stability property of the discrete operator $A_Q$ in a suitable residual norm. Next, we introduce interpolation operators in space and in time, and we bound the consistency error. We combine these bounds with the abstract conditional estimate from Lemma 2 to derive error estimates in the shifted energy norm. Finally, we establish the convergence rates for the method.

In what follows, we use the convention $A \lesssim B$ to abbreviate the inequality $A \le CB$ for positive real numbers $A$ and $B$, where the constant $C > 0$ does not depend on $h$, $\tau$, and the solution of the continuous and discrete problems.

*Remark 3.1. (Domain regularity)* The convergence analysis is conducted under the twofold assumption that (i) $\Omega$ is a Lipschitz polyhedron so that it can be meshed exactly by an affine mesh whose cells have a piecewise constant outward unit normal; (ii) the stability estimate from Lemma 2 holds true for such a domain $\Omega$. The numerical tests are also performed in this setting with $\Omega$ taken as a square domain. We notice that treating a curved domain would be quite technical due to the hybrid nature of the discretization method (see for instance [14]), and the modifications to the implementation would be non-trivial as well. Establishing Lemma 2 on a polygonal domain is a key issue that we postpone to future work.

## 3.1 Residual stability

Recall the broken operator $L_Q$ defined in (11) and the notation $G_h^{\text{int}} := G_h^{\text{t,int}} \cup G_h^{\text{x,int}}$ together with the unit normal vector $\tilde{n}_q \in \mathbb{R}^{d+1}$. For all $(\hat{v}_Q, \hat{\zeta}_Q) \in \hat{U}_Q \times \hat{U}_{Q,0}$, we define the residual norm

$$ \|\hat{v}_Q, \hat{\zeta}_Q\|^2 := \|v_Q\|_{\mathbb{R}}^2 + \|v_Q\|_{J \times \varpi}^2 + s_Q(\hat{v}_Q, \hat{v}_Q) + \sigma_Q(\hat{\zeta}_Q, \hat{\zeta}_Q), \quad (23) $$

with

$$ \|v_Q\|_{\mathbb{R}}^2 := \|hL_Q(v_Q)\|_Q^2 + \|h^{\frac{1}{2}}[\tilde{\nabla}_Q v_Q]_{G^{\text{int}}} \cdot \tilde{n}_{G^{\text{int}}}\|_{G^{\text{int}}}^2, \quad (24) $$

$$ \text{and } \|h^{\frac{1}{2}}[\tilde{\nabla}_Q v_Q]_{G^{\text{int}}} \cdot \tilde{n}_{G^{\text{int}}}\|_{G^{\text{int}}}^2 := \sum_{g \in G_h^{\text{int}}} h \|[\tilde{\nabla}_Q v_Q]_g \cdot \tilde{n}_g\|_g^2. $$

**Lemma 3 (Norm).** $\|\cdot\|$ defines a norm on $\hat{U}_Q \times \hat{U}_{Q,0}$.
Proof. Let $(\hat{v}_Q, \hat{\zeta}_Q) \in \hat{U}_Q \times \hat{U}_{Q,0}$ be such that $\|\hat{v}_Q, \hat{\zeta}_Q\| = 0$.

(1) The dual component $\hat{\zeta}_Q := (\zeta_Q, \zeta_G) \in \hat{U}_{Q,0}$ satisfies $\sigma_Q(\hat{\zeta}_Q, \hat{\zeta}_Q) = 0$. Hence, $\zeta_Q$ is piecewise constant in $Q$ and $\zeta_G$ is equal to the trace of $\zeta_Q$ on all faces. Hence, $\zeta_Q$ is globally constant, and it vanishes since $\zeta_G$ vanishes on $G_h^{x,ext}$ and $G_h^{t,ext}$. In conclusion, $\hat{\zeta}_Q = (0,0)$.

(2) The primal component $\hat{v}_Q := (v_Q, v_G) \in \hat{U}_Q$ satisfies $\|v_Q\|_R = \|v_Q\|_{J \times \bar{\omega}} = s_Q(\hat{v}_Q, \hat{v}_Q) = 0$. Since $L_Q(v_Q) = 0$, $v_Q$ is a piecewise affine polynomial in space and in time. Moreover, $s_Q(\hat{v}_Q, \hat{v}_Q) = 0$ implies that $v_Q$ is globally continuous and $v_G$ is obtained from the traces of $v_Q$ on the faces of $G_h$. Combined with the fact that $\|h^{\frac{1}{2}}[\tilde{\nabla}_Q v_Q]\|_{G^{\text{int}}} \cdot \tilde{n}_{G^{\text{int}}}\|_{G^{\text{int}}} = 0$, we infer that $v_Q$ is globally affine in $Q$. Finally, invoking that $\bar{\omega}$ contains at least one space-time cell or that the face components vanish on $G_h^{x,ext}$, we conclude that $\hat{v}_Q = (0,0)$. $\square$

Our key stability result is the following inf-sup condition.

**Lemma 4 (Inf-sup condition).** The following holds for all $(\hat{v}_Q, \hat{\zeta}_Q) \in \hat{U}_Q \times \hat{U}_{Q,0}$,

$$ \|\hat{v}_Q, \hat{\zeta}_Q\| \lesssim \sup_{(\hat{w}_Q, \hat{\eta}_Q) \in \hat{U}_Q \times \hat{U}_{Q,0} \setminus \{(0,0)\}} \frac{A_Q((\hat{v}_Q, \hat{\zeta}_Q), (\hat{w}_Q, \hat{\eta}_Q))}{\|\hat{w}_Q, \hat{\eta}_Q\|} \quad (25) $$

Proof. Let us denote by $S$ the right-hand side of (25).

(1) Stabilizations. We first use the test functions $\hat{w}_Q := \hat{v}_Q$ and $\hat{\eta}_Q := -\hat{\zeta}_Q$ to get

$$ \|v_Q\|_{J \times \bar{\omega}}^2 + s_Q(\hat{v}_Q, \hat{v}_Q) + \sigma_Q(\hat{\zeta}_Q, \hat{\zeta}_Q) = A_Q((\hat{v}_Q, \hat{\zeta}_Q), (\hat{v}_Q, -\hat{\zeta}_Q)) \le S \|\hat{v}_Q, \hat{\zeta}_Q\|. $$

(2) Gradient jumps. Let $\hat{\eta}_Q := (0, \eta_G) \in \hat{U}_{Q,0}$ with $\eta_g := h[\|B\tilde{\nabla}_Q v_Q\|_g \cdot \tilde{n}_g$ for all $g \in G_h^{\text{int}}$ and $\eta_g := 0$ otherwise. Since $\|0, \hat{\eta}_Q\| = \sigma_Q(\hat{\eta}_Q, \hat{\eta}_Q)^{\frac{1}{2}}$, integration by parts in time and in space gives

$$ \begin{aligned} \|h^{\frac{1}{2}}[\tilde{\nabla}_Q v_Q]\|_{G^{\text{int}}} \cdot \tilde{n}_{G^{\text{int}}} \|_{G^{\text{int}}}^2 &= a_Q(\hat{v}_Q, \hat{\eta}_Q) \\ &= A_Q((\hat{v}_Q, \hat{\zeta}_Q), (0, \hat{\eta}_Q)) + \sigma_Q(\hat{\zeta}_Q, \hat{\eta}_Q) \\ &\le S \sigma_Q(\hat{\eta}_Q, \hat{\eta}_Q)^{\frac{1}{2}} + \sigma_Q(\hat{\zeta}_Q, \hat{\eta}_Q) \\ &\le (S + \sigma_Q(\hat{\zeta}_Q, \hat{\zeta}_Q)^{\frac{1}{2}}) \sigma_Q(\hat{\eta}_Q, \hat{\eta}_Q)^{\frac{1}{2}}. \end{aligned} $$

Moreover, we also have

$$ \sigma_Q(\hat{\eta}_Q, \hat{\eta}_Q) = \sum_{q \in Q_h} h^{-1} \|\eta_{\partial q}\|_{\partial q}^2 \lesssim \|h^{\frac{1}{2}}[\tilde{\nabla}_Q v_Q]\|_{G^{\text{int}}} \cdot \tilde{n}_{G^{\text{int}}} \|_{G^{\text{int}}}^2. $$

This implies that

$$ \|h^{\frac{1}{2}}[\tilde{\nabla}_Q v_Q]\|_{G^{\text{int}}} \cdot \tilde{n}_{G^{\text{int}}} \|_{G^{\text{int}}}^2 \lesssim S^2 + \sigma_Q(\hat{\zeta}_Q, \hat{\zeta}_Q). $$

(3) Cell residuals. We now consider $\hat{\eta}_Q := (\eta_Q, 0)$ with $\eta_Q := h^2 L_Q(v_Q)$. Integration by parts in time and in space gives

$$ \begin{aligned} a_Q(\hat{v}_Q, \hat{\eta}_Q) &= \sum_{q \in Q_h} \{ (B \tilde{\nabla} v_q, \tilde{\nabla} \eta_q)_q - (B \tilde{\nabla} v_q \cdot \tilde{n}_q, \eta_q)_{\partial q} - (v_q - v_{\partial q}, B \tilde{\nabla} \eta_q \cdot \tilde{n}_q)_{\partial q} \} \\ &= (L_Q(v_Q), \eta_Q)_Q - \sum_{q \in Q_h} (v_q - v_{\partial q}, B \tilde{\nabla} \eta_q \cdot \tilde{n}_q)_{\partial q}. \end{aligned} $$
Using the definition of $\hat{\eta}_Q$, we have

$$
\begin{align*}
\|hL_Q(v_Q)\|_Q^2 &= a_Q(\hat{v}_Q, \hat{\eta}_Q) + \sum_{q \in Q_h} (v_q - v_{\partial q}, B\tilde{\nabla}\eta_q \cdot \tilde{n}_q)_{\partial q} \\
&= A_Q((\hat{v}_Q, \hat{\zeta}_Q), (0, \hat{\eta}_Q)) + \sigma_Q(\hat{\zeta}_Q, \hat{\eta}_Q) + \sum_{q \in Q_h} (v_q - v_{\partial q}, B\tilde{\nabla}\eta_q \cdot \tilde{n}_q)_{\partial q} \\
&\lesssim S\|0, \hat{\eta}_Q\| + (\sigma_Q(\hat{\zeta}_Q, \hat{\zeta}_Q) + s_Q(\hat{v}_Q, \hat{v}_Q))^{\frac{1}{2}} \sigma_Q(\hat{\eta}_Q, \hat{\eta}_Q)^{\frac{1}{2}},
\end{align*}
$$

since a discrete trace inverse inequality implies that

$$
\sum_{q \in Q_h} h \| \tilde{\nabla} \eta_q \cdot \tilde{n}_q \|_{\partial q}^2 \lesssim \| \tilde{\nabla}_Q \eta_Q \|_Q^2 \le \sigma_Q(\hat{\eta}_Q, \hat{\eta}_Q).
$$

Furthermore, invoking again inverse inequalities, we have

$$
\sigma_Q(\hat{\eta}_Q, \hat{\eta}_Q) = \sum_{q \in Q_h} \left( \|\tilde{\nabla} \eta_q\|_q^2 + h^{-1} \|\eta_q\|_{\partial q}^2 \right) \lesssim \|h^{-1} \eta_Q\|_Q^2 = \|hL_Q(v_Q)\|_Q^2.
$$

Using $\|0, \hat{\eta}_Q\|^2 = \sigma_Q(\hat{\eta}_Q, \hat{\eta}_Q)$ and Young's inequality, this implies that

$$
\|hL_Q(v_Q)\|_Q^2 \lesssim S^2 + s_Q(\hat{v}_Q, \hat{v}_Q) + \sigma_Q(\hat{\zeta}_Q, \hat{\zeta}_Q).
$$

Gathering the previous estimates leads to the expected inf-sup condition.

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

**Corollary 5 (Well-posedness).** The discrete problems (18) and (21) are well-posed.

*Proof.* Since both (18) and (21) amount to square linear systems, it suffices to prove that, for each system, the only solution with zero right-hand side is zero. For (18), this follows from the inf-sup condition and Lemma 3. For (21), one needs to adapt the above proofs. First, adapting the proof of Lemma 3, one observes that the norm

$$
\|v_Q, \zeta_Q\|^2 := \|hL_Q(v_Q)\|_Q^2 + \|v_Q\|_{J \times \omega}^2 + s_Q((v_Q, 0), (v_Q, 0)) + \sigma_Q((\zeta_Q, 0), (\zeta_Q, 0)),
$$

defines a norm on $U_Q \times U_Q$. Then, adapting the proof of Lemma 4, one observes that an inf-sup condition holds true for the above norm when considering the bilinear form $A_Q(((u_Q, 0), (\xi_Q, 0)), ((w_Q, 0), (\eta_Q, 0)))$. $\square$

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

3.2 Interpolation operator

In this section, we define the space-time interpolation operator used in the error analysis.
Let $q \in Q_h$. For all $v \in L^2(q)$, we define $I_q(v)$ to be the $L^2$-orthogonal projection of
$v$ onto $U_q$. For all $v \in H^1(q)$, we define $I_{\partial q}(v)$ as the the $L^2$-orthogonal projection of
$v|_{\partial q}$ onto $U_{\partial q}$. We also define $\hat{I}_q(v) := (I_q(v), I_{\partial q}(v))$. Similarly, for all $v \in H^1(Q)$, we
define the interpolation on the whole discrete space by $\hat{I}_Q(v) := (I_Q(v), I_G(v))$, where
$I_Q(v)|_q := I_q(v)$ for all $q \in Q_h$ and $I_G(v)|_g := I_{\partial q}(v)|_g$ for all $g \in G_h$ with $g \subset \partial q$. Notice
that $\hat{I}_Q(v) \in \hat{U}_Q$ if $v=0$ on $J \times \partial \Omega$, and $\hat{I}_Q(v) \in \hat{U}_{Q,0}$ if additionally $v(0) = v(T_f) = 0$.

The above interpolation operator enjoys the following convergence properties (the proof is omitted since it is standard).
**Lemma 6 (Approximation).** Let $k \ge 1$ and $\ell \ge 1$. The following holds for all $q := T^t \times T^x \in Q_h$, $T^t \in \mathcal{T}_\tau^t$, $T^x \in \mathcal{T}_h^x$, all $v \in H^{\ell+1}(T^t; H^2(T^x)) \cap H^2(T^t; H^{k+1}(T^x))$,

$$ \|I_q(v) - v\|_q \lesssim h^{\ell+1} \|v\|_{H^{\ell+1}(T^t; L^2(T^x))} + h^{k+1} \|v\|_{L^2(T^t; H^{k+1}(T^x))}, \quad (26a) $$

$$ \|\tilde{\nabla}(I_q(v) - v)\|_q \lesssim h^\ell \|v\|_{H^{\ell+1}(T^t; H^1(T^x))} + h^k \|v\|_{H^1(T^t; H^{k+1}(T^x))}, \quad (26b) $$

$$ h^{\frac{1}{2}} \|\tilde{\nabla}(I_q(v) - v)\|_{\partial q} \lesssim h^\ell \|v\|_{H^{\ell+1}(T^t; H^2(T^x))} + h^k \|v\|_{H^2(T^t; H^{k+1}(T^x))}, \quad (26c) $$

$$ h \|L(I_q(v) - v)\|_q \lesssim h^\ell \|v\|_{H^{\ell+1}(T^t; H^2(T^x))} + h^k \|v\|_{H^2(T^t; H^{k+1}(T^x))}. \quad (26d) $$

Moreover, we have

$$ \begin{aligned} h^{-\frac{1}{2}} \|I_q(v) - v\|_{\partial q} &+ h^{-\frac{1}{2}} \|I_{\partial q}(v) - v\|_{\partial q} \\ &\lesssim h^\ell \|v\|_{H^{\ell+1}(T^t; H^1(T^x))} + h^k \|v\|_{H^1(T^t; H^{k+1}(T^x))}. \end{aligned} \quad (27) $$

## 3.3 Error decomposition, consistency and a priori residual bound

Recall that $u$ denotes the solution to the exact problem (1) and that $(\hat{u}_Q, \hat{\xi}_Q)$ denotes the solution to the discrete problem (18). The discrete and interpolation errors on the primal unknown are defined as

$$ \hat{e}_Q := \hat{u}_Q - \hat{I}_Q(u), \quad \hat{\theta}_Q := (\theta_Q, \theta_G), \quad (28) $$

where $\theta_q := u|_q - I_q(u)$ for all $q \in Q_h$ and $\theta_g := u|_g - I_G(u)|_g$ for all $g \in G_h$.

We now bound the consistency error. To this purpose, we consider the norm

$$ \|\hat{\theta}_Q\|_\# := \|\tilde{\nabla}_Q \theta_Q\|_Q^2 + \|\theta_Q\|_Q^2 \\ + \sum_{q \in Q_h} \{ \|hL(\theta_q)\|_q^2 + \|h^{\frac{1}{2}} \tilde{\nabla} \theta_q\|_{\partial q}^2 + \|h^{-\frac{1}{2}} \theta_q\|_{\partial q}^2 + \|h^{-\frac{1}{2}} \theta_{\partial q}\|_{\partial q}^2 \}. \quad (29) $$

Recall that the stability norm $\|\cdot\|$ is defined in (23)-(24). In what follows, we assume that

$$ u \in H^2(J; L^2(\Omega)) \cap L^2(J; H^2(\Omega)). \quad (30) $$

**Lemma 7 (Consistency).** Let $(\hat{u}_Q, \hat{\xi}_Q)$ denote the solution to the discrete problem (18). Let $\hat{e}_Q$ and $\hat{\theta}_Q$ be the discrete and interpolation errors on the primal unknown defined in (28). Under the regularity assumption (30), we have, for all $(\hat{w}_Q, \hat{\eta}_Q) \in \hat{U}_Q \times \hat{U}_{Q,0}$,

$$ |A_Q((\hat{e}_Q, \hat{\xi}_Q), (\hat{w}_Q, 0))| \lesssim (\|\hat{\theta}_Q\|_\# + \|\delta\|_{J \times \bar{\omega}}) \|\hat{w}_Q, 0\|, \quad (31a) $$

$$ |A_Q((\hat{e}_Q, \hat{\xi}_Q), (0, \hat{\eta}_Q))| \lesssim \|\hat{\theta}_Q\|_\# \|0, \hat{\eta}_Q\|. \quad (31b) $$

*Proof.* Proof of (31a). Since $(\hat{e}_Q, \hat{\xi}_Q) = (\hat{u}_Q, \hat{\xi}_Q) - (\hat{I}_Q(u), 0)$, the equation (18a) in the discrete problem, the definition (19) of $A_Q$, and $g_\delta = g + \delta$ imply that, for all $\hat{w}_Q \in \hat{U}_Q$,

$$ \begin{aligned} A_Q((\hat{e}_Q, \hat{\xi}_Q), (\hat{w}_Q, 0)) &= (g_\delta, w_Q)_{J \times \bar{\omega}} - (I_Q(u), w_Q)_{J \times \bar{\omega}} - s_Q(\hat{I}_Q(u), \hat{w}_Q) \\ &= (\delta, w_Q)_{J \times \bar{\omega}} + (\theta_Q, w_Q)_{J \times \bar{\omega}} - s_Q(\hat{I}_Q(u), \hat{w}_Q). \end{aligned} $$
The Cauchy–Schwarz inequality yields

$$|(\delta, w_Q)_{J \times \bar{\varpi}}| + |(\theta_Q, w_Q)_{J \times \bar{\varpi}}| \lesssim (\|\delta\|_{J \times \bar{\varpi}} + \|\theta_Q\|_{J \times \bar{\varpi}}) \|w_Q\|_{J \times \bar{\varpi}}, \\ |s_Q(\hat{I}_Q(u), \hat{w}_Q)| \lesssim s_Q(\hat{I}_Q(u), \hat{I}_Q(u))^{\frac{1}{2}} s_Q(\hat{w}_Q, \hat{w}_Q)^{\frac{1}{2}}.$$

Moreover, we have

$$s_Q(\hat{I}_Q(u), \hat{I}_Q(u)) = \sum_{q \in Q_h} h^{-1} \|I_q(u) - I_{\partial q}(u)\|_{\partial q}^2 \lesssim \|\hat{\theta}_Q\|_{\#}^2.$$

The bound (31a) follows by gathering the above estimates.

Proof of (31b). Proceeding as above and using now the equation (18b) and $f = L(u)$, we have, for all $\hat{\eta}_Q \in \hat{U}_{Q,0}$,

$$A_Q((\hat{e}_Q, \hat{\xi}_Q), (0, \hat{\eta}_Q)) = (L(u), \eta_Q)_Q - a_Q(\hat{I}_Q(u), \hat{\eta}_Q).$$

Integrating by parts the first term on the right-hand side gives

$$(L(u), \eta_Q)_Q = \sum_{q \in Q_h} \{(B \tilde{\nabla} u, \tilde{\nabla} \eta_q)_q - (B \tilde{\nabla} u \cdot \tilde{n}_q, \eta_q - \eta_{\partial q})_{\partial q}\},$$

where we used the regularity assumption on $u$ and the fact $\eta_{\partial q}$ vanishes on all the boundary faces in space and in time. Recalling the definition (12)-(14) of $a_Q$, we infer that

$$A_Q((\hat{e}_Q, \hat{\xi}_Q), (0, \hat{\eta}_Q)) = \sum_{q \in Q_h} \{(B \tilde{\nabla} \theta_q, \tilde{\nabla} \eta_q)_q - (B \tilde{\nabla} \theta_q \cdot \tilde{n}_q, \eta_q - \eta_{\partial q})_{\partial q} \\ - (\theta_q - \theta_{\partial q}, B \tilde{\nabla} \eta_q \cdot \tilde{n}_q)_{\partial q}\}.$$

Using Cauchy–Schwarz and inverse inequalities, we readily obtain (31b). □

**Lemma 8 (A priori residual bound).** Under the regularity assumption (30), we have

$$\|\hat{e}_Q, \hat{\xi}_Q\| + s_Q(\hat{u}_Q, \hat{w}_Q)^{\frac{1}{2}} \lesssim \|\hat{\theta}_Q\|_{\#} + \|\delta\|_{J \times \bar{\varpi}}. \qquad (32)$$

*Proof.* Using inf-sup stability (Lemma 4), we have

$$\|\hat{e}_Q, \hat{\xi}_Q\| \lesssim \sup_{(\hat{w}_Q, \hat{\eta}_Q) \in \hat{U}_Q \times \hat{U}_{Q,0} \setminus \{(0,0)\}} \frac{A_Q((\hat{e}_Q, \hat{\xi}_Q), (\hat{w}_Q, \hat{\eta}_Q))}{\|\hat{w}_Q, \hat{\eta}_Q\|}.$$

Moreover, owing to consistency (Lemma 7), we have, for all $(\hat{w}_Q, \hat{\eta}_Q) \in \hat{U}_Q \times \hat{U}_{Q,0}$,

$$|A_Q((\hat{e}_Q, \hat{\xi}_Q), (\hat{w}_Q, \hat{\eta}_Q))| \lesssim (\|\hat{\theta}_Q\|_{\#} + \|\delta\|_{J \times \bar{\varpi}}) \|\hat{w}_Q, \hat{\eta}_Q\|.$$

Combining the two bounds gives

$$\|\hat{e}_Q, \hat{\xi}_Q\| \lesssim \|\hat{\theta}_Q\|_{\#} + \|\delta\|_{J \times \bar{\varpi}}.$$

Finally, observing that

$$s_Q(\hat{u}_Q, \hat{w}_Q)^{\frac{1}{2}} \le s_Q(\hat{I}_Q(u), \hat{I}_Q(u))^{\frac{1}{2}} + s_Q(\hat{e}_Q, \hat{e}_Q)^{\frac{1}{2}} \le \|\hat{\theta}_Q\|_{\#} + \|\hat{e}_Q, \hat{\xi}_Q\|,$$

concludes the proof. □
### 3.4 Bound on dual error norm

In this section, we prove another important result bounding some dual norm of the error. Recall the operator $L := \partial_{tt}^2 - \Delta$ so that, for all $v \in H^1(Q)$ and all $\eta \in H_0^1(Q)$, we have $\langle L(v), \eta \rangle_{H^{-1}, H_0^1} = (B\tilde{\nabla}v, \tilde{\nabla}\eta)_Q$, where $\langle \cdot, \cdot \rangle_{H^{-1}, H_0^1}$ stands for the duality product between $H^{-1}(Q)$ and $H_0^1(Q)$. Observing that the broken operator $L_Q$ acting on discrete functions $v_Q \in U_Q$ satisfies $\langle L_Q(v_Q), \eta \rangle_{H^{-1}, H_0^1} = (B\tilde{\nabla}_Q v_Q, \tilde{\nabla}\eta)_Q$, we can consider the extended operator $L_Q: H^1(Q) + U_Q \to H^{-1}(Q)$ (we use the same notation for simplicity). For all $v \in H^1(Q) + U_Q$, the corresponding dual norm is

$$ \|L_Q(v)\|_{H^{-1}(Q)} := \sup_{\substack{\eta \in H_0^1(Q) \\ \|\tilde{\nabla}\eta\|_Q = 1}} (B\tilde{\nabla}_Q v, \tilde{\nabla}\eta)_Q. $$

**Lemma 9 (Bound on dual error norm).** Under the regularity assumption (30), we have

$$ \|L_Q(u_Q - u)\|_{H^{-1}(Q)} \lesssim \|\hat{\theta}_Q\|_{\#} + \|\delta\|_{J \times \bar{\omega}}. $$

*Proof.* Let $\eta \in H_0^1(Q)$ with $\|\tilde{\nabla}\eta\|_Q = 1$. We have

$$ \langle L_Q(u_Q - u), \eta \rangle_{H^{-1}, H_0^1} = (B\tilde{\nabla}_Q u_Q, \tilde{\nabla}\eta)_Q - (L(u), \eta)_Q. $$

Let $\hat{\eta}_Q \in \hat{U}_{Q,0}$ be such that $\hat{\eta}_Q := \hat{I}_Q(\eta)$. Invoking equation (18b) gives

$$ a_Q(\hat{u}_Q, \hat{\eta}_Q) - \sigma_Q(\hat{\xi}_Q, \hat{\eta}_Q) = (L(u), \eta)_Q. $$

Subtracting this equation from the above expression, and re-arranging the terms gives

$$ \langle L_Q(u_Q - u), \eta \rangle_{H^{-1}, H_0^1} = A_1 + A_2 + A_3, $$

with

$$ A_1 := (L_Q(u_Q - u), \eta - \eta_Q)_Q, $$

$$ A_2 := (B\tilde{\nabla}_Q u_Q, \tilde{\nabla}\eta)_Q - a_Q(\hat{u}_Q, \hat{\eta}_Q) - (L_Q(u_Q), \eta - \eta_Q)_Q, $$

$$ A_3 := \sigma_Q(\hat{\xi}_Q, \hat{\eta}_Q). $$

We have

$$ \begin{aligned} & \sum_{q \in Q_h} \{(B\tilde{\nabla}u_q, \tilde{\nabla}\eta)_q - a_q(\hat{u}_q, \hat{\eta}_q)\} \\ &= \sum_{q \in Q_h} \{(B\tilde{\nabla}u_q, \tilde{\nabla}(\eta - \eta_q))_q + (B\tilde{\nabla}u_q \cdot \tilde{n}_q, \eta_q - \eta_{\partial q})_{\partial q} + (u_q - u_{\partial q}, B\tilde{\nabla}\eta_q \cdot \tilde{n}_q)_{\partial q}\} \\ &= \sum_{q \in Q_h} \{(L(u_q), \eta - \eta_q)_q + (B\tilde{\nabla}u_q \cdot \tilde{n}_q, \eta - \eta_{\partial q})_{\partial q} + (u_q - u_{\partial q}, B\tilde{\nabla}\eta_q \cdot \tilde{n}_q)_{\partial q}\} \\ &= (L_Q(u_Q), \eta - \eta_Q)_Q + \sum_{q \in Q_h} (u_q - u_{\partial q}, B\tilde{\nabla}\eta_q \cdot \tilde{n}_q)_{\partial q}, \end{aligned} $$
where we used the definition (12) of the bilinear form $a_q$, the fact that $(B\tilde{\nabla}u_q \cdot \tilde{n}_q)|_{\partial q} \in U_{\partial q}$ and the definition of $\eta_{\partial q}$ so that $(B\tilde{\nabla}u_q \cdot \tilde{n}_q, \eta - \eta_{\partial q})_{\partial q} = 0$. We can then rewrite $A_2$ as

$$A_2 = \sum_{q \in Q_h} (u_q - u_{\partial q}, B\tilde{\nabla}\eta_q \cdot \tilde{n}_q)_{\partial q}.$$

It still remains to bound $A_1$, $A_2$, $A_3$.

Bound on $A_1$. Since $\theta_Q = u - I_Q(u)$, we have

$$L_Q(u_Q - u) = L_Q(u_Q - I_Q(u)) - L_Q(\theta_Q).$$

Recalling the definition (24) of the $\|\cdot\|_R$-norm, we have

$$|(L_Q(u_Q - I_Q(u)), \eta - \eta_Q)_Q| \lesssim h^{-1} \|e_Q\|_R h \|\tilde{\nabla}\eta\|_Q,$$

where we used that $\|\eta - \eta_Q\|_Q \lesssim h\|\tilde{\nabla}\eta\|_Q$. Owing to Lemma 8, we have $\|e_Q\|_R \lesssim \|\hat{\theta}_Q\|_{\#} + \|\delta\|_{J \times \bar{\omega}}$. Moreover, the above bound on $\|\eta - \eta_Q\|_Q$ also gives

$$|(L_Q(\theta_Q), \eta - \eta_Q)_Q| \lesssim h \|L_Q(\theta_Q)\|_Q \|\tilde{\nabla}\eta\|_Q \le \|\hat{\theta}_Q\|_{\#} \|\tilde{\nabla}\eta\|_Q.$$

Gathering the above estimates yields

$$|A_1| \lesssim (\|\hat{\theta}_Q\|_{\#} + \|\delta\|_{J \times \bar{\omega}}) \|\tilde{\nabla}\eta\|_Q.$$

Bound on $A_2$. Using Cauchy-Schwarz and inverse trace inequalities and interpolation stability, we have

$$|A_2| = \left| \sum_{q \in Q_h} (u_q - u_{\partial q}, B\tilde{\nabla}\eta_q \cdot \tilde{n}_q)_{\partial q} \right| \lesssim s_Q(\hat{u}_Q, \hat{u}_Q)^{\frac{1}{2}} \|\tilde{\nabla}_Q \eta_Q\|_Q \\ \lesssim (\|\hat{\theta}_Q\|_{\#} + \|\delta\|_{J \times \bar{\omega}}) \|\tilde{\nabla}\eta\|_Q,$$

where we used Lemma 8 and $\|\tilde{\nabla}_Q \eta_Q\|_Q \lesssim \|\tilde{\nabla}\eta\|_Q$ in the last bound.

Bound on $A_3$. Recalling that $\sigma_Q(\hat{\eta}_Q, \hat{\eta}_Q) = \|\tilde{\nabla}_Q \eta_Q\|_Q^2 + s_Q(\hat{\eta}_Q, \hat{\eta}_Q)$, we infer that $\sigma_Q(\hat{\eta}_Q, \hat{\eta}_Q)^{\frac{1}{2}} \lesssim \|\tilde{\nabla}_Q \eta_Q\|_Q \lesssim \|\tilde{\nabla}\eta\|_Q$. This implies that

$$|\sigma_Q(\hat{\xi}_Q, \hat{\eta}_Q)| \lesssim \sigma_Q(\hat{\xi}_Q, \hat{\xi}_Q)^{\frac{1}{2}} \sigma_Q(\hat{\eta}_Q, \hat{\eta}_Q)^{\frac{1}{2}} \lesssim \sigma_Q(\hat{\xi}_Q, \hat{\xi}_Q)^{\frac{1}{2}} \|\tilde{\nabla}\eta\|_Q,$$

and $\sigma_Q(\hat{\xi}_Q, \hat{\xi}_Q)^{\frac{1}{2}}$ is bounded in Lemma 8. Hence, we have

$$|A_3| \lesssim (\|\hat{\theta}_Q\|_{\#} + \|\delta\|_{J \times \bar{\omega}}) \|\tilde{\nabla}\eta\|_Q.$$

Combining the above estimates and recalling that $\|\tilde{\nabla}\eta\|_Q = 1$ gives the expected bound. $\square$
### 3.5 Main result: error estimate in shifted energy norm

We are now ready to derive our main error estimate. The idea of the proof is to combine the results of Sections 3.3 and 3.4 with the conditional stability estimate from Lemma 2. We use the shifted energy norm which we extend to $H^1(Q) + U_Q$ by setting

$$ \|v\|_{\text{sft}} := \|v\|_{L^{\infty}(J;L^2(\Omega))} + \|\partial_{t,Q}v\|_{L^2(J;H^{-1}(\Omega))} + \|\nabla_Q v\|_{H^{-1}(J;L^2(\Omega))}. \quad (33) $$

We shall invoke the following classical result on finite element interpolation, here written in a space-time setting using the space-time mesh regularity assumption (i.e. $h \lesssim \tau \lesssim h$). It is established first on pairs of the form $(v_Q, 0)$ (see for instance [27, 6, 24, 26]) and then extended to pairs $\hat{v}_Q \in \hat{U}_Q$ by invoking a triangle inequality.

**Lemma 10.** For all $\hat{v}_Q \in \hat{U}_Q$, there exists $\tilde{v}_Q \in U_Q \cap H^1(Q)$ such that

$$ \|h^{-1}(\tilde{v}_Q - v_Q)\|_Q + \|\tilde{\nabla}_Q(\tilde{v}_Q - v_Q)\|_Q \lesssim s_Q(\hat{v}_Q, \hat{v}_Q)^{\frac{1}{2}}. $$

Recall the constant $C_{\text{stb}} > 0$ introduced in Lemma 2. Recall that the $\|\cdot\|_{\#}$-norm is defined in (29). We have the following convergence result.

**Theorem 11** (Error estimate in shifted energy norm). Under the regularity assumption (30) and Assumption 1, we have

$$ \|u - u_Q\|_{\text{sft}} \lesssim (1 + C_{\text{stb}})(\|\hat{\theta}_Q\|_{\#} + \|\delta\|_{J \times \bar{\omega}}). \quad (34) $$

*Proof.* (1) $H^1$-conforming approximation. Let $\tilde{u}_Q$ be the $H^1$-conforming approximation to $\hat{u}_Q$ given by Lemma 10. We have

$$ \|h^{-1}(\tilde{u}_Q - u_Q)\|_Q + \|\tilde{\nabla}_Q(\tilde{u}_Q - u_Q)\|_Q \lesssim s_Q(\hat{u}_Q, \hat{u}_Q)^{\frac{1}{2}}. \quad (35) $$

Invoking the triangle inequality yields

$$ \|u - u_Q\|_{\text{sft}} \le \|u_Q - \tilde{u}_Q\|_{\text{sft}} + \|u - \tilde{u}_Q\|_{\text{sft}}, $$

and we are left with bounding the two terms on the right-hand side.

(2) Bound on $\|u_Q - \tilde{u}_Q\|_{\text{sft}}$. We observe that

$$ \begin{aligned} \|\nabla_Q(u_Q - \tilde{u}_Q)\|_{H^{-1}(J;L^2(\Omega))} &+ \|\partial_{t,Q}(u_Q - \tilde{u}_Q)\|_{L^2(J;H^{-1}(\Omega))} \\ &\lesssim \|\tilde{\nabla}_Q(u_Q - \tilde{u}_Q)\|_{L^2(Q)} \lesssim s_Q(\hat{u}_Q, \hat{u}_Q)^{\frac{1}{2}}, \\ \|u_Q - \tilde{u}_Q\|_{L^{\infty}(J;L^2(\Omega))} &\lesssim h^{-\frac{1}{2}} \|u_Q - \tilde{u}_Q\|_{L^2(Q)} \lesssim h^{\frac{1}{2}} s_Q(\hat{u}_Q, \hat{u}_Q)^{\frac{1}{2}}, \end{aligned} $$

where we used (35), and an inverse trace inequality for the last estimate. As a consequence, owing to Lemma 8 and since $h \lesssim 1$, we infer that

$$ \|u_Q - \tilde{u}_Q\|_{\text{sft}} \lesssim s_Q(\hat{u}_Q, \hat{u}_Q)^{\frac{1}{2}} \lesssim \|\hat{\theta}_Q\|_{\#} + \|\delta\|_{J \times \bar{\omega}}. $$

(3) Bound on $\|u - \tilde{u}_Q\|_{\text{sft}}$. Since $u - \tilde{u}_Q \in H^1(Q)$, recalling Assumption 1, we can invoke the conditional stability estimate from Lemma 2. This gives

$$ \|u - \tilde{u}_Q\|_{\text{sft}} \lesssim C_{\text{stb}}(\|u - \tilde{u}_Q\|_{J \times \bar{\omega}} + \|L(u - \tilde{u}_Q)\|_{H^{-1}(Q)}). $$
It remains to bound $\|u - \tilde{u}_Q\|_{J \times \bar{\omega}}$ and $\|L(u - \tilde{u}_Q)\|_{H^{-1}(Q)}$.

(3.a) Bound on $\|u - \tilde{u}_Q\|_{J \times \bar{\omega}}$. We have

$$ \|u - \tilde{u}_Q\|_{J \times \bar{\omega}} \le \|u - I_Q(u)\|_{J \times \bar{\omega}} + \|I_Q(u) - u_Q\|_{J \times \bar{\omega}} + \|u_Q - \tilde{u}_Q\|_{J \times \bar{\omega}}, $$

where $I_Q(u)$ is defined in Section 3.2. We have

$$ \|u - I_Q(u)\|_{J \times \bar{\omega}} = \|\theta_Q\|_{J \times \bar{\omega}} \le \|\theta_Q\|_Q, $$

$$ \|I_Q(u) - u_Q\|_{J \times \bar{\omega}} = \|e_Q\|_{J \times \bar{\omega}} \le \|\hat{e}_Q, \hat{\xi}_Q\|. $$

Furthermore, the bound (35) and $h \lesssim 1$ imply that

$$ \|u_Q - \tilde{u}_Q\|_{J \times \bar{\omega}} \lesssim h s_Q(\hat{u}_Q, \hat{u}_Q)^{\frac{1}{2}} \lesssim s_Q(\hat{u}_Q, \hat{u}_Q)^{\frac{1}{2}}. $$

Gathering the above estimates and using Lemma 8 gives

$$ \|u - \tilde{u}_Q\|_{J \times \bar{\omega}} \lesssim \|\hat{\theta}_Q\|_{\#} + \|\delta\|_{J \times \bar{\omega}}. $$

(3.b) Bound on $\|L(u - \tilde{u}_Q)\|_{H^{-1}(Q)}$. Lemma 9 gives

$$ \|L_Q(u - u_Q)\|_{H^{-1}(Q)} \lesssim \|\hat{\theta}_Q\|_{\#} + \|\delta\|_{J \times \bar{\omega}}. $$

Moreover, using the Cauchy-Schwarz inequality and (35), we have

$$ \begin{aligned} \|L_Q(u_Q - \tilde{u}_Q)\|_{H^{-1}(Q)} &= \sup_{\substack{\eta \in H_0^1(Q) \\ \|\tilde{\nabla}\eta\|_Q = 1}} (B\tilde{\nabla}_Q(u_Q - \tilde{u}_Q), \tilde{\nabla}\eta)_Q \\ &\le \|B\tilde{\nabla}_Q(u_Q - \tilde{u}_Q)\|_Q \lesssim s_Q(\hat{u}_Q, \hat{u}_Q)^{\frac{1}{2}}. \end{aligned} $$

(4) Combining the bounds from the above steps proves (34). □

## 3.6 Convergence

We consider the functional space

$$ U_* := H^{\ell+1}(J; H^2(\Omega)) \cap H^2(J; H^{k+1}(\Omega)), \quad (36) $$

equipped with its natural norm. An important consequence of Lemma 6 is that, under the assumption $u \in U_*$, we have

$$ \|\hat{\theta}_Q\|_{\#} \lesssim h^{\ell} \|u\|_{H^{\ell+1}(J; H^2(\Omega))} + h^k \|u\|_{H^2(J; H^{k+1}(\Omega))}, \quad (37) $$

where we used that $h \lesssim 1$ to simplify the expression.

**Theorem 12 (Decay rates).** Assume that $u \in U_*$. We have the following decay rates on the errors:

$$ \|\hat{e}_Q, \hat{\xi}_Q\| + s_Q(\hat{u}_Q, \hat{u}_Q)^{\frac{1}{2}} \lesssim h^{\ell} \|u\|_{H^{\ell+1}(J; H^2(\Omega))} + h^k \|u\|_{H^2(J; H^{k+1}(\Omega))} + \|\delta\|_{J \times \bar{\omega}}, \quad (38a) $$

$$ \|u - u_Q\|_{\text{sft}} \lesssim (1 + C_{\text{stb}}) (h^{\ell} \|u\|_{H^{\ell+1}(J; H^2(\Omega))} + h^k \|u\|_{H^2(J; H^{k+1}(\Omega))} + \|\delta\|_{J \times \bar{\omega}}). \quad (38b) $$

*Proof.* This is a consequence of plugging (37) into (32) and (34). □
# 4 Numerical results

In this section, we present numerical experiments to check the convergence rates established in Theorem 12. We also study the influence of noise in the measurements, and we illustrate the benefits of using a high-order discretization.

The space domain is $\Omega := (0, 1)^d$ with $d \in \{1, 2\}$. All the errors are computed as the difference between the numerical solution and the $L^2(Q)$-orthogonal projection of the exact solution onto the cell discretization space. These errors are measured in the norm $\| \cdot \|_{\text{sft}} := (\| \cdot \|_{\ell^\infty(J;L^2(\Omega))}^2 + \|\nabla \cdot \|_{L^2(Q)}^2 + \|\partial_t \cdot \|_{L^2(Q)}^2)^{\frac{1}{2}}$, which is an approximation of the shifted energy norm defined in (2), whereby the $H^{-1}$-norms (which are not accessible to computation) are replaced by $L^2$-norms and the $\ell^\infty(J; L^2(\Omega))$-norm is computed as the maximum value over all the Gauss points of each time interval.

Some noise can be added to the measurements in the following way: (i) the space-time domain is divided into $10^{d+1}$ subdomains; (ii) a random noise level $\delta := a * \text{rand()}$ is assigned to each subdomain, with $a \ge 0$ the noise amplitude and $\text{rand()}$ is a C++ function returning a random number in $[-1, 1]$. Thus, every subdomain has the same noise during the whole mesh-refinement process.

The time mesh is considered to be a uniform mesh composed of $N$ time cells.

All the tests are run with the `DisSk++` library [20], and all the linear systems are solved using the Pardiso solver from the MKL library. During the matrix assembly process, we compute the value of the coefficients for the first time interval. Then, since the time step is constant, the coefficients associated with the following time intervals are the same and therefore do not need to be recomputed.

## 4.1 One-dimensional test cases

In this section, we consider the following one-dimensional setting:

$$ \Omega := (0, 1), \quad \varpi := \left(\frac{1}{4}, \frac{3}{4}\right), \quad T_f := 2, \quad u(t, x) := \cos(\pi t) \sin(2\pi x). $$

We notice that this geometry fulfills Assumption 1, since every wave has to cross the set $J \times \varpi$ at some point. We also notice that the above solution corresponds to a nonzero source term.

We use a uniform mesh in space and in time ($N$ cells in time and $M$ cells in space). Four levels of refinement are considered in space ($M \in \{16, 32, 64, 128\}$) and in time ($N \in \{10, 20, 40, 80\}$). We run two convergence studies. At first, for a good precision in time ($N = 128$, $\ell = 3$), we study the error for several successive space mesh sizes $M \in \{16, 32, 64, 128\}$ and several polynomial orders $k \in \{1, 2, 3\}$. Then, for a good precision in space ($M = 256$, $k = 3$), we study the error for several successive time mesh sizes $N \in \{10, 20, 40, 80\}$ and several polynomial orders $\ell \in \{1, 2, 3\}$. The results without noise are reported in Figure 2. We observe optimal space convergence at rate $k$ for $k \in \{2, 3\}$. Instead, we obtain superconvergence at rate 2 for $k = 1$. Moreover, we also observe the expected time convergence at rate $\ell$, for $\ell \in \{1, 2, 3\}$. (Recall that the errors are computed as the difference between the numerical solution and a projection of the exact solution.) Another relevant observation is that a higher-order scheme is more efficient in terms of degrees of freedom.
<figure><img src="image_4.png" /><figcaption>Figure 2: $\|\cdot\|_{\text{sft}}$-norm errors for the 1d test case without noise (with respect to the number of degrees of freedom). Left: space convergence ($M \in \{16, 32, 64, 128\}$, $N = 128$, $\ell = 3$). Right: time convergence ($N \in \{10, 20, 40, 80\}$, $M = 256$, $k = 3$).</figcaption></figure>

<table><caption>Table 1: Number of degrees of freedom (1d test case) for space refinement with $\ell = 3$ and $N = 128$. Left: without static condensation. Right: with static condensation.</caption><thead><tr><th colspan="4">Without static condensation</th><th colspan="4">With static condensation</th></tr><tr><th><i>M</i></th><th><i>k</i> = 1</th><th><i>k</i> = 2</th><th><i>k</i> = 3</th><th><i>M</i></th><th><i>k</i> = 1</th><th><i>k</i> = 2</th><th><i>k</i> = 3</th></tr></thead><tbody><tr><td>16</td><td>56320</td><td>76800</td><td>97280</td><td>16</td><td>23552</td><td>27648</td><td>31744</td></tr><tr><td>32</td><td>113664</td><td>154624</td><td>195584</td><td>32</td><td>48128</td><td>56320</td><td>64512</td></tr><tr><td>64</td><td>228352</td><td>310272</td><td>392192</td><td>64</td><td>97280</td><td>113664</td><td>130048</td></tr><tr><td>128</td><td>457728</td><td>621568</td><td>785408</td><td>128</td><td>195584</td><td>228352</td><td>261120</td></tr></tbody></table>

Notice that the numerical error has also been computed using the $L^2(Q)$-norm and the $H^1(\Omega)$-seminorm at initial time (using the associated face element). We observed for these errors the same convergence rate as the one for the $\|\cdot\|_{\text{sft}}$-norm used above. For the sake of clarity, these errors have not been drawn on Figure 2.

We report in Table 1 the number of degrees of freedom in the two situations of whether or not the static condensation procedure is used (see Section 2.4). We can observe a drastic reduction in the number of degrees of freedom offered by the static condensation procedure.

## 4.2 Two-dimensional test case

In this section, we consider the following two-dimensional setting:

$$ \Omega := (0, 1)^2, \quad \varpi := \{(b, 1) \times (0, 1)\} \setminus \{(b, \frac{7}{8}) \times (\frac{1}{8}, \frac{7}{8})\}, \quad T_f := 2, \\ u(t, x, y) := \cos(\sqrt{2\pi\omega t}) \sin(\pi\omega x) \sin(\pi\omega y), $$

where $b \in [0, \frac{7}{8})$ is a parameter that may vary and will be used to consider several geometric configurations, and $\omega := 3$.

One can verify that with $T_f := 2$ and $b := 0$, the geometry fulfills the geometric control condition (Assumption 1). However, this assumption is not fulfilled for $b > 0$.
This latter configuration will enable us to evaluate the impact of this assumption on the approximation results.

<figure><img src="image_5.png" /><figcaption>Figure 3: First, third, fifth, and seventh triangulations used for the 2d test case ($b = \frac{1}{2}$); all the triangulations are fitted to $\varpi$.</figcaption></figure>

The errors in the $\|\cdot\|_{\text{sft}'}$-norm are presented in Figure 4 in the absence of noise on the measurements. On the left panel, we provide a convergence study with respect to the space discretization. Seven successive meshes are considered (see Figure 3) and the space polynomial degrees are $k \in \{1, 2, 3\}$. For this test, the time discretization remains constant ($\ell = 3, N = 40$). Keeping in mind that the x-axis of this figure is the number of degrees of freedom (that scales essentially as $h^{-2}$ when we refine the space discretization), we retrieve the expected space convergence rates for the various values $k \in \{2, 3\}$. We also observe that the high-order approach is much more efficient than the low-order one, which is a known behavior for the wave equation.

<figure><img src="image_6.png" /><figcaption>Figure 4: $\|\cdot\|_{\text{sft}'}$-norm errors for the 2d test case without noise (with respect to the number of degrees of freedom). Left: space convergence. Right: time convergence.</figcaption></figure>

On the right panel, we provide a study of the time convergence for $\ell \in \{1, 2, 3\}$ and for $N \in \{5, 10, 15, 20, 30, 40\}$. As stated before for the space convergence study, we observe here that the high-order approach is much more efficient. Surprisingly, the low-order approach ($k = 1$ or $\ell = 1$) does not seem to converge for the considered meshes.

In both of these tests, we observe that the configuration with $b = \frac{1}{2}$ is less favorable than the configuration with $b = 0$. However, convergence at optimal rate is still observed in the configuration with $b = \frac{1}{2}$. The only visible drawback of that configuration is that the
<table><caption>Table 2: 2d test case for space refinement with $\ell = 3$ and $N = 40$ ($b = \frac{1}{2}$) and without noise. Left: computational time (in sec). Right: memory usage tip (in GB).</caption><thead><tr><th colspan="4">Computational time (in sec)</th><th colspan="4">Memory usage (in GB)</th></tr><tr><th>mesh</th><th>k = 1</th><th>k = 2</th><th>k = 3</th><th>mesh</th><th>k = 1</th><th>k = 2</th><th>k = 3</th></tr></thead><tbody><tr><td>1</td><td>8.6</td><td>39</td><td>100</td><td>1</td><td>10<sup>-3</sup></td><td>4.1</td><td>8.8</td></tr><tr><td>3</td><td>38</td><td>170</td><td>440</td><td>3</td><td>4.3</td><td>11.6</td><td>23.7</td></tr><tr><td>5</td><td>610</td><td>2500</td><td>6300</td><td>5</td><td>20.9</td><td>54.5</td><td>107</td></tr><tr><td>7</td><td>13000</td><td>54000</td><td>—</td><td>7</td><td>135</td><td>345</td><td>—</td></tr></tbody></table>

scaling of the errors, but not their decay rate, seems to be a bit worse. We also considered adding some noise to the right-hand side $f$, but this resulted only in a slight degradation of the solution. The results are not reported for the sake of conciseness.

The same setting is then considered in the presence of noise ($a = 10^{-2}$). The results are reported in Figure 5. We observe that the main consequence of noise is that convergence can be obstructed since the error stagnates at a certain level which depends on the geometry ($b = 0$ or $b = \frac{1}{2}$).

<figure><img src="image_7.png" /><figcaption>Figure 5: $\|\cdot\|_{\text{sft}}$-norm errors for the 2d test case with noise $a = 10^{-2}$ (with respect to the number of degrees of freedom). Left: space convergence. Right: time convergence.</figcaption></figure>

Computational time and memory usage for the configuration $b = \frac{1}{2}$ without noise are reported in Table 2. We observe that, since a direct solver is used, we attain the limits of the method in terms of memory usage. More refined meshes or 3d-space configurations need the development of an iterative solver.

**Acknowledgement:** GD was supported by the Emergence grant of Sorbonne University. EB was supported by EPSRC grants EP/T033126/1 and EP/V050400/1.

## Appendix A Technical proof

Here we prove Lemma 2. The proof is obtained by combining Proposition A.1 and Theorem A.4 from [18] where $\Omega$ is assumed to be smooth. Proposition A.1 yields that for $v = 0$ on
$J \times \partial\Omega$ we have

$$
\begin{align*}
\|v\|_{L^{\infty}(J;L^2(\Omega))} + \|\partial_t v\|_{L^2(J;H^{-1}(\Omega))} \\
&\lesssim \|v|_{t=0}\|_{L^2(\Omega)} + \|\partial_t v|_{t=0}\|_{H^{-1}(\Omega)} + \|L(v)\|_{H^{-1}(J \times \Omega)}.
\end{align*}
$$

Moreover, following the proof of Theorem A.4 from [18], we can show that under the
geometric control condition (Assumption 1), we have

$$
\|v|_{t=0}\|_{L^2(\Omega)} + \|\partial_t v|_{t=0}\|_{H^{-1}(\Omega)} \lesssim \|L(v)\|_{H^{-1}(J \times \Omega)} + \|v\|_{L^2(J \times \bar{\omega})}.
$$

Combining these two relations, we have

$$
\|v\|_{L^{\infty}(J;L^2(\Omega))} + \|\partial_t v\|_{L^2(J;H^{-1}(\Omega))} \lesssim \|v\|_{L^2(J \times \bar{\omega})} + \|L(v)\|_{H^{-1}(Q)}.
$$

To estimate the missing term, we consider

$$
\|\nabla v\|_{H^{-1}(J;L^2(\Omega))} \lesssim \|\Delta v\|_{H^{-1}(Q)} \le \|L(v)\|_{H^{-1}(Q)} + \|\partial_{tt}v\|_{H^{-1}(Q)},
$$

where, in the first bound, we used $\|\nabla\phi\|_{\Omega} \lesssim \|\Delta\phi\|_{H^{-1}(\Omega)}$ for all $\phi \in H_0^1(\Omega)$. Moreover, for $\phi \in H_0^1(Q)$ we have

$$
\begin{align*}
\langle \partial_{tt} v, \phi \rangle_{H^{-1}(Q), H_0^1(Q)} &= - \langle \partial_t v, \partial_t \phi \rangle_{L^2(J; H^{-1}(\Omega)), L^2(J; H_0^1(\Omega))} \\
&\le \|\partial_t v\|_{L^2(J; H^{-1}(\Omega))} \|\partial_t \phi\|_{L^2(J; H_0^1(\Omega))}.
\end{align*}
$$

Hence, $\|\partial_{tt}v\|_{H^{-1}(Q)} \le \|\partial_t v\|_{L^2(J;H^{-1}(\Omega))}$. This ends the proof.

References

[1] C. Bardos, G. Lebeau, and J. Rauch. Sharp sufficient conditions for the observation, control, and stabilization of waves from the boundary. *SIAM J. Control Optim.*, 30 (5):1024–1065, 1992.

[2] L. Baudouin, M. de Buhan, S. Ervedoza, and A. Osses. Carleman-based reconstruc-
tion algorithm for waves. *SIAM J. Numer. Anal.*, 59(2):998–1039, 2021.

[3] E. Bécache, L. Bourgeois, L. Franceschini, and J. Dardé. Application of mixed formu-
lations of quasi-reversibility to solve ill-posed problems for heat and wave equations:
the 1D case. *Inverse Probl. Imaging*, 9(4):971–1002, 2015.

[4] E. Burman. Stabilized finite element methods for nonsymmetric, noncoercive, and ill-posed problems. Part I: Elliptic equations. *SIAM J. Sci. Comput.*, 35(6):A2752–A2780, 2013.

[5] E. Burman. Error estimates for stabilized finite element methods applied to ill-posed problems. *C. R. Math. Acad. Sci. Paris*, 352(7-8):655–659, 2014.

[6] E. Burman and A. Ern. Continuous interior penalty hp-finite element methods for advection and advection-diffusion equations. *Math. Comp.*, 76(259):1119–1140, 2007.
[7] E. Burman and P. Hansbo. Stabilized nonconforming finite element methods for data assimilation in incompressible flows. *Math. Comp.*, 87(311):1029–1050, 2018.

[8] E. Burman and L. Oksanen. Data assimilation for the heat equation using stabilized finite element methods. *Numer. Math.*, 139(3):505–528, 2018.

[9] E. Burman and J. Preuss. Unique continuation for the wave equation based on a discontinuous Galerkin time discretization. preprint, 2024. URL https://arxiv.org/pdf/2405.04615. submitted.

[10] E. Burman, J. Ish-Horowicz, and L. Oksanen. Fully discrete finite element data assimilation method for the heat equation. *ESAIM Math. Model. Numer. Anal.*, 52 (5):2065–2082, 2018.

[11] E. Burman, M. Nechita, and L. Oksanen. Unique continuation for the Helmholtz equation using stabilized finite element methods. *J. Math. Pures Appl. (9)*, 129:1–22, 2019.

[12] E. Burman, A. Feizmohammadi, and L. Oksanen. A finite element data assimilation method for the wave equation. *Math. Comp.*, 89(324):1681–1709, 2020.

[13] E. Burman, A. Feizmohammadi, and L. Oksanen. A fully discrete numerical control method for the wave equation. *SIAM J. Control Optim.*, 58(3):1519–1546, 2020.

[14] E. Burman, M. Cicuttin, G. Delay, and A. Ern. An unfitted hybrid high-order method with cell agglomeration for elliptic interface problems. *SIAM J. Sci. Comput.*, 43(2): A859–A882, 2021.

[15] E. Burman, G. Delay, and A. Ern. A hybridized high-order method for unique continuation subject to the Helmholtz equation. *SIAM J. Numer. Anal.*, 59(5):2368–2392, 2021.

[16] E. Burman, A. Feizmohammadi, A. Münch, and L. Oksanen. Space time stabilized finite element methods for a unique continuation problem subject to the wave equation. *ESAIM Math. Model. Numer. Anal.*, 55(suppl.):S969–S991, 2021.

[17] E. Burman, G. Delay, and A. Ern. The unique continuation problem for the heat equation discretized with a high-order space-time nonconforming method. *SIAM J. Numer. Anal.*, 61(5):2534–2557, 2023.

[18] E. Burman, A. Feizmohammadi, A. Münch, and L. Oksanen. Spacetime finite element methods for control problems subject to the wave equation. *ESAIM Control Optim. Calc. Var.*, 29:Paper No. 41, 40, 2023.

[19] N. Burq and P. Gérard. Condition nécessaire et suffisante pour la contrôlabilité exacte des ondes. *C. R. Acad. Sci. Paris Sér. I Math.*, 325(7):749–752, 1997.

[20] M. Cicuttin, D. A. Di Pietro, and A. Ern. Implementation of Discontinuous Skeletal methods on arbitrary-dimensional, polytopal meshes using generic programming. *J. Comput. Appl. Math.*, 2017.
[21] N. Cîndea and A. Münch. Inverse problems for linear hyperbolic equations using mixed formulations. *Inverse Problems*, 31, 2015.

[22] N. Cîndea and A. Münch. A mixed formulation for the direct approximation of the control of minimal $L^2$-norm for linear type wave equations. *Calcolo*, 52(3):245-288, 2015.

[23] W. Dahmen, H. Monsuur, and R. Stevenson. Least squares solvers for ill-posed PDEs that are conditionally stable. *ESAIM Math. Model. Numer. Anal.*, 57(4):2227-2255, 2023.

[24] A. Ern and J.-L. Guermond. Finite element quasi-interpolation and best approximation. *ESAIM Math. Model. Numer. Anal.*, 51(4):1367-1385, 2017. ISSN 2822-7840.

[25] S. Ervedoza and E. Zuazua. *Numerical approximation of exact controls for waves*. SpringerBriefs in Mathematics. Springer, New York, 2013.

[26] P. Houston, D. Schötzau, and T. P. Wihler. Energy norm a posteriori error estimation of hp-adaptive discontinuous Galerkin methods for elliptic problems. *Math. Models Methods Appl. Sci.*, 17(1):33-62, 2007. ISSN 0218-2025.

[27] O. A. Karakashian and F. Pascal. A posteriori error estimates for a discontinuous Galerkin approximation of second-order elliptic problems. *SIAM J. Numer. Anal.*, 41(6):2374-2399, 2003. ISSN 0036-1429.

[28] J. Le Rousseau, G. Lebeau, P. Terpolilli, and E. Trélat. Geometric control condition for the wave equation with a time-dependent observation domain. *Anal. PDE*, 10(4): 983-1015, 2017.