# Damped energy-norm a posteriori error estimates using C2-reconstructions for the fully discrete wave equation with the leapfrog scheme

Théophile Chaumont-Frelet, Alexandre Ern

► To cite this version:

Théophile Chaumont-Frelet, Alexandre Ern. Damped energy-norm a posteriori error estimates using C2-reconstructions for the fully discrete wave equation with the leapfrog scheme. ESAIM: Mathematical Modelling and Numerical Analysis, In press. hal-04511867

HAL Id: hal-04511867

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

Submitted on 19 Mar 2024

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Damped energy-norm a posteriori error estimates for fully discrete approximations of the wave equation using $C^2$-reconstructions*

T. Chaumont-Frelet†, A. Ern‡

March 19, 2024

## Abstract

We derive a posteriori error estimates for the scalar wave equation discretized in space by continuous finite elements and in time by the explicit leapfrog scheme. Our analysis combines the idea of invoking extra time-regularity for the right-hand side, as previously introduced in the space semi-discrete setting, with a novel, piecewise quartic, globally twice-differentiable time-reconstruction of the fully discrete solution. Our main results show that the proposed estimator is reliable and efficient in a damped energy norm. These properties are illustrated in a series of numerical examples.

**Keywords:** a posteriori error estimates, finite element method, leapfrog scheme, time-integration, wave equation

## 1 Introduction

Given an open, bounded, Lipschitz polyhedron $\Omega \subset \mathbb{R}^d$, $d \ge 1$, with boundary $\partial\Omega$, the time interval $J := [0, +\infty)$, and a source term $f: J \times \Omega \to \mathbb{R}$, the scalar wave equation consists in finding $u: J \times \Omega \to \mathbb{R}$ such that

$$ \ddot{u} - \Delta u = f \quad \text{in } J \times \Omega, \qquad (1.1a) $$

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

$$ u|_{t=0} = \dot{u}|_{t=0} = 0 \quad \text{in } \Omega. \qquad (1.1c) $$

The homogeneous Dirichlet condition (1.1b) is considered for the sake of simplicity, and non-homogeneous coefficients in space could be considered in (1.1a). Instead, the zero initial conditions (1.1c) play a role in our analysis. Moreover, the source term $f$ is assumed to be smooth in time and supported away from zero. Notice that the time-smoothness assumption on $f$ does not preclude dealing with minimal space-regularity in domains generating corner or edge singularities.

The model problem (1.1) is of relevance in many engineering applications, so that its numerical discretization has been extensively developed and analyzed. Here, we shall focus on the method of lines in its simplest form, where continuous finite elements are employed for the space discretization, combined with the (explicit) leapfrog scheme as time-marching scheme. This is one of the most frequently used methods to discretize (1.1) owing to its computational efficiency with appropriate mass lumping techniques. Recall, in particular, that stability of explicit time-integrators for the wave equation is typically achieved under a mild CFL condition

<footnote>†Inria Univ. Lille and Laboratoire Paul Painlevé, 59655 Villeneuve-d'Ascq, France</footnote>

<footnote>‡CERMICS, Ecole des Ponts, 77455 Marne-la-Vallée Cedex 2, France and Inria Paris, 75589 Paris, France</footnote>

<footnote>*This work was supported by the ANR JCJC project APOWA (research grant ANR-23-CE40-0019-01).</footnote>
of the form $\tau \lesssim h$, which is often required for accuracy reasons anyway (here, $h$ denotes the mesh size and $\tau$ the time step).

In this work, we are interested in rigorously estimating the discretization error using an a posteriori error estimator. We aim at deriving both upper and lower bounds on the error by the estimator, i.e., reliability and efficiency properties. Perhaps surprisingly, only few works address this question in the literature, as compared to the vast number of references dealing with elliptic and parabolic problems (see, e.g., [9, 10, 18, 19, 20] and the references therein). The main obstruction for deriving reliable and efficient a posteriori error estimates for the wave equation is the lack of an inf-sup stability framework in natural norms [21, Theorem 4.2.23]. Advances towards inf-sup stable variational formulations of (1.1) have been recently reported in [3, 11, 21]. Moreover, in the context of boundary integral equations, we refer the reader to [16] for a least-squares approach, and to [14] for residual a posteriori error estimates. However, all the above formulations lead to fully coupled space-time discretizations, whereas, in the present work, our objective is instead to cover the method of lines.

The a posteriori error analysis of the wave equation discretized with the method of lines has been previously addressed in [2, 12, 13, 15]. On the one hand, an implicit time discretization is considered in [2, 12, 15], either using a second-order backward differentiation formula for the second-order time derivative or a Newmark-type scheme. The error measure is the $H^1(J; H^{-1}(\Omega)) \cap L^2(J; L^2(\Omega))$-norm in [2], the $L^\infty(J; L^2(\Omega))$-norm in [12], and the energy-norm $(H^1(J; L^2(\Omega)) \cap L^2(J; H_0^1(\Omega)))$ in [15]. Among these three works, only [2] also derives error lower bounds, but the efficiency result is somewhat polluted by the presence of additional terms involving the error energy-norm. On the other hand, the (explicit) leapfrog scheme is considered in [13], but in an infinite-dimensional setting for the space variable. This makes it difficult to consider a finite element discretization since the time-discrete solutions in [13] must belong to the domain of iterated powers of the Laplacian operator.

The above discussion clearly indicates that a gap remains in the literature concerning the a posteriori error analysis of the fully discrete wave equation using the (explicit) leapfrog scheme. The goal of the present work is to partly fill this gap. To this purpose, we rely on the approach recently introduced in [5] in the space semi-discrete setting to bypass the lack of inf-sup stability of the wave equation. There are two key ideas in [5]. The first one is to (abstractly) work with the Laplace transform and distinguish low- and high-frequency components of the error. The former can be controlled by invoking a duality argument, and the latter by invoking the time smoothness of the source term $f$. The second idea is to bound the space semi-discrete error, $e$, using the following damped energy-norm:

$$ \mathcal{E}_{\rho}^2(e) := \int_0^{+\infty} \{ \|\dot{e}(t)\|_{\Omega}^2 + \|\nabla e(t)\|_{\Omega}^2 \} e^{-2\rho t} dt, \quad (1.2) $$

where the damping parameter $\rho > 0$ scales as the reciprocal of a time and can be chosen as small as desired. In practice, one is typically interested in the solution up to some finite time-horizon $T_*$, and one then sets $\rho$ to be, e.g., the reciprocal of some multiple of $T_*$.

The main step forward accomplished herein is to extend [5] to the fully discrete setting. This step is by no means straightforward. The key idea is to introduce suitable time-reconstructed functions from the sequence of fully discrete solutions produced by the leapfrog scheme at the discrete time nodes. One crucial difficulty is that reformulating the leapfrog scheme in a time-functional setting requires introducing two time reconstructions, one which is piecewise quadratic in time and globally of class $C^0$, say $u_{h\tau}$, and the other which is piecewise quartic in time and globally of class $C^2$, say $w_{h\tau}$. The idea of introducing a time reconstruction from values produced at the discrete time nodes by a time-marching scheme is already known in the context of a posteriori error analysis [1]. In the context of the wave equation for instance, a piecewise cubic reconstruction, which is globally of class $C^1$ in time, is introduced in [12]. However, the $C^2$ time reconstruction introduced herein seems, to our knowledge, a novel idea in the analysis of the wave equation.

Using the above two reconstructions $u_{h\tau}$ and $w_{h\tau}$, the leapfrog scheme can be rewritten as
follows: For all $t \in J$ and all $v_h \in V_h$,

$$ (\ddot{w}_{h\tau}(t), v_h)_\Omega + (\nabla u_{h\tau}(t), \nabla v_h)_\Omega = (f_\tau(t), v_h)_\Omega. \quad (1.3) $$

Here, $V_h \subset V := H_0^1(\Omega)$ is a finite element space and $f_\tau$ a suitable approximation in time of the source term $f$ (precise notation is introduced in Section 2). Defining the fully discrete error

$$ e := u - w_{h\tau}, \quad (1.4) $$

the main consequence of (1.3) is that the error equation takes the following form: For all $t \in J$ and all $v \in V$,

$$ (\ddot{e}(t), v)_\Omega + (\nabla e(t), \nabla v)_\Omega = (\eta_f(t), v)_\Omega + (\nabla \delta_{h\tau}(t), \nabla v)_\Omega, \quad (1.5) $$

where $\eta_f := f - f_\tau$ represents a data oscillation term (which can be shown to decay at higher order in time than the error itself) and

$$ \delta_{h\tau} := u_{h\tau} - w_{h\tau}. \quad (1.6) $$

The last term on the right-hand side of (1.5) leads to a lack of Galerkin orthogonality in the error equation. The consequence of this fact is that norms of $\delta_{h\tau}$ appear in both upper and lower bounds on the error. Notice that $\delta_{h\tau}$ is fully computable and can be viewed as a time discretization error estimator. As confirmed by our numerical experiments, the contribution of $\delta_{h\tau}$ to the a posteriori error estimator can be made small by decreasing the time step. We also emphasize that the lack of Galerkin orthogonality is also present when dealing with the heat equation discretized using continuous finite elements and an implicit time-scheme (e.g., backward Euler), as already highlighted in [9].

Let us briefly outline our main results. Let $\mathcal{E}_\rho^2(e)$ be the damped energy norm of the fully discrete error $e = u - w_{h\tau}$ (see (1.2)). Our main result in Theorem 4.3 below states that, up to higher-order terms (h.o.t.),

$$ \mathcal{E}_\rho^2(e) \le \Lambda_\rho^2 := \int_0^{+\infty} \left\{ \eta_h^2(t) + \frac{10}{\rho^2} \|\nabla \dot{\delta}_{h\tau}(t)\|_{\Omega}^2 \right\} e^{-2\rho t} dt + \text{h.o.t.} \quad (1.7) $$

with

$$ \eta_h(t) := \|\ddot{w}_{h\tau}(t) - \Delta u_{h\tau}(t) - f_\tau(t)\|_{H^{-1}(\Omega)}. \quad (1.8) $$

Notice that the error estimator, $\Lambda_\rho$, defined in (1.7) is the sum of two terms, which may be linked, respectively, to the space and the time discretization errors. Moreover, since $\ddot{w}_{h\tau} - \Delta u_{h\tau} - f_\tau$ enjoys a Galerkin orthogonality property at all times $t \in J$ owing to (1.3), the a posteriori estimator $\eta_h(t)$ can be bounded by any technique available in the elliptic context, e.g., of residual-type or based on flux equilibration. A converse bound on $\eta_h$ is established in Theorem 5.1 under a CFL constraint on the time step, namely (up to a constant independent of the mesh size $h$, time step $\tau$, and damping parameter $\rho$)

$$ \int_0^{+\infty} \eta_h^2(t) e^{-2\rho t} dt \lesssim \mathcal{E}_\rho^2(e) + \int_0^{+\infty} \left\{ \|\nabla \delta_{h\tau}(t)\|^2 + \frac{1}{\rho^2} \|\nabla \dot{\delta}_{h\tau}(t)\|_{\Omega}^2 \right\} e^{-2\rho t} dt + \text{h.o.t.} \quad (1.9) $$

Thus, if the term

$$ \mathcal{D}_\rho^2 := \int_0^{+\infty} \left\{ \|\nabla \delta_{h\tau}(t)\|^2 + \frac{1}{\rho^2} \|\nabla \dot{\delta}_{h\tau}(t)\|_{\Omega}^2 \right\} e^{-2\rho t} dt \quad (1.10) $$

is added to the error measure and to the estimator by setting $\tilde{\mathcal{E}}_\rho^2(e) := \mathcal{E}_\rho^2(e) + \mathcal{D}_\rho^2$ and $\tilde{\Lambda}_\rho^2 := \Lambda_\rho^2 + \mathcal{D}_\rho^2$, one obtains the following reliability and efficiency result:

$$ \tilde{\Lambda}_\rho^2 - \text{h.o.t.} \lesssim \tilde{\mathcal{E}}_\rho^2(e) \lesssim \tilde{\Lambda}_\rho^2 + \text{h.o.t.} \quad (1.11) $$
Incorporating a fully computable term like $\mathcal{D}_\rho$ in the error and the estimator is standard in other contexts, such as discontinuous Galerkin discretizations of elliptic problems [17]. In our numerical experiments on the wave equation, we observe that $\Lambda_\rho$ controls $\mathcal{E}_\rho$ by a factor of at most 10.

The remainder of the paper is organized as follows. We make the continuous and discrete settings precise in Sections 2. We present the time reconstructions in Section 3, where we also investigate their accuracy. Sections 4 and 5 are, respectively, dedicated to establishing the upper and lower bounds on the error. We present numerical examples in Section 6. In Appendix 7, we derive some stability results on the leapfrog scheme in the damped energy norm. These results are useful in establishing the error lower bounds, but we believe they are of independent interest. Finally, in Appendix 8, we derive, for completeness, an a priori error estimate on the exact solution to the wave equation in the damped energy norm.

# 2 Continuous and discrete settings

In this section, we present the continuous problem and state its basic stability properties in the damped energy norm. Next, we present the discrete setting based on continuous finite elements for the space discretization and the leapfrog scheme for the time discretization. Finally, we recall the classical CFL stability condition for the leapfrog scheme, and we introduce an approximation factor quantifying how well the finite element space approximates some dual solution which is used to establish the error upper bound.

## 2.1 Time and frequency domains

We use standard notation for the Lebesgue, Sobolev, and Bochner-Sobolev spaces. In particular, $(\cdot, \cdot)_\Omega$ denotes the inner product in $L^2(\Omega)$ and $\|\cdot\|_\Omega$ the corresponding norm, and we employ the same notation for vector-valued functions with all components in $L^2(\Omega)$. For any Banach space $Y$ composed of functions defined over $\Omega$, and for any integer $r \ge 0$, we set $C_b^r(J; Y) := \{v \in C^r(J; Y) | f^{(r)} \in L^\infty(J; Y)\}$ and $C_{0,0}^2(J; Y) := \{v \in C^2(J; Y) | v(0) = \dot{v}(0) = 0\}$. Moreover, for every function $\phi \in H^1(J; L^2(\Omega)) \cap L^2(J; V)$ with $V := H_0^1(\Omega)$, we define its damped energy norm as

$$ \mathcal{E}_{\rho}^{2}(\phi) := \int_{0}^{+\infty} \left\{ \|\dot{\phi}(t)\|_{\Omega}^{2} + \|\nabla\phi(t)\|_{\Omega}^{2} \right\} e^{-2\rho t} dt. \quad (2.1) $$

We assume that the source term satisfies $f \in C_b^0(J; L^2(\Omega))$ and that it is supported away from zero. It is then natural to seek for a strong solution $u \in C_{0,0}^2(J; L^2(\Omega)) \cap C^0(J; V)$ such that

$$ (\ddot{u}(t), v)_\Omega + (\nabla u(t), \nabla v)_\Omega = (f(t), v)_\Omega, \quad \forall t \in J, \forall v \in V. \quad (2.2) $$

Notice that the boundary condition (1.1b) is encoded in the fact that $u \in C^0(J; V)$, and the initial conditions (1.1c) are encoded in the fact that $u \in C_{0,0}^2(J; L^2(\Omega))$.

A convenient way to handle the damped energy-norm introduced in (1.2) is to work in the frequency domain. To this purpose, we consider the Laplace transform

$$ \hat{v}(s) := \int_{0}^{+\infty} v(t)e^{-st}dt, \quad s := \rho + i\nu \in \mathbb{C}, \quad \Re(s) := \rho > 0, \Im(s) := \nu, \quad (2.3) $$

for any function $v \in L^\infty(J; Y)$ so that the integral is properly defined. We observe that the damping parameter $\rho > 0$ introduced in the error measure determines the real part of the complex frequency $s$. A key property of the Laplace transform is that

