We consider a stable unique continuation problem for the wave equation where the initial data is lacking and the solution is reconstructed using measurements in some subset of the bulk domain. Typically fairly sophisticated space-time methods have been used in previous work to obtain stable and accurate solutions to this reconstruction problem. Here we propose to solve the problem using a standard discontinuous Galerkin method for the temporal discretization and continuous finite elements for the space discretization. Error estimates are established under a geometric control condition. We also investigate two preconditioning strategies which can be used to solve the arising globally coupled space-time system by means of simple time-stepping procedures. Our numerical experiments test the performance of these strategies and highlight the importance of the geometric control condition for reconstructing the solution beyond the data domain.

Unique continuation, data assimilation, wave equation, finite element method, discontinuous Galerkin, preconditioning, geometric control condition

65M32, 35R30, 35L05, 65F08

§ INTRODUCTION

In this article we consider a Lipschitz-stable unique continuation problem for the wave equation.
For a bounded domain $\Omega \subset \mathbb{R}^d$, where $d \in \mathbb{N}$, with smooth strictly convex boundary $\partial \Omega$, and $T > 0$ we define the space-time
cylinder $\STdom = (0,T) \times \Omega$ with lateral boundary $\Sigma := (0,T) \times \partial \Omega$.
Analogously, for a subset $\omega \subset \Omega$ we define the data domain $ \STdata := (0,T) \times \omega$.
We are now concerned with the following data assimilation problem:
Find $u:\STdom \rightarrow \mathbb{R}$ such that

\begin{align}
(\Box u, u)  &= (0,0) \text{ in } \STdom \times \Sigma,    \label{eq:PDE+lateralBc} \\
u  &= u_{\omega} \text{ in }  \STdata, \label{eq:data-constraint}
\end{align}

where $\Box u :=\partial_t^2 u - \triangle u$ and $u_{\omega}$ is the known data. Note that only the restriction of
$u$ to the data set is given, whereas the wave-displacement and velocity at time $t=0$ in $\Omega$ are

Problem <ref>-<ref> enjoys Lipschitz stability if the
set $ \STdata \subset \STdom$ fulfills the so called geometric control condition (GCC), see [3, 4] for a precise definition thereof.

Assume that $ \STdata \subset \STdom$ fulfills the GCC.
Then there exists a constant $C > 0$ such that for any $\phi \in H^1(\STdom)$ the following stability estimate holds:
ϕ_L^∞(0,T;L^2(Ω)) + ∂_t ϕ_L^2(0,T;H^-1(Ω))
≤C (  ϕ_ L^2()  + ϕ_ L^2(Σ) +  ϕ_ H^-1 ( )      ).

See <cit.>.

Unless otherwise stated, below we assume that the GGC holds for the chosen geometry.

§.§ Literature review

Recently, at least five different methods [13, 9, 7, 19, 14] for the numerical solution of problem <ref>-<ref> have been proposed.
Before discussing these approaches let us briefly explain a common principle that drives the design of all numerical methods based on the stability estimate <ref>. To this end, let $u_h \in H^1(\STdom)$ be some discrete approximation of the solution $u$ of problem <ref>-<ref>.
The idea is then to apply <Ref> to $\phi = u - u_h$, which leads to an error estimate provided that the terms on the right hand side of <ref>, in particular the PDE residual $\norm{ \Box (u - u_h ) }_{ H^{-1} ( \STdom ) }$, can be controlled. Let us distinguish the recently proposed methods based on how they strive to
achieve this aim.

* The first method proposed in [13] in fact preceeded the stability estimate in <ref>. Here the discretization was required to fulfill
$\Box u_h \in L^2(0,T,H^{-1}(\Omega))$ to repect the PDE constraint. To this end, $C^1$-conforming (space-time) elements have been used. Verifying discrete inf-sup stability in
this setting is an intricate problem. Alternatively, the authors propose a stabilized formulation which is shown to convergence under an additional regularity assumption on
an auxiliary variable, see <cit.>.

* Subsequently, a class of methods [9, 7] based on a primal-dual stabilized finite element framework introduced in [5] were proposed. Here, the problem is formulated as a PDE constrained optimization problem and discretized. Then, (weakly)-consistent stabilization terms are added to the discretization which allow to control the PDE residual without having to evaluate dual norms or to use exotic discretization spaces.
Based on this strategy, in reference [9] a piecewise affine continuous finite element method in space was coupled with backward differences in time to yield a provably convergent
first order method for problem <ref>-<ref>.
We point out that even though backward differences in time were used, this is not a time-stepping method, as the solution of a globally coupled space-time system is required. Since the method of [9] ultimately leads to space-time systems anyhow, it seems natural to treat problem <ref>-<ref>
directly using a full space-time discretization. In fact, such an approach has been considered in reference [7]. Therein the authors proposed an arbitrary order $H^1(\STdom)$-conforming
stabilized finite element method and derived optimal error estimates. An extension of the method to a closely related control problem has been presented in [8].

* Whereas the methods introduced in [9, 7] avoid a direct evaluation of the PDE residual $\norm{ \Box (u - u_h ) }_{ H^{-1} ( \STdom ) }$, in reference [14] a least-squares finite element method has been introduced which minimizes this quantity directly. The authors deal with the arising challenges of having to prove technical inf-sup conditions and the need to replace the negative index Sobolev norms by computable quantities to be able to solve the problem numerically. We remark that the framework introduced in this reference is fairly general which allows the additional treatment of data assimilation for the heat equation and the Cauchy problem for Poisson's equation in this publication.

* Finally, physics informed neural networks have been proposed in reference [19] to solve problem <ref>-<ref> and similar data assimilation problems. Here, the PDE residual is minimized in the $L^2$ instead of the $H^{-1}$-norm on a set of training points and  estimates on the generalization error of these networks have been derived based on <Ref>. Let us mention that, even though this approach in its present form can yield accurate results for problems enjoying robust stability, it has recently been shown that it is not suitable yet to solve Helmholtz problems posed in geometrical configurations with poor stability, see reference [20].

§.§ Motivation for and novelty of this paper

The methods discussed in <ref> share the common feature of treating problem <ref>-<ref> in a space-time fashion.
This is reasonable since initial data $u(0,\cdot), u_t(0,\cdot)$ for starting a time-marching procedure is not available.
Also, in order to apply <Ref> to the error, one needs to ensure that $u_h$ is in $H^1(\STdom)$, which renders the use of conforming space-time discretizations very natural and appealing.
However, space-time methods also come with some considerable challenges. Applying for example the method from [7] to a domain $\Omega \subset \mathbb{R}^3$ would
require to generate a non-trivial (recall that $\partial \Omega$ is smooth) four-dimensional space-time mesh and a corresponding discretization which can be challenging to realize with common finite element software. Whereas this algorithmic challenge may be overcome, solving the arising extremely large and ill-conditioned linear systems appears like
a more daunting obstacle. Therefore, we want to investigate in this article to what extent problem <ref>-<ref> is amenable to the solution with computationally less-demanding time-marching methods. To this end, we consider a fairly flexible time discretization provided by a discontinuous Galerkin method (dG), which is very well-established in the literature, see e.g. [16, 17, 15]. The dG method is known to allow for time-marching procedures provided that only weak continuity in time is imposed by the variational formulation. The analysis presented in this article shows however, that some stronger constraints which couple the problem globally in time cannot be avoided if one wants to preserve the optimal error estimates from reference [7]. The contribution of this manuscript is to identify these couplings explictly and then to introduce small perturbations thereof for which time-marching is possible.
These approximate solvers can then be used as preconditioners to solve the original globally-coupled space-time problem iteratively.
The advantage of this strategy is that it allows to solve problem <ref>-<ref> using a standard time-marching procedure and that it avoids space-time discretizations. Now we can compute solutions in spatially-three dimensional domains ($d=3$) up to moderate levels of precision. The current limitations of our approach are that the preconditioners lack robustness with respect refinement of the discretization and that adaptive refinement strategies have not been considered yet.
We remark that a dG method in time for an ill-posed heat equation has already been analyzed in reference [6]. However, the authors did not make use of the dG in time structure for solving the linear systems. This idea appears novel to us.

§.§ Structure of this paper

To show the generality of the approach with respect to time discretisation we consider both a semi-discretisation in space, with time continuous and a fully-discrete formulation with a discontinuous Galerkin approximation in time.
We introduce these methods in <ref> and establish some auxiliary results valid for both approaches. The semi-discrete method in which time is still considered to be a continuous variable is analyzed in <ref>, whereas analysis of the fully-discrete method can be found in <ref>. The separate treatment of these two approaches allows us to identify very clearly the mechanism which couples provably convergent approximations of <ref>-<ref> globally in time at the discrete level.
Based on these insights, we propose two preconditioning strategies in <ref> which can be used for the iterative solution of the arising linear systems. In <ref> we investigate their performance in numerical experiments and explore further aspects which go beyond our analytical results. For instance, we investigate how the numerical method behaves when the GCC (which is crucial for the analysis) is violated. This article is accompanied by a <ref> in which interpolation into tensor-product space-time finite element spaces is discussed.

§ A SEMI-DISCRETE AND A FULLY-DISCRETE METHOD FOR UNIQUE CONTINUATION

In this section we introduce two levels of discretization for problem <ref>-<ref>.
We start in <ref> by defining a geometric partition of the the space-time cylinder $Q$
into time-slabs, which are tensor products of subintervals of $(0,T)$ with a simplicial mesh covering the spatial
domain $\Omega$. In <ref> we then introduce a preliminary method to solve problem <ref>-<ref>, which
features an ansatz space that is discrete in the spatial variable yet still continuous in time.
In <ref> we drop the latter assumption by introducing a
discontinuous Galerkin time discretization. This requires a slight enrichment of the method introduced in <ref> to impose sufficient temporal continuity
by means of the variational formulation. In <ref> and  <ref> we establish some theoretical results valid for both discretizations which are required for the error analysis of the respective methods presented in <ref> and <ref>.

§.§ Discretization of the space-time cylinder

§.§.§ Spatial mesh

For <ref> to hold we require that $\Omega$ has a smooth boundary. To account for this in a simple manner in the spatial disretization we use the same strategy as discussed in <cit.>. Essentially, we fit the boundary using curved elements and incorporate the boundary conditions weakly using Nitsche's method, similar as in <cit.>. Let $\hat{\mathcal{T}}_h$ be a quasi-uniform triangulation such that $\hat{\Omega}:= \cup_{K \in \hat{\mathcal{T}}_h} K $ contains $\Omega$, i.e. $\Omega \subset \hat{\Omega}$. We then set $
\mathcal{T}_h := \{ K \cap \Omega \mid K \in  \hat{\mathcal{T}}_h \}.$
Note that by construction $\Omega = \cup_{K \in \mathcal{T}_h} K$ holds true.
As explained in <cit.>, it is possible to construct the family $\{ \hat{\mathcal{T}}_h \mid h > 0 \}$ such that for sufficiently small $h$ the  continuous trace inequality
v_ [L^2(∂K)]^d  ≤C ( h^-1/2 v_ [L^2(K)]^d  + h^1/2  ∇v  _ [L^2(K)]^d    ),   ∀v ∈[H^1(K)]^d,   K ∈𝒯_h
holds, where here as in the remainder of this article $C$ denotes a generic constant independent of $h$. We refer to <cit.> for a proof.
We will in the following also assume that the triangulation $\mathcal{T}_h$ fits the data domain $\omega$.
Denote then by $\mathcal{F}_i$ the set of all interior facets of $\mathcal{T}_h$ and let

\[
V_h^k := \{ u \in H^1(\Omega) \mid u|_K \in \mathcal{P}_k(K) \; \forall K \in \mathcal{T}_h \},
\]

where $\mathcal{P}_k(K)$ denotes the spaces of polynomials of order at most $k \in \mathbb{N}$ on $K$.

§.§.§ Partition of the time axis

We introduce a partition $0 = t_0 < t_1 < \ldots < t_N=T$ of the time axis into subintervals $I_n := (t_n,t_{n+1}), n=0,\ldots,N-1$.
For simplicity all time intervals are assumed to be of equal length $\Delta t = \abs{t_{n+1} - t_n}$.
The discretization of the time axis implies a corresponding partition of the space-time cylinder into time slabs
Q^n := I_n ×Ω,   Σ^n = I_n ×Σ,   n=0,…,N-1,   Q := ∪_n=0^N-1 Q^n,   Σ:= ∪_n=0^N-1 Σ^n.
The space-time integrals on the slabs are denoted by

\[
(u,v)_{Q^n} := \int\limits_{I_n} \int\limits_{\Omega} u v \; \dX \; \dT,
\;
a(u,v)_{Q^n} := \! \int\limits_{I_n} \int\limits_{\Omega} \nabla u \cdot \nabla v \; \dX \; \dT,
\;
(u,v)_{ \Sigma^n } = \int\limits_{I_n}  \int\limits_{ \partial \Omega} u v \; \dS \; \dT
\]

and we set $\norm{v}^2_{Q^n} := (v,v)_{Q^n}$ and $ \norm{v}^2_{\Sigma^n} = (v,v)_{\Sigma^n}$.
The data domain and integrals thereon are defined similarly

\[
\omega^n := I_n \times \omega, \quad (u,v)_{\omega^n} := \int\limits_{I_n} \int\limits_{\omega} u v \; \dX \; \dT,   \quad \norm{ u }_{ \STdata }^2 = \sum\limits_{n=0}^{N-1}(u,v)_{\omega^n}.
\]

§.§ A semi-discrete conforming FEM

Let us first ignore the partition of the time axis in our ansatz space and define
k := H^1 ( 0,T; V_h^k ).
For tuples in the product space $ \ProdSemiDiscrSpace{k} = \SemiDiscrSpace{k} \times \SemiDiscrSpace{k} $ it is convenient to use the notation
𝐔_𝐡 := ( u_1,u_2) ∈k.

Let us now proceed to define a variational method for solving problem <ref>-<ref> in these semi-discrete spaces. It is convenient to define this method already in a time-slab wise manner, since we will later use a slightly augmented version thereof once the time variable has been discretized by a discontinuous Galerkin method.
Define the bilinear form representing a mixed formulation of the wave equation as

\begin{equation}\label{eq:bil-form-def-fully-discrete}
A[ \underline{\mathbf{U}}_h, \underline{\mathbf{Y}}_h] = \! \sum\limits_{n=0}^{N-1} \! \left\{ (\partial_t \underline{u}_2, \underline{y}_1)_{Q^n} + a(\underline{u}_1,\underline{y}_1)_{Q^n}
+ (\partial_t \underline{u}_1 - \underline{u}_2, \underline{y}_2)_{Q^n}
- ( \nabla \underline{u}_1 \! \cdot \! \mathbf{n}, \underline{y}_1)_{ \Sigma^n }
\right\}
\end{equation}

for $\underline{\mathbf{U_h}} \in \ProdSemiDiscrSpace{k}$ and $\underline{\mathbf{Y}}_h \in $ for $\ProdSemiDiscrSpace{k_{\ast}}$ with $k_{\ast} \in \mathbb{N}$. Here $\mathbf{n}$ denotes the outer normal vector of $\Omega$.
Note that for a sufficiently regular solution $u$ of <ref> we obtain for $ \mathbf{U} = (u, \partial_t u)$ via integration by parts in space that $A[ \mathbf{U} , \underline{\mathbf{Y}}_h] = 0$.
We will search for approximate solutions of problem <ref>-<ref> as stationary points of the Lagrangian
ℒ_h( 𝐔_h,𝐙_h) :=
1/2  u_1 - u_ω_ω_T^2 +  A[ 𝐔_h, 𝐙_h ] + 1/2 S_h( 𝐔_h, 𝐔_h)  - 1/2 S_h^∗( 𝐙_h, 𝐙_h),
for $(\underline{\mathbf{U}}_h,\underline{\mathbf{Z}}_h) \in \ProdSemiDiscrSpace{k} \times \ProdSemiDiscrSpace{k_{\ast}}$.
Here, the first and second term of the Lagrangian incorporate the data and PDE constraints, respectively. The contributions $ S_h(\cdot,\cdot)$ and $S_h^{\ast}(\cdot,\cdot)$ represent stabilization terms which are defined as follows. We define a spatial continuous interior penalty term $J(\cdot,\cdot)$, which penalizes the jump of the gradient $\jump{ \nabla \underline{u}_1 }$ over interior facets, a Galerkin-least squares term $G(\cdot,\cdot)$ which enforces the PDE at the element level, a term $I_{0}(\cdot,\cdot)$ which ensures that $\underline{u}_2 = \partial_t \underline{u}_1$ and a term for boundary stability:

\begin{align*}
& J( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h ) := \sum\limits_{n=0}^{N-1} \int\limits_{I_n} \sum\limits_{F \in \mathcal{F}_i} h  ( \jump{ \nabla \underline{u}_1 } , \jump{ \nabla \underline{w}_1 }  )_{F} \; \dT,
\\
& G( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h ) := \sum\limits_{n=0}^{N-1} \int\limits_{I_n} \sum\limits_{K \in \mathcal{T}_h }  h^2 ( \partial_t \underline{u}_2 - \Delta \underline{u}_1 , \partial_t \underline{w}_2 - \Delta \underline{w}_1 )_{K}  \; \dT, \quad
\\
& I_0( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h ) := \sum\limits_{n=0}^{N-1}(\underline{u}_2 - \partial_t \underline{u}_1, \underline{w}_2 - \partial_t \underline{w}_1)_{Q^n},
\quad  R( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h ) := \sum\limits_{n=0}^{N-1} h^{-1} ( \underline{u}_1, \underline{w}_1)_{ \Sigma^n }.
\end{align*}

The complete primal stabilization is given as a sum of these terms:

\begin{equation}\label{eq:Def-Sh}
S_{h}( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h ) :=
J( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h )
+ G( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h )
+ R( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h )
+ I_0( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h ).
\end{equation}

To stabilize the dual variable $ \underline{\mathbf{Z}}_h$ we will use:

\begin{equation}\label{eq:dual_stab_fully_disc}
S^{\ast}_h( \underline{\mathbf{Y}}_h, \underline{\mathbf{Z}}_h) :=
\sum\limits_{n=0}^{N-1} \big\{ (\underline{y}_1,\underline{z}_1)_{Q^n} + a(\underline{y}_1,\underline{z}_1)_{Q^n} + (\underline{y}_2,\underline{z}_2)_{Q^n} + h^{-1} (\underline{y}_1 ,\underline{z}_1 )_{ \Sigma^n }  \big\}.
\end{equation}

These stabilizations give rise to semi-norms, respectively norms:

\begin{equation}\label{eq:def-h-semi-norms}
\abs{ \underline{\mathbf{U}}_h }_{S_h} := S_h( \underline{\mathbf{U}}_h, \underline{\mathbf{U}}_h)^{1/2}, \quad
\norm{ \underline{\mathbf{Z}}_h }_{S_h^{\ast}} := S_h^{\ast}( \underline{\mathbf{Z}}_h, \underline{\mathbf{Z}}_h)^{1/2}.
\end{equation}

We will later see that the following expression defines a norm on $\ProdSemiDiscrSpace{k} \times \ProdSemiDiscrSpace{k_{\ast}}$:
( 𝐔_h, 𝐙_h) ^2_h  :=   𝐔_h _S_h^2 +   u_1 ^2_ +   𝐙_h _S_h^∗^2.
The first order optimality conditions characterizing critical points of the Lagrangian take the form:
Find $ (\underline{\mathbf{U}}_h ,\underline{\mathbf{Z}}_h) \in  \ProdSemiDiscrSpace{k} \times \ProdSemiDiscrSpace{k_{\ast}}  $ such that

