# 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, 2025, 59 (4), pp.1937-1972. 10.1051/m2an/2025027. hal-05169099

HAL Id: hal-05169099

https://hal.science/hal-05169099v1

Submitted on 18 Jul 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# DAMPED ENERGY-NORM A POSTERIORI ERROR ESTIMATES USING $C^2$-RECONSTRUCTIONS FOR THE FULLY DISCRETE WAVE EQUATION WITH THE LEAPFROG SCHEME

THÉOPHILE CHAUMONT-FRELET<fnref n="1" />,* AND ALEXANDRE ERN<fnref n="2" />,<fnref n="3" />

**Abstract.** We derive a posteriori error estimates for the 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.

**Mathematics Subject Classification.** 35L05, 65M15, 65M20, 65M60.

Received March 20, 2024. Accepted March 31, 2025.

## 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,

<footnote>Keywords and phrases. A posteriori error estimates, finite element method, leapfrog scheme, time-integration, wave equation.</footnote>

<footnote n="1">¹ Univ. Lille, CNRS, Inria, UMR 8524 - Laboratoire Paul Painlevé, F-59000 Lille, France.</footnote>

<footnote n="2">² CERMICS, ENPC, Institut Polytechnique de Paris, 6 & 8 avenue B. Pascal, 77455 Marne-la-Vallée, France.</footnote>

<footnote n="3">³ Inria, 48 Rue Barrault, 75647 Paris, France.</footnote>

<footnote>*Corresponding author: theophile.chaumont@inria.fr</footnote>
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 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., [8, 9, 20, 22, 23] and the references therein). One important challenge 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 ([24], Thm. 4.2.23). Advances towards inf-sup stable variational formulations of (1.1) have been recently reported in [4, 10, 24]. Moreover, in the context of boundary integral equations, we refer the reader to [17] 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 [3, 12, 15] using an implicit time discretization, based on either 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 [3], 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 [3] 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 only in a time semi-discrete setting. However, explicit time semi-discrete schemes cannot propagate compactly supported initial data beyond the initial support. We also notice that explicit time semi-discrete solutions must sit in the domain of iterated powers of the Laplacian operator, an assumption that is readily met at the space discrete level. Finally, we mention the recent preprint [16] where an *a posteriori* error upper bound is derived for the leapfrog scheme with mesh changes.

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 fact that extra time-regularity is required is somehow related to the lack of inf-sup stability). 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} \left\{ \|\dot{e}(t)\|_{\Omega}^{2} + \|\nabla e(t)\|_{\Omega}^{2} \right\} e^{-2\rho t} \mathrm{d}t, \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_*$. Two important remarks are in order. First, assuming that the source term satisfies $f \in L^2(J; L^2(\Omega))$, classical arguments (briefly recalled in Rem. 2.2) show that the undamped energy of the exact solution, $E_u(t) := \frac{1}{2}\|\dot{u}(t)\|_{\Omega}^2 + \frac{1}{2}\|\nabla u(t)\|_{\Omega}^2$, grows at most linearly with $t$. Under the slightly tighter assumption $f \in H^1(J; L^2(\Omega))$ and the above CFL restriction on the time step, we verify in Remark 3.7 that the same property holds true for the undamped energy of the (time reconstructed) discrete solution (see below for its precise definition). This justifies neglecting the tail of the time integral in (1.2) in practical computations. Second, although the leapfrog scheme enjoys uniform-in-time stability, this does not lead to uniform-in-time error estimates that decay to zero with the mesh size (even in the classical *a priori* error
analysis based on undamped energy arguments). Thus, the introduction of a time horizon, here by means of the damping parameter $\rho$ in (1.2), is a natural ingredient in the a priori and a posteriori error analysis.

The main step forward accomplished herein is to extend [5] to the fully discrete setting using the leapfrog scheme in time with a constant time step and a fixed space discretization. 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 Sect. 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.) as the mesh size and the time step tend to zero,

$$ \mathcal{E}_{\rho}^2(e) \le \int_0^{+\infty} \left\{ \eta_h^2(t) + \frac{20}{\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 w_{h\tau}(t) - f_{\tau}(t) \|_{H^{-1}(\Omega)}. \quad (1.8) $$

Notice that the right-hand side of (1.7) is the sum of two terms, which may be linked, respectively, to the space and the time discretization errors, plus higher-order terms further discussed in Remark 4.5. We observe that all the higher-order terms are either explicitly computable or admit a computable upper bound (which is not of higher-order, but may be useful to certify the error). 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 a triangle inequality producing the term $\delta_{h\tau}$ and 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) \mathbf{e}^{-2\rho t} dt \lesssim \mathcal{E}_\rho^2(e) + \int_0^{+\infty} \|\nabla \delta_{h\tau}(t)\|^2 \mathbf{e}^{-2\rho t} dt + \text{h.o.t.}, \quad (1.9) $$

with higher-order terms further discussed in Remark 5.2. Thus, if the term

$$ 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\} \mathbf{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) + D_{\rho}^{2}$ and $\tilde{\Lambda}_{\rho}^{2} := \int_{0}^{+\infty} \eta_{h}^{2}(t) \mathbf{e}^{-2\rho t} dt + 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 $D_{\rho}$ in the error and the estimator is standard in other contexts, such as discontinuous Galerkin discretizations of elliptic problems [19]. Notice that we prove in Lemma 3.9 below that $D_{\rho}$ decreases optimally in time (i.e., as $\tau^2$). In our numerical experiments on the wave equation, we observe that the right-hand side of (1.7) controls $\mathcal{E}_{\rho}$ by a factor of at most 10.

We notice that the present limitation of a constant time step and a fixed space discretization is not easy to lift. This is not specific to the present a posteriori analysis since the obstruction already arises in the a priori error analysis. For instance, stability issues may arise if the time step varies, as discussed in [21]. Moreover, mesh changes can, in principle, be taken into account in the error upper bound by including in the estimator the energy variations induced by the mesh changes, but the question of incorporating these changes in an error lower bound remains open.

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 Section 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. For completeness, we also derive an a priori error estimate on the exact solution in the damped energy norm. For a priori error estimates in the classical energy norm for the leapfrog scheme, we refer the reader, e.g., to [18] and the references therein.