$$ \int_{0}^{+\infty} \|v(t)\|_{Y}^{2} e^{-2\rho t} dt = \int_{\rho-i\infty}^{\rho+i\infty} \|\hat{v}(s)\|_{Y}^{2} ds. \quad (2.4) $$
For all $s \in \mathbb{C}$, we define the sesquilinear form $b_s$ on $\hat{V} \times \hat{V}$ with $\hat{V} := H_0^1(\Omega; \mathbb{C})$ such that

$$b_s(\hat{\phi}, \hat{v}) := s^2(\hat{\phi}, \hat{v})_{\Omega} + (\nabla \hat{\phi}, \nabla \hat{v})_{\Omega}. \qquad (2.5)$$

In the frequency domain, (2.2) can be rewritten as

$$b_s(\hat{u}(s), \hat{v}) = (\hat{f}(s), \hat{v})_{\Omega} \quad \forall s \in \mathbb{C}, \forall \hat{v} \in \hat{V}. \qquad (2.6)$$

## 2.2 A priori estimates on the exact solution

Defining the “frequency-domain energy norm” as

$$\|\hat{v}\|^2 := |s|^2 \|\hat{v}\|_{\Omega}^2 + \|\nabla \hat{v}\|_{\Omega}^2 \quad \forall \hat{v} \in \hat{V}, \qquad (2.7)$$

the key stability property of the sesquilinear form $b_s$ defined in (2.5) is

$$\rho \|\hat{v}\|^2 = \Re(b_s(\hat{v}, s\hat{v})) \quad \forall \hat{v} \in \hat{V}. \qquad (2.8)$$

**Lemma 2.1** (A priori estimate). Assume that $u \in C_b^1(J; L^2(\Omega)) \cap C_b^0(J; V)$. The following holds:

$$\mathcal{E}_{\rho}^{2}(u) \leq \frac{1}{\rho^{2}} \int_{0}^{+\infty} \|f(t)\|_{\Omega}^{2} e^{-2\rho t} dt. \qquad (2.9)$$

*Proof.* Owing to the stability property (2.8), we infer that, for all $s \in \mathbb{C}$,

$$\rho \|\hat{u}(s)\|^2 = \Re(b_s(\hat{u}(s), s\hat{u}(s))) = \Re(\hat{f}(s), s\hat{u}(s))_{\Omega} \leq \|\hat{f}(s)\|_{\Omega} \|\hat{u}(s)\|_{\Omega}.$$

Hence, we have

$$\|\hat{u}(s)\| \leq \frac{1}{\rho} \|\hat{f}(s)\|_{\Omega} \quad \forall s \in \mathbb{C}. \qquad (2.10)$$

The assertion follows from $s\hat{u}(s) = \hat{u}(s)$ and (2.4). $\square$

## 2.3 Discrete problem

The space discretization is realized by means of a conforming finite element method (FEM). We assume that $\Omega$ is a polyhedron and consider an affine simplicial mesh $\mathcal{T}_h$ that covers $\Omega$ exactly. The subscript $h$ refers to the mesh size. Let $V_h \subset V$ be the continuous FEM space built using Lagrange finite elements of degree $k \geq 1$. The time discretization is realized by considering an increasing sequence of discrete time nodes $(t^n)_{n \in \mathbb{N}}$ (conventionally, $0 \in \mathbb{N}$). For simplicity, we assume a constant time-step $\tau$ so that $t^n = n\tau$ for all $n \in \mathbb{N}$. Recalling that the damping parameter $\rho$ is such that $\rho T \geq 1$, where $T$ is some observation time-scale, we make in what follows the following mild assumption:

$$\rho\tau \leq 1. \qquad (2.11)$$

The second-order time derivative in the scalar wave equation is discretized by means of the leapfrog (central finite difference) scheme. The fully discrete wave equation consists in finding the sequence $(U^n)_{n \in \mathbb{N}} \subset V_h^N$, i.e., $U^n \in V_h$ for all $n \in \mathbb{N}$, such that

$$\frac{1}{\tau^2}(U^{n+1} - 2U^n + U^{n-1}, v_h)_{\Omega} + (\nabla U^n, \nabla v_h)_{\Omega} = (f^n, v_h)_{\Omega} \quad \forall n \geq 1, \forall v_h \in V_h, \qquad (2.12)$$

with $f^n := f(t^n)$ and the initial conditions $U^0 = U^1 = 0$. (Recall that $f(0) = 0$ by assumption.)

**Remark 2.2** (Time-horizon). In practice, one fixes a finite time-horizon $T_* = N\tau$ for some positive integer $N$ and computes only the first $N$ steps of the scheme (2.12).
## 2.4 CFL condition

It is well-known that the leapfrog scheme is conditionally stable under a CFL condition. This condition is not needed in our analysis to derive the error upper bound, but it is invoked in the error lower bound. To state the CFL restriction, one introduces on $V_h \times V_h$ the bilinear form

$$m_{h\tau}(v_h, w_h) := (v_h, w_h)_{\Omega} - \tau^2(\nabla v_h, \nabla w_h)_{\Omega}. \qquad (2.13)$$

Then, for all $\mu_0 \in (0, 1)$, there exists $\tau^*(\mu_0) > 0$ such that, whenever $\tau \in (0, \tau^*(\mu_0)]$, the following holds:

$$\mu_0(v_h, w_h)_{\Omega} \le m_{h\tau}(v_h, w_h) \le (v_h, w_h)_{\Omega} \quad \forall v_h, w_h \in V_h. \qquad (2.14)$$

(Notice that the second bound is trivial.) One can take $\mu_0 = \frac{1}{2}$ in what follows to fix the ideas. Invoking an inverse inequality, one readily shows that $\tau^*(\mu_0) \le Ch_{\min}$, where $h_{\min}$ is the smallest diameter of a mesh cell and the constant $C$ depends on $\mu_0$ and the shape-regularity parameter of the mesh $T_h$. We leave the dependence on the latter parameter implicit and write the CFL condition in the form

$$\tau \le C(\mu_0)h_{\min}, \qquad (2.15)$$

for some positive constant $C(\mu_0)$ such that (2.14) holds true for some $\mu_0 \in (0, 1)$.

## 2.5 Approximation factor

As in [5], the derivation of the error upper bound involves a duality argument. To this purpose, for all $\hat{g} \in L^2(\Omega; \mathbb{C})$, we consider the solution $\hat{\chi}_g \in \hat{V}$ to the adjoint problem $b_s(\hat{w}, \hat{\chi}_g) = |s|^2(\hat{w}, \hat{g})_{\Omega}$ for all $\hat{w} \in \hat{V}$. We define the approximation factor

$$\gamma_s(h) := \sup_{\substack{\hat{g} \in L^2(\Omega; \mathbb{C}) \\ |s|\|\hat{g}\|_{\Omega}=1}} \min_{\hat{v}_h \in \hat{V}_h} \|\nabla(\hat{\chi}_g - \hat{v}_h)\|_{\Omega}, \qquad (2.16)$$

where $\hat{V}_h$ is the complex-valued version of the finite element space $V_h$ introduced above.

**Lemma 2.3** (Bound on approximation factor). *The following holds:*

$$\chi_s(h) \le \min \left\{ \frac{|s|}{\rho}, C_{\text{app}} C_{\text{ell}} h^\theta \ell_\Omega^{1-\theta} |s| \left(1 + \frac{|s|}{\rho}\right) \right\}, \qquad (2.17)$$

where $\ell_\Omega := \text{diam}(\Omega)$ is a global length scale introduced for dimensional consistency, $C_{\text{app}}$ is related to the approximation properties of the finite element space $\hat{V}_h$, and $\theta \in (\frac{1}{2}, 1]$ and $C_{\text{ell}}$ are related to the elliptic regularity shift in $\Omega$.

*Proof.* Let $\hat{g} \in L^2(\Omega; \mathbb{C})$ be s.t. $|s|\|\hat{g}\|_{\Omega} = 1$. Invoking (2.8) and taking $\hat{v}_h = 0$ readily shows that

$$\gamma_s(h) \le \|\nabla \hat{\chi}_g\|_{\Omega} \le \|\hat{\chi}_g\| \le \frac{|s|}{\rho}.$$

Moreover, invoking elliptic regularity in Lipschitz polyhedra (see [7, p. 158]), we infer that there are $\theta \in (\frac{1}{2}, 1]$ and $C_{\text{ell}} > 0$ such that $\|\hat{\chi}_g\|_{H^{1+\theta}(\Omega)} \le C_{\text{ell}} \ell_\Omega^2 \|\Delta \hat{\chi}_g\|_{\Omega}$. Using the approximation properties of the finite element space $\hat{V}_h$ to bound the minimum over $\hat{v}_h$ in (2.16), this gives

$$\|\nabla(\hat{\chi}_g - \hat{v}_h)\|_{\Omega} \le C_{\text{app}} h^\theta \ell_\Omega^{1-\theta} \|\hat{\chi}_g\|_{H^{1+\theta}(\Omega)} \le C_{\text{app}} C_{\text{ell}} h^\theta \ell_\Omega^{1-\theta} \|\Delta \hat{\chi}_g\|_{\Omega} \\ \le C_{\text{app}} C_{\text{ell}} h^\theta \ell_\Omega^{1-\theta} |s| \left(1 + \frac{|s|}{\rho}\right),$$

where the last bound follows from $\Delta \hat{\chi}_g = s^2 \hat{\chi}_g - |s|^2 \hat{g}$, $|s|\|\hat{g}\|_{\Omega} = 1$, and $|s|\|\hat{\chi}_g\|_{\Omega} \le \|\hat{\chi}_g\| \le \frac{|s|}{\rho}$. $\square$
For any cutoff frequency $\omega > 0$ and any damping parameter $\rho > 0$, we set

$$ \gamma_{\rho,\omega}(h) := \max_{\substack{s=\rho+i\nu \\ |\nu| \le \omega}} \gamma_s(h). \qquad (2.18) $$

Bounds on $\gamma_{\rho,\omega}(h)$ are readily derived from Lemma 2.3. For instance, $\gamma_{\rho,\omega}(h) \le (1 + (\frac{\omega}{\rho})^2)^{\frac{1}{2}}$. More importantly, fixing $\omega > 0$ and $\rho > 0$, we infer that $\gamma_{\rho,\omega}(h) \to 0$ as $h \to 0$.

# 3 Time reconstructions for the leapfrog scheme

In this section, we introduce two time reconstructions defined on sequences in $V_h$. The two reconstructions satisfy an important commuting property with the second-order (discrete) time derivative and allow us to rewrite the leapfrog scheme in a time-functional setting.

For all $n \in \mathbb{N}$, we consider the time interval $J_n := [t^n, t^{n+1})$. Let $Y$ be a Banach space composed of functions defined over $\Omega$; typical examples include $Y = V_h$, $Y = V$ or $Y = L^2(\Omega)$. Given a polynomial degree $\ell \ge 0$, we define the following broken polynomial space in time:

$$ P^{\ell}(J_{\tau}; Y) := \{v_{h\tau} \in L^{\infty}(J; Y) \mid v_{h\tau}|_{J_n} \in \mathbb{P}^{\ell}(J_n; Y) \ \forall n \in \mathbb{N}\}, \qquad (3.1) $$

where $\mathbb{P}^{\ell}(J_n; Y)$ is composed of $Y$-valued time-polynomials of order at most $\ell$ restricted to $J_n$.

## 3.1 Definitions and key properties

Consider a sequence $V := (V^n)_{n \in \mathbb{N}}$ in $Y$. Henceforth, any such sequence is extended by zero for negative indices, i.e., we conventionally set $V^{-1} = V^{-2} = \dots := 0$. The first time reconstruction we consider is the function

$$ R(V) \in P^2(J_{\tau}; Y) \cap C^0(J; Y), \qquad (3.2) $$

which is defined such that, for all $n \in \mathbb{N}$ and all $t \in J_n$,

$$ R(V)(t) := V^n + \frac{V^{n+1} - V^{n-1}}{2\tau}(t - t^n) + \frac{V^{n+1} - 2V^n + V^{n-1}}{\tau^2}\frac{1}{2}(t - t^n)^2. \qquad (3.3) $$

Notice that $R(V)$ is the restriction to $J_n$ of the Lagrange interpolate at the discrete time nodes $t^{n-1}$, $t^n$, $t^{n+1}$. While the fact that $R(V) \in P^2(J_{\tau}; Y)$ is obvious by construction, we now justify the claim $R(V) \in C^0(J; Y)$.

**Lemma 3.1 (Time-reconstruction).** Let $R(V)$ be defined by (3.3). Then, $R(V) \in C^0(J; Y)$.

*Proof.* For all $n \in \mathbb{N}$, we observe that

$$ R(V)(t_{-}^{n+1}) = V^n + \frac{V^{n+1} - V^{n-1}}{2} + \frac{V^{n+1} - 2V^n + V^{n-1}}{2} = V^{n+1} = R(V)(t^{n+1}). $$

This completes the proof. $\square$

The second time reconstruction provides $C^2$-smoothness in time, and its construction hinges on piecewise quartic time-polynomials. Specifically, we consider the time-reconstructed function

$$ L(V) \in P^4(J_{\tau}; Y) \cap C^2(J; Y), \qquad (3.4) $$

which is defined such that, for all $n \in \mathbb{N}$ and all $t \in J_n$,

$$ L(V)(t) := \alpha^n + \beta^n(t - t^n) + \gamma^n \frac{1}{2}(t - t^n)^2 + \delta^n \frac{1}{6}(t - t^n)^3 + \epsilon^n \frac{1}{24}(t - t^n)^4, \qquad (3.5) $$
with the coefficients

$$ \alpha^n = \frac{3\mathbf{V}^{n+1} + 17\mathbf{V}^n + 5\mathbf{V}^{n-1} - \mathbf{V}^{n-2}}{24}, $$

$$ \gamma^n = \frac{\mathbf{V}^{n+1} - 2\mathbf{V}^n + \mathbf{V}^{n-1}}{\tau^2}, $$

$$ \epsilon^n = \frac{\mathbf{V}^{n+2} - 4\mathbf{V}^{n+1} + 6\mathbf{V}^n - 4\mathbf{V}^{n-1} + \mathbf{V}^{n-2}}{\tau^4}. $$

$$ \beta^n = \frac{5\mathbf{V}^{n+1} + 3\mathbf{V}^n - 9\mathbf{V}^{n-1} + \mathbf{V}^{n-2}}{12\tau}, \quad (3.6a) $$

$$ \delta^n = \frac{\mathbf{V}^{n+2} - 2\mathbf{V}^{n+1} + 2\mathbf{V}^{n-1} - \mathbf{V}^{n-2}}{2\tau^3}, \quad (3.6b) $$

$$ (3.6c) $$

As above, the fact that $L(\mathbf{V}) \in P^4(J_\tau; Y)$ is obvious by construction, and we now justify the claim $L(\mathbf{V}) \in C^2(J; Y)$.

**Lemma 3.2** (Smooth time-reconstruction). Let $L(\mathbf{V})$ be defined by (3.5)-(3.6). Then, $L(\mathbf{V}) \in C^2(J; Y)$.

*Proof.* Let us set $v := L(\mathbf{V})$. For all $n \in \mathbb{N}$, elementary manipulations show that

$$ v(t_{-}^{n+1}) = \alpha^n + \beta^n \tau + \gamma^n \frac{1}{2} \tau^2 + \delta^n \frac{1}{6} \tau^3 + \epsilon^n \frac{1}{24} \tau^4 = \alpha^{n+1} = v(t^{n+1}), $$

$$ \dot{v}(t_{-}^{n+1}) = \beta^n + \gamma^n \tau + \delta^n \frac{1}{2} \tau^2 + \epsilon^n \frac{1}{6} \tau^3 = \beta^{n+1} = \dot{v}(t^{n+1}), $$

$$ \ddot{v}(t_{-}^{n+1}) = \gamma^n + \delta^n \tau + \epsilon^n \frac{1}{2} \tau^2 = \gamma^{n+1} = \ddot{v}(t^{n+1}). $$

This concludes the proof. $\square$