\begin{align}\label{eq:optimality_cond_semi}
(\underline{u}_1,\underline{w}_1)_{\STdata} + A[ \underline{\mathbf{W}}_h, \underline{\mathbf{Z}}_h ] +  S_h(  \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h ) &= (u_{\omega},\underline{w}_1)_{ \STdata}, \quad \forall \; \underline{\mathbf{W}}_h \in  \ProdSemiDiscrSpace{k}, \nonumber \\
A[ \underline{\mathbf{U}}_h, \underline{\mathbf{Y}}_h ] - S_h^{\ast}( \underline{ \mathbf{Y}}_h, \underline{\mathbf{Z}}_h)  &= 0, \quad \forall \;  \underline{\mathbf{Y}}_h \in \ProdSemiDiscrSpace{k_{\ast}}.
\end{align}

This can be written in a more compact way:
Find $ (\underline{\mathbf{U}}_h ,\underline{\mathbf{Z}}_h) \in  \ProdSemiDiscrSpace{k} \times \ProdSemiDiscrSpace{k_{\ast}} $ such that

\begin{equation}\label{eq:opt_compact_semi}
B_h[ ( \underline{\mathbf{U}}_h,\underline{\mathbf{Z}}_h), ( \underline{\mathbf{W}}_h ,\underline{\mathbf{Y}}_h ) ] = (u_{\omega_T}, \underline{w}_1)_{ \STdata},
\end{equation}

\begin{align}\label{eq:complete_bfi_def_semi}
B_h[ ( \underline{\mathbf{U}}_h,\underline{\mathbf{Z}}_h), ( \underline{\mathbf{W}}_h ,\underline{\mathbf{Y}}_h ) ]  := & (\underline{u}_1,\underline{w}_1)_{\STdata} + A[\underline{\mathbf{W}}_h, \underline{\mathbf{Z}}_h ] + S_h( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h ) \nonumber \\
& + A[\underline{\mathbf{U}}_h, \underline{\mathbf{Y}}_h ] - S_h^{\ast}( \underline{\mathbf{Y}}_h, \underline{\mathbf{Z}}_h).
\end{align}

In <ref> we will show
that the first component $\underline{u}_1$ of the solution of <ref> converges to the exact solution $u$ of <ref>-<ref>
in the norm defined by the left hand side of <ref> at a rate of $h^k$. Note, however, that the method <ref> is not yet suitable for implementation since the time variable is still assumed to be continuous.

§.§ A fully-discrete FEM using dG in time

To discretize time, let $\mathcal{P}^q(I_n)$ denote the space of polynomials up to degree $q \in \mathbb{N}_0$ on $I_n$ and define a discontinuous (in-time) finite element space:

\[ \FullyDiscrSpace{k}{q} :=  \otimes_{n=0}^{N-1} \mathcal{P}^q(I_n) \otimes V_h^k,
\quad \ProdFullyDiscrSpace{k}{q} := \FullyDiscrSpace{k}{q} \times \FullyDiscrSpace{k}{q},
\quad q \in \mathbb{N}_{0}, \; k \in \mathbb{N}. \]

Note that functions in $ \ProdFullyDiscrSpace{k}{q} $ are allowed to be discontinuous between time-slabs.
Hence, we require a notation for jumps in time.
v^n_±(x) := lim_s →0^+ v(x,t_n ±s),  and   v^n  := v_+^n - v_-^n.
Since continuity in time is now no longer ensured by the spaces, we will need to augment the semi-discrete Lagrangian $\mathcal{L}_h(\cdot,\cdot)$ from <ref> by a stabilization term $S^{\uparrow \downarrow}_{\Delta t}$ which imposes sufficient temporal regularity of the primal variable, i.e.
ℒ( 𝐔_h,𝐙_h) :=
ℒ_h( 𝐔_h,𝐙_h) +
1/2 S^↑↓_Δt( 𝐔_h, 𝐔_h ),
(𝐔_h,𝐙_h) ∈kq ×k_∗q_∗,

\begin{equation}\label{eq:def-jump-stab}
S^{\uparrow \downarrow}_{\Delta t}( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h )
:= \underline{I}_1^{\uparrow \downarrow}( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h ) +
\underline{I}_2^{\uparrow \downarrow}( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h )
\end{equation}

\begin{align*}
\underline{I}_1^{\uparrow \downarrow}( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h ) &:= \sum\limits_{n=1}^{N-1} \left\{ \frac{1}{\Delta t} ( \jump{ \underline{u}_1^n } , \jump{ \underline{w}_1^n } )_{ \Omega}
+ \Delta t ( \jump{ \nabla \underline{u}_1^n } , \jump{ \nabla \underline{w}_1^n } )_{ \Omega} \right\}, \\
\underline{I}_2^{\uparrow \downarrow}( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h ) &:= \sum\limits_{n=1}^{N-1} \frac{1}{\Delta t} ( \jump{ \underline{u}_2^n } , \jump{ \underline{w}_2^n } )_{ \Omega}.
\end{align*}

The first order optimality conditions then take the form:
Find $(\underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h) \in \ProdFullyDiscrSpace{k}{q} \times \ProdFullyDiscrSpace{ k_{\ast} }{ q_{\ast} } $ such that
B[ ( 𝐔_h,𝐙_h), ( 𝐖_h ,𝐘_h ) ] = (u_ω, w_1)_
for all $( \underline{\mathbf{W}}_h ,\underline{\mathbf{Y}}_h ) \in \ProdFullyDiscrSpace{k}{q} \times \ProdFullyDiscrSpace{ k_{\ast} }{ q_{\ast} } $,

\begin{equation}\label{eq:complete_bfi_def_fully_discr}
B[ ( \underline{\mathbf{U}}_h,\underline{\mathbf{Z}}_h), ( \underline{\mathbf{W}}_h ,\underline{\mathbf{Y}}_h ) ] =  B_h[ ( \underline{\mathbf{U}}_h,\underline{\mathbf{Z}}_h), ( \underline{\mathbf{W}}_h ,\underline{\mathbf{Y}}_h ) ]
+ S^{\uparrow \downarrow}_{\Delta t}( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h )
\end{equation}

with $B_h$ as defined in <ref> and we assume that $k,k_{\ast},q \geq 1$ and $q_{\ast} \geq 0$.   Since $\tnorm{(\cdot,\cdot)}_h$ no longer defines a norm on the fully discrete spaces $\ProdFullyDiscrSpace{k}{q} \times \ProdFullyDiscrSpace{ k_{\ast} }{ q_{\ast} } $, it has to be augmeneted as well. To this end, we define

\[
\tnorm{ (\underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h) }^2  :=  \tnorm{ (\underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h) }^2_h +\abs{ \underline{\mathbf{U}}_h   }_{ \uparrow \downarrow }^2, \qquad
\abs{ \underline{\mathbf{U}}_h   }_{ \uparrow \downarrow } := S^{\uparrow \downarrow}_{\Delta t}( \underline{\mathbf{U}}_h, \underline{\mathbf{U}}_h )^{1/2}.
\]

§.§ Stability, continuity and interpolation results

In this subsection we derive some important auxiliary results which are later used in the error analysis of the semi-discrete and the fully-discrete method. The first aim is to establish inf-sup stability of the bilinear forms $B_h$ and $B$ on the respective spaces. First we have to verify that the expression $\tnorm{(\cdot,\cdot)}_h$, respectively $\tnorm{(\cdot,\cdot)}$, are indeed norms. This is established in the following lemma and its corollary.

* For $\underline{\mathbf{U}}_h \in \ProdSemiDiscrSpace{k} + \ProdFullyDiscrSpace{k}{q}$
and $y \in H^1_0(\STdom)$ it holds that:

\begin{align}
& \int\limits_{\STdom} \left\{ -(\partial_t \underline{u}_1)  \partial_t y + \nabla \underline{u}_1 \nabla y \right\}
=  \sum\limits_{n=0}^{N-1} \int\limits_{I_n} \sum\limits_{K \in \mathcal{T}_h }  ( \partial_t \underline{u}_2 - \Delta \underline{u}_1 , y )_{K} \; \dT \label{eq:IP-H-1} \\
& + \sum\limits_{n=0}^{N-1} ( \underline{u}_2 - \partial_t u_1,  \partial_t y )_{Q^n}  + \sum\limits_{n=0}^{N-1} \int\limits_{I_n} \sum\limits_{F \in \mathcal{F}_i}  ( \jump{ \nabla \underline{u}_1 \cdot \mathbf{n} } , y  )_{F} \; \dT
+ \sum\limits_{n=1}^{N-1}(  \jump{ \underline{u}_2^n } , y )_{ \Omega}, \nonumber \end{align}

where $ \jump{ \nabla \underline{u}_1 \cdot \mathbf{n} } $ denotes the jump of the normal derivative over interior facets.

* If for $(\underline{\mathbf{U}}_h,\underline{\mathbf{Z}}_h) \in (\ProdSemiDiscrSpace{k} + \ProdFullyDiscrSpace{k}{q}) \times  (\ProdSemiDiscrSpace{k_{\ast} } + \ProdFullyDiscrSpace{ k_{\ast} }{ q_{\ast} }) $ it holds $\tnorm{ (\mathbf{U}_h, \mathbf{Z}_h) }_h = 0$, then

\begin{align}
& \underline{\mathbf{Z}}_h = 0, \quad \underline{u}_1|_{\Sigma} = 0, \quad \underline{u}_1|_{ \STdata } = 0, \quad \partial_t \underline{u}_1 = \underline{u}_2 \label{eq:h-norm-zero-impl} \\& \text{ and } \norm{\Box \underline{u}_1}_{ H^{-1}(\STdom)} = \sup_{\substack{  y \in H^1_0(\STdom), \\ \norm{y}_{H^1(\STdom) } = 1  }} \sum\limits_{n=1}^{N-1}(  \jump{ \underline{u}_2^n } , y )_{ \Omega}.
\label{eq:h-norm-zero-res-H-1}
\end{align}

* This follows immediately from integration by parts.

* In view of definitions <ref>-<ref>, the claims in <ref> follow directly from $\tnorm{ (\underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h) }_h = 0$. The claim in <ref> is a consequence of <ref>, since by definition

\[
\norm{\Box \underline{u}_1}_{ H^{-1}(\STdom)} = \sup_{\substack{  y \in H^1_0(\STdom), \\ \norm{y}_{H^1(\STdom) } = 1  }}  \int\limits_{\STdom} \left\{ -(\partial_t \underline{u}_1)  \partial_t y + \nabla \underline{u}_1 \nabla y \right\}
\]

and all terms on the right hand side of <ref>, except for the jump over the time slice boundaries, vanish due to $\abs{ \underline{\mathbf{U}}_h }_{S_h} = 0$.

* It holds that $\tnorm{ (\cdot, \cdot) }_h$ is a norm on  $ \ProdSemiDiscrSpace{k} \times \ProdSemiDiscrSpace{k_{\ast}} $.

* The stronger expression $\tnorm{ (\cdot, \cdot) }$ is a norm on  $\ProdFullyDiscrSpace{k}{q} \times \ProdFullyDiscrSpace{ k_{\ast} }{ q_{\ast} }$.

* Let  $\tnorm{ (\underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h) }_h = 0$. Since $\underline{u}_2 \in H^1(\STdom)$, we have $\Box \underline{u}_1 = 0 $ in $H^{-1}(\STdom)$ from <ref>. As $\underline{u}_1 \in  H^1(\STdom)$, we are allowed to apply the stability estimate from <ref> to $\phi = \underline{u}_1$. Thanks to <ref> we then obtain as required that

\[ \underline{u}_1 = 0, \quad \underline{u}_2 = \partial_t \underline{u}_1 = 0, \quad \text{and }  \underline{\mathbf{Z}}_h = 0.
\]

* Since $ \abs{\underline{\mathbf{U}}_h }_{ \uparrow \downarrow } = 0$ implies $(\underline{u}_1, \underline{u}_2) \in [H^1(\STdom)]^2$, the same argument as in (a) applies.

Due to

\[
B_h[ (\underline{\mathbf{U}}_h,\underline{\mathbf{Z}}_h), ( \underline{\mathbf{U}}_h ,-\underline{\mathbf{Z}}_h ) ]  = \tnorm{ (\underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h) }_h \tnorm{ (\underline{\mathbf{U}}_h, -\underline{\mathbf{Z}}_h) }_h,
\quad
S^{\uparrow \downarrow}_{\Delta t}( \underline{\mathbf{U}}_h, \underline{\mathbf{U}}_h ) = \abs{ \underline{\mathbf{U}}_h   }_{ \uparrow \downarrow }^2
\]

and <ref> we then obtain inf-sup stability on the respective spaces:

\begin{align}
\label{eq:inf-sup}
\sup_{ (\underline{\mathbf{W}}_h,\underline{\mathbf{Y}}_h) \in  \ProdSemiDiscrSpace{k} \times \ProdSemiDiscrSpace{ k_{\ast} }    } \!\!\!\!\!\!\!\!\!\!\!\! \frac{ B_h[ (\underline{\mathbf{U}}_h,\underline{\mathbf{Z}}_h), ( \underline{\mathbf{W}}_h , \underline{\mathbf{Y}}_h ) ] }{  \tnorm{ (\underline{\mathbf{W}}_h, \underline{\mathbf{Y}}_h) }_h } &\geq C \tnorm{ ( \underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h) }_h, \\
\sup_{ (\underline{\mathbf{W}}_h,\underline{\mathbf{Y}}_h) \in  \ProdFullyDiscrSpace{k}{q} \times \ProdFullyDiscrSpace{ k_{\ast} }{ q_{\ast} }    } \!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \frac{ B[ (\underline{\mathbf{U}}_h,\underline{\mathbf{Z}}_h), ( \underline{\mathbf{W}}_h , \underline{\mathbf{Y}}_h ) ] }{  \tnorm{ (\underline{\mathbf{W}}_h, \underline{\mathbf{Y}}_h) } } &\geq C \tnorm{ ( \underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h) }. \nonumber
\end{align}

In the next lemma we record some continuity results for the bilinear form $A[\cdot,\cdot]$.

* For $\underline{\mathbf{U}}_h \in \ProdSemiDiscrSpace{k} + \ProdFullyDiscrSpace{k}{q}$ and $ \mathbf{Y} \in [H^1(\STdom)]^2 + \ProdFullyDiscrSpace{q_{\ast} }{k_{\ast} } $ we have

\[ A[\underline{\mathbf{U}_h}, \mathbf{Y} ] \leq C \abs{ \underline{\mathbf{U}}_h }_{\underline{S}_h} \left\{ \sum\limits_{n=0}^{N-1} h^{-2} \norm{ y_1 }_{Q^n}^2 +  \norm{ \nabla y_1 }_{Q^n}^2 + \norm{ y_2 }_{Q^n}^2  \right\}^{1/2}. \]

* For $\mathbf{U}|_{Q^n} \in \left[ H^{1}(Q^n) \cap L^2(0,T;H^2(\mathcal{T}_h)) \right] \times H^{1}(Q^n) $ for all $n=0,\ldots,N-1$ and $\underline{\mathbf{Y}}_h \in \ProdSemiDiscrSpace{k_{\ast}} + \ProdFullyDiscrSpace{ k_{\ast } }{ q_{\ast} }$ it holds that

\begin{align*}
A[ \mathbf{U}, \underline{\mathbf{Y}}_h ] \leq C & \bigg( \sum\limits_{n=0}^{N-1} \big\{ \norm{\partial_t u_2}_{Q^n}^2  + \norm{ \nabla u_1 }_{Q^n}^2
+ \norm{\partial_t u_1}_{Q^n}^2 + \norm{ u_2}_{Q^n}^2  \\
& + \int\limits_{I_n} \sum\limits_{K \in \mathcal{T}_h } h^2 \norm{u_1}_{H^2(K)}^2 \; \dT \big\}  \bigg)^{1/2}  \norm{  \underline{\mathbf{Y}}_h }_{ S^{\ast}_h  }.
\end{align*}

Let us first prove (a). Performing integration by parts in space yields

\begin{align*}
& \sum\limits_{n=0}^{N-1} \left\{ (\partial_t \underline{u}_2, y_1)_{Q^n} + a( \underline{u}_1,y_1)_{Q^n} - ( \nabla \underline{u}_1 \cdot \mathbf{n}, y_1)_{ \Sigma^n }
\right\} \\
&= \sum\limits_{n=0}^{N-1} \left\{  \int\limits_{I_n} \sum\limits_{K \in \mathcal{T}_h } ( \partial_t \underline{u}_2 - \Delta \underline{u}_1, y_1 )_{K} \; \dT
+ \int\limits_{I_n} \sum\limits_{F \in \mathcal{F}_i}   ( \jump{ \nabla \underline{u}_1 \cdot \mathbf{n} } , y_1   )_{F} \; \dT
\right\} \\
&\leq C \left[  G( \underline{\mathbf{U}}_h, \underline{\mathbf{U}}_h )^{1/2}  \left( \sum\limits_{n=0}^{N-1} \frac{1}{h^2} \norm{ y_1 }_{Q^n}^2   \right)^{1/2} \!\!\!\!\!\!\!\!
+ J( \underline{\mathbf{U}}_h, \underline{\mathbf{U}}_h )^{1/2} \left(  \sum\limits_{n=0}^{N-1} \sum\limits_{F \in \mathcal{F}_i} \frac{1}{h} \norm{y_1}_{F}^2 \; \dT  \right)^{1/2}
\right] \\
&\leq C \abs{ \underline{\mathbf{U}}_h }_{ S_h } \left( \sum\limits_{n=0}^{N-1} h^{-2} \norm{ y_1 }_{Q^n}^2 +  \norm{ \nabla y_1 }_{Q^n}^2  \right)^{1/2},
\end{align*}

where the trace inequality <ref> has been used in the last step.

\[
\sum\limits_{n=0}^{N-1} ( \partial_t \underline{u}_1 - \underline{u}_2, y_2  )_{ Q^n }
\leq  \sum\limits_{n=0}^{N-1} \norm{ \partial_t \underline{u}_1 - \underline{u}_2 }_{ Q^n } \norm{ y_2 }_{ Q^n }
\leq C \abs{ \underline{\mathbf{U}}_h }_{ S_h } \left( \sum\limits_{n=0}^{N-1} \norm{ y_2 }_{ Q^n }^2 \right)^{1/2}.
\]

This shows (a). Now to establish (b) we note that

\begin{align*}
\sum\limits_{n=0}^{N-1} \left\{ (\partial_t u_2, \underline{y}_1)_{Q_n} + a( u_1,\underline{y}_1)_{Q^n} \right\}  \leq C \left( \sum\limits_{n=0}^{N-1} \norm{\partial_t u_2}_{Q^n}^2  + \norm{ \nabla u_1 }_{Q^n}^2  \right)^{1/2} \norm{  \underline{\mathbf{Y}}_h }_{ S^{\ast}_h  }.
\end{align*}

Further, by using the trace inequality <ref> we obtain

\begin{align*}
\sum\limits_{n=0}^{N-1} ( \nabla \underline{u}_1 \!\! \cdot \mathbf{n}, \underline{y}_1)_{ \Sigma^n }
& \leq \left( \sum\limits_{n=0}^{N-1} \int\limits_{I_n} \! \sum\limits_{K \in \mathcal{T}_h^n }\!\!\! h \norm{ \nabla u_1 }_{ \partial \Omega \cap \partial K }^2 \! \right)^{\frac{1}{2} } \!\!\!
\left(  \sum\limits_{n=0}^{N-1} \int\limits_{I_n} \! \sum\limits_{ K \in \mathcal{T}_h^n } \frac{1}{h} \norm{ \underline{y}_1}_{ \partial \Omega \cap \partial K }^2 \right)^{\frac{1}{2} } \\
& \leq C \left(  \sum\limits_{n=0}^{N-1} \big\{  \norm{ \nabla u_1 }_{Q^n}^2  + \int\limits_{I_n} \sum\limits_{ K \in \mathcal{T}_h } h^2 \norm{u_1}_{H^2(K)}^2 \; \dT \big\}  \right)^{1/2}  \norm{  \underline{\mathbf{Y}}_h }_{S^{\ast}_h }.
\end{align*}

An application of the Cauchy-Schwarz inequality then covers the remaining term

