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

Erik Burman, Guillaume Delay, Alexandre Ern

► To cite this version:

Erik Burman, Guillaume Delay, Alexandre Ern. The unique continuation problem for the heat equation discretized with a high-order space-time nonconforming method. SIAM Journal on Numerical Analysis, 2023, 61 (5), pp.2534-2557. 10.1137/22M1508637. hal-03720960v4

HAL Id: hal-03720960

https://hal.science/hal-03720960v4

Submitted on 16 Oct 2023

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

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

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

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# The unique continuation problem for the heat equation discretized with a high-order space-time nonconforming method*

Erik Burman†

Guillaume Delay‡

Alexandre Ern§

October 16, 2023

## Abstract

We are interested in solving the unique continuation problem for the heat equation, i.e., we want to reconstruct the solution of the heat equation in a target space-time subdomain given its (noised) value in a subset of the computational domain. Both initial and boundary data can be unknown. We discretize this problem using a space-time discontinuous Galerkin method (including hybrid variables in space) and look for the solution that minimizes a discrete Lagrangian. We establish discrete inf-sup stability and bound the consistency error, leading to a priori estimates on the residual. Owing to the ill-posed nature of the problem, an additional estimate on the residual dual norm is needed to prove the convergence of the discrete solution to the exact solution in the energy norm in the target space-time subdomain. This is achieved by combining the above results with a conditional stability estimate at the continuous level. The rate of convergence depends on the conditional stability, the approximation order in space and in time, and the size of the perturbations in data. Quite importantly, the weight of the regularization term depends on the time step and the mesh size, and we show how to choose it to preserve the best possible decay rates on the error. Finally, we run numerical simulations to assess the performance of the method in practice.

## 1 Introduction

In the present work, we are interested in solving numerically a data assimilation problem subject to the heat equation. In this problem, neither the boundary conditions nor the initial data are known. In order to compensate for the lack of initial and boundary data, we use the knowledge of the solution in a subdomain. We also investigate the influence of noise on this additional datum. Specifically, we consider a bounded Lipschitz domain $\Omega \subset \mathbb{R}^d$, $d \in \{1, 2, 3\}$, a subset $\omega \subset \Omega$, and a time interval $J := (0, T_f)$ with final time

<footnote>*Draft version, October 16, 2023</footnote>

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

<footnote>‡Sorbonne Université, CNRS, Université Paris Cité, LJLL, F-75005 Paris, France.</footnote>

<footnote>§CERMICS, École des Ponts, 77455 Marne-la-Vallée cedex 2, and INRIA, Paris, France.</footnote>
$T_f > 0$. Our goal is to approximate the function $u: J \times \Omega \to \mathbb{R}$ that satisfies

$$L(u) := \frac{\partial u}{\partial t} - \Delta u = f \quad \text{in } J \times \Omega, \qquad (1)$$

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

