# A priori and a posteriori error estimates of a $C^0$-in-time method for the wave equation in second order formulation

Zhaonan Dong, Lorenzo Mascotto, Zuodong Wang

## ► To cite this version:

Zhaonan Dong, Lorenzo Mascotto, Zuodong Wang. A priori and a posteriori error estimates of a $C^0$-in-time method for the wave equation in second order formulation. 2025. hal-04768144v2

HAL Id: hal-04768144

https://inria.hal.science/hal-04768144v2

Preprint submitted on 27 Oct 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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 priori and a posteriori error estimates of a $C^0$-in-time method for the wave equation in second order formulation

Z. Dong,\* L. Mascotto,\d{z} Z. Wang\*

**Abstract**

We establish fully-discrete a priori and semi-discrete in time a posteriori error estimates for a discontinuous-continuous Galerkin discretization of the wave equation in second order formulation; the resulting method is a Petrov-Galerkin scheme based on piecewise and piecewise continuous polynomial in time test and trial spaces, respectively. Crucial tools in the a priori analysis for the fully-discrete formulation are the design of suitable projection and interpolation operators extending those used in the parabolic setting, and stability estimates based on a nonstandard choice of the test function; a priori estimates are shown, which are measured in $L^\infty$-type norms in time. For the semi-discrete in time formulation, we exhibit reliable a posteriori error estimates for the error measured in the $L^\infty(L^2)$ norm with fully explicit constants; to this aim, we design a reconstruction operator into $C^1$ piecewise polynomials over the time grid with optimal approximation properties in terms of the polynomial degree distribution and the time steps. Numerical examples illustrate the theoretical findings.

**AMS subject classification:** 65M50; 65M60; 65J10.

**Keywords:** wave equation; discontinuous Galerkin method; adaptive algorithm; hp-analysis; a posteriori error analysis.

# 1 Introduction

We establish fully-discrete a priori and semi-discrete in time a posteriori error estimates for a $C^0$-in-time method, see, e.g., [39], approximating solutions to the wave equation in second order formulation, which are explicit in the spatial mesh size, the time steps, and the polynomial degrees.

**Formulation, a priori estimates, and minimal literature.** The $C^0$-in-time method we are interested in is based on a second order formulation of the wave equation. Compared to several references where first order systems are considered [13,14,23,25], the dimension of the corresponding discrete spaces is smaller for fixed time steps.

The method lies in between a fully DG and a $C^1$ schemes in time: it employs piecewise and piecewise continuous polynomial in time for the test and trial spaces, respectively; upwind terms involving first time derivatives in time of the trial functions are included in the formulation. The polynomial degree in time of the trial functions is larger by 1 than that for the test functions, leading to square systems for each time interval, which are solved sequentially as a time marching scheme.

A key tool in deriving stability and a priori error estimates is the choice of an appropriate test function. For instance, in [39], a higher order fully-discrete version of the test function in [4] is used, leading to stability estimates in the $L^\infty$-type norms in time. Nonlinear problems are also discretized with this approach [18].

<footnote>\*Inria, 48 rue Barrault, 75013 Paris, France; CERMICS, ENPC, Institut Polytechnique de Paris, 77455 Marne-la-Vallee Cedex 2, France zhaonan.dong@inria.fr, zuodong.wang@inria.fr</footnote>

<footnote>†Dipartimento 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>
**A posteriori estimates and minimal literature.** A posteriori error estimates are well established for elliptic problems and a considerable amount of work is available for parabolic problems as well; on the other hand, hyperbolic (and in particular wave) problems are less explored.

A posteriori error estimates for wave problems in second order formulation are studied in [1] and later rigorously analysed in [5]; first order systems are instead the topic of [25, 27, 35]; a posteriori error estimates measured in the $L^\infty(L^2)$ norm are investigated for different time stepping schemes in [15, 16, 19]. A posteriori error estimates that are reliable and efficient in the same norm are instead investigated in [7, 8].

**Features of the $C^0$-in-time method.** Compared to fully DG schemes in time, the $C^0$-in-time method involves fewer unknowns; in view of the a posteriori error analysis for the semi-discrete in time version, one can employ tools from the parabolic setting and deduce a posteriori bounds that are fully explicit with respect to the polynomial degree. On the other hand, a modification of the scheme seems to be suited for dynamic mesh changes in space [9, 26].

**First main contribution of the manuscript.** We modify the analysis in [39] so as to have explicit dependence on the polynomial degree for the a priori analysis of the fully-discrete scheme. We consider uniform polynomial degree in the spatial discretization and possibly nonuniform polynomial degree in the time discretization. Static meshes in space are considered.

Our analysis hinges upon deriving stability estimates for the scheme, which are explicit in the polynomial degrees; see Section 2.1. Based on such stability estimates and the properties of an integrated Thomée operator discussed in Section 2.2, a priori error estimates are proven in Section 2.3. Error estimates are obtained in Section 2.4, which are explicit in the spatial mesh size, the time steps, and the polynomial degrees. Optimal estimates for the error measured with respect to $L^\infty$-type norms in time is shown for sufficiently regular solutions and data for fixed polynomial degrees in time and space.

**Second main contribution of the manuscript.** For the semi-discrete in time method, we design an error estimator satisfying a posteriori error estimates, with explicit dependence on the polynomial degree distribution in time. The crucial tool in a posteriori error estimates for time-dependent problems is the derivation of a reconstruction operator into smoother spaces. The original idea in the context of parabolic problems traces back to Makridakis and Nochetto [27]; the corresponding $p$-version analysis is detailed in [33] and later in [22]. We design a related operator in the wave equation setting and derive $p$-optimal approximation estimates in several norms in Section 3.2; such an operator is instrumental for designing an error estimator that is reliable for the error measured in the $L^\infty(L^2)$ norm; see Section 3.3; a posteriori error estimates in that norm are important as they are particularly suited to dynamic mesh changes as extensively explored in our recent work [9]. The upper bound is explicit in the polynomial degree distribution and the time steps, without unknown constants. Since the test and trial test functions have different polynomial degrees in time, the a posteriori error bounds involve extra oscillation terms compared to the parabolic setting.

To the best of our knowledge, we provide for the first time in the literature fully explicit, reliable a posteriori error estimates for a semi-discrete in time method for the approximation of solutions to the wave equation in second order formulation, which are explicit in the polynomial degree distribution in time and optimal in the time steps; implicit time-stepping schemes are essential in view of future analysis of coupled wave and parabolic problems, e.g. fluid-structure interaction problems. The proposed analysis does not hinge upon any CFL condition, which is greatly advantageous for adaptivity whilst compared to methods based on explicit time stepping; there, for each mesh refinement, one has to check whether the resulting spatial mesh size is sufficiently small compared to the corresponding time step and the polynomial degree [21]. This improvement is even more effective for wave problems involving higher order spatial elliptic operators [30], where the CFL condition poses even stricter constraints on the spatial mesh size.

**List of the main results and advances.** For the reader's convenience, we list here the main results of the manuscript ([APRI] = a priori analysis; [APOS] = a posteriori analysis):
* [APRI] Theorem 2.1 is concerned with fully explicit stability estimates for the $C^0$-in-time formulation, which are explicit with respect to the polynomial degrees in time and space employed throughout;

* [APRI] Proposition 2.5 analyzes the approximation properties of a novel integrated Thomée operator;

* [APRI] Theorem 2.8 discusses a priori estimates that are explicit with respect to the spatial mesh size, the time step distribution, the spatial polynomial degree, and the polynomial degree distribution in time;

* [APOS] we define an error estimator $\eta$ in (70) for the semi-discrete in time formulation;

* [APOS] corresponding reliability estimates with respect to the $L^\infty(L^2)$ norm of the error are displayed in Proposition 3.6.

**Functional setting.** Standard notation is used throughout for Sobolev and Bochner spaces. Let $D$ be a Lipschitz domain in $\mathbb{R}^d$, $d = 1, 2$, and $3$, with boundary $\partial D$. The space of Lebesgue measurable and square integrable functions over $D$ is $L^2(D)$. The Sobolev space of positive integer order $s$ is $H^s(D)$. We endow $H^s(D)$ with the inner product, seminorm, and norm

$$ (\cdot, \cdot)_{s,D}, \quad |\cdot|_{s,D}, \quad \|\cdot\|_{s,D}. $$

Interpolation theory is used to construct Sobolev spaces of positive noninteger order; duality is used to define negative order Sobolev spaces. We shall be particularly using the space $H^{-1}(D)$, which is the dual of $H_0^1(D)$; the duality pairing between the two spaces is $\langle \cdot, \cdot \rangle$. The space of polynomials of nonnegative degree $p$ over $D$ is $\mathbb{P}_p(D)$.

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

$$ \|v\|_{L^t(I; \mathcal{X})} := \begin{cases} \left( \int_I \|v\|_{\mathcal{X}}^t \right)^{\frac{1}{t}} dt & \text{for } 1 \le t < \infty \\ \ess \sup_{t \in I} \|v\|_{\mathcal{X}} & \text{for } t = \infty. \end{cases} $$

For $s$ in $\mathbb{N}$, the space $H^s(I; \mathcal{X})$ is the space of measurable functions $v$ whose derivatives in time up to order $s$ belong to $L^2(I; \mathcal{X})$. For any real number $s$ larger than or equal to 0, the space $H^s(I; \mathcal{X})$ is constructed using interpolation of integer order Bochner spaces. Bochner inner products are denoted by $(\cdot, \cdot)_{L^t(I; \mathcal{X})}$ and $(\cdot, \cdot)_{H^s(I; \mathcal{X})}$.

To avoid confusion, the seminorm symbol $|\cdot|$ is only used to denote spatial seminorms. Seminorms in time are rather displayed as $L^2$ norms of a suitable time derivative. The first and second derivative symbols are $'$ and $''$; time derivatives of order $s$ larger than 2 are displayed as $^{(s)}$.

**The continuous problem.** Let $\Omega$ be a polytopic, Lipschitz domain in $\mathbb{R}^d$, $d = 1, 2, 3$; $T$ a positive final time; $Q_T := (0, T] \times \Omega$ the space-time cylinder; $u_0$ in $H_0^1(\Omega)$; $u_1$ in $L^2(\Omega)$; $f$ in $L^2(0, T; L^2(\Omega))$.

Given $\Delta_x$ and $\nabla_x$ the spatial Laplace and gradient operator, we consider the following problem: find $u: Q_T \to \mathbb{R}$ such that