\[
\sum\limits_{n=0}^{N-1} ( \partial_t u_1 - u_2, \underline{y}_2  )_{ Q^n }
\leq C \left( \sum\limits_{n=0}^{N-1} \norm{ \partial_t u_1 }_{ Q^n }^2 + \norm{ u_2 }_{ Q^n }^2 \right)^{1/2} \left( \sum\limits_{n=0}^{N-1} \norm{ \underline{y}_2 }_{ Q^n }^2 \right)^{1/2}.
\]

We also require some interpolation operators into the semi and fully-discrete spaces.
For the proof of the following lemma we refer to the appendix.

There exist interpolation operators $\Pi_h^k$ into $\SemiDiscrSpace{k}$ and $\Pi_{h, \Delta t}^{k,q}$ into $\FullyDiscrSpace{k}{q}$ such that
for $\Pi_h \in \{ \Pi_h^k,  \Pi_{h, \Delta t}^{k,q}\}$ under the assumption $\Delta t = C h$ the following interpolation estimates are valid for

* $\sum\limits_{n=0}^{N-1} \{  h^{-1} \norm{u - \Pi_h u }_{ L^2(Q^n) } + \norm{ \nabla (u - \Pi_h u) }_{ L^2(Q^n) } \} \leq C \norm{u}_{H^1(Q)}  $.

* $\sum\limits_{n=0}^{N-1}  \{ h^{-1} \norm{u - \Pi_h u }_{ L^2(Q^n) } + \norm{ u - \Pi_h u }_{ H^1(Q^n) } \} \leq C h^s \norm{u}_{H^{m-1}(Q)}$,

* For $u \in H^m(Q^n) \cap C^0(I_n, H^2(\Omega))$ we have:

\[ \left(\sum\limits_{n=0}^{N-1} \int\limits_{I_n} \sum\limits_{K \in \mathcal{T}_h } h^2 \norm{ u - \Pi_h  u }_{ H^2(K) }^2 \; \dT \right)^{1/2} \leq C h^s \norm{u}_{H^{m}(Q)}, \]

where we have defined

\begin{equation}\label{def:s-m}
(s,m) =
\begin{cases}
(k,k+2), &  \text{ when } \Pi_h = \Pi_h^k, \\
( \min\{k,q\}, \max\{k,q\}+3) , & \text{ when } \Pi_h = \Pi_{h, \Delta t}^{k,q}.
\end{cases}
\end{equation}

In the following we continue to use the generic notation $\Pi_h$ to prove some results which are valid for both interpolation operators. Interpolation into the spaces
containing the dual variable is denoted by $\Pi_h^{\ast} \in \{ \Pi_h^{k_{\ast}},  \Pi_{h, \Delta t}^{k_{\ast},q_{\ast} }\}$.
These operators are extended to tuples $\mathbf{W} = (w_1,w_2)$ in a componentwise-manner, i.e.$\mathbf{\Pi}_h \mathbf{W} :=  ( \Pi_h w_1 , \Pi_h w_2 )$.
We also keep the assumption $\Delta t = C h$ from now on.

Let $u \in H^{m}(Q)$ solve <ref> and set $ \mathbf{U} = (u, \partial_t u)$. Let $(s,m)$ be as defined as in <ref>.

* For  $u \in H^{m }(Q) \cap C^{0}([0,T],H^2(\Omega))$ we have

\[ \abs{ \mathbf{\Pi}_h \mathbf{U} }_{ S_h} = \abs{ \mathbf{\Pi}_h \mathbf{U} - \mathbf{U} }_{ S_h}   \leq C h^{ s }  \norm{u}_{ H^{m}(\STdom)  }.  \]

* For  $ \underline{\mathbf{U}}_h \in \ProdSemiDiscrSpace{k} + \ProdFullyDiscrSpace{k}{q}$ it holds that

\[
\norm{u-\underline{u}_1}_{\Sigma}  + \norm{u-\underline{u}_1}_{\STdata}
\leq C \left(  h^{ s+1/2 }  \norm{u}_{ H^{m}(\STdom)} + \tnorm{ ( \underline{\mathbf{U}}_h - \mathbf{\Pi}_h \mathbf{U},0) }_h \right).
\]

* For $\mathbf{Y} = (y_1,y_2) \in [H_0^1(\STdom )]^2$ we have

\[ \norm{ \mathbf{\Pi}_h^{\ast}  \mathbf{Y} }_{ S^{\ast}_h } \leq C \left( \norm{y_1}_{H^1(\STdom)} + \norm{y_2}_{H^1(\STdom)} \right).
\]

* We estimate the terms in $\abs{ \mathbf{\Pi}_h \mathbf{U} }_{ S_h}$ separately. By using smoothness of $u$, the trace inequality <ref> and the interpolation estimates from <ref> we obtain

\begin{align*}
& J( \mathbf{\Pi}_h \mathbf{U} , \mathbf{\Pi}_h \mathbf{U} ) = \sum\limits_{n=0}^{N-1} \int\limits_{I_n} \sum\limits_{F \in \mathcal{F}_i} h \norm{\jump{ \nabla ( \Pi_h u  - u)}}_{F}^2 \; \dT \\
&\leq C \sum\limits_{n=0}^{N-1} \int\limits_{I_n} \sum\limits_{ K \in \mathcal{T}_h }  \left(   \norm{  \nabla ( \Pi_h u - u ) }^2_{ L^2(K)  } + h^2 \norm{  \Pi_h u - u  }^2_{ H^2(K)  }  \right) \dT \\
& \leq C h^{2 s }  \norm{u}_{ H^m(\STdom)  }^2.
\end{align*}

Inserting $\partial_t^2 u - \Delta u = 0$ and applying interpolation bounds gives

\begin{align*}
G( \mathbf{\Pi}_h \mathbf{U} ,\mathbf{\Pi}_h \mathbf{U} ) &= \sum\limits_{n=0}^{N-1} \int\limits_{I_n} \sum\limits_{K \in \mathcal{T}_h^n } h^2 \norm{ \partial_t \left( \Pi_h \partial_t u - \partial_t u \right) + \Delta ( u - \Pi_h  u) }_{K}^2  \; \dT \\
&  \leq C h^{2s} \left( \norm{\partial_t u }_{H^{m-1}(\STdom) }^2  + \norm{u}_{H^m(\STdom)}^2  \right) \leq C h^{2s} \norm{u}_{H^m(\STdom)}^2.
\end{align*}

Similarly, we estimate

\begin{align*}
I_0( \mathbf{\Pi}_h  \mathbf{U} , \mathbf{\Pi}_h  \mathbf{U} ) & =  \sum\limits_{n=0}^{N-1} \norm{  ( \Pi_h  \partial_t u - \partial_t u) +  \partial_t(u - \Pi_h  u)    }^2_{ Q^n } \leq C h^{2s} \norm{u}_{H^m(\STdom)}^2.
\end{align*}

The boundary term is treated by using the trace inequality and then appealing to the interpolation results:

\begin{align*}
&  R( \mathbf{\Pi}_h \mathbf{U} , \mathbf{\Pi}_h \mathbf{U} )  =  \sum\limits_{n=0}^{N-1} \int\limits_{I_n} \sum\limits_{K \in \mathcal{T}_h } h^{-1} \norm{ \Pi_h u - u  }_{ L^2(\partial \Omega \cap \partial K )  }^2    \dT \\
\leq &  C  \sum\limits_{n=0}^{N-1} \int\limits_{I_n} \sum\limits_{K \in \mathcal{T}_h } \left( h^{-2}  \norm{ \Pi_h u - u}_{ L^2(K ) }^2 + \norm{ \nabla \left( \Pi_h  u - u \right) }_{ L^2(K ) }^2    \right)  \dT \\
& \leq C h^{2s} \norm{u}_{H^{m-1} (\STdom)}^2.
\end{align*}

* From interpolation estimates and the definition <ref> of $\tnorm{\cdot}_h$ we obtain

\begin{align*}
\norm{u-\underline{u}_1}_{\STdata} & \leq C \left( \norm{u-\Pi_h u }_{\STdata} + \tnorm{ ( \mathbf{\Pi}_h \mathbf{U} - \underline{ \mathbf{U}}_h , 0) }_h  \right)  \\
& \leq C \left(  h^{ s+1 }  \norm{u}_{ H^{m}(\STdom)} + \tnorm{ ( \underline{  \mathbf{U}}_h  - \mathbf{\Pi}_h \mathbf{U}  , 0) }_h  \right).
\end{align*}

We proceed similarly for the boundary term using the trace inequality:

\begin{align*}
& \norm{u-\underline{u}_1}_{\Sigma}^2
\leq C  \sum\limits_{n=0}^{N-1} \int\limits_{I_n} \sum\limits_{K \in \mathcal{T}_h }  \norm{ u - \Pi_h u }_{ L^2(\partial \Omega \cap \partial K )  }^2 +  \frac{1}{h} \norm{ \Pi_h u - \underline{u}_1  }_{ L^2(\partial \Omega \cap \partial K )  }^2   \\
& \leq  C  \left(  \sum\limits_{n=0}^{N-1} \left[ \frac{1}{h} \norm{  u - \Pi_h u }^2_{Q^n} +  h \norm{ \nabla (  u - \Pi_h u ) }^2_{Q^n}   \right] +  \tnorm{ ( \underline{  \mathbf{U}}_h  - \mathbf{\Pi}_h \mathbf{U}  , 0) }_h^2  \right) \\
& \leq C \left(  h^{ 2(s+1/2) }  \norm{u}_{ H^{m}(\STdom)}^2  + \tnorm{ ( \underline{  \mathbf{U}}_h  - \mathbf{\Pi}_h \mathbf{U}  , 0) }_h^2 \right).
\end{align*}

* By applying <ref> (a) we obtain

\[
\norm{ \Pi_h^{\ast} y_1  }_{ Q^n } +  \norm{ \Pi_h^{\ast} y_2  }_{ Q^n } + \norm{ \nabla \Pi_h^{\ast}  y_1  }_{ Q^n }
\leq C  \left( \norm{y_1}_{H^1(Q^n)} + \norm{y_2}_{H^1(Q^n)}  \right).
\]

Using that $y_1$, vanishes on the lateral boundary we obtain by means of the trace inequality and <ref> (a):

\begin{align*}
& h^{-1} \norm{ \Pi_h^{\ast} y_1}_{\Sigma^n}^2 = h^{-1} \norm{ \Pi_h^{\ast} y_1 - y_1}_{\Sigma^n}^2 \\
& \leq C \int\limits_{I_n} \left( h^{-2} \norm{ \Pi_h^{\ast} y_1 - y_1  }_{ L^2(\Omega ) }^2 + \norm{ \Pi_h^{\ast} y_1 - y_1  }_{ H^1(\Omega ) }^2  \right)  \; \dT \leq C \norm{y_1}_{H^1(Q^n)}^2.
\end{align*}

§.§ Consistency and residual bounds

In this subsection we establish two additional results which are crucial for the analysis of the semi-discrete and the fully-discrete method. The first lemma mainly follows from consistency of the stabilization $S_h$ and the bilinear form $A$.

Let $u$ be a solution of <ref> and set $ \mathbf{U} = (u, \partial_t u)$. Then for any
$ \underline{\mathbf{W}}_h \in \ProdSemiDiscrSpace{k} + \ProdFullyDiscrSpace{k}{q}$ and $ \underline{\mathbf{Y}}_h \in \ProdSemiDiscrSpace{k_{\ast}} + \ProdFullyDiscrSpace{ k_{\ast } }{ q_{\ast} }$ it holds that

\begin{align*}
(u -  \Pi_h u,\underline{w}_1)_{\STdata} -  S_h(\mathbf{\Pi}_h \mathbf{U}, \underline{\mathbf{W}}_h) -
A[ \mathbf{\Pi}_h \mathbf{U} , \underline{\mathbf{Y}}_h ]
\leq C h^{s} \norm{u}_{H^m(\STdom)} \tnorm{ (\underline{\mathbf{W}}_h, \underline{\mathbf{Y}}_h) }_h,
\end{align*}

where $s$ is as defined in <ref>.

From <ref> (a) and interpolation results we immediately have

\begin{align*}
(u -  \Pi_h u,\underline{w}_1)_{\STdata} - S_h(\mathbf{\Pi}_h \mathbf{U}, \underline{\mathbf{W}}_h)
& \leq C \left( \norm{ u -  \Pi_h u}_{ \STdata} + \abs{ \mathbf{\Pi}_h \mathbf{U}}_{S_h}  \right)  \tnorm{ (\underline{\mathbf{W}}_h, 0) }_h \\
& \leq C h^{s} \norm{u}_{H^m(\STdom)}  \tnorm{ (\underline{\mathbf{W}}_h, 0) }_h.
\end{align*}

For the remaining term we use $ A[ \mathbf{U} , \underline{\mathbf{Y}}_h ] = 0 $, <ref> (b) and interpolation estimates from
<ref> (b)-(c) to obtain

\begin{align*}
& A[ \mathbf{\Pi}_h \mathbf{U} , \underline{\mathbf{Y}}_h ] = A[ \mathbf{\Pi}_h \mathbf{U} - \mathbf{U}, \underline{\mathbf{Y}}_h ]  \\
& \leq C \left( \sum\limits_{n=0}^{N-1}  \sum\limits_{i=1,2} \left\{ \norm{ \Pi_h u_i - u_i }_{ H^1(Q^n) }^2 \right\} + \int\limits_{I_n} \sum\limits_{K \in \mathcal{T}_h } h^2 \norm{  \Pi_h  u_1 - u_1  }_{ H^2(K) }^2   \right)^{1/2}
\norm{  \underline{\mathbf{Y}}_h }_{ S^{\ast}_h  } \\ %\leq C h^{s} \norm{u}_{H^m(\STdom)}  \norm{  \underline{\mathbf{Y}}_h }_{ S^{\ast}_h  }.
& \leq C h^{s} \norm{u}_{H^m(\STdom)}  \norm{  \underline{\mathbf{Y}}_h }_{ S^{\ast}_h  }.
\end{align*}

The second lemma will later allow us to control the PDE residual $\norm{\Box \underline{u}_1 }_{H^1(\STdom)}$.

Let $ ( \underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h, \mathbf{\Pi}_h^{\ast} ) \in   \{( \ProdSemiDiscrSpace{k},\ProdSemiDiscrSpace{k_{\ast}},\mathbf{\Pi}_h^{k_{\ast}}) ,(\ProdFullyDiscrSpace{k}{q}, \ProdFullyDiscrSpace{ k_{\ast } }{ q_{\ast} }, \mathbf{\Pi}_{h, \Delta t}^{k_{\ast},q_{\ast}}) \}$ such that

\begin{equation}\label{eq:optimality_dual}
A[ \underline{\mathbf{U}}_h, \mathbf{\Pi}_h^{\ast} \mathbf{W} ] = S_h^{\ast}(  \mathbf{\Pi}_h^{\ast} \mathbf{W}  , \underline{\mathbf{Z}}_h)
\end{equation}

for any $ \mathbf{W} = (w_1,w_2) = (w,0)$ with $w \in H^1_0(\STdom)$ holds true. Then for  $ \mathbf{U} = (u, \partial_t u)$ with $u$ being the solution of <ref> we have:

\begin{align*}
& \sum\limits_{n=0}^{N-1} \{ (\partial_t \underline{u}_1, \partial_t w_1)_{ Q^n } - a(\underline{u}_1,w_1)_{Q^n} \} \\
& \leq \left\vert \sum\limits_{n=1}^{N-1} ( \jump{\underline{u}_2^n}, w_{1,+}^n  )_{\Omega} \right\vert + C \norm{w}_{H^1(\STdom)} \left(  \tnorm{ ( \underline{\mathbf{U}}_h - \mathbf{\Pi}_h \mathbf{U} ,\underline{\mathbf{Z}}_h )   }_h  +  \abs{ \mathbf{\Pi}_h \mathbf{U} }_{ S_h}  \right).
%&  + \leq C \norm{w}_{H^1(\STdom)} \left(  \tnorm{ \underline{\mathbf{U}}_h - \mathbf{\Pi}_h \mathbf{U}   }_h +  \abs{ \mathbf{\Pi}_h \mathbf{U} }_{ S_h}  \right),
\end{align*}

Integration by parts in time yields

\begin{align*}
& \sum\limits_{n=0}^{N-1} \{ (\partial_t \underline{u}_1, \partial_t w_1)_{ Q^n } - a(\underline{u}_1,w_1)_{Q^n} \}
= \sum\limits_{n=0}^{N-1} (\partial_t \underline{u}_1 - \underline{u}_2, \partial_t w_1)_{ Q^n } -  \sum\limits_{n=1}^{N-1} ( \jump{\underline{u}_2^n}, w_{1,+}^n  )_{\Omega} \\
& - \sum\limits_{n=0}^{N-1} \left\{  (\partial_t \underline{u}_2 , w_1)_{ Q^n } +  a(\underline{u}_1,w_1)_{Q^n} +  (\partial_t \underline{u}_1 - \underline{u}_2 , \underbrace{ w_2 }_{=0} )_{ Q^n } - (\nabla \underline{u}_1 \cdot \mathbf{n},  \underbrace{w_1}_{ = 0  } )_{ \Sigma^n } \right\}  \\
& = \sum\limits_{n=0}^{N-1} (\partial_t \underline{u}_1 - \underline{u}_2, \partial_t w_1)_{ Q^n } -  \sum\limits_{n=1}^{N-1} ( \jump{\underline{u}_2^n}, w_{1,+}^n  )_{\Omega}
- A[\underline{\mathbf{U}_h},  \mathbf{W} ].
\end{align*}

We treat the terms separately. Clearly,

\[
\sum\limits_{n=0}^{N-1} (\partial_t \underline{u}_1 - \underline{u}_2, \partial_t w_1)_{ Q^n }
\leq \underline{I}_0( \underline{\mathbf{U}}_h, \underline{\mathbf{U}}_h )^{1/2} \norm{\partial_t w_1  }_{\STdom}
\leq C \abs{  \underline{\mathbf{U}}_h   }_{ S_h } \norm{w}_{H^1(\STdom)}.
\]

To cope with the remaining term we insert $ \mathbf{\Pi}_h^{\ast}  \mathbf{W} = (  \Pi_h^{\ast} w,0) $ so that
A[𝐔_h, 𝐖 ] = A[𝐔_h, 𝐖 -  Π_h^∗ 𝐖 ] + A[ 𝐔_h, Π_h^∗ 𝐖].
From <ref> (a) and <ref> (a)  we obtain

\begin{align*}
\underline{A}[\underline{\mathbf{U}}_h, \mathbf{W} - \mathbf{\Pi}_h^{\ast} \mathbf{W} ]
& \leq C \abs{ \underline{\mathbf{U}}_h }_{ S_h } \bigg\{ \sum\limits_{n=0}^{N-1} \frac{1}{h^2} \norm{ w_1 - \Pi_h^{\ast} w_1 }_{Q^n}^2 + \norm{ \nabla (w_1 - \Pi_h^{\ast} w_1 )}_{Q^n}^2  \bigg\}^{1/2} \\
& \leq C  \abs{ \underline{\mathbf{U}}_h }_{ S_h } \norm{w}_{H^1(\STdom)}.
\end{align*}

For the other contribution we obtain from <ref> and <ref> (c) that

\[
A[\underline{\mathbf{U}}_h,  \mathbf{\Pi}_h^{\ast} \mathbf{W} ] = \underline{S}^{\ast} ( \mathbf{\Pi}_h^{\ast} \mathbf{W}, \underline{\mathbf{Z}}_h)
\leq C \norm{ \mathbf{\Pi}_h^{\ast} \mathbf{W}}_{ S^{\ast}_h }  \norm{ \underline{\mathbf{Z}}_h }_{ S^{\ast}_h }
\leq C  \norm{ \underline{\mathbf{Z}}_h }_{ S^{\ast}_h } \norm{w}_{H^1(\STdom)}.
\]

The claim follows in view of