**Remark 3.3** (Consistency). Elementary manipulations of Taylor polynomials show that in the case of a sequence $V$ such that $\mathbf{V}^n = v(t^n)$ for some smooth function $v \in C^6(J; Y)$ supported away from zero, the coefficients defining the time-reconstructed function $L(\mathbf{V})$ in (3.6) are such that $\alpha^n = v(t^n) + O(\tau^2)$, $\beta^n = \dot{v}(t^n) + O(\tau^2)$, $\gamma^n = \ddot{v}(t^n) + O(\tau^2)$, $\delta^n = \dddot{v}(t^n) + O(\tau^2)$, and $\epsilon^n = \text{````v}(t^n) + O(\tau^2)$.

Consider a sequence $V := (\mathbf{V}^n)_{n \in \mathbb{N}}$ in $Y$ such that $V^0 = 0$. Define the sequence $(\delta_{\tau}^2 \mathbf{V}^n)_{n \in \mathbb{N}}$ such that $\delta_{\tau}^2 \mathbf{V}^n := \frac{1}{\tau^2}(\mathbf{V}^{n+1} - 2\mathbf{V}^n + \mathbf{V}^{n-1})$, for all $n \in \mathbb{N}$. (It is legitimate to set $\delta_{\tau}^2 \mathbf{V}^{-1} = 0$ since $V^0 = 0$.)

**Lemma 3.4** (Commuting with second-order time derivative). The following holds:

$$ \frac{d^2}{dt^2} L(\mathbf{V})(t) = R(\delta_{\tau}^2 \mathbf{V})(t) \quad \forall t \in J. \quad (3.7) $$

*Proof.* Since $L(\mathbf{V})$ is of class $C^2$ and $R(\delta_{\tau}^2 \mathbf{V})$ is of class $C^0$ owing to Lemma 3.2 and Lemma 3.1, respectively, it suffices to establish the identity (3.7) for all $t \in J_n$ and all $n \in \mathbb{N}$. Recalling the definitions in (3.6), elementary manipulations show that the coefficients $\gamma^n$, $\delta^n$, and $\epsilon^n$ defining $L(\mathbf{V})$ satisfy

$$ \gamma^n = \delta_{\tau}^2 \mathbf{V}^n, \quad \delta^n = \frac{\delta_{\tau}^2 \mathbf{V}^{n+1} - \delta_{\tau}^2 \mathbf{V}^{n-1}}{2\tau}, \quad \epsilon^n = \frac{\delta_{\tau}^2 \mathbf{V}^{n+1} - 2\delta_{\tau}^2 \mathbf{V}^n + \delta_{\tau}^2 \mathbf{V}^{n-1}}{\tau^2}. $$

As a consequence, we have

$$ \frac{d^2}{dt^2} L(\mathbf{V})(t) = \gamma^n + \delta^n(t - t^n) + \epsilon^n \frac{1}{2}(t - t^n)^2 = R(\delta_{\tau}^2 \mathbf{V})(t) \quad \forall t \in J_n, \forall n \in \mathbb{N}. $$

This concludes the proof. $\square$
## 3.2 Rewriting of leapfrog scheme

The commuting property established in Lemma 3.4 allows us to rewrite the leapfrog scheme in a time-functional setting. Let us set

$$f_{\tau} := R(\mathbf{F}) \in P^2(J_{\tau}; L^2(\Omega)) \cap C^0(J; L^2(\Omega)) \quad \text{with} \quad \mathbf{F}^n := f(t^n) \quad \forall n \in \mathbb{N}. \quad (3.8)$$

Notice that the sequence $(\mathbf{F}^n)_{n \in \mathbb{N}}$ can indeed be extended by zero for negative indices since $f$ is supported away from zero. We also set

$$u_{h\tau} := R(\mathbf{U}) \in P^2(J_{\tau}; V_h) \cap C^0(J; V_h), \quad w_{h\tau} := L(\mathbf{U}) \in P^4(J_{\tau}; V_h) \cap C^2(J; V_h), \quad (3.9)$$

where the sequence $\mathbf{U} := (\mathbf{U}^n)_{n \in \mathbb{N}}$ solves the fully discrete scalar wave equation (2.12). The sequence of accelerations is defined as $\mathbf{A} := (\mathbf{A}^n)_{n \in \mathbb{N}}$ with

$$\mathbf{A}^n := \delta_{\tau}^2 \mathbf{U}^n = \frac{1}{\tau^2} (\mathbf{U}^{n+1} - 2\mathbf{U}^n + \mathbf{U}^{n-1}) \quad \forall n \in \mathbb{N}. \quad (3.10)$$

Notice that $\mathbf{A}^0 = 0$ by definition and that it is legitimate to set $\mathbf{A}^{-1} = \mathbf{A}^{-2} = \dots = 0$.

**Lemma 3.5 (Time-reconstructed wave equation).** *The following holds:*

$$(\ddot{w}_{h\tau}(t), v_h)_{\Omega} + (\nabla u_{h\tau}(t), \nabla v_h)_{\Omega} = (f_{\tau}(t), v_h)_{\Omega} \quad \forall t \in J, \forall v_h \in V_h. \quad (3.11)$$

*Proof.* Since the fully discrete wave equation can be rewritten as

$$(\mathbf{A}^n, v_h)_{\Omega} + (\nabla \mathbf{U}^n, \nabla v_h)_{\Omega} = (\mathbf{F}^n, v_h)_{\Omega} \quad \forall n \in \mathbb{N}, \forall v_h \in V_h,$$

the claim follows by applying the time-reconstruction operator $\mathbf{R}$ to this equation and invoking Lemma 3.4 which gives $\mathbf{R}(\mathbf{A})(t) = \frac{d^2}{dt^2} L(\mathbf{U}) = \ddot{w}_{h\tau}(t)$. $\square$

## 3.3 Stability and approximation properties

In this section, we investigate some stability and approximation properties of the time reconstruction operator $\mathbf{R}$. In what follows, for positive real numbers $A$ and $B$, we abbreviate as $A \lesssim B$ the inequality $A \le CB$ with a generic constant $C$ whose value can change at each occurrence as long as it is independent of the time step, the mesh size, the damping parameter $\rho$, and, whenever relevant, any function involved in the bound.

**Lemma 3.6 (Stability).** Let $\mathbf{V} := (\mathbf{V}^n)_{n \in \mathbb{N}}$ be a bounded sequence in $\mathbf{Y}$ and set $\dot{\mathbf{V}}^{n+\frac{1}{2}} := \frac{1}{\tau}(\mathbf{V}^{n+1} - \mathbf{V}^n)$. The following holds:

$$\int_{0}^{+\infty} \|R(\mathbf{V})\|_{\mathbf{Y}}^{2} e^{-2\rho t} dt \lesssim \sum_{n \in \mathbb{N}} \tau \|\mathbf{V}^{n}\|_{\mathbf{Y}}^{2} e^{-2\rho t^{n}}, \quad (3.12a)$$

$$\int_{0}^{+\infty} \|\frac{d}{dt} R(\mathbf{V})\|_{\mathbf{Y}}^{2} e^{-2\rho t} dt \lesssim \sum_{n \in \mathbb{N}} \tau \|\dot{\mathbf{V}}^{n+\frac{1}{2}}\|_{\mathbf{Y}}^{2} e^{-2\rho t^{n}}. \quad (3.12b)$$

*Proof.* Invoking inverse inequalities in time shows that $\|R(\mathbf{V})\|_{L^{\infty}(J_n; \mathbf{Y})} \lesssim \max_{\delta \in \{-1, 0, 1\}} \|\mathbf{V}^{n+\delta}\|_{\mathbf{Y}}$. This gives that

$$\int_{0}^{+\infty} \|R(\mathbf{V})\|_{\mathbf{Y}}^{2} e^{-2\rho t} dt \lesssim \sum_{n \in \mathbb{N}} \tau \max_{\delta \in \{-1, 0, 1\}} \|\mathbf{V}^{n+\delta}\|_{\mathbf{Y}}^{2} e^{-2\rho t^{n}},$$

whence (3.12a) follows since $\rho\tau \le 1$ (see (2.11)). The proof of (3.12b) is similar and uses that $\|\frac{d}{dt} R(\mathbf{V})\|_{L^{\infty}(J_n; \mathbf{Y})} \lesssim \max(\|\dot{\mathbf{V}}^{n-\frac{1}{2}}\|_{\mathbf{Y}}, \|\dot{\mathbf{V}}^{n+\frac{1}{2}}\|_{\mathbf{Y}})$. $\square$
For a function $v \in C^0(J; Y)$ supported away from zero, we set $V^n := v(t^n)$ for all $n \in \mathbb{N}$ (so that $V^0 = v(0) = 0$ by assumption). We extend $v$ by zero for $t \le 0$ and we set $V^{-1} = V^{-2} = \dots := 0$. We now bound the approximation error $v - R(V)$ assuming enough smoothness of $v$ in time.

**Lemma 3.7 (Approximation).** The following holds for every function $v \in C_b^3(J; Y)$ supported away from zero:

$$ \int_0^{+\infty} \|(v - R(\mathbf{V}))(t)\|_{\mathbf{Y}}^2 e^{-2\rho t} dt \lesssim \tau^6 \int_0^{+\infty} \|\ddot{v}(t)\|_{\mathbf{Y}}^2 e^{-2\rho t} dt, \quad (3.13a) $$

$$ \int_0^{+\infty} \| \frac{d}{dt} (v - R(\mathbf{V}))(t) \|_{\mathbf{Y}}^2 e^{-2\rho t} dt \lesssim \tau^4 \int_0^{+\infty} \| \ddot{v}(t) \|_{\mathbf{Y}}^2 e^{-2\rho t} dt. \quad (3.13b) $$

*Proof.* Let $v \in C_b^3(J; Y)$ be supported away from zero and set $\dot{V}^n := \dot{v}(t^n)$, $\ddot{V}^n := \ddot{v}(t^n)$ for all $n \in \mathbb{N}$.

(1) Proof of (3.13a). We observe that

$$ \int_0^{+\infty} \|(v - R(\mathbf{V}))(t)\|_{\mathbf{Y}}^2 e^{-2\rho t} dt \le \sum_{n \in \mathbb{N}} e^{-2\rho t^n} \int_{J_n} \|(v - R(\mathbf{V}))(t)\|_{\mathbf{Y}}^2 dt. \quad (3.14) $$

Let $n \in \mathbb{N}$ and $t \in J_n$. Using a third-order Taylor expansion of $v$ with exact remainder, we observe that

$$ (v - R(\mathbf{V}))(t) = \Delta_1^n (t - t^n) + \Delta_2^n \frac{1}{2} (t - t^n)^2 + \Delta_3(t), \quad (3.15) $$

with

$$ \Delta_1^n := \dot{V}^n - \frac{V^{n+1} - V^{n-1}}{2\tau}, \quad \Delta_2^n := \ddot{V}^n - \frac{V^{n+1} - 2V^n + V^{n-1}}{\tau^2}, \quad \Delta_3(t) := \frac{1}{2} \int_{t^n}^{t} (t-s)^2 \ddot{v}(s) ds. $$

Invoking first-order Taylor expansions with exact remainder gives

$$ V^{n+1} = V^n + \tau \dot{V}^n + \int_{t^n}^{t^{n+1}} (t^{n+1} - s) \ddot{v}(s) ds, \quad V^{n-1} = V^n - \tau \dot{V}^n + \int_{t^{n-1}}^{t} (s - t^{n-1}) \ddot{v}(s) ds. $$

We infer that

$$ \Delta_1^n = \frac{1}{2} \int_{t^n}^{t^{n+1}} \psi^n(s) \ddot{v}(s) ds - \frac{1}{2} \int_{t^{n-1}}^{t^n} \psi^n(s) \ddot{v}(s) ds = \frac{1}{2} \int_{t^n}^{t^{n+1}} \psi^n(s) (\ddot{v}(s) - \ddot{v}(\hat{s})) ds, $$

with $\hat{s} := 2t^n - s$ and $\psi^n$ denotes the hat basis function in time having support in $[t^{n-1}, t^{n+1}]$ and satisfying $\psi^n(t^n) = 1$. Since $\|\ddot{v}(s) - \ddot{v}(\hat{s})\|_Y \le \int_{\hat{s}}^s \|\ddot{v}(\sigma)\|_Y d\sigma \le \tau^{\frac{1}{2}} (\int_{t^{n-1}}^{t^{n+1}} \|\ddot{v}(\sigma)\|_Y^2 d\sigma)^{\frac{1}{2}}$ for all $s \in [t^n, t^{n+1}]$, we infer that

$$ \|\Delta_1^n\|_Y^2 \lesssim \tau^3 \int_{t^{n-1}}^{t^{n+1}} \|\ddot{v}(\sigma)\|_Y^2 d\sigma. $$

Moreover, using second-order Taylor expansions with exact remainder gives

$$ \Delta_2^n = \frac{1}{2} \int_{t^n}^{t^{n+1}} \psi^n(s)^2 \ddot{v}(s) ds - \frac{1}{2} \int_{t^{n-1}}^{t^n} \psi^n(s)^2 \ddot{v}(\hat{s}) ds, $$

so that

$$ \|\Delta_2^n\|_Y^2 \lesssim \tau \int_{t^{n-1}}^{t^{n+1}} \|\ddot{v}(\sigma)\|_Y^2 d\sigma. $$
(Here, it is not necessary to invoke the fourth-order derivative of $v$ to gain an extra power of $\tau$.) Finally, we have

$$ \|\Delta_3\|_{C^0(J_n; Y)}^2 \lesssim \tau^5 \int_{t^n}^{t^{n+1}} \|\ddot{v}(\sigma)\|_Y^2 d\sigma \quad \forall t \in J_n. $$

Combining (3.15) with the above bounds gives

$$ \int_{J_n} \|(v - R(\mathbf{V}))(t)\|_Y^2 dt \lesssim \tau(\tau^2 \|\Delta_1^n\|_Y^2 + \tau^4 \|\Delta_2^n\|_Y^2 + \|\Delta_3\|_{C^0(J_n; Y)}^2) \lesssim \tau^6 \int_{t^{n-1}}^{t^{n+1}} \|\ddot{v}(\sigma)\|_Y^2 d\sigma. $$

Using this estimate in (3.14) gives

$$ \int_{0}^{+\infty} \|(v - R(\mathbf{V}))(t)\|_{Y}^{2} e^{-2\rho t} dt \lesssim \tau^{6} \sum_{n \in \mathbb{N}} e^{-2\rho t} \int_{t^{n-1}}^{t^{n+1}} \|\ddot{v}(t)\|_{Y}^{2} dt, $$

and the assertion follows from

$$ e^{-2\rho t} \int_{t^{n-1}}^{t^{n+1}} \|\ddot{v}(t)\|_Y^2 dt \le \int_{J_{n-1}} \|\ddot{v}(t)\|_Y^2 e^{-2\rho t} dt + e^2 \int_{J_n} \|\ddot{v}(t)\|_Y^2 e^{-2\rho t} dt, $$

recalling that $\rho\tau \le 1$ owing to (2.11).

(2) The proof of (3.13b) is similar and is only sketched. We observe that for all $t \in J_n$,

$$ \frac{d}{dt}(v - R(\mathbf{V}))(t) = \Delta_1^n + \Delta_2^n(t - t^n) + \Delta_4(t), $$

with $\Delta_4(t) := \int_{t^n}^t (t-s)\ddot{v}(s)ds$. Using the above bounds on $\|\Delta_1^n\|_Y$ and $\|\Delta_2^n\|_Y$ together with $\|\Delta_4\|_{C^0(J_n; Y)}^2 \lesssim \tau^3 \int_{t^n}^{t^{n+1}} \|\ddot{v}(s)\|_Y^2 ds$ readily proves (3.13b). $\square$

## 3.4 A priori estimates on the time-reconstruction error

The rewriting of the leapfrog scheme in a time-functional setting naturally leads to the notion of time-reconstruction error defined as

$$ \delta_{h\tau}(t) := u_{h\tau}(t) - w_{h\tau}(t) \quad \forall t \in J. \qquad (3.16) $$

