# A posteriori error analysis and adaptivity of a space-time finite element method for the wave equation in second order formulation

Zhaonan Dong, Emmanuil H Georgoulis, Lorenzo Mascotto, Zuodong Wang

► To cite this version:

Zhaonan Dong, Emmanuil H Georgoulis, Lorenzo Mascotto, Zuodong Wang. A posteriori error analysis and adaptivity of a space-time finite element method for the wave equation in second order formulation. 2025. hal-05248777

HAL Id: hal-05248777

https://inria.hal.science/hal-05248777v1

Preprint submitted on 10 Sep 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# *A posteriori* error analysis and adaptivity of a space-time finite element method for the wave equation in second order formulation

Zhaonan Dong\*, Emmanuil H. Georgoulis\d, Lorenzo Mascotto\dd, and Zuodong Wang\sp

## Abstract

We establish rigorous *a posteriori* error bounds for a space-time finite element method of arbitrary order discretising linear wave problems in second order formulation. The method combines standard finite elements in space and continuous piecewise polynomials in time with an upwind discontinuous Galerkin-type approximation for the second temporal derivative. The proposed scheme accepts dynamic mesh modification, as required by space-time adaptive algorithms, resulting in a discontinuous temporal discretisation when mesh changes occur. We prove *a posteriori* error bounds in the $L^\infty(L^2)$-norm, using carefully designed temporal and spatial reconstructions; explicit control on the constants (including the spatial and temporal orders of the method) in those error bounds is shown. The convergence behaviour of an error estimator is verified numerically, also taking into account the effect of the mesh change. A space-time adaptive algorithm is proposed and tested numerically.

**Keywords.** Wave equation; space-time methods; finite element method; a posteriori error analysis; adaptive algorithm.

**MSC.** 65N30, 65M50; 65M60; 65J10.

## 1 Introduction

Solutions to hyperbolic equations are often characterised by highly localised spatio-temporal features, such as travelling fronts, pulses, and their interactions. Their numerical approximation over large domains and for long times is, therefore, particularly challenging: standard numerical methods, which are based on quasi-uniform time-steps and spatial mesh-sizes, should employ spatio-temporal grids/meshes that are sufficiently fine to capture localised features that may result in very high numbers of degrees of freedom. To address the computational complexity, wave simulations are typically performed using explicit time-stepping methods, possibly combined with modern computer architectures, e.g., implementation on graphics processing units, etc.

Explicit schemes require CFL-type time-step size restrictions to ensure stable discretisations, which are typically further compounded when non-linearities are present. Aiming to remove the CFL-type restrictions and to provide localised spatio-temporal resolution, implicit space-time finite element methods have been developed, since the late 1980s [36, 37, 40, 24, 55, 25] and has continued to this day in various settings [41, 34, 54, 17, 7, 45, 42, 45, 8, 39, 49, 50, 23, 22, 6]. Of particular interest to the present work is the method presented in [54], employing discrete spaces of continuous piecewise polynomials in time, tensorised

<footnote>\*Inria, 48 Rue Barrault, 75647 Paris, France and CERMICS, Ecole nationale des ponts et chaussées, IP Paris, 6 \& 8 avenue B. Pascal, 77455 Marne-la-Vallée, France (zhaonan.dong@inria.fr)</footnote>

<footnote>\dThe Maxwell Institute for Mathematical Sciences and Department of Mathematics, Heriot-Watt University, Edinburgh EH14 4AS, United Kingdom, and Department of Mathematics, School of Applied Mathematical and Physical Sciences, National Technical University of Athens, Zografou 15780, Greece, and IACM-FORTH, Greece. (E.Georgoulis@hw.ac.uk)</footnote>

<footnote>\ddDipartimento di Matematica e Applicazioni, Università degli Studi di Milano-Bicocca, 20125 Milan, Italy; IMATI-CNR, 27100, Pavia, Italy; Fakultät für Mathematik, Universität Wien, 1090 Vienna, Austria (lorenzo.mascotto@unimib.it)</footnote>

<footnote>\spCERMICS, Ecole nationale des ponts et chaussées, IP Paris, 6 & 8 avenue B. Pascal, 77455 Marne-la-Vallée, France and Inria, 48 Rue Barrault, 75647 Paris, France (zuodong.wang@inria.fr)</footnote>
with $H^1$-conforming finite element functions in space over static spatial simplicial and box-type meshes; a corresponding scheme for a nonlinear wave problem is given in [30].

From their inception, a central argument to alleviate the additional computational overhead of implicit space-time methods has been their ability to accept locally variable meshes [36, 37, 40]. However, typically the local resolution requirements are not known a priori. One way to address this issue is the derivation of a posteriori error estimates or a posteriori error bounds in various norms. Rigorous a posteriori error analysis is very well developed in the context of elliptic (see, e.g., [53] and the references therein) and parabolic problems (see [19, 43, 53] and the reference therein); fewer results are available for wave problems. In [40], Johnson proved an a posteriori error bound based on a duality argument for a space-time method employing discontinuous Galerkin time-stepping and finite elements in space. An a posteriori error indicator via super-convergence points is presented in [1]. Goal-oriented error estimation and adaptivity frameworks are presented in [3, 4, 5]. In [9] a posteriori error bounds in the $L^\infty(H^1)$- and final time $L^2$-norms for a backward Euler time-stepping with finite elements in space are proven employing energy arguments. For the same method, $L^\infty(L^2)$-norm a posteriori error bounds are proven in [27] employing appropriate space and time reconstructions. Bounds for time discretisations by the leap-frog and cosine methods are proven in [28]. In [12], a posteriori error bounds in the damped energy norm $(H^1(L^2) + L^2(H^1))$ for the leap-frog scheme combined with a finite element spatial discretization without mesh change are proven; see also [11] for earlier developments on error estimators for spatially-discrete schemes. For the Newmark scheme a posteriori error bounds in the $L^\infty(H^1)$-norm are proven in [31]. A local time-stepping strategy combined with the leap-frog method is considered in [14]; the corresponding a posteriori bounds are given in [33]. Also, related to the developments above, the a posteriori error bounds in the $L^\infty(L^2)$-norm for the semidiscrete time-stepping method in [54] is presented in [15]. Remarkably, with the exception of the early work [40], we are not aware of rigorous a posteriori error analysis for any other arbitrary order space-time Galerkin-type method.

This work is concerned with the proof of fully computable a posteriori error bounds in the $L^\infty(L^2)$-norm for a variant of the method in [54], allowing, crucially, for dynamic mesh modification at each time-step; the results extend the respective time-discrete scheme a posteriori bound from [15] to the practical space-time setting with mesh modification. The choice of the method is motivated by the desire to strike a “cost-analysis” balance: this is an arbitrary order one-field discretisation, resulting in fewer degrees of freedom than methods defined upon rewriting the problem as a first order system, and allowing for dynamic mesh modification yet retaining the possibility of deriving rigorous error control. It is known already since the 1970s [2] that the $L^\infty(H^1)$-norm error for standard Galerkin approximations of the wave equation is susceptible to initial data approximation errors. Hence, dynamic mesh modification typically creates a compounded effect, resulting in suboptimal convergence rates for the $L^\infty(H^1)$-norm error [55, 41].

The proof of the $a posteriori L^\infty(L^2)$-norm error upper bound for the method in [54] with dynamic mesh modification requires a combination of known and new reconstructions. More specifically, the mesh modification creates discontinuity across time nodes, which is handled using averaging of a Hermite interpolation basis expansion to ensure the preservation of the first temporal derivative while simultaneously lifting the temporal discontinuity; this non-standard reconstruction is instrumental in keeping the error proliferation due to mesh modification under control. The time derivative of the resulting lifted continuous-in-time field is further reconstructed using the discontinuous Galerkin time-stepping reconstruction of Makridakis and Nochetto [44], combined with a non-standard elliptic reconstruction, motivated by the developments in [29] for parabolic problems. Special care is taken to provide values (or bounds) for all constants involved in the corresponding estimates.

Further, we validate the theoretical results by detailed numerical experiments; we also investigate the energy dissipation of the method for different approximation orders to find that dissipation reduces while increasing the order in time of the method. We assess the resulting a posteriori error estimates on test cases with smooth and singular solutions, with and without mesh modification. Moreover, we propose a practical space-time adaptive algorithm employing the constituent computable error estimators of the new a posteriori $L^\infty(L^2)$-norm error bound; a series of numerical experiments of dynamically adapted meshes showcases the relevance of the proven a posteriori bound in driving mesh modification.

The remainder of this work is structured as follows. In Section 2 we present the model problem and
provide some function space notation. The method is introduced in Section 3. The definition and the properties of the reconstructions required for the proof of the main result are discussed in Section 4. The proof of the a posteriori error estimate is given in Section 5. Section 6 contains the numerical experiments and the description of the space-time adaptive algorithm. Finally in Section 7, we draw some conclusions.

## 2 Model problem

Let $D$ be a polytopic Lipschitz domain in $\mathbb{R}^d$, $d = 1, 2, 3$. Standard notation is used throughout for Sobolev and Bochner spaces: the space of Lebesgue measurable and square integrable functions over $D$ is denoted by $L^2(D)$, while the Sobolev space of positive integer order $s$ by $H^s(D)$. $H_0^1(D)$ is the space of functions in $H^1(D)$ with zero trace over the boundary of $D$. The space $H^{-1}(D)$ is the dual of $H_0^1(D)$ and the duality pairing between these two spaces is $\langle \cdot, \cdot \rangle$.

Given $\mathcal{X}$ a real Banach space with norm $\|\cdot\|_{\mathcal{X}}$, an interval $I$, and $k$ larger than or equal to 1, we define $L^k(I; \mathcal{X})$ as the Bochner space of measurable functions $w$ from $I$ to $\mathcal{X}$ such that the following norm is finite:

$$ \|w\|_{L^k(I; \mathcal{X})} := \left( \int_I \|w(t)\|_{\mathcal{X}}^k dt \right)^{\frac{1}{k}} \quad \text{for } 1 \le k < \infty; \qquad \|w\|_{L^\infty(I; \mathcal{X})} := \ess \sup_{t \in I} \|w(t)\|_{\mathcal{X}}. $$

For any positive integer $s$ and any $k \ge 1$ the space $W^{s,k}(I; \mathcal{X})$ is the space of measurable functions $w$ whose derivatives in time up to order $s$ belong to $L^k(I; \mathcal{X})$.

Let $T$ be a positive final time and $J := (0, T]$ the evolution time interval. We use $'$ and $''$ to denote the first- and second-order time derivative, respectively. For spatial derivatives, we drop the spatial dependence for brevity, e.g., $\nabla := \nabla_x$. The spatial $L^2(D)$-inner product is abbreviated to $(\cdot, \cdot)$, the corresponding norm to $\|\cdot\|$. We omit the dependence in time and space for inner products, duality pairings, and norms when they are defined on the whole spatial domain $D$ or the whole time interval $J$, e.g., $\|\cdot\|_{L^2(L^2)} := \|\cdot\|_{L^2(J; L^2(D))}$.

Let $u_0 \in H_0^1(D)$, $u_1 \in L^2(D)$, and $f \in L^2(J; L^2(D))$. We henceforth assume that the spatial domain $D$ is convex<fnref n="1" />. We consider the linear wave initial/boundary value problem: find $u: J \times D \to \mathbb{R}$, such that