\[
\abs{ \underline{\mathbf{U}}_h }_{ S_h } + \norm{ \underline{\mathbf{Z}}_h }_{ S^{\ast}_h }
\leq \tnorm{ ( \underline{\mathbf{U}}_h - \mathbf{\Pi}_h \mathbf{U} ,\underline{\mathbf{Z}}_h )   }_h + \abs{ \mathbf{\Pi}_h \mathbf{U} }_{ S_h}.
\]

§ ERROR ANALYSIS FOR THE SEMI-DISCRETE METHOD

In this section we derive an error estimate for the semi-discrete method defined in <ref>. Thanks to the auxiliary results derived in <ref>-
<ref> and the fact that no jumps over time slice boundaries occur, we easily obtain the desired result. First we deduce convergence of the discretization error in the triple norm $\tnorm{ (\cdot , \cdot ) }_h$.

Let $u$ be a sufficiently regular solution of <ref> and set $\mathbf{U} := (u,\partial_t u) $.
Let  $ (\underline{\mathbf{U}}_h ,\underline{\mathbf{Z}}_h) \in  \ProdSemiDiscrSpace{k} \times \ProdSemiDiscrSpace{k_{\ast}}  $
be the solution of <ref>. Then

\[
\tnorm{ (\underline{\mathbf{U}}_h - \mathbf{\Pi}_h \mathbf{U} , \underline{\mathbf{Z}}_h ) }_h \leq C h^{k}  \norm{u}_{ H^{k+2}(\STdom)  }.
\]

From <ref> we have

\[
B_h[ ( \underline{\mathbf{U}}_h - \mathbf{\Pi}_h \mathbf{U}, \underline{\mathbf{Z}}_h ), ( \underline{\mathbf{W}}_h ,\underline{\mathbf{Y}}_h ) ]
= (u -  \Pi_h u,\underline{w}_1)_{\STdata} -  S_h(\mathbf{\Pi}_h \mathbf{U}, \underline{\mathbf{W}}_h) -
A[ \mathbf{\Pi}_h \mathbf{U} , \underline{\mathbf{Y}}_h]
\]

for any $ (\underline{\mathbf{W}}_h ,\underline{\mathbf{Y}}_h) \in  \ProdSemiDiscrSpace{k} \times \ProdSemiDiscrSpace{k_{\ast}} $
Hence, the claim follows from <ref> and the inf-sup condition <ref> for $B_h$.

We then improve the former result to obtain convergence in a physically significant norm by utilizing the Lipschitz stability result from <ref>.

Let $u$ be a sufficiently regular solution of <ref> and $ (\underline{\mathbf{U}}_h ,\underline{\mathbf{Z}}_h) \in  \ProdSemiDiscrSpace{k} \times \ProdSemiDiscrSpace{k_{\ast}}  $ be the solution of <ref>. Then

\[
\norm{ u - \underline{u}_1 }_{L^{\infty}(0,T;L^2(\Omega))} + \norm{\partial_t (u-\underline{u}_1) }_{L^{2}(0,T;H^{-1}(\Omega))}
\leq C h^{k} \norm{u}_{ H^{k+2}(\STdom)  }.
\]

Since $e_h := u - \underline{u}_1 \in H^1(\STdom)$, we can directly apply the stability estimate <ref> to $\phi = e_h $ to obtain

\begin{align*}
\norm{ e_h }_{L^{\infty}(0,T;L^2(\Omega))} + \norm{\partial_t e_h }_{L^{2}(0,T;H^{-1}(\Omega))} \!
\leq C \! \left(  \norm{ e_h }_{ \STdata } \! + \! \norm{ e_h }_{ L^2(\Sigma)} \! + \! \norm{ \Box e_h }_{ H^{-1} ( \STdom ) } \right) \!.
\end{align*}

From <ref> (b) and <ref> we see that the first two terms on the right hand side are bounded
by a constant times $ h^{k}  \norm{u}_{ H^{k+2}(\STdom)  }$. The same bound is valid for the last term, which we deduce by writing down the definition

\begin{align*}
\norm{ \Box e_h }_{ H^{-1} ( \STdom ) } = \norm{ - \Box \underline{u}_1}_{ H^{-1} ( \STdom ) } =  \sup_{\substack{w \in H^1_0(\STdom), \\ \norm{w}_{H^1(\STdom) } = 1 }} \sum\limits_{n=0}^{N-1} \{ (\partial_t \underline{u}_1, \partial_t w_1)_{ Q^n } - a(\underline{u}_1,w_1)_{Q^n} \}
\end{align*}

and applying <ref>[Note that <ref> holds because of the first order optimality condition <ref>.], <ref>, <ref> (a) and the fact that $ \jump{\underline{u}_2^n} = 0$ for $\underline{u}_2 \in \SemiDiscrSpace{k}$ holds.

§ ERROR ANALYSIS FOR THE FULLY-DISCRETE METHOD

In this section we derive an error estimate for the fully-discrete method defined in <ref>. In view of the results that have already been established in <ref> for the semi-discrete method, this basically amounts to treating the additional coupling terms that were introduced by adding the term $S^{\uparrow \downarrow}_{\Delta t}$ on top of the semi-discrete Lagrangian.
As indicated in <ref> we will set $(s,m) = ( \min\{k,q\}, \max\{k,q\}+3)$ and hold on to the assumption $\Delta t = C h$ in this section. First we require some additional interpolation results at the boundary of the time slices whose proof can be found in the appendix.

At the top and bottom of the time slice the following approximation bounds for the traces holds true for $n=0,\ldots,N-1$:

* $ \norm{(u - \Pi_h^{\ast} u )( \tau, \cdot )}_{\Omega} \leq C h^{1/2} \norm{u}_{H^{1}(Q^n)},  \; \tau \in \{ t_n,t_{n+1}\}$,

* $ \norm{(u - \Pi_h u )( \tau, \cdot )}_{\Omega} \leq C h^{s+1/2} \norm{u}_{H^{m-2}(Q^n)},  \; \tau \in \{ t_n,t_{n+1}\}$,

* $  \norm{\nabla(u - \Pi_h u )( \tau, \cdot )}_{\Omega} \leq C h^{s-1/2} \norm{u}_{H^{m-1}(Q^n)},  \; \tau \in \{ t_n,t_{n+1}\}$.

We then deduce the counterpart of <ref> (a) for the additional stabilization term.

Let $u$ be a sufficiently regular solution of <ref> and set $\mathbf{U} := (u,\partial_t u) $. Then

\[
\abs{ \mathbf{\Pi}_h \mathbf{U} }_{ \uparrow \downarrow } \leq C  h^{s} \norm{u}_{H^{m-1}(\STdom)}.
\]

By using <ref> and smoothness of $u$ we obtain

\begin{align*}
& I_2^{\uparrow \downarrow}(  \mathbf{\Pi}_h \mathbf{U}, \mathbf{\Pi}_h \mathbf{U}) = \sum\limits_{n=1}^{N-1} \frac{1}{\Delta t} \norm{ \jump{ (\Pi_h \partial_t u - \partial_t u)^n  } }_{\Omega}^2 \\
& \leq \frac{2}{\Delta t } \left(  \sum\limits_{n=1}^{N-1} \norm{  (\Pi_h \partial_t u - \partial_t u)(t_n,\cdot)   }_{\Omega}^2 + \sum\limits_{n=0}^{N-2} \norm{  (\Pi_h \partial_t u - \partial_t u)(t_{n+1},\cdot)   }_{\Omega}^2    \right) \\
& \leq C h^{2s} \norm{u}_{H^{m-1}(\STdom)}^2.
\end{align*}

In a similar way we bound the contributions in $I_1$, e.g.
\begin{align*}
& \sum\limits_{n=1}^{N-1} \Delta t \norm{ \jump{ \nabla (\Pi_h u - u)^n  } }_{\Omega}^2  \leq 2 \Delta t \Big(  \sum\limits_{n=1}^{N-1} \norm{ \nabla (\Pi_h u - u)(t_n,\cdot)   }_{\Omega}^2 \\
& + \sum\limits_{n=0}^{N-2} \norm{ \nabla (\Pi_h u - u)(t_{n+1},\cdot)   }_{\Omega}^2  \Big)
\leq C h^{2s} \norm{u}_{H^{m-1}(\STdom)}^2.
\end{align*}

With this result, convergence of the discretization error in the strengthened triple norm $\tnorm{(\cdot,\cdot)}$ immediately follows.

Let $u$ be a sufficiently regular solution of <ref> and set $\mathbf{U} := (u,\partial_t u) $.
Let  $ (\underline{\mathbf{U}}_h ,\underline{\mathbf{Z}}_h) \in \ProdFullyDiscrSpace{k}{q} \times \ProdFullyDiscrSpace{ k_{\ast} }{ q_{\ast} } $
be the solution of <ref>. Then

\[
\tnorm{ (\underline{\mathbf{U}}_h - \mathbf{\Pi}_h \mathbf{U} , \underline{\mathbf{Z}}_h ) } \leq C h^{s}  \norm{u}_{ H^{m}(\STdom)  }.
\]

From <ref> we obtain

\begin{align*}
B[ ( \underline{\mathbf{U}}_h - \mathbf{\Pi}_h \mathbf{U}, \underline{\mathbf{Z}}_h ), ( \underline{\mathbf{W}}_h ,\underline{\mathbf{Y}}_h ) ]
= &  (u -  \Pi_h u,\underline{w}_1)_{\STdata} - S_h(\mathbf{\Pi}_h \mathbf{U}, \underline{\mathbf{W}}_h) -
A[ \mathbf{\Pi}_h \mathbf{U} , \underline{\mathbf{Y}}_h] \\
& \!\! - S^{\uparrow \downarrow}_{\Delta t}( \mathbf{\Pi}_h \mathbf{U} , \underline{\mathbf{W}}_h ).
\end{align*}

Since $S^{\uparrow \downarrow}_{\Delta t}( \mathbf{\Pi}_h \mathbf{U} , \underline{\mathbf{W}}_h ) \leq C h^{s} \norm{u}_{H^{m-1}(\STdom)} \abs{ \underline{\mathbf{W}}_h }_{ \uparrow \downarrow }  $ thanks to <ref>, the claim follows as in <ref>.

In the semi-discrete case we were able to apply the Lipschitz stability result <ref> directly to the error $u - \underline{u}_1$ to deduce our error estimate. In the fully-discrete case, this is no longer possible since $\underline{u}_1 \notin H^1(\STdom)$. To address this issue, we follow [6] and introduce a lifting operator $L_{\Delta t}:\FullyDiscrSpace{k}{q} \rightarrow C^0((0,T);V_h^k)$ as follows.
Let $\vartheta_{n}(t) := (t_{n+1} - t)/\Delta t$  for $ t \in I_n$ and set

\begin{equation}\label{eq:def_lifting}
L_{\Delta t} \underline{w}(t) = \underline{w}(t) - \jump{ \underline{w}^n }  \vartheta_{n}(t), \quad  t \in I_n, n \geq 1, \quad
L_{\Delta t} \underline{w}(t) = \underline{w}(t),  \quad t \in I_0.
\end{equation}

Note that $\vartheta_{n}(t_n) = 1$ and $ \vartheta_{n}(t_{n+1} ) =  0$,
which implies

\[
\lim_{ t \downarrow t_n} L_{\Delta t} \underline{w}(t)
= \underline{w}_{+}^{n} - \jump{ \underline{w}^n }
= \underline{w}_{-}^{n}
= \lim_{ t \uparrow t_n} L_{\Delta t} \underline{w}(t).
\]

Thus, $ L_{\Delta t} \underline{w}$ is continuous in time. For later purposes we also record that

\begin{equation}\label{eq:theta-and-deriv-L2-norm}
\norm{ \vartheta_{n}}_{ L^2(I_n) } = ( \Delta t)^{1/2}/\sqrt{3},
\qquad
\norm{ \vartheta_{n}^{\prime} }_{ L^2(I_n) } =  ( \Delta t)^{-1/2}.
\end{equation}

For the smoothed solution $L_{\Delta t} \underline{u}_1$ convergence rates follow similarly as in <ref>, where the additional jump terms induced by the lifting <ref> are covered by the stabilization term $S^{\uparrow \downarrow}_{\Delta t}  $. Note that the convergence order is independent of the polynomial order of the dual variable.

Let $u$ be a sufficiently regular solution of <ref> and $ (\underline{\mathbf{U}}_h ,\underline{\mathbf{Z}}_h) \in \ProdFullyDiscrSpace{k}{q} \times \ProdFullyDiscrSpace{ k_{\ast} }{ q_{\ast} } $ be the solution of <ref>.
Then for $(s,m) = ( \min\{k,q\}, \max\{k,q\}+3)$ it holds

\[
\norm{ u - L_{\Delta t} \underline{u}_1 }_{L^{\infty}(0,T;L^2(\Omega))} + \norm{\partial_t (u - L_{\Delta t} \underline{u}_1) }_{L^{2}(0,T;H^{-1}(\Omega))}
\leq C h^{s} \norm{u}_{ H^{m}(\STdom)  }.
\]

Since $L_{\Delta t} \underline{u}_1 \in H^1(Q)$, we are allowed to apply the continuum stability estimate <ref> to
$\tilde{e}_h := u - L_{\Delta t} \underline{u}_1$ which yields:

\begin{align}
& \norm{ \tilde{e}_h }_{L^{\infty}(0,T;L^2(\Omega))} + \norm{\partial_t \tilde{e}_h }_{L^{2}(0,T;H^{-1}(\Omega))} \nonumber \\
& \leq C \left(  \norm{ e_h }_{ \STdata } + \norm{ e_h }_{ \Sigma}
+  \norm{ \underline{u}_1 -  L_{\Delta t} \underline{u}_1 }_{ \STdata } + \norm{ \underline{u}_1 -  L_{\Delta t} \underline{u}_1 }_{ \Sigma }
+ \norm{ \Box \tilde{e}_h }_{ H^{-1} ( \STdom ) } \right), \label{eq:appl-stability-lifted-error}
\end{align}

where $e_h := u - \underline{u}_1$ as in the proof of <ref>.
We need to control the terms on the right hand side. Starting with $\norm{ \Box \tilde{e}_h }_{ H^{-1} ( \STdom ) }$, we have for any $w \in H_0^1(\STdom)$
in view of $\Box u = 0$:

\begin{align*}
(- \partial_t \tilde{e}_h, \partial_t w)_Q + a(\tilde{e}_h,w)_Q
& =   \sum\limits_{n=0}^{N-1} \{ (\partial_t \underline{u}_1, \partial_t w_1)_{ Q^n } - a(\underline{u}_1,w_1)_{Q^n} \} + \xi( \underline{u}_1, w) \\
& \leq C h^{s}  \norm{u}_{ H^{m}(\STdom)  } \norm{w}_{H^1(\STdom)} +  \xi( \underline{u}_1, w) + \eta( \underline{u}_2, w),
\end{align*}

where <ref>, <ref> and <ref> (a) have been used and

\begin{align*}
\xi( \underline{u}_1, w) :=  \sum\limits_{n=1}^{N-1} \{ a( \vartheta_{n} \jump{ \underline{u}_1^n } , w  )_{Q^n}  -( \vartheta^{\prime}_{n} \jump{ \underline{u}_1^n } , \partial_t w  )_{Q^n} \}, \;
\eta( \underline{u}_2, w) := \left\vert \sum\limits_{n=1}^{N-1} ( \jump{\underline{u}_2^n}, w_{+}^n  )_{\Omega} \right\vert.
\end{align*}

By using <ref> we obtain

\begin{align*}
& \xi( \underline{u}_1, w) \leq \left(  \sum\limits_{n=1}^{N-1} \norm{ \vartheta^{\prime}_{n}  }^2_{L^2(I_n)} \norm{ \jump{ \underline{u}_1^n } }_{ \Omega}^2
+ \norm{ \vartheta_{n}  }^2_{L^2(I_n)} \norm{ \jump{ \nabla \underline{u}_1^n } }_{ \Omega}^2  \right)^{1/2} \norm{ w }_{ H^1(\STdom) }  \\
& \leq C \left(  \sum\limits_{n=1}^{N-1}   (\Delta t)^{-1}  \norm{ \jump{ \underline{u}_1^n } }_{ \Omega}^2
+ \Delta t \norm{ \jump{ \nabla \underline{u}_1^n } }_{ \Omega}^2  \right)^{1/2} \norm{ w }_{ H^1(\STdom ) }
\leq C \abs{ \underline{\mathbf{U}}_h }_{ \uparrow \downarrow } \norm{ w }_{ H^1(\STdom) }.
\end{align*}

To treat $\eta( \underline{u}_2, w)$, we first note that

\[
\norm{w_{+}^n }_{\Omega} \leq \norm{ w_{+}^n - \Pi_h^{\ast} w_{+}^n  }_{\Omega} + \norm{  \Pi_h^{\ast} w_{+}^n  }_{\Omega}
\leq C \left( h^{1/2} + (\Delta t)^{-1/2}   \right) \norm{w}_{H^{1}(Q^n)} ,
\]

where we applied <ref> and a discrete inverse inequality in time (see <ref>) to estimate the second term.

\[
\eta( \underline{u}_2, w) \leq C \left(  \sum\limits_{n=1}^{N-1} \frac{1}{\Delta t} \norm{  \jump{\underline{u}_2^n} }_{\Omega}^2 \right)^{1/2}
\left(  \sum\limits_{n=1}^{N-1} \Delta t \norm{w_{+}^n}_{\Omega}^2 \right)^{1/2}
\leq C \abs{ \underline{\mathbf{U}}_h }_{ \uparrow \downarrow } \norm{ w }_{ H^1(\STdom) }.
\]

From <ref> and <ref> we obtain that $\abs{ \underline{\mathbf{U}}_h }_{ \uparrow \downarrow } \leq C h^{s}  \norm{u}_{ H^{m}(\STdom)  }$  and hence we conclude that the same bound holds for $\norm{ \Box \tilde{e}_h }_{ H^{-1} ( \STdom ) }$.
The terms $ \norm{ e_h }_{ L^2(\STdata) }$ and $ \norm{ e_h }_{ L^2(\Sigma)}$ obey this bound as well which follows similarly as in the proof
of <ref> for the semi-discrete case.
It remains to treat

\begin{align*}
& \norm{ \underline{u}_1 -  L_{\Delta t} \underline{u}_1 }_{ \STdata }^2 + \norm{ \underline{u}_1 -  L_{\Delta t} \underline{u}_1 }_{ \Sigma}^2
\leq C \Delta t \sum\limits_{n=1}^{N-1} \left( \norm{ \jump{ \underline{u}_1^n } }_{ \Omega}^2 + \sum\limits_{K \in \mathcal{T}_h } \norm{ \jump{ \underline{u}_1^n }  }_{ L^2(\partial \Omega \cap \partial K )  }^2   \right) \\
& \leq C \Delta t \sum\limits_{n=1}^{N-1} \! \sum\limits_{K \in \mathcal{T}_h } \left( \frac{1}{h} \norm{ \jump{ \underline{u}_1^n }  }_{K }^2 + \! h \norm{ \jump{ \nabla \underline{u}_1^n }  }_{K }^2  \right) \leq C h \abs{ \underline{\mathbf{U}}_h }_{ \uparrow \downarrow }^2 \leq C  h^{2s+1}  \norm{u}_{ H^{m}(\STdom)  }^2,
\end{align*}

where the trace inequality <ref> and $\Delta t = C h$  have been used.
Now we have bounded all terms on the right hand side of <ref> and the proof is complete.

§ PRECONDITIONING STRATEGIES

The fully-discrete method defined in <ref> is globally coupled in time. This is due to the presence of the term $S^{\uparrow \downarrow}_{\Delta t}( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h )$ defined in <ref>, which contains jump operators on both variables $\underline{\mathbf{U}}_h$ and $\underline{\mathbf{W}}_h$. Consequently, solving <ref> with a direct solver would require to assemble and factorize a $d+1$-dimensional linear system, which is prohibitively expensive for $d=3$. Hence, iterative solutions strategies must be investigated. In this setting Krylov methods such as GMRes are very attractive as they only require matrix-vector multiplications. Note here that application of the space-time matrix resulting from <ref> can be easily implemented time-slab wise which obliterates the need for assembling an enormous space-time system. Unfortunately, an unpreconditioned GMRes iteration is not suitable for solving our problem even in the simplest possible setting as <ref> demonstrates. The column `NoPre' in the table gives the GMRes iteration numbers for a simple problem on the unit interval discretized with piecewise linear elements, which are almost as high as the total number of degrees of freedom (`ndof' in the second column). Hence, a preconditioning strategy is required.  In this chapter we propose two simple strategies based on time-marching procedures.