## 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) | v^{(r)} \in L^\infty(J; Y)\}$ and $C_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 H^1(J; L^2(\Omega)) \cap L^\infty(J; L^2(\Omega))$ so that, in particular, $f \in C_b^0(J; L^2(\Omega))$. Moreover, we assume that $f$ is supported away from zero. It is then natural to seek for a strong solution $u \in C_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^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, \quad \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}. \quad (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}. \quad (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}, \quad (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}. \quad (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) \le \frac{1}{\rho^2} \int_{0}^{+\infty} \|f(t)\|_{\Omega}^{2} e^{-2\rho t} dt. \quad (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} \le \|\hat{f}(s)\|_{\Omega} \|\hat{u}(s)\|. $$

Hence, we have

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

The assertion follows from $s\hat{u}(s) = \hat{u}(s)$ and (2.4). □
**Remark 2.2** (Bound on undamped energy). The undamped energy $E_u(t) := \frac{1}{2}\|\dot{u}(t)\|_{\Omega}^2 + \frac{1}{2}\|\nabla u(t)\|_{\Omega}^2$ satisfies the identity $\dot{E}_u(t) = (f(t), \dot{u}(t))_{\Omega}$ at all times. Fix any time $t > 0$. Then $\dot{E}_u(s) \le \frac{t}{2}\|f(s)\|_{\Omega}^2 + \frac{1}{t}E_u(s)$ for all $s \in (0, t)$, and a Gronwall-like argument readily gives $E_u(t) \le \frac{e}{2}t \int_0^t \|f(s)\|_{\Omega}^2 ds$. This shows that $E_u(t)$ grows at most linearly in $t$ if $f \in L^2(J; L^2(\Omega))$.

## 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 \ge 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}$). As motivated in the introduction, 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 \ge 1$, where $T$ is some observation time-scale, we make in what follows the following mild assumption:

$$ \rho\tau \le 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^{\mathbb{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 \ge 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.3** (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 $\mathcal{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) = |\hat{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, \quad (2.16) $$

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

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

$$ \gamma_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\}, \quad (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 $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 [6], 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

$$ \begin{aligned} \|\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), \end{aligned} $$

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). \quad (2.18) $$

Bounds on $\gamma_{\rho, \omega}(h)$ are readily derived from Lemma 2.4. 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) | v_{h\tau}|_{J_n} \in \mathbb{P}^\ell(J_n; Y) \forall n \in \mathbb{N}\}, \quad (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 + \vartheta^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{3V^{n+1} + 17V^n + 5V^{n-1} - V^{n-2}}{24}, \qquad \beta^n = \frac{5V^{n+1} + 3V^n - 9V^{n-1} + V^{n-2}}{12\tau}, \qquad (3.6a)$$

$$\gamma^n = \frac{V^{n+1} - 2V^n + V^{n-1}}{\tau^2}, \qquad \vartheta^n = \frac{V^{n+2} - 2V^{n+1} + 2V^{n-1} - V^{n-2}}{2\tau^3}, \qquad (3.6b)$$

$$\epsilon^n = \frac{V^{n+2} - 4V^{n+1} + 6V^n - 4V^{n-1} + V^{n-2}}{\tau^4}. \qquad (3.6c)$$

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

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

*Proof.* Let us set $v := L(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 + \vartheta^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 + \vartheta^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 + \vartheta^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 $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(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)$, $\vartheta^n = \dddot{v}(t^n) + O(\tau^2)$, and $\epsilon^n = \ddddot{v}(t^n) + O(\tau^2)$.

Consider a sequence $V := (V^n)_{n \in \mathbb{N}}$ in $Y$ such that $V^0 = 0$. Define the sequence $(D_{\tau}^2 V)^n_{n \in \mathbb{N}}$ such that $D_{\tau}^2 V^n := \frac{1}{\tau^2}(V^{n+1} - 2V^n + V^{n-1})$, for all $n \in \mathbb{N}$. (It is legitimate to set $D_{\tau}^2 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(V)(t) = R(D_{\tau}^2 V)(t) \quad \forall t \in J. \qquad (3.7) $$

*Proof.* Since $L(V)$ is of class $C^2$ and $R(D_{\tau}^2 V)$ is of class $C^0$ owing to Lemmas 3.2 and 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$, $\vartheta^n$, and $\epsilon^n$ defining $L(V)$ satisfy

$$ \gamma^n = D_{\tau}^2 V^n, \quad \vartheta^n = \frac{D_{\tau}^2 V^{n+1} - D_{\tau}^2 V^{n-1}}{2\tau}, \quad \epsilon^n = \frac{D_{\tau}^2 V^{n+1} - 2D_{\tau}^2 V^n + D_{\tau}^2 V^{n-1}}{\tau^2}. $$

As a consequence, we have

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

This concludes the proof. $\square$

## 3.2. Rewriting of the 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(F) \in P^2(J_{\tau}; L^2(\Omega)) \cap C^0(J; L^2(\Omega)) \quad \text{with} \quad F^n := f(t^n) \quad \forall n \in \mathbb{N}. \qquad (3.8) $$

Notice that the sequence $(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(U) \in P^2(J_{\tau}; V_h) \cap C^0(J; V_h), \quad w_{h\tau} := L(U) \in P^4(J_{\tau}; V_h) \cap C^2(J; V_h), \qquad (3.9) $$

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

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

Notice that $A^0 = 0$ by definition and that it is legitimate to set $A^{-1} = 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. \qquad (3.11) $$

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

$$ (A^n, v_h)_\Omega + (\nabla U^n, \nabla v_h)_\Omega = (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 $R$ to this equation and invoking Lemma 3.4 which gives $R(A)(t) = \frac{d^2}{dt^2} L(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 $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 $V := (V^n)_{n \in \mathbb{N}}$ be a bounded sequence in $Y$ and set $\dot{V}^{n+\frac{1}{2}} := \frac{1}{\tau}(V^{n+1} - V^n)$. The following holds:

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

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

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

$$ \int_{0}^{+\infty} \|R(V)\|_{Y}^{2} e^{-2\rho t} dt \lesssim \sum_{n \in \mathbb{N}} \tau \max_{\delta \in \{-1, 0, 1\}} \|V^{n+\delta}\|_{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(V)\|_{L^{\infty}(J_n; Y)} \lesssim \max(\|\dot{V}^{n-\frac{1}{2}}\|_Y, \|\dot{V}^{n+\frac{1}{2}}\|_Y)$. $\square$

**Remark 3.7 (Bound on undamped energy for time-reconstructed solutions).** For all $n \in \mathbb{N}$, we define the midpoint state and velocity, $U^{n+\frac{1}{2}} := \frac{1}{2}(U^{n+1} + U^n)$ and $\dot{U}^{n+\frac{1}{2}} := \frac{1}{\tau}(U^{n+1} - U^n)$, as well as the discrete undamped energy $E_{U}^{n+\frac{1}{2}} := \frac{1}{2}\|\dot{U}^{n+\frac{1}{2}}\|_{\Omega}^{2} + \frac{1}{2}\|\nabla U^{n+\frac{1}{2}}\|_{\Omega}^{2}$. Classical arguments for the leapfrog scheme (see [18] and Sect. 7.1 for a global-in-time bound on the damped energy) show that, under the CFL condition (2.15), we have

$$ E_{U}^{n+\frac{1}{2}} \lesssim t^n \sum_{m \in \{0:n\}} \tau \|F^n\|_{\Omega}^2 \lesssim t^n \|f\|_{H_{\tau}^1(0,t^n;L^2(\Omega))}^2, \quad (3.13) $$

with

$$ \|f\|_{H_{\tau}^{1}(0, t^{n}; L^{2}(\Omega))}^{2} := \|f\|_{L^{2}(0, t^{n}; L^{2}(\Omega))}^{2} + \tau^{2} \|\dot{f}\|_{L^{2}(0, t^{n}; L^{2}(\Omega))}^{2} $$

and where the second bound follows from the embedding $H^1(J_n; L^2(\Omega)) \hookrightarrow C^0(\bar{J}_n; L^2(\Omega))$. Define now, for all $t \in J$, the undamped energies $\mathbb{E}_{u_{h\tau}}(t) := \frac{1}{2}\|\dot{u}_{h\tau}(t)\|_{\Omega}^2 + \frac{1}{2}\|\nabla u_{h\tau}(t)\|_{\Omega}^2$ and $\mathbb{E}_{w_{h\tau}}(t) := \frac{1}{2}\|\dot{w}_{h\tau}(t)\|_{\Omega}^2 + \frac{1}{2}\|\nabla w_{h\tau}(t)\|_{\Omega}^2$. Owing to (3.3), we infer that, for all $n \in \mathbb{N}$,

$$ \|\dot{u}_{h\tau}\|_{L^{\infty}(J_n; L^2(\Omega))} \lesssim \max(\|\dot{U}^{n+\frac{1}{2}}\|_{\Omega}, \|\dot{U}^{n-\frac{1}{2}}\|_{\Omega}), $$

and

$$ \begin{aligned} \|\nabla u_{h\tau}\|_{L^{\infty}(J_n; L^2(\Omega))} & \lesssim \max(\|\nabla U^{n+1}\|_{\Omega}, \|\nabla U^n\|_{\Omega}, \|\nabla U^{n-1}\|_{\Omega}) \\ & \lesssim \max(\|\nabla U^{n+\frac{1}{2}}\|_{\Omega}, \|\nabla U^{n-\frac{1}{2}}\|_{\Omega}) + \frac{\tau}{h} \max(\|\dot{U}^{n+\frac{1}{2}}\|_{\Omega}, \|\dot{U}^{n-\frac{1}{2}}\|_{\Omega}), \end{aligned} $$

where the second bound follows from elementary algebraic manipulations and an inverse inequality in space. Combining these bounds and invoking the CFL condition readily shows that $\|\mathbb{E}_{u_{h\tau}}\|_{L^{\infty}(J_n)} \lesssim t^n \|f\|_{H_{\tau}^1(0, t^n; L^2(\Omega))}^2$. Similar arguments using (3.5) and (3.6) (omitted for brevity) show that $\|\mathbb{E}_{w_{h\tau}}\|_{L^{\infty}(J_n)} \lesssim t^n \|f\|_{H_{\tau}^1(0, t^n; L^2(\Omega))}^2$. In conclusion, the undamped energies of the time-reconstructed solutions, $\mathbb{E}_{u_{h\tau}}(t)$ and $\mathbb{E}_{w_{h\tau}}(t)$, grow at most linearly in $t$ if $f \in H^1(J; L^2(\Omega))$.
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.8 (Approximation).** The following holds for every function $v \in C_b^3(J; Y)$ supported away from zero:

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

$$ \int_0^{+\infty} \left\| \frac{d}{dt} (v - R(V))(t) \right\|_Y^2 e^{-2\rho t} dt \lesssim \tau^4 \int_0^{+\infty} \|\ddot{v}(t)\|_Y^2 e^{-2\rho t} dt. \quad (3.14b) $$

*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.14a). We observe that

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

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(V))(t) = \Xi_1^n (t - t^n) + \Xi_2^n \frac{1}{2} (t - t^n)^2 + \Xi_3(t), \quad (3.16) $$

with

$$ \Xi_1^n := \dot{V}^n - \frac{V^{n+1} - V^{n-1}}{2\tau}, \quad \Xi_2^n := \ddot{V}^n - \frac{V^{n+1} - 2V^n + V^{n-1}}{\tau^2}, \quad \Xi_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^n} (s - t^{n-1}) \ddot{v}(s) ds. $$

We infer that

$$ \Xi_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}(\tilde{s})) ds, $$

with $\tilde{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}(\tilde{s})\|_Y \le \int_{\tilde{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

$$ \|\Xi_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

$$ \Xi_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}(\tilde{s}) ds, $$

so that

$$ \|\Xi_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

$$ \|\Xi_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.16) with the above bounds gives

$$ \int_{J_n} \|(v - R(\mathcal{V}))(t)\|_Y^2 dt \lesssim \tau \left( \tau^2 \|\Xi_1^n\|_Y^2 + \tau^4 \|\Xi_2^n\|_Y^2 + \|\Xi_3\|_{C^0(J_n; Y)}^2 \right) \lesssim \tau^6 \int_{t^{n-1}}^{t^{n+1}} \|\ddot{v}(\sigma)\|_Y^2 d\sigma. $$

Using this estimate in (3.15) gives

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

and the assertion follows from

$$ e^{-2\rho t^n} \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.14b) is similar and is only sketched. We observe that for all $t \in J_n$,

$$ \frac{d}{dt}(v - R(\mathcal{V}))(t) = \Xi_1^n + \Xi_2^n(t - t^n) + \Xi_4(t), $$

with $\Xi_4(t) := \int_{t^n}^t (t-s)\ddot{v}(s) ds$. Using the above bounds on $\|\Xi_1^n\|_Y$ and $\|\Xi_2^n\|_Y$ together with $\|\Xi_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.14b). $\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.17) $$

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.9** (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.18a) $$