where $f \in L^2(J \times \Omega)$ is a given source term and $g \in H^1(J; H^1(\bar{\omega})') \cap L^2(J; H^1(\bar{\omega}))$ is the restriction to $J \times \bar{\omega}$ of a solution to the heat equation in $J \times \Omega$. The model problem (1)-(2) is ill-posed. Indeed, whenever a solution exists (this is the case whenever $g$ satisfies the heat equation (1) in $J \times \bar{\omega}$), it is unique, but there is no a priori estimate on the solution in the usual Hadamard sense. Another difficulty is that we want to consider perturbed data $g_\delta$ instead of $g$ to account for some noise in the measurements.

Although we do not have usual stability estimates, so-called conditional stability estimates are available. This type of estimate will play a key role in our error analysis. Conditional stability estimates essentially allow one to control the norm of a function in a target subdomain by means of weaker norms in a larger domain. Our analysis hinges on the following result.

**Lemma 1 (Conditional stability estimate).** Let $B$ be a connected subset of $\Omega$ such that $\bar{B} \subset \Omega$. Let $0 < T_1 < T_2 < T_f$. There exist $C_{\text{stb}} > 0$ and $\alpha \in (0, 1]$ such that for all $v \in H^1(J; H^{-1}(\Omega)) \cap L^2(J; H^1(\Omega))$, we have

$$\|v\|_{L^2(T_1, T_2; H^1(B))} \le C_{\text{stb}} (\|v\|_{L^2(J; L^2(\bar{\omega}))} + \|L(v)\|_{L^2(J; H^{-1}(\Omega))})^\alpha \\ \times (\|v\|_{L^2(J; L^2(\Omega))} + \|L(v)\|_{L^2(J; H^{-1}(\Omega))})^{1-\alpha}. \qquad (3)$$

Lemma 1 with the additional assumption $\bar{\omega} \subset B$ corresponds to Theorem 1 in [8]. This assumption can be removed by using the arguments of the proof of Theorem 1.1 from [23]; see also [24, 17, 1]. The conditional stability estimate from Lemma 1 will be used to prove the convergence of our approximation method. The constant $\alpha$ therein has an influence on the convergence rate. For instance, if $h$ is the mesh size and $\tau$ the time step, if polynomials of degree $k \ge 1$ (resp., $\ell \ge 0$) are used for the space (resp., time) discretization, and provided the noise in the measurements is small enough, our main result establishes an error bound with decay rate $O((h^k + \tau^{\ell+\frac{1}{2}})^\alpha)$ in the target subdomain $B$ (see Theorem 9). We also mention that Lemma 1 admits some interesting extensions. First, in the case of known boundary conditions (i.e., if we add the assumption $v = 0$ on $\partial\Omega$), then Lemma 1 remains valid with $T_2 = T_f$, $B = \Omega$, and $\alpha = 1$ (see Theorem 2 from [8]). This setting will be considered in our numerical experiments as well. Moreover, in the case of known initial conditions (i.e., if we add the assumption $v(0, \cdot) = v_0$ in $\Omega$ for some $v_0 \in L^2(\Omega)$), then Lemma 1 remains valid with $T_1 = 0$ provided $\|v_0\|_{\Omega}$ is added to $\|L(v)\|_{L^2(J; H^{-1}(\Omega))}$ in (3) [13]. Other stability estimates are available in the literature, for instance using Cauchy data on the boundary instead of interior data [22]. For an overview on analysis techniques of unique continuation for parabolic equations, we refer the reader to [27].

Since the problem (1)-(2) is ill-posed, a regularization must be considered to devise a reasonable approximation method. The usual approach is to regularize the continuous problem before embarking on any discretization method. Several regularization methods
are available, for instance the quasi-reversibility method [3, 4, 16] or the Tikhonov regularization [20]. These methods have already been applied to solve the data assimilation problem subject to the heat equation in the one-dimensional case; see [26, 2]. Other numerical approaches have also been proposed and analyzed for the one-dimensional heat equation, see [21, 28], but, to the best of our knowledge, none of the above-mentioned references provides an error analysis balancing the stability of the regularization method, the conditional stability estimate, and the approximation order of the discretization in space and in time.

An alternative approach is to first discretize the ill-posed problem and then regularize it at the discrete level. The main advantage of discretizing first and then regularizing is that it makes it possible to design regularization terms that allow for a rigorous numerical analysis, leading to error estimates with rates that match the best possible rates deduced from the conditional stability estimates. The discretize-then-regularize approach has already been considered for instance in [5, 6, 9] for the stationary version of the present problem. The analysis was extended to the data assimilation problem subject to the Helmholtz equation using a high-order discretization in [11]. Data assimilation subject to non-stationary problems was considered in [12] for the wave equation and in [8, 10] for the heat equation with known Dirichlet boundary data for the numerical tests. In particular, [8] considers the semi-discretization in space using lowest-order Lagrange finite elements, discussing several stability situations, whereas [10] deals with the full discretization of the problem using the backward Euler scheme in time and still lowest-order Lagrange finite elements in space. In the latter reference, only the case of unknown initial data, but known boundary data, is considered. In [15], a method minimizing the dual norm of the PDE residual augmented by the least-squares error on data fitting was proposed, together with error estimates for the reconstruction problem. Also in this latter work, the boundary data are assumed to be known.

In the present work, we use higher-order methods for the discretization in space and in time. Specifically, we employ a discontinuous Galerkin (dG) method in time and a hybridized dG method in space (recall that such methods attach discrete unknowns to the mesh cells and to the mesh faces). The use of a dG method in time is the natural way to extend the backward Euler scheme to higher order. Furthermore, the use of a hybridized dG method in space reduces the number of stabilization terms that are needed with respect to standard finite elements (irrespective of the order of approximation). Indeed, as we shall see below (see the inf-sup condition from Lemma 3), the hybridized dG method naturally gives a control on the residuals without the need of any stabilization (which is needed for standard finite elements). Notice that, at the algebraic level, the discrete unknowns are globally coupled in time, as is anyway the case for all the methods to solve data assimilation problems.

To sum up, the added value of the present work is that we consider the data assimilation problem with unknown boundary conditions (and initial conditions) and that we devise a higher-order method in space and in time to discretize the above ill-posed problem. The salient feature of our approximation method is that the regularization term matches the decay rates of the approximation method. We also emphasize that we rigorously derive an a priori error estimate that decays at the best possible rate in view of the approximation capacities of the discretization and of the conditional stability of the continuous problem.
Finally, we mention that the analysis entails various subtleties because some of the usual arguments in the analysis of time dG methods cannot be applied here. The reason for this is that, because of the ill-posed nature of the problem, we cannot invoke a plain coercivity argument that gives control in the $L^2(J; H_0^1(\Omega))$-norm. As a consequence, we need to add a stabilization term controlling the time jumps of the discrete solution (and of its gradient) at all the discrete time nodes. Moreover, in addition to the main analysis arguments used for well-posed problems (inf-sup stability, consistency), we need an additional, non-trivial estimate on a suitable dual residual norm. This estimate is indeed key to invoke the conditional stability estimate from Lemma 1 and prove an error estimate in the energy norm.

The rest of this work is organized as follows. In Section 2, we present the high-order discretization method in space and in time. In Section 3, we perform the error analysis, which combines arguments relevant to well-posed problems (Sections 3.2 to 3.4) with arguments specific to ill-posed problems (Section 3.5). These steps prepare the stage for our main results, which are the error estimate (Theorem 7) and the devising of the regularization to achieve optimal error decay rates (Theorem 9). Finally, in Section 4, we present some numerical results that corroborate our theoretical results. Moreover, we illustrate the benefits of high-order discretization methods.

# 2 Problem discretization

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

## 2.1 Time discretization

In what follows, for two integer numbers $p \le q$, we use the notation $\{p:q\} := \{m \in \mathbb{N}, p \le m \le q\}$. We discretize the time interval $J$ using a uniform time step $\tau := T_f/N$ (for simplicity), where $N \in \mathbb{N}^*$ is the number of time steps. We then have $0 = t_0 < t_1 < \dots < t_N = T_f$ where $t_n := n\tau$ for all $n \in \{0:N\}$. We define $I_n := (t_{n-1}, t_n)$ for all $n \in \{1:N\}$.

For a piecewise smooth function $\Phi^\tau$ defined on the above partition of $J$, we use the shorthand notation $\Phi^n := \Phi^\tau|_{I_n}$ for all $n \in \{1:N\}$. We define $\Phi^\tau(t_n^-) := \Phi^n(t_n)$, $\Phi^\tau(t_n^+) := \Phi^{n+1}(t_n)$ and $[\Phi^\tau]^n := \Phi^\tau(t_n^+) - \Phi^\tau(t_n^-)$. For convenience, we also set $[\Phi^\tau]^0 := 0$; this convention is motivated by the fact that there is no initial condition to enforce in the present problem. Finally, $\partial_t^\tau \Phi^\tau$ denotes the broken time derivative of $\Phi^\tau$ such that $(\partial_t^\tau \Phi^\tau)|_{I_n} := \partial_t \Phi^n$ for all $n \in \{1:N\}$.

## 2.2 Space discretization

Let $(\mathcal{T}_h)_{h>0}$ be a family of matching meshes of $\Omega$. In principle, the meshes can have cells that are polyhedra with planar faces in $\mathbb{R}^d$, and hanging nodes are also possible. However, the analysis below requires the mesh to be such that the underlying discontinuous polynomial approximation space has a global $H^1$-conforming subspace with optimal approximation properties. For simplicity, we will therefore restrict the discussion to meshes composed of simplices (one can also readily consider meshes composed of cuboids). The mesh cells are conventionally taken to be open subsets of $\mathbb{R}^d$, and $n_T$ denotes the unit
outward normal to the generic mesh cell $T \in \mathcal{T}_h$. For a subset $S \subset \mathbb{R}^d$, $h_S$ denotes the diameter of $S$, and for a mesh $\mathcal{T}_h$, the index $h$ refers to the maximal diameter of the mesh cells.

The mesh faces are collected in the set $\mathcal{F}_h$ which is split into the set of the mesh interfaces, $\mathcal{F}_h^{\text{int}}$, and the set of the mesh boundary faces, $\mathcal{F}_h^{\partial}$. Any mesh interface $F \in \mathcal{F}_h^{\text{int}}$ is oriented by a fixed unit normal vector $\mathbf{n}_F$. Moreover, for a piecewise smooth function $v$ and any mesh interface $F \in \mathcal{F}_h^{\text{int}}$, $[v]_F$ denotes the jump of $v$ across $F$ in the direction of $\mathbf{n}_F$. We also use the broken gradient and Laplacian operators, $\nabla_{\mathcal{T}}$ and $\Delta_{\mathcal{T}}$, which are defined such that $(\nabla_{\mathcal{T}}v)|_T := \nabla(v|_T)$ and $(\Delta_{\mathcal{T}}v)|_T := \Delta(v|_T)$ for all $T \in \mathcal{T}_h$.

To avoid technicalities, we assume henceforth that the mesh family $(\mathcal{T}_h)_{h>0}$ is quasi-uniform. Therefore, we will use $h$ to measure the diameter of any mesh cell or any mesh face. Moreover, we assume that all the meshes are fitted to the subset $\varpi$.

## 2.3 Discrete spaces and bilinear forms

Let $k \ge 1$ be the polynomial degree of the hybridized dG method in space and let $\ell \ge 0$ be the polynomial degree of the dG method in time. We denote by $\mathbb{P}^k(S)$ the set of polynomials of total degree at most $k$ on the subset $S \subseteq \Omega$. Moreover, for a linear space $U$ composed of functions defined on $\Omega$, we denote by $\mathbb{P}^\ell(I; U)$ the set of $U$-valued polynomials of degree at most $\ell$ on $I \subseteq \bar{J} = [0, T_f]$.

The discrete unknowns in space are piecewise polynomials of degree $k$ attached to the mesh cells and of the same degree $k$ attached to the mesh faces. We define the discrete spaces

$$ \hat{U}_h^k := U_T^k \times U_F^k, \quad U_T^k := \chi_{T \in \mathcal{T}_h} \mathbb{P}^k(T), \quad U_F^k := \chi_{F \in \mathcal{F}_h} \mathbb{P}^k(F). \quad (4) $$

For a generic pair $\hat{v}_h \in \hat{U}_h^k$, we write $\hat{v}_h := (v_{\mathcal{T}}, v_{\mathcal{F}})$ with $v_{\mathcal{T}} := (v_T)_{T \in \mathcal{T}_h} \in U_T^k$ and $v_{\mathcal{F}} := (v_F)_{F \in \mathcal{F}_h} \in U_F^k$. We denote by $\hat{U}_{h0}^k$ the linear subspace of $\hat{U}_h^k$ in which all the degrees of freedom attached to the mesh boundary faces are null. For a generic pair $\hat{v}_h \in \hat{U}_h^k$, its degrees of freedom associated with a generic mesh cell $T \in \mathcal{T}_h$ are denoted by

$$ \hat{v}_T := (v_T, v_{\partial T} := (v_F)_{F \in \mathcal{F}_{\partial T}}) \in \hat{U}_T^k := \mathbb{P}^k(T) \times \mathbb{P}^k(\mathcal{F}_{\partial T}), \quad (5) $$

where $\mathbb{P}^k(\mathcal{F}_{\partial T}) := \chi_{F \in \mathcal{F}_{\partial T}} \mathbb{P}^k(F)$ and $\mathcal{F}_{\partial T} := \{F \in \mathcal{F}_h \mid F \subset \partial T\}$ collects the mesh faces composing the boundary of $T$. We introduce the space-time discrete spaces

$$ \tilde{U}_h^\tau := \{\hat{v}_h^\tau \in L^2(J; \hat{U}_h^k) \mid \hat{v}_h^\tau|_{I_n} \in \mathbb{P}^\ell(I_n; \hat{U}_h^k), \forall n \in \{1:N\}\}, \quad (6) $$

$$ \tilde{U}_{h0}^\tau := \{\hat{v}_h^\tau \in \tilde{U}_h^\tau \mid \hat{v}_h^\tau|_{I_n} \in \mathbb{P}^\ell(I_n; \tilde{U}_{h0}^k), \forall n \in \{1:N\}\}. \quad (7) $$

Notice that we have $\tilde{U}_h^\tau = U_T^\tau \times U_F^\tau$ with

$$ U_T^\tau := \{v_T^\tau \in L^2(J; U_T^k) \mid v_T^n := v_T^\tau|_{I_n} \in \mathbb{P}^\ell(I_n; U_T^k), \forall n \in \{1:N\}\}, \quad (8) $$

$$ U_F^\tau := \{v_F^\tau \in L^2(J; U_F^k) \mid v_F^n := v_F^\tau|_{I_n} \in \mathbb{P}^\ell(I_n; U_F^k), \forall n \in \{1:N\}\}. \quad (9) $$

For a generic function $\hat{v}_h^\tau \in \tilde{U}_h^\tau$, we write $\hat{v}_h^\tau := (v_T^\tau, v_F^\tau)$ with $v_T^\tau \in U_T^\tau$ and $v_F^\tau \in U_F^\tau$.

We can now introduce the various bilinear forms needed to formulate the discrete problem. Let $\hat{v}_h^\tau, \hat{w}_h^\tau$ be generic functions in $\tilde{U}_h^\tau$ (primal variables) and let $\hat{\zeta}_h^\tau, \hat{\eta}_h^\tau$ be
generic functions in $\tilde{U}_{h0}^{\tau}$ (dual variables). We use the subscript $h$ to indicate the bilinear forms related to the space discretization, and the superscript $\tau$ to indicate those related to the time discretization. The two bilinear forms associated with the discretization of the heat equation are

$$a_h(\hat{v}_h^\tau, \hat{\eta}_h^\tau) := \sum_{T \in \mathcal{T}_h} \left\{ \begin{aligned} & (\nabla v_T^\tau, \nabla \eta_T^\tau)_{J \times T} - (\nabla v_T^\tau \cdot \mathbf{n}_T, \eta_T^\tau - \eta_{\partial T}^\tau)_{J \times \partial T} \\ & - (v_T^\tau - v_{\partial T}^\tau, \nabla \eta_T^\tau \cdot \mathbf{n}_T)_{J \times \partial T} \end{aligned} \right\}, \quad (10)$$

$$b^{\tau}(v_{\tau}^{\tau}, \eta_{\tau}^{\tau}) := \sum_{n \in \{1:N\}} \left\{ (\partial_t v_{\tau}^{n}, \eta_{\tau}^{n})_{I_n \times \Omega} + (\llbracket v_{\tau}^{\tau} \rrbracket^{n-1}, \eta_{\tau}^{n}(t_{n-1}^{+}))_{\Omega} \right\}. \quad (11)$$

The stabilization bilinear forms read as follows:

$$s_h^\tau(\hat{v}_h^\tau, \hat{w}_h^\tau) := d_h(\hat{v}_h^\tau, \hat{w}_h^\tau) + d^\tau(v_T^\tau, w_T^\tau), \quad (12)$$

$$\sigma_h(\hat{\zeta}_h^\tau, \hat{\eta}_h^\tau) := (\nabla \tau \zeta_T^\tau, \nabla \tau \eta_T^\tau)_{J \times \Omega} + d_h(\hat{\zeta}_h^\tau, \hat{\eta}_h^\tau), \quad (13)$$

with

$$d_h(\hat{v}_h^\tau, \hat{w}_h^\tau) := \sum_{T \in \mathcal{T}_h} h^{-1}(v_T^\tau - v_{\partial T}^\tau, w_T^\tau - w_{\partial T}^\tau)_{J \times \partial T}, \quad (14)$$

$$d^{\tau}(v_{\tau}^{\tau}, w_{\tau}^{\tau}) := \sum_{n \in \{1:N-1\}} \left\{ (\llbracket v_{\tau}^{\tau} \rrbracket^{n}, \llbracket w_{\tau}^{\tau} \rrbracket^{n})_{\Omega} + (\llbracket \nabla_{\tau} v_{\tau}^{\tau} \rrbracket^{n}, \llbracket \nabla_{\tau} w_{\tau}^{\tau} \rrbracket^{n})_{\Omega} \right\}. \quad (15)$$

(Notice that the stabilization bilinear form for the dual variable is only related to the space discretization.) Finally, the bilinear forms associated with the regularization and the measurements are

$$t_h^\tau(v_T^\tau, w_T^\tau) := \gamma c(\tau, h)^2 (v_T^\tau, w_T^\tau)_{J \times \Omega}, \quad (16)$$

$$m_{\varpi}(v_T^\tau, w_T^\tau) := (v_T^\tau, w_T^\tau)_{J \times \varpi}, \quad (17)$$

where $\gamma > 0$ and the value of $c(\tau, h)$ will result from the error analysis (we shall obtain $c(\tau, h) = h^k + \tau^{\ell+\frac{1}{2}}$, see (56)). For later use, we define

$$\|w_T^\tau\|_{J \times \varpi}^2 := m_{\varpi}(w_T^\tau, w_T^\tau). \quad (18)$$

## 2.4 Lagrangian and discrete problem

We want to find the saddle-point of the Lagrangian defined for all $(\hat{v}_h^\tau, \hat{\zeta}_h^\tau) \in \tilde{U}_h^\tau \times \tilde{U}_{h0}^\tau$ by

$$\mathcal{L}_h^\tau(\hat{v}_h^\tau, \hat{\zeta}_h^\tau) := \frac{1}{2}\|v_T^\tau - g_\delta\|_{J \times \varpi}^2 + \frac{1}{2}t_h^\tau(v_T^\tau, v_T^\tau) + \frac{1}{2}s_h^\tau(\hat{v}_h^\tau, \hat{v}_h^\tau) - \frac{1}{2}\sigma_h(\hat{\zeta}_h^\tau, \hat{\zeta}_h^\tau) \\ + a_h(\hat{v}_h^\tau, \hat{\zeta}_h^\tau) + b^\tau(v_T^\tau, \zeta_T^\tau) - (f, \zeta_T^\tau)_{J \times \Omega}, \quad (19)$$

where $g_\delta := g + \delta$ denotes the available perturbed measurement of $g$. Notice that there is no boundary condition on the primal variable, whereas there is a Dirichlet boundary condition on the dual variable. The motivation for looking for a saddle-point of the discrete Lagrangian is to minimize the discrepancy with respect to data under the constraint of
the partial differential equation. This idea is classical, but the key point here is that the Lagrangian is defined at the discrete level. This makes it simpler to enhance the convexity/concavity of the discrete Lagrangian by suitable consistent stabilization terms; see [9] for further discussion.

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

$$m_{\varpi}(u_{\tau}^{\tau}, w_{\tau}^{\tau}) + t_{h}^{\tau}(u_{\tau}^{\tau}, w_{\tau}^{\tau}) + s_{h}^{\tau}(\hat{u}_{h}^{\tau}, \hat{w}_{h}^{\tau}) + a_{h}(\hat{w}_{h}^{\tau}, \hat{\xi}_{h}^{\tau}) + b^{\tau}(w_{\tau}^{\tau}, \xi_{\tau}^{\tau}) = m_{\varpi}(g_{\delta}, w_{\tau}^{\tau}), \quad (20)$$

$$a_h(\hat{u}_h^\tau, \hat{\eta}_h^\tau) + b^\tau(u_\tau^\tau, \eta_\tau^\tau) - \sigma_h(\hat{\xi}_h^\tau, \hat{\eta}_h^\tau) = (f, \eta_\tau^\tau)_{J \times \Omega}, \quad (21)$$

where the first equation holds for all $\hat{w}_h^\tau \in \tilde{U}_h^\tau$ and the second for all $\hat{\eta}_h^\tau \in \tilde{U}_{h0}^\tau$. For all $(\hat{v}_h^\tau, \hat{\zeta}_h^\tau)$ and $(\hat{w}_h^\tau, \hat{\eta}_h^\tau)$ in $\tilde{U}_h^\tau \times \tilde{U}_{h0}^\tau$, we define the bilinear form

$$A_h^\tau((\hat{v}_h^\tau, \hat{\zeta}_h^\tau), (\hat{w}_h^\tau, \hat{\eta}_h^\tau)) := m_\varpi(v_\tau^\tau, w_\tau^\tau) + t_h^\tau(v_\tau^\tau, w_\tau^\tau) + s_h^\tau(\hat{v}_h^\tau, \hat{w}_h^\tau) + a_h(\hat{w}_h^\tau, \hat{\zeta}_h^\tau) \\ + b^\tau(w_\tau^\tau, \zeta_\tau^\tau) + b^\tau(v_\tau^\tau, \eta_\tau^\tau) + a_h(\hat{v}_h^\tau, \hat{\eta}_h^\tau) - \sigma_h(\hat{\zeta}_h^\tau, \hat{\eta}_h^\tau). \quad (22)$$

The discrete problem (20)-(21) can be rewritten as follows: Find $(\hat{u}_h^\tau, \hat{\xi}_h^\tau) \in \tilde{U}_h^\tau \times \tilde{U}_{h0}^\tau$ such that

$$A_h^\tau((\hat{u}_h^\tau, \hat{\xi}_h^\tau), (\hat{w}_h^\tau, \hat{\eta}_h^\tau)) = m_\varpi(g_\delta, w_\tau^\tau) + (f, \eta_\tau^\tau)_{J \times \Omega}, \quad \forall (\hat{w}_h^\tau, \hat{\eta}_h^\tau) \in \tilde{U}_h^\tau \times \tilde{U}_{h0}^\tau.$$

# 3 Analysis

This section is organized as follows. We first introduce a time reconstruction operator to rewrite the bilinear form $b^\tau$. This operator is classical in the context of dG methods in time; see, e.g., [19, Section 69.2.3] or [25, Section 2.3] and the references therein. Then, we study the stability properties of $A_h^\tau$ in a suitable residual norm, we introduce interpolation operators in space and in time, and we bound the consistency error. The second important step, specific to ill-posed problems, is to bound a suitable dual residual norm. We combine these bounds with the abstract conditional estimate from Lemma 1 to derive error estimates in the target subdomain B introduced therein. Finally, we state the approximation properties of the interpolation operator, we tune the size of the stabilization parameter and we establish the convergence rates for the method.

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

## 3.1 Time reconstruction operator

For all $v_\tau^\tau \in U_\tau^\tau$, its time reconstruction $R^\tau(v_\tau^\tau) \in C^0(J; U_\tau^k)$ is defined such that $R^\tau(v_\tau^\tau)(t_0^+) := v_\tau^\tau(t_0^+)$ and such that, for all $n \in \{1:N\}$, $R^n(v_\tau^\tau) := R^\tau(v_\tau^\tau)|_{I_n} \in \mathbb{P}^{\ell+1}(I_n; U_\tau^k)$ satisfies

$$(\partial_t R^n(v_\tau^\tau), q_\tau^n)_{I_n \times \Omega} := (\partial_t v_\tau^n, q_\tau^n)_{I_n \times \Omega} + (\llbracket v_\tau^\tau \rrbracket^{n-1}, q_\tau^n(t_{n-1}^+))_{\Omega}, \quad (23)$$
for all $q_T^n \in \mathbb{P}^l(I_n; U_T^k)$. Since the function $R^\tau(v_T^\tau)$ is continuous in time, its time derivative is well defined in $L^2(J; U_T^k)$. The main consequence of (23) is that the bilinear form $b^\tau$ can be rewritten as

$$b^\tau(v_T^\tau, \eta_T^\tau) = (\partial_t R^\tau(v_T^\tau), \eta_T^\tau)_{J \times \Omega}, \quad \forall (v_T^\tau, \eta_T^\tau) \in U_T^\tau \times U_T^\tau. \quad (24)$$

Moreover, it is well-known that the reconstruction operator can be rewritten as

$$R^n(v_T^\tau)(t, \boldsymbol{x}) = v_T^n(t, \boldsymbol{x}) - \left[\left[v_T\right]\right]^{n-1}(\boldsymbol{x}) \frac{(-1)^\ell}{2} (L_\ell - L_{\ell+1}) \circ T_n^{-1}(t), \quad (25)$$

for all $(t, \boldsymbol{x}) \in I_n \times \Omega$, where $L_\ell$ is the Legendre polynomial of degree $\ell$ defined on $(-1, 1)$ and $T_n$ is the affine mapping from $(-1, 1)$ to $I_n$. A consequence of (25) is that

$$R^n(v_T^\tau)(t_n^-) = v_T^n(t_n^-), \quad \forall n \in \{1:N\}. \quad (26)$$

The following stability properties are useful for the present analysis (the proof is outlined in Section 3.8.1).

**Lemma 2 (Stability properties of $R^\tau$).** For all $v_T^\tau \in U_T^\tau$, we have

$$\|\partial_t^\tau (R^\tau(v_T^\tau) - v_T^\tau)\|_{J \times \Omega} \lesssim \tau^{-\frac{1}{2}} d^\tau(v_T^\tau, v_T^\tau)^{\frac{1}{2}}, \quad (27)$$

$$\|R^\tau(v_T^\tau) - v_T^\tau\|_{J \times \Omega} + \|\nabla \tau (R^\tau(v_T^\tau) - v_T^\tau)\|_{J \times \Omega} \lesssim \tau^{\frac{1}{2}} d^\tau(v_T^\tau, v_T^\tau)^{\frac{1}{2}}. \quad (28)$$

Moreover, for all $\hat{v}_h^\tau = (v_T^\tau, v_F^\tau) \in \tilde{U}_h^\tau$, we have

$$\|h^{-\frac{1}{2}}[R^\tau(v_T^\tau)]_{\mathcal{F}_h^{\text{int}}}\|_{J \times \mathcal{F}_h^{\text{int}}} := \left( \sum_{F \in \mathcal{F}_h^{\text{int}}} h^{-1} \| [R^\tau(v_T^\tau)]_F \|_{J \times F}^2 \right)^{\frac{1}{2}} \lesssim d_h(\hat{v}_h^\tau, \hat{v}_h^\tau)^{\frac{1}{2}}. \quad (29)$$

## 3.2 Residual stability

For all $(\hat{v}_h^\tau, \hat{\zeta}_h^\tau) \in \tilde{U}_h^\tau \times \tilde{U}_{h0}^\tau$, we define the residual norm

$$\|\hat{v}_h^\tau, \hat{\zeta}_h^\tau\|^2 := \|v_T^\tau\|_{\mathrm{R}}^2 + \|v_T^\tau\|_{J \times \varpi}^2 + t_h^\tau(v_T^\tau, v_T^\tau) + s_h^\tau(\hat{v}_h^\tau, \hat{v}_h^\tau) + \sigma_h(\hat{\zeta}_h^\tau, \hat{\zeta}_h^\tau), \quad (30)$$

with

$$\|v_T^\tau\|_{\mathrm{R}}^2 := \|h(\partial_t R^\tau(v_T^\tau) - \Delta_T v_T^\tau)\|_{J \times \Omega}^2 + \|h^{\frac{1}{2}}[\nabla v_T^\tau]_{\mathcal{F}_h^{\text{int}} \cdot \mathbf{n}_{\mathcal{F}_h^{\text{int}}}}\|_{J \times \mathcal{F}_h^{\text{int}}}^2, \quad (31)$$

and $\|h^{\frac{1}{2}}[\nabla v_T^\tau]_{\mathcal{F}_h^{\text{int}} \cdot \mathbf{n}_{\mathcal{F}_h^{\text{int}}}}\|_{J \times \mathcal{F}_h^{\text{int}}}^2 := \sum_{F \in \mathcal{F}_h^{\text{int}}} h \|[\nabla v_T^\tau]_F \cdot \mathbf{n}_F\|_{J \times F}^2$. We observe that $\|\cdot\|$ defines a norm on $\tilde{U}_h^\tau \times \tilde{U}_{h0}^\tau$. Our key stability result is the following inf-sup condition. Notice that this result implies the existence and uniqueness of the solution to the discrete problem.

**Lemma 3 (Inf-sup condition).** The following holds for all $(\hat{v}_h^\tau, \hat{\zeta}_h^\tau) \in \tilde{U}_h^\tau \times \tilde{U}_{h0}^\tau$,

$$\|\hat{v}_h^\tau, \hat{\zeta}_h^\tau\| \lesssim \sup_{(\hat{w}_h^\tau, \hat{\eta}_h^\tau) \in \tilde{U}_h^\tau \times \tilde{U}_{h0}^\tau \setminus \{(0,0)\}} \frac{A_h^\tau((\hat{v}_h^\tau, \hat{\zeta}_h^\tau), (\hat{w}_h^\tau, \hat{\eta}_h^\tau))}{\|\hat{w}_h^\tau, \hat{\eta}_h^\tau\|}. \quad (32)$$
Proof. Let us denote by $S$ the right-hand side of (32).

(i) We first use the test functions $\hat{w}_h^\tau := \hat{v}_h^\tau$ and $\hat{\eta}_h^\tau := -\hat{\zeta}_h^\tau$ to get

$$ \begin{aligned} \|v_T^\tau\|_{J \times \omega}^2 + t_h^\tau(v_T^\tau, v_T^\tau) + s_h^\tau(\hat{v}_h^\tau, \hat{v}_h^\tau) + \sigma_h(\hat{\zeta}_h^\tau, \hat{\zeta}_h^\tau) &= A_h^\tau((\hat{v}_h^\tau, \hat{\zeta}_h^\tau), (\hat{v}_h^\tau, -\hat{\zeta}_h^\tau)) \\ &\le S \|\hat{v}_h^\tau, \hat{\zeta}_h^\tau\|. \end{aligned} $$

(ii) Let $\hat{\eta}_h^\tau := (0, (\eta_F^\tau)_{F \in \mathcal{F}_h})$ with $\eta_F^\tau := h[\nabla v_T^\tau]_F \cdot \mathbf{n}_F$ for all $F \in \mathcal{F}_h^{\text{int}}$ and $\eta_F^\tau := 0$ for all $F \in \mathcal{F}_h^\partial$. Since $\|0, \hat{\eta}_h^\tau\| = \sigma_h(\hat{\eta}_h^\tau, \hat{\eta}_h^\tau)^{\frac{1}{2}}$, we have

$$ \begin{aligned} \|h^{\frac{1}{2}}[\nabla v_T^\tau]_{\mathcal{F}_h^{\text{int}} \cdot \mathbf{n}_{\mathcal{F}_h^{\text{int}}}}\|_{J \times \mathcal{F}_h^{\text{int}}}^2 &= a_h(\hat{v}_h^\tau, \hat{\eta}_h^\tau) \\ &= A_h^\tau((\hat{v}_h^\tau, \hat{\zeta}_h^\tau), (0, \hat{\eta}_h^\tau)) + \sigma_h(\hat{\zeta}_h^\tau, \hat{\eta}_h^\tau) \\ &\le S\sigma_h(\hat{\eta}_h^\tau, \hat{\eta}_h^\tau)^{\frac{1}{2}} + \sigma_h(\hat{\zeta}_h^\tau, \hat{\eta}_h^\tau) \\ &\le (S + \sigma_h(\hat{\zeta}_h^\tau, \hat{\zeta}_h^\tau)^{\frac{1}{2}})\sigma_h(\hat{\eta}_h^\tau, \hat{\eta}_h^\tau)^{\frac{1}{2}}. \end{aligned} $$

Moreover, we also have

$$ \sigma_h(\hat{\eta}_h^\tau, \hat{\eta}_h^\tau) = d_h(\hat{\eta}_h^\tau, \hat{\eta}_h^\tau) = \sum_{T \in \mathcal{T}_h} h^{-1} \|\eta_{\partial T}^\tau\|_{J \times \partial T}^2 \lesssim \|h^{\frac{1}{2}}[\nabla v_T^\tau]_{\mathcal{F}_h^{\text{int}} \cdot \mathbf{n}_{\mathcal{F}_h^{\text{int}}}}\|_{J \times \mathcal{F}_h^{\text{int}}}^2. $$

This implies that

$$ \|h^{\frac{1}{2}}[\nabla v_T^\tau]_{\mathcal{F}_h^{\text{int}} \cdot \mathbf{n}_{\mathcal{F}_h^{\text{int}}}}\|_{J \times \mathcal{F}_h^{\text{int}}}^2 \lesssim S^2 + \sigma_h(\hat{\zeta}_h^\tau, \hat{\zeta}_h^\tau). $$

(iii) We now consider $\hat{\eta}_h^\tau := (\eta_T^\tau, 0)$ with $\eta_T^\tau := h^2(\partial_t R^\tau(v_T^\tau) - \Delta_T v_T^\tau)$. We have

$$ \begin{aligned} & (\partial_t R^\tau(v_T^\tau), \eta_T^\tau)_{J \times \Omega} + a_h(\hat{v}_h^\tau, \hat{\eta}_h^\tau) \\ &= (\partial_t R^\tau(v_T^\tau), \eta_T^\tau)_{J \times \Omega} + \sum_{T \in \mathcal{T}_h} \{ (\nabla v_T^\tau, \nabla \eta_T^\tau)_{J \times T} - (\nabla v_T^\tau \cdot \mathbf{n}_T, \eta_T^\tau)_{J \times \partial T} \\ & \qquad - (v_T^\tau - v_{\partial T}^\tau, \nabla \eta_T^\tau \cdot \mathbf{n}_T)_{J \times \partial T} \} \\ &= (\partial_t R^\tau(v_T^\tau) - \Delta_T v_T^\tau, \eta_T^\tau)_{J \times \Omega} - \sum_{T \in \mathcal{T}_h} (v_T^\tau - v_{\partial T}^\tau, \nabla \eta_T^\tau \cdot \mathbf{n}_T)_{J \times \partial T}. \end{aligned} $$

Using the definition of $\hat{\eta}_h^\tau$ and recalling the rewriting (24) of $b^\tau$, we have

$$ \begin{aligned} & \|h(\partial_t R^\tau(v_T^\tau) - \Delta_T v_T^\tau)\|_{J \times \Omega}^2 \\ &= a_h(\hat{v}_h^\tau, \hat{\eta}_h^\tau) + b^\tau(v_T^\tau, \eta_T^\tau) + \sum_{T \in \mathcal{T}_h} (v_T^\tau - v_{\partial T}^\tau, \nabla \eta_T^\tau \cdot \mathbf{n}_T)_{J \times \partial T} \\ &= A_h^\tau((\hat{v}_h^\tau, \hat{\zeta}_h^\tau), (0, \hat{\eta}_h^\tau)) + \sigma_h(\hat{\zeta}_h^\tau, \hat{\eta}_h^\tau) + \sum_{T \in \mathcal{T}_h} (v_T^\tau - v_{\partial T}^\tau, \nabla \eta_T^\tau \cdot \mathbf{n}_T)_{J \times \partial T} \\ &\lesssim S \|0, \hat{\eta}_h^\tau\| + (\sigma_h(\hat{\zeta}_h^\tau, \hat{\zeta}_h^\tau) + d_h(\hat{v}_h^\tau, \hat{v}_h^\tau))^{\frac{1}{2}} \sigma_h(\hat{\eta}_h^\tau, \hat{\eta}_h^\tau)^{\frac{1}{2}}, \end{aligned} $$

since a discrete trace inverse inequality implies that

$$ \sum_{T \in \mathcal{T}_h} h \|\nabla \eta_T^\tau \cdot \mathbf{n}_T\|_{J \times \partial T}^2 \lesssim \|\nabla \eta_T^\tau\|_{J \times \Omega}^2 \le \sigma_h(\hat{\eta}_h^\tau, \hat{\eta}_h^\tau). $$
Furthermore, invoking again inverse inequalities, we have

$$
\begin{align*}
\sigma_h(\hat{\eta}_h^\tau, \hat{\eta}_h^\tau) &= \sum_{T \in \mathcal{T}_h} \left( \|\nabla \eta_T^\tau\|_{J \times T}^2 + h^{-1} \|\eta_T^\tau\|_{J \times \partial T}^2 \right) \\
&\lesssim \|h^{-1} \eta_T^\tau\|_{J \times \Omega}^2 = \|h(\partial_t R^\tau(v_T^\tau) - \Delta_T v_T^\tau)\|_{J \times \Omega}^2.
\end{align*}
$$

This implies that

$$
\|h(\partial_t R^\tau(v_T^\tau) - \Delta_T v_T^\tau)\|_{J \times \Omega}^2 \lesssim S^2 + d_h(\hat{v}_h^\tau, \hat{v}_h^\tau) + \sigma_h(\hat{\zeta}_h^\tau, \hat{\zeta}_h^\tau).
$$

(iv) Gathering the previous estimates leads to

$$
\|v_T^\tau\|_{J \times \bar{\omega}}^2 + t_h^\tau(v_T^\tau, v_T^\tau) + s_h^\tau(\hat{v}_h^\tau, \hat{v}_h^\tau) + \sigma_h(\hat{\zeta}_h^\tau, \hat{\zeta}_h^\tau) + \|h^{\frac{1}{2}}[\nabla v_T^\tau]\|_{\mathcal{F}_h^{\text{int}} \cdot \mathbf{n}_{\mathcal{F}_h^{\text{int}}}}^2 \|v_{J \times \mathcal{F}_h^{\text{int}}}\| \\
+ \|h(\partial_t R^\tau(v_T^\tau) - \Delta_T v_T^\tau)\|_{J \times \Omega}^2 \lesssim S \|\hat{v}_h^\tau, \hat{\zeta}_h^\tau\| + S^2.
$$

Recalling the definition (30) of the triple norm and using Young’s inequality gives the
expected inf-sup condition.
□

3.3 Interpolation operator and error decomposition

In this section, we define the space-time interpolation operator used in the error analysis.
Its definition is motivated by orthogonality properties. To facilitate the reading, the
approximation properties of this interpolation operator are discussed later in Section 3.7.

Let us first consider the approximation in time. For all $v \in H^1(J; H^1(\Omega))$ and all $n \in \{1:N\}$, we define $\check{I}_n^\ell(v) \in \mathbb{P}^\ell(I_n; H^1(\Omega))$ by

$$
\check{I}_n^\ell(v)(t_n^-) := v(t_n), \quad (33)
$$

$$
\int_{I_n} (\check{I}_n^\ell(v) - v, q^n)_\Omega := 0, \quad \forall q^n \in \mathbb{P}^{\ell-1}(I_n; H^1(\Omega)). \quad (34)
$$

Furthermore, the approximation in space is realized by using $L^2$-orthogonal projections.
Let $y \in H^1(\Omega)$. For all $T \in \mathcal{T}_h$, we define $\Pi_T^k(y)$ as the $L^2$-orthogonal projection of $y|_T$
onto $\mathbb{P}^k(T)$, i.e., $(\Pi_T^k(y), q_T)_T := (y, q_T)_T$ for all $q_T \in \mathbb{P}^k(T)$. For all $F \in \mathcal{F}_h$, we define
$\Pi_F^k(y)$ as the $L^2$-orthogonal projection of $y|_F$ onto $\mathbb{P}^k(F)$, i.e., $(\Pi_F^k(y), q_F)_F := (y, q_F)_F$
for all $q_F \in \mathbb{P}^k(F)$.

We can now define our space-time interpolation operator by setting, for all $v \in H^1(J; H^1(\Omega))$,

$$
\hat{I}_h^\tau(v) := (I_T^\tau(v), I_F^\tau(v)) \in \tilde{U}_h^\tau,
$$

where $I_T^\tau(v) \in U_T^\tau$ and $I_F^\tau(v) \in U_F^\tau$ are such that, for all $n \in \{1:N\}$, all $T \in \mathcal{T}_h$, and all $F \in \mathcal{F}_h$,

$$
I_T^n(v) := I_T^\tau(v)|_{I_n} := \Pi_T^k(\check{I}_n^\ell(v)), \quad I_F^n(v) := I_F^\tau(v)|_{I_n} := \Pi_F^k(\check{I}_n^\ell(v)). \quad (35)
$$

We notice that the approximation operators in space and in time commute. Moreover,
proceeding as in [19, Lemma 69.16], we derive the following useful orthogonality property:
For all $v \in H^1(J; H^1(\Omega))$ and all $w_T^\tau \in \mathbb{P}^\ell(J; U_T^k)$,

$$
(\partial_t v - \partial_t R^\tau(I_T^\tau(v)), w_T^\tau)_{J \times \Omega} = 0. \quad (36)
$$
Let us finally define the discrete and interpolation errors on the primal unknown. Recall that $u$ denotes the solution to the exact problem (1)-(2) and that $(\hat{u}_h^\tau, \hat{\xi}_h^\tau)$ denotes the solution to the discrete problem (20)-(21). The discrete and interpolation errors on the primal unknown are then defined as

$$ \hat{e}_h^\tau := \hat{u}_h^\tau - \hat{I}_h^\tau(u), \quad \hat{\theta}_h^\tau := (u, u|_{\mathcal{F}_h}) - \hat{I}_h^\tau(u), \qquad (37) $$

so that we have $\hat{\theta}_h^\tau := ((\theta_T^\tau)_{T \in \mathcal{T}_h}, (\theta_F^\tau)_{F \in \mathcal{F}_h})$ with $\theta_T^\tau := u|_T - I_T^\tau(u)$ for all $T \in \mathcal{T}_h$, and $\theta_F^\tau := u|_F - I_F^\tau(u)$ for all $F \in \mathcal{F}_h$.

### 3.4 Consistency and a priori residual bound

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

$$ \begin{aligned} \|\hat{\theta}_h^\tau\|_{\#}^2 := & \|\nabla_T \theta_T^\tau\|_{J \times \Omega}^2 + \|h \Delta_T \theta_T^\tau\|_{J \times \Omega}^2 + \|\theta_T^\tau\|_{J \times \Omega}^2 \\ & + \sum_{T \in \mathcal{T}_h} \{\|h^{\frac{1}{2}} \nabla_T \theta_T^\tau\|_{J \times \partial T}^2 + \|h^{-\frac{1}{2}} \theta_T^\tau\|_{J \times \partial T}^2 + \|h^{-\frac{1}{2}} \theta_T^\tau\|_{J \times \partial T}^2\} \\ & + \sum_{n \in \{1:N-1\}} \{\|\llbracket \theta_T^\tau \rrbracket^n\|_{\Omega}^2 + \|\llbracket \nabla_T \theta_T^\tau \rrbracket^n\|_{\Omega}^2\} + \tau \|\partial_t^\tau \theta_T^\tau\|_{J \times \Omega}^2. \end{aligned} \qquad (38) $$

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

$$ u \in H^1(J; H^1(\Omega)) \cap L^2(J; H^2(\Omega)). \qquad (39) $$

**Lemma 4 (Consistency and boundedness).** Let $(\hat{u}_h^\tau, \hat{\xi}_h^\tau)$ denote the solution to the discrete problem (20)-(21). Let $\hat{e}_h^\tau$ and $\hat{\theta}_h^\tau$ be the discrete and interpolation errors on the primal unknown defined in (37). Under the regularity assumption (39), we have, for all $(\hat{w}_h^\tau, \hat{\eta}_h^\tau) \in \tilde{U}_h^\tau \times \tilde{U}_{h0}^\tau$,

$$ |A_h^\tau((\hat{e}_h^\tau, \hat{\xi}_h^\tau), (\hat{w}_h^\tau, 0))| \lesssim (\|\hat{\theta}_h^\tau\|_{\#} + \|\delta\|_{J \times \varpi} + t_h^\tau(I_T^\tau(u), I_T^\tau(u))^{\frac{1}{2}}) \|\hat{w}_h^\tau, 0\|, \qquad (40) $$

$$ |A_h^\tau((\hat{e}_h^\tau, \hat{\xi}_h^\tau), (0, \hat{\eta}_h^\tau))| \lesssim \|\hat{\theta}_h^\tau\|_{\#} \|0, \hat{\eta}_h^\tau\|. \qquad (41) $$

*Proof.* (i) Proof of (40). Using $(\hat{e}_h^\tau, \hat{\xi}_h^\tau) = (\hat{u}_h^\tau, \hat{\xi}_h^\tau) - (\hat{I}_h^\tau(u), 0)$, the equation (20) in the discrete problem, the definition (22) of $A_h^\tau$, and $g_\delta = g + \delta$, we have, for all $\hat{w}_h^\tau \in \tilde{U}_h^\tau$,

$$ \begin{aligned} A_h^\tau((\hat{e}_h^\tau, \hat{\xi}_h^\tau), (\hat{w}_h^\tau, 0)) &= m_\varpi(g_\delta, w_\tau^\tau) - m_\varpi(I_T^\tau(u), w_\tau^\tau) - t_h^\tau(I_T^\tau(u), w_\tau^\tau) - s_h^\tau(\hat{I}_h^\tau(u), \hat{w}_h^\tau) \\ &= m_\varpi(\delta, w_\tau^\tau) + m_\varpi(\theta_T^\tau, w_\tau^\tau) - t_h^\tau(I_T^\tau(u), w_\tau^\tau) - s_h^\tau(\hat{I}_h^\tau(u), \hat{w}_h^\tau). \end{aligned} $$

The Cauchy–Schwarz inequality yields

$$ \begin{aligned} |m_\varpi(\delta, w_\tau^\tau)| + |m_\varpi(\theta_T^\tau, w_\tau^\tau)| &\lesssim (\|\delta\|_{J \times \varpi} + \|\theta_T^\tau\|_{J \times \varpi}) \|w_\tau^\tau\|_{J \times \varpi}, \\ |t_h^\tau(I_T^\tau(u), w_\tau^\tau)| &\lesssim t_h^\tau(I_T^\tau(u), I_T^\tau(u))^{\frac{1}{2}} t_h^\tau(w_\tau^\tau, w_\tau^\tau)^{\frac{1}{2}}, \\ |s_h^\tau(\hat{I}_h^\tau(u), \hat{w}_h^\tau)| &\lesssim s_h^\tau(\hat{I}_h^\tau(u), \hat{I}_h^\tau(u))^{\frac{1}{2}} s_h^\tau(\hat{w}_h^\tau, \hat{w}_h^\tau)^{\frac{1}{2}}. \end{aligned} $$
Moreover, we have

$$
\begin{aligned}
s_h^T(\hat{I}_h^T(u), \hat{I}_h^T(u)) = & \sum_{T \in \mathcal{T}_h} h^{-1} \|I_T^T(u) - I_{\partial T}^T(u)\|_{J \times \partial T}^2 \\
& + \sum_{n \in \{1:N-1\}} \{ \|\llbracket \theta_T^T \rrbracket^n\|_{\Omega}^2 + \|\llbracket \nabla_T \theta_T^T \rrbracket^n\|_{\Omega}^2 \} \lesssim \|\hat{\theta}_h^T\|_{\#}^2.
\end{aligned}
$$

The bound (40) follows by gathering the above estimates.

(ii) Proof of (41). Proceeding as above and using now the equation (21) in the discrete problem, we have, for all $\hat{\eta}_h^T \in \tilde{U}_{h0}^T$,

$$
\begin{aligned}
A_h^T((\hat{e}_h^T, \hat{\xi}_h^T), (0, \hat{\eta}_h^T)) &= (f, \eta_T^T)_{J \times \Omega} - b^T(I_T^T(u), \eta_T^T) - a_h(\hat{I}_h^T(u), \hat{\eta}_h^T) \\
&= -a_h(\hat{I}_h^T(u), \hat{\eta}_h^T) - (\Delta u, \eta_T^T)_{J \times \Omega},
\end{aligned}
$$

where we used that $f = \partial_t u - \Delta u$ in $J \times \Omega$ and $(\partial_t u, \eta_T^T)_{J \times \Omega} = b^T(I_T^T(u), \eta_T^T)$ owing to (36) and (24). Recalling the definition (10) of $a_h$, we infer that

$$
\begin{aligned}
A_h^T((\hat{e}_h^T, \hat{\xi}_h^T), (0, \hat{\eta}_h^T)) &= \sum_{T \in \mathcal{T}_h} \left\{ -(\nabla I_T^T(u), \nabla \eta_T^T)_{J \times T} + (\nabla I_T^T(u) \cdot n_T, \eta_T^T - \eta_{\partial T}^T)_{J \times \partial T} \right. \\
&\qquad \left. + (I_T^T(u) - I_{\partial T}^T(u), \nabla \eta_T^T \cdot n_T)_{J \times \partial T} - (\Delta u, \eta_T^T)_{J \times T} \right\} \\
&= \sum_{T \in \mathcal{T}_h} \left\{ (\nabla \theta_T^T, \nabla \eta_T^T)_{J \times T} - (\nabla \theta_T^T \cdot n_T, \eta_T^T - \eta_{\partial T}^T)_{J \times \partial T} \right. \\
&\qquad \left. + (\theta_{\partial T}^T - \theta_T^T, \nabla \eta_T^T \cdot n_T)_{J \times \partial T} \right\},
\end{aligned}
$$

where we used the assumed regularity of $u$ and the fact that $\eta_F^T = 0$ for all $F \in \mathcal{F}_h^\partial$ to infer that $\sum_{T \in \mathcal{T}_h} (\nabla u \cdot n_T, \eta_{\partial T}^T)_{J \times \partial T} = 0$. Using Cauchy-Schwarz and inverse inequalities, we readily obtain (41). $\square$

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

$$
\|\hat{e}_h^T, \hat{\xi}_h^T\| + t_h^T(u_T^T, u_T^T)^{\frac{1}{2}} + s_h^T(\hat{u}_h^T, \hat{u}_h^T)^{\frac{1}{2}} \lesssim \|\hat{\theta}_h^T\|_{\#} + \|\delta\|_{J \times \varpi} + t_h^T(I_T^T(u), I_T^T(u))^{\frac{1}{2}}.
$$

*Proof.* Using inf-sup stability (Lemma 3) yields

$$
\|\hat{e}_h^T, \hat{\xi}_h^T\| \lesssim \sup_{(\hat{w}_h^T, \hat{\eta}_h^T) \in \tilde{U}_h^T \times \tilde{U}_{h0}^T \setminus \{(0,0)\}} \frac{A_h^T((\hat{e}_h^T, \hat{\xi}_h^T), (\hat{w}_h^T, \hat{\eta}_h^T))}{\|\hat{w}_h^T, \hat{\eta}_h^T\|}.
$$

Moreover, owing to consistency (Lemma 4), we have, for all $(\hat{w}_h^T, \hat{\eta}_h^T) \in \tilde{U}_h^T \times \tilde{U}_{h0}^T$,

$$
|A_h^T((\hat{e}_h^T, \hat{\xi}_h^T), (\hat{w}_h^T, \hat{\eta}_h^T))| \lesssim (\|\hat{\theta}_h^T\|_{\#} + \|\delta\|_{J \times \varpi} + t_h^T(I_T^T(u), I_T^T(u))^{\frac{1}{2}}) \|\hat{w}_h^T, \hat{\eta}_h^T\|.
$$

Combining the two bounds gives

$$
\|\hat{e}_h^T, \hat{\xi}_h^T\| \lesssim \|\hat{\theta}_h^T\|_{\#} + \|\delta\|_{J \times \varpi} + t_h^T(I_T^T(u), I_T^T(u))^{\frac{1}{2}}.
$$

Finally, observing that

$$
\begin{aligned}
t_h^T(u_T^T, u_T^T)^{\frac{1}{2}} &\le t_h^T(I_T^T(u), I_T^T(u))^{\frac{1}{2}} + t_h^T(e_T^T, e_T^T)^{\frac{1}{2}} \le t_h^T(I_T^T(u), I_T^T(u))^{\frac{1}{2}} + \|\hat{e}_h^T, \hat{\xi}_h^T\|, \\
s_h^T(\hat{u}_h^T, \hat{u}_h^T)^{\frac{1}{2}} &\le s_h^T(\hat{I}_h^T(u), \hat{I}_h^T(u))^{\frac{1}{2}} + s_h^T(\hat{e}_h^T, \hat{e}_h^T)^{\frac{1}{2}} \le \|\hat{\theta}_h^T\|_{\#} + \|\hat{e}_h^T, \hat{\xi}_h^T\|,
\end{aligned}
$$

concludes the proof. $\square$

<figure><img src="image_2.png" /></figure>
Remark 3.1. (Lemma 5) The a priori residual bound involves three terms, one depending on the interpolation error of the exact solution measured in the $||\cdot||_{\#}$-norm, one resulting from the presence of noise in the measurements, and one resulting from the regularization.

## 3.5 Bound on dual error norm

In this section, we prove another important result bounding some dual norm of the error. Recall the operator $L := \partial_t - \Delta$ so that, for all $v \in H^1(J; L^2(\Omega)) \cap L^2(J; H^1(\Omega))$ and all $\eta \in L^2(J; H_0^1(\Omega))$, we have $\langle L(v), \eta \rangle_{L^2(H^{-1}), L^2(H_0^1)} = (\partial_t v, \eta)_{J \times \Omega} + (\nabla v, \nabla \eta)_{J \times \Omega}$, where $\langle \cdot, \cdot \rangle_{L^2(H^{-1}), L^2(H_0^1)}$ stands for the duality product between $L^2(J; H^{-1}(\Omega))$ and $L^2(J; H_0^1(\Omega))$. The operator $L$ can be extended to space-discrete functions $v_T^\tau \in H^1(J; U_T^k)$ by setting

$$ \langle L_T(v_T^\tau), \eta \rangle_{L^2(H^{-1}), L^2(H_0^1)} := (\partial_t v_T^\tau, \eta)_{J \times \Omega} + (\nabla_T v_T^\tau, \nabla \eta)_{J \times \Omega}. \quad (42) $$

The corresponding dual norm is

$$ \|L_T(v)\|_{L^2(J; H^{-1}(\Omega))} = \sup_{\substack{\eta \in L^2(J; H_0^1(\Omega)) \\ \|\nabla \eta\|_{J \times \Omega} = 1}} \{ (\partial_t v, \eta)_{J \times \Omega} + (\nabla_T v, \nabla \eta)_{J \times \Omega} \}. $$

It is also useful to introduce the data oscillation term

$$ \Theta_f := f - \Pi^\tau(f), \quad (43) $$

where $\Pi^\tau(f)|_{I_n} := \Pi_n^\ell(f|_{I_n})$ for all $n \in \{1:N\}$, and $\Pi_n^\ell$ denotes the $L^2$-orthogonal projection from $L^2(I_n; L^2(\Omega))$ onto $\mathbb{P}^\ell(I_n; L^2(\Omega))$.

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

$$ \begin{aligned} & \|L_T(R^\tau(u_T^\tau) - u)\|_{L^2(J; H^{-1}(\Omega))} \\ & \lesssim (1 + \tau^{-\frac{1}{2}}h)\|\hat{\theta}_h^\tau\|_{\#} + \|\delta\|_{J \times \varpi} + t_h^\tau(I_T^\tau(u), I_T^\tau(u))^{\frac{1}{2}} + \|\Theta_f\|_{J \times \Omega}. \end{aligned} $$

*Proof.* Let $\eta \in L^2(J; H_0^1(\Omega))$ with $\|\nabla \eta\|_{J \times \Omega} = 1$. Recalling the definition of $\Theta_f$, we have

$$ \begin{aligned} & \langle L_T(R^\tau(u_T^\tau) - u), \eta \rangle_{L^2(H^{-1}), L^2(H_0^1)} \\ & = (\partial_t R^\tau(u_T^\tau), \eta)_{J \times \Omega} + (\nabla_T R^\tau(u_T^\tau), \nabla \eta)_{J \times \Omega} - (\Pi^\tau(f), \eta)_{J \times \Omega} - (\Theta_f, \eta)_{J \times \Omega} \\ & = (\partial_t R^\tau(u_T^\tau), \Pi^\tau(\eta))_{J \times \Omega} + (\nabla_T R^\tau(u_T^\tau), \nabla \eta)_{J \times \Omega} - (f, \Pi^\tau(\eta))_{J \times \Omega} - (\Theta_f, \eta)_{J \times \Omega}, \end{aligned} $$

where we used the $L^2$-orthogonality properties of $\Pi^\tau$. Let $\hat{\eta}_h^\tau \in \tilde{U}_{h0}^\tau$ be such that

$$ \hat{\eta}_h^\tau|_{I_n} := ((\Pi_n^\ell \Pi_T^k(\eta|_T))_{T \in \mathcal{T}_h}, (\Pi_n^\ell \Pi_F^k(\eta|_F))_{F \in \mathcal{F}_h}) \in \mathbb{P}^\ell(I_n; \hat{U}_{h0}^k), $$

for all $n \in \{1:N\}$. Invoking the second equation (21) in the discrete problem and recalling the rewriting (24) of $b^\tau$ gives

$$ (\partial_t R^\tau(u_T^\tau), \eta_T^\tau)_{J \times \Omega} + a_h(\hat{u}_h^\tau, \hat{\eta}_h^\tau) - \sigma_h(\hat{\xi}_h^\tau, \hat{\eta}_h^\tau) = (f, \eta_T^\tau)_{J \times \Omega}. $$
Subtracting this equation from the above expression, recalling that $f = \partial_t u - \Delta u$, and re-arranging the terms leads to

$$ \langle L_{\mathcal{T}}(R^{\tau}(u_{\mathcal{T}}^{\tau}) - u), \eta \rangle_{L^2(H^{-1}), L^2(H_0^1)} = A_1 + A_2 + A_3, $$

with

$$ A_1 := (\partial_t(R^\tau(u_T^\tau) - u) - \Delta_T(u_T^\tau - u), \Pi^\tau(\eta) - \eta_T^\tau)_{J \times \Omega}, $$

$$ A_2 := (\nabla_{\mathcal{T}} R^{\tau}(u_{\mathcal{T}}^{\tau}), \nabla \eta)_{J \times \Omega} - a_h(\hat{u}_h^{\tau}, \hat{\eta}_h^{\tau}) + (\Delta_{\mathcal{T}} u_{\mathcal{T}}^{\tau}, \Pi^{\tau}(\eta) - \eta_{\mathcal{T}}^{\tau})_{J \times \Omega}, $$

$$ A_3 := \sigma_h(\hat{\xi}_h^{\tau}, \hat{\eta}_h^{\tau}) - (\Theta_f, \eta)_{J \times \Omega}. $$

A straightforward calculation using the definition (10) of the bilinear form $a_h$ shows that

$$ A_2 = (\nabla_{\mathcal{T}}(R^{\tau}(u_{\mathcal{T}}^{\tau}) - u_{\mathcal{T}}^{\tau}), \nabla \eta)_{J \times \Omega} + \sum_{T \in \mathcal{T}_h} (u_T^{\tau} - u_{\partial T}^{\tau}, \nabla \eta_T^{\tau} \cdot n_T)_{J \times \partial T}, $$

since we have

$$ \begin{aligned} & (\nabla_{\mathcal{T}} u_{\mathcal{T}}^{\tau}, \nabla \eta)_{J \times \Omega} - a_h(\hat{u}_h^{\tau}, \hat{\eta}_h^{\tau}) \\ &= (\nabla_{\mathcal{T}} u_{\mathcal{T}}^{\tau}, \nabla_{\mathcal{T}} (\Pi^{\tau}(\eta) - \eta_{\mathcal{T}}^{\tau}))_{J \times \Omega} \\ &\quad + \sum_{T \in \mathcal{T}_h} \{ (\nabla u_T^{\tau} \cdot n_T, \eta_T^{\tau} - \eta_{\partial T}^{\tau})_{J \times \partial T} + (u_T^{\tau} - u_{\partial T}^{\tau}, \nabla \eta_T^{\tau} \cdot n_T)_{J \times \partial T} \} \\ &= (-\Delta_{\mathcal{T}} u_{\mathcal{T}}^{\tau}, \Pi^{\tau}(\eta) - \eta_{\mathcal{T}}^{\tau})_{J \times \Omega} \\ &\quad + \sum_{T \in \mathcal{T}_h} \{ (\nabla u_T^{\tau} \cdot n_T, \Pi^{\tau}(\eta) - \eta_{\partial T}^{\tau})_{J \times \partial T} + (u_T^{\tau} - u_{\partial T}^{\tau}, \nabla \eta_T^{\tau} \cdot n_T)_{J \times \partial T} \} \\ &= (-\Delta_{\mathcal{T}} u_{\mathcal{T}}^{\tau}, \Pi^{\tau}(\eta) - \eta_{\mathcal{T}}^{\tau})_{J \times \Omega} + \sum_{T \in \mathcal{T}_h} (u_T^{\tau} - u_{\partial T}^{\tau}, \nabla \eta_T^{\tau} \cdot n_T)_{J \times \partial T}, \end{aligned} $$

where we used that $(\nabla_{\mathcal{T}} u_{\mathcal{T}}^{\tau}, \nabla(\eta - \Pi^{\tau}(\eta)))_{J \times \Omega} = 0$ in the first equality and the fact that $(\nabla u_{\mathcal{T}}^{\tau} \cdot n_T)|_{I_n \times \partial T} \in \mathbb{P}^{\ell}(I_n; \mathbb{P}^k(\mathcal{F}_{\partial T}))$ and the definition of $\eta_{\partial T}^{\tau}$ in the third equality. It remains to bound $A_1, A_2, A_3$.

(i) Bound on $A_1$. This is the most delicate term. We observe that

$$ \begin{aligned} \partial_t (R^\tau (u_T^\tau) - u) - \Delta_T (u_T^\tau - u) &= \partial_t R^\tau (u_T^\tau - I_T^\tau(u)) - \Delta_T (u_T^\tau - I_T^\tau(u)) \\ &\quad + \partial_t (R^\tau (I_T^\tau(u)) - u) - \Delta_T (I_T^\tau(u) - u). \end{aligned} $$

Recalling the definition (31) of the $\|\cdot\|_{\mathbb{R}}$-norm, we have

$$ |(\partial_t R^\tau (u_T^\tau - I_T^\tau(u)) - \Delta_T (u_T^\tau - I_T^\tau(u)), \Pi^\tau(\eta) - \eta_T^\tau)_{J \times \Omega}| \lesssim \|e_T^\tau\|_{\mathbb{R}} \|\nabla \eta\|_{J \times \Omega}, $$

where we used that

$$ \|\Pi^{\tau}(\eta) - \eta_{\mathcal{T}}^{\tau}\|_{J \times \Omega} \lesssim h \|\nabla \eta\|_{J \times \Omega}. $$

Owing to Lemma 5, we have $\|e_T^\tau\|_{\mathbb{R}} \lesssim \|\hat{\theta}_h^\tau\|_{\#} + \|\delta\|_{J \times \varpi} + t_h^\tau(I_T^\tau(u), I_T^\tau(u))^{\frac{1}{2}}$. Moreover, the triangle inequality, the estimate (27), and the above bound on $\|\Pi^\tau(\eta) - \eta_T^\tau\|_{J \times \Omega}$ give

$$ \begin{aligned} & |(\partial_t (R^\tau(I_T^\tau(u)) - u), \Pi^\tau(\eta) - \eta_T^\tau)_{J \times \Omega}| \\ &\lesssim (\|\partial_t^\tau(R^\tau(I_T^\tau(u)) - I_T^\tau(u))\|_{J \times \Omega} + \|\partial_t^\tau(u - I_T^\tau(u))\|_{J \times \Omega})h \|\nabla \eta\|_{J \times \Omega} \\ &\lesssim (d^\tau(I_T^\tau(u), I_T^\tau(u))^{\frac{1}{2}} + \tau^{\frac{1}{2}} \|\partial_t^\tau \theta_T^\tau\|_{J \times \Omega}) \tau^{-\frac{1}{2}} h \|\nabla \eta\|_{J \times \Omega}. \end{aligned} $$
Finally, we have

$$|(\Delta_{\tau}(I_{\tau}^{T}(u) - u), \Pi^{\tau}(\eta) - \eta_{\tau})_{J \times \Omega}| \lesssim h \|\Delta_{\tau}\theta_{\tau}^{T}\|_{J \times \Omega} \|\nabla\eta\|_{J \times \Omega}.$$

Gathering the above estimates gives

$$|A_1| \lesssim ((1 + \tau^{-\frac{1}{2}}h)\|\hat{\theta}_h^T\|_{\#} + \|\delta\|_{J \times \varpi} + t_h^T(I_{\tau}^T(u), I_{\tau}^T(u))^{\frac{1}{2}})\|\nabla\eta\|_{J \times \Omega}.$$

(ii) Bound on $A_2$. Owing to the Cauchy-Schwarz inequality and (28), we have

$$|(\nabla_{\tau}(R^{\tau}(u_{\tau}^{T}) - u_{\tau}^{T}), \nabla\eta)_{J \times \Omega}| \lesssim \tau^{\frac{1}{2}} d^{\tau}(u_{\tau}^{T}, u_{\tau}^{T})^{\frac{1}{2}} \|\nabla\eta\|_{J \times \Omega}.$$

Moreover, using Cauchy-Schwarz and inverse trace inequalities, we have

$$\left| \sum_{T \in \mathcal{T}_h} (u_T^T - u_{\partial T}^T, \nabla \eta_T^T \cdot n_T)_{J \times \partial T} \right| \lesssim d_h(\hat{u}_h^T, \hat{u}_h^T)^{\frac{1}{2}} \|\nabla \eta\|_{J \times \Omega}.$$

Thus, we have

$$|A_2| \lesssim (\tau^{\frac{1}{2}} d^{\tau}(u_{\tau}^{T}, u_{\tau}^{T})^{\frac{1}{2}} + d_h(\hat{u}_h^{T}, \hat{u}_h^{T})^{\frac{1}{2}}) \|\nabla \eta\|_{J \times \Omega}.$$

Using $\tau \lesssim 1$ and Lemma 5 yields

$$|A_2| \lesssim s_h^T(\hat{u}_h^T, \hat{u}_h^T)^{\frac{1}{2}} \|\nabla\eta\|_{J \times \Omega} \lesssim (\|\hat{\theta}_h^T\|_{\#} + \|\delta\|_{J \times \varpi} + t_h^T(I_{\tau}^T(u), I_{\tau}^T(u))^{\frac{1}{2}}) \|\nabla\eta\|_{J \times \Omega}.$$

(iii) Bound on $A_3$. Since $\sigma_h(\hat{\eta}_h^T, \hat{\eta}_h^T) = \|\nabla_{\tau}\eta_T^T\|_{J \times \Omega}^2 + d_h(\hat{\eta}_h^T, \hat{\eta}_h^T)$, we infer that $\sigma_h(\hat{\eta}_h^T, \hat{\eta}_h^T)^{\frac{1}{2}} \lesssim \|\nabla\eta\|_{J \times \Omega}$. This implies that $|\sigma_h(\hat{\xi}_h^T, \hat{\eta}_h^T)| \lesssim \sigma_h(\hat{\xi}_h^T, \hat{\xi}_h^T)^{\frac{1}{2}} \sigma_h(\hat{\eta}_h^T, \hat{\eta}_h^T)^{\frac{1}{2}} \lesssim \sigma_h(\hat{\xi}_h^T, \hat{\xi}_h^T)^{\frac{1}{2}} \|\nabla\eta\|_{J \times \Omega}$ and $\sigma_h(\hat{\xi}_h^T, \hat{\xi}_h^T)^{\frac{1}{2}}$ is bounded in Lemma 5. Moreover, the bound $|(\Theta_f, \eta)_{J \times \Omega}| \lesssim \|\Theta_f\|_{J \times \Omega} \|\nabla\eta\|_{J \times \Omega}$ results from the Cauchy-Schwarz inequality combined and the global Poincaré inequality in $\Omega$. Hence, we have

$$|A_3| \lesssim (\|\hat{\theta}_h^T\|_{\#} + \|\delta\|_{J \times \varpi} + t_h^T(I_{\tau}^T(u), I_{\tau}^T(u))^{\frac{1}{2}} + \|\Theta_f\|_{J \times \Omega}) \|\nabla\eta\|_{J \times \Omega}.$$

(iv) Combining the above estimates and recalling that $\|\nabla\eta\|_{J \times \Omega} = 1$ gives the expected bound. $\square$

## 3.6 Main result: error estimate in the target subdomain

We are now ready to derive our main error estimate. The idea of the proof is to combine the results of Sections 3.4 and 3.5 with the conditional stability estimate from Lemma 1. We use the target (semi)norm

$$\|v\|_{\text{trg}} := \|\nabla_{\tau} v\|_{L^2(T_1, T_2; L^2(B))},$$

where $T_1, T_2$ and $B$ are defined in Lemma 1. Recall the quantities $C_{\text{stb}} > 0$ and $\alpha \in (0, 1]$ introduced in Lemma 1. Recall that the $\|\cdot\|_{\#}$-norm is defined in (38), the data oscillation term $\Theta_f$ in (43), that $\hat{\theta}_h^T$ is the interpolation error defined in (37), and that the coefficient $c(\tau, h)$ is used to weight the regularization bilinear form $t_h^T$ defined in (16).
**Theorem 7 (Error estimate in target subdomain).** Under the regularity assumption (39), we have

$$ \|u - u_{\mathcal{T}}^{\tau}\|_{\text{trg}} \lesssim C_{\text{stb}}(1 + \tau^{-1}h)^{\alpha}(1 + \tau^{-1}h + \gamma^{-\frac{1}{2}}c(\tau, h)^{-1})^{1-\alpha} \\ \times (\|\hat{\theta}_{h}^{\tau}\|_{\#} + \|\delta\|_{J \times \bar{\omega}} + t_{h}^{\tau}(I_{\mathcal{T}}^{\tau}(u), I_{\mathcal{T}}^{\tau}(u))^{\frac{1}{2}} + \|\Theta_{f}\|_{J \times \Omega}). \quad (44) $$

*Proof.* (i) Using an averaging operator in space (see [7, Lemmas 3.2 and 5.3] or [18, Chap. 22] and the references therein), we can build a piecewise polynomial function $\tilde{u}_{\mathcal{T}}^{\tau} \in \{v_{\mathcal{T}}^{\tau} \in H^{1}(J; H^{1}(\Omega)) | v_{\mathcal{T}}^{\tau}|_{I_{n} \times T} \in \mathbb{P}^{l+1}(I_{n}; \mathbb{P}^{k}(T)), n \in \{1:N\}, T \in \mathcal{T}_{h}\}$ such that

$$ \|h^{-1}(\tilde{u}_{\mathcal{T}}^{\tau} - R^{\tau}(u_{\mathcal{T}}^{\tau}))\|_{J \times \Omega} + \|\nabla_{\mathcal{T}}(\tilde{u}_{\mathcal{T}}^{\tau} - R^{\tau}(u_{\mathcal{T}}^{\tau}))\|_{J \times \Omega} \lesssim \|h^{-\frac{1}{2}}[R^{\tau}(u_{\mathcal{T}}^{\tau})]\|_{\mathcal{F}_{h}^{\text{int}}} \|_{J \times \mathcal{F}_{h}^{\text{int}}} \\ \lesssim d_{h}(\hat{u}_{h}^{\tau}, \hat{u}_{h}^{\tau})^{\frac{1}{2}}, \quad (45) $$

where the first bound results from the approximation properties in space of the averaging operator and the second bound from (29) (see Lemma 2). Invoking the triangle inequality yields

$$ \|u - u_{\mathcal{T}}^{\tau}\|_{\text{trg}} \le \|u_{\mathcal{T}}^{\tau} - \tilde{u}_{\mathcal{T}}^{\tau}\|_{\text{trg}} + \|u - \tilde{u}_{\mathcal{T}}^{\tau}\|_{\text{trg}}, $$

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

(ii) Bound on $\|u_{\mathcal{T}}^{\tau} - \tilde{u}_{\mathcal{T}}^{\tau}\|_{\text{trg}}$. We observe that

$$ \|u_{\mathcal{T}}^{\tau} - \tilde{u}_{\mathcal{T}}^{\tau}\|_{\text{trg}} \le \|\nabla_{\mathcal{T}}(u_{\mathcal{T}}^{\tau} - R^{\tau}(u_{\mathcal{T}}^{\tau}))\|_{J \times \Omega} + \|\nabla_{\mathcal{T}}(\tilde{u}_{\mathcal{T}}^{\tau} - R^{\tau}(u_{\mathcal{T}}^{\tau}))\|_{J \times \Omega}. $$

The first term on the right-hand side is estimated by invoking (28) (see Lemma 2) and the second term by invoking (45). Owing to Lemma 5 and since $\tau \lesssim 1$, we infer that

$$ \|u_{\mathcal{T}}^{\tau} - \tilde{u}_{\mathcal{T}}^{\tau}\|_{\text{trg}} \le d^{\tau}(u_{\mathcal{T}}^{\tau}, u_{\mathcal{T}}^{\tau})^{\frac{1}{2}} + d_{h}(\hat{u}_{h}^{\tau}, \hat{u}_{h}^{\tau})^{\frac{1}{2}} \lesssim s_{h}^{\tau}(\hat{u}_{h}^{\tau}, \hat{u}_{h}^{\tau})^{\frac{1}{2}} \\ \lesssim \|\hat{\theta}_{h}^{\tau}\|_{\#} + \|\delta\|_{J \times \bar{\omega}} + t_{h}^{\tau}(I_{\mathcal{T}}^{\tau}(u), I_{\mathcal{T}}^{\tau}(u))^{\frac{1}{2}}. $$

(iii) Bound on $\|u - \tilde{u}_{\mathcal{T}}^{\tau}\|_{\text{trg}}$. Since $\tilde{u}_{\mathcal{T}}^{\tau} - u \in H^{1}(J; H^{1}(\Omega)) \subset H^{1}(J; H^{-1}(\Omega)) \cap L^{2}(J; H^{1}(\Omega))$, we can invoke the conditional stability estimate from Lemma 1. This gives

$$ \|\tilde{u}_{\mathcal{T}}^{\tau} - u\|_{\text{trg}} \lesssim C_{\text{stb}} \left( \|\tilde{u}_{\mathcal{T}}^{\tau} - u\|_{J \times \bar{\omega}} + \|L(\tilde{u}_{\mathcal{T}}^{\tau} - u)\|_{L^{2}(J; H^{-1}(\Omega))} \right)^{\alpha} \\ \times \left( \|\tilde{u}_{\mathcal{T}}^{\tau} - u\|_{J \times \Omega} + \|L(\tilde{u}_{\mathcal{T}}^{\tau} - u)\|_{L^{2}(J; H^{-1}(\Omega))} \right)^{1-\alpha}. $$

It remains to bound $\|\tilde{u}_{\mathcal{T}}^{\tau} - u\|_{J \times \bar{\omega}}$, $\|\tilde{u}_{\mathcal{T}}^{\tau} - u\|_{J \times \Omega}$, and $\|L(\tilde{u}_{\mathcal{T}}^{\tau} - u)\|_{L^{2}(J; H^{-1}(\Omega))}$.

(iii.a) Bound on $\|\tilde{u}_{\mathcal{T}}^{\tau} - u\|_{J \times \bar{\omega}}$. We have

$$ \|\tilde{u}_{\mathcal{T}}^{\tau} - u\|_{J \times \bar{\omega}} \le \|u - I_{\mathcal{T}}^{\tau}(u)\|_{J \times \bar{\omega}} + \|u_{\mathcal{T}}^{\tau} - I_{\mathcal{T}}^{\tau}(u)\|_{J \times \bar{\omega}} + \|u_{\mathcal{T}}^{\tau} - \tilde{u}_{\mathcal{T}}^{\tau}\|_{J \times \bar{\omega}}, $$

where $I_{\mathcal{T}}^{\tau}(u)$ is defined in Section 3.3. We have

$$ \|u - I_{\mathcal{T}}^{\tau}(u)\|_{J \times \bar{\omega}} = \|\theta_{\mathcal{T}}^{\tau}\|_{J \times \bar{\omega}} \le \|\theta_{\mathcal{T}}^{\tau}\|_{J \times \Omega}, \\ \|u_{\mathcal{T}}^{\tau} - I_{\mathcal{T}}^{\tau}(u)\|_{J \times \bar{\omega}} = \|e_{\mathcal{T}}^{\tau}\|_{J \times \bar{\omega}} \le \|\hat{e}_{h}^{\tau}, \hat{\xi}_{h}^{\tau}\|. $$
Furthermore, the bound (28) (see Lemma 2) and the bound (45) imply that

$$ \begin{aligned} \|u_{\tau}^{\tau} - \tilde{u}_{\tau}^{\tau}\|_{J \times \bar{\omega}} &\le \|u_{\tau}^{\tau} - R^{\tau}(u_{\tau}^{\tau})\|_{J \times \bar{\omega}} + \|R^{\tau}(u_{\tau}^{\tau}) - \tilde{u}_{\tau}^{\tau}\|_{J \times \bar{\omega}} \\ &\lesssim \tau^{\frac{1}{2}} d^{\tau}(u_{\tau}^{\tau}, u_{\tau}^{\tau})^{\frac{1}{2}} + h d_h(\hat{u}_h^{\tau}, \hat{u}_h^{\tau})^{\frac{1}{2}} \lesssim s_h^{\tau}(\hat{u}_h^{\tau}, \hat{u}_h^{\tau})^{\frac{1}{2}}, \end{aligned} $$

where we used the definition (12) of $s_h^\tau$ and $\tau \lesssim 1$, $h \lesssim 1$. Gathering the above estimates and using Lemma 5 gives

$$ \|\tilde{u}_{\tau}^{\tau} - u\|_{J \times \bar{\omega}} \lesssim \|\hat{\theta}_{h}^{\tau}\|_{\#} + \|\delta\|_{J \times \bar{\omega}} + t_{h}^{\tau}(I_{\tau}^{\tau}(u), I_{\tau}^{\tau}(u))^{\frac{1}{2}}. $$

(iii.b) Bound on $\|\tilde{u}_{\tau}^{\tau} - u\|_{J \times \Omega}$. Using the same decomposition, we have

$$ \|\tilde{u}_{\tau}^{\tau} - u\|_{J \times \Omega} \le \|u - I_{\tau}^{\tau}(u)\|_{J \times \Omega} + \|u_{\tau}^{\tau} - I_{\tau}^{\tau}(u)\|_{J \times \Omega} + \|u_{\tau}^{\tau} - \tilde{u}_{\tau}^{\tau}\|_{J \times \Omega}. $$

The terms $\|u - I_{\tau}^{\tau}(u)\|_{J \times \Omega} + \|u_{\tau}^{\tau} - \tilde{u}_{\tau}^{\tau}\|_{J \times \Omega}$ are estimated by using the same arguments as in Step (iii.a). The bound on $\|u_{\tau}^{\tau} - I_{\tau}^{\tau}(u)\|_{J \times \Omega}$ uses, however, a different argument since we can only invoke here the regularization. Recalling the definition (16) of $t_h^\tau$, we have

$$ \|u_{\tau}^{\tau} - I_{\tau}^{\tau}(u)\|_{J \times \Omega} = \|e_{\tau}^{\tau}\|_{J \times \Omega} = \gamma^{-\frac{1}{2}} c(\tau, h)^{-1} t_{h}^{\tau}(e_{\tau}^{\tau}, e_{\tau}^{\tau})^{\frac{1}{2}}. $$

Invoking Lemma 5 gives

$$ \|\tilde{u}_{\tau}^{\tau} - u\|_{J \times \Omega} \lesssim (1 + \gamma^{-\frac{1}{2}} c(\tau, h)^{-1}) (\|\hat{\theta}_{h}^{\tau}\|_{\#} + \|\delta\|_{J \times \bar{\omega}} + t_{h}^{\tau}(I_{\tau}^{\tau}(u), I_{\tau}^{\tau}(u))^{\frac{1}{2}}). $$

(iii.c) Bound on $\|L(\tilde{u}_{\tau}^{\tau} - u)\|_{L^2(J; H^{-1}(\Omega))}$. Recalling that $L_{\tau}$ denotes the extension of $L$ to $H^1(J; U_{\tau}^k)$, the triangle inequality gives

$$ \|L(\tilde{u}_{\tau}^{\tau} - u)\|_{L^2(J; H^{-1}(\Omega))} \le \|L_{\tau}(\tilde{u}_{\tau}^{\tau} - R^{\tau}(u_{\tau}^{\tau}))\|_{L^2(J; H^{-1}(\Omega))} \\ + \|L_{\tau}(R^{\tau}(u_{\tau}^{\tau}) - u)\|_{L^2(J; H^{-1}(\Omega))}. $$

For the second term on the right-hand side, we invoke Lemma 6 to infer

$$ \begin{aligned} \|L_{\tau}(R^{\tau}(u_{\tau}^{\tau}) - u)\|_{L^2(J; H^{-1}(\Omega))} &\lesssim (1 + \tau^{-\frac{1}{2}} h) \|\hat{\theta}_{h}^{\tau}\|_{\#} \\ &\quad + \|\delta\|_{J \times \bar{\omega}} + t_{h}^{\tau}(I_{\tau}^{\tau}(u), I_{\tau}^{\tau}(u))^{\frac{1}{2}} + \|\Theta_{f}\|_{J \times \Omega}. \end{aligned} $$

To estimate the first term on the right-hand side, we bound the dual norm by considering an arbitrary test function $\eta \in L^2(J; H_0^1(\Omega))$ with $\|\nabla \eta\|_{J \times \Omega} = 1$. On the one hand, we have

$$ \begin{aligned} |(\partial_t(\tilde{u}_{\tau}^{\tau} - R^{\tau}(u_{\tau}^{\tau})), \eta)_{J \times \Omega}| &\lesssim \|\partial_t(\tilde{u}_{\tau}^{\tau} - R^{\tau}(u_{\tau}^{\tau}))\|_{J \times \Omega} \|\eta\|_{J \times \Omega} \\ &\lesssim \tau^{-1} h d_h(\hat{u}_h^{\tau}, \hat{u}_h^{\tau})^{\frac{1}{2}} \|\nabla \eta\|_{J \times \Omega}, \end{aligned} $$

where we used an inverse inequality in time, the estimate (45), and a global Poincaré inequality in $\Omega$ for $\eta$. On the other hand, using again (45) gives

$$ |(\nabla \tau(\tilde{u}_{\tau}^{\tau} - R^{\tau}(u_{\tau}^{\tau})), \nabla \eta)_{J \times \Omega}| \lesssim d_h(\hat{u}_h^{\tau}, \hat{u}_h^{\tau})^{\frac{1}{2}} \|\nabla \eta\|_{J \times \Omega}. $$
Combining the last two bounds and invoking Lemma 5, we infer that

$$ \begin{aligned} \|L_T(\tilde{u}_T^\tau - R^\tau(u_T^\tau))\|_{L^2(J; H^{-1}(\Omega))} &\lesssim (1 + \tau^{-1}h) d_h(\hat{u}_h^\tau, \hat{u}_h^\tau)^{\frac{1}{2}} \\ &\lesssim (1 + \tau^{-1}h)(\|\hat{\theta}_h^\tau\|_{\#} + \|\delta\|_{J \times \bar{\omega}} + t_h^\tau(I_T^\tau(u), I_T^\tau(u))^{\frac{1}{2}}). \end{aligned} $$

Altogether, this gives

$$ \begin{aligned} &\|L(\tilde{u}_T^\tau - u)\|_{L^2(J; H^{-1}(\Omega))} \\ &\lesssim (1 + \tau^{-1}h)(\|\hat{\theta}_h^\tau\|_{\#} + \|\delta\|_{J \times \bar{\omega}} + t_h^\tau(I_T^\tau(u), I_T^\tau(u))^{\frac{1}{2}}) + \|\Theta_f\|_{J \times \Omega}. \end{aligned} $$

(iv) Combining the bounds from the above steps proves (44). □

## 3.7 Weighting the regularization

The last step in our error analysis is to identify the weighting coefficient $c(\tau, h)$ in the regularization bilinear form. The interpolation operator defined in Section 3.3 fulfills the following convergence properties (the proof is outlined in Section 3.8.2).

**Lemma 8 (Approximation).** The following holds for all $n \in \{1:N\}$, all $T \in \mathcal{T}_h$, and all $v \in H^{\ell+1}(I_n; H^2(T)) \cap H^1(I_n; H^{k+1}(T))$,

$$ \|I_T^n(v) - v\|_{I_n \times T} \lesssim \tau^{\ell+1} \|v\|_{H^{\ell+1}(I_n; L^2(T))} + h^{k+1} \|v\|_{H^1(I_n; H^{k+1}(T))}, \quad (46) $$

$$ \|\nabla(I_T^n(v) - v)\|_{I_n \times T} \lesssim \tau^{\ell+1} \|v\|_{H^{\ell+1}(I_n; H^1(T))} + h^k \|v\|_{H^1(I_n; H^{k+1}(T))}, \quad (47) $$

$$ h^{\frac{1}{2}} \|\nabla(I_T^n(v) - v)\|_{I_n \times \partial T} \lesssim h^{\frac{1}{2}} \tau^{\ell+1} \|v\|_{H^{\ell+1}(I_n; H^2(T))} + h^k \|v\|_{H^1(I_n; H^{k+1}(T))}, \quad (48) $$

$$ h \|\Delta(I_T^n(v) - v)\|_{I_n \times T} \lesssim h \tau^{\ell+1} \|v\|_{H^{\ell+1}(I_n; H^2(T))} + h^k \|v\|_{H^1(I_n; H^{k+1}(T))}, \quad (49) $$

$$ \|\partial_t(I_T^n(v) - v)\|_{I_n \times T} \lesssim \tau^\ell \|v\|_{H^{\ell+1}(I_n; L^2(T))} + h^{k+1} \|v\|_{H^1(I_n; H^{k+1}(T))}. \quad (50) $$

Moreover, we also have for all $s \in \bar{I}_n$,

$$ \|(I_T^n(v) - v)(s)\|_T \lesssim \tau^{\ell+\frac{1}{2}} \|v\|_{H^{\ell+1}(I_n; L^2(T))} + h^{k+1} \|v\|_{H^1(I_n; H^{k+1}(T))}, \quad (51) $$

$$ \|\nabla(I_T^n(v) - v)(s)\|_T \lesssim \tau^{\ell+\frac{1}{2}} \|v\|_{H^{\ell+1}(I_n; H^1(T))} + h^k \|v\|_{H^1(I_n; H^{k+1}(T))}. \quad (52) $$

Finally, we have

$$ \begin{aligned} &h^{-\frac{1}{2}} \|I_T^n(v) - v\|_{I_n \times \partial T} + h^{-\frac{1}{2}} \|I_{\partial T}^n(v) - v\|_{I_n \times \partial T} \\ &\lesssim h^{-\frac{1}{2}} \tau^{\ell+1} \|v\|_{H^{\ell+1}(I_n; H^1(T))} + h^k \|v\|_{H^1(I_n; H^{k+1}(T))}. \end{aligned} \quad (53) $$

We consider the functional space

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

equipped with its natural norm denoted by $\|\cdot\|_*$.

$$ \|\hat{\theta}_h^\tau\|_{\#} + t_h^\tau(I_T^\tau(u), I_T^\tau(u))^{\frac{1}{2}} + \|\Theta_f\|_{J \times \bar{\omega}} \lesssim ((1 + h^{-1}\tau)^{\frac{1}{2}}\tau^{\ell+\frac{1}{2}} + h^k + \gamma^{\frac{1}{2}}c(\tau, h))\|u\|_*, \quad (55) $$

where we used that $\tau \lesssim 1$ and $h \lesssim 1$ to simplify the expression. Notice also that $u \in U_*$ implies that $f \in H^{\ell+1}(J; L^2(\Omega))$.
**Theorem 9 (Decay rates).** Assume that $u \in U_*$ and that $\tau \lesssim h, h \lesssim \tau$. Set

$$c(\tau, h) := \tau^{\ell+\frac{1}{2}} + h^k. \qquad (56)$$

We have the following decay rates on the errors:

$$\|\hat{e}_h^\tau, \hat{\xi}_h^\tau\| + t_h^\tau(u_\tau^\tau, u_\tau^\tau)^{\frac{1}{2}} + s_h^\tau(\hat{u}_h^\tau, \hat{u}_h^\tau)^{\frac{1}{2}} \lesssim (\tau^{\ell+\frac{1}{2}} + h^k)\|u\|_* + \|\delta\|_{J \times \bar{\omega}}, \qquad (57)$$

$$\|u - u_{\tau}\|_{\mathrm{trg}} \lesssim C_{\mathrm{stb}}(\tau^{\ell+\frac{1}{2}} + h^k)^{\alpha} (\|u\|_* + (\tau^{\ell+\frac{1}{2}} + h^k)^{-1} \|\delta\|_{J \times \bar{\omega}}). \qquad (58)$$

where the hidden constant scales as $\max(\gamma, 1)^{\frac{1}{2}}$ in (57) and $\max(\gamma, \gamma^{\alpha-1})^{\frac{1}{2}}$ in (58).

*Proof.* Plug (56) into (55) and invoke Lemma 5 for (57) and Theorem 7 for (58). $\square$

**Remark 3.2. (Choosing discretization parameters)** The estimates of Theorem 9 indicate that space and time refinements have to be stopped when $(\tau^{\ell+\frac{1}{2}} + h^k)\|u\|_* \simeq \|\delta\|_{J \times \bar{\omega}}$. Hence, as the noise level diminishes, finer discretizations can be employed.

## 3.8 Technical proofs

In this section, we outline the proofs of Lemma 2 and Lemma 8. The proofs use standard arguments from finite element approximation theory.

### 3.8.1 Proof of Lemma 2

Let $\hat{v}_h^\tau := (v_\tau^\tau, v_\tau^\tau) \in \tilde{U}_h^\tau$.

(i) Proof of (27). For all $n \in \{1:N\}$ and all $q_n^\tau \in \mathbb{P}^\ell(I_n; U_\tau^k)$, we have

$$(\partial_t(R^n(v_\tau^\tau) - v_\tau^n), q_\tau^n)_{I_n \times \Omega} = (\llbracket v_\tau^\tau \rrbracket^{n-1}, q_\tau^n(t_{n-1}^+))_{\Omega}.$$

Invoking a discrete trace inequality in time implies

$$\|\partial_t(R^n(v_\tau^\tau) - v_\tau^n)\|_{I_n \times \Omega} \lesssim \tau^{-\frac{1}{2}} \|\llbracket v_\tau^\tau \rrbracket^{n-1}\|_{\Omega}. \qquad (59)$$

Summing over $n \in \{1:N\}$ proves (27) since $\sum_{n \in \{1:N\}} \|\llbracket v_\tau^\tau \rrbracket^{n-1}\|_{\Omega}^2 \le d^\tau(v_\tau^\tau, v_\tau^\tau)$ (recall that $\llbracket \cdot \rrbracket^0 = 0$ by convention).

(ii) Proof of (28). Since we have $(R^n(v_\tau^\tau) - v_\tau^n)(t_n^-) = 0$ (see (26)), a Poincaré inequality in time over $I_n$ and the bound (59) yield

$$\|R^n(v_\tau^\tau) - v_\tau^n\|_{I_n \times \Omega} \lesssim \tau \|\partial_t(R^n(v_\tau^\tau) - v_\tau^n)\|_{I_n \times \Omega} \lesssim \tau^{\frac{1}{2}} \|\llbracket v_\tau^\tau \rrbracket^{n-1}\|_{\Omega}.$$

Summing over $n \in \{1:N\}$ proves as above that

$$\|R^\tau(v_\tau^\tau) - v_\tau^\tau\|_{J \times \Omega} \lesssim \tau^{\frac{1}{2}} d^\tau(v_\tau^\tau, v_\tau^\tau)^{\frac{1}{2}}. \qquad (60)$$

Moreover, (25) shows that the operators $R^\tau$ and $\nabla_\tau$ commute. Hence, using the same arguments as above shows that

$$\|\nabla_\tau(R^n(v_\tau^\tau) - v_\tau^n)\|_{I_n \times \Omega} \lesssim \tau^{\frac{1}{2}} \|\llbracket \nabla_\tau v_\tau^\tau \rrbracket^{n-1}\|_{\Omega}.$$
(ii) Proof of (47), (48), (49), and (50). We can use similar arguments to those invoked in Step (i) since we have

$$ \|\nabla(I_T^n(v) - v)\|_{I_n \times T} \le \|\nabla(\Pi_T^k(\tilde{I}_n^\ell(v)) - \tilde{I}_n^\ell(v))\|_{I_n \times T} + \|\tilde{I}_n^\ell(\nabla v) - \nabla v\|_{I_n \times T}, $$

$$ \|\partial_t(I_T^n(v) - v)\|_{I_n \times T} \le \|\Pi_T^k(\partial_t \tilde{I}_n^\ell(v)) - \partial_t \tilde{I}_n^\ell(v)\|_{I_n \times T} + \|\partial_t(\tilde{I}_n^\ell(v) - v)\|_{I_n \times T}, $$

because the operators $\tilde{I}_n^\ell$ and $\nabla$ commute, as well as the operators $\Pi_T^k$ and $\partial_t$. A similar argument can also be used when replacing $\nabla$ by $\Delta$.

(iii) The proofs of (51), (52), and (53) employ similar arguments.

# 4 Numerical results

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

The space domain is $\Omega := (0, 1)^d$ with $d \in \{1, 2\}$. All the errors are computed as the difference between the numerical solution and the $L^2(J; L^2(\Omega))$-orthogonal projection of the exact solution. These errors are measured in the $L^2(T_1, T_2; H^1(B))$-seminorm, with the subset $(T_1, T_2) \times B$ specified below for each test case. Some noise is added to the measurements in the following way: (i) the time-space domain is divided in $10^{d+1}$ subdomains; (ii) a random noise level $\delta := a * \text{rand()}$ is assigned to each subdomain, with $a \ge 0$ the noise amplitude and $\text{rand()}$ is a C++ function returning a random number in $[-1, 1]$. Thus, every subdomain has the same noise during the whole mesh-refinement process. Moreover, the Tikhonov regularization coefficient is set to $\gamma := 10^{-3}$.

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

## 4.1 One-dimensional test cases

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

$$ \Omega := (0, 1), \quad \varpi := (0.25, 0.75), \quad T_f := 2, \quad u(t, x) := \cos(\pi t) \sin(\pi x). $$

All the errors are estimated in the subdomain

$$ (T_1, T_2) \times B := (0.2, 1.8) \times (0.125, 0.875). $$

We consider two settings. First, we address the case with unknown initial data but with known Dirichlet boundary data on $J \times \partial\Omega$. In this situation, the boundary conditions can be enforced in a strong way on the boundary face degrees of freedom by searching $\hat{u}_h^\tau$ in $\tilde{U}_{h0}^\tau$ and not in $\tilde{U}_h^\tau$. Notice that in this case, we have $\alpha = 1$ in Lemma 1 (see Theorem 2 from [8]). We then expect the same convergence rates as for a well-posed problem, i.e., $k$ in space and $\ell + \frac{1}{2}$ in time. Then, we consider the case where both initial and boundary data are unknown.
<figure><img src="image_3.png" /><figcaption>Figure 1: $L^2(T_1, T_2; H^1(B))$-seminorm errors with respect to the number of dofs for the 1d test case without any noise. Both cases with known or unknown boundary data are considered. Left: space convergence ($M \in \{16, 32, 64, 128\}$, $N = 128$, $\ell = 3$). Right: time convergence ($N \in \{10, 20, 40, 80\}$, $M = 256$, $k = 3$).</figcaption></figure>

We use a uniform mesh in space and in time ($N$ cells in time and $M$ cells in space). Four levels of refinement are considered in space ($M \in \{16, 32, 64, 128\}$) and in time ($N \in \{10, 20, 40, 80\}$). We run two convergence studies. At first, for a good precision in time ($N = 128$, $\ell = 3$), we study the error for several successive space mesh sizes $M \in \{16, 32, 64, 128\}$ and several polynomial orders $k \in \{1, 2, 3\}$. Then, for a good precision in space ($M = 256$, $k = 3$), we study the error for several successive time mesh sizes $N \in \{10, 20, 40, 80\}$ and several polynomial orders $\ell \in \{0, 1, 2\}$. The results without any noise are reported in Figure 1. We observe optimal space convergence at rate $k$ for $k \in \{2, 3\}$. Instead, we obtain superconvergence at rate 2 for $k = 1$. Moreover, we observe time convergence at rate $\ell + 1$, for $\ell \in \{1, 2\}$, which is slightly better than the expected rate $\ell + \frac{1}{2}$ from Theorem 9. (Recall that the errors are computed as the difference between the numerical solution and a projection of the exact solution.) The errors in Figure 1 are reported as a function of the total number of space-time degrees of freedom (dofs). This gives a perspective on the efficiency of the method with respect to its cost. In particular, we see that, for a given number of dofs, we can reach a better precision by employing a high-order method. We also notice that the errors with known boundary data are always smaller than those with unknown boundary data. For completeness, the number of dofs is reported in Table 1. The numbers correspond to the case of known boundary data; those for unknown boundary data are slightly higher due to the use of additional boundary unknowns.

Let us now consider some noise in the measurements. In this situation, we obtain the same results as above until the discrete error becomes small and then the convergence stops. To illustrate this point, let us focus on the highest polynomial degree (which generates the lowest error). The results are reported in Figure 2 for $\ell = 2$ and $k = 3$. As expected, the higher the noise, the larger the value at which the error stagnates. More precisely, in the case of unknown boundary data, the error stagnates at about twice the level of noise ($2 \cdot 10^{-5}$ for $a = 10^{-5}$ and $2 \cdot 10^{-3}$ for $a = 10^{-3}$), whereas it stagnates at
<table><caption>Table 1: Number of degrees of freedom (1d test case with known boundary data). Left: space refinement for $l = 3$ and $N = 128$. Right: time refinement for $M = 256$ and $k = 3$.</caption><thead><tr><th>M</th><th>k = 1</th><th>k = 2</th><th>k = 3</th><th>N</th><th>l = 0</th><th>l = 1</th><th>l = 2</th></tr></thead><tbody><tr><td>16</td><td>4.8·10<sup>4</sup></td><td>6.5·10<sup>4</sup></td><td>8.1·10<sup>4</sup></td><td>10</td><td>2.6·10<sup>4</sup></td><td>5.1·10<sup>4</sup></td><td>7.7·10<sup>4</sup></td></tr><tr><td>32</td><td>9.7·10<sup>4</sup></td><td>1.3·10<sup>5</sup></td><td>1.6·10<sup>5</sup></td><td>20</td><td>5.1·10<sup>4</sup></td><td>1.0·10<sup>5</sup></td><td>1.5·10<sup>5</sup></td></tr><tr><td>64</td><td>2.0·10<sup>5</sup></td><td>2.6·10<sup>5</sup></td><td>3.3·10<sup>5</sup></td><td>40</td><td>1.0·10<sup>5</sup></td><td>2.0·10<sup>5</sup></td><td>3.1·10<sup>5</sup></td></tr><tr><td>128</td><td>3.9·10<sup>5</sup></td><td>5.2·10<sup>5</sup></td><td>6.5·10<sup>5</sup></td><td>80</td><td>2.0·10<sup>5</sup></td><td>4.1·10<sup>5</sup></td><td>6.1·10<sup>5</sup></td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Figure 2: $L^2(T_1, T_2; H^1(B))$-seminorm errors with respect to the number of dofs for the 1d test case and various noise levels. Left: space convergence ($M \in \{16, 32, 64, 128\}$, $k = 3$, $N = 128$, $l = 3$). Right: time convergence ($N \in \{10, 20, 40, 80\}$, $l = 2$, $M = 256$, $k = 3$).</figcaption></figure>

$8 \cdot 10^{-5}$ for $a = 10^{-3}$ in the case of known boundary data. This means that the test case with unknown boundary data is more sensitive to the presence of noise.

## 4.2 Two-dimensional test case

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

$$ \Omega := (0, 1)^2, \quad \varpi := \Omega \setminus (0, 0.875) \times (0.125, 0.875), \quad T_f := 2, \\ (T_1, T_2) \times B := (0.2, 1.8) \times (0.125, 0.875)^2, \quad u(t, x, y) := \cos(\pi t) \sin(\pi x) \sin(\pi y). $$

Notice that $\bar{B} \cap \partial\Omega = \emptyset$ and that $\varpi \not\subset B$. Both cases with and without boundary data are considered.

The space convergence is studied on a sequence of seven triangulations of increasing refinement ($h \in \{0.25, 0.18, 0.125, 0.09, 0.0625, 0.045, 0.03125\}$) and $k \in \{1, 2\}$. The first, third, fifth, and seventh triangulations of the sequence are shown in Figure 3. Notice that all the triangulations are fitted to the sets $\varpi$ and $B$. We use $N = 40$ time intervals and the time degree $l = 2$. The time convergence is studied using $N \in \{5, 10, 15, 20, 25, 30\}$, $l \in \{0, 1\}$, and the space discretization using the finest time mesh and $k = 2$.

The $L^2(T_1, T_2; H^1(B))$-seminorm errors obtained without any noise are reported in Figure 4 as a function of the total number of space-time dofs. The convergence in space
<figure><img src="image_5.png" /><figcaption>Figure 3: First, third, fifth, and seventh triangulations used for the 2d test case; all the triangulations are fitted to $\bar{\omega}$ and $B$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: $L^2(T_1, T_2; H^1(B))$-seminorm errors for the 2d test case without any noise. Left: space convergence ($N = 40$ except the curve with red bullets for which $N = 80$, $\ell = 1$). Right: time convergence (finest mesh, $k = 2$).</figcaption></figure>

(left panel) indicates that the scheme has a convergence rate of $k$ with and without the knowledge of boundary data. (Since the number of dofs scales as $h^{-d}$, the slopes in the left panel of Figure 4 have to be multiplied by two to get the convergence rate.) This convergence stops in the case of unknown boundary data when we reach the lowest error enabled by the time refinement (this error can be extrapolated using the curve corresponding to $\ell = 1$ on the right panel). When this lowest error is reached, continuing refining the mesh increases the error. To confirm that the obstruction comes from the time discretization error, we also report in the left panel of Figure 4 the space convergence error obtained with a finer time discretization ($N = 80$). We also notice that, since the convergence rate for $k = 1$ is indeed one in the case of known boundary data, we conjecture that the superconvergence observed in the 1d test case was due to the fact that the domain was one-dimensional. The convergence in time (right panel of Figure 4) seems to correspond to a convergence at rate $\ell + 1/2$, even if the boundary data are not known. (Recall that for the one-dimensional case, we got $\ell + 1$ for $\ell \in \{1, 2\}$.) An interesting observation is that, in the case of known boundary data, space refinement has to be finer since the error stops converging on the right panel for $\ell = 1$ (which indicates a lack of space refinement), whereas, for unknown boundary data, convergence stops on the left panel for
<table><caption>Table 2: Number of degrees of freedom (2d test case with known boundary value). Left: space refinement for $ℓ = 1$ and $N = 40$. Right: time refinement for $h = 0.03125$ and $k = 2$.</caption><thead><tr><th>h</th><th>k = 1</th><th>k = 2</th><th>N</th><th>ℓ = 0</th><th>ℓ = 1</th></tr></thead><tbody><tr><td>2.5·10<sup>-1</sup></td><td>6.0·10<sup>4</sup></td><td>1.1·10<sup>5</sup></td><td>5</td><td>2.6·10<sup>5</sup></td><td>5.2·10<sup>5</sup></td></tr><tr><td>1.8·10<sup>-1</sup></td><td>1.2·10<sup>5</sup></td><td>2.0·10<sup>5</sup></td><td>10</td><td>5.2·10<sup>5</sup></td><td>1.0·10<sup>6</sup></td></tr><tr><td>1.3·10<sup>-1</sup></td><td>1.9·10<sup>5</sup></td><td>3.3·10<sup>5</sup></td><td>15</td><td>7.8·10<sup>5</sup></td><td>1.6·10<sup>6</sup></td></tr><tr><td>9.0·10<sup>-2</sup></td><td>3.6·10<sup>5</sup></td><td>6.2·10<sup>5</sup></td><td>20</td><td>1.0·10<sup>6</sup></td><td>2.1·10<sup>6</sup></td></tr><tr><td>6.3·10<sup>-2</sup></td><td>5.5·10<sup>5</sup></td><td>9.7·10<sup>5</sup></td><td>25</td><td>1.3·10<sup>6</sup></td><td>2.6·10<sup>6</sup></td></tr><tr><td>4.5·10<sup>-2</sup></td><td>1.3·10<sup>6</sup></td><td>2.2·10<sup>6</sup></td><td>30</td><td>1.6·10<sup>6</sup></td><td>3.1·10<sup>6</sup></td></tr><tr><td>3.1·10<sup>-2</sup></td><td>2.4·10<sup>6</sup></td><td>4.2·10<sup>6</sup></td><td>.</td><td>.</td><td>.</td></tr></tbody></table>

$k = 2$, (which indicates a lack of time refinement). For completeness, the number of dofs is reported in Table 2. We see that we reach several millions of dofs, which is high for the use of a direct solver like Pardiso. In fact, these numerical simulations use more than 100GB of RAM.

<figure><img src="image_7.png" /><figcaption>Figure 5: $L^2(T_1, T_2; H^1(B))$-seminorm errors for the 2d test case with noised data ($a \in \{0, 0.05, 0.1\}$). Left: space convergence ($k = 2$ and $ℓ = 1$). Right: time convergence (finest mesh from Figure 3, $k = 2$ and $ℓ = 1$).</figcaption></figure>

Finally, we consider some noise at level $a \in \{0, 0.05, 0.1\}$. The results are reported in Figure 5 along with the solution without any noise for $k = 2$ and $ℓ = 1$. In the left panel, $N = 40$ (resp. $N = 80$) is considered for known (resp. unknown) boundary data. We observe once more that the test case without the knowledge of boundary data is much more sensitive to the presence of noise. Indeed, the presence of noise changes the result only in this case.

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

[1] G. Alessandrini, L. Rondi, E. Rosset, and S. Vessella. The stability for the Cauchy problem for elliptic equations. *Inverse Problems*, 25(12):123004, 47, 2009.

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

[3] L. Bourgeois. A mixed formulation of quasi-reversibility to solve the Cauchy problem for Laplace's equation. *Inverse Problems*, 21(3):1087–1104, 2005.

[4] L. Bourgeois. Convergence rates for the quasi-reversibility method to solve the Cauchy problem for Laplace's equation. *Inverse Problems*, 22(2):413–430, 2006.

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

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

[7] E. Burman and A. Ern. Continuous interior penalty *hp*-finite element methods for advection and advection-diffusion equations. *Math. Comp.*, 76(259):1119–1140, 2007.

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

[9] E. Burman, P. Hansbo, and M. G. Larson. Solving ill-posed control problems by stabilized finite element methods: an alternative to Tikhonov regularization. *Inverse Problems*, 34(3):035004, 36, 2018.

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

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

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

[13] E. Burman, G. Delay, A. Ern, and L. Oksanen. A stability estimate for data assimilation subject to the heat equation with initial datum, 2022. To appear in *Comptes Rendus Math. (CRAS)*.

[14] M. Cicuttin, D. A. Di Pietro, and A. Ern. Implementation of Discontinuous Skeletal methods on arbitrary-dimensional, polytopal meshes using generic programming. *J. Comput. Appl. Math.*, 2017.
[15] W. Dahmen, R. Stevenson, and J. Westerdiep. Accuracy controlled data assimilation for parabolic problems. *Math. Comp.*, 91(334):557-595, 2022.

[16] J. Dardé, A. Hannukainen, and N. Hyvönen. An $H_{div}$-based mixed quasi-reversibility method for solving elliptic Cauchy problems. *SIAM J. Numer. Anal.*, 51(4):2123-2148, 2013.

[17] O. Y. Èmanuilov. Controllability of parabolic equations. *Mat. Sb.*, 186(6):109-132, 1995.

[18] A. Ern and J.-L. Guermond. *Finite elements. I Approximation and Interpolation*, volume 72 of *Texts in Applied Mathematics*. Springer, Cham, 2021.

[19] A. Ern and J.-L. Guermond. *Finite elements III. First-order and time-dependent PDEs*, volume 74 of *Texts in Applied Mathematics*. Springer, Cham, 2021.

[20] K. Ito and B. Jin. *Inverse problems. Tikhonov theory and algorithms*, volume 22 of *Series on Applied Mathematics*. World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2015. ISBN 978-981-4596-19-0.

[21] S. E. Jenkins, C. J. Budd, M. A. Freitag, and N. D. Smith. The effect of numerical model error on data assimilation. *J. Comput. Appl. Math.*, 290:567-588, 2015.

[22] M. V. Klibanov. Estimates of initial conditions of parabolic equations and inequalities via lateral Cauchy data. *Inverse Problems*, 22(2):495-514, 2006.

[23] J. Korevaar and J. L. H. Meyers. Logarithmic convexity for supremum norms of harmonic functions. *Bull. London Math. Soc.*, 26(4):353-362, 1994.

[24] L. Robbiano. Théorème d'unicité adapté au contrôle des solutions des problèmes hyperboliques. *Comm. Partial Differential Equations*, 16(4-5):789-800, 1991.

[25] I. Smears. Robust and efficient preconditioners for the discontinuous Galerkin time-stepping method. *IMA J. Numer. Anal.*, 37(4):1961-1985, 2017.

[26] M. Tadi, M. V. Klibanov, and W. Cai. An inversion method for parabolic equations based on quasireversibility. *Comput. Math. Appl.*, 43(8-9):927-941, 2002.

[27] S. Vessella. Unique continuation properties and quantitative estimates of unique continuation for parabolic equations. In *Handbook of differential equations: evolutionary equations. Vol. V*, Handb. Differ. Equ., pages 421-500. Elsevier/North-Holland, Amsterdam, 2009.

[28] Y. B. Wang, J. Cheng, J. Nakagawa, and M. Yamamoto. A numerical method for solving the inverse heat conduction problem without initial value. *Inverse Probl. Sci. Eng.*, 18(5):655-671, 2010.