§.§ Monolithic time-marching

A simple preconditioner can be obtained by relaxing the bi-directional coupling $\uparrow \downarrow$ in $S^{\uparrow \downarrow}$ to a forward coupling $\uparrow$ in time only. This is achieved by replacing $S^{\uparrow \downarrow}$ in <ref> by

\begin{equation}\label{eq:def-jump-stab-pre}
S^{\uparrow }_{\Delta t}( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h )
:= \underline{I}_1^{\uparrow}( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h ) +
\underline{I}_2^{\uparrow}( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h ),
\end{equation}

\begin{align*}
\underline{I}_1^{\uparrow}( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h ) :=& \sum\limits_{n=1}^{N-1} \left\{ \frac{1}{\Delta t} ( \jump{ \underline{u}_1^n } ,  \underline{w}_{1,+}^n  )_{ \Omega}
+ \Delta t ( \jump{ \nabla \underline{u}_1^n } , \nabla \underline{w}_{1,+}^n  )_{ \Omega} \right\}, \\
\underline{I}_2^{\uparrow}( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h ) :=& \sum\limits_{n=1}^{N-1} \frac{1}{\Delta t} ( \jump{ \underline{u}_2^n } ,  \underline{w}_{2,+}^n  )_{ \Omega}.
\end{align*}

This allows to solve for $(\underline{\mathbf{U}}_h,\underline{\mathbf{Z}}_h)$ time-slab by time-slab, i.e. starting with $Q^0$ and progressing one slab $Q^n$ at a time. Indeed, if the solution on $Q^{n-1}$ has been determined then in particular $\underline{u}_{j,-}^n$ for $j=1,2$ are known and can be shifted to the right hand side of the linear system as is common practice in the dG time-stepping method, see e.g. the monograph of Thomée [25] for details. This means that the cost of one application of the preconditioner is the same as the solution of a wave equation with a dG time-stepping method.

Since the dual stabilization and the slab-local part $S_h(\cdot,\cdot)$ of the primal stabilization remain fully active in the preconditioner, it remains possible to choose the polynomial degrees of primal and dual variable independent of each other.

* In the following we will refer to the strategy using the full dual order, i.e. $k_{\ast}=k$ and $ q_{\ast} = q$, as `M-f'.

* The choice using minimal dual order $k_{\ast}=1, q_{\ast}=0$ will be denoted by `M-l'.

§.§ Decoupled forward-backward solve

The strategy proposed in the previous subsection treats primal and dual variable in a monolithic fashion. Alternatively, we can (approximately) solve for these variables separately by exploiting the fact that $\underline{\mathbf{Z}}_h$ converges to zero, see <ref>. Hence, we may set $\underline{\mathbf{Z}}_h = 0$ in <ref>
and solve for $\underline{\mathbf{U}}_h$ and subsequently determine a corrected  $\underline{\mathbf{Z}}_h$ based on $\underline{\mathbf{U}}_h$. This requires, however, that $A[ \cdot, \cdot]$ on its own provides a consistent and stable discretization of the wave equation which is currently not the case. Consequently, we have to replace $A[ \cdot, \cdot]$ by an enriched bilinear form $\tilde{A}[ \cdot, \cdot]$ defined as follows:

\begin{align}
\tilde{A}[ \underline{\mathbf{U}}_h, \underline{\mathbf{Y}}_h] &=  A[ \underline{\mathbf{U}}_h, \underline{\mathbf{Y}}_h]  \nonumber \\
& + ( \underline{u}_1 , \underline{y}_1 )_{\STdata } + \frac{\lambda}{h} \sum\limits_{n=0}^{N-1} ( \underline{u}_1 , \underline{y}_1  )_{ \Sigma^n } + \sum\limits_{n=1}^{N-1}  \left\{  ( \jump{ \underline{u}_1^n } ,  \underline{y}_{2,+}^n  )_{\Omega} +  ( \jump{ \underline{u}_2^n } ,  \underline{y}_{1,+}^n  )_{\Omega}  \right\} \label{eq:tildeA-def}
% \end{aligned}
\end{align}

with the Nitsche parameter $\lambda > 0$ sufficiently large as usual. We remark that this enriched bilinear form is a minor modification of a method originally proposed and analyzed in [17]. In order to account for the presence of the jumps terms in  $\tilde{A}[ \cdot, \cdot]$, the dual stabilization has to be adapted as follows:

\[
\tilde{S}^{\ast}_h( \underline{\mathbf{Y}}_h, \underline{\mathbf{Z}}_h) = S^{\ast}_h( \underline{\mathbf{Y}}_h, \underline{\mathbf{Z}}_h) + \Delta t \sum\limits_{n=1}^{N-1}  \left\{  ( \underline{y}_{1,+}^n, \underline{z}_{1,+}^n  )_{\Omega} +  ( \underline{y}_{2,+}^n, \underline{z}_{2,+}^n  )_{\Omega}  \right\}.
\]

In $\tilde{A}$ we have also added a Luenberger observer type term on $\STdata$. It is known that this can improve the stability of the corrected dynamics (see [11, 12]), which may improve the convergence of the iterations. Since the perturbation is consistent it is straightforward to include in the error analysis above.
The first-order optimality conditions then become:
Find $(\underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h) \in \ProdFullyDiscrSpace{k}{q} \times \ProdFullyDiscrSpace{ k_{\ast} }{ q_{\ast} } $ such that

\begin{align}
(\underline{u}_1,\underline{w}_1)_{\STdata} + \tilde{A}[ \underline{\mathbf{W}}_h, \underline{\mathbf{Z}}_h ] + \tilde{S}_h(  \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h )  &= (u_{\omega},\underline{w}_1)_{ \STdata}, \quad \forall \; \underline{\mathbf{W}}_h \in \ProdFullyDiscrSpace{k}{q}, \label{eq:optimality_cond_fb_dual} \\
\tilde{A}[ \underline{\mathbf{U}}_h, \underline{\mathbf{Y}}_h ] -  \tilde{S}_h^{\ast}( \underline{ \mathbf{Y}}_h, \underline{\mathbf{Z}}_h)  &=  (u_{\omega},\underline{y}_1)_{ \STdata} , \quad \forall \;  \underline{\mathbf{Y}}_h \in \ProdFullyDiscrSpace{ k_{\ast} }{ q_{\ast} },  \label{eq:optimality_cond_fb_primal}
\end{align}

for $ \tilde{S}_h(  \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h ) := S_h(  \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h ) + S^{\uparrow \downarrow}_{\Delta t}( \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h ) $.
Assuming $k=k_{\ast}$ and $q = q_{\ast}$ we propose the following preconditioning strategy, which we will refer to as `DFB'.

1. Set  $\underline{\mathbf{Z}}_h = 0$ in <ref> and solve

\begin{equation}\label{eq:tildeA-fb-Uh}
\tilde{A}[ \underline{\mathbf{U}}_h, \underline{\mathbf{Y}}_h ]  =  (u_{\omega},\underline{y}_1)_{ \STdata} , \quad \forall \;  \underline{\mathbf{Y}}_h \in \ProdFullyDiscrSpace{ k_{\ast} }{ q_{\ast} }
\end{equation}

by marching forward in time to obtain $\underline{\mathbf{U}}_h$.

2. Then insert $\underline{\mathbf{U}}_h$ in <ref> and solve

\begin{equation}\label{eq:tildeA-fb-Zh}
\tilde{A}[ \underline{\mathbf{W}}_h, \underline{\mathbf{Z}}_h ] = (u_{\omega} - \underline{u}_1,\underline{w}_1)_{ \STdata} -  \tilde{S}_h(  \underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h ), \quad \forall \; \underline{\mathbf{W}}_h \in \ProdFullyDiscrSpace{k}{q},
\end{equation}

backward in time to obtain $\underline{\mathbf{Z}}_h$.

We point out that we integrate by parts in <ref> to write

\begin{align*}
\begin{aligned}
& \tilde{A}[ \underline{\mathbf{W}}_h, \underline{\mathbf{Z}}_h] = \sum\limits_{n=0}^{N-1} \big\{ - (\underline{w}_2,\partial_t  \underline{z}_1 + \underline{z}_2 )_{Q^n} + a(\underline{w}_1,\underline{z}_1)_{Q^n}
- ( \underline{w}_1 , \partial_t \underline{z}_2)_{Q^n}  \\
& - ( \nabla \underline{w}_1 \cdot \mathbf{n}, \underline{z}_1)_{ \Sigma^n } +  \frac{\lambda}{h} ( \underline{w}_1 , \underline{z}_1  )_{ \Sigma^n }  \big\} + ( \underline{w}_1 ,  \underline{z}_1 )_{ \STdata} -  ( \underline{w}_{2,+}^0 , \underline{z}_{1,+}^0  )_{\Omega} -  ( \underline{w}_{1,+}^0 , \underline{z}_{2,+}^0  )_{\Omega} \\
&  +  ( \underline{w}_{2,-}^N , \underline{z}_{1,-}^N  )_{\Omega} +  ( \underline{w}_{1,-}^N , \underline{z}_{2,-}^N  )_{\Omega} - \sum\limits_{n=1}^{N-2} \left\{ (  \underline{w}_{2,-}^{n+1} , \jump{ \underline{z}_1^{n+1} }  )_{\Omega} + (  \underline{w}_{1,-}^{n+1} , \jump{ \underline{z}_2^{n+1} }  )_{\Omega}   \right\},
\end{aligned}
\end{align*}

which then indeed allows to obtain $\underline{\mathbf{Z}}_h$ time-slab wise by starting at $n=N-1$ and progressing backwards in time.
We also mention that this strategy is limited to the setting in which primal and dual polynomial degrees coincide as we need to make sure that the linear systems in <ref> and <ref> are both well-posed.

§.§ Comparison

Let us briefly compare the proposed strategies. Both approaches only require ingredients available in standard dG time-stepping methods. In particular, only linear systems on the time slabs $Q^n$ need to be solved. When the full polynomial degree for the dual variables is used, i.e.  $k_{\ast}=k$ and $ q_{\ast} = q$, then the linear systems on the slabs to be inverted in the decoupled method (`DFB') are half the size of the ones in the monolithic method (`M-f'). This is because of the monolithic treatment of primal and dual variable in the latter. However, if the variant `M-l' of the monolithic method with minimal dual polynomial order is used, then the difference in size between the linear systems quickly diminishes as the primal polynomial order increases. Therefore, it seems fair to focus in the numerical experiments presented in <ref> on a comparison between the `DFB' and `M-l' methods.

§ NUMERICAL EXPERIMENTS

In this section we present some numerical experiments for the fully-discrete method defined in <ref> and a comparison of the preconditioning strategies proposed in <ref>. Our implementation uses the  library [2, 23, 22, 1]. Reproduction material for the numerical experiments is available at Zenodo [10].

§.§ Preconditioning experiments for $d=1$

We start by considering a simple experiment on the unit interval $\Omega = [0,1]$. The data domain is given by $\omega = [0,1/4] \cup [3/4,1]$ and the exact solution is $u = \cos(\pi t) \sin(\pi x)$. We measure up to $T=1/2$ so that the GCC is fulfilled.

columns/N/.style=  column name=$N$, column type/.add=> [gray].9   ,
columns/ndof/.style=  column name=ndof, column type/.add=> [gray].9     ,
columns/vanilla/.style=  column name=NoPre, string replace=3000,empty cells with=-,  ,
columns/block/.style=  column name=Block, string replace=3000,empty cells with=-,   ,
columns/MTM-lo/.style=  column name=M-l , string replace=3000,empty cells with=-  ,
columns/DFB/.style=  column name=DFB, string replace=3000,empty cells with=-  ,
columns/MTM-full/.style=  column name=M-f , string replace=3000,empty cells with=-  ,
] precond-1d-iters-order1.dat

[scale = 1.0]
height = 4.25cm,
width = 5.5cm,
every axis plot/.append style=thick,
axis y line*=left,

label style=font=,
tick label style=font=,

xlabel=   $\Delta t$  ,
y tick label style= xshift=2.0em,yshift=-0.4em ,
x tick label style= xshift=1.2em,yshift=1.1em ,
x label style=at=(axis description cs:0.35,+0.25),anchor=east,
title =    $\norm{\partial_t \tilde{e}_h }_{L^{2}(0,T;L^2(\Omega))}$    ,
legend style=at=(0.75, 0.4),anchor=north,
title style=at=(0.55,0.935),anchor=north,
[red,very thick,mark=o,forget plot]
table[x=deltat,y=DFB] precond-1d-L2L2ut-order1.dat;
[blue,very thick,mark=square,forget plot]
table[x=deltat,y=MTM-lo] precond-1d-L2L2ut-order1.dat;
[green!70!black,very thick,mark=triangle,forget plot]
table[x=deltat,y=MTM-full] precond-1d-L2L2ut-order1.dat;
[lightgray,dashed,ultra thick]
table[mark=none,x=deltat,y expr =.75*0] precond-1d-L2L2ut-order1.dat;
$ \mathcal{O}(\Delta t) $

columns/N/.style=  column name=$N$, column type/.add=> [gray].9   ,
columns/ndof/.style=  column name=ndof, column type/.add=> [gray].9     ,
columns/MTM-lo/.style=  column name=M-l, string replace=3000,empty cells with=-  ,
columns/DFB/.style=  column name=DFB, string replace=3000,empty cells with=-  ,
columns/MTM-full/.style=  column name=M-f, string replace=3000,empty cells with=-  ,
] precond-1d-iters-order2.dat

[scale = 1.0]
group style=
group size=2 by 1,
horizontal sep=4.0pt,
vertical sep=40pt,
label style=font=,
tick label style=font=,
height = 4.25cm,
width = 5.0cm,
every axis plot/.append style=thick,
axis y line*=left,
ymode = log,
xlabel=  iter  ,
y tick label style= xshift=0.25em,yshift=-0.4em ,
x tick label style= xshift=-0.1em,yshift=-0.2em ,
ytick = 1e7,1e3,1e-1,1e-5 ,
x label style=at=(axis description cs:0.525,+0.25),anchor=east,
title =  Residual ,
legend style=at=(0.6, 0.8),anchor=north,
title style=at=(0.55,0.935),anchor=north,
xmax = 177,

[red,only marks, mark =*,mark options=scale=0.5,each nth point=4 ]
table[x=iter,y=res] precond-1d-DFB-order2-residuals-ref-lvl3.dat;
[ red,only marks, mark =o,mark options=scale=0.7,each nth point=4 ]
table[x=iter,y=res] precond-1d-DFB-order2-residuals-ref-lvl4.dat;

[ blue,only marks, mark = square*,mark options=scale=0.7,forget plot,each nth point=4  ]
table[x=iter,y=res] precond-1d-MTM-lo-order2-residuals-ref-lvl3.dat;
[ blue,only marks, mark =square,mark options=scale=0.7,forget plot,each nth point=4  ]
table[x=iter,y=res] precond-1d-MTM-lo-order2-residuals-ref-lvl4.dat;
$N=8$, $N=16$

xtick = 1e-1,
axis y line*=left,

xlabel=   $\Delta t$  ,
y tick label style= xshift=-0.1em,yshift= 0.0em ,
x tick label style= xshift=-0.7em,yshift=1.4em ,
x label style=at=(axis description cs:0.75,0.1),anchor=east,
title =    $\norm{\partial_t \tilde{e}_h }_{L^{2}(0,T;L^2(\Omega))}$    ,
legend style=at=(0.35, 0.75),anchor=north,
title style=at=(0.45,0.935),anchor=north,
[blue,very thick,mark=square,forget plot]
table[x=deltat,y=MTM-lo] precond-1d-L2L2ut-order2.dat;
[green!70!black,very thick,mark=triangle,forget plot]
table[x=deltat,y=MTM-full] precond-1d-L2L2ut-order2.dat;
[red,very thick,mark=o,forget plot]
table[x=deltat,y=DFB] precond-1d-L2L2ut-order2.dat;
[lightgray,dashed,ultra thick]
table[mark=none,x=deltat,y expr =.75*0*0] precond-1d-L2L2ut-order2.dat;
$ \mathcal{O}(\Delta t^2) $

Results for one-dimensional experiment of <ref> using $h = \Delta t$.
The upper panel refers to $q=k=1$ while the lower panel uses $q=k=2$. The tables contain iteration numbers for GMRes preconditioned by the methods proposed in <ref> using tolerance $10^{-7}$. The column `ndof' gives the total number of degrees of freedom for the `M-f'-method.
The two plots on the right show convergence of $\norm{\partial_t \tilde{e}_h }_{L^{2}(0,T;L^2(\Omega))}$ with $\tilde{e}_h := u - L_{\Delta t} \underline{u}_1$. In the history of the residual only every fourth iteration is shown as a data point.

The results for piecewise affine and quadratic approximation are displayed in <ref>.
In addition to the preconditioning approaches introduced in <ref> we also show the iteration numbers for a `Block'-Jacobi preconditioner in which the coupling terms $S^{\uparrow \downarrow}_{\Delta t}( \cdot, \cdot )$ have been completely removed. The poor performance of this approach shows that including some form of communication between time-slabs is crucial. Concerning the approaches from  <ref> we make the follwing observations:

* The iteration numbers for all the methods increase at least linearly with $N$. For the monolithic method this is evident as the term $S^{\uparrow \downarrow}_{\Delta t}(\cdot,\cdot)$ which is only approximately represented by the preconditioner scales with $1/\Delta t \sim N$.

* The method `M-f' has the lowest iteration numbers. However, it is also the most expensive.
Fortunately, the iteration numbers for the `M-l' method are nearly the same as for the  `M-f' method when $q=k >1$ and its cost is similar to that of the `DFB' method. Hence, the `M-l' method is a promising candidate to be used for higher-order discretizations.

* For $q=k = 1$ the  `M-l' method cannot compete with the  `DFB' method due to the poorly resolved dual variable. The plot of the error in the upper panel of <ref> also shows that the latter enjoys some super-convergence properties for piecewise affine discretizations. This is courtesy of the additional jumps terms which have been integrated into the modified bilinear form $\tilde{A}[ \cdot, \cdot]$ in <ref>. These terms cannot be accommodated in the method `M-l' as they are incompatable with the monolithic forward sweep of <ref>.
A minor drawback of the `DFB' method is that the corresponding residuals are very high in the first GMRes iterations and they also increase with $N$ as <ref> shows. Here, the  `M-l' method behaves more favourably. Nevertheless, the `DFB' method is attractive for  $q=k = 1$.

§.§ Unit cube geometry with data given all around

Let us now consider a computationally more demanding experiment on the unit cube $\Omega = [0,1]^3$ with data domain $\omega = \Omega \setminus [1/4,3/4]^3$ and final time $T= 1/2$.  The data domain is sketched in the center of <ref>. The reference solution is given by

\[
u(t,x,y,z) = \cos(\sqrt{3} \pi t) \sin(\pi x) \sin(\pi y) \sin(\pi z).
\]

We set $\Delta t = h$ for this experiment and report the errors in the $\norm{ \cdot }_{L^{\infty}(0,T;L^2(\Omega))}$ and $\norm{ \partial_t( \cdot) }_{L^{2}(0,T;L^2(\Omega))}$-norms[For simplicity we measured the $L^{2}(0,T;L^2(\Omega))$ instead of the $L^{2}(0,T;H^{-1}(\Omega))$ norm.] in terms of the polynomial order $q=k$ in <ref>.  For the wave displacement some superconvergence can be observed, as has also been reported in previous publications [9, 7] for the spatially one-dimensional case. The velocity marginally superconverges for the
`DFB' method when piecewise affine approximation is used as already observed in the previous experiment.
Overall, the numerical convergence rates are in line with the theoretical result from <ref>.