$$ \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.18b) $$

$$ \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.18c) $$

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.18d) $$
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!} \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.19) $$

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.20a) $$

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

$$ \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.20c) $$

We can now invoke Corollary 7.3 (see Sect. 7) to conclude the proof. Indeed, (3.18a) and (3.18d) readily follow from (7.7a), (3.19), and (3.20a). The estimate (7.7a) can also be used to bound the first term on the right-hand side of (3.20b). To bound the second term on the right-hand side of (3.20b) as well as the right-hand side of (3.20c), 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}{\tau}(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}=1} |\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.17) the time-reconstruction error defined as $\delta_{h\tau}(t) := u_{h\tau}(t) - w_{h\tau}(t)$. In the proofs below, it is useful to consider the modified residual $\mathcal{R}^{\dagger}(t) \in V'$ such that $\mathcal{R}^{\dagger}(t) := \mathcal{R}(t) + \Delta\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 (1 + 40\gamma_s(h)^2) \|\hat{\mathcal{R}}(s)\|_{V'}^2 + \frac{18}{\rho^2} \|\nabla \hat{\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 \hat{\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} \left| (\hat{\eta}_f(s), \hat{\chi}_e(s))_{\Omega} + (\nabla \hat{\delta}_{h\tau}(s), \nabla \hat{\chi}_e(s))_{\Omega} \right| &\le (\|\hat{\eta}_f(s)\|_{\Omega}^2 + \|\nabla \hat{\delta}_{h\tau}(s)\|_{\Omega}^2)^{\frac{1}{2}} \frac{1}{|s|} \|\hat{\chi}_e(s)\| \\ &\le (\|\hat{\eta}_f(s)\|_{\Omega}^2 + \|\nabla \hat{\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} \left| \langle \hat{\mathcal{R}}^{\dagger}(s), \hat{\chi}_e(s) - \hat{v}_h(s) \rangle \right| \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'}. $$

Since $\|\hat{\mathcal{R}}^{\dagger}(s)\|_{V'} \le \|\hat{\mathcal{R}}(s)\|_{V'} + \|\nabla \hat{\delta}_{h\tau}(s)\|_{\Omega}$ by the triangle inequality and since $\gamma_s(h) \le \frac{|s|}{\rho}$ owing to Lemma 2.4, we infer that

$$ |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}}(s)\|_{V'} + \frac{1}{\rho} \|\nabla \hat{\delta}_{h\tau}(s)\|_{\Omega}^2 =: \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(5\gamma_s(h)^2 \|\hat{\mathcal{R}}(s)\|_{V'}^2 + \frac{5}{4}\frac{1}{\rho^2} \|\nabla \hat{\delta}_{h\tau}(s)\|_{\Omega}^2)$ and re-arranging the terms. $\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^r(J; L^2(\Omega))$:

$$ \|\hat{e}(s)\| \le \frac{1}{\rho} \left( \|\nabla \hat{\delta}_{h\tau}(s)\|_{\Omega}^2 + \|\hat{\eta}_f(s)\|_{\Omega}^2 \right)^{\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 \left( (\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} \right) \\ &\le \left\{ \left( \|\hat{f}(s)\|_{\Omega} + \|\hat{\eta}_f(s)\|_{\Omega} \right)^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} \|\widehat{f^{(r)}}(s)\|_{\Omega} + \frac{1}{\rho} \left( \|\hat{\eta}_f(s)\|_{\Omega}^2 + \|\nabla \hat{\delta}_{h\tau}(s)\|_{\Omega}^2 \right)^{\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.17), respectively, and let the residual $\mathcal{R}(t) \in V'$ be defined in (4.3). 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:

$$ \mathcal{E}_{\rho}^{2}(e) \le \int_{0}^{+\infty} \left\{ \|\mathcal{R}(t)\|_{V'}^{2} + \frac{20}{\rho^{2}} \|\nabla \hat{\delta}_{h\tau}(t)\|_{\Omega}^{2} \right. \\ \left. \qquad + 40\gamma_{\rho,\omega}(h)^{2} \|\mathcal{R}(t)\|_{V'}^{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. \quad (4.11) $$

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

$$ \mathcal{E}_{\rho}^{2}(e) = \int_{0}^{+\infty} \left\{ \|\dot{e}(t)\|_{\Omega}^{2} + \|\nabla e(t)\|_{\Omega}^{2} \right\} 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\{ (1 + 40\gamma_{\rho,\omega}(h)^{2}) \|\hat{\mathcal{R}}(s)\|_{V'}^{2} + \frac{18}{\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}} \|\widehat{f^{(r)}}(s)\|_{\Omega}^2 \right\} ds. $$

Putting the above two bounds together proves the assertion.

**Remark 4.4** (Thm. 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$.

**Remark 4.5** (Higher-order terms). The error upper bound (4.11) fits the form (1.7) if one sets

$$ \text{h.o.t.} := \int_{0}^{+\infty} \left\{ 40\gamma_{\rho,\omega}(h)^2 \|\mathcal{R}(t)\|_{V'}^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. $$

Let us motivate that the three terms on the right-hand side can be considered as higher-order terms. First, we observe that the data time-oscillation term $\eta_f$ converges to third-order in $\tau$ owing to Lemma 3.8 provided $f \in C_b^3(J; L^2(\Omega))$. Concerning the other two terms on the right-hand side, we adapt the arguments of Corollary 5.2 from [5]. We first set the cutoff frequency as

$$ \omega^2 = \left( \left( \frac{\ell_{\Omega}}{h} \right)^{\frac{\theta}{2}} (\rho \ell_{\Omega})^{\beta} - 1 \right) \rho^2, $$

where the exponent $\beta \ge 0$ will be chosen later on (the above right-hand side is positive if $h$ is small enough). Owing to the bound (2.17) on $\gamma_s(h)$, we obtain

$$ \gamma_{\rho,\omega}(h) \le 2C_{\text{app}}C_{\text{ell}}(\rho\ell_{\Omega})^{1-\beta} \left(\frac{h}{\ell_{\Omega}}\right)^{\frac{\theta}{2}}. \quad (4.12) $$

Assuming (to fix the ideas) that $h \le \ell_{\Omega}/16$, $\rho\ell_{\Omega} \le 1$, and since $\theta > \frac{1}{2}$, we infer that $\omega^2 \ge \frac{1}{2}\rho^2(\frac{\ell_{\Omega}}{h})^{\frac{\theta}{2}}(\rho\ell_{\Omega})^{-\beta}$. Observing that

$$ \frac{8}{\rho^2} \frac{1}{\omega^{2r}} \|f^{(r)}(t)\|_{\Omega}^2 \le \frac{2^{3+r}}{\rho^{2r+2}} (\rho\ell_{\Omega})^{\beta r} \left(\frac{h}{\ell_{\Omega}}\right)^{\frac{r\theta}{2}} \|f^{(r)}(t)\|_{\Omega}^2, \quad (4.13) $$

we select $r$ so that $r \ge \frac{2k+1}{\theta} \ge 4k+2$. Finally, we can choose $\beta = \frac{3+2r}{1+r} \approx 2+\frac{1}{r} \approx 2$ to balance the powers of $\rho$ in the bounds (4.12) and (4.13).

## 5. ERROR LOWER BOUND

In this section, we establish an error lower bound. We need to consider only the first two terms on the right-hand side of (4.11) since the other three terms can be considered as higher-order terms (see Rem. 4.5). However, the second term on the right-hand side, involving $\nabla \hat{\delta}_{h\tau}$, is, at the same time, an error indicator and an error. Therefore, we focus here on bounding the first term on the right-hand side of (4.11) involving the dual norm of the residual.

We consider the $L^2$-orthogonal projection $\pi_h: L^2(\Omega) \to V_h$ defined by $(\pi_h(\phi), v_h)_\Omega = (\phi, v_h)_\Omega$ for all $\phi \in L^2(\Omega)$ and all $v_h \in V_h$. Under mild assumptions on the mesh grading, which accommodate newest vertex bisection in two dimensions (see [2, 11] and also [7], Rem. 22.23), we have, for all $\phi \in V$,

$$ \|\phi - \pi_h(\phi)\|_{\Omega} \lesssim h \|\nabla \phi\|_{\Omega}, \quad \|\nabla \pi_h(\phi)\|_{\Omega} \lesssim \|\nabla \phi\|_{\Omega}. \quad (5.1) $$
**Theorem 5.1** (Error lower bound). Assume $f \in C_b^r(J; L^2(\Omega))$. Assume the CFL condition (2.15). Assume that the mesh is graded so that (5.1) is satisfied. The following holds for all $r \ge 2$:

$$ \int_0^{+\infty} \|\mathcal{R}(t)\|_{V'}^2 e^{-2\rho t} \, dt \lesssim \int_0^{+\infty} \left\{ \|\dot{e}(t)\|_{\Omega}^2 + \|\nabla e(t)\|_{\Omega}^2 \right\} e^{-2\rho t} \, dt + \int_0^{+\infty} \|\nabla \delta_{h\tau}(t)\|_{\Omega}^2 e^{-2\rho t} \, dt \\ + h^2 \int_0^{+\infty} \|\eta_f(t)\|_{\Omega}^2 e^{-2\rho t} \, dt + \frac{h^{2r}}{\rho^2} \int_0^{+\infty} \|f^{(r)}(t)\|_{\Omega}^2 e^{-2\rho t} \, dt. \quad (5.2) $$

*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). To bound the norm of $\mathcal{R}^\dagger(t)$ in $V'$, we pick an arbitrary $v \in V$ with $\|\nabla v\|_{\Omega} = 1$. Invoking the Galerkin orthogonality property (4.6) gives

$$ \begin{aligned} \langle \mathcal{R}^\dagger(t), v \rangle &= \langle \mathcal{R}^\dagger(t), (I - \pi_h)(v) \rangle \\ &= (\ddot{e}(t), (I - \pi_h)(v))_{\Omega} + (\nabla e(t), \nabla(I - \pi_h)(v))_{\Omega}, + (\nabla \delta_{h\tau}(t), \nabla(I - \pi_h)(v))_{\Omega} - (\eta_f(t), (I - \pi_h)(v))_{\Omega} \\ &= ((I - \pi_h)\ddot{e}(t), (I - \pi_h)(v))_{\Omega} + (\nabla e(t), \nabla(I - \pi_h)(v))_{\Omega}, + (\nabla \delta_{h\tau}(t), \nabla(I - \pi_h)(v))_{\Omega} \\ &\quad - (\eta_f(t), (I - \pi_h)(v))_{\Omega}, \end{aligned} $$

where we employed the $L^2$-orthogonality property of $\pi_h$. Invoking the properties of $\pi_h$ in (5.1), we infer that

$$ \|\mathcal{R}^\dagger(t)\|_{V'} \lesssim h \|(I - \pi_h)\ddot{e}(t)\|_{\Omega} + \|\nabla e(t)\|_{\Omega} + \|\nabla \delta_{h\tau}(t)\|_{\Omega} + h \|\eta_f(t)\|_{\Omega}, $$

for all $t \ge 0$. Owing to the triangle inequality, we infer that

$$ \|\mathcal{R}(t)\|_{V'} \lesssim h \|(I - \pi_h)\ddot{e}(t)\|_{\Omega} + \|\nabla e(t)\|_{\Omega} + \|\nabla \delta_{h\tau}(t)\|_{\Omega} + h \|\eta_f(t)\|_{\Omega}, \quad (5.3) $$

and it remains to bound the first term on the right-hand side.

(2) Owing to the identity (2.4), we infer that

$$ h^2 \int_0^{+\infty} \|(I - \pi_h)(\ddot{e}(t))\|_{\Omega}^2 e^{-2\rho t} \, dt = h^2 \int_{\rho-i\infty}^{\rho+i\infty} |s|^2 \|(I - \pi_h)\hat{e}(s)\|_{\Omega}^2 \, ds. $$

We split the integral on the right-hand side depending on whether $|s| \le h^{-1}$ or $|s| \ge h^{-1}$. On the one hand, we have

$$ h^2 \int_{\rho-i\infty}^{\rho+i\infty} |s|^2 \|(I - \pi_h)\hat{e}(s)\|_{\Omega}^2 1_{|s| \le h^{-1}} \, ds \le \int_{\rho-i\infty}^{\rho+i\infty} \|(I - \pi_h)\hat{e}(s)\|_{\Omega}^2 \, ds = \int_0^{+\infty} \|(I - \pi_h)\dot{e}(t)\|_{\Omega}^2 e^{-2\rho t} \, dt. $$

On the other hand, we have

$$ \begin{aligned} h^2 \int_{\rho-i\infty}^{\rho+i\infty} |s|^2 \|(I - \pi_h)\hat{e}(s)\|_{\Omega}^2 1_{|s| > h^{-1}} \, ds &\le h^2 \int_{\rho-i\infty}^{\rho+i\infty} |s|^2 \|\hat{u}(s)\|_{\Omega}^2 1_{|s| > h^{-1}} \, ds \\ &\le \frac{h^{2r}}{\rho^2} \int_0^{+\infty} \|f^{(r)}(t)\|_{\Omega}^2 e^{-2\rho t} \, dt. \end{aligned} $$

Combining these bounds proves that

$$ h^2 \int_0^{+\infty} \|(I - \pi_h)(\ddot{e}(t))\|_{\Omega}^2 e^{-2\rho t} \, dt \le \int_0^{+\infty} \|(I - \pi_h)\dot{e}(t)\|_{\Omega}^2 e^{-2\rho t} \, dt + \frac{h^{2r}}{\rho^2} \int_0^{+\infty} \|f^{(r)}(t)\|_{\Omega}^2 e^{-2\rho t} \, dt. \quad (5.4) $$

Since $\|(I - \pi_h)\dot{e}(t)\|_{\Omega} \le \|\dot{e}(t)\|_{\Omega}$, the combination of (5.3) and (5.4) concludes the proof. $\square$
**Remark 5.2** (Higher-order terms). The error lower bound (5.2) fits the form (1.9) if one sets

$$ \text{h.o.t.} := h^2 \int_0^{+\infty} \|\eta_f(t)\|_{\Omega}^2 e^{-2\rho t} dt + \frac{h^{2r}}{\rho^2} \int_0^{+\infty} \|f^{(r)}(t)\|_{\Omega}^2 e^{-2\rho t} dt. $$

We refer the reader to Remark 4.5 for the discussion on the data time-oscillation term $\eta_f$. Moreover, we can take here $r \ge k + 1$.

**Remark 5.3** (Alternative error lower bound). It is also possible to establish an error lower bound without invoking the $H^1$-stability of the $L^2$-orthogonal projection, but this leads to the additional term $\frac{h^4}{\rho^2} \int_0^{+\infty} \|\hat{f}(t)\|_{\Omega}^2 e^{-2\rho t} dt$ on the right-hand side of (5.2). This term is not fully satisfactory as it is not of higher-order. To establish the claim, the proof proceeds again in two steps. The first step is similar to the above one (but invokes any $H^1$-stable quasi-interpolation operator instead of the $L^2$-orthogonal projection), leading to (compare with (5.3))

$$ \|\mathcal{R}(t)\|_{V'} \lesssim h \|\ddot{e}(t)\|_{\Omega} + \|\nabla e(t)\|_{\Omega} + \|\nabla \delta_{h\tau}(t)\|_{\Omega} + h \|\eta_f(t)\|_{\Omega}. $$

To bound $h\|\ddot{e}(t)\|_{\Omega}$, we consider as above the Laplace transform. In the low-frequency regime ($|s| \le h^{-1}$), we have

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

In the high-frequency regime ($|s| \ge h^{-1}$), we first invoke the triangle inequality so that $h\|\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 bounded as

$$ 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^r}{\rho} \|\widehat{f^{(r)}(s)}\|_{\Omega}. $$

On the other hand, we have $h|s|^2 \|\hat{w}_{h\tau}(s)\|_{\Omega} \le \|\hat{w}_{h\tau}(s)\|_{\Omega}$, and we invoke the bound (3.18d) from Lemma 3.9 together with the CFL condition. Altogether, this gives

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

whence the claim.

# 6. NUMERICAL RESULTS

In this section, we present numerical results to assess the 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, \qquad e^2(U) := \tau \sum_{n=0}^{+\infty} \|u(t^n) - U^n\|_{\Omega}^2 e^{-2\rho t^n}, \qquad (6.1a)$$

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

$$e_t^2(z) := \int_0^{+\infty} \|\dot{u}(t) - \dot{z}_{h\tau}(t)\|_{\Omega}^2 e^{-2\rho t} \, dt, \qquad 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}. \qquad (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). \qquad (6.2)$$

In practice, all the integrals and sums are computed up to $t = T_\star$. As argued above, $T_\star$ 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. \qquad (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}. \qquad (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, \qquad (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, \qquad (6.5b)$$

and using the following shorthand notation for the estimator:

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

the error upper bound (4.11) rewrites

$$\mathcal{E}_{\rho}^2 \le \Lambda_{\rho}^2 + \text{h. o. t.} \qquad (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)), \qquad (6.8)$$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Time profiles in the standing wave example.</figcaption></figure>

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. \qquad (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. \qquad (6.10)$$

### 6.2.2. Propagating wave

We set

$$f(t, x) := -2(t - t_0)e^{-(t-t_0)^2}e^{-x^2}, \qquad (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)}. \qquad (6.12)$$

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, \qquad (6.13)$$
which stems from the Green function $\frac{1}{2}1_{t>\vert x\vert}$ 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 \left\{ u_{\infty}(t, 2nL + (-1)^n x) + u_{\infty}(t, -2nL + (-1)^n x) \right\}. \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 \left\{ f(t, 2nL + (-1)^n x) + f(t, -2nL + (-1)^n x) \right\}, \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 $(\mathbb{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–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.

## 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$. From the results reported in Figure 7, we again 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$. Results are reported in Figure 8. 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$. Results are reported in Figure 9. We obtain asymptotically a guaranteed error upper bound with effectivity indices of about 10, except for $k = 1$ where the asymptotic regime is not yet reached.
<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.5. Time-step variation

Finally, we investigate the effect of reducing the time step on fixed meshes for $k \in \{2, 3\}$. The goal here is to numerically illustrate that for a small time step, the proposed estimator for the fully discrete scheme approaches the estimator of [5] when the time-discretization error becomes negligible. 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.

## 7. DAMPED STABILITY AND A PRIORI ERROR ESTIMATES

This section collects some stability and error estimates for the leapfrog scheme in the damped energy norm. The stability 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 0, \forall v_h \in V_h, \quad (7.1) $$
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Standing wave: reconstruction errors on the space derivative.</figcaption></figure>

with the initial condition $X^0 = X^{-1} = 0$. Notice that here we do not assume that $G^0 = 0$, and consequently we can have $X^1 \neq 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), \qquad (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, \qquad (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}} \ge \mu_0 \|\dot{X}^{n+\frac{1}{2}}\|_{\Omega}^2 + \|\nabla X^{n+\frac{1}{2}}\|_{\Omega}^2 \qquad (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) with $X^0 = X^{-1} = 0$, 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} \le 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}. \qquad (7.5)$$
<figure><img src="image_6.png" /><figcaption>FIGURE 5. Standing wave: reconstruction errors on the time derivative.</figcaption></figure>

*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| & \le \tau \|G^n\|_{\Omega} \|\dot{X}^{n+\frac{1}{2}}\|_{\Omega} + \tau \|G^n\|_{\Omega} \|\dot{X}^{n-\frac{1}{2}}\|_{\Omega} \\ & \le 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 \\ & \le 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}} \le 2\gamma\mu_0^{-1}\frac{\tau}{\rho}\|G^n\|_{\Omega}^2 + \left(1 + \frac{1}{4\gamma}\theta\right) 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} \le 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}}.$$
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Standing wave example with $\rho = 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>

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

$$E_X^{n+\frac{1}{2}} e^{-\rho t^n} \le c_1 \left(\frac{13}{24}, \theta\right) \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{align*}
\sum_{n \in \mathbb{N}} \tau E_X^{n+\frac{1}{2}} e^{-2\rho t^n} &\le c_1 \left(\frac{13}{24}, \theta\right) \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 \left(\frac{13}{24}, \theta\right) \mu_0^{-1} \frac{\tau^2}{\rho} \sum_{m \ge 0} \left( \sum_{n \ge m} e^{-\rho t^n} \right) \|G^m\|_{\Omega}^2 e^{-\rho t^m} \\
&\le c_1 \left(\frac{13}{24}, \theta\right) \theta (1 - e^{-\theta})^{-1} \mu_0^{-1} \frac{1}{\rho^2} \sum_{m \ge 0} \tau \|G^m\|_{\Omega}^2 e^{-2\rho t^m},
\end{align*}$$

since $\sum_{n \ge m} e^{-\rho t^n} = e^{-\rho t^m} (1 - e^{-\theta})^{-1}$ and $\theta = \rho\tau$. $\square$

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

**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: stability 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. \qquad (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. \qquad (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 \left( \mu_0 \| \dot{A}^{n+\frac{1}{2}} \|_{\Omega}^2 + \| \nabla A^{n+\frac{1}{2}} \|_{\Omega}^2 \right) 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, \qquad (7.7a)$$

$$\sum_{n \in \mathbb{N}} \tau \left( \mu_0 \| \dot{B}^{n+\frac{1}{2}} \|_{\Omega}^2 + \| \nabla B^{n+\frac{1}{2}} \|_{\Omega}^2 \right) 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, \qquad (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 $A^0 = A^{-1} = 0$ and 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

$$\begin{aligned} \sum_{n \in \mathbb{N}} \tau \left( \mu_0 \| \dot{A}^{n+\frac{1}{2}} \|_{\Omega}^2 + \| \nabla A^{n+\frac{1}{2}} \|_{\Omega}^2 \right) e^{-2\rho t^n} &\le \sum_{n \in \mathbb{N}} \tau E_A^{n+\frac{1}{2}} e^{-2\rho t^n} \\ &\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. \end{aligned}$$

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 $B^0 = B^{-1} = 0$ and 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}(\tilde{t})) dt \\ &= \int_{t^{n-\frac{1}{2}}}^{t^{n+1}} \phi^n(t) \left( \int_{\tilde{t}}^{t} \ddot{f}(s) ds \right) dt \\ &= \int_{t^{n-2}}^{t^{n+1}} \dddot{f}(s) \xi(s) ds, \quad \xi(s) := \int_{\tilde{s}}^{t^{n+1}} \phi_n(t) dt, \end{aligned} $$

with $\tilde{t} := 2t^{n-\frac{1}{2}} - t$ for all $t \in [t^{n-\frac{1}{2}}, t^{n+1}]$ and $\tilde{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$

### 7.3. 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 seminorm $|\phi|_{H^m(\mathcal{T}_h)}^2 := \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 seminorm $|\cdot|_Y$ equipping the space $Y$ composed of functions defined over $\Omega$, we use the shorthand notation $|v|_{L^2_{\rho}(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 \kappa h^k |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)}$ (a sharper $L^2$-estimate is not needed). 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 7.4** (Damped energy-norm a priori error estimate on $(u-u_{h\tau})$). Assume the CFL condition (2.15). Assume that $C_{\text{tim}}(u) := \sum_{r \in \{3:4\}} \kappa^{r-2} |u^{(r)}|_{L^2_{\rho}(V)}$ and $C_{\text{spa}}(u) := \sum_{r \in \{0:3\}} \kappa^r |u^{(r)}|_{L^2_{\rho}(H^{k+1}(\mathcal{T}_h))}$ are bounded. The following holds:

$$ \mathcal{E}_{\rho}^2(u - u_{h\tau}) \lesssim \tau^4 C_{\text{tim}}(u)^2 + h^{2k} C_{\text{spa}}(u)^2. \quad (7.8) $$
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 \left( \kappa |\dot{u}|_{L_p^2(H^{k+1}(T_h))} + |u|_{L_p^2(H^{k+1}(T_h))} \right).$$

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

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

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

(3) We observe that $X := u_h - U$ solves the leapfrog scheme with $X^0 = X^{-1} = 0$ and right-hand side

$$G^n := G_1^n + G_2^n := \pi_h^E((D_\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 $(D_\tau^2 u)^n := \frac{1}{\tau^2}(u(t^{n+1}) - 2u(t^n) + u(t^{n-1}))$. Combining the results of Lemmas 3.6 and 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.8 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)|_V^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}(T_h)}^2 dt + h^{2k} \ell_{\Omega}^2 \tau^2 \int_{J_n} |u^{(3)}(t)|_{H^{k+1}(T_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_p^2(V)} + \kappa^2 h^k \left( |u^{(2)}|_{L_p^2(H^{k+1}(T_h))} + \kappa |u^{(3)}|_{L_p^2(H^{k+1}(T_h))} \right).$$

This completes the proof.

**Corollary 7.5** (Damped energy-norm a priori error estimate on $(u - w_{h\tau})$). Under the assumptions of Theorem 7.4, and provided $C_{\text{tim}}(f) := \sum_{r \in \{2,3\}} \kappa^{r-2} \|f^{(r)}\|_{L_p^2(L^2)}$ is bounded, the following holds:

$$\mathcal{E}_{\rho}^{2}(u - w_{h\tau}) \lesssim \tau^{4}(C_{\text{tim}}(u)^{2} + \kappa^{2}C_{\text{tim}}(f)^{2}) + h^{2k}C_{\text{spa}}(u)^{2}. \quad (7.9)$$

Proof. Combine Theorem 7.4 with Lemma 3.9.

**Remark 7.6** (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$).

## FUNDING

This work was supported by the ANR JCJC project APOWA (research grant ANR-23-CE40-0019-01).
DATA AVAILABILITY STATEMENT

The research data associated with this article are included in the article.

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) 133–160.

[2] R.E. Bank and H. Yserentant, On the $H^1$-stability of the $L_2$-projection onto finite element spaces. *Numer. Math.* **126** (2014) 361–381.

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

[4] P. Bignardi and A. Moiola, A space-time continuous and coercive formulation for the wave equation. Preprint arXiv:2312.07268 (2023).

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

[6] M. Dauge, Elliptic Boundary Value Problems on Corner Domains. *Lecture Notes in Mathematics*. Vol. 1341. Springer-Verlag, Berlin (1988).

[7] A. Ern and J.-L. Guermond, Finite Elements I: Approximation and Interpolation. *Texts in Applied Mathematics*. Vol. 72. Springer Nature, Cham, Switzerland (2021).

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

[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) 2811–2834.

[10] T. Führer, R. González and M. Karkulik, Well-posedness of first-order acoustic wave equations and space-time finite element approximation. Preprint arXiv:2311.10536 (2023).

[11] F.D. Gaspoz, C.-J. Heine and K.G. Siebert, Optimal grading of the newest vertex bisection and $H^1$-stability of the $L_2$-projection. *IMA J. Numer. Anal.* **36** (2016) 1217–1241.

[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) 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) 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] M. Grote, O. Lakkis and C. Santos, A posteriori error estimates for the wave equation with mesh change in the leapfrog method. Preprint arXiv:2411.16933 (2024).

[17] 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).

[18] P. Joly, Variational methods for time-dependent wave propagation problems, in Topics in Computational Wave Propagation. *Lect. Notes Comput. Sci. Eng.* Vol. 31. Springer, Berlin (2003) 201–264.

[19] 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) 2374–2399.

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

[21] R.D. Skeel, Variable step size destabilizes the Störmer/leapfrog/Verlet method. *BIT* **33** (1993) 172–175.

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

[23] R. Verfürth, A Posteriori Error Estimation Techniques for Finite Element Methods. Oxford Science Publications (2013).

[24] M. Zank, *Inf-sup stable space-time methods for time-dependent partial differential equations*. Ph.D. thesis, Technische Universität Graz (2019).
<figure><img src="image_13.png" /></figure>

Please help to maintain this journal in open access!

This journal is currently published in open access under the **Subscribe to Open model** (S2O). We are thankful to our subscribers and supporters for making it possible to publish this journal in open access in the current year, free of charge for authors and readers.

Check with your library that it subscribes to the journal, or consider making a personal donation to the S2O programme by contacting subscribers@edpsciences.org.

More information, including a list of supporters and financial transparency reports, is available at https://edpsciences.org/en/subscribe-to-open-s2o.