The goal of this section is to derive some a priori estimates on the time-reconstruction error. It is natural to expect that this quantity is second-order accurate in $\tau$. We now establish a more precise result using the damped energy norm. Recall that $\delta_{h\tau} \in C^0(J; V_h)$ and $w_{h\tau} \in C^2(J; V_h)$, so that the weak time-derivatives $\dot{\delta}_{h\tau}$ and $\ddot{w}_{h\tau}$ can be evaluated by computing locally the time derivative(s) in each time interval $J_n$.

**Lemma 3.8** (Bound on time-reconstruction error). Assume that the sequence $(U^n)_{n \in \mathbb{N}}$ solves the leapfrog scheme (2.12) with the initial conditions $U^0 = U^1 = 0$. Let $u_{h\tau}$ and $w_{h\tau}$ be defined in (3.9). Assume that $f \in C_b^3(J; L^2(\Omega))$. Assume the CFL condition (2.15). The following holds:

$$ \int_{0}^{+\infty} \|\dot{\delta}_{h\tau}(t)\|_{\Omega}^{2} e^{-2\rho t} dt \lesssim \frac{\tau^4}{\rho^2} \int_{0}^{+\infty} \|\ddot{f}(t)\|_{\Omega}^{2} e^{-2\rho t} dt, \qquad (3.17a) $$

$$ \int_{0}^{+\infty} \|\nabla \delta_{h\tau}(t)\|_{\Omega}^{2} e^{-2\rho t} dt \lesssim \frac{\tau^4}{\rho^2} \int_{0}^{+\infty} \left\{ \|\ddot{f}(t)\|_{\Omega}^{2} + \tau^2 \|\dddot{f}(t)\|_{\Omega}^{2} \right\} e^{-2\rho t} dt, \qquad (3.17b) $$

$$ \int_{0}^{+\infty} \|\nabla \dot{\delta}_{h\tau}(t)\|_{\Omega}^{2} e^{-2\rho t} dt \lesssim \frac{\tau^4}{\rho^2} \int_{0}^{+\infty} \|\dddot{f}(t)\|_{\Omega}^{2} e^{-2\rho t} dt. \qquad (3.17c) $$

In addition, we have

$$ \int_{0}^{+\infty} \|\ddot{w}_{h\tau}(t)\|_{\Omega}^{2} e^{-2\rho t} dt \lesssim \frac{1}{\rho^2} \int_{0}^{+\infty} \|\ddot{f}(t)\|_{\Omega}^{2} e^{-2\rho t} dt. \qquad (3.17d) $$
Proof. Recall the definition (3.10) of $A^n$ and set $A^{n+\frac{1}{2}} := \frac{1}{2}(A^{n+1}+A^n)$, $\dot{A}^{n+\frac{1}{2}} := \frac{1}{\tau}(A^{n+1}-A^n)$ for all $n \ge -1$. A direct computation shows that for all $t \in J_n$ and all $n \in \mathbb{N}$,

$$ \tau^{-2}\delta_{h\tau}(t) = \frac{1}{12}(A^{n-\frac{1}{2}} + \tau\dot{A}^{n-\frac{1}{2}}) - \frac{\tau}{12}\dot{A}^{n-\frac{1}{2}}\zeta_1^n(t) + \frac{\tau}{2}(\dot{A}^{n+\frac{1}{2}} + \dot{A}^{n-\frac{1}{2}})\zeta_3^n(t) \\ + \tau(\dot{A}^{n+\frac{1}{2}} - \dot{A}^{n-\frac{1}{2}})\zeta_4^n(t), $$

$$ \tau^{-2}\dot{\delta}_{h\tau}(t) = -\frac{1}{12}\dot{A}^{n-\frac{1}{2}} + \frac{1}{2}(\dot{A}^{n+\frac{1}{2}} + \dot{A}^{n-\frac{1}{2}})\zeta_2^n(t) + (\dot{A}^{n+\frac{1}{2}} - \dot{A}^{n-\frac{1}{2}})\zeta_3^n(t), $$

$$ \ddot{w}_{h\tau}(t) = \frac{1}{2}(\dot{A}^{n+\frac{1}{2}} + \dot{A}^{n-\frac{1}{2}}) + (\dot{A}^{n+\frac{1}{2}} - \dot{A}^{n-\frac{1}{2}})\zeta_1^n(t), $$

where we used the shorthand notation $\zeta_m^n(t) := \frac{1}{m!}\frac{1}{\tau^m}(t-t^n)^m$ for all $m \in \{0, \dots, 4\}$. Notice that all of these polynomials are bounded on $J_n$. Hence, we have

$$ \begin{aligned} \int_0^{+\infty} \|\dot{\delta}_{h\tau}(t)\|_{\Omega}^2 e^{-2\rho t} dt &\lesssim \tau^4 \sum_{n \in \mathbb{N}} (\|\dot{A}^{n-\frac{1}{2}}\|_{\Omega}^2 + \|\dot{A}^{n+\frac{1}{2}}\|_{\Omega}^2) \int_{J_n} e^{-2\rho t} dt \\ &\le \tau^4 \sum_{n \in \mathbb{N}} \tau (\|\dot{A}^{n-\frac{1}{2}}\|_{\Omega}^2 + \|\dot{A}^{n+\frac{1}{2}}\|_{\Omega}^2) e^{-2\rho t^n} \\ &\lesssim \tau^4 \sum_{n \in \mathbb{N}} \tau \|\dot{A}^{n+\frac{1}{2}}\|_{\Omega}^2 e^{-2\rho t^n}, \end{aligned} \quad (3.18) $$

where we used that $\dot{A}^{-\frac{1}{2}} = 0$ and $\rho\tau \le 1$ (see (2.11)). Similarly, we have

$$ \int_0^{+\infty} \|\ddot{w}_{h\tau}(t)\|_{\Omega}^2 e^{-2\rho t} dt \lesssim \sum_{n \in \mathbb{N}} \tau \|\dot{A}^{n+\frac{1}{2}}\|_{\Omega}^2 e^{-2\rho t^n}, \quad (3.19a) $$

$$ \int_0^{+\infty} \|\nabla \delta_{h\tau}(t)\|_{\Omega}^2 e^{-2\rho t} dt \lesssim \tau^4 \sum_{n \in \mathbb{N}} \tau (\|\nabla A^{n+\frac{1}{2}}\|_{\Omega}^2 + \tau^2 \|\nabla \dot{A}^{n+\frac{1}{2}}\|_{\Omega}^2) e^{-2\rho t^n}, \quad (3.19b) $$

$$ \int_0^{+\infty} \|\nabla \dot{\delta}_{h\tau}(t)\|_{\Omega}^2 e^{-2\rho t} dt \lesssim \tau^4 \sum_{n \in \mathbb{N}} \tau \|\nabla \dot{A}^{n+\frac{1}{2}}\|_{\Omega}^2 e^{-2\rho t^n}. \quad (3.19c) $$

We can now invoke Lemma 7.3 (see Appendix 7) to conclude the proof. Indeed, (3.17a) and (3.17d) readily follow from (7.7a), (3.18), and (3.19a). The estimate (7.7a) can also be used to bound the first term on the right-hand side of (3.19b). To bound the second term on the right-hand side of (3.19b) as well as the right-hand side of (3.19c), we invoke (7.7b), recalling that $B^n := \frac{1}{\tau^3}(U^{n+1}-3U^n+3U^{n-1}-U^{n-2})$, $B^{n+\frac{1}{2}} := \frac{1}{2}(B^{n+1}+B^n)$, $\dot{B}^{n+\frac{1}{2}} := \frac{1}{\tau}(B^{n+1}-B^n)$, so that $B^{n+1} = \dot{A}^{n+\frac{1}{2}}$ (notice that $B^n$ is meant to approximate the third-order time-derivative). To conclude, we observe that

$$ \|\nabla \dot{A}^{n+\frac{1}{2}}\|_{\Omega} = \|\nabla B^{n+1}\|_{\Omega} \le \|\nabla B^{n+\frac{1}{2}}\|_{\Omega} + \frac{1}{2}\tau \|\nabla \dot{B}^{n+\frac{1}{2}}\|_{\Omega} \le \|\nabla B^{n+\frac{1}{2}}\|_{\Omega} + \frac{1}{2}(1-\mu_0)^{\frac{1}{2}} \|\dot{B}^{n+\frac{1}{2}}\|_{\Omega}, $$

where we used (2.14) owing to the CFL condition. $\square$

# 4 Asymptotically constant-free error upper bound

Our goal is to bound the error $e = u - w_{h\tau}$ (see (1.4)), where $u$ is the exact solution and $w_{h\tau}$ is the $C^2$-reconstruction of the fully discrete solution defined in (3.9). The evolutionary PDE governing the error can be written as

$$ (\ddot{e}(t), v)_{\Omega} + (\nabla e(t), \nabla v)_{\Omega} = (\eta_f(t), v)_{\Omega} + \langle \mathcal{R}(t), v \rangle \quad \forall t \in J, \forall v \in V, \quad (4.1) $$

where we introduced the data time-oscillation term $\eta_f \in C^0(J; L^2(\Omega))$ such that

$$ \eta_f(t) := f(t) - f_{\tau}(t), \quad (4.2) $$
and the residual $\mathcal{R}(t) \in V'$ associated with (3.11) such that

$$ \langle \mathcal{R}(t), v \rangle := (f_{\tau}(t), v)_{\Omega} - (\ddot{w}_{h\tau}(t), v)_{\Omega} - (\nabla w_{h\tau}(t), \nabla v)_{\Omega} \quad \forall v \in V, \quad (4.3) $$