[scale = 1.0]
[ ]

height = 4.75cm,
width = 5.7cm,
ymax = 2e-1,
axis y line*=left,
xlabel=  $\Delta t = h$,
x label style=at=(axis description cs:0.35,+0.225),anchor=east,
legend style =  column sep = 10pt, legend columns = 1, legend to name = grouplegend,,
title style=at=(0.5,1.0735),anchor=north,
title =   $\norm{ u - \mathcal{L}_{\Delta t} \underline{u}_1 }_{L^{\infty}(0,T;L^2(\Omega))}$ ,
legend style=at=(0.5,-0.1),anchor=north,

[blue,only marks,mark=square,mark options=scale=1.25,forget plot]
table[x=deltat,y=MTM-lo] precond-3d-GCC-LinftyL2u-order1.dat;
[blue,very thick]
table[x=deltat,y=MTM-lo] precond-3d-GCC-LinftyL2u-order1.dat; $q=1$
[red,very thick,only marks, mark=o,mark options=scale=1.25,forget plot]
table[x=deltat,y=DFB] precond-3d-GCC-LinftyL2u-order1.dat;

[blue,only marks,mark=square,mark options=scale=1.25,forget plot]
table[x=deltat,y=MTM-lo] precond-3d-GCC-LinftyL2u-order2.dat;
[blue,very thick,dashed]
table[x=deltat,y=MTM-lo] precond-3d-GCC-LinftyL2u-order2.dat; $q=2$
[red,very thick,only marks, mark=o,mark options=scale=1.25,forget plot]
table[x=deltat,y=DFB] precond-3d-GCC-LinftyL2u-order2.dat;

[blue,only marks,mark=square,mark options=scale=1.25,forget plot]
table[x=deltat,y=MTM-lo] precond-3d-GCC-LinftyL2u-order3.dat;
[blue,very thick,densely dotted]
table[x=deltat,y=MTM-lo] precond-3d-GCC-LinftyL2u-order3.dat; $q=3$
[red,very thick,only marks, mark=o,mark options=scale=1.25,forget plot]
table[x=deltat,y=DFB] precond-3d-GCC-LinftyL2u-order3.dat;

[lightgray,dashed,ultra thick]
table[mark=none,x=deltat,y expr =.15*0] precond-3d-GCC-LinftyL2u-order1.dat;
[lightgray,dotted,ultra thick]
table[mark=none,x=deltat,y expr =.15*0*0 ] precond-3d-GCC-LinftyL2u-order2.dat;
[lightgray,dashdotted,ultra thick]
table[mark=none,x=deltat,y expr =.05*0*0*0] precond-3d-GCC-LinftyL2u-order3.dat;

at ($(displ) + (-0.0cm,-2.55cm)$) <ref>;

(geom) at (1.0,2.25) ;
[draw,fill=white] (la) at (1.15,3.25)  $\omega $ ;
[] (Z2) at (-2.25,0.55) ;
(err) at (1.25,-0.4) ;
[draw,fill=white] (la) at (1.05,0.875) $\vert u(0,\cdot) - \mathcal{L}_{\Delta t} \underline{u}_1(0,\cdot) \vert$;
[lightgray,ultra thick, ->] (-0.3,-0.55 ) – (-3.145, 1.695);

height = 4.75cm,
width = 5.7cm,
ymax = 1e-0,
axis y line*=left,
xlabel=  $\Delta t = h$,
legend style =  column sep = 10pt, legend columns = 1, legend to name = grouplegendv,,
x label style=at=(axis description cs:0.35,+0.225),anchor=east,
title style=at=(0.5,1.0735),anchor=north,
title =   $\norm{ \partial_t( u - \mathcal{L}_{\Delta t} \underline{u}_1 ) }_{L^{2}(0,T;L^2(\Omega))}$ ,
legend style=at=(0.5,-0.1),anchor=north,

[blue,only marks,mark=square,mark options=scale=1.25,forget plot]
table[x=deltat,y=MTM-lo] precond-3d-GCC-L2L2ut-order1.dat;
[blue,very thick,forget plot]
table[x=deltat,y=MTM-lo] precond-3d-GCC-L2L2ut-order1.dat;
[red,very thick,only marks, mark=o,mark options=scale=1.25,forget plot]
table[x=deltat,y=DFB] precond-3d-GCC-L2L2ut-order1.dat;

[blue,only marks,mark=square,mark options=scale=1.25,forget plot]
table[x=deltat,y=MTM-lo] precond-3d-GCC-L2L2ut-order2.dat;
[blue,very thick,dashed,forget plot]
table[x=deltat,y=MTM-lo] precond-3d-GCC-L2L2ut-order2.dat;
[red,very thick,only marks, mark=o,mark options=scale=1.25,forget plot]
table[x=deltat,y=DFB] precond-3d-GCC-L2L2ut-order2.dat;

[blue,only marks,mark=square,mark options=scale=1.25,forget plot]
table[x=deltat,y=MTM-lo] precond-3d-GCC-L2L2ut-order3.dat;
[blue,very thick,densely dotted,forget plot]
table[x=deltat,y=MTM-lo] precond-3d-GCC-L2L2ut-order3.dat;
[red,very thick,only marks, mark=o,mark options=scale=1.25,forget plot]
table[x=deltat,y=DFB] precond-3d-GCC-L2L2ut-order3.dat;

[lightgray,dashed,ultra thick]
table[mark=none,x=deltat,y expr =1.3*0] precond-3d-GCC-L2L2ut-order1.dat;  $ \mathcal{O}(\Delta t) $
[lightgray,dotted,ultra thick]
table[mark=none,x=deltat,y expr =1.4*0*0 ] precond-3d-GCC-L2L2ut-order2.dat;  $ \mathcal{O}((\Delta t)^2) $
[lightgray,dashdotted,ultra thick]
table[mark=none,x=deltat,y expr =.5*0*0*0] precond-3d-GCC-L2L2ut-order3.dat;  $ \mathcal{O}((\Delta t)^3) $

at ($(vel) + (-0.0cm,-2.55cm)$) <ref>;

Convergence plots for the unit cube with data given in  $\omega = \Omega \setminus [1/4,3/4]^3$. Results are given for the `M-l'  and `DFB' methods. The central figure displays the absolute error on the indicated refinement level.

every head row/.style=
before row=
3cDFB     3cM-l

after row=,
columns/N/.style=  column name=$N$, column type/.add=> [gray].9   ,
columns/ndof-tot-DFB/.style=  column name=ndof-tot    ,
columns/ndof-inv-DFB/.style=  column name=ndof-inv    ,
columns/DFB-iter/.style=  column name=iter ,
columns/ndof-tot-MTM-lo/.style=  column name=ndof-tot    ,
columns/ndof-inv-MTM-lo/.style=  column name=ndof-inv    ,
columns/MTM-lo-iter/.style=  column name=iter ,
] precond-3d-GCC-iters-order1.dat

every head row/.style=
before row=
3cDFB     3cM-l

after row=,
columns/N/.style=  column name=$N$, column type/.add=> [gray].9   ,
columns/ndof-tot-DFB/.style=  column name=ndof-tot    ,
columns/ndof-inv-DFB/.style=  column name=ndof-inv    ,
columns/DFB-iter/.style=  column name=iter ,
columns/ndof-tot-MTM-lo/.style=  column name=ndof-tot    ,
columns/ndof-inv-MTM-lo/.style=  column name=ndof-inv    ,
columns/MTM-lo-iter/.style=  column name=iter ,
] precond-3d-GCC-iters-order2.dat

every head row/.style=
before row=
3cDFB     3cM-l

after row=,
columns/N/.style=  column name=$N$, column type/.add=> [gray].9   ,
columns/ndof-tot-DFB/.style=  column name=ndof-tot    ,
columns/ndof-inv-DFB/.style=  column name=ndof-inv    ,
columns/DFB-iter/.style=  column name=iter ,
columns/ndof-tot-MTM-lo/.style=  column name=ndof-tot    ,
columns/ndof-inv-MTM-lo/.style=  column name=ndof-inv    ,
columns/MTM-lo-iter/.style=  column name=iter ,
] precond-3d-GCC-iters-order3.dat

GMRes iteration numbers for the unit cube with data given in  $\omega = \Omega \setminus [1/4,3/4]^3$ using tolerance $10^{-5}$. The left half of the tables displays the results for the `DFB' method and the right half for the `M-l' method. The column `ndof-tot' gives the total number of degrees of freedom of the entire space-time system for the respective method, whereas the column `ndof-inv' indicates the size of the linear system which needs to be inverted to apply the preconditioner.

The GMRes iteration numbers displayed in <ref> confirm our observations from the one-dimensional example: The `DFB' method is preferred for  $q=k = 1$, whereas the iteration numbers for the `M-l' method are far lower for $q=k > 1$. Since the degrees of freedom in the linear system which needs to be inverted to apply the preconditioner (column `ndof-inv' in the table) are only marginally larger for the `M-l' method, this is clearly the method of choice for higher order discretizations.
Even though this preconditioning strategy is not robust with respect to refinement of the discretization, it still allows us to solve a globally-coupled unique continuation problem posed on a four dimensional space-time geometry up to a reasonable level of accuracy. Solving the $d+1$ dimensional space-time systems directly would require to invert a linear system which involves a factor[The ratio of the columns `ndof-tot' and  `ndof-inv' in <ref>.] $N$ (respectively $2N$ for the `DFB' method) more degrees of freedom than the linear system to be factorized to apply the preconditioner.

§.§ Unit cube geometry without GCC

[scale = 1.0]
[ ]

height = 4.75cm,
width = 5.7cm,
ymax = 7e-1,
axis y line*=left,
xlabel=  $\Delta t = h$,
x label style=at=(axis description cs:0.55,+0.21),anchor=east,
legend style =  column sep = 10pt, legend columns = 1, legend to name = grouplegendn,,
title =   $\norm{ u - \mathcal{L}_{\Delta t} \underline{u}_1 }_{L^{\infty}(0,T;L^2(\Omega))}$ ,
title style=at=(0.5,1.0735),anchor=north,
legend style=at=(0.5,-0.1),anchor=north,

[blue,only marks,mark=square,mark options=scale=1.25,forget plot]
table[x=deltat,y=MTM-lo] precond-3d-noGCC-LinftyL2u-order1.dat;
[blue,very thick]
table[x=deltat,y=MTM-lo] precond-3d-noGCC-LinftyL2u-order1.dat; $q=1$
[red,very thick,forget plot]
table[x=deltat,y=DFB] precond-3d-noGCC-LinftyL2u-order1.dat;
[red,very thick,only marks, mark=o,mark options=scale=1.25,forget plot]
table[x=deltat,y=DFB] precond-3d-noGCC-LinftyL2u-order1.dat;

[blue,only marks,mark=square,mark options=scale=1.25,forget plot]
table[x=deltat,y=MTM-lo] precond-3d-noGCC-LinftyL2u-order2.dat;
[blue,very thick,dashed]
table[x=deltat,y=MTM-lo] precond-3d-noGCC-LinftyL2u-order2.dat; $q=2$
[red,very thick,only marks, mark=o,mark options=scale=1.25,forget plot]
table[x=deltat,y=DFB] precond-3d-noGCC-LinftyL2u-order2.dat;
[red,very thick,dashed,forget plot]
table[x=deltat,y=DFB] precond-3d-noGCC-LinftyL2u-order2.dat;

[blue,only marks,mark=square,mark options=scale=1.25,forget plot]
table[x=deltat,y=MTM-lo] precond-3d-noGCC-LinftyL2u-order3.dat;
[blue,very thick,densely dotted]
table[x=deltat,y=MTM-lo] precond-3d-noGCC-LinftyL2u-order3.dat; $q=3$
[red,very thick,only marks, mark=o,mark options=scale=1.25,forget plot]
table[x=deltat,y=DFB] precond-3d-noGCC-LinftyL2u-order3.dat;
[red,very thick,densely dotted]
table[x=deltat,y=DFB] precond-3d-noGCC-LinftyL2u-order3.dat;

[lightgray,dashed,ultra thick]
table[mark=none,x=deltat,y expr =.9*0] precond-3d-noGCC-LinftyL2u-order1.dat;
[lightgray,dotted,ultra thick]
table[mark=none,x=deltat,y expr =.35*0*0 ] precond-3d-noGCC-LinftyL2u-order2.dat;
[lightgray,dashdotted,ultra thick]
table[mark=none,x=deltat,y expr =.25*0*0*0] precond-3d-noGCC-LinftyL2u-order3.dat;

at ($(displ) + (-0.0cm,-2.55cm)$) <ref>;

(geom) at (1.15,2.25) ;
[draw,fill=white] (la) at (0.46,3.15)  $\omega $ ;
[] (Z2) at (-2.25,0.55) ;
(err) at (1.0,-0.4) ;
[draw,fill=white] (la) at (1.05,0.78) $\vert u(0,\cdot) - \mathcal{L}_{\Delta t} \underline{u}_1(0,\cdot) \vert$;
[lightgray,ultra thick, ->] (-0.3,-0.55 ) – (-3.145, 2.195);

height = 4.75cm,
width = 5.7cm,
ymax = 1e-0,
xlabel=  $\Delta t = h$,
axis y line*=left,
legend style =  column sep = 10pt, legend columns = 1, legend to name = grouplegendnv,,
x label style=at=(axis description cs:0.55,+0.21),anchor=east,
title =   $\norm{ \partial_t( u - \mathcal{L}_{\Delta t} \underline{u}_1 ) }_{L^{2}(0,T;L^2(\Omega))}$ ,
legend style=at=(0.5,-0.1),anchor=north,
title style=at=(0.5,1.0735),anchor=north,

[lightgray,dashed,ultra thick]
table[mark=none,x=deltat,y expr =1.3*0] precond-3d-noGCC-L2L2ut-order1.dat;  $ \mathcal{O}(\Delta t) $
[lightgray,dotted,ultra thick]
table[mark=none,x=deltat,y expr =1.4*0*0 ] precond-3d-noGCC-L2L2ut-order2.dat;  $ \mathcal{O}((\Delta t)^2) $
[lightgray,dashdotted,ultra thick]
table[mark=none,x=deltat,y expr =.5*0*0*0] precond-3d-noGCC-L2L2ut-order3.dat;  $ \mathcal{O}((\Delta t)^3) $

[blue,only marks,mark=square,mark options=scale=1.25,forget plot]
table[x=deltat,y=MTM-lo] precond-3d-noGCC-L2L2ut-order1.dat;
[blue,very thick,forget plot]
table[x=deltat,y=MTM-lo] precond-3d-noGCC-L2L2ut-order1.dat; $q=1$
[red,very thick,forget plot]
table[x=deltat,y=DFB,forget plot] precond-3d-noGCC-L2L2ut-order1.dat;
[red,very thick,only marks, mark=o,mark options=scale=1.25,forget plot]
table[x=deltat,y=DFB] precond-3d-noGCC-L2L2ut-order1.dat;

[blue,only marks,mark=square,mark options=scale=1.25,forget plot]
table[x=deltat,y=MTM-lo] precond-3d-noGCC-L2L2ut-order2.dat;
[blue,very thick,dashed,forget plot]
table[x=deltat,y=MTM-lo] precond-3d-noGCC-L2L2ut-order2.dat; $q=2$
[red,very thick,only marks, mark=o,mark options=scale=1.25,forget plot]
table[x=deltat,y=DFB] precond-3d-noGCC-L2L2ut-order2.dat;
[red,very thick,dashed,forget plot]
table[x=deltat,y=DFB] precond-3d-noGCC-L2L2ut-order2.dat;

[blue,only marks,mark=square,mark options=scale=1.25,forget plot]
table[x=deltat,y=MTM-lo] precond-3d-noGCC-L2L2ut-order3.dat;
[blue,very thick,densely dotted,forget plot]
table[x=deltat,y=MTM-lo] precond-3d-noGCC-L2L2ut-order3.dat; $q=3$
[red,very thick,only marks, mark=o,mark options=scale=1.25,forget plot]
table[x=deltat,y=DFB] precond-3d-noGCC-L2L2ut-order3.dat;
[red,very thick,densely dotted,forget plot]
table[x=deltat,y=DFB] precond-3d-noGCC-L2L2ut-order3.dat;

at ($(vel) + (-0.0cm,-2.55cm)$) <ref>;

every head row/.style=
before row=
3cDFB     3cM-l

after row=,
columns/N/.style=  column name=$N$, column type/.add=> [gray].9   ,
columns/ndof-tot-DFB/.style=  column name=ndof-tot    ,
columns/ndof-inv-DFB/.style=  column name=ndof-inv    ,
columns/DFB-iter/.style=  column name=iter ,
columns/ndof-tot-MTM-lo/.style=  column name=ndof-tot    ,
columns/ndof-inv-MTM-lo/.style=  column name=ndof-inv    ,
columns/MTM-lo-iter/.style=  column name=iter ,
] precond-3d-noGCC-iters-order2.dat

On top: Convergence plots for the unit cube with data given in merely in $\omega = [0,1/4] \times [0,1]^2$. Results are given for the `M-l'  and `DFB' methods. Bottom: GMRes iteration numbers for $q=k=2$.