$$ \left\{ \begin{array}{ll} u'' - \Delta_x u = f & \text{in } Q_T \\ u(t, \cdot) = 0 & \text{on } (0, T] \times \partial \Omega \\ u(0, \cdot) = u_0(\cdot) & \text{on } \Omega \\ u'(0, \cdot) = u_1(\cdot) & \text{on } \Omega. \end{array} \right. \qquad (1) $$

Introduce the spaces

$$ X := H^2(0, T; H^{-1}(\Omega)) \cap L^2(0, T; H_0^1(\Omega)) \cap H^1(0, T; L^2(\Omega)), \quad Y := L^2(0, T; H_0^1(\Omega)), $$

and the bilinear form on $H_0^1(\Omega) \times H_0^1(\Omega)$

$$ a(u, v) := (\nabla_x u, \nabla_x v)_{0,\Omega}. $$
We consider the following weak formulation of problem (1):

$$ egin{cases} 
\text{find } u \in X \text{ such that} \\
\int_0^T [\langle u'', v \rangle + a(u, v)] dt = \int_0^T (f, v)_{0,\Omega} dt & \forall v \in Y \\
u(0, \cdot) = u_0(\cdot) \text{ in } H_0^1(\Omega), & u'(0, \cdot) = u_1(\cdot) \text{ in } L^2(\Omega). 
\end{cases} 
\quad (2) $$

Problem (1) is well posed; see, e.g., [31, Chapter 8].

In [39], inhomogeneous Dirichlet and inhomogeneous Neumann boundary conditions are considered; This results in further complication in the analysis below, cf. [39, Example 3.2]. For this reason, we prefer to stick to the setting in (2).

**Spatial meshes, time grids, and polynomial degree distributions**. We consider either a simplicial or tensor-product conforming mesh $\mathcal{T}_h$ of $\Omega$ and a corresponding $H^1$-conforming Lagrangian finite element space $V_h$ of uniform order $p_x$. We assume the existence of a constant $\gamma$ in $(0, 1)$ such that each element of $\mathcal{T}_h$ is star-shaped with respect to a ball of radius larger than or equal to the diameter of that element; moreover, we assume quasi-uniformity of the mesh, i.e., given $h_1$ and $h_2$ the diameters of two arbitrary elements with $h_1 \le h_2$, one has $h_2 \le \gamma h_1$. $^1$.

We further consider a decomposition $0 = t_0 < t_1 < \cdots < t_N = T$ of $[0, T]$ and introduce $\tau_n := t_n - t_{n-1}$ for all $n = 1, \dots, N$. With each time interval $I_n := (t_{n-1}, t_n]$, we associate a local polynomial degree $p_n^t$; we collect such polynomial orders (in time) in the vector $\mathbf{p}^t$ in $\mathbb{N}_+^N$, and set $\mathbf{p}_n^t := p_n^t$. For $k \in \mathbb{Z}$, $\mathbf{p}^t + k$ is the vector of entries $p_n^t + k$.

For all $n = 1, \dots, N$, we set $(v')^+(t_{n-1}, \cdot) := v'_{|_{I_n}}(t_{n-1}, \cdot)$ and $(v')^-(t_{n-1}, \cdot) := v'_{|_{I_{n-1}}}(t_{n-1}, \cdot)$; for all $v$ piecewise continuous in time, we set $v^+(t_{n-1}, \cdot) := v_{|_{I_n}}(t_{n-1}, \cdot)$ and $v^-(t_{n-1}, \cdot) := v_{|_{I_{n-1}}}(t_{n-1}, \cdot)$.

We also define the tensor product space

$$ \mathbb{P}_{p_n^t}(I_n; V_h) := \{W_h \in L^2(I_n; V_h) \mid W_h = a(x)b(t), a \in V_h, b \in \mathbb{P}_{p_n^t}(I_n)\}. $$

Throughout we assume that

$$ p_n^t \ge 2 \qquad \forall n = 1, \dots, N. $$

**The fully-discrete $C^0$-in-time method.** Let $u_{0,h}$ and $u_{1,h}$ be approximations of $u_0$ and $u_1$ in $V_h$. Throughout, we pick $u_{0,h}$ as the elliptic projection of $u_0$ defined in display (21) below; $u_{1,h}$ as the $L^2$ projection of $u_1$ onto $V_h$. Other variants are possible, but are omitted here.

We define

$$ X_{h, \tau} := \{ U_h \in C^0(0, T; V_h) \mid U_h(0, \cdot) = u_{0,h}, U_h|_{I_n} \in \mathbb{P}_{p_n^t}(I_n; V_h) \forall n = 1, \dots, N \} $$

and the upwind jump operator for the time derivative as

$$ [[U_h']](t_{n-1}, \cdot) = \begin{cases} U_h'_{|I_1}(0, \cdot) - u_{1,h}(\cdot) & \text{if } n = 1 \\ U_h'_{|I_n}(t_{n-1}, \cdot) - U_h'_{|I_{n-1}}(t_{n-1}, \cdot) & \text{if } n = 2, \dots, N. \end{cases} $$

The $C^0$-in-time method, see, e.g., [39], reads as follows: find $U_h$ in $X_{h, \tau}$ such that

$$ \int_{I_n} [(U_h'', W_h)_{0,\Omega} + a(U_h, W_h)] dt + ([[U_h']](t_{n-1}, \cdot), W_h^+(t_{n-1}, \cdot))_{0,\Omega} \\ = (f, W_h)_{L^2(I_n; L^2(\Omega))} \qquad \forall W_h \in \mathbb{P}_{p_n^t-1}(I_n; V_h), \quad \forall n = 1, \dots, N. \quad (3) $$

The initial condition $u_{0,h}$ is imposed strongly in $X_{h, \tau}$; the initial condition $u_{1,h}$ is imposed weakly through the upwind term at the initial time.

Method (3) is solved time-slab by time-slab as a time marching scheme. Initial conditions on each time slab are assigned taking the values of the solution at the final time of the previous slab and upwinding the first time derivative.

<footnote n="1">¹In the a priori analysis, the only place where quasi-uniformity is used in lieu of local quasi-uniformity is while deriving estimates of the elliptic projector in the $L^2$ norm; see displays (21) and (22), and Remark 2 below</footnote>
The existence and uniqueness of a solution and the data of method (3) follow, e.g., assuming sufficient smoothness of the solution to problem (2), showing stability estimates as those in Theorem 2.1 below (which imply uniqueness), and using the fact that on each time slab the linear system to solve is square (which entails that existence is equivalent to uniqueness).

**The semi-discrete in time $C^0$-in-time method.** Define the space

$$X_{\tau} := \{U \in C^0(0, T; H_0^1(\Omega)) \mid U(0, \cdot) = u_0, U|_{I_n} \in \mathbb{P}_{p_n^t}(I_n; H_0^1(\Omega))\}.$$

In Section 3 below, we prove fully explicit, reliable a posteriori error estimates for the time semi-discrete in time version of (3). More precisely, we look for $U$ in $X_{\tau}$ such that

$$\int_{I_n} [(U'', V)_{0,\Omega} + a(U, V)]dt + (\llbracket U' \rrbracket (t_{n-1}, \cdot), V^+(t_{n-1}, \cdot))_{0,\Omega} = (f, V)_{L^2(0,T;L^2(\Omega))} \quad (4) \\ \forall V \in \mathbb{P}_{p_n^t-1}(I_n; H_0^1(\Omega)), \quad \forall n = 1, \dots, N,$$

with $u_0$ imposed strongly in $X_{\tau}$ and $u_1$ imposed weakly through upwinding.

Let $\Pi_{p_n^t-1}^0$ denote the piecewise $L^2$ projector onto $\mathbb{P}_{p_n^t-1}(I_n; L^2(\Omega))$ for all $n = 1, \dots, N$. We can replace the right-hand side of (4) with the following expression:

$$(f, V)_{L^2(I_n; L^2(\Omega))} = (\Pi_{p_n^t-1}^0 f, V)_{L^2(I_n; L^2(\Omega))} \quad \forall n = 1, \dots, N. \quad (5)$$

The well-posedness of (4) follows from standard arguments of DG time-stepping schemes; cf. [36, Ch. 12].

**Structure of the remainder of the paper.** We discuss stability and error estimates of the fully-discrete method (3) in Section 2, which are explicit in the spatial mesh size, the time steps, and the polynomial degrees. A posteriori error estimates for the semi-discrete in time method (4) are derived in Section 3. We assess the numerical findings with numerical experiments in Section 4, and draw some conclusions in Section 5.

# 2 A priori error analysis

This section is concerned with proving stability and a priori estimates for method (3): in Section 2.1, we show stability estimates following the analysis in [39] by tracking the explicit dependence on the polynomial degree distributions; in Section 2.2, we introduce and show the approximation properties of an integrated Thomée-type operator in terms of the spatial mesh size, the time steps, and the polynomial degrees, which are instrumental in deriving abstract error estimates in Section 2.3; standard polynomial approximation results yield error estimates in Section 2.4.

## 2.1 Stability estimates

The stability of method (3) in certain norms is investigated in [39, Theorem 4.5]. The main idea behind the derivation of stability estimates is to take a suitable test function, namely the $L^2$ projection onto the correct test space of the time derivative of the discrete solution times a weight mimicking an exponential function; this idea traces back to [13]; cf. [7, 8] for more recent similar results. The reason for this is that testing only with the time derivative of the discrete solution would yield to stability estimates at the time grid point only, i.e., no global stability estimates would be available; see [39, eq. (4.1)].

The constants in the stability estimates in [39], depend implicitly on the distribution $p^t$ of polynomial degrees in time; for this reason, we revisit that proof so as to carry out an explicit analysis in terms of the polynomial distribution in time.

To this aim, given a generic element $g$ in the space $X_{h,\tau}$, let $m = m(g)$ so that $I_m$ is the interval where

$$\|g'\|_{L^\infty(I_m; L^2(\Omega))}^2 + |g|_{L^\infty(I_m; H^1(\Omega))}^2 = \max_{n=1}^N (\|g'\|_{L^\infty(I_n; L^2(\Omega))}^2 + |g|_{L^\infty(I_n; H^1(\Omega))}^2), \quad (6)$$
The index $m$ depends on the choice of $g$ and can be expected on most occasions to be attained at the final interval. We also define

$$ \mu_n := \frac{1}{1024(p_n^t)^2(2p_n^t + 1)} \quad \forall n = 1, \dots, N. \quad (7) $$

**Theorem 2.1.** Let $U_h$ be the solution to (3) and $f$ be the source term in (1). The following stability estimate holds true<fnref n="2" />:

$$ \mu_m \left( \|U'_h\|_{L^\infty(I_m; L^2(\Omega))}^2 + |U_h|_{L^\infty(I_m; H^1(\Omega))}^2 \right) + \frac{1}{4} \sum_{n=1}^{m} \|[[U'_h]](t_{n-1}, \cdot)\|_{0,\Omega}^2 \\ \le \frac{1}{2} (|u_0|_{1,\Omega}^2 + \|u_1\|_{0,\Omega}^2) + \frac{t_m}{\mu_m} \|f\|_{L^2(0, t_m; L^2(\Omega))}^2. \quad (8) $$

*Proof.* Define

$$ \lambda_n := \frac{1}{4(2p_n^t + 1)\tau_n} \quad \forall n = 1, \dots, N. $$

Restrict the (piecewise in time) $L^2$ projector $\Pi_{\mathbf{p}^{t-1}}^0$ to $\mathbb{P}_{p_n^{t-1}}(I_n; V_h)$ and consider the following test function already used in [39, Theorem 4.5]:

$$ W_{h|I_n} := \Pi_{\mathbf{p}^{t-1}}^0([1 - \lambda_n(t - t_{n-1})]U'_h) \quad \forall n = 1, \dots, N. $$

Picking $W_h$ as above in (3) and proceeding as in the proof of [39, Theorem 4.5], more precisely see [39, eq. (4.2)], yield, for all $n = 1, \dots, N$,

$$ (1 - \lambda_n \tau_n) \frac{1}{2} (|U_h(t_n, \cdot)|_{1,\Omega}^2 + \|(U_h^-)'(t_n, \cdot)\|_{0,\Omega}^2) + \frac{\lambda_n}{2} \left( \|U'_h\|_{L^2(I_n; L^2(\Omega))}^2 + |U_h|_{L^2(I_n; H^1(\Omega))}^2 \right) \\ + \frac{1}{2} \|[[U'_h]](t_{n-1}, \cdot)\|_{0,\Omega}^2 \le \frac{1}{2} (|U_h(t_{n-1}, \cdot)|_{1,\Omega}^2 + \|(U_h^-)'(t_{n-1}, \cdot)\|_{0,\Omega}^2) \\ + ([[[U'_h]](t_{n-1}, \cdot), (U_h^+)'(t_{n-1}, \cdot) - W_h^+(t_{n-1}, \cdot))_{0,\Omega} + \int_{I_n} [1 - \lambda_n(\cdot - t_{n-1})](\Pi_{\mathbf{p}^{t-1}}^0 f, U'_h)_{0,\Omega} dt. $$

We estimate the last two terms on the right-hand side separately: one involving the jump of the first derivative at $t_{n-1}$; the other involving the source term $f$.

As for the “jump” term, we invoke [39, Corollary 4.4] and the definition of $\lambda_n$, use standard manipulations, and get

$$ \begin{aligned} & ([[U'_h]](t_{n-1}, \cdot), (U_h^+)'(t_{n-1}, \cdot) - W_h^+(t_{n-1}, \cdot))_{0,\Omega} \\ & \le \|[[U'_h]](t_{n-1}, \cdot)\|_{0,\Omega} \|(U_h^+)'(t_{n-1}, \cdot) - W_h^+(t_{n-1}, \cdot)\|_{0,\Omega} \\ & \le \|[[U'_h]](t_{n-1}, \cdot)\|_{0,\Omega} \lambda_n \sqrt{(2p_n^t + 1)\tau_n} \|U'_h\|_{L^2(I_n; L^2(\Omega))} \\ & \le \sqrt{(2p_n^t + 1)\tau_n \lambda_n} \left( \frac{1}{2} \|[[U'_h]](t_{n-1}, \cdot)\|_{0,\Omega}^2 + \frac{\lambda_n}{2} \|U'_h\|_{L^2(I_n; L^2(\Omega))}^2 \right) \\ & = \frac{1}{2} \left( \frac{1}{2} \|[[U'_h]](t_{n-1}, \cdot)\|_{0,\Omega}^2 + \frac{\lambda_n}{2} \|U'_h\|_{L^2(I_n; L^2(\Omega))}^2 \right) \quad \forall n = 1, \dots, N. \end{aligned} $$

As for the “source” term, we write

$$ \int_{I_n} [1 - \lambda_n(\cdot - t_{n-1})](\Pi_{\mathbf{p}^{t-1}}^0 f, U'_h)_{0,\Omega} dt \le \|\Pi_{\mathbf{p}^{t-1}}^0 f\|_{L^1(I_n; L^2(\Omega))} \|U'_h\|_{L^\infty(I_n; L^2(\Omega))}. $$

We combine the three displays above: for all $n = 1, \dots, N$,

$$ \begin{aligned} & (1 - \lambda_n \tau_n) \frac{1}{2} (|U_h(t_n, \cdot)|_{1,\Omega}^2 + \|(U_h^-)'(t_n, \cdot)\|_{0,\Omega}^2) \\ & + \frac{\lambda_n}{4} \left( \|U'_h\|_{L^2(I_n; L^2(\Omega))}^2 + |U_h|_{L^2(I_n; H^1(\Omega))}^2 \right) + \frac{1}{4} \|[[U'_h]](t_{n-1}, \cdot)\|_{0,\Omega}^2 \\ & \le \frac{1}{2} (|U_h(t_{n-1}, \cdot)|_{1,\Omega}^2 + \|(U_h^-)'(t_{n-1}, \cdot)\|_{0,\Omega}^2) + \|\Pi_{\mathbf{p}^{t-1}}^0 f\|_{L^1(I_n; L^2(\Omega))} \|U'_h\|_{L^\infty(I_n; L^2(\Omega))}. \end{aligned} $$

<footnote n="2">²The norms of $U'_h$ increase cubically in $p$ with respect to the norm of the initial conditions.</footnote>
Using that $\lambda_n \tau_n < 1/4$, and $1 - \lambda_n \tau_n \ge 0$, we simplify this inequality:

$$ 
\begin{aligned}
& \frac{\lambda_n}{4} \left( \|U'_h\|_{L^2(I_n; L^2(\Omega))}^2 + |U_h|_{L^2(I_n; H^1(\Omega))}^2 \right) + \frac{1}{4} \|[[U'_h]](t_{n-1}, \cdot)\|_{0,\Omega}^2 \\
& \le \frac{1}{2} \left( |U_h(t_{n-1}, \cdot)|_{1,\Omega}^2 + \|(U_h^-)'(t_{n-1}, \cdot)\|_{0,\Omega}^2 \right) + \|\Pi_{\mathbf{p}'}^0 (-1) f\|_{L^1(I_n; L^2(\Omega))} \|U'_h\|_{L^\infty(I_n; L^2(\Omega))}.
\end{aligned}
\quad (9) $$

In [39, eq. (4.1)], a stability estimate is proven taking $U'_h$ as a test function, which gives, for all $n = 2, \dots, N$,

$$ 
\begin{aligned}
& \frac{1}{2} \left( |U_h(t_{n-1}, \cdot)|_{1,\Omega}^2 + \|(U_h^-)'(t_{n-1}, \cdot)\|_{0,\Omega}^2 \right) + \frac{1}{4} \sum_{\ell=1}^{n-1} \|[[U'_h]](t_{\ell-1}, \cdot)\|_{0,\Omega}^2 \\
& \le \frac{1}{2} \left( |U_h(0, \cdot)|_{1,\Omega}^2 + \|(U_h^-)'(0, \cdot)\|_{0,\Omega}^2 \right) \\
& \qquad + \|\Pi_{\mathbf{p}'}^0 (-1) f\|_{L^1(0, t_{n-1}; L^2(\Omega))} \left( \frac{1}{2} \|U'_h\|_{L^\infty(0, t_{n-1}; L^2(\Omega))}^2 + \frac{1}{2} |U_h|_{L^\infty(0, t_{n-1}; H^1(\Omega))}^2 \right)^{\frac{1}{2}}.
\end{aligned}
$$

Inserting the bound

$$ \|U'_h\|_{L^\infty(I_n; L^2(\Omega))} \le \left( \|U'_h\|_{L^\infty(I_n; L^2(\Omega))}^2 + |U_h|_{L^\infty(I_n; H^1(\Omega))}^2 \right)^{\frac{1}{2}} \quad \forall n = 1, \dots, N $$

in (9), selecting $n$ to be $m = m(U_h)$ as in (6), and combining the three displays above give

$$ 
\begin{aligned}
& \frac{\lambda_m}{4} \left( \|U'_h\|_{L^2(I_m; L^2(\Omega))}^2 + |U_h|_{L^2(I_m; H^1(\Omega))}^2 \right) + \frac{1}{4} \sum_{n=1}^{m} \|[[U'_h]](t_{n-1}, \cdot)\|_{0,\Omega}^2 \\
& \le \frac{1}{2} \left( |U_h(0, \cdot)|_{1,\Omega}^2 + \|(U_h^-)'(0, \cdot)\|_{0,\Omega}^2 \right) \\
& \qquad + 2 \left( \sum_{n=1}^{m} \|\Pi_{\mathbf{p}'}^0 (-1) f\|_{L^1(I_m; L^2(\Omega))} \right) \left( \|U'_h\|_{L^\infty(I_m; L^2(\Omega))}^2 + |U_h|_{L^\infty(I_m; H^1(\Omega))}^2 \right)^{\frac{1}{2}}.
\end{aligned}
$$

We recall the one dimensional $L^\infty$ to $L^2$-norm polynomial inverse inequality in [34, eq. (3.6.4)]:

$$ |U_h|_{L^\infty(I_m; H^1(\Omega))}^2 \le \frac{32(p_m^t)^2}{\tau_n} |U_h|_{L^2(I_m; H^1(\Omega))}^2. $$

Using that $2\mu_m = 1/(512(p_m^t)^2(2p_m^t + 1))$, we deduce

$$ 
\begin{aligned}
& \frac{1}{512(p_m^t)^2(2p_m^t + 1)} \left( \|U'_h\|_{L^\infty(I_m; L^2(\Omega))}^2 + |U_h|_{L^\infty(I_m; H^1(\Omega))}^2 \right) + \frac{1}{4} \sum_{n=1}^{m} \|[[U'_h]](t_{n-1}, \cdot)\|_{0,\Omega}^2 \\
& \le \frac{1}{2} \left( |U_h(0, \cdot)|_{1,\Omega}^2 + \|(U_h^-)'(0, \cdot)\|_{0,\Omega}^2 \right) \\
& \qquad + 2 \left( \sum_{n=1}^{m} \|\Pi_{\mathbf{p}'}^0 (-1) f\|_{L^1(I_m; L^2(\Omega))} \right) \left( \|U'_h\|_{L^\infty(I_m; L^2(\Omega))}^2 + |U_h|_{L^\infty(I_m; H^1(\Omega))}^2 \right)^{\frac{1}{2}}.
\end{aligned}
$$

Using Young's inequality $2ab \le \varepsilon a^2 + \frac{b^2}{\varepsilon}$, with $\varepsilon = 1024(p_m^t)^2(2p_m^t + 1) =: \mu_m^{-1}$ for the last term in the above relation, we infer

$$ 
\begin{aligned}
& \mu_m \left( \|U'_h\|_{L^\infty(I_m; L^2(\Omega))}^2 + |U_h|_{L^\infty(I_m; H^1(\Omega))}^2 \right) + \frac{1}{4} \sum_{n=1}^{m} \|[[U'_h]](t_{n-1}, \cdot)\|_{0,\Omega}^2 \\
& \le \frac{1}{2} \left( |U_h(0, \cdot)|_{1,\Omega}^2 + \|(U_h^-)'(0, \cdot)\|_{0,\Omega}^2 \right) + \frac{1}{\mu_m} \left( \sum_{n=1}^{m} \|\Pi_{\mathbf{p}'}^0 (-1) f\|_{L^1(I_m; L^2(\Omega))} \right)^2.
\end{aligned}
$$

We use the Hölder inequality, the stability of the $L^2$ projector in the $L^2(I_n)$ norm, and the $\ell^2$ Cauchy-Schwarz inequality in the last term on the right-hand side above:

$$ \sum_{n=1}^{m} \|\Pi_{\mathbf{p}'}^0 (-1) f\|_{L^1(I_n; L^2(\Omega))} \le \sum_{n=1}^{m} \tau_n^{\frac{1}{2}} \|f\|_{L^2(I_n; L^2(\Omega))} \le t_m^{\frac{1}{2}} \|f\|_{L^2(0, t_m; L^2(\Omega))}. \quad (10) $$

The assertion follows by using the stability of the discrete initial conditions $u_{0,h}$ and $u_{1,h}$ with respect to the $H_0^1(\Omega)$ and $L^2(\Omega)$ norms. $\square$
Remark 1. The norms appearing on the left-hand side of (8) are of the same sort of those in [39, Theorem 4.5]. Instead, that on the right-hand side involving the source term $f$ differs a bit: here, we employ an $L^2(L^2)$-type norm weighted with $t_m^{\frac{1}{2}}/\mu_m^{\frac{1}{2}}$, which scales exactly as the $L^1(L^2)$ norm used by Walkington. The reason for this change resides in inequality (10); by proceeding as in [39], we would end up with further suboptimality by two orders in the polynomial degree due to the $L^1$ stability of the $L^2$ projector.

## 2.2 An integrated Thomée-type operator

Given a Hilbert space $\mathcal{X}$ with inner product $(\cdot, \cdot)_{\mathcal{X}}$ and induced norm $\|\cdot\|_{\mathcal{X}}$, we introduce an operator $\mathcal{P}_{\mathbf{p}^t}$ mapping $\mathcal{C}^1(0, T; \mathcal{X})$ into the space of piecewise polynomials with polynomial distribution $\mathbf{p}^t$ over the time grid. This operator is the integrated version of an operator in [36] in the framework of parabolic problems, see also [2], and whose $p$-approximation properties are derived in [32].

The operator $\mathcal{P}_{\mathbf{p}^t}$ is piecewise defined fixing the following conditions: for all $n = 1, \dots, N$,

$$ \begin{cases} (w' - \mathcal{P}_{\mathbf{p}^t}(w)', q_{p_n^t-2}^{L^2(I_n; \mathcal{X})}) = 0 & \forall n = 1, \dots, N, \quad \forall q_{p_n^t-2} \in \mathbb{P}_{p_n^t-2}^{L^2(I_n; \mathcal{X})}; \\ \mathcal{P}_{\mathbf{p}^t}(w)'(t_n^-, \cdot) = w'(t_n, \cdot); & \mathcal{P}_{\mathbf{p}^t}(w)(t_{n-1}, \cdot) = w(t_{n-1}, \cdot) & \text{in } \mathcal{X}. \end{cases} $$

As shown in [39, Section 5.1], this operator is well defined and the above definition is equivalent to

$$ \begin{cases} (w - \mathcal{P}_{\mathbf{p}^t}(w), q_{p_n^t-3}^{L^2(I_n; \mathcal{X})}) = 0 & \forall n = 1, \dots, N, \quad \forall q_{p_n^t-3} \in \mathbb{P}_{p_n^t-3}^{L^2(I_n; \mathcal{X})}, \\ \mathcal{P}_{\mathbf{p}^t}(w)'(t_n^-, \cdot) = w'(t_n, \cdot), & \mathcal{P}_{\mathbf{p}^t}(w)(t_n, \cdot) = w(t_n, \cdot), & \mathcal{P}_{\mathbf{p}^t}(w)(t_{n-1}, \cdot) = w(t_{n-1}, \cdot) & \text{in } \mathcal{X}, \end{cases} \quad (11) $$

where the orthogonality condition in the first line of (11) is not taken into account if $p_n^t = 2$.

The function $\mathcal{P}_{\mathbf{p}^t}(w)$ is globally continuous but not necessarily $\mathcal{C}^1$ in time and satisfies the following property.

**Lemma 2.2.** For all $u$ in $H^2(I_n; \mathcal{X})$ and all $q_{p_n^t-1}$ in $\mathbb{P}_{p_n^t-1}(I_n; \mathcal{X})$, the following identity holds true:

$$ ((u - \mathcal{P}_{\mathbf{p}^t} u)'', q_{p_n^t-1}^{L^2(I_n; \mathcal{X})}) = (\| (\mathcal{P}_{\mathbf{p}^t} u)' \| (t_{n-1}, \cdot), q_{p_n^t-1}^{L^2(I_n; \mathcal{X})}) \quad \forall n = 1, \dots, N. \quad (12) $$

*Proof.* We have

$$ \begin{align*} & ((u - \mathcal{P}_{\mathbf{p}^t} u)'', q_{p_n^t-1}^{L^2(I_n; \mathcal{X})}) \\ &= -((u - \mathcal{P}_{\mathbf{p}^t} u)', q_{p_n^t-1}^{L^2(I_n; \mathcal{X})}) + ((u - \mathcal{P}_{\mathbf{p}^t} u)'(t_n, \cdot), q_{p_n^t-1}^{L^2(I_n; \mathcal{X})}) \chi \\ &\quad - ((u - \mathcal{P}_{\mathbf{p}^t} u)'(t_{n-1}, \cdot), q_{p_n^t-1}^{L^2(I_n; \mathcal{X})}) \chi \\ &\stackrel{(11)}{=} -((u - \mathcal{P}_{\mathbf{p}^t} u)'(t_{n-1}, \cdot), q_{p_n^t-1}^{L^2(I_n; \mathcal{X})}) \chi \stackrel{(11)}{=} (\| (\mathcal{P}_{\mathbf{p}^t} u)' \| (t_{n-1}, \cdot), q_{p_n^t-1}^{L^2(I_n; \mathcal{X})}) \chi, \end{align*} $$

which is the assertion. $\square$

The scope of this section is showing the approximation properties of $\mathcal{P}_{\mathbf{p}^t}$. We recall the properties of the original Thomée operator [36]. Introduce $\tilde{\mathcal{P}}_{\mathbf{p}^t}$ mapping $\mathcal{C}^0(0, T; \mathcal{X})$ into the space of piecewise polynomials with degree distribution $\mathbf{p}^t$ over the time grid (in particular, discontinuous functions in time). The operator $\tilde{\mathcal{P}}_{\mathbf{p}^t}$ is defined fixing the following conditions:

$$ \begin{cases} (w - \tilde{\mathcal{P}}_{\mathbf{p}^t}(w), q_{p_n^t-1}^{L^2(I_n; \mathcal{X})}) = 0 & \forall n = 1, \dots, N, \quad \forall q_{p_n^t-1} \in \mathbb{P}_{p_n^t-1}^{L^2(I_n; \mathcal{X})}; \\ \tilde{\mathcal{P}}_{\mathbf{p}^t}(w)(t_n^-, \cdot) = w(t_n, \cdot) & \text{in } \mathcal{X}. \end{cases} \quad (13) $$

The following approximation results can be found in [32, Lemmas 3.6–3.8, Theorem 3.10] and [6, Lemma 3.3].
**Lemma 2.3.** Let $\tilde{\mathcal{P}}_{\mathbf{p}^t}$ be the operator in (13). Then, for all $w$ in $H^{s+1}(I_n; \mathcal{X})$ with $s \ge 0$, the following inequalities hold true for all $n = 1, \dots, N$:

$$ \|w - \tilde{\mathcal{P}}_{\mathbf{p}^t} w\|_{L^2(I_n; \mathcal{X})} \lesssim \|w - \Pi_{\mathbf{p}^t}^0 w\|_{L^2(I_n; \mathcal{X})} + \frac{\tau_n}{p_n^t} \|w'\|_{L^2(I_n; \mathcal{X})}, \quad (14a) $$

$$ \begin{aligned} \|w - \tilde{\mathcal{P}}_{\mathbf{p}^t} w\|_{L^2(I_n; \mathcal{X})} & \lesssim \inf_{q_{p_n^t} \in \mathbb{P}_{p_n^t}(I_n; \mathcal{X})} \left( \|w - q_{p_n^t}\|_{L^2(I_n; \mathcal{X})} + \frac{\tau_n}{p_n^t} \|(w - q_{p_n^t})'\|_{L^2(I_n; \mathcal{X})} \right) \\ & \lesssim \left( \frac{\tau_n}{p_n^t} \right)^{\min(s, p_n^t) + 1} \|w^{(s+1)}\|_{L^2(I_n; \mathcal{X})}, \end{aligned} \quad (14b) $$

$$ \|\tilde{\mathcal{P}}_{\mathbf{p}^t} w\|_{L^2(I_n; \mathcal{X})} \lesssim \|w\|_{L^2(I_n; \mathcal{X})} + \frac{\tau_n}{p_n^t} \|w'\|_{L^2(I_n; \mathcal{X})}, \quad (14c) $$

$$ \begin{aligned} \|(w - \tilde{\mathcal{P}}_{\mathbf{p}^t} w)(t_{n-1}, \cdot)\|_{\mathcal{X}} & \lesssim \left(\frac{\tau_n}{p_n^t}\right)^{\frac{1}{2}} \inf_{q_{p_n^t} \in \mathbb{P}_{p_n^t}(I_n; \mathcal{X})} \|(w - q_{p_n^t})'\|_{L^2(I_n; \mathcal{X})} \\ & \lesssim \left(\frac{\tau_n}{p_n^t}\right)^{\min(s, p_n^t) + \frac{1}{2}} \|w^{(s+1)}\|_{L^2(I_n; \mathcal{X})}. \end{aligned} \quad (14d) $$

On the reference interval $\hat{I} = [-1, 1]^3$, the bounds in Lemma 2.3 are a consequence of the properties of the Legendre polynomials, the identity

$$ (w - \tilde{\mathcal{P}}_{\mathbf{p}^t} w)(t, \cdot) = (w - \Pi_{\mathbf{p}^t}^0 w)(t, \cdot) + (w - \Pi_{\mathbf{p}^t}^0 w)(1, \cdot) L_p(t) \quad \text{in } \mathcal{X}, \quad (15) $$

and sharp estimates [32, Lemmas 3.5 and 3.6] for the second term on the right-hand side of (15).

The operators $\mathcal{P}_{\mathbf{p}^t}$ and $\tilde{\mathcal{P}}_{\mathbf{p}^t}$ in (11) and (13) are related by the following identity: for all $t$ in $I_n$,

$$ \mathcal{P}_{\mathbf{p}^t}(w)(t, \cdot) = \int_{t_{n-1}}^{t} (\tilde{\mathcal{P}}_{\mathbf{p}^t-1}(w'))(s, \cdot) \, ds + w(t_{n-1}, \cdot) \quad \text{in } \mathcal{X} \quad \forall n = 1, \dots, N. \quad (16) $$

Combining (16) with Lemma 2.3, and using [32, Theorem 3.10], we get the following result.

**Lemma 2.4.** Let $\mathcal{P}_{\mathbf{p}^t}$ be the operator in (11). Then, for all $w$ in $H^{s+1}(I_n; \mathcal{X})$ with $s \ge 1$, the following inequalities hold true for all $n = 1, \dots, N$:

$$ \|(w - \mathcal{P}_{\mathbf{p}^t} w)'\|_{L^2(I_n; \mathcal{X})} \lesssim \|w' - \Pi_{\mathbf{p}^t-1}^0 w'\|_{L^2(I_n; \mathcal{X})} + \frac{\tau_n}{p_n^t} \|w''\|_{L^2(I_n; \mathcal{X})}, \quad (17a) $$

$$ \begin{aligned} \|(w - \mathcal{P}_{\mathbf{p}^t} w)'\|_{L^2(I_n; \mathcal{X})} & \lesssim \inf_{q_{p_n^t} \in \mathbb{P}_{p_n^t}(I_n; \mathcal{X})} \left( \|(w - q_{p_n^t})'\|_{L^2(I_n; \mathcal{X})} + \frac{\tau_n}{p_n^t} \|(w - q_{p_n^t})''\|_{L^2(I_n; \mathcal{X})} \right) \\ & \lesssim \left( \frac{\tau_n}{p_n^t} \right)^{\min(s, p_n^t)} \|w^{(s+1)}\|_{L^2(I_n; \mathcal{X})}, \end{aligned} \quad (17b) $$

$$ \|(\mathcal{P}_{\mathbf{p}^t} w)'\|_{L^2(I_n; \mathcal{X})} \lesssim \|w'\|_{L^2(I_n; \mathcal{X})} + \frac{\tau_n}{p_n^t} \|w''\|_{L^2(I_n; \mathcal{X})}. \quad (17c) $$

*Proof.* Inequality (17a) follows from (14a) and (16). Inequality (17c) follows from (14c) and (16). Inequality (17b) follows from (17c), the triangle inequality, and the fact that $\mathcal{P}_{\mathbf{p}^t}$ preserves polynomials of degree $p_n^t$ in time for all $n = 1, \dots, N$. $\square$

We derive bounds for the operator $\mathcal{P}_{\mathbf{p}^t}$ in $L^2$-type norms. To this aim, we introduce the piecewise $H^1$ projector $\Pi_{\mathbf{p}^t-1}^1$ onto polynomials in time as follows:

$$ \left\{ \begin{array}{ll} \left( (w - \Pi_{\mathbf{p}^t-1}^1 w)' , q_{p_n^t-1}' \right)_{0, I_n} = 0 & \forall w \in H^1(I_n; \mathcal{X}), \quad q_{p_n^t-1} \in \mathbb{P}_{p_n^t-1}(I_n; \mathcal{X}); \\ (w - \Pi_{\mathbf{p}^t-1}^1 w)(t_{n-1}, \cdot) = 0 & \text{in } \mathcal{X} \quad \forall n = 1, \dots, N. \end{array} \right. \quad (18) $$

On the other hand, the operator $\Pi_{\mathbf{p}^t-1}^0$ denotes the piecewise $L^2$ projector onto polynomials in time.

Some properties of the operator $\mathcal{P}_{\mathbf{p}^t}$ are detailed in the next result.

<footnote n="3">³With an abuse of notation, we use the same notation for the operators on the reference interval.</footnote>
**Proposition 2.5.** Let $\mathcal{P}_{\mathbf{p}^t}$ and $\Pi_{\mathbf{p}^t-1}^1$ be the operators in (11) and (18). Then, the following inequality holds true: for all $w$ in $H^{s+1}(I_n; \mathcal{X})$ with $s \ge 1$

$$ \|w - \mathcal{P}_{\mathbf{p}^t} w\|_{L^2(I_n; \mathcal{X})} \lesssim \|w - \Pi_{\mathbf{p}^t}^1 w\|_{L^2(I_n; \mathcal{X})} + \frac{\tau_n^2}{(p_n^t)^2} \|w''\|_{L^2(I_n; \mathcal{X})} \quad \forall n = 1, \dots, N. \quad (19) $$

Since $\mathcal{P}_{\mathbf{p}^t}$ preserves polynomials in time of order $p_n^t$, we also have

$$ \begin{aligned} \|w - \mathcal{P}_{\mathbf{p}^t} w\|_{L^2(I_n; \mathcal{X})} &\lesssim \inf_{q_{p_n^t} \in \mathbb{F}_{\mathbf{p}_n^t}(I_n; \mathcal{X})} \left( \frac{\tau_n}{p_n^t} \|(w - q_{p_n^t})'\|_{L^2(I_n; \mathcal{X})} + \frac{\tau_n^2}{(p_n^t)^2} \|(w - q_{p_n^t})''\|_{L^2(I_n; \mathcal{X})} \right) \\ &\lesssim \left( \frac{\tau_n}{p_n^t} \right)^{\min(s, p_n^t)+1} \|w^{(s+1)}\|_{L^2(I_n; \mathcal{X})} \quad \forall n = 1, \dots, N. \end{aligned} \quad (20) $$

*Proof.* We prove the assertion on the reference interval $\hat{I} = (-1, 1)$; the general bound follows from a scaling argument.

Identity (15) implies

$$ (w - \mathcal{P}_{\mathbf{p}^t} w)'(t, \cdot) \stackrel{(16)}{=} (w' - \tilde{\mathcal{P}}_{\mathbf{p}^t-1}(w'))(t, \cdot) = (w' - \Pi_{\mathbf{p}^t-1}^0 w')(t, \cdot) + (w' - \Pi_{\mathbf{p}^t-1}^0 w')(1, \cdot) L_{p_n^t-1}(t) \quad \text{in } \mathcal{X}. $$

Since $(w - \mathcal{P}_{\mathbf{p}^t} w)(-1, \cdot) = 0$, an integration by parts gives

$$ \begin{aligned} (w - \mathcal{P}_{\mathbf{p}^t} w)(t, \cdot) &= \int_{-1}^{t} (w - \mathcal{P}_{\mathbf{p}^t} w)'(s, \cdot) ds \\ &= \int_{-1}^{t} (w' - \Pi_{\mathbf{p}^t-1}^{0} w')(s, \cdot) ds + (w' - \Pi_{\mathbf{p}^t-1}^{0} w')(1, \cdot) \int_{-1}^{t} L_{p_n^t-1}(s) ds =: T_1 + T_2. \end{aligned} $$

We estimate the two terms on the right-hand side separately. The properties of the $L^2$ and $H^1$ projectors imply $\Pi_{\mathbf{p}^t-1}^0 w' = (\Pi_{\mathbf{p}^t}^1 w)'$. Therefore, recalling from (18) that $w(-1, \cdot) = \Pi_{\mathbf{p}^t}^1 w(-1, \cdot)$, we write

$$ T_1 = \int_{-1}^{t} (w' - (\Pi_{\mathbf{p}^t}^1 w)')(s, \cdot) ds = (w - \Pi_{\mathbf{p}^t}^1 w)(t, \cdot) \quad \text{in } \mathcal{X}. $$

This identity leads to the estimate on the first term on the right-hand side of (19).

As for the term $T_2$, we expand $w'$ with respect to Legendre polynomials:

$$ w'(t, \cdot) = \sum_{j=0}^{+\infty} w_j'(\cdot) L_j(t) \quad \text{in } \mathcal{X}, \qquad w_j' \in \mathcal{X}. $$

Standard properties of the $L^2$ projector imply

$$ T_2 = \left( \sum_{j=p_n^t}^{+\infty} w_j'(\cdot) \right) \int_{-1}^{t} L_{p_n^t-1}(s) ds \quad \text{in } \mathcal{X}. $$

From [32, eq. (3.5)], we have

$$ \left| \sum_{j=p_n^t}^{+\infty} w_j'(\cdot) \right| \lesssim (p_n^t)^{-\frac{1}{2}} \|w''(\cdot)\|_{0, \hat{I}} \quad \text{in } \mathcal{X}. $$

Using that $\int_{-1}^{t} L_{p_n^t-1}(s) ds = (2p_n^t - 1)^{-1}(L_{p_n^t} - L_{p_n^t-2})$ and $\|L_{p_n^t}\|_{0, \hat{I}} \approx (p_n^t)^{-\frac{1}{2}}$ yields

$$ \left\| \int_{-1}^{t} L_{p_n^t-1}(s) ds \right\|_{0, \hat{I}} \lesssim (p_n^t)^{-1} \left( \|L_{p_n^t}\|_{0, \hat{I}} + \|L_{p_n^t-2}\|_{0, \hat{I}} \right) \lesssim (p_n^t)^{-\frac{3}{2}}. $$

Collecting the two displays above implies

$$ T_2 \lesssim (p_n^t)^{-2} \|w''(\cdot)\|_{0, \hat{I}} \quad \text{in } \mathcal{X}. $$

This concludes the proof of (19). As for (20), we add and subtract $q_{p_n^t}$, employ (19), note that $\mathcal{P}_{\mathbf{p}^t}$ piecewise preserves piecewise polynomials of degree $p_n^t$ in time, and use the properties of the projector $\Pi_{\mathbf{p}^t}$ as in [34, Corollary 3.15]. $\square$
## 2.3 Abstract error analysis

Let $\Pi_h^\xi: H_0^1(\Omega) \to V_h$ denote the elliptic projector defined as

$$a(u - \Pi_h^\xi u, v_h) = 0 \quad \forall v_h \in V_h. \qquad (21)$$

We have the following stability estimate for the elliptic projector $\Pi_h^\xi$ in (21), which is a consequence of the Aubin-Nitsche technique: for any nonnegative $s$ larger than or equal to 1, and possibly equal to $\infty$, there exists $\alpha$ in $(0, 1]$ depending on $\Omega$ such that

$$\|\Pi_h^\xi v\|_{L^s(I_n; L^2(\Omega))} \lesssim \|v\|_{L^s(I_n; L^2(\Omega))} + \frac{h^\alpha}{p_x^\alpha} |v|_{L^s(I_n; H^1(\Omega))} \quad \forall n = 1, \dots, N. \qquad (22)$$

Estimate (22) holds true under the assumption of quasi-uniformity of the sequence of spatial meshes. This is the only occasion in the manuscript where local quasi-uniformity cannot be used.

**Remark 2.** If local quasi-uniformity is assumed, then we can still derive optimal error estimates in $h$; however, optimal error estimates in $p$ are still an open problem to the best of our knowledge. For convex domains, [3] (1D) and [28] (2D and 3D) proved optimal error estimates in the $L^2$ norm under a suitable graded mesh assumption; however, their proof hinges upon proving a discrete inf-sup condition in mesh dependent norms based on employing several inverse inequalities, whence the estimates cannot be optimal in the polynomial degree. In [29], analogous results are derived based on Aubin-Nitsche tricks in 2D polygonal domains based on weighted norms.

Recall that $\mathcal{P}_{\mathbf{p}^t}$ is given in (11) and let

$$\tilde{U}_h = \mathcal{P}_{\mathbf{p}^t} \Pi_h^\xi u.$$

We analyze a priori estimates of the following error quantity, which we split into two contributions:

$$e_h = u - U_h = (u - \tilde{U}_h) + (\tilde{U}_h - U_h) =: \rho + \xi_h. \qquad (23)$$

We show an upper bound on a suitable norm of the two terms on the right-hand side of (23). We begin discussing the term $\xi_h$.

**Lemma 2.6.** Consider $m = m(\xi_h)$, $\mu_m$, and $\xi_h$ as in (6), (7), and (23). Let $u$ and $U_h$ be the solutions to (2) and (3). Assume that $\Delta_x u$ and $u''$ belong to $H^2(I_n; L^2(\Omega))$ and $L^2(I_n; H_0^1(\Omega))$ for all $n = 1, \dots, m$. Recall that their initial conditions are $u_0$ and $u_1$, and $u_{0,h}$ and $u_{1,h}$, respectively. Then, the following estimate holds true:

$$\begin{aligned} & \mu_m \left( \|\xi_h'\|_{L^\infty(I_m; L^2(\Omega))}^2 + |\xi_h|_{L^\infty(I_m; H^1(\Omega))}^2 \right) + \frac{1}{4} \sum_{n=1}^m \| [[\xi_h']](t_{n-1}, \cdot) \|_{0,\Omega}^2 \\ & \le \frac{1}{2} \left[ \begin{aligned}[t] & (u_0 - u_{0,h})(\cdot)|_{1,\Omega}^2 + \|(u_1 - u_{1,h})(\cdot)\|_{0,\Omega}^2 + 2\|u_1(\cdot) - (\mathcal{P}_{\mathbf{p}^t} u)'(0, \cdot)\|_{0,\Omega}^2 \\ & \quad + 2\|(\mathcal{P}_{\mathbf{p}^t}(I - \Pi_h^\xi)u)'(0, \cdot)\|_{0,\Omega}^2 \end{aligned} \right. \\ & \left. \qquad + 2 \frac{t_m}{\mu_m} \left( \sum_{n=1}^m \|(I - \mathcal{P}_{\mathbf{p}^t}) \Delta_x u\|_{L^2(I_n; L^2(\Omega))}^2 + \sum_{n=1}^m \|(I - \Pi_h^\xi) u''\|_{L^2(I_n; L^2(\Omega))}^2 \right) \right]. \end{aligned}$$

The first term on the right-hand side above vanishes if we discretize the initial condition $u_0$ with the elliptic projector in (21).

*Proof.* Let $e_h$ be as in (23). Subtracting (2) and (3), for all $n = 1, \dots, N$, we arrive at

$$\int_{I_n} [(e_h'', W_h)_{0,\Omega} + a(e_h, W_h)]dt + ([[e_h']](t_{n-1}, \cdot), W_h(t_{n-1}, \cdot))_{0,\Omega} = 0 \quad \forall W_h \in \mathbb{P}_{p_n'-1}^t(I_n; V_h).$$

Splitting (23) yields

$$\begin{aligned} & \int_{I_n} [(\xi_h'', W_h)_{0,\Omega} + a(\xi_h, W_h)]dt + ([[\xi_h']](t_{n-1}, \cdot), W_h(t_{n-1}, \cdot))_{0,\Omega} \\ & = - \int_{I_n} (\rho'', W_h)_{0,\Omega}dt - \int_{I_n} a(\rho, W_h)dt - ([[\rho']](t_{n-1}, \cdot), W_h(t_{n-1}, \cdot))_{0,\Omega} \quad \forall W_h \in \mathbb{P}_{p_n'-1}^t(I_n; V_h). \end{aligned} \qquad (24)$$
We rewrite the right-hand side of (24) as $(f_{\xi_h}, V)_{L^2(0,T;L^2(\Omega))}$, i.e., as the right-hand side of method (3), for a suitable $f_{\xi_h}$, which we fix in (27), so as to exploit the stability estimates (8) and deduce the assertion.

We focus on the second term on the right-hand side of (24). Using the fact that $\mathcal{P}_{\mathbf{p}^t}$ and $\Pi_h^\xi$ commute, definition (21) of $\Pi_h^\xi$, the fact that $\mathcal{P}_{\mathbf{p}^t}$ preserves polynomials in time, and an integration by parts, we arrive at

$$ 
\begin{aligned}
\int_{I_n} a(\rho, W_h)dt &= \int_{I_n} a((I - \mathcal{P}_{\mathbf{p}^t} \Pi_h^\xi)u, W_h)dt \\
&= \int_{I_n} a((I - \mathcal{P}_{\mathbf{p}^t})u, W_h)dt = -(\Delta_{\mathbf{x}}(I - \mathcal{P}_{\mathbf{p}^t})u, W_h)_{L^2(I_n; L^2(\Omega))}.
\end{aligned} 
\quad (25) $$

As for the first and third terms on the right-hand side of (24), we use the continuity in time of $u'$ and (12), and get

$$ 
\begin{aligned}
& (\rho'', W_h)_{L^2(I_n; L^2(\Omega))} + (\llbracket \rho' \rrbracket(t_{n-1}, \cdot), W_h(t_{n-1}, \cdot))_{0,\Omega} \\
&= (u'', W_h)_{L^2(I_n; L^2(\Omega))} + (\llbracket u' \rrbracket(t_{n-1}, \cdot), W_h(t_{n-1}, \cdot))_{0,\Omega} \\
&\quad - ((\mathcal{P}_{\mathbf{p}^t} \Pi_h^\xi u)'', W_h)_{L^2(I_n; L^2(\Omega))} - (\llbracket (\mathcal{P}_{\mathbf{p}^t} \Pi_h^\xi u)' \rrbracket(t_{n-1}, \cdot), W_h(t_{n-1}, \cdot))_{0,\Omega} \\
&= (u'', W_h)_{L^2(I_n; L^2(\Omega))} - ((\mathcal{P}_{\mathbf{p}^t} \Pi_h^\xi u)'', W_h)_{L^2(I_n; L^2(\Omega))} - (\llbracket (\mathcal{P}_{\mathbf{p}^t} \Pi_h^\xi u)' \rrbracket(t_{n-1}, \cdot), W_h(t_{n-1}, \cdot))_{0,\Omega} \\
&\stackrel{(12)}{=} (u'', W_h)_{L^2(I_n; L^2(\Omega))} - ((\Pi_h^\xi u)'', W_h)_{L^2(I_n; L^2(\Omega))}.
\end{aligned} 
$$

For all $n = 1, \dots, N$, we deduce

$$ [(\rho'', W_h)_{L^2(I_n; L^2(\Omega))} + (\llbracket \rho' \rrbracket(t_{n-1}, \cdot), W_h(t_{n-1}, \cdot))_{0,\Omega} = ((I - \Pi_h^\xi)u'', W_h)_{L^2(I_n; L^2(\Omega))}. \quad (26) $$

We insert (25) and (26) in (24) and arrive at

$$ 
\begin{aligned}
& \int_{I_n} [(\xi_h'', W_h)_{0,\Omega} + a(\xi_h, W_h)]dt + (\llbracket \xi_h' \rrbracket(t_{n-1}, \cdot), W_h(t_{n-1}, \cdot))_{0,\Omega} \\
&= ((I - \mathcal{P}_{\mathbf{p}^t})\Delta_{\mathbf{x}}u, W_h)_{L^2(I_n; L^2(\Omega))} - ((I - \Pi_h^\xi)u'', W_h)_{L^2(I_n; L^2(\Omega))} \\
&= ((I - \mathcal{P}_{\mathbf{p}^t})\Delta_{\mathbf{x}}u - (I - \Pi_h^\xi)u'', W_h)_{L^2(I_n; L^2(\Omega))} \quad \forall W_h \in \mathbb{P}_{p_n-1}^t(I_n; V_h), \quad \forall n = 1, \dots, N.
\end{aligned} 
$$

In words, $\xi_h$ is the solution to method (3) with right-hand side given by $(f_{\xi_h}, V)_{L^2(0,T;L^2(\Omega))}$ with

$$ f_{\xi_h} := (I - \mathcal{P}_{\mathbf{p}^t})\Delta_{\mathbf{x}}u - (I - \Pi_h^\xi)u''. \quad (27) $$

We are now in a position to apply the discrete stability estimates (8). Let $m = m(\xi_h)$ be defined in (6) and get

$$ 
\begin{aligned}
& \mu_m \left( \|\xi_h'\|_{L^\infty(I_m; L^2(\Omega))}^2 + |\xi_h|_{L^\infty(I_m; H^1(\Omega))}^2 \right) + \frac{1}{4} \sum_{n=1}^m \|\llbracket \xi_h \rrbracket(t_n, \cdot)\|_{0,\Omega}^2 \\
&\le \frac{1}{2} \left( |\xi_h(0, \cdot)|_{1,\Omega}^2 + \|(\xi_h^-)'(0, \cdot)\|_{0,\Omega}^2 \right) + \frac{t_m}{\mu_m} \|f_{\xi_h}\|_{L^2(0, t_m; L^2(\Omega))}^2 =: T_1 + T_2.
\end{aligned} 
$$

We estimate the two terms on the right-hand side. The triangle inequality implies

$$ T_1 \le \frac{1}{2} \left( |\xi_h(0, \cdot)|_{1,\Omega}^2 + \|(\xi_h^-)'(0, \cdot)\|_{0,\Omega}^2 \right) = \frac{1}{2} [T_{1,1} + T_{1,2}]. $$

Using (11) and the stability of the elliptic projection $\Pi_h^\xi$, we deduce

$$ T_{1,1} = |\xi_h(0, \cdot)|_{1,\Omega}^2 = |U_h(0, \cdot) - \Pi_h^\xi u(0, \cdot)|_{1,\Omega}^2 \le |(u_0 - u_{0,h})(\cdot)|_{1,\Omega}^2. $$

The term $T_{1,1}$ vanishes if the initial condition $u_0$ is discretized with the elliptic projector in (21). We further have

$$ 
\begin{aligned}
T_{1,2} &= \|(\xi_h^-)'(0, \cdot)\|_{0,\Omega}^2 = \|(U_h^-)'(0, \cdot) - (\mathcal{P}_{\mathbf{p}^t} \Pi_h^\xi u)'(0, \cdot)\|_{0,\Omega}^2 \\
&\le 2\|(u_1 - u_{1,h})(\cdot)\|_{0,\Omega}^2 + 2\|u_1(\cdot) - (\mathcal{P}_{\mathbf{p}^t} \Pi_h^\xi u)'(0, \cdot)\|_{0,\Omega}^2 \\
&\le 2\|(u_1 - u_{1,h})(\cdot)\|_{0,\Omega}^2 + 4\|u_1(\cdot) - (\mathcal{P}_{\mathbf{p}^t} u)'(0, \cdot)\|_{0,\Omega}^2 + 4\|(\mathcal{P}_{\mathbf{p}^t}(I - \Pi_h^\xi)u)'(0, \cdot)\|_{0,\Omega}^2.
\end{aligned} 
$$
Next, we deal with the term $T_2$:

$$T_2 = \frac{t_m}{\mu_m} \|f_{\xi_h}\|_{L^2(0,t_m;L^2(\Omega))}^2 := \frac{t_m}{\mu_m} \sum_{n=1}^{m} \|(I - \mathcal{P}_{\mathbf{p}^t}) \Delta_{\mathbf{x}} u - (I - \Pi_h^{\mathcal{E}}) u''\|_{L^2(I_n;L^2(\Omega))}^2 := \frac{t_m}{\mu_m} \sum_{n=1}^{m} T_{2,n}.$$

We have

$$T_{2,n} \le 2 \|(I - \mathcal{P}_{\mathbf{p}^t}) \Delta_{\mathbf{x}} u\|_{L^2(I_n; L^2(\Omega))}^2 + 2 \|(I - \Pi_h^{\mathcal{E}}) u''\|_{L^2(I_n; L^2(\Omega))}^2.$$

Collecting the bounds on $T_1$ and $T_2$ yields the assertion.

□

Introduce $k$ in $1, \dots, N$ such that

$$\left( \|e_h'\|_{L^\infty(I_k; L^2(\Omega))}^2 + |e_h|_{L^\infty(I_k; H^1(\Omega))}^2 \right) := \max_{n=1}^N \left( \|e_h'\|_{L^\infty(I_n; L^2(\Omega))}^2 + |e_h|_{L^\infty(I_n; H^1(\Omega))}^2 \right). \quad (28)$$

The index $k$ is defined similarly to the index $m$ in (6). However, on the one hand, with $k$ we maximize positive functionals in Bochner spaces and not in space-time finite element spaces; on the other hand, using a different nomenclature of the indices allows us to improve the readability of the estimates in Theorems 2.7 and 2.8 below.

**Theorem 2.7.** Consider $m = m(\xi_h)$, $\mu_m$, $k$, and $\xi_h$ and $\rho$ as in (6), (7), (28), and (23). Let $u$ and $U_h$ be the solutions to (2) and (3). Assume that $\Delta_{\mathbf{x}}u$ and $u''$ belong to $H^2(I_n; L^2(\Omega))$ and $L^2(I_n; H_0^1(\Omega))$ for all $n = 1, \dots, m$, respectively. Recall that their initial conditions are $u_0$ and $u_1$, and $u_{0,h}$ and $u_{1,h}$, respectively. Then, the following estimate holds true:

$$\begin{aligned} \max_{n=1}^{N} \left( \|e_h'\|_{L^\infty(I_n; L^2(\Omega))}^2 + |e_h|_{L^\infty(I_n; H^1(\Omega))}^2 \right) \le \frac{1}{\mu_m} \left( \left| (u_0 - u_{0,h})(\cdot) \right|_{1,\Omega}^2 \right. \\ \left. + 2 \|(u_1 - u_{1,h})(\cdot)\|_{0,\Omega}^2 + 4 \|u_1(\cdot) - (\mathcal{P}_{\mathbf{p}^t} u)'(0, \cdot)\|_{0,\Omega}^2 + 4 \|(\mathcal{P}_{\mathbf{p}^t} (I - \Pi_h^{\mathcal{E}}) u)'(0, \cdot)\|_{0,\Omega}^2 \right) \\ + \frac{4t_m}{\mu_m^2} \sum_{n=1}^{m} \|(I - \mathcal{P}_{\mathbf{p}^t}) \Delta_{\mathbf{x}} u\|_{L^2(I_n; L^2(\Omega))}^2 + \frac{4t_m}{\mu_m^2} \sum_{n=1}^{m} \|(I - \Pi_h^{\mathcal{E}}) u''\|_{L^2(I_n; L^2(\Omega))}^2 \\ + 2 \|\rho'\|_{L^\infty(I_k; L^2(\Omega))}^2 + 2 |\rho|_{L^\infty(I_k; H^1(\Omega))}^2 =: \sum_{j=1}^{8} \mathfrak{J}_j^2.
\end{aligned} \quad (29)$$

*Proof.* Using the triangle inequality and the definition of $m = m(\xi_h)$ in (6), we deduce

$$\begin{aligned} \max_{n=1}^{N} \left( \|e_h'\|_{L^\infty(I_n; L^2(\Omega))}^2 + |e_h|_{L^\infty(I_n; H^1(\Omega))}^2 \right) & \\ \le 2 \left( \|\rho'\|_{L^\infty(I_k; L^2(\Omega))}^2 + |\rho|_{L^\infty(I_k; H^1(\Omega))}^2 \right) + 2 \left( \|\xi_h'\|_{L^\infty(I_k; L^2(\Omega))}^2 + |\xi_h|_{L^\infty(I_k; H^1(\Omega))}^2 \right) & \\ \le 2 \left( \|\rho'\|_{L^\infty(I_k; L^2(\Omega))}^2 + |\rho|_{L^\infty(I_k; H^1(\Omega))}^2 \right) + 2 \left( \|\xi_h'\|_{L^\infty(I_m; L^2(\Omega))}^2 + |\xi_h|_{L^\infty(I_m; H^1(\Omega))}^2 \right). & \end{aligned}$$

The assertion follows using Lemma 2.6.

□

## 2.4 Error estimates

For $m = m(\xi_h)$ and $\mu_m$ as in (6) and (7), we derive error estimates for method (3), which are explicit in the spatial mesh size, the time steps, and the polynomial degrees in space and time, with respect to the norm

$$\max_{n=1}^{N} \left( \|\cdot'\|_{L^\infty(I_n; L^2(\Omega))}^2 + |\cdot|_{L^\infty(I_n; H^1(\Omega))}^2 \right). \quad (30)$$

To this aim, we give explicit bounds on the terms $\mathfrak{J}_j$, $j = 1, \dots, 8$, appearing on the right-hand side of (29), and collect the resulting estimates in the following result.

**Theorem 2.8.** Let $u$ and $U_h$ be the solutions to (2) and (3), and $e_h$ be as in (23); $u_0$, $u_1$, and $u$ be sufficiently smooth; $m = m(\xi_h)$, $\mu_m$, and $k$ be as in (6), (7), and (28); $\alpha$ be the elliptic regularity parameter as in (22). Then, the following a priori estimate holds true:

$$\max_{n=1}^{N} \left( \|e_h'\|_{L^\infty(I_n; L^2(\Omega))}^2 + |e_h|_{L^\infty(I_n; H^1(\Omega))}^2 \right) \lesssim \tau_1 + \tau_2 + \tau_3,$$
where we have set

$$
\begin{align*}
\mathfrak{J}_1 := & \frac{1}{\mu_m} \left( \frac{h^{2 \min(p_x, s-\frac{1}{2})}}{p_x^{2s-1}} \|u_0(\cdot)\|_{s+\frac{1}{2}, \Omega}^2 \right. \\
& \left. + \frac{h^{2 \min(p_x, 1+s-\frac{1}{2})}}{p_x^{2s-1}} \|u_1(\cdot)\|_{s-\frac{1}{2}, \Omega}^2 + \left(\frac{\tau_1}{p_1^t}\right)^{2 \min(p_1^t, s)-1} \|u^{(s+1)}\|_{L^2(I_1, L^2(\Omega))}^2 \right. \\
& \left. + \left(\frac{\tau_1}{p_t^t}\right)^{h^{2 \min(p_x+1, s-1)-2(1-\alpha)}} \|u''\|_{L^2(I_1; H^{s-1}(\Omega))}^2 + \frac{h^{2 \min(p_x+1, s-\frac{1}{2})-2(1-\alpha)}}{p_x^{2s-2-2(1-\alpha)}} \|u_1(\cdot)\|_{s-\frac{1}{2}, \Omega}^2 \right),
\end{align*}
$$

$$
\begin{align*}
\mathfrak{J}_2 := & \frac{t_m}{\mu_m^2} \sum_{n=1}^{m} \left(\frac{\tau_n}{p_n^t}\right)^{2 \min(s-1, p_n^t+1)} \left\|\Delta_x u^{(s-1)}\right\|_{L^2(I_n; L^2(\Omega))}^2 \\
& + \frac{t_m}{\mu_m^2} \sum_{n=1}^{m} \frac{h^{2 \min(p_x+1, s-1)-2(1-\alpha)}}{p_x^{2s-2-2(1-\alpha)}} \left\|u''\right\|_{L^2(I_n; H^{s-1}(\Omega))}^2,
\end{align*}
$$

and

$$
\begin{align*}
\mathfrak{J}_3 := & \left( \frac{h^{2 \min(s, p_k^t)+2 \tau_k}}{(p_k^t)^{2s+2}} \left( 1 + \frac{h^{-2+2\alpha}}{p_x^{-2+2\alpha}} \right) \left( \tau_k^{-1} \|u'\|_{L^\infty(I_k; H^{s+1}(\Omega))}^2 \right) \right. \\
& \left. + \frac{\tau_k^{2 \min(s, p_k^t)-1}}{(p_k^t)^{2 \min(s, p_k^t)-2}} \left( \|u^{(s+1)}\|_{L^2(I_k; L^2(\Omega))}^2 + \frac{h^{2\alpha-2}}{p_x^{2\alpha}} h^2 \|u^{(s+1)}\|_{L^2(I_k; H^1(\Omega))}^2 \right) \right) \\
& + \left( \frac{h^{2 \min(p_x, s)} \tau_k^{-1}}{p_x^{2s}} \left( \tau_k \|u\|_{L^\infty(I_k; H^{s+1}(\Omega))} \right)^2 + \left( \frac{\tau_k}{p_k^t} \right)^{2 \min(s, p_k^t+1)-1} \|u^{(s)}\|_{L^2(I_k; H^1(\Omega))}^2 \right).
\end{align*}
$$

Proof. Estimates on $\mathfrak{J}_1$. We recall that $u_{0,h}$ is the elliptic projection of $u_0$ as in (21) onto $V_h$. Standard polynomial approximation estimates imply

$$
\mu_m \mathfrak{J}_1^2 = \|(u_0 - u_{0,h})(\cdot)\|_{1,\Omega}^2 \lesssim \frac{h^{2 \min(p_x, s-\frac{1}{2})}}{p_x^{2s-1}} \|u_0(\cdot)\|_{s+\frac{1}{2}, \Omega}^2. \quad (31)
$$

**Estimates on $\mathfrak{J}_2$.** We recall that $u_{1,h}$ is the $L^2$ projection of $u_1$ onto $V_h$. Standard polynomial approximation estimates imply

$$
\frac{\mu_m}{2} \mathfrak{J}_2^2 = \|(u_1 - u_{1,h})(\cdot)\|_{0,\Omega}^2 \lesssim \frac{h^{2 \min(p_x+1, s-\frac{1}{2})}}{p_x^{2s-1}} \|u_1(\cdot)\|_{s-\frac{1}{2}, \Omega}^2. \quad (32)
$$

Estimates on $\mathfrak{J}_3$. Using identity (16) and the approximation estimate (14d), we infer

$$
\frac{\mu_m}{4} \mathfrak{J}_3^2 = \|(u - \mathcal{P}_{\mathbf{p}^t} u)'(0, \cdot)\|_{0,\Omega}^2 = \|(u' - \tilde{\mathcal{P}}_{\mathbf{p}^t-1} u')(0, \cdot)\|_{0,\Omega}^2 \lesssim \left(\frac{\tau_1}{p_1^t}\right)^{2 \min(p_1^t, s)-1} \|u^{(s+1)}\|_{L^2(I_1, L^2(\Omega))}^2. \quad (33)
$$

Estimates on $\mathfrak{J}_4$. Using the triangle inequality, we have

$$
\frac{\mu_m}{4} \mathfrak{J}_4^2 = \left\| (\mathcal{P}_{\mathbf{p}^t}(I - \Pi_h^\mathcal{E})u)'(0, \cdot) \right\|_{0,\Omega}^2 \le 2 \left\| \left( (I - \mathcal{P}_{\mathbf{p}^t})(I - \Pi_h^\mathcal{E})u \right)'(0, \cdot) \right\|_{0,\Omega}^2 + 2 \left\| (I - \Pi_h^\mathcal{E})u'(0, \cdot) \right\|_{0,\Omega}^2 =: \mathfrak{J}_{4,1}^2 + \mathfrak{J}_{4,2}^2.
$$

Following the proof of the bound on $\mathfrak{J}_3^2$ we get the following bound on $\mathfrak{J}_{4,1}^2$:

$$
\begin{align*}
\mathfrak{J}_{4,1}^2 & \stackrel{(14d)}{\lesssim} \left(\frac{\tau_1}{p_1^t}\right) \|(I - \Pi_h^\mathcal{E}) u''\|_{L^2(I_1; L^2(\Omega))}^2 \\
& \stackrel{(22)}{\lesssim} \left(\frac{\tau_1}{p_1^t}\right) \inf_{q_{p_x} \in L^2(I_1; V_h)} \left( \|u'' - q_{p_x}\|_{L^2(I_1; L^2(\Omega))}^2 + \|\Pi_h^\mathcal{E}(u'' - q_{p_x})\|_{L^2(I_1; L^2(\Omega))}^2 \right) \\
& \stackrel{(22)}{\lesssim} \left(\frac{\tau_1}{p_1^t}\right) \inf_{q_{p_x} \in L^2(I_1; V_h)} \left( \|u'' - q_{p_x}\|_{L^2(I_1; L^2(\Omega))}^2 + \frac{h^{2\alpha}}{p_x^{2\alpha}} |u'' - q_{p_x}|_{L^2(I_1; H^1(\Omega))}^2 \right) \\
& \stackrel{(22)}{\lesssim} \left(\frac{\tau_1}{p_1^t}\right) \left( \frac{h^{2 \min(p_x+1, s-1)}}{p_x^{2(s-1)}} \|u''\|_{L^2(I_1; H^{s-1}(\Omega))}^2 + \frac{h^{2\alpha} h^{2 \min(p_x+1, s-1)-2}}{p_x^{2(s-1)-2}} \|u''\|_{L^2(I_1; H^{s-1}(\Omega))}^2 \right) \\
& \stackrel{(22)}{\lesssim} \left(\frac{\tau_1}{p_1^t}\right) \frac{h^{2 \min(p_x+1, s-1)-2(1-\alpha)}}{p_x^{2s-2-2(1-\alpha)}} \|u''\|_{L^2(I_1; H^{s-1}(\Omega))}^2.
\end{align*}
$$
As for the term $\mathfrak{J}_{4,2}$, we proceed similarly:

$$
\begin{aligned}
\mathfrak{J}_{4,2}^2 &= \left\| (I - \Pi_h^\mathcal{E}) u_1(\cdot) \right\|_{0,\Omega}^2 \lesssim \inf_{q_{p_x} \in L^2(I_1; V_h)} \left( \left\| (u_1 - q_{p_x})(\cdot) \right\|_{0,\Omega}^2 + \left\| \Pi_h^\mathcal{E} (u_1 - q_{p_x})(\cdot) \right\|_{0,\Omega}^2 \right) \\
&\lesssim \left( \frac{h^{2 \min(p_x+1, s-\frac{1}{2})}}{p_x^{2s-1}} \|u_1(\cdot)\|_{s-\frac{1}{2}, \Omega}^2 + \frac{h^{2\alpha} h^{2 \min(p_x+1, s-\frac{1}{2})-2}}{p_x^{2\alpha} p_x^{2s-3}} \|u_1(\cdot)\|_{s-\frac{1}{2}, \Omega}^2 \right) \\
&\lesssim \frac{h^{2 \min(p_x+1, s-\frac{1}{2})-2(1-\alpha)}}{p_x^{2s-1-2(1-\alpha)}} \|u_1(\cdot)\|_{s-\frac{1}{2}, \Omega}^2.
\end{aligned}
$$

Collecting the two displays above gives

$$
\mu_m \mathfrak{J}_4^2 \lesssim \left( \frac{\tau_1}{p_1^t} \right) \frac{h^{2 \min(p_x+1, s-1)-2(1-\alpha)}}{p_x^{2s-2-2(1-\alpha)}} \|u''\|_{L^2(I_1; H^{s-1}(\Omega))}^2 + \frac{h^{2 \min(p_x+1, s-\frac{1}{2})-2(1-\alpha)}}{p_x^{2s-1-2(1-\alpha)}} \|u_1(\cdot)\|_{s-\frac{1}{2}, \Omega}^2. \quad (34)
$$

**Estimates on $\mathfrak{J}_5$.** Using Proposition 2.5, we get

$$
\frac{\mathfrak{J}_5^2}{4} = \frac{t_m}{\mu_m^2} \sum_{n=1}^{m} \|(I - \mathcal{P}_{\mathbf{p}^t}) \Delta_{\mathbf{x}} u\|_{L^2(I_n; L^2(\Omega))}^2 \lesssim \frac{t_m}{\mu_m^2} \sum_{n=1}^{m} \left(\frac{\tau_n}{p_n^t}\right)^{2 \min(s-1, p_n^t+1)} \|\Delta_{\mathbf{x}} u^{(s-1)}\|_{L^2(I_n; L^2(\Omega))}^2. \quad (35)
$$

**Estimates on $\mathfrak{J}_6$.** Using (22) and the fact that $\Pi_h^\mathcal{E}$ preserves spatial polynomials of degree $p_x$, we get

$$
\begin{aligned}
\frac{\mathfrak{J}_6^2}{4} &= \frac{t_m}{\mu_m^2} \sum_{n=1}^{m} \left\| (I - \Pi_h^\mathcal{E}) u'' \right\|_{L^2(I_n; L^2(\Omega))}^2 \\
&\lesssim \frac{t_m}{\mu_m^2} \sum_{n=1}^{m} \inf_{q_{p_x} \in L^2(I_n; V_h)} \left( \|u'' - q_{p_x}\|_{L^2(I_n; L^2(\Omega))}^2 + \|\Pi_h^\mathcal{E}(u'' - q_{p_x})\|_{L^2(I_n; L^2(\Omega))}^2 \right) \\
&\lesssim \frac{t_m}{\mu_m^2} \sum_{n=1}^{m} \inf_{q_{p_x} \in L^2(I_n; V_h)} \left( \|u'' - q_{p_x}\|_{L^2(I_n; L^2(\Omega))}^2 + \frac{h^{2\alpha}}{p_x^{2\alpha}} |u'' - q_{p_x}|_{L^2(I_n; H^1(\Omega))}^2 \right) \quad (36) \\
&\lesssim \frac{t_m}{\mu_m^2} \sum_{n=1}^{m} \left( \frac{h^{2 \min(p_x+1, s-1)}}{p_x^{2(s-1)}} \|u''\|_{L^2(I_n; H^{s-1}(\Omega))}^2 + \frac{h^{2\alpha} h^{2 \min(p_x+1, s-1)-2}}{p_x^{2\alpha} p_x^{2(s-1)-2}} \|u''\|_{L^2(I_n; H^{s-1}(\Omega))}^2 \right) \\
&\lesssim \frac{t_m}{\mu_m^2} \sum_{n=1}^{m} \frac{h^{2 \min(p_x+1, s-1)-2(1-\alpha)}}{p_x^{2s-2-2(1-\alpha)}} \|u''\|_{L^2(I_n; H^{s-1}(\Omega))}^2.
\end{aligned}
$$

**Estimates on $\mathfrak{J}_7$.** Using (22) and the definition of $\rho$ in (23) gives

$$
\begin{aligned}
\frac{\mathfrak{J}_7^2}{2} &= \|\rho'\|_{L^\infty(I_k; L^2(\Omega))}^2 = \|(u - \mathcal{P}_{\mathbf{p}^t} \Pi_h^\mathcal{E} u)'\|_{L^\infty(I_k; L^2(\Omega))}^2 \\
&\lesssim \left( \|(u - \Pi_h^\mathcal{E} u)'\|_{L^\infty(I_k; L^2(\Omega))}^2 + \|(u - \mathcal{P}_{\mathbf{p}^t} u)'\|_{L^\infty(I_k; L^2(\Omega))}^2 + \left(\frac{h}{p_x}\right)^{2\alpha} |(u - \mathcal{P}_{\mathbf{p}^t} u)'|_{L^\infty(I_k; H^1(\Omega))}^2 \right).
\end{aligned}
\quad (37)
$$

First, we focus on the first term on the right-hand side: for any $q_{p_x}$ in $W^{1,\infty}(I_k; V_h)$,

$$
\begin{aligned}
\|(u - \Pi_h^\mathcal{E} u)'\|_{L^\infty(I_k; L^2(\Omega))}^2 &\lesssim \|(u - q_{p_x})'\|_{L^\infty(I_k; L^2(\Omega))}^2 + \|\Pi_h^\mathcal{E}(u - q_{p_x})'\|_{L^\infty(I_k; L^2(\Omega))}^2 \\
&\lesssim \left\| (u - q_{p_x})' \right\|_{L^\infty(I_k; L^2(\Omega))}^2 + \frac{h^{2\alpha}}{p^{2\alpha}} |(u - q_{p_x})'|_{L^\infty(I_k; H^1(\Omega))}^2 \\
&\lesssim \frac{h^{2 \min(s, p_k^t)+2} \tau_k}{(p_k^t)^{2s+2}} (\tau_k^{-1} \|u'\|_{L^\infty(I_k; H^{s+1}(\Omega))}^2) + \frac{h^{2 \min(s, p_k^t)+2\alpha} \tau_k}{(p_k^t)^{2s+2\alpha}} (\tau_k^{-1} \|u'\|_{L^\infty(I_k; H^{s+1}(\Omega))}^2) \\
&= \frac{h^{2 \min(s, p_k^t)+2} \tau_k}{(p_k^t)^{2s+2}} \left(1 + \frac{h^{-2+2\alpha}}{p_x^{-2+2\alpha}}\right) (\tau_k^{-1} \|u'\|_{L^\infty(I_k; H^{s+1}(\Omega))}^2).
\end{aligned}
$$

Next, we deal with the second term in the parenthesis on the right-hand side: using a polynomial
inverse inequality as in [34, eq. (3.6.4)], for any $q_{p_k^t}$ in $\mathbb{P}_{p_k^t}(I_k; L^2(\Omega))$,

$$
\begin{align*}
\|(u - \mathcal{P}_{\mathbf{p}^t} u)'\|_{L^\infty(I_k; L^2(\Omega))}^2 &\lesssim \|(u - q_{p_k^t})'\|_{L^\infty(I_k; L^2(\Omega))}^2 + \|\mathcal{P}_{\mathbf{p}^t}'(u - q_{p_k^t})\|_{L^\infty(I_k; L^2(\Omega))}^2 \\
&\lesssim \|(u - q_{p_k^t})'\|_{L^\infty(I_k; L^2(\Omega))}^2 + \frac{(p_k^t)^2}{\tau_k} \|\mathcal{P}_{\mathbf{p}^t}'(u - q_{p_k^t})\|_{L^2(I_k; L^2(\Omega))}^2 \\
&\lesssim \|(u - q_{p_k^t})'\|_{L^\infty(I_k; L^2(\Omega))}^2 + \frac{(p_k^t)^2}{\tau_k} \|(u - q_{p_k^t})'\|_{L^2(I_k; L^2(\Omega))}^2 + \tau_k \|(u - q_{p_k^t})''\|_{L^2(I_k; L^2(\Omega))}^2.
\end{align*}
$$

Standard polynomial approximation properties in 1D give

$$
\begin{align*}
\|(u - q_{p_k^t})'\|_{L^\infty(I_k; L^2(\Omega))}^2 &\lesssim \frac{\tau_k^{2 \min(s, p_k^t) - 1}}{(p_k^t)^{2 \min(s, p_k^t) - 1}} \|u^{(s+1)}\|_{L^2(I_k; L^2(\Omega))}^2, \\
\frac{(p_k^t)^2}{\tau_k} \|(u - q_{p_k^t})'\|_{L^2(I_k; L^2(\Omega))}^2 &\lesssim \frac{\tau_k^{2 \min(s, p_k^t) - 1}}{(p_k^t)^{2 \min(s, p_k^t) - 2}} \|u^{(s+1)}\|_{L^2(I_k; L^2(\Omega))}^2,
\end{align*}
$$

and

$$
\tau_k \|(u - q_{p_k^t})''\|_{L^2(I_k; L^2(\Omega))}^2 \lesssim \frac{\tau_k^{2 \min(s, p_k^t) - 1}}{(p_k^t)^{2 \min(s, p_k^t) - 2}} \|u^{(s+1)}\|_{L^2(I_k; L^2(\Omega))}^2.
$$

We collect the four displays above and get

$$
\| (u - \mathcal{P}_{\mathbf{p}^t} u)' \|_{L^\infty(I_k; L^2(\Omega))}^2 \lesssim \frac{\tau_k^{2 \min(s, p_k^t) - 1}}{(p_k^t)^{2 \min(s, p_k^t) - 2}} \| u^{(s+1)} \|_{L^2(I_k; L^2(\Omega))}^2 .
$$

We proceed similarly for the third term on the right-hand side of (37):

$$
\begin{align*}
\left(\frac{h}{p_x}\right)^{2\alpha} \left\|\left(u - \mathcal{P}_{\mathbf{p}^t} u\right)'\right\|_{L^\infty(I_k; H^1(\Omega))}^2 &\lesssim \left(\frac{h}{p_x}\right)^{2\alpha} \frac{\tau_k^{2 \min(s, p_k^t) - 1}}{(p_k^t)^{2 \min(s, p_k^t) - 2}} \left\|u^{(s+1)}\right\|_{L^2(I_k; H^1(\Omega))}^2 \\
&= \frac{h^{2\alpha-2}}{p_x^{2\alpha}} \frac{\tau_k^{2 \min(s, p_k^t) - 1}}{(p_k^t)^{2 \min(s, p_k^t) - 2}} h^2 \left\|u^{(s+1)}\right\|_{L^2(I_k; H^1(\Omega))}^2.
\end{align*}
$$

Recalling (7), we arrive at

$$
\mathfrak{J}_7^2 \lesssim \left[ \frac{h^{2 \min(s, p_k^t) + 2} \tau_k}{(p_k^t)^{2s+2}} \left( 1 + \frac{h^{-2+2\alpha}}{p_x^{-2+2\alpha}} \right) \tau_k^{-1} \|u\|_{W^{1,\infty}(I_k; H^{s+1}(\Omega))}^2 \right. \\
\left. + \frac{\tau_k^{2 \min(s, p_k^t) - 1}}{(p_k^t)^{2 \min(s, p_k^t) - 2}} \left( \|u^{(s+1)}\|_{L^2(I_k; L^2(\Omega))}^2 + \frac{h^{2\alpha-2}}{p_x^{2\alpha}} h^2 \|u^{(s+1)}\|_{L^2(I_k; H^1(\Omega))}^2 \right) \right]. \quad (38)
$$

Estimates on $\mathfrak{J}_8$. Using the stability of the elliptic projector $\Pi_h^\epsilon$ in $H^1(\Omega)$ and the 1D Sobolev embedding in [24, eq. (1.3)], we write

$$
\begin{align}
\frac{\mathfrak{J}_8^2}{2} &= |\rho|^2_{L^\infty(I_k; H^1(\Omega))} = |u - \mathcal{P}_{\mathbf{p}^t} \Pi_h^\epsilon u|^2_{L^\infty(I_k; H^1(\Omega))} \nonumber \\
&\lesssim \left( |u - \Pi_h^\epsilon u|^2_{L^\infty(I_k; H^1(\Omega))} + |\Pi_h^\epsilon (\mathrm{Id} - \mathcal{P}_{\mathbf{p}^t}) u|^2_{L^\infty(I_k; H^1(\Omega))} \right) \tag{39} \\
&\le \left( |u - \Pi_h^\epsilon u|^2_{L^\infty(I_k; H^1(\Omega))} + |u - \mathcal{P}_{\mathbf{p}^t} u|^2_{L^\infty(I_k; H^1(\Omega))} \right) \nonumber \\
&\le \left( |u - \Pi_h^\epsilon u|^2_{L^\infty(I_k; H^1(\Omega))} + |u - \mathcal{P}_{\mathbf{p}^t} u|^2_{L^2(I_k; H^1(\Omega))} \left| (u - \mathcal{P}_{\mathbf{p}^t} u)' \right|_{L^2(I_k; H^1(\Omega))} \right). \nonumber
\end{align}
$$

As for the first term on the right-hand side, polynomial approximation properties (in space) give

$$
|u - \Pi_h^\epsilon u|^2_{L^\infty(I_k; H^1(\Omega))} \le |u - q_{p_x}|^2_{L^\infty(I_k; H^1(\Omega))} \lesssim \frac{h^{2 \min(p_x, s)} \tau_k^{-1}}{p_x^{2s}} (\tau_k \|u\|^2_{L^\infty(I_k; H^{s+1}(\Omega))}).
$$
As for the second term on the right-hand side of (39), we use (20) and (17b), and obtain

$$|u - \mathcal{P}_{\mathbf{p}^t} u|_{L^2(I_k; H^1(\Omega))} |(u - \mathcal{P}_{\mathbf{p}^t} u)'|_{L^2(I_k; H^1(\Omega))} \lesssim \left(\frac{\tau_k}{p_k^t}\right)^{2 \min(s, p_k^t+1)-1} \|u^{(s)}\|_{L^2(I_k; H^1(\Omega))}^2.$$

Combining the above displays and recalling (7) entail

$$\mathbb{J}_8^2 \lesssim \left( \frac{h^{2 \min(p_{\mathbf{x}}, s)} \tau_k^{-1}}{p_{\mathbf{x}}^{2s}} (\tau_k \|u\|_{L^{\infty}(I_k; H^{s+1}(\Omega))}^2) + \left( \frac{\tau_k}{p_k^t} \right)^{2 \min(s, p_k^t+1)-1} \|u^{(s)}\|_{L^2(I_k; H^1(\Omega))}^2 \right). \quad (40)$$

**Error estimates.** The assertion follows combining (31), (32), (33), (34), (35), (36), (38), and (40). □

The a priori estimates in Theorem 2.8 can be simplified under the assumption of elliptic regularity, i.e., assuming the parameter $\alpha$ in (22) to be 1 (this happens for instance if $\Omega$ is convex); under more regularity on the solution; requiring the isotropy of the spatial and time meshes; fixing the polynomial degrees in space and time (in what follows, $p$ denotes the polynomial degree of the scheme).

Notably, the next result contains two estimates: the first one holds true for smooth exact solutions and the rate in the time step has optimal order $p$; the second one holds true for exact solutions with finite total Sobolev regularity indices $s$ smaller than or equal to $p$ and is optimal in terms of the Sobolev scaling.

**Corollary 2.9.** Let $u$ and $U_h$ be the solutions to (2) and (3), and $e_h$ be as in (23). We assume that the initial conditions $u_0$ and $u_1$ in (2) and $u$ are smooth in the sense $s$ is larger than or equal to $p_n^t + 2$ on each $I_n$ for all $n = 1, \dots, N$. Assume that the parameter $\alpha$ in (22) is equal to 1, i.e., elliptic regularity holds for the domain $\Omega$. Let $\tau_n = \tau$ for all $n = 1, \dots, N$, and choose $h = O(\tau)$. For given $p$ in $\mathbb{N}$, we further demand that $p_{\mathbf{x}} = p$ and $p_n^t = p$ for all $n = 1, \dots, N$. Recall that $m = m(\xi_h)$ is defined in (6). Then, the following a priori error estimate is valid:

$$\max_{n=1}^{N} \left( \|e_h'\|_{L^{\infty}(I_n; L^2(\Omega))}^2 + |e_h|_{L^{\infty}(I_n; H^1(\Omega))}^2 \right)^{\frac{1}{2}} \lesssim \tau^p.$$

On the other hand, if $s$ is smaller than or equal to $p$, then the following a priori error estimate is valid:

$$\max_{n=1}^{N} \left( \|e_h'\|_{L^{\infty}(I_n; L^2(\Omega))}^2 + |e_h|_{L^{\infty}(I_n; H^1(\Omega))}^2 \right)^{\frac{1}{2}} \lesssim \tau^{s-\frac{1}{2}}. \quad (41)$$

**Proof.** The proof boils down to using Theorem 2.8 and note that, for the estimates of the terms $\mathbb{J}_3$ and $\mathbb{J}_7$ in (29), we can use

$$\|u^{(p+1)}\|_{L^2(I_k; L^2(\Omega))} \le \tau_k^{\frac{1}{2}} \|u^{(p+1)}\|_{L^\infty(I_k; L^2(\Omega))}.$$

All other $\mathbb{J}_j$ terms, $j = 1, 2, 4, 5, 6, 8$, already gave $O(\tau^p)$ rates. □

**Remark 3.** The error measure in (30) does not involve the jumps at the time nodes of the first time derivative. Some comments about this feature are in order.

* Lemma 2.6 involves an error measure containing the sum of jump terms up to $m$, $m = m(\xi_h)$ as in (6); a modification of the proofs of Theorems 2.7 and 2.8 would lead to estimates for norms of the error containing jumps, but only up to the time node $t_m$ and leading to an optimal convergence rate $O(\tau^{p-\frac{1}{2}})$, which is however suboptimal for the first term in the error measure (30).

* Walkington's strategy [39] involves the use of a special test function, leading to errors measured in $L^\infty$-type norms in time. Such norms are different from $L^2$-type norms in time, which are typically obtained by testing with "more standard" functions and typically come together with jumps at the time nodes.
* Optimal convergence rates for a norm involving jumps in our setting may be derived by using stability estimates obtained by using the test function $W_h = U_h'$. This is shown, for instance, in [39, eq. (4.1)] and would give an optimal convergence rate for the jump terms.

* Even though we did not prove convergence rates for the jumps on the theoretical level, in Section 4 below, we shall investigate their practical behaviour.

# 3 A posteriori error estimates for the semi-discrete in time method

This section is concerned with introducing an error estimator for the semi-discrete in time method (4), and prove fully explicit, reliable a posteriori estimates for the error measured in the $L^\infty(L^2)$ norm under extra assumptions discussed in Section 3.1. In Section 3.2, we introduce a novel reconstruction operator and exhibit its approximation properties, which are instrumental in the a posteriori error estimates given in Section 3.3.

## 3.1 Data assumptions

Throughout, we make the following assumption, which is instrumental in deriving the a posteriori error estimates: for $u$ and $U$ solutions to (2) to (4),

$$ \Delta_x u \in L^1(0, T; L^2(\Omega)) \quad (\text{which implies } \Delta_x U \text{ belongs to } L^1(0, T; L^2(\Omega))). \quad (42) $$

Assumption (42) can be proven for instance under certain conditions on the data, as detailed in the next result.

**Proposition 3.1.** Assume that

$$ \begin{array}{l} \textit{the spatial domain } \Omega \textit{ is convex;} \\ u_0 \textit{ and } u_1 \textit{ belong to } H^2(\Omega) \cap H_0^1(\Omega) \textit{ and } H_0^1(\Omega); \\ \textit{the right-hand side } f \textit{ belongs to } H^1(0, T; L^2(\Omega)). \end{array} \quad (43) $$

Then, the solution $u$ to (2) is such that

$$ u \in L^\infty(0, T; H^2(\Omega)) \cap W^{1,\infty}(0, T; H_0^1(\Omega)) \cap W^{2,\infty}(0, T; L^2(\Omega)) \cap W^{3,\infty}(0, T; H^{-1}(\Omega)). $$

In particular, property (42) holds true.

*Proof.* A proof for domains with sufficiently smooth boundary can be found in [10, Theorem 5, Chapter 7.2] and is based on the Faedo-Galerkin technique, based on taking the limits of expansions into eigenfunctions. In turns, the regularity in space only depends on the regularity of the eigenfunctions of $\Delta_x$, which on convex domains (43) is $H^2(\Omega)$; see, e.g., [20]. □

## 3.2 A reconstruction operator

In view of deriving a posteriori error estimates in Section 3.3 below, we discuss here the properties of a generalization of the reconstruction operator introduced in [27] for parabolic problems; see also [22, 33] for a proof of p-approximation properties of the operator in [27].

Given a Hilbert space $\mathcal{X}$ with inner product $(\cdot, \cdot)_\mathcal{X}$ and induced norm $\|\cdot\|_{\mathcal{X}}$, and $V$ in $C^0(0, T; \mathcal{X})$ with $V|_{I_n}$ in $P_{p_n^t}(I_n; \mathcal{X})$ and $(V^-)'(0, \cdot)$ is an element in $\mathcal{X}$, let $\hat{V}$ be piecewise defined for all $n = 1, \dots, N$ as

$$ \left\{ \begin{array}{l} (\hat{V}'', q_{p_n^t-1})_{L^2(I_n; \mathcal{X})} = (V'', q_{p_n^t-1})_{L^2(I_n; \mathcal{X})} + (\llbracket V \rrbracket (t_{n-1}, \cdot), q_{p_n^t-1}(t_{n-1}, \cdot))_\mathcal{X} \\ \hat{V}(t_{n-1}, \cdot) = V(t_{n-1}, \cdot), \quad \hat{V}'(t_{n-1}, \cdot) = (V^-)'(t_{n-1}, \cdot) \text{ in } \mathcal{X} \quad \forall q_{p_n^t-1} \in P_{p_n^t-1}(I_n; \mathcal{X}). \end{array} \right. \quad (44) $$

The corresponding operator for parabolic problems also allows for spatial mesh changes, see, e.g., [17], a topic that is still open for the wave equation in second order formulation and that we shall investigate in the future.

We begin by proving the following property of $\hat{V}$.
**Proposition 3.2.** The function $\hat{V}$ is a $C^1$ piecewise polynomial in time reconstruction of a $C^0$ piecewise polynomial in time.

*Proof.* The assertion follows using the last condition in (44) and proving that

$$ \hat{V}'(t_n, \cdot) = (V^{-})'(t_n, \cdot), \quad \hat{V}(t_n, \cdot) = V(t_n, \cdot) \quad \text{in } \mathcal{X}. \qquad (45) $$

**Proving the first identity in (45).** Taking $q_{p_n^t-1} = c_X$ to be fixed in $\mathcal{X}$ below and independent of time in the first line of (44) for all $n = 1, \dots, N$, and integrating by parts on both sides lead to

$$ \begin{aligned} & (\hat{V}'(t_n, \cdot) - \hat{V}'(t_{n-1}, \cdot), c_X)_\mathcal{X} \\ & = ((V^{-})'(t_n, \cdot) - \underbrace{(V^{+})'(t_{n-1}, \cdot) + (V^{+})'(t_{n-1}, \cdot)}_{=0} - (V^{-})'(t_{n-1}, \cdot), c_X)_\mathcal{X}. \end{aligned} $$

Using that $\hat{V}'(t_{n-1}, \cdot) = (V^{-})'(t_{n-1}, \cdot)$ in $\mathcal{X}$ and taking $c_X$ equal to $\hat{V}'(t_n, \cdot) - (V^{-})'(t_n, \cdot)$, we deduce the first identity in (45).

**Proving the second identity in (45).** We take $q_{p_n^t-1} = (t - t_{n-1})\tilde{c}_X$ with $\tilde{c}_X$ to be fixed in $\mathcal{X}$ below and independent of time in the first condition of (44) for all $n = 1, \dots, N$, integrate by parts, and get

$$ (\hat{V}'(t_n, \cdot), \tau_n \tilde{c}_X)_\mathcal{X} - (\hat{V}', \tilde{c}_X)_{L^2(I_n; \mathcal{X})} = ((V^{-})'(t_n, \cdot), \tau_n \tilde{c}_X)_\mathcal{X} - (V', \tilde{c}_X)_{L^2(I_n; \mathcal{X})}. $$

Using $\hat{V}'(t_n, \cdot) = (V^{-})'(t_n, \cdot)$ entails

$$ (\hat{V}', \tilde{c}_X)_{L^2(I_n; \mathcal{X})} = (V', \tilde{c}_X)_{L^2(I_n; \mathcal{X})}. $$

Integrating by parts again gives

$$ (\hat{V}(t_n, \cdot) - \hat{V}(t_{n-1}, \cdot), \tilde{c}_X)_\mathcal{X} = (V(t_n, \cdot) - V(t_{n-1}, \cdot), \tilde{c}_X)_\mathcal{X}. $$

Using that $\hat{V}(t_{n-1}, \cdot) = V(t_{n-1}, \cdot)$ in $\mathcal{X}$, which is the first initial condition in (44), and choosing $\tilde{c}_X$ equal to $\hat{V}(t_n, \cdot) - V(t_n, \cdot)$ yield the second identity in (45). $\square$

We have additional properties on the operator $\hat{V}$.

**Lemma 3.3.** Consider $V$ in $C^0(0, T; \mathcal{X})$ with $V|_{I_n}$ in $\mathbb{P}_{p_n^t}(I_n; \mathcal{X})$ and $\hat{V}$ as in (44). For all $n = 1, \dots, N$, the following identities hold true:

$$ \|(V - \hat{V})'\|_{L^2(I_n; \mathcal{X})}^2 = \tau_n c_1 (p_n^t)^2 \|(V')'\|_{(I_n; \mathcal{X})}^2, \qquad (46a) $$

$$ \|(V - \hat{V})'\|_{L^\infty(I_n; \mathcal{X})}^2 = \|(V')'\|_{(I_n; \mathcal{X})}^2, \qquad (46b) $$

$$ \|V - \hat{V}\|_{L^2(I_n; \mathcal{X})}^2 \le \tau_n^3 c_2 (p_n^t)^2 \|(V')'\|_{(I_n; \mathcal{X})}^2, \qquad (46c) $$

where

$$ c_1(p_n^t)^2 := \frac{p_n^t}{(2p_n^t - 1)(2p_n^t + 1)}, \qquad c_1(p_n^t) \approx (p_n^t)^{-\frac{1}{2}}, \qquad (47a) $$

$$ c_2(p_n^t)^2 := \begin{cases} \frac{p_n^t}{\frac{1}{4}(p_n^t - 2)(p_n^t - 1)(2p_n^t - 1)(2p_n^t + 1)} & \text{if } p_n^t \ge 3, \\ \frac{2}{15\pi^2} & \text{if } p_n^t = 2, \end{cases} \qquad c_2(p_n^t) \approx (p_n^t)^{-\frac{3}{2}}. \qquad (47b) $$

*Proof.* Identity (46a) is proven in [33, Theorem 2]. Identity (46b) is essentially proven in [22, Lemma 1]. As for inequality (46c), we first observe that definition (44), an integration by parts, the first identity in (45), and the smoothness of $\hat{V}$ imply, for all $n = 1, \dots, n$,

$$ \begin{aligned} & - ((V - \hat{V})'(t_{n-1}, \cdot), q_{p_n^t-1}(t_{n-1}, \cdot)) \chi \stackrel{(44)}{=} ((V - \hat{V})'', q_{p_n^t-1})_{L^2(I_n; \mathcal{X})} \\ & \stackrel{(IBP)}{=} - ((V - \hat{V})', q'_{p_n^t-1})_{L^2(I_n; \mathcal{X})} \\ & \qquad + ((V - \hat{V})'(t_n, \cdot), q_{p_n^t-1}(t_n, \cdot)) \chi - ((V - \hat{V})'(t_{n-1}, \cdot), q_{p_n^t-1}(t_{n-1}, \cdot)) \chi \\ & \stackrel{(45)}{=} - ((V - \hat{V})', q'_{p_n^t-1})_{L^2(I_n; \mathcal{X})} - ((V')'(t_{n-1}, \cdot), q_{p_n^t-1}(t_{n-1}, \cdot)) \chi \quad \forall q_{p_n^t-1} \in \mathbb{P}_{p_n^t-1}(I_n; \mathcal{X}). \end{aligned} $$
**Proving (46c) for $p_n^t$ larger than 2.** A further integration by parts and the smoothness of $V$ and $\hat{V}$ at the endpoints of each $I_n$ entail

$$ (V - \hat{V}, q''_{p_n^t-1})_{L^2(I_n; \mathcal{X})} = 0 \quad \forall q_{p_n^t-1} \in \mathbb{P}_{p_n^t-1}(I_n; \mathcal{X}). \quad (48) $$

For all $t$ in $I_n$, we pick

$$ q_{p_n^t-1}(t, \cdot) := \int_{t_{n-1}}^t \int_{t_{n-1}}^s \Pi_{\mathbf{p}^t-3}^0 (V - \hat{V})(r, \cdot) dr ds \quad \text{in } \mathcal{X} $$

and deduce

$$ (V - \hat{V}, \Pi_{\mathbf{p}^t-3}^0 (V - \hat{V}))_{L^2(I_n; \mathcal{X})} = 0. \quad (49) $$

Using [34, Theorem 3.11] entails

$$ \|(I - \Pi_{\mathbf{p}^t-3}^0)(V - \hat{V})\|_{L^2(I_n; \mathcal{X})} \le \frac{1}{2} \frac{\tau_n}{\sqrt{(p_n^t - 2)(p_n^t - 1)}} \|(V - \hat{V})'\|_{L^2(I_n; \mathcal{X})}. \quad (50) $$

Approximation properties as in [34, Theorem 3.11] imply

$$ \begin{aligned} \|V - \hat{V}\|_{L^2(I_n; \mathcal{X})}^2 &= (V - \hat{V}, V - \hat{V})_{L^2(I_n; \mathcal{X})} \\ &\stackrel{(49)}{=} (V - \hat{V}, (I - \Pi_{\mathbf{p}^t-3}^0)(V - \hat{V}))_{L^2(I_n; \mathcal{X})} \\ &\stackrel{(50)}{\le} \frac{1}{2} \frac{\tau_n}{\sqrt{(p_n^t - 2)(p_n^t - 1)}} \|V - \hat{V}\|_{L^2(I_n; \mathcal{X})} \|(V - \hat{V})'\|_{L^2(I_n; \mathcal{X})}. \end{aligned} $$

The first line in inequality (46c) follows recalling identity (46a).

**Proving (46c) for $p_n^t$ equal to 2.** The last display above modifies as follows: the 1D Poincaré inequality in time implies

$$ \|V - \hat{V}\|_{L^2(I_n; \mathcal{X})}^2 \le \frac{\tau_n^2}{\pi^2} \|(V - \hat{V})'\|_{L^2(I_n; \mathcal{X})}^2. $$

The second line in inequality (46c) follows recalling identity (46a) with $p_n^t = 2$. $\square$

## 3.3 An error estimator for the semi-discrete in time method

We construct a computable error estimator and prove fully explicit, reliable a posteriori error estimates. We define $\xi$ and $m$ such that

$$ \|(u - \hat{U})(\xi, \cdot)\|_{0,\Omega} = \|u - \hat{U}\|_{L^{\infty}(0,T;L^2(\Omega))} \quad \text{with } \xi \in I_m. \quad (51) $$

We introduce

$$ v_B(t, \cdot) = \int_t^\xi (u - \hat{U})(s, \cdot) ds \quad \text{in } \mathcal{X}. \quad (52) $$

The function $v_B$ belongs to $C^1(0, T; L^2(\Omega))$. In definition (44), we set

$$ \hat{U}'(0, \cdot) = u_1(\cdot) \quad \text{in } L^2(\Omega). \quad (53) $$

We state an auxiliary technical result.

**Lemma 3.4.** Given $v_B$ and $\xi$ as in (52) and (51), the following estimates hold true:

$$ \inf_{q_{p_n^t} \in \mathbb{P}_{p_n^t}(I_n; L^2(\Omega))} \|v_B - q_{p_n^t}\|_{L^\infty(I_n; L^2(\Omega))} \le \frac{\tau_n}{\pi} c_3(p_n^t) \|v_B'\|_{L^\infty(I_n; L^2(\Omega))}, \quad (54a) $$

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

where we have set, also for future convenience,

$$ c_3(p_n^t) := \begin{cases} \sqrt{\pi} & \text{if } p_n^t = 0, 1, 2 \\ \frac{1}{p_n^t - 2} & \text{if } p_n^t \ge 3, \end{cases} \quad c_3(p_n^t) \approx (p_n^t)^{-1}. \quad (55) $$
*Proof.* Inequality (54a) in the case $p_n^t$ larger than 2 is the Bochner version of [37, Theorem 7.2] up to a scaling argument. Inequality (54a) in the case $p_n^t$ equal to 0, 1, and 2 follows from [24, eq. (1.3)] and the 1D Poincaré inequality in time. Inequality (54b) follows from the following computations:

$$
\begin{align*}
\|v_B\|_{L^\infty(t_{m-1}, \xi; L^2(\Omega))} &= \sup_{t \in (t_{m-1}, \xi)} \|v_B(t, \cdot)\|_{L^2(\Omega)} = \sup_{t \in (t_{m-1}, \xi)} \left( \int_{\Omega} |v_B(t, \mathbf{x})|^2 d\mathbf{x} \right)^{\frac{1}{2}} \\
&= \sup_{t \in (t_{m-1}, \xi)} \left( \int_{\Omega} |\int_t^\xi (u - \hat{U})(s, \mathbf{x}) ds|^2 d\mathbf{x} \right)^{\frac{1}{2}} \le \sup_{t \in (t_{m-1}, \xi)} \left( \int_{\Omega} |\xi - t| \int_t^\xi |(u - \hat{U})(s, \mathbf{x})|^2 ds \, d\mathbf{x} \right)^{\frac{1}{2}} \\
&\le \tau_m^{\frac{1}{2}} \sup_{t \in (t_{m-1}, \xi)} \left( \int_{\Omega} \int_t^\xi |(u - \hat{U})(s, \mathbf{x})|^2 ds \, d\mathbf{x} \right)^{\frac{1}{2}} = \tau_m^{\frac{1}{2}} \sup_{t \in (t_{m-1}, \xi)} \left( \int_{\Omega} \int_t^\xi |v_B'(s, \mathbf{x})|^2 ds \, d\mathbf{x} \right)^{\frac{1}{2}} \\
&\le \tau_m \|v_B'\|_{L^\infty(t_{m-1}, \xi; L^2(\Omega))}.
\end{align*}
$$

□

**Proposition 3.5.** Let assumption (42) hold true. Consider $u$ and $U$ the solutions to (2) and (4), and the operator $\hat{U}$ in (11). Given $v_B$, $\xi$ and $m$, $c_2(p_n^t)$, and $c_3(p_n^t)$ as in (52), (51), (47b), and (55), we have

$$
\begin{align}
& \|u - \hat{U}\|_{L^{\infty}(0,T;L^2(\Omega))} \nonumber \\
& \le 2 \left( \begin{aligned}
& \left( \sum_{n=1}^{m-1} \frac{\tau_n}{\pi} c_3(p_n^t - 1) \|f - \Pi_{\mathbf{p}^t-1}^0 f\|_{L^1(I_n;L^2(\Omega))} + \tau_m \|f - \Pi_{\mathbf{p}^t-1}^0 f\|_{L^1(I_m;L^2(\Omega))} \right) \\
& + \left( \sum_{n=1}^{m-1} \frac{\tau_n}{\pi} c_3(p_n^t - 1) \|\Delta_{\mathbf{x}}(U - \Pi_{\mathbf{p}^t-1}^0 U)\|_{L^1(I_n;L^2(\Omega))} + \tau_m \|\Delta_{\mathbf{x}}(U - \Pi_{\mathbf{p}^t-1}^0 U)\|_{L^1(I_m;L^2(\Omega))} \right) \\
& + \left( \sum_{n=1}^{m-1} \frac{\tau_n^3}{\pi} c_2(p_n^t) c_4(p_n^t - 3) \|[[\Delta_{\mathbf{x}}U']](t_{n-1}, \cdot)\|_{0,\Omega} + \tau_m^3 c_2(p_m^t) \|[[\Delta_{\mathbf{x}}U']](t_{m-1}, \cdot)\|_{0,\Omega} \right)
\end{aligned} \right), \tag{56}
\end{align}
$$

with $c_4(p_n^t - 3)$ defined as

$$
c_4(p_n^t - 3) := \begin{cases} \pi^{\frac{|t_m - t_{n-1}|}{\tau_n}} & \text{if } p_n^t = 2 \\ c_3(p_n^t - 3) & \text{if } p_n^t \ge 3. \end{cases} \quad (57)
$$

*Proof.* Let $\hat{U}$ be the reconstruction operator in (44). Using (4), we have the identity

$$
\begin{align*}
& (\hat{U}''(V), V)_{L^2(I_n; L^2(\Omega))} - (\Delta_{\mathbf{x}}U, V)_{L^2(I_n; L^2(\Omega))} \\
& = (f, V)_{L^2(I_n; L^2(\Omega))} \quad \forall V \in \mathbb{P}_{p_n^t-1}^0(I_n; H_0^1(\Omega)) \quad \forall n = 1, \dots, N.
\end{align*}
$$

Recall the time semi-discrete right-hand side in (5). Due to (42), for all $n = 1, \dots, N$, we have the following identity in $L^2(\Omega)$ inside each time interval:

$$
\hat{U}'' - \Pi_{\mathbf{p}^t-1}^0 \Delta_{\mathbf{x}} U = \Pi_{\mathbf{p}^t-1}^0 f. \quad (58)
$$

Using $\Delta_{\mathbf{x}}\Pi_{\mathbf{p}^t-1}^0 U = \Pi_{\mathbf{p}^t-1}^0 \Delta_{\mathbf{x}} U$, and subtracting (58) to (1), we get the following identity in $L^2(I_n; L^2(\Omega))$:

$$
(u - \hat{U})'' - \Delta_{\mathbf{x}}(u - \hat{U}) = (f - \Pi_{\mathbf{p}^t-1}^0 f) + \Delta_{\mathbf{x}}(\hat{U} - \Pi_{\mathbf{p}^t-1}^0 U). \quad (59)
$$

Next, we derive an error equation testing (59) with a particular function mimicking that proposed in [4, eq. (3.7)], i.e., the function in (52). It is immediate to check

$$
-v_B'(t, \cdot) = (u - \hat{U})(t, \cdot) \tag{60}
$$

and

$$
v_B(\xi, \cdot) = 0. \tag{61}
$$
Multiplying (59) with the function in (52), integrating in space, and integrating by parts, we arrive at the following identity in $L^2(0, T)$:

$$((u - \hat{U})'', v_B)_{0,\Omega} + a(u - \hat{U}, v_B) = (f - \Pi_{\mathbf{p}^t-1}^0 f, v_B)_{0,\Omega} + (\Delta_{\mathbf{x}}(\hat{U} - \Pi_{\mathbf{p}^t-1}^0 U), v_B)_{0,\Omega}.$$

Using (60), we readily deduce the following identity in $L^2(0, T)$:

$$((u - \hat{U})'', v_B)_{0,\Omega} - a(v_B', v_B) = (f - \Pi_{\mathbf{p}^t-1}^0 f, v_B)_{0,\Omega} + (\Delta_{\mathbf{x}}(\hat{U} - \Pi_{\mathbf{p}^t-1}^0 U), v_B)_{0,\Omega}.$$

The following trivial but crucial identity holds true: given sufficiently smooth in time functions $\aleph$ and $\backslash$, we have

$$(\aleph'', \backslash)_{0,\Omega} = (\aleph', \backslash)'_{0,\Omega} - (\aleph', \backslash)'_{0,\Omega} \quad \text{in } L^2(0, T).$$

This identity and (60) imply

$$((u - \hat{U})', v_B)'_{0,\Omega} + ((u - \hat{U})', u - \hat{U})_{0,\Omega} - a(v_B', v_B) = (f - \Pi_{\mathbf{p}^t-1}^0 f, v_B)_{0,\Omega} + (\Delta_{\mathbf{x}}(\hat{U} - \Pi_{\mathbf{p}^t-1}^0 U), v_B)_{0,\Omega} \quad \text{in } L^1(0, T).$$

Equivalently, we write

$$((u - \hat{U})', v_B)'_{0,\Omega} + \frac{1}{2}(\|u - \hat{U}\|_{0,\Omega}^2)' - \frac{1}{2}a(v_B, v_B)' = (f - \Pi_{\mathbf{p}^t-1}^0 f, v_B)_{0,\Omega} + (\Delta_{\mathbf{x}}(\hat{U} - \Pi_{\mathbf{p}^t-1}^0 U), v_B)_{0,\Omega} \quad \text{in } L^1(0, T).$$

We integrate in time over $(0, \xi)$, $\xi$ as in the choice of the test function in (52). We arrive at

$$\frac{1}{2} \|(u - \hat{U})(\xi, \cdot)\|_{0,\Omega}^2 - \frac{1}{2} a(v_B(\xi, \cdot), v_B(\xi, \cdot)) = \frac{1}{2} \|(u - \hat{U})(0, \cdot)\|_{0,\Omega}^2 - \frac{1}{2} a(v_B(0, \cdot), v_B(0, \cdot)) \\ + \int_0^\xi \left( (f - \Pi_{\mathbf{p}^t-1}^0 f, v_B)_{0,\Omega} + (\Delta_{\mathbf{x}}(\hat{U} - \Pi_{\mathbf{p}^t-1}^0 U), v_B)_{0,\Omega} - ((u - \hat{U})', v_B)'_{0,\Omega} \right) dt.$$

Using (44) and (53), we have $(u - \hat{U})'(0, \cdot) = 0$, whence we write

$$\begin{aligned} & - \int_0^\xi ((u - \hat{U})', v_B)'_{0,\Omega} dt \\ &= -((u - \hat{U})'(\xi, \cdot), v_B(\xi, \cdot))_{0,\Omega} + ((u - \hat{U})'(0, \cdot), v_B(0, \cdot))_{0,\Omega} \\ &\stackrel{(61)}{=} ((u - \hat{U})'(0, \cdot), v_B(0, \cdot))_{0,\Omega} = 0. \end{aligned}$$

Combining the two above displays again with (61), adding and subtracting $U$ in the last term on the right-hand side, and using that $(u - \hat{U})(0, \cdot) = 0$ yield

$$\frac{1}{2} \|(u - \hat{U})(\xi, \cdot)\|_{0,\Omega}^2 + \frac{1}{2} a(v_B(0, \cdot), v_B(0, \cdot)) = \int_0^\xi \left( (f - \Pi_{\mathbf{p}^t-1}^0 f, v_B)_{0,\Omega} + (\Delta_{\mathbf{x}}(U - \Pi_{\mathbf{p}^t-1}^0 U), v_B)_{0,\Omega} + (\Delta_{\mathbf{x}}(\hat{U} - U), v_B)_{0,\Omega} \right) dt.$$

Let $q_{p_n^t-1}$ and $q_{p_n^t-3}$ realize (54a) of degree $p_n^t - 1$ and $p_n^t - 3$ (with $q_{p_n^t-3} = 0$ for $p_n^t = 2$), respectively. Using the properties of $\Pi_{\mathbf{p}^t-1}^0$, the orthogonal property (48), and the choice of $\xi$ in (51) (including the fact that $\xi$ lies in the time interval $I_m$), we end up with

$$\begin{aligned} & \frac{1}{2} \|u - \hat{U}\|_{L^\infty(0,T;L^2(\Omega))}^2 + \frac{1}{2} a(v_B(0, \cdot), v_B(0, \cdot)) \\ &\le \sum_{n=1}^{m-1} \int_{I_n} (f - \Pi_{\mathbf{p}^t-1}^0 f, v_B - q_{p_n^t-1})_{0,\Omega} dt + \int_{t_{m-1}}^\xi (f - \Pi_{\mathbf{p}^t-1}^0 f, v_B)_{0,\Omega} dt \\ &\quad + \sum_{n=1}^{m-1} \int_{I_n} (\Delta_{\mathbf{x}}(U - \Pi_{\mathbf{p}^t-1}^0 U), v_B - q_{p_n^t-1})_{0,\Omega} dt + \int_{t_{m-1}}^\xi (\Delta_{\mathbf{x}}(U - \Pi_{\mathbf{p}^t-1}^0 U), v_B)_{0,\Omega} dt \\ &\quad + \sum_{n=1}^{m-1} \int_{I_n} (\Delta_{\mathbf{x}}(\hat{U} - U), v_B - q_{p_n^t-3})_{0,\Omega} dt + \int_{t_{m-1}}^\xi (\Delta_{\mathbf{x}}(\hat{U} - U), v_B)_{0,\Omega} dt = \sum_{i=1}^6 T_i. \end{aligned} \quad (62)$$
As for the term $T_1$, we use Hölder's inequality, inequality (54a) (with constant $c_3(p_n^t - 1)$), and the fact that $\|v_B'\|_{L^\infty(I_n; L^2(\Omega))}$ is smaller than $\|v_B'\|_{L^\infty(0, T; L^2(\Omega))}$, and end up with

$$ 
\begin{aligned}
T_1 &\le \sum_{n=1}^{m-1} \|f - \Pi_{\mathbf{p}^t-1}^0 f\|_{L^1(I_n; L^2(\Omega))} \|v_B - q_{p_n^t-1}\|_{L^\infty(I_n; L^2(\Omega))} \\
&\le \sum_{n=1}^{m-1} \frac{\tau_n}{\pi} c_3(p_n^t - 1) \|f - \Pi_{\mathbf{p}^t-1}^0 f\|_{L^1(I_n; L^2(\Omega))} \|v_B'\|_{L^\infty(I_n; L^2(\Omega))} \quad (63) \\
&\le \left( \sum_{n=1}^{m-1} \frac{\tau_n}{\pi} c_3(p_n^t - 1) \|f - \Pi_{\mathbf{p}^t-1}^0 f\|_{L^1(I_n; L^2(\Omega))} \right) \|v_B'\|_{L^\infty(0, T; L^2(\Omega))}.
\end{aligned}
$$

Similarly, using (54b) and the fact that $\|v_B'\|_{L^\infty(t_{m-1}, \xi; L^2(\Omega))}$ is smaller than $\|v_B'\|_{L^\infty(0, T; L^2(\Omega))}$, we also infer the following bound on the term $T_2$:

$$ 
\begin{aligned}
T_2 &\le \|f - \Pi_{\mathbf{p}^t-1}^0 f\|_{L^1(t_{m-1}, \xi; L^2(\Omega))} \|v_B\|_{L^\infty(t_{m-1}, \xi; L^2(\Omega))} \\
&\le \tau_m \|f - \Pi_{\mathbf{p}^t-1}^0 f\|_{L^1(I_m; L^2(\Omega))} \|v_B'\|_{L^\infty(t_{m-1}, \xi; L^2(\Omega))} \quad (64) \\
&\le \tau_m \|f - \Pi_{\mathbf{p}^t-1}^0 f\|_{L^1(I_m; L^2(\Omega))} \|v_B'\|_{L^\infty(0, T; L^2(\Omega))}.
\end{aligned}
$$

Next, we bound the term $T_3$ using again (54a) with constant $c_3(p_n^t - 1)$:

$$ 
\begin{aligned}
T_3 &\le \sum_{n=1}^{m-1} \|\Delta_{\mathbf{x}}(U - \Pi_{\mathbf{p}^t-1}^0 U)\|_{L^1(I_n; L^2(\Omega))} \|v_B - q_{p_n^t-1}\|_{L^\infty(I_n; L^2(\Omega))} \\
&\le \sum_{n=1}^{m-1} \frac{\tau_n}{\pi} c_3(p_n^t - 1) \|\Delta_{\mathbf{x}}(U - \Pi_{\mathbf{p}^t-1}^0 U)\|_{L^1(I_n; L^2(\Omega))} \|v_B'\|_{L^\infty(I_n; L^2(\Omega))} \quad (65) \\
&\le \left( \sum_{n=1}^{m-1} \frac{\tau_n}{\pi} c_3(p_n^t - 1) \|\Delta_{\mathbf{x}}(U - \Pi_{\mathbf{p}^t-1}^0 U)\|_{L^1(I_n; L^2(\Omega))} \right) \|v_B'\|_{L^\infty(0, T; L^2(\Omega))}.
\end{aligned}
$$

We bound the term $T_4$ based on (54b) and proceeding as in the bound of $T_3$:

$$ 
\begin{aligned}
T_4 &\le \|\Delta_{\mathbf{x}}(U - \Pi_{\mathbf{p}^t-1}^0 U)\|_{L^1(t_{m-1}, \xi; L^2(\Omega))} \|v_B\|_{L^\infty(t_{m-1}, \xi; L^2(\Omega))} \\
&\le \tau_m \|\Delta_{\mathbf{x}}(U - \Pi_{\mathbf{p}^t-1}^0 U)\|_{L^1(t_{m-1}, \xi; L^2(\Omega))} \|v_B'\|_{L^\infty(t_{m-1}, \xi; L^2(\Omega))} \quad (66) \\
&= \tau_m \|\Delta_{\mathbf{x}}(U - \Pi_{\mathbf{p}^t-1}^0 U)\|_{L^1(I_m; L^2(\Omega))} \|v_B'\|_{L^\infty(0, T; L^2(\Omega))}.
\end{aligned}
$$

Using estimates (54a) with constant $c_3(p_n^t - 3)$ and (46c) with constant $c_2(p_n^t)$, we show an upper bound on the term $T_5$ for $p_n^t \ge 3$:

$$ 
\begin{aligned}
T_5 &\le \sum_{n=1}^{m-1} \tau_n^{\frac{1}{2}} \|\Delta_{\mathbf{x}}(\hat{U} - U)\|_{L^2(I_n; L^2(\Omega))} \|v_B - q_{p_n^t-3}\|_{L^\infty(I_n; L^2(\Omega))} \\
&\le \sum_{n=1}^{m-1} \frac{\tau_n^{\frac{3}{2}}}{\pi} c_3(p_n^t - 3) \|\Delta_{\mathbf{x}}(\hat{U} - U)\|_{L^2(I_n; L^2(\Omega))} \|v_B'\|_{L^\infty(I_n; L^2(\Omega))} \quad (67) \\
&\le \left( \sum_{n=1}^{m-1} \frac{\tau_n^3}{\pi} c_2(p_n^t) c_3(p_n^t - 3) \|\Delta_{\mathbf{x}}U'\|_{(t_{n-1}, \cdot)} \right) \|v_B'\|_{L^\infty(0, T; L^2(\Omega))}.
\end{aligned}
$$

For $p_n^t = 2$, using the definition of $v_B$ and $\xi$ in (52) and (51), and the fact that $|\xi - t_{n-1}|$ is larger than $|\xi - t_n|$, we have

$$ \|v_B\|_{L^\infty(I_n; L^2(\Omega))} \le |\xi - t_{n-1}| \|v_B'\|_{L^\infty(I_n; L^2(\Omega))} \le |t_m - t_{n-1}| \|v_B'\|_{L^\infty(I_n; L^2(\Omega))}. $$
Using the definition of $c_4(p_n^t - 3)$ in (57), we have the following bound on the term $T_5$ for all $p_n^t$ larger than or equal to 2:

$$
\begin{aligned}
T_5 &\le \sum_{n=1}^{m-1} \tau_n^{\frac{1}{2}} \|\Delta_x(\hat{U} - U)\|_{L^2(I_n; L^2(\Omega))} \|v_B\|_{L^\infty(I_n; L^2(\Omega))} \\
&\le \sum_{n=1}^{m-1} |t_m - t_{n-1}| \tau_n^{\frac{1}{2}} \|\Delta_x(\hat{U} - U)\|_{L^2(I_n; L^2(\Omega))} \|v_B'\|_{L^\infty(I_n; L^2(\Omega))} \\
&\le \left( \sum_{n=1}^{m-1} |t_m - t_{n-1}| \tau_n^2 c_2(p_n^t) \|\llbracket \Delta_x U' \rrbracket (t_{n-1}, \cdot)\|_{0,\Omega} \right) \|v_B'\|_{L^\infty(0,T; L^2(\Omega))}.
\end{aligned}
$$

We bound the term $T_6$ using estimates (46c) and (54b):

$$
\begin{aligned}
T_6 &\le \tau_m^{\frac{1}{2}} \|\Delta_x(\hat{U} - U)\|_{L^2(t_{m-1}, \xi; L^2(\Omega))} \|v_B\|_{L^\infty(t_{m-1}, \xi; L^2(\Omega))} \\
&\le \tau_m^{\frac{3}{2}} \|\Delta_x(\hat{U} - U)\|_{L^2(I_m; L^2(\Omega))} \|v_B'\|_{L^\infty(0, T; L^2(\Omega))} \\
&\le \tau_m^3 c_2(p_m^t) \|\llbracket \Delta_x U' \rrbracket (t_{m-1}, \cdot)\|_{0, \Omega} \|v_B'\|_{L^\infty(0, T; L^2(\Omega))}.
\end{aligned}
\tag{68}
$$

Finally, using the fact that $\|v_B'\|_{L^\infty(0,T;L^2(\Omega))}$ is equal to $\|u - \hat{U}\|_{L^\infty(0,T;L^2(\Omega))}$, the assertion follows inserting (63), (64), (65), (66), (67), and (68) in (62). $\square$

Remark 4. With the notation as in the proof of Proposition 3.5, for the case $p_n^t = 2$, the parameter $c_4(p_n^t - 3)$ in (57) scales as $\tau_n^{-1}$, which reduces by one order the convergence rate of the term $T_5$ compared to the case $p_n^t \ge 3$. On the other hand, the convergence rate for the term $T_5$ for $p_n^t \ge 3$ is one order higher than that of the error $\|u - U\|_{L^\infty(0,T;L^2(\Omega))}$. Therefore, for $p_n^t = 2$, the term $T_5$ converges with the same rate as that of the error $\|u - U\|_{L^\infty(0,T;L^2(\Omega))}$.

The error measure on the left-hand side of (56) contains the reconstruction operator $\hat{U}$, which we do not want to compute in practice. For this reason, we elaborate further. Recall that $m$ is as in (51). Consider the time steps and time polynomial degree distributions $\tau$ and $p^t$ as discussed in Section 1; $c_1(p_n^t)$, $c_2(p_n^t)$, $c_3(p_n^t)$, and $c_4(p_n^t)$ as defined in (47a), (47b), (55), and (57). Define

$$ \eta_1 := \max_{n=1,...,N} \tau_n (c_1(p_n^t) c_2(p_n^t))^{\frac{1}{2}} \|\llbracket U' \rrbracket (t_{n-1}, \cdot)\|_{0,\Omega}. $$

Given $m$ as in (51), for all $n = 1, \dots, m-1$, we also define

$$
\eta_{2,n} := \begin{cases} \frac{2}{\pi} (\tau_n c_3(p_n^t - 1) \|\Delta_x(U - \Pi_{\mathbf{p}^t-1}^0 U)\|_{L^1(I_n; L^2(\Omega))} \\ + \tau_n^3 c_2(p_n^t) c_4(p_n^t - 3) \|\llbracket \Delta_x U' \rrbracket (t_{n-1}, \cdot)\|_{0,\Omega}) & \text{if } n = 1, \dots, m-1 \\ 2(\tau_m \|\Delta_x(U - \Pi_{\mathbf{p}^t-1}^0 U)\|_{L^1(I_m; L^2(\Omega))} + c_2(p_n^t) \tau_m^3 \|\llbracket \Delta_x U' \rrbracket (t_{m-1}, \cdot)\|_{0,\Omega}) & \text{if } n = m, \end{cases}
$$

and the data oscillation terms

$$
\mathrm{osc}_n(f) := \begin{cases} \frac{2\tau_n}{\pi} c_3(p_n^t - 1) \|f - \Pi_{\mathbf{p}^t-1}^0 f\|_{L^1(I_n; L^2(\Omega))} & \text{if } n = 1, \dots, m-1 \\ 2\tau_m \|f - \Pi_{\mathbf{p}^t-1}^0 f\|_{L^1(I_m; L^2(\Omega))} & \text{if } n = m. \end{cases}
$$

With this at hand, for $m$ as in (51), we introduce

$$ \eta_2 := \sum_{n=1}^{m} \eta_{2,n}, \qquad \mathrm{osc}(f) := \sum_{n=1}^{m} \mathrm{osc}_n(f). \tag{69} $$

We are now in a position to derive a fully explicit, reliable a posteriori upper bound for the $L^\infty(L^2)$ error with respect to the error estimator

$$ \eta := \eta_1 + \eta_2. \tag{70} $$
**Proposition 3.6.** Let $u$ and $U$ be the solutions to (2) and (4); $\eta$ as in (70), and $\mathrm{osc}(f)$ as in (69). Then, the following bound holds true:

$$ \|u - U\|_{L^{\infty}(0,T;L^2(\Omega))} \leq \eta + \mathrm{osc}(f). \qquad (71) $$

*Proof.* The triangle inequality implies

$$ \|u - U\|_{L^{\infty}(0,T;L^2(\Omega))} \leq \|u - \hat{U}\|_{L^{\infty}(0,T;L^2(\Omega))} + \|U - \hat{U}\|_{L^{\infty}(0,T;L^2(\Omega))}. $$

Next, we pick $\tilde{m} = 1, \dots, N$ and $\tilde{\xi}$ such that

$$ \|(U - \hat{U})(\tilde{\xi}, \cdot)\|_{0,\Omega} = \|U - \hat{U}\|_{L^{\infty}(0,T;L^2(\Omega))} \quad \text{with } \tilde{\xi} \in I_{\tilde{m}}. $$

An upper bound for the first term in the maximum is a consequence of (56). As for the second term, we use the (1D in time) Sobolev embedding $H^1 \hookrightarrow L^\infty$, see, e.g., [24, eq. (1.3)], and observe

$$ \begin{aligned} \|U - \hat{U}\|_{L^{\infty}(0,T;L^2(\Omega))}^2 &\leq \|U - \hat{U}\|_{L^2(I_{\tilde{m}};L^2(\Omega))} \|(U - \hat{U})'\|_{L^2(I_{\tilde{m}};L^2(\Omega))} \\ &\leq \tau_{\tilde{m}}^2 c_1(p_{\tilde{m}}^t)c_2(p_{\tilde{m}}^t) \|[U']\|_{(t_{\tilde{m}-1}, \cdot)}\|_{0,\Omega}^2. \end{aligned} $$

As $\tilde{m}$ and $\tilde{\xi}$ cannot be determined a priori, we further elaborate the above bound as follows:

$$ \|U - \hat{U}\|_{L^{\infty}(0,T;L^2(\Omega))}^2 \leq \max_{n=1,\dots,N} \tau_n^2 c_1(p_n^t)c_2(p_n^t) \|[U']\|_{(t_{n-1}, \cdot)}\|_{0,\Omega}^2. $$

The assertion follows combining the above displays.

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

*Remark 5.* Property (42) is crucial in deriving bounds (71). In principle, we may also assume the weaker requirement that $\Delta_x U$ belongs to $L^2(0, T; H^{-1}(\Omega))$ but end up with an error estimator involving a negative norm, which is harder to realize in practice than the $L^2(0, T; L^2(\Omega))$ norm.

# 4 Numerical experiments

We assess numerically the a priori and a posteriori error estimates proven in Theorem 2.8 and Corollary 2.9, and Proposition 3.6. The numerical experiments are conducted with the `Gridap.jl` library [38] in the `Julia` programming language.

**Spatial and time discretization.** In what follows, we always consider the spatial domain $\Omega = (-1, 1)^2$ partitioned into uniform tensor-product meshes; we use tensor Lagrangian (equidistributed) nodal basis functions. As for the time discretization, we fix $T = 1$ and take the Lagrangian basis functions in time for simplicity.

**Error measures.** We are interested in investigating different error measures for the a priori and a posteriori error estimates.

Given $u$ and $U_h$ the solutions to (2) and (3), and $e_h := u - U_h$, we consider the two error measures appearing on the left-hand side of (30) separately:

$$ \max_{n=1}^{N} \|e_h'\|_{L^{\infty}(I_n; L^2(\Omega))}, \qquad \max_{n=1}^{N} |e_h|_{L^{\infty}(I_n; H^1(\Omega))}. \qquad (72) $$

We also consider the errors

$$ |e_h|_{L^2(0,T;H^1(\Omega))}, \qquad \|e_h'\|_{L^2(0,T;L^2(\Omega))}, \qquad \|e_h\|_{L^\infty(0,T;L^2(\Omega))}, \qquad (73) $$

and the jump error

$$ \left( \sum_{n=1}^{N} \| [e_h'] (t_{n-1}, \cdot) \|_{0,\Omega}^2 \right)^{\frac{1}{2}}. \qquad (74) $$
The L²-type norms in time are computed with Gauß-Legendre quadrature formulas of order 2p<sub>n</sub><sup>t</sup>+3; the L<sup>∞</sup>-type norms in time are evaluated at 2p<sub>n</sub><sup>t</sup> + 3 equally distributed nodes in time for each time interval; the spatial norms are computed by using tensor product Gauß-Legendre quadrature formulas of order 2p<sub>x</sub> + 3.

**Imposing Dirichlet boundary conditions.** We consider homogeneous Dirichlet boundary conditions, which are imposed strongly at the boundary degrees of freedom. The case of inhomogeneous boundary conditions is discussed in [39, Sections 3 and 5.3].

## 4.1 Uniform refinements

We investigate the convergence rate of the fully-discrete method (3) under uniform time step and polynomial degree in time refinements in various norms for three tests.

The initial and boundary conditions, and the right-hand side of (1) are computed accordingly to the explicit formula of the different solutions we pick below.

### 4.1.1 Uniform refinements: test case 1

We consider the analytic solution

$$u(x, y, t) := (1 - x^2)(1 - y^2) \cos(4t). \qquad (75)$$

As for the spatial discretization, we fix $p_x = 2$ and $h = 0.4$. Since the exact solution is a quadratic polynomial in space and we use a quadratic nodal tensor product basis, up to machine precision, the spatial error is zero.

**Uniform time step refinements.** For $p_n^t$ in $\{2, 3\}$, we pick $\tau$ in $\{2 \times 10^{-1}, 10^{-1}, 5 \times 10^{-2}, 2.5 \times 10^{-2}, 1.25 \times 10^{-2}\}$; for $p_n^t$ in $\{4, 5\}$, we pick $\tau$ in $\{2 \times 10^{-1}, 1.25 \times 10^{-1}, 9.09 \times 10^{-2}, 7.15 \times 10^{-2}, 5.88 \times 10^{-2}\}$. We display the errors in (72), (73), and (74) in Figure 1.

<figure><img src="image_3.png" /><figcaption>Figure 1: Exact solution as in (75), uniform $\tau$-refinement.</figcaption></figure>
The expected optimal convergence rates as discussed in Corollary 2.9 are observed for the errors in (72); the same convergence rate is achieved by the other error measures, which is not covered by the theoretical results from Section 2. In particular, the error in (74) has order $O(\tau^{p_n^t - \frac{1}{2}})$; see also Remark 3 for further comments on this point.

**Uniform polynomial degree in time refinements.** We pick $p_n^t$ in $\{2, 3, 4, 5, 6\}$ and $\tau$ in $\{2 \times 10^{-1}, 10^{-1}\}$. We display the errors in (72) (73), and (74) in Figure 2.

<figure><img src="image_4.png" /><figcaption>Figure 2: Exact solution as in (75), uniform $p_n^t$-refinement.</figcaption></figure>

We observe exponential convergence rate for all the errors. Even though this is not covered by the results in Section 2, we can expect this behaviour from the smoothness of the function in (75) and standard p-FEM techniques [34].

### 4.1.2 Uniform refinements: test case 2

We consider analytic in space solutions

$$u(x, y, t) := (1 - x^2)(1 - y^2)t^\alpha, \quad \alpha > 1.5, \qquad (76)$$

which belong to $H^{\alpha+\frac{1}{2}}(0, T; C^\infty(\Omega))$. We fix $p_x = 2$ and $h = 0.4$.

**Uniform time step refinements.** We pick $p_n^t = 2$ and $\tau$ in $\{2 \times 10^{-1}, 10^{-1}, 5 \times 10^{-2}, 2.5 \times 10^{-2}, 1.25 \times 10^{-2}, 6.13 \times 10^{-3}, 3.06 \times 10^{-3}, 1.53 \times 10^{-3}\}$ for $\alpha = 1.75$ and $\tau$ in $\{2 \times 10^{-1}, 10^{-1}, 5 \times 10^{-2}, 2.5 \times 10^{-2}, 1.25 \times 10^{-2}\}$ for $\alpha = 2.5$. We do not consider higher polynomial degrees in time, since they deliver the same convergence rates. The results are exhibited in Figure 3.

<figure><img src="image_5.png" /><figcaption>Figure 3: Exact solution as in (76), uniform $\tau$-refinement.</figcaption></figure>

We observe optimal convergence rates for the errors in (72) as dictated by Corollary 2.9; similar rates are achieved by the error measures in (73). The error measured in the $W^{1,\infty}(L^2)$-seminorm confirms estimate (41); that seminorm converges with the same rate of the error in (74).
**Uniform polynomial degree in time refinements.** We pick $p_n^t$ in $\{2, ..., 10\}$ and $\tau = 0.2$. The results are displayed in Figure 4.

<figure><img src="image_6.png" /><figcaption>**Figure 4:** Exact solution as in (76), $p_n^t$-refinement.</figcaption></figure>

We observe doubling order convergence rate in $p_n^t$, which is standard in p-FEM while approximating functions with growth of $t^\alpha$ type [34, Section 3.3.5], for the jump of the $H^1(L^2)$-seminorm, whereas other quantities display a super-convergence phenomenon.

### 4.1.3 Simultaneous space-time uniform refinements: test case 3

We consider the analytic solution

$$u(x, y, t) := \sin(\pi nx) \sin(\pi my) \cos(\omega\pi t), \quad (77)$$

with $m = n = 1$, $\omega = \sqrt{2}$, and pick $p_x = p_n^t + 1$ and $h = \tau$, i.e., we are interested in simultaneous space and time refinements. On the other hand, for $p_n^t \in \{2, 3\}$, we choose $\tau$ in $\{2 \times 10^{-1}, 10^{-1}, 5 \times 10^{-2}, 2.5 \times 10^{-2}, 1.25 \times 10^{-2}\}$; for $p_n^t = 4$, we choose $\tau$ in $\{2 \times 10^{-1}, 10^{-1}, 6.67 \times 10^{-2}, 5 \times 10^{-2}, 4 \times 10^{-2}\}$; for $p_n^t = 5$, we choose $\tau$ in $\{5 \times 10^{-1}, 3.33 \times 10^{-1}, 2.5 \times 10^{-1}, 2 \times 10^{-1}, 1.67 \times 10^{-1}\}$. The results are displayed in Figure 5.

We observe optimal convergence rates for the errors (72) and (73).

### 4.1.4 Long time behaviour: test case 3

Here, we investigate the influence of the final time $T$ on the errors and estimators for the test case 3 with exact solution in (77) with $m = n = 1$ and $\omega = \sqrt{2}$. In particular, we take $T = 6, 8, 10$, $p_x = p_n^t = 2, 3$, $h = 2\tau = 0.4$. In Table 1, we report the errors and estimators in the selected norms. We present only a few norms compared to the previous tests, since all the $L^\infty$-in-time type norms have similar behaviour to the $L^\infty(L^2)$-norm, and all $L^2$-in-time type norms have similar behaviour to the $L^2(H^1)$-norm.

<table><caption>**Table 1** Exact solution as in (77) with $m = n = 1$ and $\omega = \sqrt{2}$, long time behaviour.</caption><thead><tr><th colspan="5">error (rate)</th></tr><tr><th colspan="5">p = q = 2</th></tr><tr><th>T</th><th>L<sup>∞</sup>(L<sup>2</sup>)</th><th>L<sup>2</sup>(H<sup>1</sup>)</th><th>jump of H<sup>1</sup>(L<sup>2</sup>)</th><th>η</th></tr></thead><tbody><tr><td>6</td><td>1.12e-1</td><td>5.41e-1</td><td>2.06e0</td><td>3.16e0</td></tr><tr><td>8</td><td>1.46e-1 (0.93)</td><td>8.24e-1 (1.46)</td><td>2.32e0 (0.42)</td><td>4.18e0 (0.97)</td></tr><tr><td>10</td><td>1.84e-1 (1.01)</td><td>1.14e0 (1.47)</td><td>2.58e0 (0.47)</td><td>5.23e0 (1.00)</td></tr></tbody></table>

<table><thead><tr><th colspan="5">error (rate)</th></tr><tr><th colspan="5">p = q = 3</th></tr><tr><th>T</th><th>L<sup>∞</sup>(L<sup>2</sup>)</th><th>L<sup>2</sup>(H<sup>1</sup>)</th><th>jump of H<sup>1</sup>(L<sup>2</sup>)</th><th>η</th></tr></thead><tbody><tr><td>6</td><td>1.32e-3</td><td>7.96e-3</td><td>1.97e-1</td><td>5.79e-2</td></tr><tr><td>8</td><td>1.69e-3 (0.85)</td><td>1.03e-2 (0.89)</td><td>2.29e-1 (0.52)</td><td>7.78e-2 (1.03)</td></tr><tr><td>10</td><td>2.07e-3 (0.89)</td><td>1.30e-2 (1.04)</td><td>2.58e-1 (0.52)</td><td>9.94e-2 (1.09)</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Figure 5: Exact solution as in (77), $\tau$-refinement.</figcaption></figure>

From Table 1, we observe a linear dependence on the final time for the $L^\infty(L^2)$-norm and the estimator $\eta$; the $L^2$-in-time type norms have slightly worse dependence for $p = q = 2$; but the same behaviour for $p = q = 3$; the jump error has half an order dependence.

### 4.1.5 Higher oscillating modes: test case 3

Here we investigate the influence on the performance of the scheme in presence higher oscillating modes in the exact solution. We still consider the test case 3 with exact solution in (77) with parameters $m = n = 10$ and $\omega = 10\sqrt{2}$. We fix $T = 1$, $p_n^t = p_x$ in $\{2, 3, 4, 5\}$, with $h = \tau$ in $\{2 \times 10^{-1}, 10^{-1}, 6.67 \times 10^{-2}, 5 \times 10^{-2}, 4 \times 10^{-2}, 3.33 \times 10^{-2}, 2.86 \times 10^{-2}\}$. In Figure 6, we display the errors in the $L^\infty(L^2)$-norm; other norms mentioned in the above sections have also been tested and omitted here for brevity since they have a similar behaviour to the $L^\infty(L^2)$-norm.

<figure><img src="image_8.png" /><figcaption>Figure 6: Exact solution as in (77) with parameters $m = n = 10$ and $\omega = 10\sqrt{2}$: polynomial degrees VS errors.</figcaption></figure>

From Figure 6, we observe that higher polynomial degrees reduce the pre-asymptotic regime.
## 4.2 Efficiency of the error estimator for the semi-discrete in time scheme

Henceforth, the semi-discrete formulation (4) is considered. We check the efficiency of the error estimator $\eta$ in (70) (with $m = N$ and $\xi = t_N$ in (69)), compared to the error measured in the $L^\infty(0, T; L^2(\Omega))$ norm. In particular, we investigate the behaviour of the effectivity index

$$ \kappa := \frac{\eta}{\|e_h\|_{L^\infty(0,T;L^2(\Omega))}}. \qquad (78) $$

We focus on the test cases with exact solutions in (75) and (76), and consider here uniform time steps and polynomial in time refinements.

In Figure 7, we present the results we obtained under uniform time steps refinements. For the test case with exact solution as in (75), we pick $p_n^t$ in $\{2, 3, 4\}$; for the test case with exact solution as in (76), we pick $p_n^t = 2$, $\alpha = 1.75$. For the test case with exact solution as in (77), we pick $m = n = 10$, $\omega = 10\sqrt{2}$, $p_x = 4$, $h = 1.33 \times 10^{-1}$ and $p_n^t$ in $\{2, 3, 4\}$.

<figure><img src="image_9.png" /><figcaption>Figure 7: Exact solution as in (75), (76), and (77) uniform $\tau$-refinement.</figcaption></figure>

The estimator has the optimal convergence rate as the error measured in the $L^\infty(0, T; L^2(\Omega))$ norm. Notably, the effectivity index in (78) seems stable with respect to $\tau$, i.e., is uniformly bounded by a constant with respect to $\tau$.

Then, in Figure 8, the same tests are investigated with $N = 5$, i.e., for a fixed $\tau = 0.2$ under $p_n^t$-refinements in time; for the test case with exact solution in (76) we only consider $\alpha = 1.75$.
<figure><img src="image_10.png" /><figcaption>Figure 8: Exact solution as in (75), (76) and (77), $p_n^t$-refinement.</figcaption></figure>

Also in this case, the estimator has the same convergence rate as that of the $L^\infty(0, T; L^2(\Omega))$ norm of the error. For the test case with exact solution in (75) and (77), the effectivity index $\kappa$ is uniformly bounded in terms of $p_n^t$; for the test case with exact solution in (76), $\kappa$ increases with rate $1/2$ in terms of $p_n^t$.

## 4.3 Adaptive refinements in time

Since the data oscillation terms in the a posteriori error estimates (71) are not dominant, we omit them for simplicity. We consider here an adaptive algorithm with the usual structure

$$ \text{SOLVE} \Rightarrow \text{ESTIMATE} \Rightarrow \text{MARK} \Rightarrow \text{REFINE}. $$

The **ESTIMATE** step is driven by using the error estimator $\eta$; to this aim, we propose an algorithm for the localization of $\eta$, notably to determine a practical value of $m$ in (69). As for the **MARK** step, we use Dörfler's marking (with $\ell^1$-type summation) with a given threshold $\theta$ in (0, 1]. The **REFINE** step is realized by the bisection of the marked time intervals.

In what follows, we denote the Kronecker delta function by $\delta_{i,j}$.

**Adaptive algorithm for the localization of each ESTIMATE step.**
1. Find the interval index $m$ where $\eta_1$ attains the maximum.

2. Compute

$$ \eta_{2,n} := \begin{cases} \frac{2}{\pi} (\tau_n c_3(p_n^t - 1) \|\Delta_x(U - \Pi_{p^t-1}^0 U)\|_{L^1(I_n; L^2(\Omega))} \\ + \tau_n^3 c_2(p_n^t) c_4(p_n^t - 3) \|\|\Delta_x U'\|\|_{(t_{n-1}, \cdot)\|_{0,\Omega}}) & \text{for all } n = 1, \dots, m \\ 0 & \text{otherwise,} \end{cases} $$

and $\eta_2$ by (69) with the above given $m$.

3. Denote $\eta_n = \delta_{m,n}\eta_1 + \eta_{2,n}$ the local error estimator on the time interval $I_n$ for $n = 1, \dots, N$.

### 4.3.1 Numerical results: the adaptive algorithm

We consider the test case with exact solution as in (76) and $\alpha = 1.75$; $p_n^t$ in $\{2, 3, 4\}$; $p_x = 2$; 5 nodes in each space direction; Dörfler's marking parameter $\theta = 0.5$. We define $DoFs = N \times p_n^t \times \text{card}(V_h)$. In Figure 9, we display the $L^\infty(L^2)$ uniform and adaptive errors and estimators $\eta$ in the left panels; the effectivity indices $\kappa$ are presented in the right panels.

<figure><img src="image_11.png" /><figcaption>Figure 9: Exact solution as in (76).</figcaption></figure>

Some remarks for this test case are in order:
<figure><img src="image_12.png" /><figcaption>Figure 10: Exact solution as in (76), time mesh visualization.</figcaption></figure>

• the adaptive algorithm delivers optimal convergence rate in terms of the number *DoFs* of the method;

• the effectivity index is uniformly bounded for fixed $p_n^t$, and the adaptive algorithm asymptotically returns smaller effectivity indices.

In Figure 10, we illustrate the final time meshes produced with $p_n^t = 2$ and 4 from the adaptive and uniform refinement algorithms.

The adaptive algorithm generates a time mesh with a strong grading towards the initial time. Such temporal meshes are more graded in the case $p_n^t = 4$ compared to the case $p_n^t = 2$.

# 5 Conclusions

For a $C^0$-in-time discretization of the wave equation in second order formulation

• we derived a priori estimates for the fully-discrete method, which are explicit in the spatial mesh size, the time steps, and the polynomial degree distributions in space and time (the errors were measured in $L^\infty$-type norms in time);

• we derived fully explicit, reliable a posteriori estimates for the $L^\infty(L^2)$ error in terms of a novel error estimator, which involves jumps of the time derivatives and the spatial Laplacian at the time nodes.

A campaign of numerical experiments revealed that

• for sufficiently smooth solutions, the order of convergence under uniform refinements is optimal also for $L^2$-type errors in time;

• the proposed error estimator is efficient under time step refinements and may be inefficient under *p*-refinements in time.

More recently [9], we analyzed and assessed the performance of a fully discrete scheme for the wave equation in 2nd order formulation with dynamic mesh change. Crucial open questions remain the proof of a localized lower bound for the error estimator in (70), which, to the best of our knowledge, has been an open problem for at least the last two decades. Another challenge is the proof of convergence and optimality of an adaptive algorithm driven by that error estimator; recent results in [11, 12] show that this can be accomplished for parabolic problems without mesh change in the energy norm, while no similar results are currently available for the wave equation.

**Acknowledgements.** ZD has been partially funded by the French National Research Agency (ANR, STEERS, project number ANR-24-CE56-0127-01). 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).

## 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] A. K. Aziz and P. Monk. Continuous finite elements in space and time for the heat equation. *Math. Comp.*, 52(186):255-274, 1989.

[3] I. Babuška and J. Osborn. Analysis of finite element methods for second order boundary value problems using mesh dependent norms. *Numer. Math.*, 34(1):41-62, 1980.

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

[5] 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.

[6] P. Castillo, B. Cockburn, D. Schötzau, and Ch. Schwab. Optimal a priori error estimates for the hp-version of the local discontinuous Galerkin method for convection-diffusion problems. *Math. Comp.*, 71(238):455-478, 2002.

[7] Th. 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.

[8] Th. Chaumont-Frelet and A. Ern. Damped energy-norm a posteriori error estimates for fully discrete approximations of the wave equation using $C^2$-reconstructions. *ESAIM Math. Model. Numer. Anal.*, 59(4):1937-1972, 2025.

[9] Z. Dong, E. H. Georgoulis, L. Mascotto, and Z. Wang. A posteriori error analysis and adaptivity of a space-time finite element method for the wave equation in second order formulation. http://arxiv.org/abs/2509.08537, 2025.

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

[11] M. Feischl. Inf-sup stability implies quasi-orthogonality. *Math. Comp.*, 91(337):2059-2094, 2022.

[12] M. Feischl and D. Niederkofler. Optimal adaptive implicit time stepping. https://arxiv.org/abs/2506.18809, 2025.

[13] D. A. French. A space-time finite element method for the wave equation. *Comput. Methods Appl. Mech. Engrg.*, 107(1-2):145-157, 1993.

[14] 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.

[15] 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.

[16] 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.

[17] E. H. Georgoulis, O. Lakkis, and Th. P. Wihler. A posteriori error bounds for fully-discrete hp-discontinuous Galerkin timestepping methods for parabolic problems. *Numer. Math.*, 148:363-386, 2021.

[18] S. Gómez and V. Nikolič. Combined DG-CG finite element method for the Westervelt equation. *IMA J. Numer. Anal.*, page draf080, 09 2025.

[19] 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.

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

[21] M. Grote, S. Michel, and S. Sauter. Stabilized leapfrog based local time-stepping method for the wave equation. *Math. Comp.*, 90(332):2603-2643, 2021.

[22] 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.

[23] 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.

[24] 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.

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

[26] 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.

[27] 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.
[28] Ch. G. Makridakis. On the Babuška-Osborn approach to finite element analysis: $L^2$ estimates for unstructured meshes. *Numer. Math.*, 139(4):831–844, 2018.

[29] F. Müller, D. Schötzau, and Ch. Schwab. Symmetric interior penalty discontinuous Galerkin methods for elliptic problems in polygons. *SIAM J. Numer. Anal.*, 55(5):2490–2521, 2017.

[30] N. Nataraj, R. Ruiz-Baier, and A. Yousuf. Semi and fully-discrete analysis of lowest-order nonstandard finite element methods for the biharmonic wave problem. *Comput. Methods Appl. Math.*, 25(4):921–948, 2025.

[31] 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. [Collection of Applied Mathematics for the Master's Degree]. Masson, Paris, 1983.

[32] D. Schötzau and Ch. Schwab. Time discretization of parabolic problems by the hp-version of the discontinuous Galerkin finite element method. *SIAM J. Numer. Anal.*, 38(3):837–875, 2000.

[33] 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.

[34] 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.

[35] E. Süli. A posteriori error analysis and adaptivity for finite element approximations of hyperbolic problems. In *An Introduction to Recent Developments in Theory and Numerics for Conservation Laws: Proceedings of the International School on Theory and Numerics for Conservation Laws*, Freiburg/Littenweiler, October 20–24, 1997, pages 123–194. Springer, 1999.

[36] V. Thomée. Galerkin Finite Element Methods for Parabolic Problems, volume 25. Springer Science & Business Media, 2007.

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

[38] 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.

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