with the brackets denoting the duality pairing between $V'$ and $V$. The $\|\cdot\|_{V'}$-norm is defined by equipping $V$ with the $H^1$-seminorm, i.e., we set $\|\mathcal{R}(t)\|_{V'} := \sup_{v \in V}, \|\nabla v\|_{\Omega} = \langle \mathcal{R}(t), v \rangle$. The key consistency property we shall use for the residual is the following perturbed Galerkin orthogonality:

$$ \langle \mathcal{R}(t), v_h \rangle = (\nabla \delta_{h\tau}(t), \nabla v_h)_{\Omega} \quad \forall v_h \in V_h, \quad (4.4) $$

recalling from (3.16) the time-reconstruction error defined as $\delta_{h\tau}(t) := u_{h\tau}(t) - w_{h\tau}(t)$. We also introduce the modified residual $\mathcal{R}^{\dagger}(t) \in V'$ such that $\mathcal{R}^{\dagger}(t) := \mathcal{R}(t) - \nabla \delta_{h\tau}(t)$, i.e.,

$$ \langle \mathcal{R}^{\dagger}(t), v \rangle := (f_{\tau}(t), v)_{\Omega} - (\ddot{w}_{h\tau}(t), v)_{\Omega} - (\nabla u_{h\tau}(t), \nabla v)_{\Omega} \quad \forall v \in V, \quad (4.5) $$

which satisfies the exact Galerkin orthogonality

$$ \langle \mathcal{R}^{\dagger}(t), v_h \rangle = 0 \quad \forall v_h \in V_h. \quad (4.6) $$

In the frequency domain, using obvious notation, the error equation (4.1) becomes

$$ b_s(\hat{e}(s), \hat{v}) = (\hat{\eta}_f(s), \hat{v})_{\Omega} + \langle \hat{\mathcal{R}}(s), \hat{v} \rangle \quad \forall s \in \mathbb{C}, \forall \hat{v} \in \hat{V}. \quad (4.7) $$

We now derive two bounds on $\|\hat{e}(s)\|$, respectively called low-frequency and high-frequency bounds because in our final error estimate, the first bound will be used for $|s| \le \omega$ and the second bound for $|s| \ge \omega$, where $\omega > 0$ is a cutoff frequency.

**Lemma 4.1 (Low-frequency bound).** The following holds for all $s \in \mathbb{C}$:

$$ \|\hat{e}(s)\|^2 \le \|\hat{\mathcal{R}}(s)\|_{V'}^2 + 8\gamma_s(h)^2 \|\hat{\mathcal{R}}^{\dagger}(s)\|_{V'}^2 + \frac{8}{\rho^2} \|\nabla \delta_{h\tau}(s)\|_{\Omega}^2 + \frac{9}{\rho^2} \|\hat{\eta}_f(s)\|_{\Omega}^2, \quad (4.8) $$

with the approximation factor $\gamma_s(h)$ defined in (2.16).

*Proof.* We need to bound $|s|\|\hat{e}(s)\|_{\Omega}$ and $\|\nabla \hat{e}(s)\|_{\Omega}$ for all $s \in \mathbb{C}$.

(1) Bound on $|s|\|\hat{e}(s)\|_{\Omega}$. Let $\hat{\chi}_e(s) \in \hat{V}$ solve the adjoint problem $b_s(\hat{w}, \hat{\chi}_e(s)) = |s|^2(\hat{w}, \hat{e}(s))_{\Omega}$ for all $\hat{w} \in \hat{V}$. The stability property (2.8) readily implies that

$$ \rho \|\hat{\chi}_e(s)\|^2 = \Re(b_s(\hat{\chi}_e(s), s\hat{\chi}_e(s))) = \Re(\bar{s}|s|^2(\hat{\chi}_e(s), \hat{e}(s))_{\Omega}) \le |s|^2 \|\hat{\chi}_e(s)\| \|\hat{e}(s)\|_{\Omega}, $$

so that

$$ \|\hat{\chi}_e(s)\| \le \frac{|s|}{\rho} |s| \|\hat{e}(s)\|_{\Omega}. $$

Moreover, testing (4.7) with $\hat{v} := \hat{\chi}_e(s)$ gives

$$ \begin{aligned} |s|^2 \|\hat{e}(s)\|_{\Omega}^2 &= b_s(\hat{e}(s), \hat{\chi}_e(s)) = (\hat{\eta}_f(s), \hat{\chi}_e(s))_{\Omega} + \langle \hat{\mathcal{R}}(s), \hat{\chi}_e(s) \rangle \\ &= (\hat{\eta}_f(s), \hat{\chi}_e(s))_{\Omega} + (\nabla \delta_{h\tau}(s), \nabla \hat{\chi}_e(s))_{\Omega} + \langle \hat{\mathcal{R}}^{\dagger}(s), \hat{\chi}_e(s) - \hat{v}_h(s) \rangle, \end{aligned} $$

where we introduced the Laplace-transformed modified residual, $\hat{\mathcal{R}}^{\dagger}(s)$, and we exploited the exact Galerkin orthogonality property (4.6) to introduce in the rightmost term an arbitrary discrete function $\hat{v}_h(s) \in \hat{V}_h$ for all $s \in \mathbb{C}$. Owing to the Cauchy–Schwarz inequality and the above bound on $\|\hat{\chi}_e(s)\|$, we infer that

$$ \begin{aligned} |(\hat{\eta}_f(s), \hat{\chi}_e(s))_{\Omega} + (\nabla \delta_{h\tau}(s), \nabla \hat{\chi}_e(s))_{\Omega}| &\le (\|\hat{\eta}_f(s)\|_{\Omega}^2 + \|\nabla \delta_{h\tau}(s)\|_{\Omega}^2)^{\frac{1}{2}} \frac{1}{|s|} \|\hat{\chi}_e(s)\| \\ &\le (\|\hat{\eta}_f(s)\|_{\Omega}^2 + \|\nabla \delta_{h\tau}(s)\|_{\Omega}^2)^{\frac{1}{2}} \frac{1}{\rho} |s| \|\hat{e}(s)\|_{\Omega}. \end{aligned} $$
Here, we used that $\hat{\delta}_{h\tau}(s) = s\hat{\delta}_{h\tau}(s)$. Moreover, recalling the approximation factor $\gamma_s(h)$ defined in (2.16) and since $\hat{v}_h(s)$ is arbitrary in $\hat{V}_h$, we have

$$ \inf_{\hat{v}_h(s) \in \hat{V}_h} |\langle \hat{\mathcal{R}}^{\dagger}(s), \hat{\chi}_e(s) - \hat{v}_h(s) \rangle| \le \gamma_s(h) \|\hat{\mathcal{R}}^{\dagger}(s)\|_{V'} |s| \|\hat{e}(s)\|_{\Omega}. $$

Putting the above bounds together gives

$$ |s| \|\hat{e}(s)\|_{\Omega} \le \frac{1}{\rho} (\|\hat{\eta}_f(s)\|_{\Omega}^2 + \|\nabla \hat{\delta}_{h\tau}(s)\|_{\Omega}^2)^{\frac{1}{2}} + \gamma_s(h) \|\hat{\mathcal{R}}^{\dagger}(s)\|_{V'} =: \Upsilon. $$

(2) Bound on $\|\nabla \hat{e}(s)\|_{\Omega}$. Since $\|\nabla \hat{v}\|_{\Omega}^2 = b_s(\hat{v}, \hat{v}) - s^2 \|\hat{v}\|_{\Omega}^2$ for all $\hat{v} \in \hat{V}$, we have

$$ \|\nabla \hat{e}(s)\|_{\Omega}^2 \le |b_s(\hat{e}(s), \hat{e}(s))| + \Upsilon^2, $$

and it remains to bound the first term on the right-hand side. Testing (4.7) with $\hat{v} := \hat{e}(s)$ gives

$$ b_s(\hat{e}(s), \hat{e}(s)) = (\hat{\eta}_f(s), \hat{e}(s))_{\Omega} + \langle \hat{\mathcal{R}}(s), \hat{e}(s) \rangle. $$

Using the Cauchy–Schwarz inequality and $\rho \|\hat{e}(s)\|_{\Omega} \le |s| \|\hat{e}(s)\|_{\Omega} \le \Upsilon$, this implies that

$$ |b_s(\hat{e}(s), \hat{e}(s))| \le \|\hat{\mathcal{R}}(s)\|_{V'} \|\nabla \hat{e}(s)\|_{\Omega} + \frac{1}{\rho} \|\hat{\eta}_f(s)\|_{\Omega} \Upsilon. $$

Putting the above bounds together gives

$$ \|\nabla \hat{e}(s)\|_{\Omega}^2 \le \|\hat{\mathcal{R}}(s)\|_{V'} \|\nabla \hat{e}(s)\|_{\Omega} + \frac{1}{\rho} \|\hat{\eta}_f(s)\|_{\Omega} \Upsilon + \Upsilon^2. $$

Using Young's inequality and re-arranging the terms, we conclude that

$$ \|\nabla \hat{e}(s)\|_{\Omega}^2 \le \|\hat{\mathcal{R}}(s)\|_{V'}^2 + \frac{1}{\rho^2} \|\hat{\eta}_f(s)\|_{\Omega}^2 + 3\Upsilon^2. $$

(3) Bound on $\|\hat{e}(s)\|^2$. Combining the bounds from Steps (1) and (2) yields

$$ \|\hat{e}(s)\|^2 \le \|\hat{\mathcal{R}}(s)\|_{V'}^2 + \frac{1}{\rho^2} \|\hat{\eta}_f(s)\|_{\Omega}^2 + 4\Upsilon^2. $$

The claim follows from $\Upsilon^2 \le \frac{2}{\rho^2} \|\hat{\eta}_f(s)\|_{\Omega}^2 + \frac{2}{\rho^2} \|\nabla \hat{\delta}_{h\tau}(s)\|_{\Omega}^2 + 2\gamma_s(h)^2 \|\hat{\mathcal{R}}^{\dagger}(s)\|_{V'}^2$. $\square$

**Lemma 4.2 (High-frequency bound).** The following holds for all $s \in \mathbb{C}$ and for all $r \ge 0$ such that $f \in C_b^T(J; L^2(\Omega))$:

$$ \|\hat{e}(s)\| \le \frac{1}{\rho} (\|\nabla \hat{\delta}_{h\tau}(s)\|_{\Omega}^2 + \|\hat{\eta}_f(s)\|_{\Omega}^2)^{\frac{1}{2}} + \frac{2}{\rho} \frac{1}{|s|^r} \|\widehat{f^{(r)}}(s)\|_{\Omega}. \quad (4.9) $$

*Proof.* The triangle inequality gives

$$ \|\hat{e}(s)\| \le \|\hat{u}(s)\| + \|\hat{w}_{h\tau}(s)\|, $$

and we bound the two terms on the right-hand side.

(1) Owing to (2.10), we infer that

$$ \|\hat{u}(s)\| \le \frac{1}{\rho} \|\hat{f}(s)\|_{\Omega} = \frac{1}{\rho} \frac{1}{|s|^r} \|\widehat{f^{(r)}}(s)\|_{\Omega}, $$

where the last equality follows by invoking the smoothness of the source term $f$ in time.
(2) Taking the Laplace transform of (3.11) and re-organizing the terms gives

$$b_s(\hat{w}_{h\tau}(s), \hat{v}) = (\hat{f}(s), \hat{v})_{\Omega} - (\hat{\eta}_f(s), \hat{v})_{\Omega} - (\nabla \hat{\delta}_{h\tau}(s), \nabla \hat{v})_{\Omega} \quad \forall s \in \mathbb{C}, \forall \hat{v} \in \hat{V}.$$

Owing to the stability property (2.8), we infer that

$$\begin{aligned} \rho \|\hat{w}_{h\tau}(s)\|^2 &= \Re((\hat{f}(s), s\hat{w}_{h\tau}(s))_{\Omega} - (\hat{\eta}_f(s), s\hat{w}_{h\tau}(s))_{\Omega} - (\nabla \hat{\delta}_{h\tau}(s), s\nabla \hat{w}_{h\tau}(s))_{\Omega}) \\ &\le \left\{ (\|\hat{f}(s)\|_{\Omega} + \|\hat{\eta}_f(s)\|_{\Omega})^2 + \|\nabla \hat{\delta}_{h\tau}(s)\|_{\Omega}^2 \right\}^{\frac{1}{2}} \|\hat{w}_{h\tau}(s)\|. \end{aligned} \quad (4.10)$$

Bounding $\|\hat{f}(s)\|_{\Omega}$ as in Step (1) and using that $((a+b)^2+c^2)^{\frac{1}{2}} \le a+(b^2+c^2)^{\frac{1}{2}}$ for nonnegative real numbers $a, b, c$, we infer that

$$\|\hat{w}_{h\tau}(s)\| \le \frac{1}{\rho} \frac{1}{|s|^r} \|\hat{f}^{(r)}(s)\|_{\Omega} + \frac{1}{\rho} (\|\hat{\eta}_f(s)\|_{\Omega}^2 + \|\nabla \hat{\delta}_{h\tau}(s)\|_{\Omega}^2)^{\frac{1}{2}}.$$

Putting the above two bounds together proves the claim. □

**Theorem 4.3 (Error upper bound).** Let the error $e \in C^2(J; V)$ be defined in (1.4) and recall the definition (2.1) of the damped energy norm. Let the functions $\eta_f(t) \in L^2(\Omega)$ and $\delta_{h\tau}(t) \in V_h$ be defined in (4.2) and (3.16), respectively, and let the residuals $\mathcal{R}(t)$, $\mathcal{R}^{\dagger}(t) \in V'$ be defined in (4.3) and (4.5), respectively. Let the cutoff frequency $\omega > 0$ and the parameter $\rho > 0$ be fixed, and let the approximation factor $\gamma_{\rho,\omega}(h)$ be defined in (2.18). Let $r \ge 0$ be such that $f \in C_b^r(J; L^2(\Omega))$. The following holds:

$$\begin{aligned} \mathcal{E}_{\rho}^{2}(e) \le \int_{0}^{+\infty} & \left\{ \|\mathcal{R}(t)\|_{V'}^{2} + 8\gamma_{\rho,\omega}(h)^{2} \|\mathcal{R}^{\dagger}(t)\|_{V'}^{2} \right. \\ & \left. + \frac{10}{\rho^{2}} \|\nabla \hat{\delta}_{h\tau}(t)\|_{\Omega}^{2} + \frac{11}{\rho^{2}} \|\eta_{f}(t)\|_{\Omega}^{2} + \frac{8}{\rho^{2}} \frac{1}{\omega^{2r}} \|f^{(r)}(t)\|_{\Omega}^{2} \right\} e^{-2\rho t} dt. \end{aligned} \quad (4.11)$$

*Proof.* Owing to the identity (2.4), we infer that

$$\mathcal{E}_{\rho}^{2}(e) = \int_{0}^{+\infty} \{\|\dot{e}(t)\|_{\Omega}^{2} + \|\nabla e(t)\|_{\Omega}^{2}\} e^{-2\rho t} dt = \int_{\rho-i\infty}^{\rho+i\infty} \|\hat{e}(s)\|_{\Omega}^{2} ds.$$

We split the integral on the right-hand side depending on whether $|s| \le \omega$ or $|s| \ge \omega$. Owing to Lemma 4.1 and the definition of $\gamma_{\rho,\omega}(h)$, we have

$$\int_{\rho-i\infty}^{\rho+i\infty} \|\hat{e}(s)\|_{\Omega}^{2} 1_{\{|s| \le \omega\}} ds \le \int_{\rho-i\infty}^{\rho+i\infty} \left\{ \|\hat{\mathcal{R}}(s)\|_{V'}^{2} + 8\gamma_{\rho,\omega}(h)^{2} \|\hat{\mathcal{R}}^{\dagger}(s)\|_{V'}^{2} \right. \\ \left. + \frac{8}{\rho^{2}} \|\nabla \hat{\delta}_{h\tau}(s)\|_{\Omega}^{2} + \frac{9}{\rho^{2}} \|\hat{\eta}_{f}(s)\|_{\Omega}^{2} \right\} ds,$$

where $1_A$ denotes the characteristic function of the subset $A \subset \mathbb{C}$. Moreover, invoking Lemma 4.2, we infer that

$$\int_{\rho-i\infty}^{\rho+i\infty} \|\hat{e}(s)\|_{\Omega}^{2} 1_{\{|s| \ge \omega\}} ds \le \int_{\rho-i\infty}^{\rho+i\infty} \left\{ \frac{2}{\rho^2} \|\nabla \hat{\delta}_{h\tau}(s)\|_{\Omega}^2 + \frac{2}{\rho^2} \|\hat{\eta}_f(s)\|_{\Omega}^2 + \frac{8}{\rho^2} \frac{1}{\omega^{2r}} \|\hat{f}^{(r)}(s)\|_{\Omega}^2 \right\} ds.$$

Putting the above two bounds together proves the assertion. □

**Remark 4.4 (Theorem 4.3).** The estimate (4.11) bounds the damped energy norm of the error (recall that the damping parameter $\rho$ is typically proportional to the reciprocal of the simulation time $T_*$) in terms of the dual norm of the residual, $\mathcal{R}$ (representative of the space discretization
error), the time-reconstruction error, $\delta_{h\tau}$ (representative of the time discretization error), the data time-oscillation term, $\eta_f$, and a term depending on higher-order time-derivatives of $f$ which can be made as small as desired by choosing the cutoff frequency $\omega$. Once the cutoff frequency and the damping parameter are fixed, the approximation factor $\gamma_{\rho,\omega}(h)$ can be made as small as desired by taking $h$ small enough (see (2.17)). Finally, we observe that the data time-oscillation term converges to third-order in $\tau$ owing to Lemma 3.7 provided $f \in C_b^3(J; L^2(\Omega))$.

# 5 Error lower bound

Let $V_h$ be the collection of the mesh vertices, which we split into the collection of interior vertices, $V_h^i$, and the collection of boundary vertices, $V_h^b$. For all $a \in V_h$, let $\omega_a$ denote the (interior of the) support of the hat basis function $\psi_a$ attached to $a$. Let us define the functional space

$$ V_a := \begin{cases} \{v \in H^1(\omega_a), (v, 1)_{\omega_a} = 0\} & \text{if } a \in V_h^i, \\ \{v \in H^1(\omega_a), v|_{\Gamma_a} = 0\} & \text{if } a \in V_h^b. \end{cases} \qquad (5.1) $$

where $\Gamma_a := \{x \in \partial\omega_a | \psi_a(x) \neq 0\}$. We localize the residual $\mathcal{R}(t) \in V'$ by introducing the linear form $\mathcal{R}_a(t) \in V_a'$ such that $\langle \mathcal{R}_a(t), v \rangle := \langle \mathcal{R}(t), \psi_a v \rangle$ for all $v \in V_a$, where on the right-hand side $\psi_a v$ is understood to be extended by zero outside $\omega_a$. The dual norm of the local residual is defined such that

$$ \|\mathcal{R}_a(t)\|_{V_a'} := \sup_{v \in V_a, \|\nabla v\|_{\omega_a} = 1} |\langle \mathcal{R}_a(t), v \rangle|. \qquad (5.2) $$

A well-known fact exploiting a local Poincaré inequality in $V_a$ is that (see, e.g., [10, §3.2])

$$ h_a^{-1} \|\psi_a v\|_{\omega_a} + \|\nabla(\psi_a v)\|_{\omega_a} \le C \|\nabla v\|_{\omega_a} \quad \forall v \in V_a, \qquad (5.3) $$

where $h_a$ denotes the diameter of $\omega_a$ and the constant $C$ only depends on the mesh shape-regularity. Observe that $h_a \lesssim h$.

**Theorem 5.1 (Error lower bound).** Assume $f \in C_b^{r+1}(J; L^2(\Omega))$. Assume the CFL condition (2.15). The following holds:

$$ \begin{aligned} \int_0^{+\infty} \|\mathcal{R}(t)\|_{V'}^2 e^{-2\rho t} dt & \lesssim \int_0^{+\infty} \sum_{a \in V_h} \{h_a^2 \|\eta_f(t)\|_{\omega_a}^2 + \|\nabla \delta_{h\tau}(t)\|_{\omega_a}^2\} e^{-2\rho t} dt \\ & + \int_0^{+\infty} \sum_{a \in V_h} \|\nabla e(t)\|_{\omega_a}^2 e^{-2\rho t} dt + \omega^2 h^2 \int_0^{+\infty} \|\dot{e}(t)\|_{\Omega}^2 e^{-2\rho t} dt \\ & + \int_0^{+\infty} \left\{ \frac{h^2}{\rho^2 \omega^{2r}} \|f^{(r+1)}(t)\|_{\Omega}^2 + \frac{h^2}{\rho^2 \omega^2} \|\ddot{f}(t)\|_{\Omega}^2 \right\} e^{-2\rho t} dt. \end{aligned} \qquad (5.4) $$

*Proof.* (1) The triangle inequality gives $\|\mathcal{R}(t)\|_{V'} \le \|\mathcal{R}^\dagger(t)\|_{V'} + \|\nabla \delta_{h\tau}(t)\|_{\Omega}$, where $\mathcal{R}^\dagger(t) \in V'$ is the modified residual defined in (4.5). Owing to the partition-of-unit property of the basis functions, we have $\|\nabla \delta_{h\tau}(t)\|_{\Omega}^2 \le \sum_{a \in V_h} \|\nabla \delta_{h\tau}(t)\|_{\omega_a}^2$, and we are left with bounding $\|\mathcal{R}^\dagger(t)\|_{V'}$.

(2) A classical argument (see [4, 6] and [8, Prop. 34.7]) shows that it is possible to bound $\|\mathcal{R}^\dagger(t)\|_{V'}$ by the Hilbertian sum of the local residual norms $\|\mathcal{R}_a^\dagger(t)\|_{V_a'}$ over the mesh vertices. Let us briefly recall the argument for completeness. Let $v \in V$ with $\|\nabla v\|_{\Omega} = 1$. For all $a \in V_h$, set $\bar{v}_a := \frac{1}{|\omega_a|}(v, 1)_{\omega_a}$ if $a \in V_h^i$, and $\bar{v}_a := 0$ if $a \in V_h^b$. Then, $v|_{\omega_a} - \bar{v}_a \in V_a$ for all $a \in V_h$. Exploiting the partition of unity property $\sum_{a \in V_h} \psi_a = 1$ and the Galerkin
orthogonality property (4.6) which gives $\langle \mathcal{R}^\dagger(t), \psi_a \bar{v}_a \rangle = 0$ for all $a \in \mathcal{V}_h$, we infer that

$$ \begin{aligned} \langle \mathcal{R}^\dagger(t), v \rangle &= \sum_{a \in \mathcal{V}_h} \langle \mathcal{R}^\dagger(t), \psi_a v \rangle = \sum_{a \in \mathcal{V}_h} \langle \mathcal{R}^\dagger(t), \psi_a (v - \bar{v}_a) \rangle \\ &= \sum_{a \in \mathcal{V}_h} \langle \mathcal{R}_a^\dagger(t), v - \bar{v}_a \rangle \le \sum_{a \in \mathcal{V}_h} \| \mathcal{R}_a^\dagger(t) \|_{V'_a} \| \nabla v \|_{\omega_a} \le (d+1)^{\frac{1}{2}} \left\{ \sum_{a \in \mathcal{V}_h} \| \mathcal{R}_a^\dagger(t) \|_{V'_a}^2 \right\}^{\frac{1}{2}}, \end{aligned} $$

where the last bound follows from the Cauchy–Schwarz inequality and the fact that every mesh cell is covered by $(d+1)$ sets $\omega_a$. To bound the local residual norm $\| \mathcal{R}_a^\dagger(t) \|_{V'_a}$, we observe that for all $v \in V_a$,

$$ \langle \mathcal{R}_a^\dagger(t), v \rangle = (\ddot{e}(t), \psi_a v)_{\Omega} + (\nabla e(t), \nabla(\psi_a v))_{\Omega} - (\eta_f(t), \psi_a v)_{\Omega} - (\nabla \delta_{h\tau}(t), \nabla(\psi_a v))_{\Omega}. $$

Using (5.3) readily yields

$$ \| \mathcal{R}_a^\dagger(t) \|_{V'_a} \lesssim h_a \| \ddot{e}(t) \|_{\omega_a} + \| \nabla e(t) \|_{\omega_a} + h_a \| \eta_f(t) \|_{\omega_a} + \| \nabla \delta_{h\tau}(t) \|_{\omega_a}. $$

This implies that

$$ \sum_{a \in \mathcal{V}_h} \| \mathcal{R}_a^\dagger(t) \|_{V'_a}^2 \lesssim h^2 \| \ddot{e}(t) \|_{\Omega}^2 + \sum_{a \in \mathcal{V}_h} \{ \| \nabla e(t) \|_{\omega_a}^2 + h_a^2 \| \eta_f(t) \|_{\omega_a}^2 + \| \nabla \delta_{h\tau}(t) \|_{\omega_a}^2 \}. $$

It only remains to estimate the first term on the right-hand side (notice that we do not localize this term in space).

(3) To bound $h\|\ddot{e}(t)\|_{\Omega}$, we consider as above the Laplace transform, i.e., we bound $h\|\hat{e}(s)\|_{\Omega}$ for all $s \in \mathbb{C}$. In the low-frequency regime ($|s| \le \omega$), we have

$$ h \| \hat{e}(s) \|_{\Omega} = h |s| \| \hat{e}(s) \|_{\Omega} \le h \omega \| \hat{e}(s) \|_{\Omega}. $$

In the high-frequency regime ($|s| \ge \omega$), we first invoke the triangle inequality so that

$$ h \| \hat{e}(s) \|_{\Omega} = h |s|^2 \| \hat{e}(s) \|_{\Omega} \le h |s|^2 \| \hat{u}(s) \|_{\Omega} + h |s|^2 \| \hat{w}_{h\tau}(s) \|_{\Omega}. $$

The first term on the right-hand side is readily bounded since

$$ h|s|^2 \|\hat{u}(s)\|_{\Omega} \le h|s|\|\|\hat{u}(s)\|\| \le \frac{h}{\rho}|s|\|\hat{f}(s)\|_{\Omega} \le \frac{h}{\rho\omega^r}\|\widehat{f^{(r+1)}}(s)\|_{\Omega}, $$

where we used that $|s| \ge \omega$ in the last bound. On the other hand, we have

$$ h|s|^2 \|\hat{w}_{h\tau}(s)\|_{\Omega} \le h \frac{1}{\omega} \|\hat{w}_{h\tau}(s)\|_{\Omega}, $$

and we invoke the bound (3.17a) from Lemma 3.8. This gives

$$ \int_{0}^{+\infty} h^2 \|\ddot{e}(t)\|_{\Omega} e^{-2\rho t} dt \lesssim \int_{0}^{+\infty} \left\{ h^2 \omega^2 \|\dot{e}(t)\|_{\Omega}^2 + \frac{h^2}{\rho^2} \frac{1}{\omega^{2r}} \|f^{(r+1)}(t)\|_{\Omega}^2 + \frac{h^2}{\rho^2 \omega^2} \|\ddot{f}(t)\|_{\Omega}^2 \right\} e^{-2\rho t} dt. $$

(4) Putting everything together proves the claim. $\square$

**Remark 5.2 (Theorem 5.1).** The error lower bound from Theorem 5.1 is not completely satisfactory. First, the term $\int_0^{+\infty} \sum_{a \in \mathcal{V}_h} \|\nabla \delta_{h\tau}(t)\|_{\omega_a}^2 e^{-2\rho t} dt$ is not present in the error upper bound and thus represents an additional contribution of the time-reconstruction error. Fortunately, Lemma 3.8 (see (3.17b)) shows that this term converges to second-order with the time step $\tau$. Second, the last term on the right-hand side of (5.4) is limited to the second-order time derivative of $f$ (one cannot go beyond since the reconstructed function $w_{h\tau}$ is only of class $C^2$ in time). Fortunately, this term can be made as small as desired by choosing a large enough cutoff frequency $\omega$. Finally, the error lower bound requires the CFL condition (2.15) to make use of the stability of the leapfrog scheme.
# 6 Numerical results

In this section, we present numerical results to assess a posteriori error estimates derived in the previous sections.

## 6.1 Setting

Here, we describe the discretization parameters, the error measures, and how we estimate the space discretization error. In all cases, we consider the one-dimensional domain $\Omega := (-L, L)$ with $L := 10$, and set the computational time to $T_* := 1000$. For the values of $\rho$ employed hereafter, this final time is such that $e^{-\rho T_*} \le 5 \cdot 10^{-6}$ in all the simulations, so that stopping the simulation at $t = T_*$ has little effect.

### 6.1.1 Discretization parameters

We use uniform meshes, and continuous finite elements of degree $k \in \{1, 2, 3\}$. We employ mesh sizes ranging from $h = 10$ to $h = 0.039$, i.e., the interval $\Omega$ is divided into at least 2 and at most 512 elements. For each $k$, we denote by $\alpha_k$ the largest value for which the leapfrog scheme is stable with the time step $\tau = \alpha_k h$. We empirically found these values to be $\alpha_1 = 0.59$, $\alpha_2 = 0.26$, and $\alpha_3 = 0.15$. In our examples, given $h$ and $k$, we fix the time step by setting $\tau := r\alpha_k h$ with $r \in (0, 1)$. Unless explicitly specified, we use $r := 0.9$.

When $k = 3$, we also employ finer time-step values so that the order of the time discretization error matches that of the space discretization error. To do so, we fix $\tau_0 = r_0\alpha_3 h_0$ on the coarsest mesh (i.e., $h_0 = 10$), and then adjust $\tau$ so that $\tau^2/h^3$ remains constant for all the finer meshes. This is indicated by the notation "$\tau^2 \sim h^3$" in the graphs below.

### 6.1.2 Error measures

For convenience, we use the following shorthand notation for the errors, where $z$ stands either for $u$ or for $w$:

$$e^2(z) := \int_0^{+\infty} \|u(t) - z_{h\tau}(t)\|_{\Omega}^2 e^{-2\rho t} dt, \quad e^2(U) := \tau \sum_{n=0}^{+\infty} \|u(t^n) - U^n\|_{\Omega}^2 e^{-2\rho t^n}, \quad (6.1a)$$

$$e_x^2(z) := \int_0^{+\infty} \|\nabla(u(t) - z_{h\tau}(t))\|_{\Omega}^2 e^{-2\rho t} dt, \quad e_x^2(U) := \tau \sum_{n=0}^{+\infty} \|\nabla(u(t^n) - U^n)\|_{\Omega}^2 e^{-2\rho t^n}, \quad (6.1b)$$

$$e_t^2(z) := \int_0^{+\infty} \|\dot{u}(t) - \dot{z}_{h\tau}(t)\|_{\Omega}^2 e^{-2\rho t} dt, \quad e_t^2(U) := \tau \sum_{n=0}^{+\infty} \|\dot{u}(t^n) - \frac{U^{n+1} - U^{n-1}}{\tau}\|_{\Omega}^2 e^{-2\rho t^n}. \quad (6.1c)$$

For the error measured in the damped energy-norm, we use the shorthand notation

$$\mathcal{E}_{\rho}^2 := \mathcal{E}_{\rho}^2(u - w_{h\tau}) = e_t^2(w) + e_x^2(w). \quad (6.2)$$

In practice, all the integrals and sums are computed up to $t = T_*$. As argued above, $T_*$ has been chosen in such a way that the tail of the integrals and series is negligible.

### 6.1.3 Estimators

To estimate the space discretization error, we construct a continuous piecewise polynomial function $\sigma_{h\tau}$ with vanishing spatial mean value such that

$$\partial_x \sigma_{h\tau}(t) = f_{\tau}(t) - \ddot{w}_{h\tau}(t) \quad \forall t \in J. \quad (6.3)$$
A simple integration by parts reveals that

$$ \|\mathcal{R}(t)\|_{V'} \le \|\partial_x w_{h\tau}(t) + \sigma_{h\tau}(t)\|_{\Omega}. \quad (6.4) $$

As a result, introducing the quantities

$$ R^2 := \int_{0}^{+\infty} \|\partial_x w_{h\tau}(t) + \sigma_{h\tau}(t)\|_{\Omega}^2 e^{-2\rho t} dt, \quad (6.5a) $$

$$ M^2 := \int_{0}^{+\infty} \frac{1}{\rho^2} \|\partial_x (\dot{u}_{h\tau}(t) - \dot{w}_{h\tau}(t))\|_{\Omega}^2 e^{-2\rho t} dt, \quad (6.5b) $$

and using the following shorthand notation for the estimator:

$$ \Lambda_{\rho}^{2} := R^{2} + 10M^{2}, \quad (6.6) $$

the error upper bound (4.11) rewrites

$$ \mathcal{E}_{\rho} \le \Lambda_{\rho} + \text{h. o. t.} \quad (6.7) $$

## 6.2 Benchmark solutions

We consider two different analytical solutions corresponding to “standing” and “propagating” waves.

### 6.2.1 Standing wave

We set

$$ f(t, x) := -2(t - t_0)e^{-(t-t_0)^2} \sin(a(x - L)), \quad (6.8) $$

where $t_0 := 4$ and $a := m\pi/(2L)$ with $m := 5$. Although $f$ does not vanish at $t = 0$, we do have $|f(0, \cdot)| \le 8 \cdot e^{-16} \le 0.9 \cdot 10^{-6}$, $|\dot{f}(0, \cdot)| \le 62 \cdot e^{-16} \le 7.0 \cdot 10^{-6}$ and $|\ddot{f}(0, \cdot)| \le 464 \cdot e^{-16} \le 5.3 \cdot 10^{-5}$, whereas $\max_{t,x} |f(t, x)| \ge 0.5$. The corresponding solution reads

$$ u(t, x) = \Re(\psi(t - t_0)) \sin(a(x - L)), \quad \psi(t) := e^{iat} \int_{-\infty}^{t} e^{-\theta^2 - ia\theta} d\theta. \quad (6.9) $$

To see this, observe that $\ddot{\psi}(t) = (-2t + ia)e^{-t^2} - a^2\psi(t)$ so that $\Re(\ddot{\psi}(t) + a^2\psi(t)) = -2te^{-t^2}$. We observe that $|u(0, \cdot)| \le 1.1 \cdot 10^{-8}$ and $|\dot{u}(0, \cdot)| \le 8.3 \cdot 10^{-8}$, whereas $\max_{t,x} |u(t, x)| \ge 4 \cdot 10^{-2}$, so that setting the initial conditions to zero produces an error that is (much) smaller than the discretization error. The time profiles of $f$ and $u$ are shown in Figure 1. In practice, we compute $\psi$ with the erf function of the Faddeeva software package by observing that

$$ \psi(t) = e^{-\frac{a^2}{4} + iat} \int_{-\infty}^{t} e^{-(\theta + i\frac{a}{2})^2} d\theta = e^{-\frac{a^2}{4} + iat} \int_{-\infty}^{t + i\frac{a}{2}} e^{-z^2} dz. \quad (6.10) $$

### 6.2.2 Propagating wave

We set

$$ f(t, x) := -2(t - t_0)e^{-(t-t_0)^2} e^{-x^2}, \quad (6.11) $$

where, again, $t_0 := 4$. The solution in full space reads

$$ u_{\infty}(t, x) := \frac{1}{4} \left( e^{\frac{1}{2}(t-t_0+x)^2} \int_{-\infty}^{t-t_0+x} e^{-\frac{1}{2}\tau^2} d\tau + e^{\frac{1}{2}(t-t_0-x)^2} \int_{-\infty}^{t-t_0-x} e^{-\frac{1}{2}\tau^2} d\tau \right) e^{-(x^2+(t-t_0)^2)}. \quad (6.12) $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Time profiles in the standing wave example</figcaption></figure>

This expression is obtained from the representation

$$u_{\infty}(t, x) = \frac{1}{2} \int_{y \in \mathbb{R}} \int_{\theta > |y|} f(t - \theta, x - y) d\theta dy, \quad (6.13)$$

which stems from the Green function $\frac{1}{2}1_{t>|x|}$ of the one-dimensional wave equation in free space. Notice that $u_{\infty}$ does not satisfy the Dirichlet boundary conditions at $x = \pm L$. However, we can obtain the correct solution by using the mirror image principle, and setting

$$u(t, x) := u_{\infty}(t, x) + \sum_{n=1}^{+\infty} (-1)^n \{u_{\infty}(t, 2nL + (-1)^n x) + u_{\infty}(t, -2nL + (-1)^n x)\}. \quad (6.14)$$

One readily checks that $u$ satisfies the boundary conditions at $x \pm L$. As above, we do not have $u(0, \cdot) = 0$, but the initial values are small enough so that they do not affect the numerical experiments. Notice that $u$ in fact solves the wave equation over $J \times \mathbb{R}$ with right-hand side

$$\tilde{f}(t, x) := f(t, x) + \sum_{n=1}^{+\infty} (-1)^n \{f(t, 2nL + (-1)^n x) + f(t, -2nL + (-1)^n x)\}, \quad (6.15)$$

but the difference $\tilde{f} - f$ is small enough over $\Omega$ so that it does not impact the accuracy of our numerical experiments. For any finite time $t \ge 0$ and any $x \in \Omega$, the terms in the series defining $u(t, x)$ decay exponentially. For our simulations, we simply compute the sum over $n \in \{1:50\}$. As above, we use the Faddeeva software package to evaluate the integral of the Gaussian numerically. Space profiles of $u$ at different times are shown in Figure 2.

## 6.3 Accuracy of the reconstructions

Our first goal is to verify that the reconstructions $u_{h\tau}$ and $w_{h\tau}$ obtained from the time-step values $(U^n)_{n \in \mathbb{N}}$ have the expected accuracy. We focus on the standing wave example with $\rho := 0.02$. As can be seen from Figures 3, 4 and 5, all the errors are very close, as anticipated. We also observe the expected convergence rates. We obtained very similar results with different values of $\rho$ as well as in the propagating wave example. For the sake of shortness, we do not reproduce all the curves here.
<figure><img src="image_3.png" /><figcaption>Figure 2: Propagating wave: space profiles of the solution at various times</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Standing wave: reconstruction errors on the solution</figcaption></figure>

## 6.4 Properties of the error estimator

We now focus on the properties of the error estimator, and we start with the standing wave example. We first study short-time error control by setting $\rho := 1$. The corresponding results are reported in Figure 6. The first observation is that the data-oscillation term $\eta_f$ always superconverges. Besides, asymptotically, we indeed obtain a guaranteed error upper bound, and the effectivity index is at most 3 in all cases. Moreover, for $k = 1$, the estimator seems to be asymptotically exact. Since the leapfrog scheme provides a time discretization error decaying as $\tau^2 \sim h^2$, this is to be expected, and it is in agreement with the results observed in the space semi-discrete case in [5]. Similarly, we see that, for $k = 1$, the error estimator $M$ superconverges.

Still for the standing wave example, we now address long-time error control by setting $\rho := 0.02$. We still observe a superconvergence of the data oscillation term $\eta_f$ for all the values of $k$. Asymptotically, the bound given by $\eta$ is indeed guaranteed, and the effectivity index is about 10 at worst. For $k = 1$, we see that $M$ is in fact the dominant part of the estimator, meaning that the asymptotic regime where the time discretization error becomes negligible has not yet been reached. This is in agreement with the above observation regarding the convergence rates.

Interestingly, the cases where $k \ge 2$ show that the $M$ component of the estimator is indeed required to obtain a guaranteed error upper bound. Indeed, in those cases, $R/\mathcal{E}_\rho \ll 1$, even asymptotically. The question of whether the factor $\sqrt{10}$ in the bound is sharp remains open.

We now consider the propagating wave example. Addressing first short-time error control, we set $\rho := 0.2$. Asymptotically, we obtain guaranteed error bounds with effectivity indices of 10 at worst.

Finally, to study long-time error control, we set $\rho := 0.01$. We obtain asymptotically a
<figure><img src="image_5.png" /><figcaption>Figure 4: Standing wave: reconstruction errors on the space derivative</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Standing wave: reconstruction errors on the time derivative</figcaption></figure>

guaranteed error upper bound with effectivity indices of about 10, except for $k = 1$ where the asymptotic regime is not yet reached.

## 6.5 Time-step variation

Finally, we investigate the effect of reducing the time step on fixed meshes for $k \in \{2, 3\}$. To do so, for quadratic elements, we fix $\tau$ using the CFL condition with $r = 0.9, 0.8, 0.5, 0.2$, and $0.1$. For cubic elements, we select $\tau^2 \sim h^3$ as described above with $r_0 = 0.9, 0.8, 0.5, 0.2$, and $0.1$. We set $\rho := 0.05$. As can be seen in Figure 10, reducing the time step has little effect on the actual error itself. However, it does improve the efficiency of the estimator. In particular, we obtain asymptotically constant-free error upper bounds when $\tau$ is selected ten times smaller than the CFL constraint.

## References

[1] G. Akrivis, C. Makridakis, and R. H. Nochetto, *Optimal order a posteriori error estimates for a class of Runge-Kutta and Galerkin methods*, Numer. Math. **114** (2009), no. 1, 133–160.

[2] C. Bernardi and E. Süli, *Time and space adaptivity for the second-order wave equation*, Math. Models Meth. Appl. Sci. **12** (2005), no. 2, 199–225.

[3] P. Bignardi and A. Moiola, *A space-time continuous and coercive formulation for the wave equation*, arXiv:2312.07268, 2023.
<figure><img src="image_7.png" /><figcaption>Figure 6: Standing wave example with $ρ = 1$</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Standing wave example with $\rho = 0.02$</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Propagating wave example with $ρ = 0.2$</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Propagating wave example with $\rho = 0.01$</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Time-step variation in the propagating wave example, $\rho = 0.05$</figcaption></figure>

[4] C. Carstensen and S. A. Funken, *Fully reliable localized error control in the FEM*, SIAM J. Sci. Comput. **21** (1999/00), no. 4, 1465–1484.

[5] T. Chaumont-Frelet, *Asymptotically constant-free and polynomial-degree-robust a posteriori estimates for space discretizations of the wave equation*, 2022, HAL preprint, https://hal.inria.fr/hal-03632468.

[6] P. Ciarlet, Jr. and M. Vohralík, *Localization of global norms and robust a posteriori error control for transmission problems with sign-changing coefficients*, ESAIM Math. Model. Numer. Anal. **52** (2018), no. 5, 2037–2064.

[7] M. Dauge, *Elliptic boundary value problems on corner domains*, Lecture Notes in Mathematics, vol. 1341, Springer-Verlag, Berlin, 1988.

[8] A. Ern and J.-L. Guermond, *Finite elements II: Galerkin approximation, elliptic and mixed pdes*, Texts in Applied Mathematics, vol. 73, Springer Nature, Cham, Switzerland, 2021.

[9] A. Ern, I. Smears, and M. Vohralík, *Guaranteed, locally space-time efficient, and polynomial-degree robust a posteriori error estimates for high-order discretizations of parabolic problems*, SIAM J. Numer. Anal. **55** (2017), no. 6, 2811–2834.

[10] A. Ern and M. Vohralík, *Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations*, SIAM J. Numer. Anal. **53** (2015), no. 2, 1058–1081.

[11] T. Führer, R. González, and M. Karkulik, *Well-posedness of first-order acoustic wave equations and space-time finite element approximation*, arXiv:2311.10536, 2023.
[12] E. H. Georgoulis, O. Lakkis, and C. Makridakis, A posteriori $L^\infty(L^2)$-error bounds for finite element approximations to the wave equation, IMA J. Numer. Anal. 33 (2013), 1245-1264.

[13] E. H. Georgoulis, O. Lakkis, C. 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 (2016), no. 1, 120-136.

[14] H. Gimperlein, C. Özdemir, D. Stark, and E. P. Stephan, A residual a posteriori error estimate for the time-domain boundary element method, Numer. Math. 146 (2020), no. 2, 239-280.

[15] 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 (2019), 1672-1705.

[16] D. Hoonhout, R. Löscher, O. Steinbach, and C. Urzúa-Torres, *Stable least-squares space-time boundary element methods for the wave equation*, arXiv:2312.12547, 2023.

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

[18] C. Makridakis and R. H. Nochetto, *Elliptic reconstruction and a posteriori error estimates for parabolic problems*, SIAM J. Numer. Anal. 41 (2003), no. 4, 1585-1594.

[19] R. Verfürth, A posteriori error estimates for finite element discretizations of the heat equation, Calcolo 40 (2003), 195-212.

[20] ______, A posteriori error estimation techniques for finite element methods, Oxford science publications, 2013.

[21] M. Zank, Inf-sup stable space-time methods for time-dependent partial differential equations, Ph.D. thesis, Technische Universität Graz, 2019.

# 7 Damped stability estimates for leapfrog scheme

This section collects some stability estimates for the leapfrog scheme in the damped energy norm. These estimates are essentially a discrete counterpart of the a priori estimate established in Lemma 2.1 at the continuous level. They are established by adapting known arguments devised for leapfrog stability estimates in the classical energy norm.

## 7.1 Abstract estimate

Let $(G^n)_{n \in \mathbb{N}}$ be a sequence in $L^2(\Omega)$. Let $(X^n)_{n \in \mathbb{N}} \subset V_h^N$ be such that

$$ \frac{1}{\tau^2}(X^{n+1} - 2X^n + X^{n-1}, v_h)_\Omega + (\nabla X^n, \nabla v_h)_\Omega = (G^n, v_h)_\Omega \quad \forall n \ge 1, \forall v_h \in V_h, \quad (7.1) $$

with the initial condition $X^0 = 0$. Notice that here we do not assume that $G^0 = 0$, and consequently we can have $X^1 \ne 0$. For all $n \ge -1$, we define the following midpoint state and velocity:

$$ X^{n+\frac{1}{2}} := \frac{1}{2}(X^{n+1} + X^n), \quad \dot{X}^{n+\frac{1}{2}} := \frac{1}{\tau}(X^{n+1} - X^n), \quad (7.2) $$

as well as the discrete energy functional

$$ E_X^{n+\frac{1}{2}} := m_{h\tau}(\dot{X}^{n+\frac{1}{2}}, \dot{X}^{n+\frac{1}{2}}) + \|\nabla X^{n+\frac{1}{2}}\|_{\Omega}^2, \quad (7.3) $$
with the bilinear form $m_{h\tau}$ defined in (2.13). Recall that the discrete energy functional defines a quadratic form on $\dot{X}^{n+\frac{1}{2}}$ and $\nabla X^{n+\frac{1}{2}}$ under the CFL condition (2.15) since (2.14) implies that

$$E_X^{n+\frac{1}{2}} \geq \mu_0 \|\dot{X}^{n+\frac{1}{2}}\|_{\Omega}^2 + \|\nabla X^{n+\frac{1}{2}}\|_{\Omega}^2. \quad (7.4)$$

with $\mu_0 \in (0, 1)$.

**Lemma 7.1** (Damped stability for leapfrog scheme). Assume the CFL condition (2.15). Let $(X^n)_{n \in \mathbb{N}}$ solve (7.1) and let the discrete energy functional $E_X^{n+\frac{1}{2}}$ be defined in (7.3). Let $\theta := \rho\tau \in (0, 1]$ and set $C_X(\theta) := \frac{169}{12}(13 - 6\theta)^{-1}\theta(1 - e^{-\theta})^{-1}$. The following holds:

$$\sum_{n \in \mathbb{N}} \tau E_X^{n+\frac{1}{2}} e^{-2\rho t^n} \leq C_X(\theta) \mu_0^{-1} \frac{1}{\rho^2} \sum_{n \in \mathbb{N}} \tau \|G^n\|_{\Omega}^2 e^{-2\rho t^n}. \quad (7.5)$$

*Proof.* The starting point is the following well-known energy identity for the leapfrog scheme:

$$E_X^{n+\frac{1}{2}} - E_X^{n-\frac{1}{2}} = (G^n, X^{n+1} - X^{n-1})_{\Omega} = \tau(G^n, \dot{X}^{n+\frac{1}{2}} + \dot{X}^{n-\frac{1}{2}})_{\Omega} \quad \forall n \in \mathbb{N}.$$

Invoking the Cauchy–Schwarz and Young’s inequality (with parameter $\gamma > \frac{1}{4}$) together with (7.4) gives

$$\begin{aligned} \tau|(G^n, \dot{X}^{n+\frac{1}{2}} + \dot{X}^{n-\frac{1}{2}})_\Omega| &\leq \tau \|G^n\|_{\Omega} \|\dot{X}^{n+\frac{1}{2}}\|_{\Omega} + \tau \|G^n\|_{\Omega} \|\dot{X}^{n-\frac{1}{2}}\|_{\Omega} \\ &\leq 2\gamma \mu_0^{-1} \frac{\tau}{\rho} \|G^n\|_{\Omega}^2 + \frac{1}{4\gamma} \rho \tau \mu_0 \|\dot{X}^{n+\frac{1}{2}}\|_{\Omega}^2 + \frac{1}{4\gamma} \rho \tau \mu_0 \|\dot{X}^{n-\frac{1}{2}}\|_{\Omega}^2 \\ &\leq 2\gamma \mu_0^{-1} \frac{\tau}{\rho} \|G^n\|_{\Omega}^2 + \frac{1}{4\gamma} \rho \tau E_X^{n+\frac{1}{2}} + \frac{1}{4\gamma} \rho \tau E_X^{n-\frac{1}{2}}. \end{aligned}$$

Recalling that $\theta := \rho\tau \in (0, 1]$ and re-arranging the terms gives

$$(1 - \frac{1}{4\gamma}\theta)E_X^{n+\frac{1}{2}} \leq 2\gamma\mu_0^{-1}\frac{\tau}{\rho}\|G^n\|_{\Omega}^2 + (1 + \frac{1}{4\gamma}\theta)E_X^{n-\frac{1}{2}}.$$

Setting $c_1(\gamma, \theta) := \frac{2\gamma}{1 - \frac{1}{4\gamma}\theta}$ and $c_2(\gamma, \theta) := \frac{1 + \frac{1}{4\gamma}\theta}{1 - \frac{1}{4\gamma}\theta}e^{-\theta}$, this implies that

$$E_X^{n+\frac{1}{2}} e^{-\rho t^n} \leq c_1(\gamma, \theta) \mu_0^{-1} \frac{\tau}{\rho} \|G^n\|_{\Omega}^2 e^{-\rho t^n} + c_2(\gamma, \theta) E_X^{n-\frac{1}{2}} e^{-\rho t^{n-1}}.$$

The ideal choice of $\gamma$ would be to minimize $c_1(\gamma, 1)$ while ensuring that $c_2(\gamma, \theta) \leq 1$ for all $\theta \in (0, 1]$. A fairly optimal choice is $\gamma = \frac{13}{24}$ (see Remark 7.2 below for some further discussion) giving

$$E_X^{n+\frac{1}{2}} e^{-\rho t^n} \leq c_1(\frac{13}{24}, \theta) \mu_0^{-1} \frac{\tau}{\rho} \sum_{m \in \{0:n\}} \|G^m\|_{\Omega}^2 e^{-\rho t^m},$$

since $E_X^{-\frac{1}{2}} = 0$. Summing over all $n \in \mathbb{N}$ and exchanging the summations on the right-hand side, we infer that

$$\begin{aligned} \sum_{n \in \mathbb{N}} \tau E_X^{n+\frac{1}{2}} e^{-2\rho t^n} &\leq c_1(\frac{13}{24}, \theta) \mu_0^{-1} \frac{\tau^2}{\rho} \sum_{n \in \mathbb{N}} e^{-\rho t^n} \sum_{m \in \{0:n\}} \|G^m\|_{\Omega}^2 e^{-\rho t^m} \\ &= c_1(\frac{13}{24}, \theta) \mu_0^{-1} \frac{\tau^2}{\rho} \sum_{m \geq 0} \left( \sum_{n \geq m} e^{-\rho t^n} \right) \|G^m\|_{\Omega}^2 e^{-\rho t^m} \\ &\leq c_1(\frac{13}{24}, \theta) \theta (1 - e^{-\theta})^{-1} \mu_0^{-1} \frac{1}{\rho^2} \sum_{m \geq 0} \tau \|G^m\|_{\Omega}^2 e^{-2\rho t^m}, \end{aligned}$$

since $\sum_{n \geq m} e^{-\rho t^n} = e^{-\rho t^m} (1 - e^{-\theta})^{-1}$ and $\theta = \rho\tau$. $\square$
**Remark 7.2** (Constant $C_X(\theta)$). We notice that $C_X(\theta) \le \frac{169}{84} \frac{e}{e-1} \approx 3.1828$. Moreover, we observe that $C_X(\theta) \to C_X(0^+) := \frac{13}{12}$ as $\theta \to 0^+$. This limit is relevant as $\theta = \rho\tau$ and $\tau \ll T \le \rho^{-1}$. The exact counterpart of the a priori estimate (2.9) would be $C_X(0^+) = 1$. This would require taking $\gamma = \frac{1}{2}$, but then $c_2(\gamma, \theta)$ can take values larger than one; some further optimization should be possible by considering two distinct coefficients in the Young inequalities related to $\dot{X}^{n+\frac{1}{2}}$ and to $\dot{X}^{n-\frac{1}{2}}$, but this is not further explored here.

## 7.2 Application: estimates on the acceleration

Let $(U^n)_{n \in \mathbb{N}}$ solve the leapfrog scheme (2.12) with $U^1 = U^0 = 0$. Recall the definition (3.10) of the acceleration $A^n$ for all $n \in \mathbb{N}$, and set

$$A^{n+\frac{1}{2}} := \frac{1}{2}(A^{n+1} + A^n), \quad \dot{A}^{n+\frac{1}{2}} := \frac{1}{\tau}(A^{n+1} - A^n) \quad \forall n \ge -1. \quad (7.6a)$$

Similarly, let us set $B^n := \frac{1}{\tau^3}(U^{n+1} - 3U^n + 3U^{n-1} - U^{n-2})$ for all $n \ge \mathbb{N}$, and

$$B^{n+\frac{1}{2}} := \frac{1}{2}(B^{n+1} + B^n), \quad \dot{B}^{n+\frac{1}{2}} := \frac{1}{\tau}(B^{n+1} - B^n) \quad \forall n \ge -1. \quad (7.6b)$$

Notice that $B^{n+1} = \dot{A}^{n+\frac{1}{2}}$.

**Corollary 7.3** (Estimates on acceleration). Assume the CFL condition (2.15). Let $(U^n)_{n \in \mathbb{N}}$ solve (2.12) with $U^1 = U^0 = 0$. Let $A^{n+\frac{1}{2}}$, $\dot{A}^{n+\frac{1}{2}}$ be defined in (7.6a), and let $B^{n+\frac{1}{2}}$, $\dot{B}^{n+\frac{1}{2}}$ be defined in (7.6b). Recall that $\theta := \rho\tau \in (0, 1]$. The following holds:

$$\sum_{n \in \mathbb{N}} \tau (\mu_0 \|\dot{A}^{n+\frac{1}{2}}\|_{\Omega}^2 + \|\nabla A^{n+\frac{1}{2}}\|_{\Omega}^2) e^{-2\rho t^n} \le C_A(\theta) \mu_0^{-1} \frac{1}{\rho^2} \int_0^{+\infty} \|\ddot{f}(t)\|_{\Omega}^2 e^{-2\rho t} dt, \quad (7.7a)$$

$$\sum_{n \in \mathbb{N}} \tau (\mu_0 \|\dot{B}^{n+\frac{1}{2}}\|_{\Omega}^2 + \|\nabla B^{n+\frac{1}{2}}\|_{\Omega}^2) e^{-2\rho t^n} \le C_B(\theta) \mu_0^{-1} \frac{1}{\rho^2} \int_0^{+\infty} \|\dddot{f}(t)\|_{\Omega}^2 e^{-2\rho t} dt, \quad (7.7b)$$

with $C_A(\theta) := \frac{2}{3}C_X(\theta)(1 + e^{2\theta})$, $C_B(\theta) := \frac{11}{20}C_X(\theta)(2 + e^{2\theta})$, and $C_X(\theta)$ defined in Lemma 7.1.

*Proof.* (1) Proof of (7.7a). We observe that the sequence $(A^n)_{n \in \mathbb{N}}$ solves the generic leapfrog scheme (7.1) with right-hand side

$$G^n := \frac{F^{n+1} - 2F^n + F^{n-1}}{\tau^2} = \frac{1}{\tau} \int_{t^{n-1}}^{t^{n+1}} \psi^n(t) \ddot{f}(t) dt,$$

where $\psi^n(t)$ denotes the hat basis function in time having support in $[t^{n-1}, t^{n+1}]$ and satisfying $\psi^n(t^n) = 1$. Since $\int_{t^{n-1}}^{t^{n+1}} \psi_n(t)^2 dt = \frac{2}{3}\tau$, a Cauchy–Schwarz inequality shows that

$$\|G^n\|_{\Omega}^2 \le \frac{2}{3\tau} \int_{t^{n-1}}^{t^{n+1}} \|\ddot{f}(t)\|_{\Omega}^2 dt.$$

Letting $E_A^{n+\frac{1}{2}} := m_{h\tau}(\dot{A}^{n+\frac{1}{2}}, \dot{A}^{n+\frac{1}{2}}) + \|\nabla A^{n+\frac{1}{2}}\|_{\Omega}^2$ and since $m_{h\tau}(\dot{A}^{n+\frac{1}{2}}, \dot{A}^{n+\frac{1}{2}}) \ge \mu_0 \|\dot{A}^{n+\frac{1}{2}}\|_{\Omega}^2$ owing to (2.14), Lemma 7.1 gives

$$\sum_{n \in \mathbb{N}} \tau (\mu_0 \|\dot{A}^{n+\frac{1}{2}}\|_{\Omega}^2 + \|\nabla A^{n+\frac{1}{2}}\|_{\Omega}^2) e^{-2\rho t^n} \le \sum_{n \in \mathbb{N}} \tau E_A^{n+\frac{1}{2}} \le C_X(\theta) \mu_0^{-1} \frac{1}{\rho^2} \sum_{n \in \mathbb{N}} \tau \|G^n\|_{\Omega}^2 \tau e^{-2\rho t^n} \\ \le \frac{2}{3} C_X(\theta) \mu_0^{-1} \frac{1}{\rho^2} \sum_{n \in \mathbb{N}} e^{-2\rho t^n} \int_{t^{n-1}}^{t^{n+1}} \|\ddot{f}(t)\|_{\Omega}^2 dt.$$
Since $\theta = \rho\tau \le 1$, we observe that

$$e^{-2\rho t^n} \int_{t^{n-1}}^{t^{n+1}} \|\ddot{f}(t)\|_{\Omega}^2 dt \le \int_{J_{n-1}} \|\ddot{f}(t)\|_{\Omega}^2 e^{-2\rho t} dt + e^{2\theta} \int_{J_n} \|\ddot{f}(t)\|_{\Omega}^2 e^{-2\rho t} dt.$$

This completes the proof of (7.7a).

(2) The proof of (7.7b) is similar and is only sketched. We observe that the sequence $(B^n)_{n \in \mathbb{N}}$ solves the generic leapfrog scheme (7.1) with right-hand side

$$G^n := \frac{F^{n+1} - 3F^n + 3F^{n-1} - F^{n-2}}{\tau^3} = \frac{1}{\tau^2} \int_{t^{n-2}}^{t^{n+1}} \phi^n(t) \ddot{f}(t) dt,$$

with $\phi_n(t) := \psi_n(t) - \psi_{n-1}(t)$. Since $\phi_n$ is skew-symmetric around $t^{n-\frac{1}{2}} := \frac{1}{2}(t^{n-1} + t^n)$, we infer using Fubini's theorem that

$$\begin{aligned} \int_{t^{n-2}}^{t^{n+1}} \phi^n(t) \ddot{f}(t) dt &= \int_{t^{n-\frac{1}{2}}}^{t^{n+1}} \phi^n(t) (\ddot{f}(t) - \ddot{f}(\hat{t})) dt \\ &= \int_{t^{n-\frac{1}{2}}}^{t^{n+1}} \phi^n(t) \left( \int_{\hat{t}}^{t} \dddot{f}(s) ds \right) dt \\ &= \int_{t^{n-2}}^{t^{n+1}} \dddot{f}(s) \xi(s) ds, \quad \xi(s) := \int_{\hat{s}}^{t^{n+1}} \phi_n(t) dt, \end{aligned}$$

with $\hat{t} := 2t^{n-\frac{1}{2}} - t$ for all $t \in [t^{n-\frac{1}{2}}, t^{n+1}]$ and $\hat{s} := t^{n-\frac{1}{2}} + |s - t^{n-\frac{1}{2}}|$ for all $s \in [t^{n-2}, t^{n+1}]$. A direct calculation shows that

$$\xi(s) = \begin{cases} \int_s^{t^{n+1}} \frac{t^{n+1}-t}{\tau} dt = \frac{\tau}{2} \left(\frac{t^{n+1}-s}{\tau}\right)^2 & s \in [t^n, t^{n+1}], \\ \frac{\tau}{2} + \int_s^{t^n} 2 \frac{t-t^{n-\frac{1}{2}}}{\tau} dt = \frac{3\tau}{4} - \tau \left(\frac{s-t^{n-\frac{1}{2}}}{\tau}\right)^2 & s \in [t^{n-\frac{1}{2}}, t^n]. \end{cases}$$

This implies that

$$\int_{t^{n-2}}^{t^{n+1}} \xi(s)^2 ds = 2 \int_{t^{n-\frac{1}{2}}}^{t^{n+1}} \xi(s)^2 ds = 2 \left( \int_{t^{n-\frac{1}{2}}}^{t^n} \xi(s)^2 ds + \int_{t^n}^{t^{n+1}} \xi(s)^2 ds \right) = \frac{11}{20} \tau^3,$$

since $\int_{t^{n-\frac{1}{2}}}^{t^n} \xi(s)^2 ds = \tau^3 \int_0^{\frac{1}{2}} (\frac{3}{4} - u^2)^2 du = \frac{9\tau^3}{40}$ and $\int_{t^n}^{t^{n+1}} \xi(s)^2 ds = \frac{\tau^3}{4} \int_0^1 u^4 du = \frac{\tau^3}{20}$. Invoking the Cauchy–Schwarz inequality, we infer that

$$\|G^n\|_{\Omega}^2 \le \frac{1}{\tau} \frac{11}{20} \int_{t^{n-2}}^{t^{n+1}} \|\ddot{f}(t)\|_{\Omega}^2 dt.$$

The conclusion is now straightforward. $\square$

# 8 Damped energy-norm a priori error estimates

For $m \in \mathbb{N}$, we say that a function $\phi \in L^2(\Omega)$ belongs to $H^m(\mathcal{T}_h)$ if $\phi|_K \in H^m(K)$ for all $K \in \mathcal{T}_h$, and we introduce the semi-norm $|\phi|_{H^m(\mathcal{T}_h)} := \sum_{K \in \mathcal{T}_h} |\phi|_K|^2_{H^m(K)}$. For dimensional consistency, we set $\kappa := \max(\ell_{\Omega}, \rho^{-1})$. For a semi-norm $|\cdot|_Y$ equipping the space $Y$ composed of functions defined over $\Omega$, we use the shorthand notation $|v|_{L^2_Y}^2 := \int_0^{+\infty} |v(t)|_Y^2 e^{-2\rho t} dt$. Let $\pi_h^E: V \to V_h$ denote the Riesz elliptic projector such that, for all $v \in V$, the discrete function $\pi_h^E(v) \in V_h$ is uniquely defined by requiring that $(\nabla(\pi_h^E(v) - v), \nabla w_h)_\Omega = 0$ for all $w_h \in V_h$. Recall that $\|\pi_h^E(v) - v\|_{\Omega} \lesssim \ell_{\Omega} h^k \kappa |v|_{H^{k+1}(\mathcal{T}_h)}$ and $\|\nabla(\pi_h^E(v) - v)\|_{\Omega} \lesssim h^k |v|_{H^{k+1}(\mathcal{T}_h)}$. Moreover, we have the $H^1$-stability properties $\|\pi_h^E(v)\|_{\Omega} \lesssim \ell_{\Omega} |v|_V$ and $|\pi_h^E(v)|_V \le |v|_V$.
**Theorem 8.1** (Damped energy-norm a priori error estimate on $(u - u_{h\tau})$). Assume the CFL condition (2.15). Assume that $C_\tau(u) := \kappa|u^{(3)}|_{L^2_\rho(V)} + \kappa^2|u^{(4)}|_{L^2_\rho(V)}$ and $C_h(u) := |u|_{L^2_\rho(H^{k+1}(\tau_h))} + \kappa|\dot{u}|_{L^2_\rho(H^{k+1}(\tau_h))} + \kappa^2|u^{(2)}|_{L^2_\rho(H^{k+1}(\tau_h))} + \kappa^3|u^{(3)}|_{L^2_\rho(H^{k+1}(\tau_h))}$ are bounded. The following holds:

$$ \mathcal{E}_{\rho}^{2}(u - u_{h\tau}) \lesssim \tau^{4}C_{\tau}(u)^{2} + h^{2k}C_{h}(u)^{2}. \qquad (8.1) $$

*Proof.* We consider the sequence $u_h := (u_h^n)_{n \in \mathbb{N}}$ such that $u_h^n := \pi_h^E(u(t^n))$ for all $n \in \mathbb{N}$. We write

$$ u - u_{h\tau} = (u - \pi_h^E(u)) + (\pi_h^E(u) - R(u_h)) + R(u_h - U), $$

and bound the damped energy norm of the three terms on the right-hand side.

(1) The approximation properties of the elliptic projector and $\ell_\Omega \le \kappa$ give

$$ \mathcal{E}_{\rho}(u - \pi_h^E(u)) \lesssim h^k (\kappa |\dot{u}|_{L^2_\rho(H^{k+1}(\tau_h))} + |u|_{L^2_\rho(H^{k+1}(\tau_h))}). $$

(2) Invoking Lemma 3.7 and the $H^1$-stability of the elliptic projection gives

$$ \mathcal{E}_{\rho}(\pi_h^E(u) - R(u_h)) \lesssim \tau^2 (\ell_{\Omega}|u^{(3)}|_{L^2_{\rho}(V)} + \tau|u^{(3)}|_{L^2_{\rho}(V)}) \le \tau^2 \kappa |u^{(3)}|_{L^2_{\rho}(V)}, $$

since $\tau \le \rho^{-1} \le \kappa$.

(3) We observe that $X := u_h - U$ solves the leapfrog scheme with right-hand side

$$ G^n := G_1^n + G_2^n := \pi_h^E((\delta_\tau^2 u)^n - \ddot{u}(t^n)) + (\pi_h^E(\ddot{u}(t^n)) - \ddot{u}(t^n)) \quad \forall n \in \mathbb{N}, $$

with $(\delta_\tau^2 u)^n := \frac{1}{\tau^2}(u(t^{n+1}) - 2u(t^n) + u(t^{n-1}))$. Combining the results of Lemma 3.6 and Lemma 7.1 gives

$$ \mathcal{E}_{\rho}^{2}(R(u_h - U)) \lesssim \sum_{n \in \mathbb{N}} \tau E_{X}^{n+\frac{1}{2}} e^{-2\rho t^n} \lesssim \frac{1}{\rho^2} \sum_{n \in \mathbb{N}} \tau (\|G_1^n\|_{\Omega}^2 + \|G_2^n\|_{\Omega}^2) e^{-2\rho t^n}. $$

Proceeding as in the proof of Lemma 3.7 and using the $H^1$-stability of the elliptic projection, we infer that $\|G_1^n\|_{\Omega}^2 \lesssim \tau^3 \ell_{\Omega}^2 \int_{t^{n-1}}^{t^{n+1}} |u^{(4)}(t)|_Y^2 dt$. Moreover, using the standard estimate $\tau|\psi(t^n)|_Y^2 \lesssim \int_{J_n} |\psi(t)|_Y^2 dt + \tau^2 \int_{J_n} |\dot{\psi}(t)|_Y^2 dt$, we infer that

$$ \tau \|G_2^n\|_{\Omega}^2 \lesssim h^{2k} \ell_{\Omega}^2 \int_{J_n} |u^{(2)}(t)|_{H^{k+1}(\tau_h)}^2 dt + h^{2k} \ell_{\Omega}^2 \tau^2 \int_{J_n} |u^{(3)}(t)|_{H^{k+1}(\tau_h)}^2 dt. $$

Putting the above two bounds together and taking the square root gives

$$ \mathcal{E}_{\rho}(R(u_h - U)) \lesssim \kappa^2 \tau^2 |u^{(4)}|_{L^2_{\rho}(V)} + \kappa^2 h^k (|u^{(2)}|_{L^2_{\rho}(H^{k+1}(\tau_h))} + \kappa |u^{(3)}|_{L^2_{\rho}(H^{k+1}(\tau_h))}). $$

This completes the proof. $\square$

**Corollary 8.2** (Damped energy-norm a priori error estimate on $(u - w_{h\tau})$). Under the assumptions of Theorem 8.1, and provided $C_\tau(f) := \|\dot{f}\|_{L^2_\rho(L^2)} + \kappa \|\ddot{f}\|_{L^2_\rho(L^2)}$ is bounded, the following holds

$$ \mathcal{E}_{\rho}^{2}(u - w_{h\tau}) \lesssim \tau^{4}(C_{\tau}(u)^{2} + \kappa^{2}C_{\tau}(f)^{2}) + h^{2k}C_{h}(u)^{2}. \qquad (8.2) $$

*Proof.* Combine Theorem 8.1 with Lemma 3.8. $\square$

**Remark 8.3** (Sharper a priori estimates). At several places in the above proof, we used suboptimal estimates, like $\tau \le \kappa$. Sharper estimates can be derived. They lead to the same convergence rates, but with higher powers of $\tau$ multiplying some of the higher time-derivatives of $u$ (and $f$).