Let us change the setup from <ref> by removing part of the data. If we set
$\omega = [0,1/4] \times [0,1]^2$ keeping the final time $T=1/2$, then the GCC, which is crucial for our analysis, will be violated.
The numerical results shown in <ref> indicate that the error in this case can only be reduced to a certain threshold (depending on the type of discretization) at which stagnation occurs. Interestingly, the accuracies achieved with the `DFB' and `M-l' methods differ significantly. For $q=k \in \{2,3\}$ the latter is about one order of magnitude more accurate than the former. On the other hand, the GMRes iteration numbers for the `M-l' method are more adversely affected by the lack of the GCC than for the `DFM' method.
In the end, both methods are hampered by the fact that the error for the reconstruction of the wave displacement outside the data domain is fairly large and will not converge to zero (at a reasonable speed) when the discretization is refined as the central plot in <ref> shows.

[scale = 1.0]
[ ]

height = 4.15cm,
width = 5.45cm,
ymax = 7e-1,
axis y line*=left,
xlabel=  $\Delta t = h$,
x label style=at=(axis description cs:0.75,+0.31),anchor=east,
y tick label style= xshift=2.85em,yshift=0.6em ,
ytick = 1e-3,1e-6,
legend style =  column sep = 10pt, legend columns = 1, legend to name = grouplegendnoGCC,,

title =  $ u$ in $ \Omega$ (solid) vs $ B_t $ (dashed) ,
title style=at=(0.5,1.0735),anchor=north,
legend style=at=(0.5,-0.0),anchor=north,

[red,very thick,mark=*]
table[x=deltat,y=LinftyL2u-all] noGCC-restricted-1d-order1.dat; $q=k=1$
[blue,very thick,mark=triangle]
table[x=deltat,y=LinftyL2u-all] noGCC-restricted-1d-order2.dat;  $q=k=2$
[green!70!black,very thick,mark=x]
table[x=deltat,y=LinftyL2u-all] noGCC-restricted-1d-order3.dat;  $q=k=3$

[red,very thick,dashed,forget plot]
table[x=deltat,y=LinftyL2u-restrict] noGCC-restricted-1d-order1.dat;
[blue,very thick,dashed,forget plot]
table[x=deltat,y=LinftyL2u-restrict] noGCC-restricted-1d-order2.dat;
[green!70!black,very thick,dashed,forget plot]
table[x=deltat,y=LinftyL2u-restrict] noGCC-restricted-1d-order3.dat;

[red,only marks,mark=*,forget plot]
table[x=deltat,y=LinftyL2u-restrict] noGCC-restricted-1d-order1.dat;
[blue,only marks,mark=triangle,forget plot]
table[x=deltat,y=LinftyL2u-restrict] noGCC-restricted-1d-order2.dat;
[green!70!black,only marks,mark=x,forget plot]
table[x=deltat,y=LinftyL2u-restrict] noGCC-restricted-1d-order3.dat;

[lightgray,dashed,ultra thick,forget plot]
table[mark=none,x=deltat,y expr =0.07*0] noGCC-restricted-1d-order1.dat;  $ \mathcal{O}(\Delta t) $
[lightgray,dotted,ultra thick,forget plot]
table[mark=none,x=deltat,y expr =.04*0*0 ] noGCC-restricted-1d-order2.dat;  $ \mathcal{O}((\Delta t)^2) $
[lightgray,dashdotted,ultra thick,forget plot]
table[mark=none,x=deltat,y expr =.02*0*0*0] noGCC-restricted-1d-order3.dat;  $ \mathcal{O}((\Delta t)^3) $

[draw,circle,ultra thick,lightgray] (Z1) at (axis cs:0.03125,4.5e-6) ;

at ($(displ) + (-0.0cm,-2.25cm)$) <ref>;

(abserr) at (0,1) ;
[->,ultra thick] (-2.35,-0.95) – node [right,at end] $x$  (1.0,-0.95);
[->,ultra thick] (-2.35,-0.95) – node [above,at end] $t$  (-2.35,2.5);
[white,ultra thick,dashed] (-1.05,-0.5 ) – (-0.3, 1.0);
[white,ultra thick,dashed] (-0.3,1.0 ) – (-1.05, 2.4);
[] (B) at (-1.3,0.95)   $\textcolor{white}{B}$ ;
[lightgray,ultra thick, ->] (-2.0,1.15 ) – (-4.8, 0.95);

height = 4.15cm,
width = 5.45cm,
ymax = 1e-0,
axis y line*=left,
xlabel=  $\Delta t = h$,
ytick = 1e-3,1e-5,
legend style =  column sep = 10pt, legend columns = 1, legend to name = grouplegendnoGCCv,,
x label style=at=(axis description cs:0.75,+0.31),anchor=east,
y tick label style= xshift=2.85em,yshift=-0.4em ,

title =  $\partial_t u$ in $ \Omega$ (solid) vs $ B_t $ (dashed) ,
title style=at=(0.5,1.0735),anchor=north,
legend style=at=(0.5,-0.1),anchor=north,
[red,very thick,mark=*,forget plot]
table[x=deltat,y=L2L2ut-all] noGCC-restricted-1d-order1.dat; $q=k=1$
[blue,very thick,mark=triangle,forget plot]
table[x=deltat,y=L2L2ut-all] noGCC-restricted-1d-order2.dat;  $q=k=2$
[green!70!black,very thick,mark=x,forget plot]
table[x=deltat,y=L2L2ut-all] noGCC-restricted-1d-order3.dat;  $q=k=3$

[red,very thick,mark=*,dashed,forget plot]
table[x=deltat,y=L2L2ut-restrict] noGCC-restricted-1d-order1.dat;
[blue,very thick,mark=triangle,dashed,forget plot]
table[x=deltat,y=L2L2ut-restrict] noGCC-restricted-1d-order2.dat;
[green!70!black,very thick,mark=x,dashed,forget plot]
table[x=deltat,y=L2L2ut-restrict] noGCC-restricted-1d-order3.dat;

[lightgray,dashed,ultra thick]
table[mark=none,x=deltat,y expr =1.3*0] noGCC-restricted-1d-order1.dat;  $ \mathcal{O}(\Delta t) $
[lightgray,dotted,ultra thick]
table[mark=none,x=deltat,y expr =0.4*0*0 ] noGCC-restricted-1d-order2.dat;  $ \mathcal{O}((\Delta t)^2) $
[lightgray,dashdotted,ultra thick]
table[mark=none,x=deltat,y expr =.1*0*0*0] noGCC-restricted-1d-order3.dat;  $ \mathcal{O}((\Delta t)^3) $

at ($(vel) + (-0.0cm,-2.25cm)$) <ref>;

Results for $\Omega = [0,1]$ with data given in $[0,1/4]$ computed with the `M-f' method.
On the left: $\norm{ u - \mathcal{L}_{\Delta t} \underline{u}_1 }_{L^{\infty}(0,T;L^2(\Omega))}$ (solid lines) and $\norm{ u - \mathcal{L}_{\Delta t} \underline{u}_1 }_{L^{\infty}(0,T;L^2(B_t))}$ (dashed lines) for the restricted set $B_t$ defined in <ref>. Similarly on the right for the time derivative in the $L^2-L^2$-norm.

To further investigate this phenomenon, let us consider an analogous experiment in one spatial dimension. Here, we consider $\Omega = [0,1]$ with data domain $\omega = [0,1/4]$, final time $T=1/2$ and exact solution $u(t,x) = \cos(\pi t) \sin(\pi x)$. The results shown in <ref> are similar to the spatially three-dimensional case. If we measure in the full spatial domain $\Omega$ (solid lines in the plot), then the error appears to stagnate. The method is unable to reconstruct the solution far away from the data domain as the space-time plot of the absolute error in the center of <ref> demonstrates. However, if we settle for measuring the spatial error in the sets $B_t$ defined as

\begin{equation}\label{def:Bt}
B = \cup_{t \in [0,T]} B_t, \quad B_t := \{ x \leq 1/4 + t \text{ for }  t \leq 1/4 \mid x \leq 3/4 - t \text{ for }  t \geq 1/4\},
\end{equation}

then optimal convergence can be observed (see the dashed lines in <ref>). It would hence be interesting to investigate whether the Lipschitz stability result from <ref> could be adapted to cover the considered setting.

§ INTERPOLATION IN SPACE-TIME

In reference [21] an interpolation theory for unfitted space-time finite element spaces has been developed.
Here we recall the results for the fitted case and give some extensions which are
required for the applications in this paper. The final aim of this appendix is to prove <ref> and <ref>.

We first define interpolation on the extended time slices
$\hat{Q}^n = I_n \times \hat{\Omega}$ and $\hat{Q} = I \times \hat{\Omega}$
and then obtain the result on the orignal domains by means of a Sobolev extension argument.
To this end, let

\[
\hat{V}_h^k := \{ u \in H^1(\hat{\Omega}) \mid u|_K \in \mathcal{P}_k(K) \; \forall K \in \hat{\mathcal{T}}_h \}, \qquad \FullyDiscrSpaceHat{k}{q} :=  \otimes_{n=0}^{N-1} \mathcal{P}^q(I_n) \otimes \hat{V}_h^k.
\]

We start in <ref> by defining a purely temporal interpolation operator and introduce its purely spatial counterpart in <ref>.
In <ref> we then concatenate these operators to achieve space-time interpolation. To seprate spatial and temporal continuity, we define $t$-anisotropic Sobolev spaces
H^q,k(Q̂) := { u |∂_t^p D_x^αu ∈L^2(Q̂), α/k + p/q ≤1  }.
For $q=k$ the isotropic Sobolev spaces are retained: $H^{k,k}(\hat{Q}) = H^{k}(\hat{Q})$.

§.§ Interpolation in time

For $q \in \mathbb{N}_0$ define $\pi^{q}_n: H^{0,0}(\hat{Q}^n) \rightarrow \mathcal{P}^q(I_n) \otimes L^2(\hat{\Omega})  $ by

\begin{equation}\label{eq:Def temporal L2-projection}
(  \pi^{q}_n(u) - u,v)_{\hat{Q}^{n}} = 0 \quad \forall v \in  \mathcal{P}^{q}(I_n) \otimes L^2(\hat{\Omega}).
\end{equation}

On sufficiently smooth functions, $\pi^{q}_n$ coincides with a standard temporal $L^2(I_{n})$-projection $I^{q}_{n}$ defined by

\begin{equation*}
(I^{q}_{n} u,\chi)_{I_{n}} = (u,\chi)_{I_{n}} \quad \forall \; \chi \in \mathcal{P}^{q}(I_n).
\end{equation*}

For $u \in C^{0}( \hat{\Omega} ;L^2(I_{n}))$ there holds

\begin{equation*}
\pi^{q}_n u = I_n^{q}u \text{ in the } L^2(\hat{Q}^{n}) \text{ sense.}
\end{equation*}

Moreover, for $u \in H^{0,1}(\hat{Q}^{n})$ the temporal interpolation commutes with spatial derivatives

\begin{equation*}
\partial_{x_{j}} \pi^{q}_{n} u = \pi^{q}_{n} \partial_{x_{j}}u, \quad j=1,\ldots,d.
\end{equation*}

See <cit.>.

We denote the temporal interpolation error by $\mathcal{L}^{q}_n := \mathrm{id} - \pi^{q}_n $.
The stability and approximation properties of  $ \pi^{q}_{n}$ are studied below.

The following stability and approximation results hold.

* $\norm{  \pi^{q}_{n} u }_{\hat{Q}^{n}} \leq \norm{u}_{\hat{Q}^{n}} \quad \forall u \in L^2(\hat{Q}^{n})$ for $n=1,\ldots,N-1$.

* $\norm{ \partial_{t} \pi^{q}_{n}  u}_{\hat{Q}^{n}} \leq C \norm{u}_{H^{1,0}(\hat{Q}^{n})} \quad \forall u \in H^{1,0}(\hat{Q}^{n})$ for $n=1,\ldots,N-1$.

* For $n=1,\ldots,N-1 $ we have $\norm{\mathcal{L}^{q}_{n} u}_{ \hat{Q}^{n}} \leq C \Delta{t}^{l} \norm{u}_{H^{l,0}(\hat{Q}^{n})} $ $\forall u \in H^{l,0}(\hat{Q}^{n})$ and $ l \in \{0,\ldots,q+1 \}$.

* For $n=1,\ldots,N-1 $ we have $\norm{ \partial_{t} \mathcal{L}^{q}_{n} u}_{ \hat{Q}^{n}} \leq C \Delta{t}^{l-1} \norm{u}_{H^{l,0}( \hat{Q}^{n})}$ for all $ u \in H^{l,0}(\hat{Q}^{n})$ and $ l \in \{1,\ldots,q+1 \}$.

The results for $ \Pi^{q}_{n}$ have been shown in <cit.>.
Note that in this reference the results (c) and (d) are only stated for $ l \in \{1,\ldots,q+1 \}$ and  $ l \in \{2,\ldots,q+1 \}$, respectively.
However, the case $l=0$ in (c) immediately follows from (a) and the case $l=1$ in (d) follows from (b).

§.§ Interpolation in space

In this section we consider a generic interval $\tilde{I} \subset (0,T)$ and set $ \tilde{Q} = \tilde{I} \times \hat{\Omega}$. Later, we will choose $\tilde{I} = (0,T)$ to define an interpolation operator into the semi-discrete space $\SemiDiscrSpace{k}$ and $\tilde{I} = I_n$ to define a time-slab-wise interpolation into the fully-discrete space $\FullyDiscrSpace{k}{q}$.

Define the projector $ \pi^{k} : L^2(\tilde{Q}) \rightarrow L^2(\tilde{I}) \otimes  \hat{V}_{h}^k $, so that for $u \in L^2(\tilde{Q})$

\begin{equation}
(\pi^{k}(u) -u,v_{h})_{\tilde{Q}} = 0, \quad \text{ for all } v_{h} \in L^2(\tilde{I}) \otimes  \hat{V}_{h}^k.
\end{equation}

On sufficiently smooth functions $\pi^{k}$ coincides with a spatial $L^2( \hat{\Omega} )$ projector $I^{k}: L^2( \hat{\Omega} ) \rightarrow \hat{V}_{h}^k$.

Let $I^{k}: L^2( \hat{\Omega})  \rightarrow \hat{V}_{h}^k$ be the spatial $L^2( \hat{\Omega} )$ projection operator. Then there holds for $u \in C^{0}( \tilde{I},L^2( \hat{\Omega} ))$

\begin{equation}
\pi^{k}u = I^{k}u \quad \text{in the } L^2( \tilde{Q} ) \text{ sense.}
\end{equation}

Further, for $u \in H^{1,0}(\tilde{Q})$ there holds:

\begin{equation}
\pi^{k} \partial_{t} u = \partial_{t} \pi^{k} u.
\end{equation}

See <cit.>.

Let $ \mathcal{L}^{k} = \mathrm{id} - \pi^{k}$ denote the spatial interpolation error. The next lemma shows that the projection into the semi-discrete space is stable and has the expected approximation properties.

For the projector $\pi^{k}: L^2(\tilde{Q}) \rightarrow L^2(\tilde{I}) \otimes  \hat{V}_{h}^k $ there hold the following stability and approximation results:

* $\norm{\pi^{k}u}_{\tilde{Q}} \leq \norm{u}_{\tilde{Q}} \quad \forall u \in L^2(\tilde{Q})$,

* $\norm{ \nabla \pi^{k} u}_{ \tilde{Q} } \leq C \norm{ u}_{H^{0,1}(\tilde{Q})} \quad \forall u \in H^{0,1}(\tilde{Q})$,

* $\norm{\mathcal{L}^{k} u}_{\tilde{Q}} \leq C h^{s} \norm{u}_{H^{0,s}(\tilde{Q})} \quad \forall u \in H^{0,s}(\tilde{Q}), \quad s \in \{1,\ldots,k+1 \}$,

* $\norm{ \nabla \mathcal{L}^{k} u}_{\tilde{Q}} \leq C h^{s-1} \norm{u}_{H^{s,0}(\tilde{Q})} \quad \forall u \in H^{0,s}(\tilde{Q}), \quad s \in \{1,\ldots,k+1 \}$.

* Also, for all $ \forall u \in H^{0,s}(\tilde{Q}) \cap C^{0}(\tilde{I},H^2(\hat{\Omega} ))$ and $ s \in \{2,\ldots,k+1 \}$ it holds that

\[ \left( \int\limits_{\tilde{I}} \sum\limits_{K \in \hat{\mathcal{T}}_h } \norm{\mathcal{L}^{k} u}_{ H^2(K) }^2 \; \dT \right)^{1/2} \leq C h^{s-2} \norm{u}_{H^{0,s}(\tilde{Q})}. \]

We only show (e) since the results (a)-(d) have already been established in <cit.>.
Note that since $u \in C^{0}(\tilde{I},L^2( \hat{\Omega} ))$ we have $\mathcal{L}^{k}u = u - I^{k} u$.
Moreover, as $u(t,\cdot) \in H^2(\hat{\Omega} )$ and $d \leq 3$ the standard Lagrange interpolation operator $I^{L}:H^2( \hat{\Omega}  ) \rightarrow \hat{V}_{h}^k $ is well-defined
and has the approximation property
I^Lu - u_H^2(K) ≤C h^s-2 u_H^s(K),   for  s ∈{ 2, …, k+1 }.
Then we have

\[ \norm{\mathcal{L}^{k} u}_{ H^2(K) } \leq \norm{ u -  I^{L}u }_{H^2(K)} + \norm{ I^{L}u - I^{k} u }_{H^2(K)}.   \]

The integrand in the second term is a discrete function to which the usual inverse inequalities can be applied:

\begin{align*}
\norm{ I^{L}u - I^{k} u }_{H^2(K)} &\leq  h^{-2} C \norm{ I^{L}u - I^{k} u }_{L^2(K)}  \\
& \leq h^{-2} C \left(  \norm{ I^{L}u - u }_{L^2(K)} + \norm{ I^{k} u - u }_{L^2(K)}   \right).
\end{align*}

Note that $\norm{ I^{k} u - u }_{L^2( \hat{\Omega}  ) } \leq  \norm{ I^{L}u - u }_{L^2( \hat{\Omega}  )}$ by the best approximation property of the
$L^2(\hat{\Omega} )$-projection.
Therefore, summing up over $K \in \hat{\mathcal{T}}_h$ yields

\[
h^{-2} \sum\limits_{K \in \hat{\mathcal{T}}_h }  \norm{ I^{k} u - u }_{L^2(K)}^2 = h^{-2} \norm{ I^{k} u - u }_{L^2(\hat{\Omega} )}^2 \leq h^{-2} \norm{ I^{L} u - u }_{L^2(\hat{\Omega} )}^2.
\]

Combining these results it follows that

\begin{align*}
\int\limits_{\tilde{I}} \sum\limits_{K \in \hat{\mathcal{T}}_h } \norm{\mathcal{L}^{k} u}_{ H^2(K) }^2 \; \dT
&\leq C \int\limits_{ \tilde{I}}  \sum\limits_{K \in \hat{\mathcal{T}}_h  } \left\{  \norm{ I^{L}u - u }_{H^2(K)}^2 + h^{-2} \norm{ I^{L}u - u }_{L^2(K)}^2  \right\} \dT   \\
& \leq C h^{2(s-2)} \norm{u}_{H^{0,s}(\tilde{Q} )}^2, \; s \in \{2,\ldots,k+1 \}.
\end{align*}

Using the result established above it is easy to deduce the interpolation results
required for analysis of the semi-discretization. Proof of <ref> for $ \Pi_h = \Pi_h^k$:  To define an interpolation operator into $\SemiDiscrSpace{k}$, we precompose $\pi^{k}$ defined on $\tilde{I} = (0,T)$ with a continuous (and degree independent) Sobolev extension operator $\mathcal{E}: H^s(\STdom) \rightarrow H^s(\hat{Q})$, see <cit.>. That is, let $\Pi_h^k u :=  (\pi^{k} \mathcal{E}u)|_{Q}$. Then by construction  $\Pi_h^k u \in  L^2(\tilde{I}) \otimes  V_{h}^k $. Thanks to <ref> we have for sufficiently regular $u$ that

\[
\partial_t \Pi_h^k u =   (\partial_t \pi^{k} \mathcal{E}u)|_{Q}.
= ( \pi^{k} \partial_t \mathcal{E}u)|_{Q} \in  L^2(\tilde{I}) \otimes  V_{h}^k.
\]

Hence, $\Pi_h^k u  \in H^1 \left( 0,T; V_h^k \right) =  \SemiDiscrSpace{k}$.
To apply the interpolation results established in <ref> to $\Pi_h^k$, we note that for any $D \in \{ \mathrm{id},\nabla,\partial_t \}$
we have

\begin{align*}
\norm{D ( \Pi_h^ku -u ) }_{Q} = \norm{D ( \pi^{k} \mathcal{E}u  - \mathcal{E} u ) }_{Q}
\leq  \norm{D ( \pi^{k} \mathcal{E}u  - \mathcal{E} u ) }_{ \hat{Q}}.
\end{align*}

Using the interpolation results from <ref> then allows to estimate the last expression by some power of $h$
times $C \norm{  \mathcal{E}u }_{ H^s( \hat{Q}) }$ which can be bounded by $C \norm{u }_{ H^s( Q) }$ thanks to continuity of the Sobolev extension operator. Using this argument the interpolation results in <ref> follow from <ref>.

§.§ Interpolation in space-time

To define an interpolation operator into the fully-discrete spaces $\FullyDiscrSpace{k}{q}$, some additional work is necessary. If we apply the results from <ref> on $\tilde{I} = I_n$ we obtain an interpolation operator $\pi^{k}_n$ into $L^2(I_n) \otimes  \hat{V}_{h}^k$ for $n=0,\ldots,N-1$. The space-time interpolation operator $\pi^{q,k}_n$ is then obtained by concatenating the purely spatial and temporal operators introduced in the previous two subsections. We define

\begin{equation}\label{eq:interp_time_slab}
\pi_{n}^{q,k}: H^{0,0}(\hat{Q}^n) \rightarrow \mathcal{P}^q(I_n) \otimes \hat{V}_h^k, \quad \pi^{q,k}_{n} = \pi^{q}_{n} \circ \pi^{k}_n
\end{equation}

and obtain the global operator by means of restriction to the time-slabs:

\begin{equation}\label{eq:Definition space-time interpolation operator_DG}
\Pi^{q,k}_{n}: H^{0,0}(\hat{Q}) \rightarrow \FullyDiscrSpaceHat{q}{k}, \quad \Pi^{q,k}_{n}u =  \pi^{q,k}_{n} u|_{ \hat{Q}^n }.
\end{equation}