$$ \left\{ \begin{array}{ll} u'' - \Delta u = f & \text{in } J \times D \\ u = 0 & \text{on } J \times \partial D \\ u(0) = u_0, \quad u'(0) = u_1 & \text{in } D, \end{array} \right. \qquad (2.1) $$

in weak form: find $u \in H^2(J; H^{-1}(D)) \cap L^2(J; H_0^1(D)) \cap H^1(J; L^2(D))$, such that

$$ \left\{ \begin{array}{ll} \displaystyle \int_J (\langle u'', w \rangle + (\nabla u, \nabla w)) \, dt = \int_J (f, w) \, dt & \text{for all } w \in L^2(J; H_0^1(D)) \\ u(0) = u_0, \quad u'(0) = u_1. & \end{array} \right. \qquad (2.2) $$

Problem (2.2) is well-posed [46, Theorem 8.2-2]. With minor modifications to the analysis below, variable coefficients, other linear spatial differential operators, and different boundary conditions are possible to be treated.

## 3 A space-time finite element method

For $N \in \mathbb{N}$, we consider a grid $0 = t_0 < t_1 < \dots < t_N = T$ of $J$ and introduce the time-step $\tau_n := t_n - t_{n-1}$ for all $n \in \mathcal{N} := \{1, 2, \dots, N\}$. With each time interval $I_n := (t_{n-1}, t_n]$, we associate a local polynomial degree $q_n$; we collect these in the polynomial degree vector $q := (q_1, q_2, \dots, q_N)$. For $k \in \mathbb{Z}$, $q+k$ is the vector of entries $q_n+k$, $n \in \mathcal{N}$. For all $n \in \mathcal{N} \cup \{0\}$, we consider sequences of simplicial conforming meshes $T_h^n$ over the spatial domain $D$, set $T_h^0 := T_h^1$, assume that neighbouring elements have uniformly comparable size, and

<footnote n="1">¹This assumption is needed in Section 5.</footnote>
(for $d > 1$) the diameter of each element is uniformly comparable to that of each of its facets. For all cell $K$ in $\mathcal{T}_h^n$, we define the mesh-size $h_K$ as its diameter; $h_n$ is the piecewise constant function $h_{n|K}(\mathbf{x}) := h_K$ for all $n \in \mathcal{N}$, $K \in \mathcal{T}_h^n$, and $\mathbf{x}$ in the interior of $K$. The set of all interior facets of a mesh $\mathcal{T}_h^n$ is $\mathcal{F}_I^n$; the set of the facets of a cell $K$ is $\mathcal{F}_K$. We assume that two consecutive meshes are constructed by a uniformly bounded finite number of mesh refinements and coarsening steps.

With each cell $K$ and each facet $F$, we associate the outward normal vector $\mathbf{n}_K$ and the (fixed once-and-for-all amongst the two possible ones) normal vector $\mathbf{n}_F$. Given the broken Sobolev space

$$H^1(\mathcal{T}_h^n) := \{w \in L^2(D) \mid w|_K \in H^1(K) \quad \forall K \in \mathcal{T}_h^n\},$$

we define the spatial jump operator $[\cdot]$ as the identity on each boundary facet and, on each $F$ in $\mathcal{F}_I^n$ by

$$[w] := (w|_{K_1 \mathbf{n}_{K_1}} \cdot \mathbf{n}_F)|_F + (w|_{K_2 \mathbf{n}_{K_2}} \cdot \mathbf{n}_F)|_F \quad \forall w \in H^1(\mathcal{T}_h^n), \quad F = K_1 \cap K_2.$$

Moreover, for all $n \in \mathcal{N} \cup \{0\}$, $V_h^n$ is an $H^1$-conforming finite element space of uniform degree $p$ over a mesh $\mathcal{T}_h^n$. $\Pi_p^n$ is the orthogonal $L^2$-projector from $L^2(D)$ to $V_h^n$; its global counterpart $\Pi_p: L^2(J; L^2(D)) \to L^2(J; H_0^1(D))$ is given by $\Pi_p|_{I_n} := \Pi_p^n$. In addition, we consider its orthogonal complement projector $\Pi_p^\perp|_{I_n} := \Pi_p^{n,\perp} := I_d - \Pi_p^n$. Then, we introduce the intersection and union spaces

$$V_h^{n, \oplus} := V_h^n \cap V_h^{n-1} \quad \text{and} \quad V_h^{n, \oplus} := V_h^n \cup V_h^{n-1}, \quad \forall n \in \mathcal{N},$$

with associated meshes $\mathcal{T}_h^{n, \oplus}$ and $\mathcal{T}_h^{n, \ominus}$, interior facet sets $\mathcal{F}_I^{n, \oplus}$ and $\mathcal{F}_I^{n, \ominus}$, and mesh-size functions $h_n^\oplus$ and $h_n^\ominus$, respectively. We also define coarsened (resp. refined) cells as the cells that only belong to one between $\mathcal{T}_h^{n-1}$ and $\mathcal{T}_h^n$, and also belong to $\mathcal{T}_h^{n, \oplus}$ (resp. $\mathcal{T}_h^{n, \ominus}$), the subdomain $D_c^n$ (resp. $D_r^n$) as the union of cells sharing at least one vertex with the coarsened (resp. refined) cells between $\mathcal{T}_h^n$ and $\mathcal{T}_h^{n-1}$. A concrete example is given in Figure 1 where two meshes from consecutive time-steps are illustrated.

<figure><img src="image_2.png" /><figcaption>Figure 1: Illustration of two meshes from two consecutive time-steps.</figcaption></figure>

For any Banach space $\mathcal{X}$, we define the local and global tensor-product spaces

$$\mathbb{P}_{q_n}(I_n; \mathcal{X}) := \{w \in L^2(I_n; \mathcal{X}) \mid w = \nu b(t), \nu \in \mathcal{X}, b \in \mathbb{P}_{q_n}(I_n)\},$$

$$\mathbb{P}_q(J; \mathcal{X}) := \{w \in L^2(J; \mathcal{X}) \mid w|_{I_n} \in \mathbb{P}_q(I_n; \mathcal{X}) \forall n \in \mathcal{N}\}.$$

For any polynomial degree $q_n$, given the $L^2$-projector $\Pi_{q_n}: L^2(I_n; \mathcal{X}) \to \mathbb{P}_{q_n}(I_n; \mathcal{X})$ in the time interval $I_n$, we define the $L^2$-projector in time $\Pi_q: L^2(J; \mathcal{X}) \to \mathbb{P}_q(J; \mathcal{X})$ as $\Pi_q|_{I_n} := \Pi_{q_n}$, its orthogonal complement projector is given by $\Pi_q^\perp|_{I_n} := \Pi_q^\perp := I_d - \Pi_{q_n}$. Given the local space-time projector $\Pi_{p,q_n} := \Pi_p^n \circ \Pi_{q_n} = \Pi_{q_n} \circ \Pi_p^n: L^2(I_n; L^2(D)) \to \mathbb{P}_{q_n}(I_n; V_h^n)$ at each time-step, we define the global space-time projector $\Pi_{p,q}: L^2(J; L^2(D)) \to \mathbb{P}_q(J; H_0^1(D))$ as $\Pi_{p,q}|_{I_n} := \Pi_{p,q_n}$.

Henceforth, we assume that $p \ge 1$ and $q_n \ge 2$ for all $n \in \mathcal{N}$. Let $u_{0,h}$ and $u_{1,h}$ be approximations of the initial data $u_0$ and $u_1$ onto $V_h^0$ with optimal approximation properties in the $L^2$-norm; for instance, in the
numerical experiments of Section 6 we pick the $L^2$-projection onto the finite element space over the initial mesh and (if both initial data $u_0$ and $u_1$ in (2.1) are defined pointwise) the Lagrangian interpolant.

We define the space-time finite element space

$$V_{p,q} := \{W \in \mathbb{P}_q(J; H_0^1(D)) \mid W|_{I_n} \in \mathbb{P}_{q_n}(I_n; V_h^n) \quad \forall n \in \mathcal{N}\}.$$

For $w(t_0^-) \in H_0^1(D)$, we define a time jump operator $[w]_n := w(t_n^+) - w(t_n^-)$ for $n = 0, 1, \dots, N-1$; for $n = N$, we set $[w]_N := 0$. Note that $[w]_n = 0$ if $w$ continuous across $t_n$. For $w = U$, we define $U(t_0^-) := u_{0,h}$ and $U'(t_0^-) := u_{1,h}$.

The space-time method for problem (2.2) we shall be concerned with reads: find $U \in V_{p,q}$, such that

$$\left\{ \begin{array}{ll} \displaystyle \int_{I_n} \{(U'', W) + (\nabla U, \nabla W)\}dt + ([U']_{n-1}, W(t_{n-1}^+)) = \int_{I_n} (f, W)dt & \text{for all } W \in \mathbb{P}_{q_{n-1}}(I_n; V_h^n) \\ U(t_{n-1}^+) := \mathcal{P}_p^n U(t_{n-1}^-) & \forall n \in \mathcal{N} \end{array} \right. \quad (3.1)$$

with $\mathcal{P}_p^n: V_h^{n-1} \to V_h^n$ denoting the orthogonal $H^1$-projection operator. In the case of static spatial meshes, the method (3.1) coincides to that in [54], resulting in a continuous approximation in time.

# 4 Reconstructions and Baker's test function

We introduce some reconstructions that will be used below in the proof of the a posteriori error bound. Also, we show some approximation results for a non-standard test function to be used by, whose use in error analysis of finite element methods for the wave equation originates in [2].

## 4.1 A time reconstruction

Given $V \in \mathbb{P}_{q_n}(I_n; L^2(D))$, for all $n \in \mathcal{N}$ with given $V'(t_0^-) \in L^2(D)$, we define $\hat{V}$ with $\hat{V}|_{I_n} \in \mathbb{P}_{q_{n+1}}(I_n; L^2(D))$ by the relations

$$\left\{ \begin{array}{ll} \displaystyle \int_{I_n} (\hat{V}'', W) dt = \int_{I_n} (V'', W) dt + ([V']_{n-1}, W(t_{n-1}^+)) & \text{for all } W \in \mathbb{P}_{q_{n-1}}(I_n; L^2(D)) \\ \hat{V}(t_{n-1}^+) = V(t_{n-1}^+), & \hat{V}'(t_{n-1}^+) = V'(t_{n-1}^-). \end{array} \right. \quad (4.1)$$

The above definition is motivated by the respective construction due to Makridakis and Nochetto [44] for discontinuous Galerkin time-stepping methods for first order evolution problems; see also [47] for an extension to the hp-version and [15] for (4.1) in the case of static meshes and semi-discrete in time setting.

The time reconstruction above satisfies certain smoothness properties. In particular, it maps $C^0$ piecewise polynomials into $C^1$ piecewise polynomials.

**Lemma 4.1** (Smoothness of the time reconstruction). $\hat{V}$ satisfies

$$\hat{V}'(t_n^-) = V'(t_n^-) \in V_h^n, \qquad \hat{V}(t_n^-) = V(t_n^-) \in V_h^n \quad \forall n \in \mathcal{N}, \quad (4.2)$$

i.e., $\hat{V}' \in C^0(J; L^2(D))$. Further, if $V \in C^0(J; L^2(D))$, i.e.,

$$V(t_{n-1}^-) = V(t_{n-1}^+) \quad \forall n \in \mathcal{N} \cup \{0\}, \quad (4.3)$$

then $\hat{V} \in C^1(J; L^2(D))$ and it coincides with the reconstruction operator in [15, Sect. 3.2].

*Proof.* The proof is an extension of that of [15, Proposition 3.2], taking now into account mesh modification.

We begin by selecting $W = c_p \equiv c_p(x) \in L^2(D)$ a constant function with respect to $t$ in (4.1). Integrating by parts in time the first line in (4.1) on both sides, and imposing the initial condition on the time derivative in (4.1), we deduce $(\hat{V}'(t_n^-) - V'(t_n^-), c_p) = 0$. Choosing $c_p = \hat{V}'(t_n^-) - V'(t_n^-)$, we infer the first identity in (4.2).
As for the second identity in (4.2), integrating by parts in time twice both sides of the first line of (4.1) and invoking the first identity in (4.2), we find

$$\int_{I_n} (\hat{V} - V, W'')dt - (\hat{V}(t_n^-) - V(t_n^-), W'(t_n^-)) + (\hat{V}(t_{n-1}^+) - V(t_{n-1}^+), W'(t_{n-1}^+)) = 0.$$

We infer the second identity in (4.2) by picking $W = (\hat{V}(t_n^-) - V(t_n^-))(t - t_{n-1})$ and using $\hat{V}(t_{n-1}^+) = V(t_{n-1}^+)$. Finally, the continuity of $V$ in (4.3) implies that $\hat{V}$ is globally continuous. $\square$

We now estimate the error between $V$ and $\hat{V}$ by the norm of the jump of the time derivative of $V$.

**Lemma 4.2** (Orthogonality and bounds for the time reconstruction). For $n \in \mathcal{N}$, we have

$$\int_{I_n} (V - \hat{V}, W'')dt = 0 \quad \forall W \in \mathbb{P}_{q_{n-1}}(I_n; L^2(D)). \quad (4.4)$$

Moreover, we have

$$\|V' - \hat{V}'\|_{L^2(I_n; L^2)}^2 = \tau_n c_1(q_n)^2 \|[V']_{n-1}\|^2 \quad (4.5a)$$

$$\|V' - \hat{V}'\|_{L^{\infty}(I_n; L^2)}^2 = \|[V']_{n-1}\|^2 \quad (4.5b)$$

$$\|V - \hat{V}\|_{L^2(I_n; L^2)}^2 \le \tau_n^3 c_2(q_n)^2 \|[V']_{n-1}\|^2 \quad (4.5c)$$

$$\|V'' - \hat{V}'\|_{L^2(I_n; L^2)}^2 \le \tau_n^{-1} q_n^2 \|[V']_{n-1}\|^2, \quad (4.5d)$$

where

$$c_1(q) := \left( \frac{q}{(2q-1)(2q+1)} \right)^{\frac{1}{2}}, \quad c_2(q) := \left( \frac{q}{4(q-2)(q-1)(2q-1)(2q+1)} \right)^{\frac{1}{2}}, \quad (4.6)$$

for $q \ge 3$, while $c_1(2) = \sqrt{2/15}$ and $c_2(2) = \sqrt{1/15}$.

*Proof.* The proof of (4.5a) and (4.5b) is given in [15, Lemma 3.3]; the corresponding bounds can also be found in [35, 47]. To show (4.4), we integrate by parts twice in time the first line in (4.1), and invoke (4.2) and the second line in (4.1), to find

$$\int_{I_n} (V - \hat{V}, W'') dt = (V(t_n^-) - \hat{V}(t_n^-), W'(t_n^-)) = 0 \quad \forall W \in \mathbb{P}_{q_{n-1}}(I_n; L^2(D)).$$

We split the proof of (4.5c) into two cases: $q_n \ge 3$ and $q_n = 2$. For the case $q_n \ge 3$, we pick

$$W(t) := \int_{t_n}^{t} \int_{t_n}^{s} \Pi_{q-3}(V - \hat{V})(r) dr ds,$$

in (4.4), to deduce $\int_{I_n} (V - \hat{V}, \Pi_{q-3}(V - \hat{V}))dt = 0$. Now, [48, Theorem 3.11] implies

$$\|(I_d - \Pi_{q-3})(V - \hat{V})\|_{L^2(I_n; L^2)} \le \frac{\tau_n}{2\sqrt{(q_n - 2)(q_n - 1)}} \|V' - \hat{V}'\|_{L^2(I_n; L^2)}.$$

Combining the three displays above gives

$$\|V - \hat{V}\|_{L^2(I_n; L^2)}^2 = \int_{I_n} (V - \hat{V}, (I_d - \Pi_{q_n-3})(V - \hat{V}))dt \\ \le \frac{\tau_n}{2\sqrt{(q_n - 2)(q_n - 1)}} \|V - \hat{V}\|_{L^2(I_n; L^2)} \|(V - \hat{V})'\|_{L^2(I_n; L^2)},$$

which, together with (4.5a) and (4.5b), implies the third bound in (4.5) for the case $q_n \ge 3$.
For $q_n = 2$, we have the Poincaré inequality

$$ \int_0^1 |w(t)|^2 dt = \int_0^1 \left( \int_0^t w'(s) ds \right)^2 dt \le \int_0^1 t dt \int_0^1 |w'(s)|^2 ds \le \frac{1}{2} \int_0^1 |w'(s)|^2 ds, $$

for $w \in H^1(0, 1)$, with $w(0) = 0$, which upon scaling implies $\|V - \hat{V}\|_{L^2(I_n; L^2)} \le 2^{-\frac{1}{2}} \tau_n \|(V - \hat{V})'\|_{L^2(I_n; L^2)}$.

In order to obtain (4.5d), we introduce the lifting operator $L_n : V_{p,q-1} \to \mathbb{P}_{q_{n-1}}(I_n; L^2(D))$ as in [47, Section 4.1] for all $n \in \mathcal{N}$:

$$ \int_{I_n} (L_n(Z), W) dt = ([Z]_{n-1}, W) \quad \forall Z \in V_{p,q-1}, w \in \mathbb{P}_{q_{n-1}}(I_n; L^2(D)). $$

It readily follows from [47, Proposition 2] that $\|L_n(V')\|_{L^2(I_n; L^2)}^2 = q_n^2 \tau_n^{-1} \|[V']_{n-1}\|^2$. Invoking the above lifting operator and testing with $V'' - \hat{V}''$ in (4.1), we find

$$ \begin{aligned} \|(V - \hat{V})''\|_{L^2(I_n; L^2)}^2 &= \int_{I_n} (L_n(V'), (V - \hat{V})'') dt \le \|L_n(V')\|_{L^2(I_n; L^2)} \|(V - \hat{V})''\|_{L^2(I_n; L^2)} \\ &= q_n \tau_n^{-\frac{1}{2}} \|[V']_{n-1}\| \|(V - \hat{V})''\|_{L^2(I_n; L^2)}. \end{aligned} $$

□

We next show computable bounds of the spatial approximation error of the time reconstruction.

**Lemma 4.3 (Spatial error of the time reconstruction).** Let $V \in V_{p,q}$ with $V(t_0^-) \in V_h^0$, and consider any non-negative, piecewise constant on each space-time element $I_n \times K$, $n \in \mathcal{N}$, $K \in \mathcal{T}_h^n$, function $\rho_0$. Then, for $\ell \in \{2, 3\}$, we have

$$ \|\rho_0(\hat{V} - \Pi_p^n \hat{V})^{(\ell)}\|_{L^\infty(I_n; L^2)} \le c_3(q_n; \ell) \tau_n^{1-\ell} \|\rho_0 \Pi_p^{n,\perp} V'(t_{n-1}^-)\|, \quad (4.7) $$

$$ \text{where } c_3(q; \ell) := 4\sqrt{2}q(q+1) (2\sqrt{3}(q+1)^2)^{\ell-2}. $$

*Proof.* Since $V = \Pi_p V$, we find $\|\rho_0(\hat{V} - \Pi_p \hat{V})^{(\ell)}\|_{L^\infty(I_n; L^2)} = \|\rho_0 \Pi_p^{n,\perp}(\hat{V} - V)^{(\ell)}\|_{L^\infty(I_n; L^2)}$. We consider the cases $\ell = 2$ and $\ell = 3$ separately. For $\ell = 2$, since $\rho_0$ in constant on $I_n \times K$, the $L^\infty$-to-$L^2$ polynomial inverse inequality, see, e.g., [48, Theorem 3.92], implies

$$ \|\rho_0(\hat{V} - \Pi_p \hat{V})''\|_{L^\infty(I_n; L^2)} \le 4\sqrt{2}(q_n + 1)\tau_n^{-\frac{1}{2}} \|\rho_0(\hat{V}(t) - \Pi_p \hat{V}(t))''\|_{L^2(I_n; L^2)}. $$

Then we proceed as in (4.5d) to get $\|\rho_0(\hat{V} - \Pi_p \hat{V})''\|_{L^2(I_n; L^2)} \le q_n \tau_n^{-\frac{1}{2}} \|\rho_0 \Pi_p^{n,\perp}[V']_{n-1}\|$. The case $\ell = 2$ is concluded by observing that $V'(t_{n-1}^+) = \Pi_p^n V'(t_{n-1}^+)$.

For $\ell = 3$, we use the $H^1$-to-$L^2$ polynomial inverse inequality, see, e.g., [48, Theorem 3.91] (which contributes $2\sqrt{3}(q_n + 1)^2 \tau_n^{-1}$), to find

$$ \|\rho_0(\hat{V} - \Pi_p \hat{V})^{(3)}\|_{L^\infty(I_n; L^2)} \le 8\sqrt{6}(q_n + 1)^3 \tau_n^{-1} \|\rho_0(\hat{V}(t) - \Pi_p \hat{V}(t))''\|_{L^2(I_n; L^2)}. $$

Then we proceed as in the case $\ell = 2$ to conclude the proof.

□

## 4.2 A Hermite-type time reconstruction

**Definition 4.4 (A Hermite-type polynomial basis).** For all $n \in \mathcal{N}$, we consider the cubic Hermite polynomials $\{\phi_i^n\}_{i=1}^4$, defined as

$$ \left\{ \begin{array}{ll} \phi_1^n(t_{n-1}) = 1, & \phi_1^n(t_n) = (\phi_1^n)'(t_{n-1}) = (\phi_1^n)'(t_n) = 0 \\ \phi_2^n(t_n) = 1, & \phi_2^n(t_{n-1}) = (\phi_2^n)'(t_{n-1}) = (\phi_2^n)'(t_n) = 0 \\ (\phi_3^n)'(t_{n-1}) = 1, & \phi_3^n(t_n) = \phi_3^n(t_{n-1}) = (\phi_3^n)'(t_n) = 0 \\ (\phi_4^n)'(t_n) = 1, & \phi_4^n(t_{n-1}) = (\phi_4^n)'(t_{n-1}) = \phi_4^n(t_n) = 0. \end{array} \right. $$
Picking any basis {$\hat{\phi}_i^n$}$_{i=1}^{q_n-2}$ spanning $\mathbb{P}_{q_n-2}(I_n; \mathbb{R})$, we define the $\mathbb{C}^1$-bubble functions

$$ \phi_i^n := \left( (t - t_{n-1})(t_n - t) \right)^2 \hat{\phi}_{i-4}^n \quad \forall i = 5, \dots, q_n + 2. $$

The set {$\hat{\phi}_i^n$}$_{i=1}^{q_n+2}$ forms a basis of $\mathbb{P}_{q_n+1}(I_n; \mathbb{R})$. In the remainder of this section, $V \in \mathbb{P}_{q+1}(J; L^2(D))$.

**Definition 4.5** (Hermite-type time reconstruction operator). Noting the expansion $V|_{I_n} = V(t_{n-1}^+)\phi_1^n + V(t_n^-)\phi_2^n + \sum_{i=3}^{q_n+2} V_i\phi_i^n$, with $V(t_N^+) := V(t_N^-)$, we define the Hermite-type time reconstruction by

$$ \mathcal{I}_{\mathcal{C}}(V)|_{I_n} := \frac{V(t_{n-1}^-) + V(t_{n-1}^+)}{2} \phi_1^n + \frac{V(t_n^-) + V(t_n^+)}{2} \phi_2^n + \sum_{i=3}^{q_n+2} V_i \phi_i^n \quad \forall n \in \mathcal{N} \setminus \{1\}, \quad (4.8a) $$

$$ \mathcal{I}_{\mathcal{C}}(V)|_{I_1} := V(t_0^+)\phi_1^1 + \frac{V(t_1^-) + V(t_1^+)}{2}\phi_2^1 + \sum_{i=3}^{q_1+2} V_i\phi_i^1. \quad (4.8b) $$

In other words, the reconstruction operator modifies the coefficient of the first (resp. second) Hermite basis function from its right (resp. left) limit to its average. The reconstruction operator in (4.8a) satisfies the following conditions at the endpoints of each time interval: for all $n \in \mathcal{N}$,

$$ \mathcal{I}_{\mathcal{C}}(V)(t_{n-1}^+) = \mathcal{I}_{\mathcal{C}}(V)(t_{n-1}^-), \quad \mathcal{I}_{\mathcal{C}}(V)'(t_{n-1}^+) = V'(t_{n-1}^+), \quad \mathcal{I}_{\mathcal{C}}(V)'(t_{n-1}^-) = V'(t_{n-1}^-), \quad (4.9) $$

i.e., it maps piecewise polynomials into $\mathbb{C}^0$ piecewise polynomials as discussed in the next result.

**Lemma 4.6** (Smoothness of $\mathcal{I}_{\mathcal{C}}$). For any $V \in \mathbb{P}_{q+1}(J; L^2(D))$, $\mathcal{I}_{\mathcal{C}}(V) \in \mathbb{C}^0(J; L^2(D))$. If $V$ is also continuous at all time nodes, then $\mathcal{I}_{\mathcal{C}}(V) = V$. In addition, if $V' \in \mathbb{C}^0(J; L^2(D))$, then $\mathcal{I}_{\mathcal{C}}(V) \in \mathbb{C}^1(J; L^2(D))$. $\square$

Next, we estimate the error between $V$ and $\mathcal{I}_{\mathcal{C}}(V)$ by the jump of $V$ at the time nodes in several norms.

**Lemma 4.7** (Computable bounds for $V - \mathcal{I}_{\mathcal{C}}(V)$). Let $V \in \mathbb{P}_{q+1}(J; L^2(D))$. Then, for all $n \in \mathcal{N}$, we have

$$ \|V - \mathcal{I}_{\mathcal{C}}(V)\|_{L^1(I_n; L^2)} \le \frac{\tau_n}{4} (\|[V]_{n-1}\| + \|[V]_n\|), \quad (4.10a) $$

$$ \|V' - \mathcal{I}_{\mathcal{C}}(V)'\|_{L^1(I_n; L^2)} \le \frac{1}{2} (\|[V]_{n-1}\| + \|[V]_n\|). \quad (4.10b) $$

$$ \|V - \mathcal{I}_{\mathcal{C}}(V)\|_{L^{\infty}(I_n; L^2)} \le \frac{1}{2} (\|[V]_{n-1}\| + \|[V]_n\|), \quad (4.10c) $$

In addition, for all $w \in L^\infty((0, \xi]; L^2(D))$, $w|_{I_n} \in W^{1,\infty}(I_n; L^2(D))$ with $w(\xi) = 0$ for a given $0 < \xi \le T$, let $m \in \mathcal{N}$ the index such that $\xi \in I_m$, and set $\tilde{\tau}_n := 2\tau_n + \tau_{n-1} + (\xi - t_n)|\tau_{n-1}/\tau_n - 1|$. We have

$$ \left| \int_0^\xi (w, (I_d - \mathcal{I}_c)(V)) dt \right| \le \|w'\|_{L^\infty((0, \xi]; L^2)} \left[ \sum_{n=1}^m \left( \frac{\tau_n \tilde{\tau}_n}{4} \|[V]_{n-1}\| \right) + \frac{\tau_m^2}{4} \|[V]_m\| \right]. \quad (4.11) $$

*Proof.* For all $n \in \mathcal{N}$, direct computations given $\phi_1^n = 2\left(\frac{t-t_{n-1}}{\tau_n}\right)^3 - 3\left(\frac{t-t_{n-1}}{\tau_n}\right)^2 + 1$ and $\phi_2^n = -2\left(\frac{t-t_{n-1}}{\tau_n}\right)^3 + 3\left(\frac{t-t_{n-1}}{\tau_n}\right)^2$. We infer that

$$ \|\phi_1^n\|_{L^1(I_n)} = \|\phi_2^n\|_{L^1(I_n)} = \frac{\tau_n}{2}, \quad \|(\phi_1^n)'\|_{L^1(I_n)} = \|(\phi_2^n)'\|_{L^1(I_n)} = 1, \quad \|\phi_1^n\|_{L^\infty(I_n)} = \|\phi_2^n\|_{L^\infty(I_n)} = 1. $$

Definition (4.8a) further implies that in $I_n$,

$$ V - \mathcal{I}_{\mathcal{C}}(V) = \frac{[V]_{n-1}}{2} \phi_1^n - \frac{[V]_n}{2} \phi_2^n. $$
We combine the two displays above for all $n \in \mathcal{N}$ to get (4.10a), (4.10b), and (4.10c).

To show (4.11), we introduce $\Phi_1^n := \frac{(t-t_{n-1})^4}{2\tau_n^3} - \frac{(t-t_{n-1})^3}{\tau_n^2} + t-t_{n-1} + b_1^n$ and $\Phi_2^n := \frac{-(t-t_{n-1})^4}{2\tau_n^3} + \frac{(t-t_{n-1})^3}{\tau_n^2} + b_2^n$ as the anti-derivatives of $\phi_1^n$ and $\phi_2^n$, respectively, with $b_1^n, b_2^n \in \mathbb{R}$ to be determined later. Let $a_n := \frac{[V]_n}{2}$. Using that $[V]_0 = 0$, $w(\xi) = 0$, $\Phi_1^n(t_n^-) = \frac{\tau_n}{2} + b_1^n$, $\Phi_1^n(t_{n-1}^+) = b_1^n$, $\Phi_2^n(t_n^-) = \frac{\tau_n}{2} + b_2^n$ and $\Phi_2^n(t_{n-1}^+) = b_2^n$, we have

$$ \begin{aligned} \int_0^\xi (w, (I_d - \mathcal{I}_c)(V))dt &= \sum_{n=1}^{m-1} \left(w, a_{n-1}\Phi_1^n - a_n\Phi_2^n\right)|_{t_{n-1}^+}^{t_n^-} + \left(w, a_{m-1}\Phi_1^m - a_m\Phi_2^m\right)|_{t_{m-1}^+}^{\xi} \\ &\quad - \sum_{n=1}^{m-1} \int_{I_n} \left(w', a_{n-1}\Phi_1^n - a_n\Phi_2^n\right)dt - \int_{t_{m-1}}^\xi \left(w', a_{m-1}\Phi_1^m - a_m\Phi_2^m\right)dt \\ &\le \|w'\|_{L^\infty((0,\xi);L^2)} \sum_{n=1}^m \left( \|a_{n-1}\Phi_1^n\|_{L^1(I_n;L^2)} + \|a_n\Phi_2^n\|_{L^1(I_n;L^2)} \right) \\ &\quad + \sum_{n=1}^{m-1} \left( w(t_n), a_{n-1}(\frac{\tau_n}{2} + b_1^n) - a_n(\frac{\tau_n}{2} + b_2^n + b_1^{n+1}) + a_{n+1}b_2^{n+1} \right). \end{aligned} $$

We fix $b_1^n := \frac{-\tau_{n-1}}{2}$ and $b_2^n := 0$. This leads to

$$ \begin{aligned} \|a_{n-1}\Phi_1^n\|_{L^1(I_n;L^2)} &\le \frac{\tau_n(\tau_n + \tau_{n-1})}{4}\|[V]_{n-1}\|, \quad \|a_n\Phi_2^n\|_{L^1(I_n;L^2)} \le \frac{\tau_n^2}{4}\|[V]_n\|, \\ \left(w(t_n), a_{n-1}(\frac{\tau_n}{2} + b_1^n)\right) &\le (\xi - t_n)\|w'\|_{L^\infty((0,\xi);L^2)} \frac{\tau_n}{4} \left|\frac{\tau_{n-1}}{\tau_n} - 1\right| \|[V]_{n-1}\|. \end{aligned} $$

Inserting the above bounds in $\int_0^\xi (w, (I_d - \mathcal{I}_c)(V))dt$, we complete the proof. $\square$

**Remark 4.8.** Bounds (4.10b) and (4.10c) are optimal in $q_n$, while bound (4.10a) is $q_n$-suboptimal by one order.

**Remark 4.9** (Restriction on the adaptive time-steps). All terms in the telescope sum in (4.11) remain $O(\tau_n^2)$ if $|\tau_{n-1} - \tau_n| = O(\tau_n^2/(\xi - t_n))$ and/or if $\|[V]_{n-1}\|$ is small enough.

## 4.3 An elliptic reconstruction

Let $\Delta_{pw}$ denote the broken Laplace operator ('pw' stands for 'piecewise'), defined as $\Delta_{pw}V|_{I_n \times K} := \Delta V$ for all $n \in \mathcal{N} \cup \{0\}$ with $I_0 = \{0\}$, $K \in \mathcal{T}_h^n$, and $V \in V_h^n$. Consider also the discrete Laplacian $\Delta_h: H^1(D) \to H^1(D)$, whose restriction $\Delta_h^n := \Delta_h|_{I_n}: V_h^n \to V_h^n$ on each time interval $I_n$, $n \in \mathcal{N}$, is defined as follows: for $W \in V_h^n$,

$$ \langle -\Delta_h^n W, Z \rangle = (-\Delta_h^n W, Z) := (\nabla W, \nabla Z), \quad \text{for all } Z \in V_h^n. \quad (4.12) $$

**Definition 4.10** (Elliptic reconstruction). For $U \in V_{p,q}$ (below $U$ will be the solution to (3.1)) and for $t \in I_n$, $n \in \mathcal{N}$, we define its elliptic reconstruction $\tilde{U}$ with $\tilde{U}|_{I_n} \in \mathbb{P}_{q_n}(I_n; H_0^1(D) \cap H^2(D))$ of $U$, by

$$ (\nabla \tilde{U}, \nabla w) := (g, w) \quad \text{for all } w \in H_0^1(D), \quad \text{where } g = g(U) := -\Delta_h U - (\Pi_{p,q-1} - \Pi_{q-1})f + \Pi_p^\perp \tilde{U}'' . \quad (4.13) $$

At $t = t_0^-$, upon setting $g(t_0^-) := -\Delta_h^0 u_{0,h}$ and $g'(t_0^-) := -\Delta_h^0 u_{1,h}$, we define $\tilde{U}(t_0^-) \in H_0^1(D)$ and $\tilde{U}'(t_0^-) \in H_0^1(D)$ to be the solutions to the elliptic equations $-\Delta \tilde{U}(t_0^-) = g(t_0^-)$ and $-\Delta \tilde{U}'(t_0^-) = g'(t_0^-)$ equipped with homogeneous Dirichlet boundary conditions.

The function $g$ in (4.13) is a piecewise polynomial in time on each time interval taking values in $L^2(D)$. An integration by parts in space on the left-hand side of (4.13) implies that $\Delta \tilde{U}$ is also a polynomial in time; moreover we have that

$$ -\Delta \tilde{U} = g \in \mathbb{P}_q(J; L^2(D)). \quad (4.14) $$
For a given positive constant $c_{L^2}$ independent of the spatial mesh-size and $p$, an arbitrary finite element space $V_h$ with associated cell set $\mathcal{T}_h$ and interior facet set $\mathcal{F}_I$, $w$ in $H_0^1(D)$ with $\Delta w|_K$ in $L^2(K)$ for all $K$ in $\mathcal{T}_h$, $s$ in $L^2(D)$, we introduce the spatial residual-type error estimator

$$ \mathcal{E}(V_h, w, s) := c_{L^2} \left( \left\| \frac{h^2}{p^2} (s + \Delta_{pw} w) \right\|^2 + \mathcal{J}^2(\mathcal{F}_I, w) \right)^{\frac{1}{2}}, \quad \mathcal{J}(\mathcal{F}_I, w)^2 := \sum_{F \in \mathcal{F}_I} \frac{h_F^3}{p^3} \left\| \mathbf{n}_F \cdot \nabla w \right\|_{L^2(F)}^2. \quad (4.15) $$

The error between $U$ to (3.1) and its elliptic reconstruction $\tilde{U}$ in (4.13) can be measured in terms of the spatial error estimator above; of course, other choices of error estimators are also possible. For the sake of brevity, we introduce the notation $\Delta_h^- := \Delta_{pw} - \Delta_h$.

**Lemma 4.11** (Computable bounds for the elliptic reconstruction operator). Let $U \in V_{p,q}$, and $\tilde{U}$ and $g$ be as in (4.13). Then, for all $n \in \mathcal{N}$, $t \in I_n$, $l = 0, 1$, we have the error bounds

$$ \begin{aligned} \|(U - \tilde{U})^{(l)}(t)\| &\le c_{L^2} \left( \left\| \frac{h_n^2}{p^2} \Delta_h^- U^{(l)}(t) \right\| + \mathcal{J}(\mathcal{F}_I^n, U^{(l)}(t)) \right. \\ &\qquad \left. + \left\| \frac{h_n^2}{p^2} (\Pi_{q-1} \Pi_p^{n,\perp} f)^{(l)}(t) \right\| + c_3(q_n; l+2) \tau_n^{-l-1} \left\| \frac{h_n^2}{p^2} \Pi_p^{n,\perp} U'(t_{n-1}^-) \right\| \right). \end{aligned} \quad (4.16) $$

For $n = 0$ and $l = 0, 1$, we have

$$ \|(U - \tilde{U})^{(l)}(t_0^-)\| \le c_{L^2} \left( \left\| \frac{h_0^2}{p^2} \Delta_h^- U^{(l)}(t_0^-) \right\| + \mathcal{J}(\mathcal{F}_I^1, U^{(l)}(t_0^-)) \right). \quad (4.17) $$

In addition, the following bounds also hold true for all $n \in \mathcal{N} \setminus \{1\}$ and $l = 0, 1$:

$$ \begin{aligned} \|[(U - \tilde{U})^{(l)}]_{n-1}\| &\le c_{L^2} \left( \left\| \frac{(h_n^{\ominus})^2}{p^2} [\Delta_h^- U^{(l)}]_{n-1} \right\| + \mathcal{J}(\mathcal{F}_I^{n,\ominus}, [U^{(l)}]_{n-1}) \right. \\ &\qquad \left. + \left\| \frac{(h_n^{\ominus})^2}{p^2} [(\Pi_{q-1} \Pi_p^{\perp} f)^{(l)}]_{n-1} \right\| + \sum_{k=n-1}^n \frac{c_3(q_k; 2)}{\tau_k^{l+1}} \left\| \frac{(h_n^{\ominus})^2}{p^2} \Pi_p^{k,\perp} U'(t_{k-1}^-) \right\| \right). \end{aligned} \quad (4.18) $$

For $n = 1$, we have

$$ \| [U - \tilde{U}]_0 \| \le c_{L^2} \left\| \frac{(h_1^{\ominus})}{p^2} [(\Pi_{q-1} \Pi_p^{\perp} f)(t_0^{+})] \right\| \quad (4.19a) $$

$$ \|[(U - \tilde{U})']_0\| \le c_{L^2} \left( \left\| \frac{(h_1^{\ominus})^2}{p^2} [\Delta_h^- U']_0 \right\| + \mathcal{J}(\mathcal{F}_I^{1,\ominus}, [U']_0) + \left\| \frac{(h_1^{\ominus})}{p^2} [(\Pi_{q-1} \Pi_p^{\perp} f)'(t_0^{+})] \right\| \right). \quad (4.19b) $$

*Proof.* By construction, (4.13) implies the Galerkin orthogonality $(\nabla(U - \tilde{U}), \nabla W) = 0$ for all $W \in V_h^n$, standard arguments, see, e.g., [53], reveal that

$$ \|(U - \tilde{U})(t)\| \le \mathcal{E}(V_h^n, U(t), g(t)), \quad (4.20) $$

for specific choices of $c_{L^2}$ in (4.15). Now, using the definition of $g$ in (4.13) and $\mathcal{E}(\cdot, \cdot, \cdot)$ in (4.15), we arrive at

$$ \|(U - \tilde{U})(t)\| \le c_{L^2} \left( \left\| \frac{h_n^2}{p^2} \Delta_h^- U(t) \right\| + \mathcal{J}(\mathcal{F}_I^n, U(t)) + \left\| \frac{h_n^2}{p^2} (\Pi_{q-1} \Pi_p^{n,\perp} f)(t) \right\| + \left\| \frac{h_n^2}{p^2} \Pi_p^{n,\perp} \hat{U}''(t) \right\| \right). $$

Applying Lemma 4.3 with $\rho_0 = \frac{h_n^2}{p^2}$, we find $\|\frac{h_n^2}{p^2} \Pi_p^{n,\perp} \hat{U}''(t)\| \le c_3(q_n; 2) \tau_n^{-1} \|\frac{h_n^2}{p^2} \Pi_p^{n,\perp} U'(t_{n-1}^-)\|$, which leads to (4.16) with $l = 0$. Inequality (4.16) with $l = 1$ and (4.17) are proved analogously.
Next, we show (4.18). Notice that $[\tilde{U}]_{n-1}$ satisfies $-\Delta[\tilde{U}]_{n-1} = [g]_{n-1}$ in $L^2(D)$, with homogeneous Dirichlet boundary condition and, thus, we have the Galerkin orthogonality $(\nabla([U - \tilde{U}]_{n-1}), \nabla W) = 0$ for all $W \in V_h^{n, \ominus}$, implying, as before,

$$ \| [U - \tilde{U}]_{n-1} \| \le \mathcal{E}(V_h^{n, \ominus}, [U]_{n-1}, [g]_{n-1}). \quad (4.21) $$

Then, (4.15) and (4.7) lead to (4.18) with $l = 0$. The bound, for $l = 1$, as well as (4.19a) and (4.19b) are proven completely analogously. □

**Remark 4.12** (Sharper bound on the elliptic reconstruction operator). In the proof of Lemma 4.11, we used the inequality

$$ \| \frac{(h_n^{\ominus})^2}{p^2} [\Pi_p^{\perp} \hat{U}^{\prime\prime}]_{n-1} \| \le c_3(q_n; 2) \tau_n^{-1} \| \frac{(h_n^{\ominus})^2}{p^2} \Pi_p^{n, \perp} U'(t_{n-1}^{-}) \| + c_3(q_{n-1}; 2) \tau_{n-1}^{-1} \| \frac{(h_n^{\ominus})^2}{p^2} \Pi_p^{n-1, \perp} U'(t_{n-2}^{-}) \|, $$

to derive an easily computable bound without computing $\hat{U}$. Since $\hat{U}$ is a polynomial in time, $\hat{U}^{\prime\prime}$ can be computed practically; this would lead to a sharper bound in the analysis at the price of extra computations on the practical level.

## 4.4 A space-time reconstruction and Baker's test function

**Definition 4.13.** We define the space-time reconstruction operator

$$ \bar{U} := \mathcal{I}_c(\hat{\tilde{U}}). \quad (4.22) $$

Due to Lemmas 4.1 and 4.6, we have $\bar{U} \in \mathcal{C}^1(J; H_0^1(D) \cap H^2(D))$: the space-time reconstruction operator maps piecewise polynomials into globally $\mathcal{C}^1$ in time, piecewise polynomials.

Let $U$ be the solution to the space-time finite element method (3.1). We define $\xi \in J$ and $m \in \mathcal{N}$ as

$$ \|(u - \bar{U})(\xi)\|_{L^2} := \|u - \bar{U}\|_{L^\infty(L^2)} \quad \text{with } \xi \in I_m. \quad (4.23) $$

Inspired by the classical construction in [2], we set $v \in \mathcal{C}^0(J; H_0^1(D)) \cap \mathcal{C}^1(J; L^2(D))$, with

$$ v(t) := \int_t^\xi (u - \bar{U})(s) \, ds. \quad (4.24) $$

Next, we recall the following technical result.

**Lemma 4.14** ([15, Lemma 3.4], Approximation and Poincaré inequality for $v$). For $m$ and $v$ as in (4.23) and (4.24), we have

$$ \|v - \Pi_q^\infty v\|_{L^\infty(I_n; L^2)} \le \frac{\tau_n}{\pi} c_4(q) \|v'\|_{L^\infty(I_n; L^2)} \quad n = 1, \dots, m-1, \quad q \in \mathbb{N}, \quad (4.25a) $$

$$ \|v\|_{L^\infty((t_{m-1}, \xi]; L^2)} \le \tau_m \|v'\|_{L^\infty((t_{m-1}, \xi]; L^2)}, \quad (4.25b) $$

with $\Pi_q^\infty$ being the $L^\infty$-projector [51, Theorem 7.2] in time of order $q$, and $c_4(q) := \sqrt{\pi}$ if $q < 3$, while $c_4(q) := (q-2)^{-1}$ if $q \ge 3$. □

## 5 A posteriori error upper bounds

We henceforth assume that

$$ \Delta u \in L^1(J; L^2(D)). \quad (5.1) $$

This implies that $v \in \mathcal{C}^0(J; H_0^1(D) \cap H^2(D)) \cap \mathcal{C}^1(J; L^2(D))$ owing to the domain convexity assumption, with $v$ defined in (4.24). Sufficient conditions for the validity of (5.1) are $u_0 \in H^2(D) \cap H_0^1(D)$, $u_1 \in H^1(D)$, and
$f \in H^1(J; L^2(D))$, following the argument in [21, Theorem 5] combined with the elliptic regularity result $H^2(D)$ on convex domains in [32].

We introduce

$$e := u - U = (u - \bar{U}) + (\bar{U} - U) =: \rho_1 + \rho_2. \qquad (5.2)$$

To prove an upper bound of $\|u - U\|_{L^\infty(L^2)}$, we estimate from above the corresponding norms of $\rho_1$ and $\rho_2$ separately.

We will frequently use the Kronecker delta function defined as $\delta_{i,j} = 1$ if $i = j$, and $\delta_{i,j} = 0$ if $i \neq j$. We also recall that $\Pi_p^\perp$ and $\Pi_q^\perp$ are the orthogonal complement projectors in space and time, respectively, and $\tilde{\tau}_n = 2\tau_n + \tau_{n-1} + (\xi - t_n)|\tau_{n-1}/\tau_n - 1|$. Also, given $m$ as in (4.23), we define

$$c_5(n) := \begin{cases} 1 & \text{if } n = m \\ \frac{t_m - t_{n-1}}{\tau_n} & \text{if } q_n = 2, n < m; \\ \frac{c_4(q_n - 3)}{\pi} & \text{if } q_n > 2, n < m, \end{cases} \quad \text{and} \quad \tilde{c}_5(n) := \begin{cases} \frac{c_4(q_n - 1)}{\pi} & \text{if } n < m; \\ 1 & \text{if } n = m. \end{cases} \qquad (5.3)$$

For convex $D$, the classical Poincaré-Steklov inequality states that, for $C_{PS} := \text{diam}(D)/\pi$, we have

$$\|w\| \le C_{PS} \|\nabla w\| \quad \text{for all } w \in H_0^1(D). \qquad (5.4)$$

We prove bounds for the $L^\infty(L^2)$ norms of $\rho_1$ and $\rho_2$ in (5.2). We begin with $\rho_1$.

**Lemma 5.1** (Upper bound on $\|\rho_1\|_{L^\infty(L^2)}$). Let $u$ and $U$ the solutions to (2.2) and (3.1) with initial conditions $u_0$ and $u_1$, and $u_{0,h}$ and $u_{1,h}$, respectively. For $c_1, c_2, c_5, \tilde{c}_5, m$, as in (4.6) (5.3), and (4.23), respectively, we have

$$\|\rho_1\|_{L^\infty(L^2)} \le \sqrt{2}\|u_0 - u_{0,h}\| + \sqrt{2}C_{PS}\|u_1 - u_{1,h}\| + \sqrt{2}C_{PS}\mathcal{E}(V_h^0, u_{1,h}, g'(t_0^-)) \\ + 2\mathcal{E}(V_h^1, U(t_0^+), g(t_0^+)) + \tau_m^2\|[g]_m\| + \mathcal{E}(V_h^{m+1, \ominus}, [U]_m, [g]_m) \\ + 2 \sum_{n=1}^m \left[ \begin{aligned} & \tau_n \tilde{c}_5(n) \|\Pi_{q-1}^\perp f\|_{L^1(I_n; L^2)} + \|[U]_{n-1}\| + \int_{I_n} \mathcal{E}(V_h^n, U'(t), g'(t)) dt \\ & + \tau_n c_1(q_n) \mathcal{E}(V_h^{n, \ominus}, [U']_{n-1}, [g']_{n-1}) + \mathcal{E}(V_h^{n, \ominus}, [U]_{n-1}, [g]_{n-1}) \\ & + \frac{\tau_n \tilde{\tau}_n}{4} \|[g]_{n-1}\| + \tau_n^3 c_2(q_n) c_5(n) \|[g']_{n-1}\| + \tau_n \tilde{c}_5(n) \|\Pi_{q-1}^\perp \Delta_h U\|_{L^1(I_n; L^2)} \end{aligned} \right].$$

*Proof.* For each $I_n, n \in \mathcal{N}$, using (4.1), we can rewrite (3.1) in pointwise strong form as

$$(\Pi_p \hat{U})'' - \Pi_{q-1} \Delta_h U = \Pi_{p,q-1} f.$$

Employing (4.13), (4.14), and (4.22), we find

$$\bar{U}'' - \Delta \bar{U} = \Pi_{q-1} f + (\bar{U} - \hat{U})'' + \Delta(\tilde{U} - \bar{U}) - \Pi_{q-1}^\perp \Delta_h U.$$

Subtracting the above identity from (2.1), we arrive at

$$\rho_1'' - \Delta \rho_1 = \Pi_{q-1}^\perp f + (\hat{U} - \bar{U})'' + \Delta(\bar{U} - \tilde{U}) + \Pi_{q-1}^\perp \Delta_h U.$$

Testing the last identity against $v$ given in (4.24), integrating in space and in the time interval $(0, \xi)$, with $\xi$ as in (4.23), and, finally, integrating by parts in space on the left-hand side, we deduce the identity

$$T_1 + T_2 = T_3 + T_4 + T_5 + T_6,$$

where

$$T_1 := \int_0^\xi (\rho_1'', v) \, dt, \quad T_2 := \int_0^\xi (\nabla \rho_1, \nabla v) \, dt, \quad T_3 := \int_0^\xi (\Pi_{q-1}^\perp f, v) \, dt,$$

$$T_4 := \int_0^\xi ((\hat{U} - \bar{U})'', v) \, dt, \quad T_5 := \int_0^\xi (\Delta(\bar{U} - \tilde{U}), v) \, dt, \quad T_6 := \int_0^\xi (\Pi_{q-1}^\perp \Delta_h U, v) \, dt.$$

We estimate $T_1 + T_2$ from below, and $T_3, \dots, T_6$ from above, and eventually collect the bounds.
**Lower bound on $T_1 + T_2$.** We will use the trivial identity $(f_1'', f_2) = -(f_1', f_2') + (f_1', f_2)$ for all sufficiently smooth functions $f_1$ and $f_2$. Since $v \in C^0(J; H_0^1(D)) \cap C^1(J; L^2(D))$ and $\bar{U} \in C^1(J; H_0^1(D))$, the above identity with $v(\xi) = 0$, and Young's inequality imply

$$ T_1 = \int_0^\xi (\rho_1', \rho_1) dt + \int_0^\xi (\rho_1', v)' dt = \frac{1}{2} \|\rho_1(\xi)\|^2 - \frac{1}{2} \|\rho_1(t_0^+)\|^2 - (\rho_1'(t_0^+), v(t_0^+)) \\ \geq \frac{(5.4)}{2} \|\rho_1(\xi)\|^2 - \frac{1}{2} \|\rho_1(t_0^+)\|^2 - C_{PS}^2 \|\rho_1'(t_0^+)\|^2 - \frac{1}{4} \|\nabla v(t_0^+)\|^2. $$

Definition (4.24) implies $v' = -\rho_1$ and $v(\xi) = 0$, whence we infer that

$$ T_2 = - \int_0^\xi (\nabla v', \nabla v) dt = \frac{1}{2} \|\nabla v(t_0^+)\|^2. $$

Since $\rho_1$ is continuously differentiable in $t$ over $J$, we observe that

$$ \|\rho_1(t_0^+)\| \stackrel{(4.9)}{=} \|u_0 - \tilde{U}(t_0^+)\| \stackrel{(4.1)}{=} \|u_0 - \tilde{U}(t_0^+)\| \le \|u_0 - u_{0,h}\| + \|[U]_0\| + \|(U - \tilde{U})(t_0^+)\|, \\ \|\rho_1'(t_0^+)\| \stackrel{(4.9)}{=} \|u_1 - \tilde{U}'(t_0^+)\| \stackrel{(4.1)}{=} \|u_1 - \tilde{U}'(t_0^-)\| \le \|u_1 - u_{1,h}\| + \|(U - \tilde{U})'(t_0^-)\|. $$

Combining the three displays above with $[U]_0 = 0$, (4.20) and (4.23) implies

$$ T_1 + T_2 \geq \frac{1}{2} \|\rho_1\|_{L^\infty(L^2)}^2 - \|u_0 - u_{0,h}\|^2 - \mathcal{E}(V_h^1, U(t_0^+), g(t_0^+))^2 - C_{PS}^2 (\|u_1 - u_{1,h}\|^2 + \mathcal{E}(V_h^0, u_{1,h}, g'(t_0^-))^2). $$

**Upper bound on $T_3$.** Set $w = \Pi_{q-1}^\infty v$ as in (4.25) to get

$$ T_3 = \sum_{n=1}^{m-1} \int_{I_n} (\Pi_{\mathbf{q}-1}^{\perp} f, v) dt + \int_{t_{m-1}}^{\xi} (\Pi_{\mathbf{q}-1}^{\perp} f, v) dt = \sum_{n=1}^{m-1} \int_{I_n} (\Pi_{\mathbf{q}-1}^{\perp} f, v - w) dt + \int_{t_{m-1}}^{\xi} (\Pi_{\mathbf{q}-1}^{\perp} f, v) dt. $$

For $1 \le n \le m - 1$, we have

$$ \int_{I_n} (\Pi_{\mathbf{q}-1}^{\perp} f, v - w) dt \stackrel{(4.25a)}{\le} \frac{\tau_n}{\pi} c_4(q_n - 1) \|\Pi_{\mathbf{q}-1}^{\perp} f\|_{L^1(I_n; L^2)} \|v'\|_{L^\infty(I_n; L^2)} \\ \le \frac{\tau_n}{\pi} c_4(q_n - 1) \|\Pi_{\mathbf{q}-1}^{\perp} f\|_{L^1(I_n; L^2)} \|v'\|_{L^\infty(L^2)}, $$

and, for $n = m$,

$$ \int_{t_{m-1}}^{\xi} (\Pi_{\mathbf{q}-1}^{\perp} f, v) dt \stackrel{(4.25b)}{\le} \|\Pi_{\mathbf{q}-1}^{\perp} f\|_{L^1(I_m; L^2)} \tau_m \|v'\|_{L^\infty(I_m; L^2)} \le \|\Pi_{\mathbf{q}-1}^{\perp} f\|_{L^1(I_m; L^2)} \tau_m \|v'\|_{L^\infty(L^2)}. $$

Combining, we find

$$ T_3 \le \sum_{n=1}^{m} \left( \tau_n \tilde{c}_5(n) \|\Pi_{\mathbf{q}-1}^{\perp} f\|_{L^1(I_n; L^2)} \right) \|\rho_1\|_{L^\infty(L^2)}. $$

**Upper bound on $T_4$.** Using that $v(\xi) \stackrel{(4.24)}{=} 0$, $\hat{U}'(t_0^+) \stackrel{(4.1)}{=} u_{1,h}$, $\Delta_h^0 u_{1,h} = \Delta \tilde{U}'(t_0^-) \stackrel{(4.22)}{=} \Delta \bar{U}'(t_0^+)$, and the fact that $\hat{U}'$ and $\bar{U}'$ belong to $C^0(J; L^2(D))$, we write

$$ T_4 = \int_0^\xi ((\hat{U} - \bar{U})', v)' dt + \int_0^\xi ((\hat{U} - \bar{U})', \rho_1) dt \\ = \int_0^\xi ((\hat{U} - \bar{U})', \rho_1) dt + ((u_{1,h} - \Delta^{-1}\Delta_h^0 u_{1,h}), v(t_0^-)) \\ \le \int_0^\xi ((\hat{U} - \bar{U})', \rho_1) dt + C_{PS}^2 \mathcal{E}(V_h^0, U'(t_0^-), g'(t_0^-))^2 + \frac{1}{4} \|\nabla v(t_0^+)\|^2. $$
The first term on the right-hand side can be rewritten as

$$ \int_0^\xi ((\hat{U} - \bar{U})', \rho_1) dt = \int_0^\xi ((\hat{U} - \tilde{\hat{U}})', \rho_1) dt + \int_0^\xi ((\tilde{\hat{U}} - \bar{U})', \rho_1) dt =: T_{4,1} + T_{4,2}. $$

As for $T_{4,1}$, we have

$$ T_{4,1} \le \|\rho_1\|_{L^\infty(L^2)} \sum_{n=1}^{m} \|(\hat{U} - \tilde{\hat{U}})\|_{L^1(I_n; L^2)}. $$

For all $n \in \{1, \dots, m\}$, we write

$$ \begin{aligned} \|(\hat{U} - \tilde{\hat{U}})'\|_{L^1(I_n; L^2)} & \stackrel{(4.5a)}{\le} \|(U - \tilde{U})'\|_{L^1(I_n; L^2)} + \tau_n c_1(q_n) \|[(U - \tilde{U})']_{n-1}\| \\ & \stackrel{(4.20), (4.21)}{\le} \int_{I_n} \mathcal{E}(V_h^n, U'(t), g'(t)) dt + \tau_n c_1(q_n) \mathcal{E}(V_h^{n, \ominus}, [U']_{n-1}, [g']_{n-1}). \end{aligned} $$

As for $T_{4,2}$, we have $T_{4,2} \le \|\rho_1\|_{L^\infty(L^2)} \sum_{n=1}^m \|(\hat{U} - \bar{U})'\|_{L^1(I_n; L^2)}$. We arrive at

$$ \begin{aligned} \sum_{n=1}^{m} \|(\hat{U} - \bar{U})'\|_{L^1(I_n; L^2)} & \stackrel{(4.10b)}{\le} \|\|\tilde{\hat{U}}\|_{n-1}\| + \frac{1}{2}\|\|\tilde{\hat{U}}\|_{m}\| \stackrel{(4.1),(4.2)}{=} \sum_{n=1}^{m} \|\|\tilde{U}\|_{n-1}\| + \frac{1}{2}\|\|\tilde{U}\|_{m}\| \\ & \le \sum_{n=1}^{m} (\|\|U\|_{n-1}\| + \|\|\tilde{U} - U\|_{n-1}\|) + \frac{1}{2}(\|\|U\|_{m}\| + \|\|\tilde{U} - U\|_{m}\|). \end{aligned} $$

Finally, combining the above displays, we get

$$ T_4 \le \|\rho_1\|_{L^\infty(L^2)} \sum_{n=1}^{m} \left[ \begin{aligned} & \int_{I_n} \mathcal{E}(V_h^n, U'(t), g'(t)) dt + \tau_n c_1(q_n) \mathcal{E}(V_h^{n, \ominus}, [U']_{n-1}, [g']_{n-1}) \\ & + \| [U]_{n-1} \| + \mathcal{E}(V_h^{n, \ominus}, [U]_{n-1}, [g]_{n-1}) + \frac{1}{2} \| [U]_m \| + \frac{1}{2} \mathcal{E}(V_h^{m+1, \ominus}, [U]_m, [g]_m) \end{aligned} \right]. $$

**Upper bound on $T_5$.** We split $T_5$ as

$$ T_5 = \int_0^\xi (\Delta(\bar{U} - \tilde{\hat{U}}), v) + \int_0^\xi (\Delta(\tilde{\hat{U}} - \bar{U}), v) =: T_{5,1} + T_{5,2}. $$

As for $T_{5,1}$, the crucial bound (4.11), (4.2), and (4.14) give

$$ T_{5,1} \le \|\rho_1\|_{L^\infty(L^2)} \sum_{n=1}^{m} \left( \frac{\tau_n \tilde{\tau}_n}{4} \| [g]_{n-1} \| \right) + \frac{\tau_m^2}{4} \| [g]_m \|. $$

For $T_{5,2}$, we distinguish two cases. For $q_n > 2$, using that $v|_{I_n} \in C^0(I_n; H^2(D) \cap H_0^1(D))$ for all $n \in 1, \dots, m-1$, invoking the orthogonality (4.4) with $w(t) = \int_{t_n}^t \int_{t_n}^s (\Delta \Pi_{q_n-3}^\infty v)(r) dr ds$ in $\mathbb{P}_{q_n-1}(I_n; L^2(D))$, we have

$$ \int_{I_n} (\Delta(\tilde{\hat{U}} - \tilde{\hat{U}}), \Pi_{q_n-3}^\infty v) dt = \int_{I_n} (\tilde{\hat{U}} - \hat{U}, \Delta \Pi_{q_n-3}^\infty v) dt = \int_{I_n} (\tilde{\hat{U}} - \hat{U}, w'') dt \stackrel{(4.4)}{=} 0. $$

Thus,

$$ T_{5,2} = \sum_{n=1}^{m-1} \int_{I_n} (\Delta(\tilde{\hat{U}} - \tilde{\hat{U}}), v - \Pi_{q_n-3}^\infty v) dt + \int_{t_{m-1}}^\xi (\Delta(\tilde{\hat{U}} - \tilde{\hat{U}}), v) dt. $$
The first terms on the right-hand side can be estimated from above as follows:

$$
\begin{align*}
\int_{I_n} (\Delta(\tilde{U} - \tilde{U}), v - \Pi_{q_n-3}^{\infty} v) dt &\stackrel{(4.25a)}{\le} \frac{\tau_n^3}{\pi} c_4(q_n-3) \|\Delta(\tilde{U} - \tilde{U})\|_{L^2(I_n; L^2)} \|\rho_1\|_{L^{\infty}(L^2)} \\
&\stackrel{(4.5b)}{=} \frac{\tau_n^3}{\pi} c_2(q_n) c_4(q_n-3) \|[\Delta\tilde{U}']_{n-1}\| \|\rho_1\|_{L^{\infty}(L^2)} \stackrel{(4.14)}{=} \frac{\tau_n^3}{\pi} c_2(q_n) c_4(q_n-3) \|[g']_{n-1}\| \|\rho_1\|_{L^{\infty}(L^2)}.
\end{align*}
$$

For $q_n = 2$, we have $\|v\|_{L^\infty(I_n;L^2)} \le (t_m - t_{n-1})\|\rho_1\|_{L^\infty(L^2)}$, whence we infer

$$
\begin{align*}
\int_{I_n} (\Delta(\tilde{U} - \tilde{U}), v) dt &\stackrel{(4.25a), (4.5)}{\le} (t_m - t_{n-1}) \tau_n^2 c_2(q_n) \|[\Delta \tilde{U}']_{n-1}\| \|\rho_1\|_{L^\infty(L^2)} \\
&\stackrel{(4.14)}{=} (t_m - t_{n-1}) \tau_n^2 c_2(q_n) \|[g']_{n-1}\| \|\rho_1\|_{L^\infty(L^2)}.
\end{align*}
$$

The second term on the right-hand side is controlled by

$$
\int_{t_{m-1}}^{\xi} (\Delta(\tilde{U} - \tilde{U}), v) dt \stackrel{(4.25b), (4.14), (4.5)}{\le} \tau_m^3 c_2(q_m) \|[g']_{m-1}\| \|\rho_1\|_{L^\infty(L^2)}.
$$

Collecting the above bounds and recalling the definition of $c_5(n)$ in (5.3), we find

$$
T_{5,2} \le \left( \sum_{n=1}^{m} \tau_n^3 c_2(q_n) c_5(n) \|[g']_{n-1} \| \right) \|\rho_1\|_{L^\infty(L^2)}.
$$

Therefore, the above bounds collect into

$$
T_5 \le \|\rho_1\|_{L^\infty(L^2)} \sum_{n=1}^{m} \left[ \frac{\tau_n \tilde{\tau}_n}{4} \|[g]_{n-1}\| + \tau_n^3 c_2(q_n) c_5(n) \|[g']_{n-1}\| \right] + \frac{\tau_m^2}{4} \|[g]_m\|.
$$

**Upper bound on $T_6$.** For $m$ as in (4.23) we set $w = \Pi_{q-1}^{\infty} v$ as in (4.25) and get

$$
\begin{align*}
T_6 &= \sum_{n=1}^{m-1} \int_{I_n} (\Pi_{\boldsymbol{q}-1}^{\perp} \Delta_h U, v-w) dt + \int_{t_{m-1}}^{\xi} (\Pi_{\boldsymbol{q}-1}^{\perp} \Delta_h U, v) dt \\
&\le \|\rho_1\|_{L^{\infty}(L^2)} \left( \sum_{n=1}^{m} \tau_n \tilde{c}_5(n) \|\Pi_{\boldsymbol{q}-1}^{\perp} \Delta_h U\|_{L^1(I_n; L^2)} \right).
\end{align*}
$$

**The final bound.** Collecting all bounds, we employ at a bound of the form $\frac{1}{2}a^2 \le b + ac$ which implies
$a \le \sqrt{2b} + 2c$ for $a, b, c \in \mathbb{R}^+$ to conclude the assertion.
$\square$

Next, we derive an upper bound for $\rho_2$ in (5.2).

**Lemma 5.2** (Upper bound on $\|\rho_2\|_{L^\infty(L^2)}$). With the notation of Lemma 5.1, we have

$$
\|\rho_2\|_{L^\infty(L^2)} \le \max_{n \in \mathcal{N}} \|[U]_{n-1}\| + \max_{n \in \mathcal{N}} \mathcal{E}(V_h^{n,\ominus}, [U]_{n-1}, [g]_{n-1}) + \max_{n \in \mathcal{N}} (c_1(q_n)c_2(q_n))^{\frac{1}{2}} \tau_n \mathcal{E}(V_h^{n,\ominus}, [U']_{n-1}, [g']_{n-1}) \\
+ \max_{n \in \mathcal{N}} (c_1(q_n)c_2(q_n))^{\frac{1}{2}} \tau_n \|[U']_{n-1}\| + \max_{n \in \mathcal{N}} \sup_{t \in I_n} \mathcal{E}(V_h^n, U(t), g(t)).
$$

*Proof.* Adding and removing terms in $\|\rho_2\|_{L^\infty(L^2)}$ leads to

$$
\|\rho_2\|_{L^\infty(L^2)} \le \|\bar{U} - \tilde{U}\|_{L^\infty(L^2)} + \|\tilde{U} - \tilde{U}'\|_{L^\infty(L^2)} + \|\tilde{U} - U\|_{L^\infty(L^2)} =: T_7 + T_8 + T_9.
$$

We estimate the three terms on the right-hand side from above.
Upper bound on $T_7$. We have

$$ 
\begin{aligned}
T_7 \stackrel{(4.10c)}{\le} \max_{n \in \mathcal{N}} \|\hat{[U]}_{n-1}\|^{(4.1),(4.2)} &\stackrel{(4.1),(4.2)}{=} \max_{n \in \mathcal{N}} \|\tilde{[U]}_{n-1}\| \le \max_{n \in \mathcal{N}} (\|[U]_{n-1}\| + \|\tilde{[U]} - [U]_{n-1}\|) \\
&\stackrel{(4.21)}{\le} \max_{n \in \mathcal{N}} (\|[U]_{n-1}\| + \mathcal{E}(V_h^{n,\ominus}, [U]_{n-1}, [g]_{n-1})).
\end{aligned}
$$

Upper bound on $T_8$. We have (4.18) to get

$$ 
\begin{aligned}
T_8 \stackrel{[38, \text{eq. (1.3)}]}{\le} & \max_{n \in \mathcal{N}} (\|\hat{U} - \tilde{U}\|_{L^2(I_n;L^2)} \|\hat{U} - \tilde{U}\|_{L^2(I_n;L^2)})^{\frac{1}{2}} \\
& \stackrel{(4.5)}{\le} \max_{n \in \mathcal{N}} (c_1(q_n)c_2(q_n))^{\frac{1}{2}} \tau_n \|\tilde{[U']}_{n-1}\| \le \max_{n \in \mathcal{N}} (c_1(q_n)c_2(q_n))^{\frac{1}{2}} \tau_n (\|\tilde{[U-U']}_{n-1}\| + \|[U']_{n-1}\|) \\
& \stackrel{(4.21)}{\le} \max_{n \in \mathcal{N}} (c_1(q_n)c_2(q_n))^{\frac{1}{2}} \tau_n (\mathcal{E}(V_h^{n,\ominus}, [U']_{n-1}, [g']_{n-1}) + \|[U']_{n-1}\|).
\end{aligned}
$$

Upper bound on $T_9$. We invoke (4.20) to find $T_9 \le \max_{n \in \mathcal{N}} \sup_{t \in I_n} \mathcal{E}(V_h^n, U(t), g(t))$. Collecting the above estimates completes the proof. □

Recalling that $\Delta_h^- = \Delta_{pw} - \Delta_h$, we define the **initial error estimator**

$$ \eta_{\text{init}} := \sqrt{2}\|u_0 - u_{0,h}\| + \sqrt{2}C_{\text{PS}}\|u_1 - u_{1,h}\| + \sqrt{2}c_{L^2}C_{\text{PS}}\left(\|\frac{h_0^2}{p^2}\Delta_h^- u_{1,h}\| + \mathcal{J}(\mathcal{F}_I^1, u_{1,h})\right); \quad (5.5) $$

the source term oscillations

$$ 
\begin{aligned}
\eta_f := 2 & \sum_{n=1}^{m} \left[ \tau_n \tilde{c}_5(n) \|\Pi_{\boldsymbol{q}-1}^{\perp} f\|_{L^1(I_n;L^2)} + c_{L^2} \left\| \frac{h_n^2}{p^2} (\Pi_{\boldsymbol{q}-1} \Pi_p^{n,\perp} f)' \right\|_{L^1(I_n;L^2)} \right. \\
& + (1 - \delta_{1,n}) \left( c_{L^2} c_1(q_n) \tau_n \left\| \frac{(h_n^{\ominus})^2}{p^2} [(\Pi_{\boldsymbol{q}-1} \Pi_p^{\perp} f)']_{n-1} \right\| + c_{L^2} \left\| \frac{(h_n^{\ominus})^2}{p^2} [\Pi_{\boldsymbol{q}-1} \Pi_p^{\perp} f]_{n-1} \right\| \right. \\
& + \left. \frac{\tau_n \tilde{\tau}_n}{4} \| [\Pi_{\boldsymbol{q}-1} \Pi_p^{\perp} f]_{n-1} \| + \tau_n^3 c_2(q_n) c_5(n) \| [(\Pi_{\boldsymbol{q}-1} \Pi_p^{\perp} f)']_{n-1} \| \right) \\
& + (2c_{L^2} c_1(q_1) + (c_1(q_1)c_2(q_1))^{\frac{1}{2}}) \tau_1 \left\| \frac{h_1^2}{p^2} (\Pi_{\boldsymbol{q}-1} \Pi_p^{\perp} f)'(t_0^+) \right\| + 3c_{L^2} \left\| \frac{h_1^2}{p^2} \Pi_{\boldsymbol{q}-1} \Pi_p^{\perp} f(t_0^+) \right\| \\
& + \frac{3\tau_1^2}{2} \| (\Pi_{\boldsymbol{q}-1} \Pi_p^{\perp} f)(t_0^+) \| + \tau_m^2 \| [\Pi_{\boldsymbol{q}-1} \Pi_p^{\perp} f]_m \| + 2\tau_1^3 c_2(q_1)c_5(1) \| (\Pi_{\boldsymbol{q}-1} \Pi_p^{\perp} f)'(t_0^+) \| \\
& + c_{L^2} \left[ \max_{n \in \mathcal{N} \setminus \{1\}} 2 \left\| \frac{(h_n^{\ominus})^2}{p^2} [\Pi_{\boldsymbol{q}-1} \Pi_p^{\perp} f]_{n-1} \right\| + \max_{n \in \mathcal{N}} (1 + \sqrt{2}\delta_{1,n}) \sup_{t \in I_n} \left\| \frac{h_n^2}{p^2} \Pi_{\boldsymbol{q}-1} \Pi_p^{\perp} f(t) \right\| \right. \\
& \left. + \max_{n \in \mathcal{N} \setminus \{1\}} (c_1(q_n)c_2(q_n))^{\frac{1}{2}} \tau_n \left\| \frac{(h_n^{\ominus})^2}{p^2} [(\Pi_{\boldsymbol{q}-1} \Pi_p^{\perp} f)']_{n-1} \right\| \right];
\end{aligned}
\quad (5.6) $$

the temporal error estimator

$$ 
\begin{aligned}
\eta_{\text{time}} := 2 & \sum_{n=1}^{m} \tau_n^3 c_2(q_n) c_5(n) \|\Delta_h U'\|_{n-1} + \max_{n \in \mathcal{N}} \tau_n (c_1(q_n)c_2(q_n))^{\frac{1}{2}} \|\tilde{[U']}_{n-1}\| \\
& + 2 \sum_{n=1}^{m} \tau_n \tilde{c}_5(n) \|\Pi_{\boldsymbol{q}-1}^{\perp} \Delta_h U\|_{L^1(I_n;L^2)};
\end{aligned}
\quad (5.7) $$
the spatial error estimator

$$
\begin{align}
\eta_{\text{space}} &:= 2c_{L^2} \sum_{n=1}^{m} \left[ \int_{I_n} \left( \|\frac{h_n^2}{p^2} \Delta_h^- U'(t)\| + \mathcal{J}(\mathcal{F}_I^n, U'(t)) \right) \mathrm{d}t \right. \nonumber \\
&\qquad \left. + c_1(q_n) \tau_n \left\| \frac{(h_n^{\ominus})^2}{p^2} [\Delta_h^- U']_{n-1} \right\| + \left\| \frac{(h_n^{\ominus})^2}{p^2} [\Delta_h^- U]_{n-1} \right\| \right. \nonumber \\
&\qquad \left. + c_1(q_n) \tau_n \mathcal{J}(\mathcal{F}_I^{n,\ominus}, [U']_{n-1}) + \mathcal{J}(\mathcal{F}_I^{n,\ominus}, [U]_{n-1}) \right] \nonumber \\
&\qquad + c_{L^2} \max_{n \in \mathcal{N}} (1 + \sqrt{2}\delta_{1,n}) \sup_{t \in I_n} \left( \|\frac{h_n^2}{p^2} \Delta_h^- U(t)\| + \mathcal{J}(\mathcal{F}_I^n, U(t)) \right) \nonumber \\
&\qquad + c_{L^2} \max_{n \in \mathcal{N}} \left\| \frac{(h_n^{\ominus})^2}{p^2} [\Delta_h^- U]_{n-1} \right\| + c_{L^2} \max_{n \in \mathcal{N}} (c_1(q_n)c_2(q_n))^{\frac{1}{2}} \tau_n \left\| \frac{(h_n^{\ominus})^2}{p^2} [\Delta_h^- U']_{n-1} \right\| \nonumber \\
&\qquad + 2c_{L^2} \max_{n \in \mathcal{N}} \mathcal{J}(\mathcal{F}_I^{n,\ominus}, [U]_{n-1}) + 2c_{L^2} \max_{n \in \mathcal{N}} (c_1(q_n)c_2(q_n))^{\frac{1}{2}} \tau_n \mathcal{J}(\mathcal{F}_I^{n,\ominus}, [U']_{n-1}); \tag{5.8}
\end{align}
$$

the mesh-change error estimator

$$
\begin{align}
\eta_{\text{mesh}} := 2 \sum_{n=1}^{m} & \left[ c_{L^2} \frac{c_3(q_n; 3)}{\tau_n} \| \frac{h_n^2}{p^2} \Pi_p^{n, \perp} U'(t_{n-1}^-) \| + \| [U]_{n-1} \| + \frac{\tau_n \tilde{\tau}_n}{4} \| [\Delta_h U]_{n-1} \| \right. \nonumber \\
& \left. + (1 - \delta_{1,n}) \left( c_{L^2} + \frac{\tau_n \tilde{\tau}_n}{4} \right) \sum_{l=n-1}^{n} \frac{c_3(q_l; 2)}{\tau_l} \| \frac{(h_n^{\ominus})^2}{p^2} \Pi_p^{l, \perp} U'(t_{l-1}^-) \| \right. \nonumber \\
& \left. + (1 - \delta_{1,n}) \left( c_{L^2} c_1(q_n) \tau_n + \tau_n^3 c_2(q_n) c_5(n) \right) \sum_{l=n-1}^{n} \frac{c_3(q_l; 3)}{\tau_l^2} \| \frac{(h_n^{\ominus})^2}{p^2} \Pi_p^{l, \perp} U'(t_{l-1}^-) \| \right] \nonumber \\
& + \tau_m^2 \left( \| [\Delta_h U]_m \| + \sum_{l=m}^{m+1} \frac{c_3(q_l; 2)}{\tau_l} \| \frac{(h_n^{\ominus})^2}{p^2} \Pi_p^{l, \perp} U'(t_{l-1}^-) \| \right) + 2 \| [U]_m \| \nonumber \\
& + c_{L^2} \left[ \max_{n \in \mathcal{N}} (1 + \sqrt{2}\delta_{1,n}) \frac{c_3(q_n; 2)}{\tau_n} \| \frac{h_n^2}{p^2} \Pi_p^{n, \perp} U'(t_{n-1}^-) \| + \max_{n \in \mathcal{N} \setminus \{1\}} \sum_{l=n-1}^{n} \frac{2c_3(q_l; 2)}{\tau_l} \| \frac{(h_n^{\ominus})^2}{p^2} \Pi_p^{l, \perp} U'(t_{l-1}^-) \| \right. \nonumber \\
& \left. + \max_{n \in \mathcal{N} \setminus \{1\}} (c_1(q_n)c_2(q_n))^{\frac{1}{2}} \tau_n \sum_{l=n-1}^{n} \frac{c_3(q_l; 3)}{\tau_l^2} \| \frac{(h_n^{\ominus})^2}{p^2} \Pi_p^{l, \perp} U'(t_{l-1}^-) \| + \max_{n \in \mathcal{N}} \| [U]_{n-1} \| \right]. \tag{5.9}
\end{align}
$$

We now have all the ingredient to prove a fully computable upper bound for the error $\|u - U\|_{L^\infty(L^2)}$, by combining the estimates for $\rho_1$ and $\rho_2$ from Lemmas 5.1 and 5.2, and the estimates for $\mathcal{E}$ from Lemma 4.3.

**Theorem 5.3 (a posteriori error bound).** With the above notation for the individual estimators, we have

$$
\|u - U\|_{L^{\infty}(L^2)} \leq \eta_{\text{init}} + \eta_f + \eta_{\text{time}} + \eta_{\text{space}} + \eta_{\text{mesh}}. \quad \square
$$

**Remark 5.4 (On the control of $\eta_{\text{mesh}}$).** The Scott-Zhang operator $\mathcal{I}_{sz}^n$ (in fact, also other classes of quasi-interpolators as in [20]) is invariant in $D \setminus D_c^n$ for functions in $V_h^{n, \oplus}$, cf. Figure 1, whence

$$
\| \Pi_p^{n, \perp} U'(t_{n-1}^-) \| \le \|(I_d - \mathcal{I}_{sz}^n) U'(t_{n-1}^-)\| = \|(I_d - \mathcal{I}_{sz}^n) U'(t_{n-1}^-)\|_{L^2(D_c^n)}.
$$

This implies that a small portions of coarsening lead to smaller $\|\Pi_p^{n,\perp} U'(t_{n-1}^-)\|$. Similarly, we have

$$
\begin{align*}
\| [U]_{n-1} \| &= \| (I_d - \mathcal{P}_p^n) U(t_{n-1}^-) \| \le C_{PS} \| \nabla (I_d - \mathcal{P}_p^n) U(t_{n-1}^-) \| \\
&\le C_{PS} \| \nabla (I_d - \mathcal{I}_{sz}^n) U(t_{n-1}^-) \| = C_{PS} \| \nabla (I_d - \mathcal{I}_{sz}^n) U(t_{n-1}^-) \|_{L^2(D_c^n)}.
\end{align*}
$$
As for the terms involving $[\Delta_h U]_{n-1}$, we have

$$ \begin{aligned} \|[\Delta_h U]_{n-1}\| &= \sup_{v \in L^2(D), \|v\|=1} ([\Delta_h U]_{n-1}, v) = \sup_{v \in L^2(D), \|v\|=1} \left( (\Delta_h U(t_{n-1}^+), \Pi_p^n v) - (\Delta_h U(t_{n-1}^-), \Pi_p^{n-1} v) \right) \\ &\stackrel{(4.12)}{=} \sup_{v \in L^2(D), \|v\|=1} \left( -(\nabla U(t_{n-1}^+), \nabla \Pi_p^n v) + (\nabla U(t_{n-1}^-), \nabla \Pi_p^{n-1} v) \right) \\ &= \sup_{v \in L^2(D), \|v\|=1} \left( \underbrace{(\nabla(I_d - \mathcal{P}_p^n)U(t_{n-1}^-), \nabla \Pi_p^n v)}_{=0} + \underbrace{(\nabla U(t_{n-1}^-), \nabla(\Pi_p^{n-1} - \Pi_p^n)v)}_{=0} \right). \end{aligned} $$

We introduce the discrete Laplacian $\Delta_h^{\oplus}$ in the space $V_h^{n, \oplus}$ by $(\Delta_h^{\oplus} U(t_{n-1}^-), W) := (\nabla U(t_{n-1}^-), \nabla W)$, for all $W \in V_h^{n, \oplus}$. We denote by $\mathcal{I}_{sz}^{\ominus, n}$ the Scott-Zhang interpolation operator on $V_h^{n, \ominus}$ and apply the $H^1$ orthogonality of $(\Pi_p^{n-1} - \Pi_p^n)v$ in $V_h^{n, \ominus}$ to find

$$ (\nabla U(t_{n-1}^-), \nabla(\Pi_p^{n-1} - \Pi_p^n)v) = (\Delta_h^{\oplus} U(t_{n-1}^-), (\Pi_p^{n-1} - \Pi_p^n)v) = ((I_d - \mathcal{I}_{sz}^{\ominus, n})\Delta_h^{\oplus} U(t_{n-1}^-), (\Pi_p^{n-1} - \Pi_p^n)v). $$

Observing that $I_d - \mathcal{I}_{sz}^{\ominus}$ is not zero only on $D_c^n \cup D_r^n$ and collecting the above displays, we get

$$ \|[\Delta_h U]_{n-1}\| \le 2 \|(I_d - \mathcal{I}_{sz}^{\ominus}) \Delta_h^{\oplus} U(t_{n-1}^-)\| = 2 \|(I_d - \mathcal{I}_{sz}^{\ominus}) \Delta_h^{\oplus} U(t_{n-1}^-)\|_{L^2(D_c^n \cup D_r^n)}. $$

Therefore, if the mesh-change is kept under control in the adaptive scheme, then the mesh-change error estimator is also kept under control.

# 6 Numerical aspects

We explore several numerical aspects of method (3.1) and the error bounds proven above. The numerical experiments are developed using the Gridap.jl library [52] in the Julia programming language.

In what follows, we consider the spatial domain $D = (-1, 1)^2$, which we partition into sequences of quasi-uniform structured and unstructured shape-regular triangular meshes, and homogeneous Dirichlet boundary conditions, which we impose strongly at the boundary degrees of freedom. We use Lagrangian uniform nodal basis functions of uniform degree $p$ in space; for the time discretization, we consider Lagrangian basis functions of degree $q$.

We consider two test cases. The first one admits a smooth solution reading

$$ u(x, y, t) := \sin(\pi x) \sin(\pi y) \cos(\sqrt{2\pi t}); \qquad (6.1) $$

the initial conditions and source term are computed accordingly. We fix the final time $T = 1$.

The second case admits an unknown in closed form solution; we solve (2.1) with data

$$ u_0(x, y) := \begin{cases} \frac{1}{c_r \epsilon} \exp(\frac{1}{r^2/\epsilon^2 - 1}) & \text{if } r < \epsilon \\ 0 & \text{otherwise,} \end{cases} \quad u_1(x, y) := \begin{cases} 0.9 - r^2 & \text{if } r < 1 \\ 0 & \text{otherwise,} \end{cases}, \quad f(x, y, t) := 0, \quad (6.2) $$

for $\epsilon = 0.1$, $r = \sqrt{x^2 + y^2}$, $c_r = 0.4439938161680794$, and $T = 0.5$. The initial condition $u_0$ is analytic and has compact support $\{r \le \epsilon\}$, while $u_1$ is discontinuous.

## 6.1 Energy dissipation

Before studying the quality of the a posteriori error bounds and their use within an adaptive algorithm below, we begin by assessing the energy dissipation properties of the method (3.1). To that end, we set $f = 0$ and introduce the total energy

$$ E(t, v) := \frac{1}{2} \|\partial_t v(t)\|^2 + \frac{1}{2} \|\nabla v(t)\|^2. $$
Owing to the dissipative nature of nonconforming methods in time, we expect energy dissipation rather than energy conservation (i.e., $E(t, u) = E(0, u)$); cf. [54, eq. (4.1)]. In Table 1, we assess the influence of $q$ on the energy dissipation; we select $h = \tau = 10^{-2}$, $p = 1$, $q = 2, \dots, 5$, a uniform static mesh, and the test case with finite Sobolev regularity exact solution with data as in (6.2).

<table><caption>Table 1: Relative numerical energy dissipation on static meshes with $h = \tau = 10^{-2}$, $p = 1$, and $q = 2, \dots, 5$.</caption><thead><tr><th>q</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>E(0,U) - E(T,U) / E(0,U)</td><td>1.79e-2</td><td>1.33e-4</td><td>5.48e-7</td><td>1.63e-9</td></tr></tbody></table>

From Table 1, we see that the energy dissipation vanishes exponentially with respect to $q$: although the scheme is formally dissipative, dissipation appears to be negligible. As we will see below, when dynamic mesh modification is considered, energy dissipation is more apparent due to the loss of information when projecting the discrete solution from a time slab to the next one, cf. [18].

## 6.2 Effectivity of the error estimators

In practice, the time instance when the error attains its maximum, cf. (4.23), is unknown. We employ an adaptive algorithm to determine a practical value of $m$ and set $\xi = t_{m-1}$, using the error indicator $\tilde{\eta}_n := \tau_n^3 c_2(q_n)c_5(n) \| [\Delta_h U']_{n-1} \| + \tau_n \tilde{c}_5(n) \| \Pi_{q-1}^\perp \Delta_h U \|_{L^1(I_n; L^2)}$. The procedure is as follows: 1) For $n = 1$, set $m = 1$, $\xi = 0$, and compute $\tilde{\eta}_1$. 2) For $n = 2, \dots, N$, compute $\tilde{\eta}_n$. If $\tilde{\eta}_n \ge \tilde{\eta}_m$, then update $m = n$; otherwise, leave $m$ unchanged.

We assess the efficiency of the temporal (5.7) and spatial (5.8) error estimators for the test case with analytic exact solution (6.1). To this end, we introduce the effectivity indices

$$ \kappa_{time} = \frac{\eta_{time}}{\|e\|_{L^{\infty}(L^2)}}, \qquad \kappa_{space} = \frac{\eta_{space}}{\|e\|_{L^{\infty}(L^2)}}. \qquad (6.3) $$

Here and below, we set $c_{L^2} = \frac{1}{10}$, $\frac{1}{30}$ and $\frac{1}{40}$ for $p = 1, 2, 3$, respectively; these choices of $c_{L^2}$ are based on our empirical experience for elliptic problems. We consider polynomial degrees $p = q = 2, 3$, and uniformly-refined sequences of spatial (either structured or quasi-uniform) meshes. We tested approximation of the initial conditions $u_0$ and $u_1$, and the source term $f$ with either their Lagrangian interpolants or $L^2$-projections. Based on the numerical results, both choices of meshes and all choices of data approximation lead to similar behaviour. For brevity, we only present the numerical results based on the structured meshes and interpolated data.

In Figure 2, we illustrate the $L^\infty(L^2)$-errors along with the estimators $\eta_{space}$ and $\eta_{time}$ in (5.7) and (5.8) together with their averaged slopes on the left-panel, and the effectivity indices in (6.3) on the right-panel.

From the results in Figure 2, it is apparent that the estimators $\eta_{space}$ and $\eta_{time}$ converge optimally. We also report (without presenting the results for brevity) that using either the Lagrangian interpolation or the $L^2$-projection of the initial conditions and the source term, the effectivity indices in (6.3) do not essentially grow or grow very moderately. The above results are in contrast with [31, Tables 7 and 8], where Lagrangian interpolation results to numerical results with effectivity indices growing considerably with mesh refinements. We conjecture that the estimator behaviour in [31, Tables 7 and 8] is due to the different a posteriori error estimator used there.

In Table 2, we investigate numerically the behaviour of the errors and estimators for the $(p, q)$-version of the method. We only consider the case of Lagrangian interpolated data, a uniform spatial mesh with 6 nodes along each direction, and set a uniform time-step equal to $h$; similar results have been obtained using the $L^2$-projection of the data and the unstructured meshes (omitted for brevity).

The $L^\infty(L^2)$-error and space/time estimators converge exponentially in terms of $p$ and $q$. We observe that the effectivity indices appear to grow fairly slowly. One reason may be that the error estimator from Section 5 is slightly subopitmal in $q$ by one order; see Remark 4.8.
<figure><img src="image_3.png" /><figcaption>Figure 2: Smooth solution in (6.1); sequences of uniformly-refined structured triangular meshes; $p = q = 2, 3$. Left-panel: $L^\infty(L^2)$-error, and estimators in (5.7) and (5.8). Right-panel: effectivity indices in (6.3).</figcaption></figure>

## 6.3 Effect of mesh modification

We investigate the influence of the mesh-change on the performance of the error and error estimators. We consider the test case with the smooth solution in (6.1). We consider two structured meshes $T_1$ and $T_2$ with $N_1 = 10$ and $N_2 = 5, 6, \dots, 30$ nodes in each direction, respectively. We select $T_1$ and $T_2$ on the time intervals $I_n$ with odd and even indices $n$, respectively; the uniform time-step is set to $\tau = \frac{1}{N_1}$. The operator $P_p^n$ in (3.1) is selected to be the Lagrange interpolant. The results are shown in Figure 3.

From Figure 3, it is apparent that the mesh-change has a significant influence on the behaviour of the $L^\infty(L^2)$-error, which matches with the theoretical findings in [41, Theorem 3.2] for the nonlinear wave equation.

Next, we assess the practical behaviour of the estimators $\eta_{space}$, $\eta_{time}$, and $\eta_{mesh}$. We pick $N_1 \in \{5, 8, 11, 14, 17\}$, $N_2 = N_1 - 1$, $\tau = \frac{1}{N_1}$, and $p = q = 2, 3, 4$, $P_p^n$ selected either as the Lagrange interpolant, $L^2$-, or $H^1$-projections. In this case, we have $T_h^{n, \ominus} = \{D\}$, $F_I^{n, \ominus} = \emptyset$ and $h_n^{\ominus} = 2\sqrt{2}$. For the sake of brevity, we do not report the numerical results and only state the conclusions here. Owing to the mesh change, the error in the $L^\infty(L^2)$-norm only converges of order $O(\tau^q)$, i.e., we lose one order of convergence, with the
<table><caption>Table 2: Smooth solution in (6.1) with $p = 2, ..., 7$, $q = 2, ..., 7$. $L^\infty(L^2)$-error (first table); estimator $\eta_{space}$ (5.8) (second table); estimator $\eta_{time}$ (5.7) (third table); effectivity indices (6.3) computed from the diagonal entries of the other three tables (fourth table).</caption><thead><tr><th></th><th>error</th><th>q = 2</th><th>q = 3</th><th>q = 4</th><th>q = 5</th><th>q = 6</th><th>q = 7</th></tr></thead><tbody><tr><td></td><td>p = 2</td><td>6.08e-2</td><td>4.99e-2</td><td>4.99e-2</td><td>4.99e-2</td><td>4.99e-2</td><td>4.99e-2</td></tr><tr><td></td><td>p = 3</td><td>2.12e-2</td><td>4.98e-3</td><td>5.19e-3</td><td>5.40e-3</td><td>5.44e-3</td><td>5.47e-3</td></tr><tr><td></td><td>p = 4</td><td>2.12e-2</td><td>9.00e-4</td><td>6.77e-4</td><td>6.84e-4</td><td>6.89e-4</td><td>6.87e-4</td></tr><tr><td></td><td>p = 5</td><td>2.12e-2</td><td>6.33e-4</td><td>7.12e-5</td><td>6.85e-5</td><td>6.99e-5</td><td>7.08e-5</td></tr><tr><td></td><td>p = 6</td><td>2.12e-2</td><td>6.22e-4</td><td>2.20e-5</td><td>6.82e-6</td><td>6.83e-6</td><td>7.15e-6</td></tr><tr><td></td><td>p = 7</td><td>2.12e-2</td><td>6.21e-4</td><td>2.13e-5</td><td>1.09e-6</td><td>5.07e-7</td><td>5.31e-7</td></tr><tr><td></td><td>η<sub>space</sub></td><td>q = 2</td><td>q = 3</td><td>q = 4</td><td>q = 5</td><td>q = 6</td><td>q = 7</td></tr><tr><td></td><td>p = 2</td><td>1.07e0</td><td>1.31e0</td><td>1.47e0</td><td>1.55e0</td><td>1.56e0</td><td>1.56e0</td></tr><tr><td></td><td>p = 3</td><td>7.79e-2</td><td>9.65e-2</td><td>1.70e-1</td><td>1.99e-1</td><td>2.27e-1</td><td>2.39e-1</td></tr><tr><td></td><td>p = 4</td><td>6.33e-2</td><td>6.68e-2</td><td>9.02e-3</td><td>1.53e-2</td><td>1.87e-2</td><td>2.18e-2</td></tr><tr><td></td><td>p = 5</td><td>6.68e-4</td><td>5.51e-4</td><td>7.03e-4</td><td>1.03e-3</td><td>1.82e-3</td><td>1.79e-3</td></tr><tr><td></td><td>p = 6</td><td>6.41e-5</td><td>4.55e-5</td><td>5.18e-5</td><td>6.74e-5</td><td>1.04e-4</td><td>1.45e-4</td></tr><tr><td></td><td>p = 7</td><td>1.09e-6</td><td>4.25e-6</td><td>3.75e-6</td><td>4.72e-6</td><td>5.92e-6</td><td>1.05e-5</td></tr><tr><td></td><td>η<sub>time</sub></td><td>q = 2</td><td>q = 3</td><td>q = 4</td><td>q = 5</td><td>q = 6</td><td>q = 7</td></tr><tr><td></td><td>p = 2</td><td>8.12e-1</td><td>1.28e-1</td><td>4.77e-2</td><td>9.28e-3</td><td>1.34e-3</td><td>1.49e-4</td></tr><tr><td></td><td>p = 3</td><td>7.04e-1</td><td>9.29e-2</td><td>5.04e-2</td><td>2.29e-2</td><td>5.79e-3</td><td>1.43e-3</td></tr><tr><td></td><td>p = 4</td><td>3.86e-1</td><td>3.07e-2</td><td>1.37e-2</td><td>8.56e-3</td><td>3.55e-3</td><td>1.17e-3</td></tr><tr><td></td><td>p = 5</td><td>3.09e-1</td><td>1.53e-2</td><td>3.26e-3</td><td>2.02e-3</td><td>1.22e-3</td><td>5.03e-4</td></tr><tr><td></td><td>p = 6</td><td>2.97e-1</td><td>1.18e-2</td><td>9.33e-4</td><td>3.74e-4</td><td>1.99e-4</td><td>1.17e-4</td></tr><tr><td></td><td>p = 7</td><td>2.95e-1</td><td>1.13e-2</td><td>4.75e-4</td><td>6.88e-5</td><td>3.01e-5</td><td>1.80e-5</td></tr><tr><td></td><td>polynomial degree</td><td>2</td><td>3</td><td>4</td><td>5</td><td>6</td><td>7</td></tr><tr><td></td><td>κ<sub>space</sub></td><td>17.59</td><td>19.37</td><td>13.32</td><td>15.03</td><td>15.22</td><td>19.74</td></tr><tr><td></td><td>κ<sub>time</sub></td><td>13.35</td><td>18.65</td><td>20.23</td><td>29.48</td><td>29.13</td><td>33.89</td></tr></tbody></table>

exception of the following terms:

$$ \eta_{space,1} := \sum_{n=1}^{m} \tau_n c_1(q_n) \left( \left\| \frac{(h_n^{\ominus})^2}{p^2} [\Delta_h^{-} U']_{n-1} \right\| + \mathcal{J}(\mathcal{F}_I^{n,\ominus}, [U']_{n-1}) \right) \quad (6.4a) $$

$$ \eta_{space,2} := \sum_{n=1}^{m} \left( \left\| \frac{(h_n^{\ominus})^2}{p^2} [\Delta_h^{-} U]_{n-1} \right\| + \mathcal{J}(\mathcal{F}_I^{n,\ominus}, [U]_{n-1}) \right) \quad (6.4b) $$

$$ \eta_{mesh,1} := \sum_{n=1}^{m} \tau_n^{-1} c_3(q_n; 3) \left\| \frac{(h_n^{\ominus})^2}{p^2} \Pi_p^{n,\perp} U'(t_{n-1}^{-}) \right\|. \quad (6.4c) $$

In Table 3, we present the values and rates for each term in (6.4), with $\mathcal{P}_p^n$ in (3.1) selected to be the Lagrange interpolant. The $L^2$- and $H^1$-projectors have also been tested and they deliver results similar to those obtained using the Lagrange interpolant (omitted for brevity).

From Table 3, we observe that $\eta_{space,i}$, $i = 1, 2$, and $\eta_{mesh,1}$ are close to $O(\tau^{q-2})$, and $\eta_{mesh,1}$ is dominant. Finally, in order to better understand the different influence of $\mathcal{P}_p^n$ on the above estimators, in Table 4 we show the values and rates of $\eta_{mesh,1}$ for three choices of $\mathcal{P}_p^n$, with $p = q = 3$. The other terms are also tested and the corresponding, analogous results are omitted for brevity. From Table 4, we observe that the values of the estimators obtained using the $H^1$-projector are smaller than those obtained using the other two projections.
<figure><img src="image_4.png" /><figcaption>Figure 3: Errors for the smooth solution in (6.1) on sequences of meshes with number of nodes $N_1$ and $N_2$ in each spatial direction for odd and even indices $n$, respectively, with fixed mesh-size and time-step, and $p = q = 2, \dots, 5$.</figcaption></figure>

## 6.4 A space-time adaptive scheme

We propose and assess a space-time adaptive algorithm for (3.1). To minimize the computational cost of computing the estimators, we omit the data oscillation term $\eta_f$, which appears to be of higher order in the examples selected, and we focus on the the dominant and “easy-to-compute” terms. In particular, the space-time adaptive algorithm does not take into account the mesh-change estimator $\eta_{\text{mesh}}^n$. We use $\eta_{\text{time}}$ to refine/coarsen in time. $\eta_{\text{space}}$ is used to refine/coarsen in time in space: elements are refined via the newest-vertex-bisection strategy, coarsening strategy follows that from [13, Sect. 3].

The adaptive algorithm, inspired by the algorithms in [13, 26, 10], employs localised error indicators associated with the error estimators in (5.5), (5.7), (5.8), and (5.9). Setting

$$ \mathcal{H}(h, p, W) := \|\frac{h^2}{p^2}(\Delta - \Delta_h)W\|_{L^2(K)} + \sum_{F \in \mathcal{F}_K} \|\frac{h^{\frac{3}{2}}}{p^{\frac{3}{2}}}[\mathbf{n}_F \cdot \nabla W]\|_{L^2(F)}, $$

for brevity, we define the local initial error indicators

$$ \eta_{0,K} := \|u_0 - u_{0,h}\|_{L^2(K)} + C_{PS}\|u_1 - u_{1,h}\|_{L^2(K)} + c_{L^2}C_{PS}\mathcal{H}(h_0, p, u_{1,h}), \quad (6.5) $$

$K \in \mathcal{T}_h^0$ and their global counterpart $\eta_0 = (\sum_{K \in \mathcal{T}_h^0} \eta_{0,K}^2)^{\frac{1}{2}}$; the time error indicators

$$ \eta_{\text{time}}^n := \tau_n(c_1(q_n)c_2(q_n))^{\frac{1}{2}} \|[U']_{n-1}\| \quad n \in \mathcal{N}; $$

the spatial error indicators

$$ \eta_{\text{space},K,1}^n := \sup_{t \in I_n} \mathcal{H}(h_n, p, U), \quad \eta_{\text{space},K,2}^n := c_1(q_n)\mathcal{H}(h_n^\ominus, p, U'), \quad \eta_{\text{space},K,3}^n := \mathcal{H}(h_n^\ominus, p, \tau_n^{-1}[U]_{n-1}), $$
<table><caption>Table 3: Smooth solution in (6.1); “switching” meshes with $N_1 \in \{5, 8, 11, 14, 17\}$ and $N_2 = N_1 - 1$ nodes in each spatial direction, respectively. Top table: $p = q = 2$. Bottom table: $p = q = 3$.</caption><thead><tr><th>value (rate) with p = q = 2</th><th>mesh 1</th><th>mesh 2</th><th>mesh 3</th><th>mesh 4</th><th>mesh 5</th></tr></thead><tbody><tr><th>η<sub>space,1</sub></th><td>8.07e0</td><td>8.19e0(0.06)</td><td>8.36e0(-0.01)</td><td>8.29e0(0.06)</td><td>8.30e0(-0.01)</td></tr><tr><th>η<sub>space,2</sub></th><td>1.23e1</td><td>1.28e1(-0.17)</td><td>1.33e1(-0.20)</td><td>1.33e1(0.01)</td><td>1.34e1(-0.05)</td></tr><tr><th>η<sub>mesh,1</sub></th><td>1.30e2</td><td>1.32e2(-0.06)</td><td>1.39e2(-0.29)</td><td>1.37e2(0.08)</td><td>1.39e2(-0.11)</td></tr></tbody></table><table><thead><tr><th>value (rate) with p = q = 3</th><th>mesh 1</th><th>mesh 2</th><th>mesh 3</th><th>mesh 4</th><th>mesh 5</th></tr></thead><tbody><tr><th>η<sub>space,1</sub></th><td>1.37e1</td><td>8.39e-1(1.10)</td><td>6.24e-1(0.96)</td><td>5.24e-1(0.74)</td><td>4.22e-1(1.13)</td></tr><tr><th>η<sub>space,2</sub></th><td>2.10e0</td><td>1.66e0(0.52)</td><td>1.35e0(0.67)</td><td>1.22e0(0.40)</td><td>1.01e0(1.02)</td></tr><tr><th>η<sub>mesh,1</sub></th><td>1.06e2</td><td>6.90e1(0.97)</td><td>3.54e1(2.17)</td><td>2.67e1(1.19)</td><td>1.95e1(1.64)</td></tr></tbody></table>

<table><caption>Table 4: Smooth solution in (6.1); “switching” meshes with $N_1 \in \{5, 8, 11, 14, 17\}$ and $N_2 = N_1 - 1$ nodes in each spatial direction and various choices of $P_p^n$, $p = q = 3$.</caption><thead><tr><th>η<sub>mesh,1</sub> with p = q = 3</th><th>mesh 1</th><th>mesh 2</th><th>mesh 3</th><th>mesh 4</th><th>mesh 5</th></tr></thead><tbody><tr><th>interpolator</th><td>1.06e2</td><td>6.90e1(0.97)</td><td>3.54e1(2.17)</td><td>2.67e1(1.19)</td><td>1.95e1(1.64)</td></tr><tr><th>L<sup>2</sup>-projector</th><td>1.21e2</td><td>8.91e1(0.69)</td><td>5.03e1(1.86)</td><td>4.08e1(0.89)</td><td>3.06e1(1.51)</td></tr><tr><th>H<sup>1</sup>-projector</th><td>6.25e1</td><td>3.91e1(1.05)</td><td>1.93e1(2.30)</td><td>1.35e1(1.51)</td><td>9.29e0(1.97)</td></tr></tbody></table>

collected in $η_{space,K}^n := η_{space,K,1}^n + η_{space,K,2}^n + η_{space,K,3}^n$, $n \in \mathcal{N}$, $K \in \mathcal{T}_h^n$, and their global counterpart

$$η_{space}^n := c_{L^2} \left( \left( \sum_{K \in \mathcal{T}_h^n} (η_{space,K,1}^n)^2 \right)^{\frac{1}{2}} + \sum_{k=1}^{n} \left( \sum_{K \in \mathcal{T}_h^n} (η_{space,K,2}^n + η_{space,K,3}^n)^2 \right)^{\frac{1}{2}} \right).$$

We employ Dörfler’s marking strategy [16] to select elements for refinement and/or coarsening, given parameters $\theta_c$ (coarsening) and $\theta_r$ (refinement). Specifically, we mark two sets of elements with minimal cardinality such that the sum of the error estimator contributions over each set exceeds $\theta_r$ or $\theta_c$ times the total error estimator, respectively. The pseudocode of the proposed adaptive procedure is given in Algorithm 1.

## 6.5 Numerical results for the adaptive scheme

We assess the performance of Algorithm 1 for the test case with finite Sobolev regularity in (6.2). We construct a reference solution by employing a static uniform spatial mesh with 309123 space-time DoFs at each time step, $\tau = h$, $p = 1$, and $q = 2$.

We set $p = 1$, $q = 2$, $\delta_0 = \delta_{time} = \delta_{space} = \delta_{mesh} = 10^{-2}$, $\theta_r = 0.1$, $\theta_c = 0.1$, $\delta_{mesh} = \delta_{space}/2$, $k_{time} = 10$, $k_{space} = 10$, and we consider the $H^1$-projection of the initial data at each time-step. In Figure 4, the numerical solution with around $4 \times 10^4$ space-time DoFs is compared with the reference solution above. In the first row, the initial mesh and the initial data are depicted; in the second row, we provide the mesh, and the numerical and reference solutions at the final time; the third row contains the number of space-time DoFs in each time interval, the size of all time-steps, and the estimators at the time nodes, respectively. We also tested the algorithm for $q = 3, 4$, which lead to similar results as those obtained with $q = 2$; in particular, $q = 2, 3, 4$ lead to time-steps of approximate sizes 0.015, 0.0225, 0.035, respectively.

From Figure 4, since the initial condition $u_1$ in (6.2) has a jump on $\{\sqrt{x^2 + y^2} = 1\}$ and $u_0$ in (6.2) has compact support located at $(0, 0)$, we deduce that the initial mesh $\mathcal{T}_h^0$ contains cells galore next to those curve and point. At the final time $T$, the reference solution displays the most variation around a circle centred at $(0, 0)$, whence several cells are located near that circular layer. The number of space-time DoFs in each time interval is decreasing for small time indices $n$ and then it increases, while the time-steps do not significantly change during the evolution. In Table 5, we provide the relative energy dissipation produced by the adaptive algorithm with dynamic mesh modification and a corresponding one without mesh modification having comparable total space-time degrees of freedom. We stress that the time steps in this example are
<caption>**Algorithm 1** Space-time adaptive algorithm</caption>

**Require:** Tolerances $\delta_0$, $\delta_{space}$, $\delta_{mesh}$, $\delta_{time} > 0$, coarsening parameter $\theta_c \in (0, 1)$, refinement parameter $\theta_r \in (0, 1)$, initial mesh $\mathcal{T}_h^0$, final time $T$, initial time-step $\tau_1$, maximal $\#$ of attempts $k_{time}, k_{space}$.

1: **while** $\eta_0 > \delta_0$ **do**
2: &nbsp;&nbsp;&nbsp;&nbsp;refine $\mathcal{T}_h^0$ based on $\eta_{0,K}$ and Dörfler's marking with $\theta_r$.
3: **end while**
4: set $t = 0, n = 1$.
5: **while** $t < T$ **do**
6: &nbsp;&nbsp;&nbsp;&nbsp;set $k \leftarrow 0, \mathcal{T}_h^n \leftarrow \mathcal{T}_h^{n-1}$.
7: &nbsp;&nbsp;&nbsp;&nbsp;update $U|_{I_n}$ with time-step $\tau_n$ and mesh $\mathcal{T}_h^n$.
8: &nbsp;&nbsp;&nbsp;&nbsp;**if** $\eta_{time}^n > \delta_{time}, k < k_{time}$ **then**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\tau_n \leftarrow \frac{\tau_n}{2}, k \leftarrow k+1$, **go to** 7.
10: &nbsp;&nbsp;&nbsp;&nbsp;**else if** $\eta_{time}^n < \frac{\delta_{time}}{2}, k < k_{time}$ **then**
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\tau_n \leftarrow \min(\frac{3\tau_n}{2}, T-t), k \leftarrow k+1$, **go to** 7.
12: &nbsp;&nbsp;&nbsp;&nbsp;**end if**
13: &nbsp;&nbsp;&nbsp;&nbsp;set $k \leftarrow 0$.
14: compute $U|_{I_n}$ with time-step $\tau_n$ and mesh $\mathcal{T}_h^n$.
15: &nbsp;&nbsp;&nbsp;&nbsp;**if** $c_{L^2}(\sum_{K \in \mathcal{T}_h^n}(\eta_{space,K}^n)^2)^{\frac{1}{2}} > \delta_{space}$ and $k < k_{space}$ **then**
16: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;refine $\mathcal{T}_h^n$, based on $\eta_{space,K}^n$ and Dörfler's marking with $\theta_r$.
17: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;compute $U|_{I_n}$ with time-step $\tau_n$ and mesh $\mathcal{T}_h^n$.
18: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$k \leftarrow k+1$.
19: &nbsp;&nbsp;&nbsp;&nbsp;**end if**
20: coarsen $\mathcal{T}_h^n$, based on $\eta_{space,K}^n$ and Dörfler's marking with $\theta_c$.
21: &nbsp;&nbsp;&nbsp;&nbsp;**if** $(\sum_{K \in \mathcal{T}_h^n}(\eta_{space,K}^n)^2)^{\frac{1}{2}} > \delta_{mesh}$ **then**
22: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\theta_c \leftarrow \frac{\theta_c}{2}$, **go to** and **redo** 20.
23: &nbsp;&nbsp;&nbsp;&nbsp;**end if**
24: compute $U|_{I_n}$ with time-step $\tau_n$ and mesh $\mathcal{T}_h^n$.
25: $t \leftarrow t + \tau_n, \tau_{n+1} \leftarrow \tau_n, n \leftarrow n+1$.
26: **end while**

$\triangleright$ control initial error

$\triangleright$ trial update

$\triangleright$ time adaptivity

$\triangleright$ initialize time parameters

$\triangleright$ initialize refinement parameters

$\triangleright$ initialize mesh refinement parameters

$\triangleright$ initialize mesh coarsening parameters

$\triangleright$ mesh refinement

$\triangleright$ mesh coarsening

$\triangleright$ control coarsening error

$\triangleright$ definitive update

$\triangleright$ go to next time-step

considerably larger than the test from Table 1, leading to more excessive energy dissipation for the same numerical test case. We note, however, that the overall impact on the energy conservation is not dramatic. For $q = 2$, the mesh modification increases energy dissipation by about 3%; for $q = 3$ by about 70%; for $q = 4$ by about 89%.

<table><caption>Table 5: Relative numerical energy dissipation on adaptive meshes and static meshes with, $p = 1$, and $q = 2, \dots, 4$.</caption><thead><tr><th>E(0,U)-E(T,U)/E(0,U)</th><th>q = 2</th><th>q = 3</th><th>q = 4</th></tr></thead><tbody><tr><td>Adaptive algorithm with mesh modification</td><td>1.99e-1</td><td>8.99e-2</td><td>7.38e-2</td></tr><tr><td>Adaptive algorithm without mesh modification</td><td>1.42e-1</td><td>5.26e-2</td><td>3.89e-2</td></tr></tbody></table>

Finally, in Figure 5, we illustrate the influence of the spatial polynomial degree $p$ in our adaptive scheme with approximately the same number of spatial elements in the initial and the final meshes: around 8000 and 5700 elements, respectively. We fix $q = 2$, $k_{space} = k_{time} = 10$, $\theta_c = \theta_r = 0.1$ and $\delta_{mesh} = 3\delta_{space}$. For $p = 1$, we choose $\delta_0 = 1.05 \times 10^{-2}$, $\delta_{time} = \delta_{space} = 3.9 \times 10^{-2}$; for $p = 2$, we choose $\delta_0 = 5 \times 10^{-3}$, $\delta_{time} = \delta_{space} = 2 \times 10^{-3}$.

From Figure 5, we observe that the case $p = 2$ captures the layers better than the case $p = 1$ with
<figure><img src="image_5.png" /><figcaption>Figure 4: Visualization of numerical solutions of the test case with exact solution with data (6.2) computed with Algorithm 1; $p = 1$; $q = 2$ (the other parameters are given in the text). (a): initial mesh $\mathcal{T}_h^0$ generated adaptively based on the indicator in (6.5). (b) and (c): Lagrangian interpolants $u_{0,h}$ and $u_{1,h}$ of $u_0$ and $u_1$ on the mesh $\mathcal{T}_h^0$. (d): mesh $\mathcal{T}_h^N$ produced by Algorithm 1. (e): adaptive solution $U(T^{-})$ on $\mathcal{T}_h^N$. (f): reference solution computed using 309123 space-time DoFs at each time-step on a fine static uniform spatial mesh. (g): evolution of the number of space-time DoFs in all time intervals using the adaptive algorithm. (h): evolution of time-steps during the adaptive process. (i): evolution of the estimators during the adaptive process.</figcaption></figure>

approximately the same number of cells. Further, we observe that different $p$ lead to almost the same time-steps, which we omit here for brevity, owing to the unconditional stability of the underlying space-time method. The evolution of the space-time DoFs has the same behaviour as in Figure 4 (g).
<figure><img src="image_6.png" /><figcaption>Figure 5: Visualization of dynamic mesh modification of the test case with exact solution with data as in (6.2) computed with Algorithm 1; $q = 2$; $p$ equals 1 and 2 for top figures and bottom figures, respectively; $T$ equals 0, 0.1, 0.3 and 0.5 from left to right.</figcaption></figure>

# 7 Conclusions

This work has been concerned with the proof of constant explicit *a posteriori* upper bounds in the $L^\infty(L^2)$-norm for the wave equation, with a particular emphasis on dynamic mesh modification. The *a posteriori* error estimates required the design and analysis of novel reconstructions, including a Hermite-type in time and an elliptic reconstructions. The efficiency of an error estimator has been validated in practice on different test cases; its proof is particularly challenging and is postponed to future investigations. The method is non-conforming in time so as to permit dynamic mesh changes; this comes at the price that it dissipates the energy. We verified practically that for given spatial meshes and time steps, higher order methods are less dissipative. The proposed error estimator consists of several distinct terms (dealing with data oscillations; spatial mesh refinement/coarsening; time refinement/coarsening), some of which were used to propose a full space-time adaptive algorithm with dynamic mesh modification; the corresponding effect within that adaptive algorithm, as well as the inclusion of local time-stepping, are theoretical challenges that deserve special attention and will be considered in the future. Upon minor modifications, the analysis presented above can be extended to other second order in time problems, e.g., wave problems with a variable diffusion coefficient in space, damped wave problems, wave problems with other linear elliptic operators, ...

**Acknowledgements.** ZD has been partially funded by the French National Research Agency (ANR, STEERS, project number ANR-24-CE56-0127-01). EHG gratefully acknowledges the financial support of EPSRC (grant number EP/W005840/2). LM has been partially funded by the European Union (ERC, NEMESIS, project number 101115663); views and opinions expressed are however those of the author only and do not necessarily reflect those of the EU or the ERC Executive Agency. LM has been partially funded by MUR (PRIN2022 research grant n. 202292JW3F). LM is also member of the Gruppo Nazionale Calcolo Scientifico-Istituto Nazionale di Alta Matematica (GNCS-INdAM). ZW would like to thank Dr. Ari Rappaport for the kind discussions concerning implementation details.
References

[1] S. Adjerid. A posteriori finite element error estimation for second-order hyperbolic problems. *Comput. Methods Appl. Mech. Engrg.*, 191(41-42):4699–4719, 2002.

[2] G. A. Baker. Error estimates for finite element methods for second order hyperbolic equations. *SIAM J. Numer. Anal.*, 13(4):564–576, 1976.

[3] W. Bangerth, M. Geiger, and R. Rannacher. Adaptive Galerkin finite element methods for the wave equation. *Comput. Methods Appl. Math.*, 10(1):3–48, 2010.

[4] W. Bangerth and R. Rannacher. Finite element approximation of the acoustic wave equation: error control and mesh adaptation. *East-West J. Numer. Math.*, 7(4):263–282, 1999.

[5] W. Bangerth and R. Rannacher. Adaptive finite element techniques for the acoustic wave equation. *J. Comput. Acoust.*, 9(2):575–591, 2001.

[6] L. Banjai, E. H. Georgoulis, and B. Hennessy. Space-time finite element methods for nonlinear wave equations via elliptic regularisation. *arXiv preprint arXiv:2507.22757*, 2025.

[7] L. Banjai, E. H. Georgoulis, and O. Lijoka. A Trefftz polynomial space-time discontinuous Galerkin method for the second order wave equation. *SIAM J. Numer. Anal.*, 55(1):63–86, 2017.

[8] P. Bansal, A. Moiola, I. Perugia, and Ch. Schwab. Space-time discontinuous Galerkin approximation of acoustic waves with point singularities. *IMA J. Numer. Anal.*, 41(3):2056–2109, 2021.

[9] Ch. Bernardi and E. Süli. Time and space adaptivity for the second-order wave equation. *Math. Models Methods Appl. Sci.*, 15(02):199–225, 2005.

[10] A. Cangiani, E. H. Georgoulis, and O. J. Sutton. Adaptive non-hierarchical Galerkin methods for parabolic problems with application to moving mesh and virtual element methods. *Math. Models Methods Appl. Sci.*, 31(4):711–751, 2021.

[11] T. Chaumont-Frelet. Asymptotically constant-free and polynomial-degree-robust a posteriori estimates for space discretizations of the wave equation. *SIAM J. Sci. Comput.*, 45(4):A1591–A1620, 2023.

[12] T. Chaumont-Frelet and A. Ern. Damped Energy-norm *a posteriori* error estimates using $C^2$-reconstructions for the fully discrete wave equation with the leapfrog scheme. *ESAIM Math. Model. Numer. Anal.*, 59(4):1937–1972, 2025.

[13] L. Chen and C. Zhang. A coarsening algorithm on adaptive grids by newest vertex bisection and its applications. *J. Comput. Math.*, 28(6):767–789, 2010.

[14] J. Diaz and M. J. Grote. Energy conserving explicit local time stepping for second-order wave equations. *SIAM J. Sci. Comput.*, 31(3):1985–2014, 2009.

[15] Z. Dong, L. Mascotto, and Z. Wang. A priori and a posteriori error estimates of a DG-CG method for the wave equation in second order formulation. http://arxiv.org/abs/2411.03264, 2024.

[16] W. Dörfler. A convergent adaptive algorithm for Poisson's equation. *SIAM J. Numer. Anal.*, 33(3):1106–1124, 1996.

[17] W. Dörfler, S. Findeisen, and Ch. Wieners. Space-time discontinuous Galerkin discretizations for linear first-order hyperbolic evolution systems. *Comput. Methods Appl. Math.*, 16(3):409–428, 2016.

[18] T. Dupont. $L^2$-estimates for Galerkin methods for second order hyperbolic equations. *SIAM J. Numer. Anal.*, 10:880–889, 1973.

[19] K. Eriksson, D. Estep, P. Hansbo, and C. Johnson. Introduction to adaptive methods for differential equations. *Acta Numer.*, pages 105–158, 1995.

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

[21] L. C. Evans. *Partial Differential Equations*, volume 19. American Mathematical Society, 2022.

[22] M. Ferrari and S. Fraschini. Stability of conforming space-time isogeometric methods for the wave equation. *published online, Math. Comp.*, 2025.

[23] M. Ferrari, S. Fraschini, G. Loli, and I. Perugia. Unconditionally stable space-time isogeometric discretization for the wave equation in hamiltonian formulation. *arXiv preprint arXiv:2411.00650*, 2024.
[24] D. A. French. A space-time finite element method for the wave equation. *Comput. Methods Appl. Mech. Engrg.*, 107(1-2):145-157, 1993.

[25] D. A. French and T. E. Peterson. A continuous space-time finite element method for the wave equation. *Math. Comp.*, 65(214):491–506, 1996.

[26] F. D. Gaspoz, K. Siebert, Ch. Kreuzer, and D. A. Ziegler. A convergent time-space adaptive dG(s) finite element method for parabolic problems motivated by equal error distribution. *IMA J. Numer. Anal.*, 39(2):650–686, 2019.

[27] E. H. Georgoulis, O. Lakkis, and Ch. Makridakis. A posteriori $L^\infty(L^2)$-error bounds for finite element approximations to the wave equation. *IMA J. Numer. Anal.*, 33(4):1245–1264, 2013.

[28] E. H. Georgoulis, O. Lakkis, Ch. Makridakis, and J. M. Virtanen. A posteriori error estimates for leap-frog and cosine methods for second order evolution problems. *SIAM J. Numer. Anal.*, 54(1):120–136, 2016.

[29] E. H. Georgoulis and Ch. G. Makridakis. Lower bounds, elliptic reconstruction and a posteriori error control of parabolic problems. *IMA J. Numer. Anal.*, 43(6):3212–3242, 2023.

[30] S. Gómez and V. Nikolič. Combined DG-CG finite element method for the Westervelt equation. https://arxiv.org/abs/2412.09095, 2024.

[31] O. Gorynina, A. Lozinski, and M. Picasso. Time and space adaptivity of the wave equation discretized in time by a second-order scheme. *IMA J. Numer. Anal.*, 39(4):1672–1705, 2019.

[32] P. Grisvard. *Elliptic Problems in Nonsmooth Domains*, volume 69 of *Classics in Applied Mathematics*. SIAM, 2011.

[33] M. Grote, O. Lakkis, and C. Santos. A posteriori error estimates for the wave equation with mesh change in the leapfrog method. http://arxiv.org/abs/2411.16933, 2024.

[34] M. J. Grote, A. Schneebeli, and D. Schötzau. Discontinuous Galerkin finite element method for the wave equation. *SIAM J. Numer. Anal.*, 44(6):2408–2431, 2006.

[35] B. Holm and Th. P. Wihler. Continuous and discontinuous Galerkin time stepping methods for nonlinear initial value problems with application to finite time blow-up. *Numer. Math.*, 138(3):767–799, 2018.

[36] Th. J. R. Hughes and G. M. Hulbert. Space-time finite element methods for elastodynamics: formulations and error estimates. *Comput. Methods Appl. Mech. Engrg.*, 66(3):339–363, 1988.

[37] G. M. Hulbert and Th. J. R. Hughes. Space-time finite element methods for second-order hyperbolic equations. *Comput. Methods Appl. Mech. Engrg.*, 84(3):327–348, 1990.

[38] A. Ilyin, A. Laptev, M. Loss, and S. Zelik. One-dimensional interpolation inequalities, Carlson–Landau inequalities, and magnetic Schrödinger operators. *Int. Math. Res. Not.* IMRN, 2016(4):1190–1222, 2016.

[39] L.-M. Imbert-Gérard, A. Moiola, and P. Stocker. A space-time quasi-Trefftz DG method for the wave equation with piecewise-smooth coefficients. *Math. Comp.*, 92(341):1211–1249, 2023.

[40] C. Johnson. Discontinuous Galerkin finite element methods for second order hyperbolic problems. *Comput. Methods Appl. Mech. Engrg.*, 107(1-2):117–129, 1993.

[41] O. Karakashian and Ch. Makridakis. Convergence of a continuous Galerkin method with mesh modification for nonlinear wave equations. *Math. Comp.*, 74(249):85–102, 2005.

[42] F. Kretzschmar, A. Moiola, I. Perugia, and S. M. Schnepp. *A priori error analysis of space-time Trefftz discontinuous Galerkin methods for wave problems*. *IMA J. Numer. Anal.*, 36(4):1599–1635, 2016.

[43] Ch. Makridakis. Space and time reconstructions in a posteriori analysis of evolution problems. In *ESAIM Proceedings*. Vol. 21, volume 21, pages 31–44. EDP Sci., Les Ulis, 2007.

[44] Ch. Makridakis and R. H. Nochetto. A posteriori error analysis for higher order dissipative methods for evolution problems. *Numer. Math.*, 104(4):489–514, 2006.

[45] A. Moiola and I. Perugia. A space-time Trefftz discontinuous Galerkin method for the acoustic wave equation in first-order formulation. *Numer. Math.*, 138(2):389–435, 2018.

[46] P.-A. Raviart and J.-M. Thomas. *Introduction à l'Analyse Numérique des Équations aux Dérivées Partielles*. Collection Mathématiques Appliquées pour la Maîtrise. Masson, Paris, 1983.

[47] D. Schötzau and Th. P. Wihler. A posteriori error estimation for hp-version time-stepping methods for parabolic partial differential equations. *Numer. Math.*, 115(3):475–509, 2010.
[48] Ch. Schwab. *p- and hp-Finite Element Methods: Theory and Applications in Solid and Fluid Mechanics*. The Clarendon Press, Oxford University Press, New York, 1998.

[49] O. Steinbach and M. Zank. A stabilized space-time finite element method for the wave equation. In *Advanced finite element methods with applications*, volume 128 of *Lect. Notes Comput. Sci. Eng.*, pages 341–370. Springer, Cham, 2019.

[50] O. Steinbach and M. Zank. Coercive space-time finite element methods for initial boundary value problems. *Electron. Trans. Numer. Anal.*, 52:154–194, 2020.

[51] L. N. Trefethen. *Approximation Theory and Approximation Practice, Extended Edition*. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2019.

[52] F. Verdugo and S. Badia. The software design of Gridap: a finite element package based on the Julia JIT compiler. *Comput. Phys. Commun.*, 276:Paper No. 108341, 24, 2022.

[53] R. Verfürth. *A Posteriori Error Estimation Techniques for Finite Element Methods*. Numerical Mathematics and Scientific Computation. Oxford University Press, Oxford, 2013.

[54] N. J. Walkington. Combined DG-CG time stepping for wave equations. *SIAM J. Numer. Anal.*, 52(3):1398–1417, 2014.

[55] D. Q. Yang. Grid modification for second-order hyperbolic problems. *Math. Comp.*, 64(212):1495–1509, 1995.