An important ingredient for establishing the approximation results of these operators is that spatial derivatives commute with the temporal interpolation and vice versa (cf. <ref> and <ref>).

For $l_{q},l_{k} \in \mathbb{N}$ define $l_{\max} = \max{\{l_{q},l_{k} \}}$.
The following stability results and approximation error bounds hold:

* $ \norm{ \pi^{q,k}_{n}  u}_{L^2(\hat{Q}^{n})} \leq \norm{u}_{L^2(\hat{Q}^{n})} $ and $ \norm{ \pi^{q,k}_{n}  u}_{H^1(\hat{Q}^{n})} \leq C  \norm{u}_{H^1(\hat{Q}^{n})}  $.

* $\norm{u -  \pi^{q,k}_{n} u}_{\hat{Q}^{n}} \leq C \left( \Delta{t}^{l_{q}}  + h^{l_{k}}  \right) \norm{u}_{H^{l_{\max} }(\hat{Q}^{n})}$,
for $l_{k} \in \{1,\ldots,k + 1 \}$ and $l_{q} \in \{1,\ldots,q+1 \}$.

* $\norm{\partial_{t} (u - \pi^{q,k}_{n} u) }_{\hat{Q}^{n}} \leq C \left( \Delta{t}^{l_{q}-1}  + h^{l_{k}}\right) \norm{u}_{H^{l_{\max}+1 }(\hat{Q}^{n})}$, for $l_{k} \in \{1,\ldots,k +1 \}$ and $l_{q} \in \{2,\ldots,q+1 \}$.

* $\norm{\nabla{ (u -  \pi^{q,k}_{n} u) }}_{\hat{Q}^{n}} \leq C \left( \Delta{t}^{l_{q}}  + h^{l_{k}-1}  \right) \norm{u}_{H^{l_{\max}+1 }(\hat{Q}^{n})}$, for $l_{k} \in \{2,\ldots,k +1 \}$ and $l_{q} \in \{1,\ldots,q+1 \}$.

* For $u \in H^{l_{\max}+2 }(\hat{Q}^{n}) \cap C^{0}(I_{n},H^2(\hat{\Omega} ))$ we have

\[  \left( \int\limits_{I_n} \sum\limits_{K \in \hat{\mathcal{T}}_h } \norm{ u - \pi^{q,k}_{n} u }_{ H^2(K) }^2 \; \dT \right)^{1/2} \leq C \left( \Delta{t}^{l_{q}}  + h^{l_{k}-2}  \right) \norm{u}_{H^{l_{\max}+2 }(\hat{Q}^{n})}  \]

for $l_{k} \in \{2,\ldots,k +1 \}$ and $l_{q} \in \{1,\ldots,q+1 \}$.

We only show (a) and (e) as the other results have been established in <cit.>.

From <ref> (a) and <ref> (a) we obtain

\[ \norm{ \pi^{q,k}_{n}  u}_{\hat{Q}^{n}} = \norm{  \pi^{q}_{n}  \pi^{k}_n  u}_{\hat{Q}^{n}} \leq \norm{  \pi^{k}_n u }_{\hat{Q}^{n}} \leq  \norm{ u }_{\hat{Q}^{n}}.   \]

By using that spatial derivatives commute with the interpolation in time (see <ref>) and applying
the stability results of <ref> (a) and <ref> (b)
we get

\[
\norm{\nabla \pi^{q,k}_{n}   u}_{\hat{Q}^{n}} =  \norm{  \pi^{q}_{n}  \nabla \pi^{k}_n u}_{\hat{Q}^{n}} \leq \norm{ \nabla \pi^{k}_n u}_{\hat{Q}^{n}} \leq C \norm{ u}_{H^{0,1}(\hat{Q}^{n})}.
\]

From the stability of the temporal interpolation, <ref> (b), we have

\[  \norm{\partial_t \pi^{q,k}_{n} u}_{Q^{n}} = \norm{ \partial_t \pi^{q}_{n}  ( \pi^{k}_n u) }_{\hat{Q}^{n}} \leq C  \norm{  \pi^{k}_n u }_{ H^{1,0}(\hat{Q}^{n}) }.   \]

Since temporal derivatives commute with the spatial interpolation (as shown in <ref>), we have

\[
\norm{  \pi^{k}_n u }_{ H^{1,0}(\hat{Q}^{n}) }^2 = \norm{ \pi^{k}_n \partial_t u }_{\hat{Q}_n}^2 + \norm{ \pi^{k}_n u }_{\hat{Q}_n}^2 \leq C  \norm{ u}_{H^{1,0}(\hat{Q}^{n})}^2
\]

by means of  <ref> (a).

We utilize the splitting $u- \pi^{q,k}_{n}  u = \mathcal{L}^q u + \pi^{q}_{n} \mathcal{L}^{k} u$.

\begin{equation*}
\int\limits_{I_n} \sum\limits_{K \in \hat{\mathcal{T}}_h } \norm{ u - \pi^{q,k}_{n} u }_{ H^2(K) }^2
\leq C  \int\limits_{I_n} \sum\limits_{K \in \mathcal{T}_h^n } \left\{  \norm{ \mathcal{L}^q u }_{ H^2(K) }^2 + \norm{ \pi^{q} \mathcal{L}^{k} u }_{ H^2(K) }^2   \right\}.
\end{equation*}

Since spatial derivatives commute with $ \mathcal{L}^q $ (see <ref>), we obtain from <ref> (c):

\[
\int\limits_{I_n} \sum\limits_{K \in \mathcal{T}_h }  \norm{ \mathcal{L}^q u }_{ H^2(K) }^2 \; \dT \leq C \Delta{t}^{2 l_q} \norm{u}_{H^{l_q,2}(\hat{Q}^{n})}^2 \quad  l_q \in \{0,\ldots,q+1 \}.
\]

For the other term we use stability of the temporal interpolation (see the <ref> (a)) and
<ref> (e) to obtain

\begin{align*}
\int\limits_{I_n} \sum\limits_{K \in \hat{\mathcal{T}}_h } \norm{ \pi^{q}_n  \mathcal{L}^{k} u }_{ H^2(K) }^2 \; \dT
& \leq
\int\limits_{I_n} \sum\limits_{K \in \hat{\mathcal{T}}_h } \norm{  \mathcal{L}^{k}u }_{ H^2(K) }^2 \; \dT \\
& \leq
C h^{2(l_k-2)} \norm{u}_{H^{0,l_k}(\hat{Q}^{n})}^2
\end{align*}

for $l_k \in \{2,\ldots,k+1 \}$ which yields the desired result.

To derive a bound on the approximation error at discrete time levels from the previous results, we need an additional discrete inverse inequality in time:

Let $M$ be a bounded Lipschitz domain.
There exists a constant $C > 0$ such that

\begin{equation}\label{eq:Discrete inverse inequality in time}
\norm{u(\tau,\cdot) }_{M} \leq C \Delta{t}^{-1/2} \norm{u}_{I_n \times M }
\end{equation}

for all $u \in  \mathcal{P}^{q}(I_n) \otimes L^2(M)$ and $\tau \in \{ t_n,t_{n+1}\}$.

See <cit.>. In this reference the case $\tau = t_n$ has been shown, however, the proof for $\tau = t_{n+1}$ is completely analogous.

At the top and bottom of the time slice the following approximation bounds for the traces hold true:

* For $u \in H^{l_{q},0}(\hat{Q}^{n}) \cap H^{0,l_{k}}(\hat{Q}^{n})$ we have

\begin{equation}\label{eq: space-time interpolation error estimate at fixed time level}
\norm{(u - \pi^{q,k}_{n} u )( \tau, \cdot )}_{\hat{\Omega } } \leq C \left( \Delta{t}^{l_{q}-1/2} \norm{u}_{H^{l_{q},0}(\hat{Q}^{n})} + \Delta{t}^{-1/2} h^{l_{k}} \norm{u}_{H^{0,l_{k}}(\hat{Q}^{n})}  \right),
\end{equation}

with $\tau \in \{ t_n,t_{n+1}\}$ for $l_{k} \in \{1,\ldots,k_{s} +1 \}$ and $l_{q} \in \{1,\ldots,q+1 \}$,

* For $u \in H^{l_{q},0}(\hat{Q}^{n}) \cap H^{0,l_{k}}(\hat{Q}^{n})$ we have

\begin{align}
& \norm{\nabla (u - \pi^{q,k}_{n} u )( \tau, \cdot )}_{\hat{\Omega} } \label{eq:time_lvl_approx_grad}
\\
& \leq C \left( \Delta{t}^{l_{q}-1/2} \norm{\nabla u}_{H^{l_{q},0}(\hat{Q}^{n})} + \Delta{t}^{-1/2} h^{l_{k}-1} \norm{u}_{H^{0,l_{k}}(\hat{Q}^{n})}  \right), \nonumber
\end{align}

with $\tau \in \{ t_n,t_{n+1}\}$, for $l_{k} \in \{1,\ldots,k_{s} +1 \}$ and $l_{q} \in \{1,\ldots,q+1 \}$.

We will show (b). The proof of (a) is analogous and can be found in <cit.>.
Recall the splitting $u- \pi^{q,k}_{n}  u = \mathcal{L}^q u + \pi^{q}_{n} \mathcal{L}^{k} u$.
An application of the triangle inequality yields

\[
\norm{\nabla (u - \pi^{q,k}_{n}  u )( \tau, \cdot )}_{ \hat{\Omega} }
\leq \norm{ \nabla ( \mathcal{L}^q u  )( \tau, \cdot )}_{ \hat{\Omega} }
+ \norm{ \nabla ( \pi^{q}_{n} \mathcal{L}^{k} u )( \tau, \cdot )}_{ \hat{\Omega} }
\]

To treat the second term, we can apply the discrete inverse inequality in time, see the <ref>, which yields

\begin{align*}
& \norm{ \nabla ( \Pi^{q}_{n} \mathcal{L}^{k} u )( \tau, \cdot )}_{ \hat{\Omega}  }
= \norm{ \Pi^{q}_{n} \nabla \mathcal{L}^{k} u ( \tau, \cdot )}_{ \hat{\Omega} } \\
& \leq C (\Delta t)^{-1/2} \norm{ \nabla \mathcal{L}^{k} u }_{ \hat{Q}^n }
\leq C (\Delta t)^{-1/2} h^{ l_{k} -1 } \norm{u}_{H^{0,l_{k}}(\hat{Q}^{n})}
\end{align*}

for $l_{k} \in \{1,\ldots,k_{s} +1 \}$,  where the spatial approximation results from <ref> have been employed in the last step.
To treat the term $  \norm{ \nabla ( \mathcal{L}^q u  )( \tau, \cdot )}_{\hat{\Omega}} = \norm{  \mathcal{L}^q g  ( \tau, \cdot )}_{\hat{\Omega}}  $ with $g=\nabla u$,
we transform the interval $I_n$ to the reference interval $\tilde{I} = (0,1]$ and consider the problem there.
Denote the corresponding transformation, which only depends on the temporal variable, by $\Phi: \tilde{Q} \rightarrow \hat{Q}^{n}$ with $\tilde{Q} = \tilde{I} \times \hat{\Omega} $.
The transformed function and operator are given by $\tilde{g} = g \circ \Phi$ and $ \tilde{\mathcal{L}}^q = \Phi^{-1} \circ \mathcal{L}^q \circ \Phi $.
Using continuity of the time trace operator in $H^{1,0}(\tilde{Q})$ and the interpolation estimates from <ref> (c) for $ \mathcal{L}^q$ yields (denoting $\tilde{\tau} = \Phi^{-1}(\tau)$)

\begin{align*}
& \norm{  \mathcal{L}^q g  ( \tau, \cdot )}_{\hat{\Omega}  }^2
=  \norm{ \tilde{\mathcal{L}}^q  \tilde{g}  ( \tilde{\tau}, \cdot )}_{\hat{\Omega} }^2
\leq C \norm{ \tilde{\mathcal{L}}^q  \tilde{g} }_{ H^{1,0}(\tilde{Q}) }^2
= C \left(  \norm{ \tilde{\mathcal{L}}^q  \tilde{g} }_{ \tilde{Q} }^2 + \norm{ \partial_t \tilde{\mathcal{L}}^q  \tilde{g} }_{ \tilde{Q} }^2   \right) \\
& \leq C  \left( \frac{1}{ \Delta t }   \norm{ \mathcal{L}^q  g }_{ \hat{Q}^{n} }^2 + \Delta t \norm{ \partial_t \mathcal{L}^q  g }_{ \hat{Q}^{n} }^2   \right) \\
& \leq C \left( \frac{ (\Delta t)^{2  l_{q} } }{ \Delta t } \norm{g}_{ H^{l_{t},0}(\hat{Q}^n) }^2 + \Delta (\Delta t)^{ 2(l_{q} -1  )  } \norm{g}_{ H^{l_{t},0}(\hat{Q}^n) }^2  \right)
\end{align*}

for $l_{q} \in \{1,\ldots,q+1 \}$.
Combining these estimates yields the claim.

To transfer the results of <ref> and <ref> from $\hat{Q}^{n}$ to $Q^{n}$ we proceed analogous to the semi-discrete case. Proof of <ref> and <ref> for $ \Pi_h = \Pi_{h, \Delta t}^{k,q} $ As before let $\mathcal{E}: H^s(\STdom) \rightarrow H^s(\hat{Q})$ be a continuous Sobolev etension operator. We define $\Pi_{h, \Delta t}^{k,q}$ time-slab wise by

\[
\Pi_{h, \Delta t}^{k,q}u|_{Q^n} =   (\pi^{q,k}_{n} \mathcal{E}u)|_{Q^n} \in \mathcal{P}^q(I_n) \otimes V_h^k.
\]

Then for $D \in \{ \mathrm{id}, \partial_t, \nabla, \partial_{x_i} \partial_{x_j} \}$ we have

\begin{align*}
\norm{D ( \Pi_{h, \Delta t}^{k,q} u -u ) }_{Q^n} = \norm{D ( \pi^{q,k}_{n} \mathcal{E}u  - \mathcal{E} u ) }_{Q^n}
\leq  \norm{D ( \pi^{q,k}_{n} \mathcal{E}u  - \mathcal{E} u ) }_{ \hat{Q}^n}.
\end{align*}

Then using <ref> and continuity of the Sobolev extension operator this can be bounded by an appropriate power of $h$ (specified in the theorem) times $C \norm{u}_{H^s(Q^n)}$.
We proceed similarly for <ref>. With this argument we obtain <ref> from <ref> and <ref> from <ref>.

[1]
M. S. Alnaes, A. Logg, K. B. Ølgaard, M. E. Rognes, and G. N. Wells,
Unified Form Language: A domain-specific language for weak
formulations of partial differential equations, ACM Trans. Math. Software,
40 (2014), <https://doi.org/10.1145/2566630>.

[2]
I. A. Baratta, J. P. Dean, J. S. Dokken, M. Habera, J. S. Hale, C. N.
Richardson, M. E. Rognes, M. W. Scroggs, N. Sime, and G. N. Wells, DOLFINx: the next generation FEniCS problem solving environment.
preprint, 2023, <https://doi.org/10.5281/zenodo.10447666>.

[3]
C. Bardos, G. Lebeau, and J. Rauch, Un exemple d'utilisation des
notions de propagation pour le contrôle et la stabilisation de
problèmes hyperboliques, Rend. Sem. Mat. Univ. Politec. Torino,  (1988),
pp. 11–31 (1989).
Nonlinear hyperbolic equations in applied sciences.

[4]
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 (1992), pp. 1024–1065,

[5]
E. Burman, Stabilized Finite Element Methods for
Nonsymmetric, Noncoercive, and Ill-Posed Problems. Part I:
Elliptic Equations, SIAM J. Sci. Comput., 35 (2013), pp. A2752–A2780,

[6]
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 (2023),
pp. 2534–2557, <https://doi.org/10.1137/22M1508637>.

[7]
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 (2021),
pp. S969–S991, <https://doi.org/10.1051/m2an/2020062>.

[8]
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 (2023), p. 41,

[9]
E. Burman, A. Feizmohammadi, and L. Oksanen, A finite element data
assimilation method for the wave equation, Math. Comp., 89 (2020),
pp. 1681–1709, <https://doi.org/10.1090/mcom/3508>.

[10]
E. Burman and J. Preuss, Reproduction material: Unique continuation
for the wave equation based on a discontinuous Galerkin time
discretization, May 2024, <https://doi.org/10.5281/zenodo.11102764>.

[11]
D. Chapelle, N. Cîndea, M. De Buhan, and P. Moireau, Exponential
convergence of an observer based on partial field measurements for the wave
equation, Math. Probl. Eng.,  (2012), pp. Art. ID 581053, 12,

[12]
D. Chapelle, N. Cîndea, and P. Moireau, Improving convergence in
numerical analysis using observers—the wave-like equation case, Math.
Models Methods Appl. Sci., 22 (2012), pp. 1250040, 35,

[13]
N. Cîndea and A. Münch, Inverse problems for linear
hyperbolic equations using mixed formulations, Inverse Problems, 31 (2015),
p. 075001, <https://doi.org/10.1088/0266-5611/31/7/075001>.

[14]
W. Dahmen, H. Monsuur, and R. Stevenson, Least squares solvers for
ill-posed pdes that are conditionally stable, ESAIM Math. Model. Numer.
Anal., 57 (2023), pp. 2227–2255, <https://doi.org/10.1051/m2an/2023050>.

[15]
D. A. French, A space-time finite element method for the wave
equation, Comput. Methods Appl. Mech. Engrg., 107 (1993), pp. 145–157,

[16]
T. J. R. Hughes and G. M. Hulbert, Space-time finite element methods
for elastodynamics: formulations and error estimates, Comput. Methods Appl.
Mech. Eng., 66 (1988), p. 339–363,

[17]
C. Johnson, Discontinuous Galerkin finite element methods for
second order hyperbolic problems, Comput. Methods Appl. Mech. Engrg., 107
(1993), pp. 117–129, <https://doi.org/10.1016/0045-7825(93)90170-3>.

[18]
C. Lehrenfeld, On a Space-Time Extended Finite Element Method for
the Solution of a Class of Two-Phase Mass Transport Problems, PhD thesis,
RWTH Aachen, February 2015,

[19]
S. Mishra and R. Molinaro, Estimates on the generalization error of
physics-informed neural networks for approximating a class of inverse
problems for PDEs, IMA J. Numer. Anal., 42 (2021), pp. 981–1022,

[20]
M. Nechita, Solving ill-posed Helmholtz problems with
physics-informed neural networks, J. Numer. Anal. Approx. Theory, 52 (2023),
p. 90–101, <https://doi.org/10.33993/jnaat521-1305>.

[21]
J. Preuß, Higher order unfitted isoparametric space-time FEM on
moving domains, master's thesis, NAM, University of Göttingen, February
2018, <https://doi.org/10.25625/UACWXS>,
<https://data.goettingen-research-online.de/file.xhtml?persistentId=doi:10.25625/UACWXS/HKWJE5   version=1.0>.

[22]
M. W. Scroggs, I. A. Baratta, C. N. Richardson, and G. N. Wells, Basix: a runtime finite element basis evaluation library, J. Open Source
Softw., 7 (2022), p. 3982, <https://doi.org/10.21105/joss.03982>.

[23]
M. W. Scroggs, J. S. Dokken, C. N. Richardson, and G. N. Wells, Construction of arbitrary order finite element degree-of-freedom maps on
polygonal and polyhedral cell meshes, ACM Trans. Math. Software, 48 (2022),
pp. 18:1–18:23, <https://doi.org/10.1145/3524456>.

[24]
E. M. Stein, Singular Integrals and Differentiability
Properties of Functions (PMS-30), Volume 30, Princeton University
Press, Princeton, 1971, <https://doi.org/doi:10.1515/9781400883882>.

[25]
V. Thomée, Galerkin finite element methods for parabolic
problems, vol. 25, Springer Science & Business Media, 2007.