# Asymptotically constant-free and polynomial-degree-robust a posteriori estimates for space discretizations of the wave equation

Théophile Chaumont-Frelet

► To cite this version:

Théophile Chaumont-Frelet. Asymptotically constant-free and polynomial-degree-robust a posteriori estimates for space discretizations of the wave equation. SIAM Journal on Scientific Computing, 2023, 45 (4), pp.A1591-A1620. 10.1137/22M1485619. hal-03632468

HAL Id: hal-03632468

https://inria.hal.science/hal-03632468

Submitted on 6 Apr 2022

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
ASYMPTOTICALLY CONSTANT-FREE AND
POLYNOMIAL-DEGREE-ROBUST A POSTERIORI ESTIMATES FOR
SPACE DISCRETIZATIONS OF THE WAVE EQUATION

T. CHAUMONT-FRELET†

ABSTRACT. We derive an equilibrated a posteriori error estimator for the space (semi) discretization of the scalar wave equation by finite elements. In the idealized setting where time discretization is ignored and the simulation time is large, we provide fully-guaranteed upper bounds that are asymptotically constant-free and show that the proposed estimator is efficient and polynomial-degree-robust, meaning that the efficiency constant does not deteriorate as the approximation order is increased. To the best of our knowledge, this work is the first to derive provably efficient error estimates for the wave equation. We also explain, without analysis, how the estimator is adapted to cover time discretization by an explicit time integration scheme. Numerical examples illustrate the theory and suggest that it is sharp.

KEY WORDS. a posteriori error estimate; equilibrated flux; finite element method; wave equation

1. INTRODUCTION

Given a domain $\Omega$, a partition of its boundary into disjoint sets $\Gamma_D$ and $\Gamma_A$ and space-
variable coefficients $\mu, \mathbf{A}$ and $\gamma$, the scalar wave equation consists in finding $u$ such that

$$
(1.1) \qquad \left\{
\begin{array}{l}
\mu\ddot{u} - \nabla \cdot (\mathbf{A}\nabla u) = \mu f \quad \text{in } \Omega, \\
\gamma\dot{u} + \mathbf{A}\nabla u \cdot \mathbf{n} = \gamma g \quad \text{on } \Gamma_A, \\
u = 0 \quad \text{on } \Gamma_D,
\end{array}
\right.
$$

$u|_{t=0} = u_0$ and $\dot{u}|_{t=0} = u_1$, where $f,g,u_0$ and $u_1$ are given right-hand sides and initial conditions. This problem models the propagation of (small) acoustics waves [34], and is a good mathematical simplification to later address more complicated physical waves arising for instance in elastodynamics and electromagnetism. As a result, the analysis of (1.1) and its numerical discretization has attracted much attention over the past decades [7, 10, 13, 16, 24, 25, 28, 29, 32, 33, 37].

Here, we are interested in finite element discretizations of (1.1), which have the advantage to easily handle complex geometries. For the sake of simplicity we focus on conforming Lagrange elements [11, 16], but discontinuous Galerkin [24, 32, 33] or hybridized [13, 29] methods could be considered as well. As compared to elliptic and parabolic problems (see, e.g, [4, 21] and [20, 36, 40], and the references therein), the a posteriori error analysis of (1.1) has not received much attention and it is the focus of the present work.

To the best of the author’s knowledge, the first references considering a posteriori error
estimation for the wave equation are [2, 3]. There, an asymptotically exact error estimator
is proposed for the space semi discretization. The analysis hinges on hierarchical polynomial

<footnote>†Inria Université Côte d'Azur, LJAD, CNRS</footnote>
basis, and strongly relies on tensor-product meshes and the smoothness of the solution. Hence, this methodology is not suited for general geometries which requires simplicial (or general hexahedral) meshes and where the solution may exhibit singularities.

In [10], an estimator is proposed for the full discretization of (1.1), where an implicit Euler scheme is used for time integration. Lower and upper bounds for the error are derived, but these bounds are not entirely satisfactory. Indeed, the norms employed to measure the error in the upper and lower bounds are different (compare [10, Corollary 5.3] and [10, Theorem 5.1]), indicating that the estimator may not be efficient. Besides, since this approach explicitly relies on the implicit Euler scheme, it is unclear whether it can be extended to space semi-discretization or to explicit time stepping. Finally, the behaviour of the constants in the estimates with respect to the problem parameters is unclear.

Another family of estimators, relying on the (finite element) solve of an elliptic problem at each time step is presented and analyzed in [25, 28, 37]. Upper bounds for the error measured in various norms, including the "$L^\infty(L^2)$" and "$L^2(H^1)$" norms, have been derived, but lower bounds are not available so far. The idea of employing an elliptic reconstruction is a technique inherited from parabolic problems [36], where implicit integration schemes are commonly used, due to restrictive CFL conditions. For wave propagation however, explicit time-stepping is possible, and almost systematically employed in large scale applications [7, 16, 24, 32, 33]. In this context, computing an elliptic reconstruction at each time step is in general very costly as compared to simply computing the solution.

Goal-oriented adaptive algorithms are also derived in [8, 9], where the focus is more on the algorithmic design and numerical tests than on theoretical analysis.

In this work, we discretize (1.1) with conforming Lagrange finite elements in space. We neglect the time-discretization, and study the resulting semi-discrete problem. In spirit, the design of our estimator is similar to the one in [10], as it is directly constructed using the "instantaneous residual" of the PDE and avoids the computation of an elliptic reconstruction. However, in contrast with [10], our analysis does not rely on the use of an implicit scheme and, although we focus on the semi-discrete case here, it could be extended to cover an explicit time scheme. Another key feature of our analysis is that we obtain constants that are much more explicit, in the upper bound and in the lower bound, than in [10]. Actually, although our methodology is general, we focus on equilibrated estimators so that the constant in our upper bound is fully computable and asymptotically equal to 1, and the constant in the lower bound is polynomial-degree-robust. Notice that similar results could be achieved with, e.g., a residual based estimator, at the price of introducing a (in general unknown) quasi-interpolation constant in the upper bound, and a $p$-dependent constant in the lower bound which is due to the use of "bubble functions" in the efficiency analysis.

Our setting is able to handle absorbing boundary conditions, which appears to be new in the context of a posteriori error estimation. We also allow for discontinuous material coefficients and (possibly) non-smooth solutions in space. However, we require that the right-hand sides $f$ and $g$ are smooth in time and for the sake of simplicity, we neglect the approximation of initial conditions (i.e. we assume that $u_0$ and $u_1$ are discrete functions). We also assume for shortness that $f$ and $g$ are piecewise polynomial in space, but we could classically alleviate this limitation by introducing the usual data oscillation terms. These assumptions are rigorously stated in Section 2.
The key idea of this work is to measure the error in a “damped energy norm”

$$ \mathcal{E}_{\rho}^{2} := \int_{0}^{+\infty} \left( \|\dot{u} - \dot{u}_{h}\|_{\mu, \Omega}^{2} + \frac{1}{\rho} \|\dot{u} - \dot{u}_{h}\|_{\gamma, \Gamma_{\Lambda}}^{2} + \|\nabla(u - u_{h})\|_{\mathbf{A}, \Omega}^{2} \right) e^{-2\rho t} dt, $$

where $\rho > 0$ is a user-defined damping parameter that can be selected as small as desired, and $u_h$ is the (semi-discrete) finite element approximation (notice that this norm is stronger than the $L^{\infty}(L^2)$ and $L^2(H^1)$ norms of [25] and [37]). We then propose an estimator $\eta(t)$ computed at each time step from the “instantaneous” residual

$$ v \to (\mu f(t), v)_{\Omega} + (\gamma g(t), v)_{\Gamma_{\Lambda}} - (\mu \ddot{u}_h(t), v) - (\gamma \dot{u}_h(t), v)_{\Gamma_{\Lambda}} - (\mathbf{A} \nabla u_h(t), \nabla v)_{\Omega}. $$

The computation of $\eta(t)$ is local in both time and space, and follows the equilibrated flux construction designed in [12, 14, 17, 21]. We emphasize that other standard techniques, such as residual-based constructions, can be employed to define $\eta(t)$, although the final constants in the estimates are less explicit in this case. Then, letting

$$ (1.2) \qquad \Lambda_{\rho}^{2} := \int_{0}^{+\infty} \eta^{2}(t) e^{-2\rho t} dt, $$

our key results read as follows. We have the upper bound

$$ (1.3) \qquad \mathcal{E}_{\rho}^{2} \le (1 + 4\gamma_{\rho,\omega,h}^{2})\Lambda_{\rho}^{2} + \left(\frac{\rho}{\omega}\right)^{2r} \mathrm{osc}_{\rho,r}^{2}, $$

where $\omega > 0$ and $r \in \mathbb{N}$ are arbitrary, $\mathrm{osc}_{\rho,r}$ is a fully computable “time-domain” data oscillation term that involves $\rho$ weighted norms of $f^{(r)}$ and $g^{(r)}$, and $\gamma_{\rho,\omega,h}$ is a constant such that

$$ (1.4) \qquad \gamma_{\rho,\omega,h} \le \sqrt{1 + \frac{\omega}{\rho}}, \quad \lim_{(h/p) \to 0} \gamma_{\rho,\omega,h} = 0. $$

The “cut-off frequency” $\omega$ can be freely selected, and in particular, we can choose it large enough so that the second term in the right-hand side of (1.3) is of “higher-order”, meaning that it converges to zero faster than the error (see Corollary 3.4 below). The estimate in (1.4) shows that all the constants in our upper bound are indeed explicitly controlled, while the limit justifies that our reliability estimate is asymptotically constant-free. On the other hand, we also show that

$$ (1.5) \qquad \Lambda_{\rho}^{2} \le C_{\mathrm{lb}}^{2} \kappa_{\mathbf{A}} \left\{ \left( \kappa_{\mathbf{A}} + \frac{\rho h_{\star}}{\vartheta_{\star}} + \left( \frac{\omega h_{\star}}{\vartheta_{\star}} \right)^{2} \right) \mathcal{E}_{\rho}^{2} + \left( \frac{\rho h_{\star}}{\vartheta_{\star}} \right)^{2} \left( \frac{\rho}{\omega} \right)^{2r} \mathrm{osc}_{\rho,r+1}^{2} \right\}, $$

where $h_{\star}/\vartheta_{\star} \simeq \max_K h_K/\vartheta_K$ with $\vartheta_K$ the wave speed in the element $K$, $\kappa_{\mathbf{A}}$ is the “contrast” in the coefficient $\mathbf{A}$, and $C_{\mathrm{lb}}$ is a generic constant solely depending on the shape-regularity of the mesh. All the constants appearing in (1.5) are independent of the polynomial degree $p$ of the finite element approximation, so that our estimator is efficient and polynomial-degree-robust. We refer the reader to Section 3 for more details.

In practice, one is often interested in understanding the behaviour of the solution $u$ of (1.1) on the time interval $(0, T)$ for some known $T > 0$. This can be accommodated in our setting by selecting $\rho \sim 1/T$. Another important comment is that in addition to ignoring time discretization, our setting also assumes that the simulation time is infinite, since the definition of the estimator $\Lambda_{\rho}$ in (1.2) includes an integral over $(0, +\infty)$. This is not an important limitation in practice as long as the simulation time is large enough since the integrand in (1.2) decays exponentially. In practical computations, the integral in (1.2) is
computed by accumulating the integrand over time steps, and the simulation can be stopped when stagnation is numerically observed. This is explained in detail in Section 7.1.

As far as the analysis is concerned, we take an approach that is fairly different from the aforementioned works on this topic. Indeed, our key idea is to perform a Laplace transform and to work in the frequency domain. We emphasize that this Laplace transform is only used as a theoretical device: it is never computed numerically. We then rely on recent results for the Helmholtz equation to control the low-frequency content of the error [14, 19] and introduce new arguments to treat the high-frequency content. We refer the reader to Section 4 for an in-depth discussion.

We provide several numerical examples where linear and quadratic Lagrange finite elements are coupled with an explicit leap-frog scheme [26]. These examples are in full agreement with our theoretical findings and suggest that they are sharp. Besides, they seem to indicate that time-discretization can be entirely neglected for linear elements (most likely, because they are coupled with a second-order integration scheme here), and only moderately affects quadratic elements when choosing a time-step close to the CFL stability limit.

The remaining of this work is organized as follows. Section 2 collects all the necessary notation to properly introduce our main results, which we state in Section 3. In Section 4, we present the strategy we employ in our analysis as well as some preliminary results in the frequency domain. Sections 5 and 6 respectively contain the reliability and efficiency proofs. We present a variety of numerical examples in Section 7 before giving some concluding remarks.

## 2. SETTINGS

This section describes notation and recalls preliminary results.

2.1. **Domain and coefficients.** $\Omega \subset \mathbb{R}^d$, $d = 2$ or $3$, is a polytopal Lischiptz domain. The boundary $\partial\Omega$ is split into two relatively open, polytopal and disjoint subsets $\Gamma_D$ and $\Gamma_A$ with Lipschitz boundaries such that $\partial\Omega = \bar{\Gamma}_D \cup \bar{\Gamma}_A$. The situations where $\Gamma_D = \emptyset$ or $\Gamma_A = \emptyset$ are allowed.

Our model wave equation features three coefficients, namely $\mu: \Omega \to \mathbb{R}$, $A: \Omega \to \mathbb{S}^d$ and $\gamma: \Gamma_A \to \mathbb{R}$. For the sake of simplicity, we assume that there exists a partition $\mathcal{P}$ of $\Omega$ into non-overlapping Lipschitz polytopal subdomains such that for all $P \in \mathcal{P}$, $\mu|_P$ and $A|_P$ are constant and, if $B := \partial P \cap \Gamma_A$ has positive surface measure, that $\gamma|_B$ is constant.

We assume that $\mu > 0$ in $\Omega$ and, for the sake of simplicity, that $\gamma > 0$ on $\Gamma_A$. We also introduce, for a.e. $x \in \Omega$

$$a_*(x) := \min_{\substack{d \in \mathbb{R}^d \\ |d|=1}} A(x) d \cdot d,$$

$$a^*(x) := \max_{\substack{d \in \mathbb{R}^d \\ |d|=1}} A(x) d \cdot d,$$

and require that $a_* > 0$ in $\Omega$.

2.2. **Mesh, hat functions and vertex patches.** The domain $\Omega$ is partitioned into a computational mesh $\mathcal{T}_h$ of (closed) simplicial elements $K$. We assume that the mesh is conforming, meaning that the intersection $K_+ \cap K_-$ of two distinct elements $K_\pm \in \mathcal{T}_h$ is either a single vertex, a full edge, or a full face of both elements. This requirement is entirely standard and does not prevent strong mesh grading (see, e.g., [11]).

For an element $K \in \mathcal{T}_h$, $h_K$ and $\rho_K$ denote the diameters of the smallest ball containing $K$ and of the largest ball contained in $K$. Then, $\beta_K := h_K/\rho_K$ is the shape regularity parameter
of $K$. If $\mathcal{T} \subset \mathcal{T}_h$ is a collection of elements $\beta_{\mathcal{T}} := \max_{K \in \mathcal{T}} \beta_K$ is the shape-regularity parameter of $\mathcal{T}$. The notations $h_{\max} := \max_{K \in \mathcal{T}_h} h_K$ and $\rho_{\min} := \min_{K \in \mathcal{T}_h}$ will also be useful.

$\mathcal{F}_h$ is the set of faces of $\mathcal{T}_h$. We assume that $\mathcal{T}_h$ conforms with the partition boundary, meaning that for each face $F \in \mathcal{F}_h$ such that $F \subset \partial\Omega$, either $F \subset \bar{\Gamma}_A$ or $F \subset \bar{\Gamma}_D$. The set $\mathcal{F}_h^A$ collects those faces contained in $\Gamma_A$.

We will also assume that the mesh fits the physical partition, meaning that for all $K \in \mathcal{T}_h$, there exists a $P \in \mathcal{P}$ such that $K \subset \bar{P}$. Equivalently, it means that $\mu|_K$ and $A|_K$ are constant for $K \in \mathcal{T}_h$, and that $\gamma|_F$ is constant for all $F \in \mathcal{F}_h^A$.

$\mathcal{V}_h$ is the set of vertices of $\mathcal{T}_h$. For $a \in \mathcal{V}_h$, $\psi_a$ denotes the associated hat function, that is, the only piecewise affine function on $\mathcal{T}_h$ such that $\psi_a(b) = \delta_{a,b}$ for all $b \in \mathcal{V}_h$. We set $\omega_a := \text{supp } \psi_a$. Then, the vertex patch $\mathcal{T}_a \subset \mathcal{T}_h$ collecting the elements $K \in \mathcal{T}_h$ sharing the vertex $a$ covers $\omega_a$.

**2.3. Key functional spaces.** If $U \subset \Omega$ is an open set, $L^2(U)$, $H^1(U)$ and $H^2(U)$ are the usual (real-valued) Lebesgue and Sobolev spaces [1]. $\mathcal{H}_{\Gamma_D}^1(\Omega)$ and $\mathcal{H}^2(\Omega)$ for the complex-valued counterparts of $H_{\Gamma_D}^1(\Omega)$ and $H^2(\Omega)$. $L_0^2(U)$ is the subset of $L^2(U)$ of functions with vanishing mean value, and $L^2(U) := [L^2(U)]^d$. The natural inner-products and norm of $L^2(U)$ and $L^2(U)$ are denoted by $(\cdot, \cdot)_U$ and $\|\cdot\|_U$. We also employ the equivalent norms $\|\cdot\|_{\mu,U}^2 := (\mu, \cdot)_U$ and $\|\cdot\|_{A,U}^2 := (A, \cdot)_U$ on $L^2(U)$ and $L^2(U)$. For a relatively open $V \subset \Gamma_A$, $L^2(V)$ is the Lebesgue space on $V$ equipped with the surface measure. $(\cdot, \cdot)_V$ and $\|\cdot\|_V$ are the norm and inner-product of $L^2(V)$, and we set $\|\cdot\|_{\gamma,V}^2 := (\gamma, \cdot)_V$. $H(\text{div}, U)$ is the subspace of $L^2(U)$ consisting of functions with weak $L^2(U)$ divergence [27]. If $w \subset \partial U$ is a relatively open subset, then $H_w^1(U)$ is the subset of $H^1(U)$ of functions with vanishing traces on $w$. Similarly, $H_w(\text{div}, U)$ is the subset $H(\text{div}, U)$ of functions with vanishing normal traces on $w$ (we refer the reader to [23] for a rigorous definition of normal traces on parts of the boundary).

When analyzing the flux construction, we will also need local spaces associated to vertex patches. If $a \in \partial\Omega$, we set $L_*^2(\omega_a) := L^2(\omega_a)$, $H_*^1(\omega_a) := H_{\partial\Omega}^1(\omega_a)$ and $H_0(\text{div}, \omega_a) := H_{\partial\omega_a \setminus \partial\Omega}(\text{div}, \omega_a)$. On the other hand, when $a \in \Omega$, we let $L_*^2(\omega_a) := L_0^2(\omega_a)$, $H_*^1(\omega_a) := H^1(\omega_a) \cap L_*^2(\omega_a)$ and $H_0(\text{div}, \omega_a) := H_{\partial\omega_a}(\text{div}, \omega_a)$.

**2.4. Local inequalities.** For all $a \in \mathcal{V}_h$, there exist constants $C_{P,a}$ and $C_{\text{tr},a}$ such that

$$ (2.1a) \qquad \|v\|_{\omega_a} \le C_{P,a} h_a \|\nabla v\|_{\omega_a}, \qquad \|v\|_{\partial\omega_a} \le C_{\text{tr},a} h_a^{1/2} \|\nabla v\|_{\omega_a} $$

for all $v \in H_*^1(\omega_a)$. Setting $C_{\text{cont},a} := 1 + C_{P,a} h_a \|\nabla \psi_a\|_{L^\infty(\omega_a)}$, we also have

$$ (2.1b) \qquad \|\nabla(\psi_a v)\|_{\omega_a} \le C_{\text{cont},a} \|\nabla v\|_{\omega_a}, $$

We then set $C_{\text{geo},a} := \max(C_{\text{cont},a}, C_{\text{tr},a})$, and $C_{\text{geo}} := 3(d+1) \max_{a \in \mathcal{V}_h} C_{\text{geo},a}$.

**2.5. Finite element spaces.** Let $K \in \mathcal{T}_h$ and $q \ge 0$. $\mathcal{P}_q(K)$ is the space of polynomial function on $K$ of degree at most $q$ and $\mathcal{P}_q(K) := [\mathcal{P}_q(K)]^d$. We also employ the notation $RT_q(K) := \mathcal{P}_q(K) + x\mathcal{P}_q(K)$ for the set of Raviart-Thomas polynomials [11]. Similarly, if $\mathcal{T} \subset \mathcal{T}_h$ and $q \ge 0$, $\mathcal{P}_q(\mathcal{T})$ and $RT_q(\mathcal{T})$ are the spaces of (discontinuous) functions whose restriction to each $K \in \mathcal{T}$ respectively belongs to $\mathcal{P}_q(K)$ and $RT_q(K)$.

In the remaining of this work, we fix a polynomial degree $p \ge 1$, and employ the notations $V_h := \mathcal{P}_p(\mathcal{T}_h) \cap H_{\Gamma_D}^1(\Omega)$ for the space Lagrange finite element of degree $p$, and $\mathcal{V}_h$ for its
complex-valued counterpart. Notice that there exists an interpolation operator $\mathcal{I}_h$ mapping $\mathcal{H}^2(\Omega) \cap \mathcal{H}_{\Gamma_D}^1(\Omega)$ into $\mathcal{V}_h$ such that

$$ (2.2) \qquad \|\nabla(v - \mathcal{I}_h v)\|_{\Omega} \le C_i h_{\max} \|\nabla^2 v\|_{\Omega} \quad \forall v \in \mathcal{H}^2(\Omega) \cap \mathcal{H}_{\Gamma_D}^1(\Omega), $$

where

$$ \|\nabla^2 v\|_{\Omega}^2 := \sum_{j\ell=1}^{d} \left\| \frac{\partial^2 v}{\partial \mathbf{x}_j \partial \mathbf{x}_{\ell}} \right\|_{\Omega}^2 $$

and $C_i$ only depends on the shape-regularity parameter of the mesh and is explicitly available [6, 35].

**2.6. Local wave speed and contrast.** The wave speed in the direction $\mathbf{d} \in \mathbb{R}^d$, $|\mathbf{d}| = 1$ at $\mathbf{x} \in \Omega$ is usually defined as $\vartheta(\mathbf{x}, \mathbf{d}) := \sqrt{\mathbf{A}(\mathbf{x})\mathbf{d} \cdot \mathbf{d}/\mu(\mathbf{x})}$. Besides, if the absorbing boundary condition on $\Gamma_A$ is designed to be exact for normally incident waves, the coefficient $\gamma$ is chosen as $\gamma(\mathbf{x}) = \vartheta(\mathbf{x}, \mathbf{n})/(\mathbf{A}(\mathbf{x})\mathbf{n} \cdot \mathbf{n})$. If $U \subset \Omega$ is an open set, this motivates the definition

$$ (2.3) \qquad \vartheta_U := \min \left( \sqrt{\frac{\inf_U a_*}{\sup_U \mu}}, \frac{\inf_U a_*}{\sup_{\partial U \cap \Gamma_A} \gamma} \right), $$

for the minimum wave speed in $U$, where we ignore the second term in the minimum if $\partial U \cap \Gamma_A = \emptyset$. We also employ the shorthand notation $\vartheta_a = \vartheta_{\omega_a}$ for $\mathbf{a} \in \mathcal{V}_h$. The quantity $\max_{\mathbf{a} \in \mathcal{V}_h} \nu h_{\mathbf{a}}/\vartheta_{\mathbf{a}}$, where $\nu > 0$ is a frequency will often appear in the analysis. To simplify it, we introduce $h_*/\vartheta_* := \max_{\mathbf{a} \in \mathcal{V}_h} h_{\mathbf{a}}/\vartheta_{\mathbf{a}}$. We also set $\vartheta_{\min} := \vartheta_{\Omega}$.

Our reliability estimate will also depend on the “contrast” in the coefficient $\mathbf{A}$ that we defined by

$$ (2.4) \qquad \kappa_{\mathbf{A},\mathbf{a}} := \frac{\sup_{\omega_{\mathbf{a}}} a^*}{\inf_{\omega_{\mathbf{a}}} a_*}, \qquad \kappa_{\mathbf{A}} := \max_{\mathbf{a} \in \mathcal{V}_h} \kappa_{\mathbf{A},\mathbf{a}}. $$

**2.7. Model problem.** The domain $\Omega$ and the partition $\{\Gamma_D, \Gamma_A\}$ of its boundary being defined, the problem is closed by specifying the right-hand sides $f$ and $g$, as well as the initial condition $u_0$ and $u_1$.

For the right-hand sides, we require that for all $t \in \mathbb{R}^+$, $f(t) \in \mathcal{P}_p(\mathcal{T}_h)$ and $g(t) \in \mathcal{P}_p(\mathcal{F}_h^A)$. The general case can be treated by adding “data oscillation” terms in the estimate, which we avoid here for the sake of shortness. We will further assume that $f \in C^\infty(\mathbb{R}^+, \mathcal{P}_p(\mathcal{T}_h))$ and $g \in C^\infty(\mathbb{R}^+, \mathcal{F}_h^A)$. Our analysis do require smoothness in time, but the $C^\infty$ assumption can be lowered to finite regularity. For the initial condition, we will for assume for the sake of simplicit that $u_0, u_1 \in \mathcal{V}_h$, so that the initial error at $t = 0$ vanish. The general case can be handled with additional terms in the error estimates.

The solution to the wave equation is the only function $u \in C^\infty(\mathbb{R}_+, H_{\Gamma_D}^1(\Omega))$ satisfying $u(0) = u_0$, $\dot{u}(0) = u_1$, and

$$ (2.5) \qquad (\mu\ddot{u}(t), v)_{\Omega} + (\gamma\dot{u}(t), v)_{\Gamma_A} + (\mathbf{A}\nabla u(t), \nabla v)_{\Omega} = (\mu f(t), v)_{\Omega} + (\gamma g(t), v)_{\Gamma_A} $$

for all $v \in H_{\Gamma_D}^1(\Omega)$ and $t \in \mathbb{R}_+$.
2.8. **Discretization.** Classically, the (semi) discrete solution is obtained by replacing the Sobolev space $H_{\Gamma_D}^1(\Omega)$ in (2.5) by its discrete counterpart $V_h$. Hence, we define the discrete solution as the unique function $u_h \in C^\infty(\mathbb{R}_+, V_h)$ such that $u_h(0) = u_0$, $\dot{u}_h(0) = u_1$, and

$$ (2.6) \quad (\mu\ddot{u}_h(t), v_h)_\Omega + (\gamma\dot{u}_h(t), v_h)_{\Gamma_A} + (\mathbf{A}\nabla u_h(t), \nabla v_h)_\Omega = (\mu f(t), v_h)_\Omega + (\gamma g(t), v_h)_{\Gamma_A} $$ 
for all $v_h \in V_h$ and $t \in \mathbb{R}_+.$

2.9. **Frequency-domain problem.** Similar to the results in [14, 19], our reliability estimate includes an “approximation factor” $\gamma_{\rho,\omega,h}$. In order to properly define it, we first need some notation from the frequency-domain.

For $s \in \mathbb{C}$ with $\rho := \text{Re } s > 0$, we introduce the sesquilinear form

$$ (2.7) \quad b_s(v, w) = s^2(\mu v, w)_\Omega + s(\gamma v, w)_{\Gamma_A} + (\mathbf{A}\nabla v, \nabla w)_\Omega \quad \forall v, w \in \mathcal{H}_{\Gamma_D}^1(\Omega), $$

which corresponds to the Laplace transform of the left-hand sides of (2.5) and (2.6). Since $\rho > 0$, $b_s$ is coercive, and if $\phi \in \mathcal{H}_{\Gamma_D}^1(\Omega)$, then there exists a unique $\mathcal{S}_s^*(\phi) \in \mathcal{H}_{\Gamma_D}^1(\Omega)$ such that

$$ (2.8) \quad b_s(w, \mathcal{S}_s^*(\phi)) = |s|^2(\mu w, \phi)_\Omega + \frac{|s|^2}{\rho}(\gamma w, \phi)_{\Gamma_A} \quad \forall v \in \mathcal{H}_{\Gamma_D}^1(\Omega). $$

With the notation $\{\cdot\}_{\mathcal{H}_s^1(\Omega)}^2 := |s|^2 \|\cdot\|_{\mu,\Omega}^2 + (|s|^2/\rho) \|\cdot\|_{\gamma,\Gamma_A}^2$, we can introduce the frequency-domain approximation factor

$$ (2.9) \quad \hat{\gamma}_{s,h} := \sup_{\substack{\phi \in \mathcal{H}_{\Gamma_D}^1(\Omega) \\ \{\phi\}_{\mathcal{H}_s^1(\Omega)} = 1}} \min_{v_h \in V_h} \|\nabla(\mathcal{S}_s^*(\phi) - v_h)\|_{\mathbf{A},\Omega}, $$

whose definition closely follows [14, 19].

2.10. **Approximation factor.** We are now in place to introduce the (time-domain) approximation factor. For any cut-off frequency $\omega > 0$ and damping parameter $\rho > 0$, it is defined by

$$ (2.10) \quad \gamma_{\rho,\omega,h} := \sup_{\substack{s=\rho+i\nu \\ |\nu|<\omega}} \hat{\gamma}_{s,h}. $$

The definition of $\gamma_{\rho,\omega,h}$ is a little bit intricate, but we can easily summarize its key properties. On the one hand, it increases when we increase $\omega$ and/or decrease $\rho$. On the other hand, it converges to zero as $(h/p) \to 0$. Proposition 2.1 states this more formally.

**Proposition 2.1 (Approximation factor).** We have

$$ (2.11) \quad \gamma_{\rho,\omega,h} \le \sqrt{1 + \frac{\omega}{\rho}}. $$

In addition, if $\Omega$ is convex, $\mathbf{A} = \mathbf{I}$ and $\Gamma_A = \emptyset$, then we have

$$ (2.12) \quad \gamma_{\rho,\omega,h} \le 2C_i \left( \frac{\rho h_{\max}}{\vartheta_{\min}} + \frac{\omega \omega h_{\max}}{\rho \vartheta_{\min}} \right). $$
2.11. **Local minimization problems.** Following [12, 14, 21], the construction of our estimator will rely on local divergence-constrained minimization problems. The following result is paramount to establish the efficiency of the estimator, its proof can be found in [12] for the 2D case, and in [22] for the 3D case.

**Proposition 2.2** (Discrete stable minimization). Let $q \in \mathbb{N}$, $\chi_q \in \mathbf{RT}_q(\mathcal{T}_a)$, $d_q \in \mathcal{P}_q(\mathcal{T}_a)$ and $b_q \in \mathcal{P}_q(\mathcal{F}_a)$. If $a \notin \bar{\Gamma}_D$, assume that

$$ (2.13) \qquad (d_q, 1)_{\omega_a} = (b_q, 1)_{\Gamma_A}. $$

We have

$$ (2.14) \qquad \min_{\substack{\tau_q \in \mathbf{RT}_q(\mathcal{T}_a) \cap \mathbf{H}_0(\text{div}, \omega_a) \\ \nabla \cdot \tau_q = d_q \text{ in } \omega_a \\ \tau_q \cdot n = b_q \text{ on } \Gamma_A}} \|\tau_q - \chi_p\|_{\omega_a} \le C_{\text{st},a} \min_{\substack{\tau \in \mathbf{H}_0(\text{div}, \omega_a) \\ \nabla \cdot \tau = d_q \text{ in } \omega_a \\ \tau \cdot n = b_q \text{ on } \Gamma_A}} \|\tau - \chi_p\|_{\omega_a} $$

where $C_{\text{st},a}$ only depends on the shape-regularity parameter of the patch and in particular, does not depend on $q$.

The constant $C_{\text{st}} := \max_{a \in \mathcal{V}_h} C_{\text{st},a}$ will be useful to state our efficiency result.

### 3. MAIN RESULTS

This section summarizes the key findings of this work.

**3.1. Equilibrated flux.** We first clarify what we mean by an equilibrated flux and propose a localized construction.

**Definition 3.1** (Equilibrated flux). An equilibrated flux is a function $\tau_h : \mathbb{R}_+ \to \mathbf{H}(\text{div}, \Omega)$ such that

$$ (3.1a) \qquad \nabla \cdot \tau_h(t) = \mu(f(t) - \ddot{u}_h(t)) \quad \text{in } \Omega $$

and

$$ (3.1b) \qquad \tau_h(t) \cdot n = \gamma(g(t) - \dot{u}_h(t)) \quad \text{on } \Gamma_A $$

for all $t \in \mathbb{R}_+$. For such $\tau_h$, we set

$$ (3.2a) \qquad \eta(t) := \|\mathbf{A}^{-1}\tau_h(t) + \nabla u_h(t)\|_{\mathbf{A}, \Omega}, $$

and

$$ (3.2b) \qquad \Lambda_{\rho}^{2} := \int_{0}^{+\infty} \eta(t)^{2} e^{-2\rho t} dt, $$

for $\rho > 0$.

Our first contribution is a construction of such a flux that is local in both space and time. This construction is standard, and follows the line of [12, 17, 21]. Specifically, for each $t \in \mathbb{R}_+$ and each vertex $a \in \mathcal{V}_h$, we set

$$ d^a := \psi_a \mu(f(t) - \ddot{u}_h(t)) - \mathbf{A} \nabla \psi_a \cdot \nabla u_h(t) \quad \text{in } \omega_a, \quad b^a := \psi_a \gamma(\dot{u}_h(t) - g(t)) \quad \text{on } \Gamma_A. $$

We then define local contributions by

$$ (3.3a) \qquad \sigma_h^a(t) := \arg \min_{\substack{\tau_h \in \mathbf{RT}_{p+1}(\mathcal{T}_a) \cap \mathbf{H}_0(\text{div}, \omega_a) \\ \nabla \cdot \tau_h = d^a(t) \text{ in } \omega_a \\ \tau_h \cdot n = b^a(t) \text{ on } \Gamma_A}} \|\mathbf{A}^{-1}\tau_h + \nabla u_h(t)\|_{\mathbf{A}, \omega_a} $$
that we assemble as

$$ (3.3b) \qquad \sigma_h(t) := \sum_{a \in V_h} \sigma_h^a(t). $$

Following the lines of [12, 14, 17, 21], we can easily show that the construction is indeed valid and provides an equilibrated flux as per (3.1). We skip the proof here for the sake of shortness.

**Proposition 3.2 (Localized flux construction).** The local mixed problems in (3.3a) are well-posed, and the construction (3.3b) provides an equilibrated flux satisfying (3.1).

**3.2. Reliability.** Our next set of results concerns the reliability of the proposed estimator. The error will be measured in the following norm

$$ \|u - u_h\|_{\rho}^2 := \int_{0}^{+\infty} \left\{ \|\dot{u} - \dot{u}_h\|_{\mu, \Omega}^2 + \frac{1}{\rho} \|\dot{u} - \dot{u}_h\|_{\gamma, \Gamma_{\Lambda}}^2 + \|\nabla(u - u_h)\|_{A, \Omega}^2 \right\} e^{-2\rho t} dt $$

where $\rho > 0$ is an arbitrary damping parameter. Our general result reads as follows.

**Theorem 3.3 (Reliability).** Assume that $\tau_h$ satisfies (3.1) and $\Lambda_\rho$ is defined by (3.2). Then, for all $\rho, \omega > 0$ and $r \in \mathbb{N}$, we have

$$ \|u - u_h\|_{\rho}^2 \le (1 + 4\gamma_{h,\omega,\rho}^2)\Lambda_{\rho}^2 + \left(\frac{\rho}{\omega}\right)^{2r}\mathrm{osc}_{\rho,r}^2 $$

where

$$ (3.4) \qquad \mathrm{osc}_{\rho,r}^2 = \frac{4}{\rho^{2r}} \left\{ \frac{1}{\rho^2} \|f^{(r)}e^{-\rho t}\|_{\mu,\Omega}^2 + \|g^{(r)}e^{-\rho t}\|_{\gamma,\Gamma_{\Lambda}}^2 \right\}. $$

Crucially, we can select $\omega$ and $r$ large enough so that the “oscillation” term converges to zero faster than the error itself. On the other hand, $\gamma_{h,\omega,\rho} \to 0$ as $(h/p) \to 0$, justifying that our upper bound is “asymptotically constant-free”. When the domain is convex and surrounded by a Dirichlet boundary condition, a simpler expression can be derived.

**Corollary 3.4 (Simplified error estimate).** Let $\tau_h$ satisfy (3.1) and define $\Lambda_\rho$ by (3.2a). Under the assumptions that $\Omega$ is convex, that $A = I$, and that $\Gamma_\Lambda = \emptyset$, we have

$$ \|u - u_h\|_{\rho}^2 \le \left( 1 + 16C_i^2 \left( \left( \frac{\rho h_{\max}}{\vartheta_{\min}} \right)^{1/2} + \frac{\rho h_{\max}}{\vartheta_{\min}} \right)^2 \right) \Lambda_{\rho}^2 + \left( \frac{\rho h_{\max}}{\vartheta_{\min}} \right)^{2(p+1)} \mathrm{osc}_{\rho,3(p+1)}^2 $$

for all $\rho > 0$.

**3.3. Efficiency.** Finally, we present our efficiency results. We start with the most general form.

**Theorem 3.5 (Efficiency).** Assume that $\tau_h$ has been constructed through the construction described in (3.3), and define $\Lambda_\rho$ with (3.2). The estimate

$$ \Lambda_{\rho}^{2} \le C_{\mathrm{lb}}^{2} \kappa_{A} \left\{ \left( \kappa_{A} + \frac{\rho h_{\star}}{\vartheta_{\star}} + \left( \frac{\omega h_{\star}}{\vartheta_{\star}} \right)^{2} \right) \|u - u_{h}\|_{\rho}^{2} + \left( \frac{\rho h_{\star}}{\vartheta_{\star}} \right)^{2} \left( \frac{\rho}{\omega} \right)^{2r} \mathrm{osc}_{\rho, r+1}^{2} \right\} $$

holds true for all $\rho, \omega > 0$ and $r \in \mathbb{N}$, with $C_{\mathrm{lb}} := C_{\mathrm{geo}}C_{\mathrm{st}}$.

Similar to the reliability results, $\omega$ and $r$ can be chosen so that the oscillation term converges to zero faster than the error. This is clearly highlighted in Corollary 3.6.
**Corollary 3.6** (Simplified efficiency estimate). Under the assumption of Theorem 3.5, for all $\rho > 0$, we have

$$ \Lambda_{\rho}^{2} \leq C_{\text{lb}}^{2} \kappa_{\mathbf{A}} \left\{ \left( \kappa_{\mathbf{A}} + 2 \frac{\rho h_{\star}}{\vartheta_{\star}} \right) \|u - u_{h}\|_{\rho}^{2} + \left( \frac{\rho h_{\star}}{\vartheta_{\star}} \right)^{2(p+1)} \mathrm{osc}_{\rho, 2(p+1)}^{2} \right\}. $$

4. PRELIMINARY RESULTS IN THE FREQUENCY DOMAIN

For the sake of simplicity, we will define $\xi_h \in C^\infty(\mathbb{R}_+, H_{\Gamma_D}^1(\Omega))$ by setting

$$ (4.1) \qquad \xi_h(t) := u(t) - u_h(t). $$

Classically, a central aspect of our analysis is to view the error $\xi_h$ as a particular solution to our PDE model, the associated right-hand side being the “residual”. To this end, we introduce for each $t \in \mathbb{R}_+$ the residual functional $\mathcal{R}(t) \in (H_{\Gamma_D}^1(\Omega))'$ by

$$ (4.2) \qquad \langle \mathcal{R}(t), v \rangle := (\mu(f(t) - \ddot{u}_h(t)), v)_{\Omega} + (\gamma(g(t) - \dot{u}_h(t)), v)_{\Gamma_A} - (\mathbf{A}\nabla u_h(t), \nabla v)_{\Omega} $$
for all $v \in H_{\Gamma_D}^1(\Omega)$, as well as its norm

$$ \|\mathcal{R}(t)\|_{-1,\Omega} := \sup_{\substack{v \in H_{\Gamma_D}^1(\Omega) \\ \|\nabla v\|_{\mathbf{A},\Omega}=1}} \langle \mathcal{R}, v \rangle. $$

Notice that then, $\xi_h(0) = \dot{\xi}_h(0) = 0$ due to our assumptions that $u_0, u_1 \in V_h$, and we have

$$ (4.3) \qquad (\mu\ddot{\xi}_h(t), v)_{\Omega} + (\gamma\dot{\xi}_h(t), v)_{\Gamma_A} + (\mathbf{A}\nabla\xi_h(t), \nabla v)_{\Omega} = \langle \mathcal{R}(t), v \rangle $$

for all $v \in H_{\Gamma_D}^1(\Omega)$ and $t \in \mathbb{R}_+$ as can be seen by subtracting (2.5) and (2.6).

When solving a steady problem with an inf-sup stable left-hand side, (4.3) can readily be employed to bound the discretization error by the residual norm (see [15, Equation (5.1)] for instance). This approach is also fruitful for parabolic problems (see, e.g., [20, Section 5]), where a suitable space-time inf-sup condition is available [39]. Unfortunately, to the best of the author's knowledge, such a framework is not available for the wave equation, making the link between the residual norm and the error harder to establish.

In this section, we develop the main idea of this work. It consists in establishing a relation between the residual norm and the error in the frequency-domain, and to treat distinctly the low-frequency and high-frequency content. For the low-frequency part, the analysis follows the lines of [14, 19], with the main difference that here, the frequency is complex-valued, with a positive imaginary part. We employ separate stability arguments to deal with the high-frequency content.

We establish all our main results in terms of the residual norm $\|\mathcal{R}(t)\|_{-1,\Omega}$, that we sometimes call the “idealized” estimator. Indeed, we believe this form is more general, since the residual norm can then be controlled by different types of estimators (including the equilibrated estimator we are focusing on here).

4.1. Laplace transform. The key tool we employ to connect the time and frequency domains is the Laplace transform. It is classically defined by

$$ \mathcal{L}\{v\}(s) = \int_{0}^{+\infty} v(t)e^{-st}dt, $$
whenever the integral is properly defined. For all $\rho > 0$, we have

$$ (4.4) \qquad \int_{0}^{+\infty} |v(t)|^2 e^{-2\rho t} dt = \int_{\rho-i\infty}^{\rho+i\infty} |\mathcal{L}\{v\}(s)|^2 ds $$

If $v$ is sufficiently regular and with $v(0) = 0$, we have

$$ (4.5) \qquad \mathcal{L}\{\dot{v}\}(s) = s\mathcal{L}\{v\}(s) $$

for all $s \in \mathbb{C}$ with $\mathrm{Re}\,s > 0$. Finally, we will use the following result to estimate high-frequency contents

$$ (4.6) \qquad \int_{\rho-i\infty}^{\rho+i\infty} |\mathcal{L}_{\rho}\{v\}(s)|^2 \chi_{|s|>\mu} ds \le \mu^{-2q} \int_{0}^{+\infty} |v^{(q)}(t)|^2 e^{-2\rho t} dt $$

In the remaining of this section we fix a complex number $s \in \mathbb{C}$ with $\rho := \mathrm{Re}\,s > 0$. For the sake of shortness, we often employ the notation $\hat{\phi} := \mathcal{L}\{\phi\}(s)$ for any function $\phi$ in the proofs.

**4.2. Frequency-domain problems.** Recalling (4.5), taking the Laplace transform of (2.5) and (2.6), we have

$$ (4.7a) \qquad b_s(\mathcal{L}\{u\}(s), v) = (\mu\mathcal{L}\{f\}(s), v)_{\Omega} + (\gamma\mathcal{L}\{g\}(s), v)_{\Gamma_A}, $$

and

$$ (4.7b) \qquad b_s(\mathcal{L}\{u_h\}(s), v_h) = (\mu\mathcal{L}\{f\}(s), v_h)_{\Omega} + (\gamma\mathcal{L}\{g\}(s), v_h)_{\Gamma_A}, $$

for all $v \in \mathcal{H}_{\Gamma_D}^1(\Omega)$ and $v_h \in \mathcal{V}_h$, where $b_s$ is the sesquilinear form defined at (2.7). Similarly, noticing that the definition of $\mathcal{R}$ naturally extends over $\mathcal{H}_{\Gamma_D}^1(\Omega)$, we can define $\mathcal{L}\{\mathcal{R}\}(s) \in (\mathcal{H}_{\Gamma_D}^1(\Omega))'$ by setting

$$ \langle \mathcal{L}\{\mathcal{R}\}(s), v \rangle = \int_{0}^{+\infty} \langle \mathcal{R}(t), v \rangle e^{-st} dt, $$

for all $v \in \mathcal{H}_{\Gamma_D}^1(\Omega)$, and we have

$$ (4.8) \qquad b_s(\mathcal{L}\{\xi_h\}(s), v) = \langle \mathcal{L}\{\mathcal{R}\}(s), v \rangle. $$

**4.3. Stability.** In the following, we equip $\mathcal{H}_{\Gamma_D}^1(\Omega)$ with norm

$$ \|v\|_{\mathcal{H}_s^1(\Omega)}^2 = |s|^2 \|v\|_{\mu, \Omega}^2 + \frac{|s|^2}{\rho} \|v\|_{\gamma, \Gamma_A}^2 + \|\nabla v\|_{\mathbf{A}, \Omega}^2 \quad \forall v \in H_{\Gamma_D}^1(\Omega). $$

Straightforward arguments then show that $b_s$ is coercive in the $\|\cdot\|_{\mathcal{H}_s^1(\Omega)}$ norm. Specifically,

$$ (4.9) \qquad \|v\|_{\mathcal{H}_s^1(\Omega)}^2 = \frac{1}{\rho} \mathrm{Re} \, b_s(v, sv) \le \frac{|s|}{\rho} |b_s(v, v)| \quad \forall v \in \mathcal{H}_{\Gamma_D}^1(\Omega). $$

As a direct consequence, we obtain a (coarse) upper bound for the frequency-domain error.

**Lemma 4.1** (Coarse frequency-domain upper bound). We have

$$ (4.10) \qquad \|\mathcal{L}\{\xi_h\}(s)\|_{\mathcal{H}_s^1(\Omega)}^2 \le 4 \left\{ \frac{1}{\rho^2} \|\mathcal{L}\{f\}(s)\|_{\mu, \Omega}^2 + \|\mathcal{L}\{g\}(s)\|_{\gamma, \Gamma_A}^2 \right\}. $$
*Proof.* On the one hand, recalling (4.9) and (4.7), we have

$$ \|\hat{u}_h\|_{\mathcal{H}_s^1(\Omega)}^2 \le \frac{|s|}{\mathrm{Re}\,s} |b_s(\hat{u}_h, \hat{u}_h)| = \frac{|s|}{\mathrm{Re}\,s} |(\mu \hat{f}, \hat{u}_h)_{\Omega} + (\gamma \hat{g}, \hat{u}_h)_{\Gamma}|. $$

On the other hand, we have

$$ \begin{aligned} |s| |(\mu \hat{f}, \hat{u}_h)_{\Omega} + (\gamma \hat{g}, \hat{u}_h)_{\Gamma}| &\le \|\hat{f}\|_{\mu, \Omega} |s| \|\hat{u}_h\|_{\mu, \Omega} + \rho \|\hat{g}\|_{\gamma, \Gamma_A} \frac{|s|}{\rho} \|\hat{u}_h\|_{\Gamma, \Gamma_A} \\ &\le \left( \|\hat{f}\|_{\mu, \Omega}^2 + \rho^2 \|\hat{g}\|_{\gamma, \Omega}^2 \right)^{1/2} \|\hat{u}_h\|_{\mathcal{H}_s^1(\Omega)}, \end{aligned} $$

and therefore

$$ (4.11) \qquad \|\hat{u}_h\|_{\mathcal{H}_s^1(\Omega)}^2 \le \frac{1}{\rho^2} \|\hat{f}\|_{\mu, \Omega}^2 + \|\hat{g}\|_{\gamma, \Gamma_A}^2. $$

Similar arguments show that (4.11) also holds for $u$, and (4.10) follows from the triangle inequality since $\hat{\xi}_h = \hat{u} - \hat{u}_h$. $\square$

**4.4. Approximation factor.** In order to refine the above error estimate, we will employ the approximation factor. To simplify the discussion below, we introduce the norm

$$ \{\phi\}_{\mathcal{H}_s^1(\Omega)}^2 := |s|^2 \|\phi\|_{\mu, \Omega}^2 + \frac{|s|^2}{\rho} \|\phi\|_{\gamma, \Gamma_A}^2 \quad \phi \in \mathcal{H}_{\Gamma_D}^1(\Omega) $$

on $\mathcal{H}_{\Gamma_D}^1(\Omega)$. Although it is not important in the forthcoming analysis, we note that this norm is not equivalent to the usual norm, and $\mathcal{H}_{\Gamma_D}^1(\Omega)$ is not an Hilbert space equipped with it. Denoting by $\Pi_h$ the orthogonal projection onto $\mathcal{V}_h$ for the $(\mathbf{A}\nabla\cdot, \cdot)_{\Omega}$ inner-product, we have

$$ (4.12) \qquad \|\nabla(\mathcal{S}_s^*(\phi) - \Pi_h \mathcal{S}_s^*(\phi))\|_{\mathbf{A}, \Omega} \le \hat{\gamma}_{s,h}\{\phi\}_{\mathcal{H}_s^1(\Omega)} \quad \forall \phi \in \mathcal{H}_s^1(\Omega). $$

Explicit upper bounds for the approximation factor are available, as we next demonstrate.

**Lemma 4.2 (Approximation factor).** *The estimate*

$$ (4.13) \qquad \hat{\gamma}_{s,h} \le \frac{|s|}{\rho} $$

holds true. In addition, if $\Omega$ is convex, $\mathbf{A} \equiv \mathbf{I}$ and $\Gamma_A = \emptyset$. Then, we have

$$ (4.14) \qquad \hat{\gamma}_{s,h} \le 2C_1 \frac{|s|}{\rho} \frac{|s|h_{\max}}{\vartheta_{\min}}. $$

*Proof.* Let $\phi \in \mathcal{H}_{\Gamma_D}^1(\Omega)$. Using (4.9), we have

$$ \begin{aligned} \|\mathcal{S}_s^*(\phi)\|_{\mathcal{H}_s^1(\Omega)}^2 &\le \frac{|s|}{\rho} |b_s(\mathcal{S}_s^*(\phi), \mathcal{S}_s^*(\phi))| \\ &= \frac{|s|}{\rho} \{\phi\}_{\mathcal{H}_s^1(\Omega)} \{\mathcal{S}_s^*(\phi)\}_{\mathcal{H}_s^1(\Omega)} \le \frac{|s|}{\rho} \{\phi\}_{\mathcal{H}_s^1(\Omega)} \|\mathcal{S}_s^*(\phi)\|_{\mathcal{H}_s^1(\Omega)}, \end{aligned} $$

so that $|s|\|\mathcal{S}_s^*(\phi)\|_{\mu, \Omega} \le \|\mathcal{S}_s^*(\phi)\|_{\mathcal{H}_s^1(\Omega)} \le (|s|/\rho)\{\phi\}_{\mathcal{H}_s^1(\Omega)}$, and (4.13) follows since

$$ \hat{\gamma}_{s,h} \le \sup_{\substack{\phi \in \mathcal{H}_{\Gamma_D}^1(\Omega) \\ \{\phi\}_{\mathcal{H}_s^1(\Omega)} = 1}} \|\nabla \mathcal{S}_s^*(\phi)\|_{\mathbf{A}, \Omega} \le \sup_{\substack{\phi \in \mathcal{H}_{\Gamma_D}^1(\Omega) \\ \{\phi\}_{\mathcal{H}_s^1(\Omega)} = 1}} \|\mathcal{S}_s^*(\phi)\|_{\mathcal{H}_s^1(\Omega)} \le \frac{|s|}{\rho}. $$
Assuming now that $\Omega$ is convex, $\mathbf{A} \equiv \mathbf{I}$ and $\Gamma_A = \emptyset$, we can apply [30, Theorem 3.2.1.2] and [31, Theorem 2.2.1], showing that

$$ \|\nabla^2 \mathcal{S}_s^*(\phi)\|_{\Omega} = \|\Delta \mathcal{S}_s^*(\phi)\|_{\Omega} = \|s^2 \mu\phi - s^2 \mu \mathcal{S}_s^*(\phi)\|_{\Omega}. $$

Since $\mathbf{A} \equiv \mathbf{I}$, we have $\mu \le \vartheta_{\min}^{-2}$, leading to

$$ \begin{aligned} \|\nabla^2 \mathcal{S}_s^*(\phi)\|_{\Omega} &\le \frac{1}{\vartheta_{\min}} (|s|^2 \|\phi\|_{\mu,\Omega} + |s|^2 \|\mathcal{S}_s^*(\phi)\|_{\mu,\Omega}) \\ &\le \frac{1}{\vartheta_{\min}} \left( |s| + \frac{|s|^2}{\rho} \right) \{\phi\}_{\mathcal{H}_s^1(\Omega)} \le \frac{2|s|^2}{\vartheta_{\min}\rho} \{\phi\}_{\mathcal{H}_s^1(\Omega)}, \end{aligned} $$

and (4.14) follows from the interpolation error estimate in (2.2) since

$$ \min_{v_h \in \mathcal{Y}_h} \|\nabla(\mathcal{S}_s^*(\phi) - v_h)\|_{\mathbf{A},\Omega} \le \|\nabla(\mathcal{S}_s^*(\phi) - \mathcal{I}_h \mathcal{S}^*(\phi))\|_{\Omega} \le C_h h_{\max} \|\nabla^2 \mathcal{S}_s^*(\phi)\|_{\Omega}. $$

□

**4.5. Reliability.** We are now ready to establish the key result of this section, which concerns the reliability of the “idealized estimator” in the frequency domain.

**Theorem 4.3 (Frequency-domain reliability).** The estimates

$$ (4.15) \qquad \|\mathcal{L}\{\xi_h\}(s)\|_{\mathcal{H}_s^1(\Omega)} \le \frac{|s|}{\rho} \|\mathcal{L}\{\mathcal{R}\}(s)\|_{-1,\Omega} $$

and

$$ (4.16) \qquad \|\mathcal{L}\{\xi_h\}(s)\|_{\mathcal{H}_s^1(\Omega)} \le (1 + 4\hat{\gamma}_{s,h}^2)^{1/2} \|\mathcal{L}\{\mathcal{R}\}(s)\|_{-1,\Omega} $$

hold true.

*Proof.* Selecting $v = \hat{\xi}_h$ in (4.8), we have

$$ (4.17) \qquad |b_s(\hat{\xi}_h, \hat{\xi}_h)| = |\langle \hat{\mathcal{R}}, \hat{\xi}_h \rangle| \le \|\hat{\mathcal{R}}\|_{-1,\Omega} \|\nabla \hat{\xi}_h\|_{\mathbf{A},\Omega}, $$

and (4.15) follows from (4.9).

The proof of (4.16) relies on an “Aubin-Nitsche trick”. Letting $\theta := \mathcal{S}_s^*(\hat{\xi}_h)$, picking $v = \hat{\xi}_h$ in (2.8) and using Galerkin orthogonality, we have

$$ \{\hat{\xi}_h\}_{\mathcal{H}_s^1(\Omega)}^2 = b_s(\hat{\xi}_h, \theta) = b_s(\hat{\xi}_h, \theta - \Pi_h \theta). $$

On the other hand, it follows from (4.8) and (4.12) that

$$ |b_s(\hat{\xi}_h, \theta - \Pi_h \theta)| \le \|\hat{\mathcal{R}}\|_{-1,\Omega} \|\nabla(\theta - \Pi_h \theta)\|_{\mathbf{A},\Omega} \le \hat{\gamma}_{s,h} \|\hat{\mathcal{R}}\|_{-1,\Omega} \{\hat{\xi}_h\}_{\mathcal{H}_s^1(\Omega)}, $$

so that

$$ (4.18) \qquad \{\hat{\xi}_h\}_{\mathcal{H}_s^1(\Omega)} \le \hat{\gamma}_{s,h} \|\hat{\mathcal{R}}\|_{-1,\Omega}. $$

Then, we have

$$ \|\nabla \xi_h\|_{\mathbf{A},\Omega}^2 = b(\hat{\xi}_h, \hat{\xi}_h) - s^2 \|\hat{\xi}\|_{\mu,\Omega}^2 - s \|\hat{\xi}\|_{\gamma,\Gamma}^2 \le |b(\hat{\xi}_h, \hat{\xi}_h)| + |s|^2 \|\hat{\xi}\|_{\mu,\Omega}^2 + |s| \|\hat{\xi}\|_{\gamma,\Gamma}^2, $$

and since $1 \le |s|/\rho$, we obtain that

$$ \|\nabla \xi_h\|_{\mathbf{A},\Omega}^2 \le |b(\hat{\xi}_h, \hat{\xi}_h)| + \{\hat{\xi}\}_{\mathcal{H}_s^1(\Omega)}^2. $$
It then follows from (4.17) and (4.18) that

$$
\begin{aligned}
\|\hat{\xi}_h\|_{\mathcal{H}_s^1(\Omega)}^2 &= \{\hat{\xi}_h\}_{\mathcal{H}_s^1(\Omega)}^2 + \|\nabla \xi_h\|_{\mathbf{A},\Omega}^2 \le |b(\hat{\xi}_h, \hat{\xi}_h)| + 2\{\hat{\xi}_h\}_{\mathcal{H}_s^1(\Omega)}^2 \\
&\le \|\hat{\mathcal{R}}\|_{-1,\Omega} \|\nabla \hat{\xi}\|_{\mathbf{A},\Omega} + 2\hat{\gamma}_{s,h}^2 \|\hat{\mathcal{R}}\|_{-1,\Omega}^2 \le \frac{1}{2} \|\hat{\xi}\|_{\mathcal{H}_s^1(\Omega)}^2 + \frac{1}{2} \|\hat{\mathcal{R}}\|_{-1,\Omega}^2 + 2\hat{\gamma}_{s,h}^2 \|\hat{\mathcal{R}}\|_{-1,\Omega}^2,
\end{aligned}
$$

from which (4.16) readily follows.

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

5. RELIABILITY

Here, we establish our reliability results.

5.1. **The damped energy norm.** The error will be measured in a damped energy norm.
For a damping parameter $\rho > 0$, we consider the norm

$$
\|v\|_{\rho}^{2} := \|\dot{v}e^{-\rho t}\|_{\mu, \Omega}^{2} + \frac{1}{\rho}\|\dot{v}e^{-\rho t}\|_{\gamma, \Gamma_{A}}^{2} + \|\nabla v e^{-\rho t}\|_{\mathbf{A}, \Omega}^{2}
$$

for all $v \in W^{1,\infty}(\mathbb{R}_+, H_{\Gamma_D}^1(\Omega))$, where we employed the notation

$$
\|v\|_{\dagger}^{2} = \int_{0}^{+\infty} \|v(t, \cdot)\|_{\dagger}^{2} dt,
$$

for any of the “space norms” $\|\cdot\|_{\dagger}$ introduced in Section 2.3. It is easily seen from (4.4) and (4.5) that actually

$$
(5.1) \qquad \|v\|_{\rho}^{2} = \int_{\rho-i\infty}^{\rho+\infty} \|\mathcal{L}\{v\}(s)\|_{\mathcal{H}_{s}^{1}(\Omega)}^{2} ds,
$$

5.2. **Approximation factor.** Recall the definition of the approximation factor $\gamma_{\rho,\omega,h}$ from (2.10). Then, Proposition 2.1 is a direct consequence of Lemma 4.2.

5.3. **Abstract reliability.** We now formulate our main reliability result in terms of the residual norm. As pointed out above, this abstract formulation permits to cover other types of estimators.

**Theorem 5.1 (Reliability).** For all $\omega > 0$, and $r \in \mathbb{N}$, we have

$$
(5.2) \qquad \|\xi_h\|_{\rho}^2 \le (1 + 4\gamma_{\rho,\omega,h}^2) \|\mathcal{R}e^{-\rho t}\|_{-1,\Omega}^2 + \left(\frac{\rho}{\omega}\right)^{2r} \mathrm{osc}_{\rho,r}^2.
$$

*Proof.* Let $\omega > 0$. Recalling (5.1), we have

$$
\|\xi_h\|_{\rho}^2 = \int_{\rho-i\infty}^{\rho+i\infty} \|\mathcal{L}\{\xi_h\}(s)\|_{\mathcal{H}_s^1(\Omega)}^2 ds = I_{\text{small}} + I_{\text{large}}
$$

with

$$
I_{\text{small}} := \int_{\rho-i\infty}^{\rho+i\infty} \|\mathcal{L}\{\xi_h\}(s)\|_{H_s^1(\Omega)}^2 \chi_{|s|<\omega} ds \quad I_{\text{large}} := \int_{\rho-i\infty}^{\rho+i\infty} \|\mathcal{L}\{\xi_h\}(s)\|_{H_s^1(\Omega)}^2 \chi_{|s|>\omega} ds.
$$

On the one hand, using (4.16), the definition of $\gamma_{\rho,\omega,h}$ in (2.10) and (4.4), we have

$$
\begin{aligned}
I_{\text{small}} &\le \int_{\rho-i\infty}^{\rho+i\infty} (1+4\hat{\gamma}_{s,h}^2) \|\mathcal{L}\{\mathcal{R}\}(s)\|_{-1,\Omega}^2 \chi_{|s|<\omega} ds \\
&\le (1+4\gamma_{\rho,\omega,h}^2) \int_{\rho-i\infty}^{\rho+i\infty} \|\mathcal{L}\{\mathcal{R}\}(s)\|_{-1,\Omega}^2 ds = (1+4\gamma_{\rho,\omega,h}^2) \|\mathcal{R}e^{-\rho t}\|_{-1,\Omega}^2.
\end{aligned}
$$
On the other hand, (4.10) and (4.6) show that

$$ I_{\text{small}} \le 4 \int_{\rho-i\infty}^{\rho+i\infty} \left\{ \frac{1}{\rho^2} \|\mathcal{L}\{f\}(s)\|_{\mu,\Omega}^2 + \|\mathcal{L}\{g\}(s)\|_{\gamma,\Gamma_A}^2 \right\} \\ \le \frac{4}{\omega^{2r}} \left\{ \frac{1}{\rho^2} \|f^{(r)}e^{-\rho t}\|_{\mu,\Omega}^2 + \|g^{(r)}e^{-\rho t}\|_{\gamma,\Gamma_A}^2 \right\} $$

and (5.2) follows recalling the definition of $\text{osc}_{\rho,r}$ in (3.4). □

We then give a simplified estimate assuming that holds when the $\Omega$ is convex, $A \equiv I$, and $\Gamma_A = \emptyset$. It follows by plugging (2.12) into (5.2), and then selecting $r = 3(p+1)$ and $\omega > 0$ such that $(\omega/\rho)^3 = \omega h_{\max}/\vartheta_{\min}$.

**Corollary 5.2 (Simplified error estimate).** Assume that $\Omega$ is convex, that $A \equiv I$ and that $\Gamma_A = \emptyset$. Then, for all $\rho > 0$, we have

$$ (5.3) \quad \|\xi_h\|_{\rho}^2 \le \left( 1 + 16C_i^2 \left( \left( \frac{\rho h_{\max}}{\vartheta_{\min}} \right)^{1/2} + \frac{\rho h_{\max}}{\vartheta_{\min}} \right)^2 \right) \|\mathcal{R}e^{-\rho t}\|_{-1,\Omega}^2 \\ \qquad + \left( \frac{\rho h_{\max}}{\vartheta_{\min}} \right)^{2(p+1)} \text{osc}_{\rho, 3p+1}^2. $$

5.4. **Application to the equilibrated estimator.** We conclude our reliability analysis by showing that the residual norm can be controlled by the equilibrated estimator using a Prager-Synge type argument [38]. We omit the proof as it is classical (see, e.g., [14, Proposition 4.1]).

**Proposition 5.3 (Control of the residual).** If $\tau_h$ satisfies (3.1) and $\eta$ and $\Lambda_\rho$ are defined by (3.2a), we have

$$ \|\mathcal{R}(t)\|_{-1,\Omega} \le \eta(t) \quad \forall t \in \mathbb{R}_+, \qquad \|\mathcal{R}e^{-\rho t}\|_{-1,\Omega} \le \Lambda_\rho \quad \forall \rho > 0. $$

At that point, Theorem 3.3 and Corollary 3.4 follow from Theorem 5.1, Corollary 5.2 and Proposition 5.3.

## 6. EFFICIENCY

In this section, we establish efficiency properties of the proposed estimator. The approach is in part similar to the analysis presented for the time-harmonic equations in [14, 19], but additional arguments are required to treat the second-time derivative (as opposed to a multiplication by $-\omega^2$).

**6.1. Localized norm.** An important ingredient of the forthcoming analysis will be the following localized residual norm

$$ (6.1) \quad \|\mathcal{R}(t)\|_{-1,\mathbf{a}} := \sup_{\substack{v \in H_\star^1(\omega_\mathbf{a}) \\ \|\nabla v\|_{A,\Omega}=1}} \langle \mathcal{R}(t), \psi_{\mathbf{a}}v \rangle $$

associated with each vertex $\mathbf{a} \in \mathcal{V}_h$. We start with a standard dual characterization showing that this norm is actually the minimum of a continuous version of the minimization problem defining the local flux contribution $\sigma_h^a$ in (3.3a). The proof is standard (see, e.g., [14, Lemma 4.3]) and omitted here for shortness.
**Lemma 6.1 (Dual characterization).** We have

$$ \|\mathcal{R}(t)\|_{-1, \mathbf{a}} = \min_{\substack{v \in H_{\Gamma_A}^{(\mathrm{div}, \omega_\mathbf{a})} \\ \nabla \cdot v = d^a \text{ in } \omega_\mathbf{a} \\ v \cdot n = b^a \text{ on } \Gamma_A}} \|\mathbf{A}^{-1}v + \psi_\mathbf{a} \nabla u_h(t)\|_{\mathbf{A}, \omega_\mathbf{a}} $$

for all vertices $\mathbf{a} \in \mathcal{V}_h$ and all $t \in \mathbb{R}_+$.

To establish a link between the localized and global norms of the residual, we observe that because the $\psi_\mathbf{a}$ form a partition of unity, we have

$$ \langle \mathcal{R}(t), v \rangle = \sum_{\mathbf{a} \in \mathcal{V}_h} \langle \mathcal{R}(t), \psi_\mathbf{a} v \rangle \le \left( \sum_{\mathbf{a} \in \mathcal{V}_h} \|\mathcal{R}(t)\|_{-1, \mathbf{a}}^2 \right)^{1/2} \left( \sum_{\mathbf{a} \in \mathcal{V}_h} \|\nabla v\|_{\mathbf{A}, \omega_\mathbf{a}}^2 \right)^{1/2} $$

Since each simplex $K$ has $d+1$ vertices, we obtain the following upper bound:

$$ (6.2) \qquad \|\mathcal{R}(t)\|_{-1, \Omega}^2 \le (d+1) \sum_{\mathbf{a} \in \mathcal{V}_h} \|\mathcal{R}(t)\|_{-1, \mathbf{a}}^2 \quad \forall t \in \mathbb{R}_+. $$

## 6.2. Abstract efficiency.
We now derive efficiency results for the idealized estimator. We first show that point-wise in time and patch-wise in space, the idealized estimator is a lower bound for a measure of the error that in addition to the desired norm, also includes the second time-derivative.

**Lemma 6.2 (Local efficiency).** For all vertices $\mathbf{a} \in \mathcal{V}_h$ and $t \in \mathbb{R}_+$, the estimate

$$ (6.3) \qquad \|\mathcal{R}(t)\|_{-1, \mathbf{a}} \le C_{\mathrm{geo}, \mathbf{a}} \left( \frac{h_{\mathbf{a}}}{\vartheta_{\mathbf{a}}} \|\ddot{\xi}_h(t)\|_{\mu, \omega_{\mathbf{a}}} + \left( \frac{\rho h_{\mathbf{a}}}{\vartheta_{\mathbf{a}}} \right)^{1/2} \frac{1}{\rho^{1/2}} \|\dot{\xi}_h(t)\|_{\gamma, \omega_{\mathbf{a}} \cap \Gamma_A} + \kappa_{\mathbf{A}}^{1/2} \|\nabla \xi_h(t)\|_{\mathbf{A}, \omega_{\mathbf{a}}} \right) $$

holds true.

*Proof.* Let $t \in \mathbb{R}_+$ and $v \in H_{\Gamma_D}^1(\Omega)$. Recalling (4.3), we have

$$ \begin{aligned} |\langle \mathcal{R}(t), \psi_{\mathbf{a}} v \rangle| &= |(\mu \ddot{\xi}_h(t), \psi_{\mathbf{a}} v)_{\Omega} + (\gamma \dot{\xi}_h(t), \psi_{\mathbf{a}} v)_{\Gamma_A} + (\mathbf{A} \nabla \xi_h(t), \nabla(\psi_{\mathbf{a}} v))_{\Omega}| \\ &\le \|\ddot{\xi}_h(t)\|_{\mu, \Omega} \|v\|_{\mu, \Omega} + \|\dot{\xi}_h(t)\|_{\gamma, \Gamma_A} \|v\|_{\gamma, \Gamma_A} + \|\nabla \xi_h(t)\|_{\mathbf{A}, \Omega} \|\nabla(\psi_{\mathbf{a}} v)\|_{\mathbf{A}, \Omega}. \end{aligned} $$

Then, (6.3) follows by applying (2.1), recalling the definitions of $\vartheta_\mathbf{a}$ and $\kappa_\mathbf{A}$ in (2.3) and (2.4) as well as the definition of the localized residual norm in (6.1). $\square$

Next, we show that globally in time and space, the second time-derivative can be removed, at the price of introducing a data oscillation term.

**Lemma 6.3 (Second time-derivative).** For all $\omega > 0$ and $r \in \mathbb{N}$, we have

$$ (6.4) \qquad \|\ddot{\xi}_h e^{-\rho t}\|_{\mu, \Omega}^2 \le \omega^2 \|\dot{\xi}_h e^{-\rho t}\|_{\mu, \Omega}^2 + \rho^2 \left(\frac{\rho}{\omega}\right)^{2r} \mathrm{osc}_{\rho, r+1}^2. $$

*Proof.* Standard identity (4.4) of the Laplace transform shows that

$$ \|\ddot{\xi}_h e^{-\rho t}\|_{\mu, \Omega}^2 = \int_{\rho-i\infty}^{\rho+i\infty} \|\mathcal{L}\{\ddot{\xi}_h\}(s)\|_{\mu, \Omega}^2 ds. $$
We split the integral in the right-hand side in two. On the one hand, using (4.5) and (4.4) again, we have

$$ \int_{\rho-i\infty}^{\rho+i\infty} \|\mathcal{L}\{\ddot{\xi}_h\}(s)\|_{\mu,\Omega}^2 \chi_{|s|<\omega} ds = \int_{\rho-i\infty}^{\rho+i\infty} |s|^2 \|\mathcal{L}\{\dot{\xi}_h\}(s)\|_{\mu,\Omega}^2 \chi_{|s|<\omega} ds. \\ \le \omega^2 \int_{\rho-i\infty}^{\rho+i\infty} |\mathcal{L}\{\dot{\xi}_h\}(s)|_{\mu,\Omega}^2 ds = \omega^2 \|\dot{\xi}_h e^{-\rho t}\|_{\mu,\Omega}^2. $$

On the other hand, using (4.5) and (4.10)

$$ \int_{\rho-i\infty}^{\rho+i\infty} \|\mathcal{L}\{\ddot{\xi}_h\}(s)\|_{\mu,\Omega}^2 \chi_{|s|>\omega} ds \le \int_{\rho-i\infty}^{\rho+i\infty} |s|^4 \|\mathcal{L}\{\dot{\xi}_h\}(s)\|_{\mu,\Omega}^2 \chi_{|s|>\omega} ds \\ \le \frac{4}{\rho^2} \int_{\rho-i\infty}^{\rho+i\infty} |s|^2 |\mathcal{L}\{f\}(s)|_{\mu,\Omega}^2 \chi_{|s|>\omega} ds + 4 \int_{\rho-i\infty}^{\rho+i\infty} |s|^2 |\mathcal{L}\{g\}(s)|_{\gamma,\Gamma_A}^2 \chi_{|s|>\omega} ds. $$

Recalling the definition of $\text{osc}_{\rho,r}$ in (3.4), we conclude the proof with (4.5) and (4.6), since

$$ \int_{\rho-i\infty}^{\rho+i\infty} |s|^2 |\mathcal{L}\{f\}(s)|_{\mu,\Omega}^2 \chi_{|s|>\omega} ds = \int_{\rho-i\infty}^{\rho+i\infty} |\mathcal{L}\{\dot{f}\}(s)|_{\mu,\Omega}^2 \chi_{|s|>\omega} ds \le \omega^{-2r} \|f^{(1+r)} e^{-\rho t}\|_{\mu,\Omega}^2, $$

and a similar estimate holds for $g$, for all $r \in \mathbb{N}$. $\square$

**Remark 6.4** (Localized treatment of the second time-derivative). It is possible to obtain a partially localized version of (6.4), with the term $\|\dot{\xi}_h e^{-\rho t}\|_{\mu,\Omega}$ and $\|\ddot{\xi}_h e^{-\rho t}\|_{\mu,\Omega}$ respectively replaced by $\|\ddot{\xi}_h e^{-\rho t}\|_{\mu,\omega_a}$ and $\|\dot{\xi}_h e^{-\rho t}\|_{\mu,\omega_a}$. However, the author does not see a way to localize the data oscillation term.

Combining Lemmas 6.2 and 6.3, we arrive at our key efficiency result.

**Theorem 6.5** (Global efficiency). The estimate

$$ (6.5) \quad \|\mathcal{R}e^{-\rho t}\|_{-1,\Omega}^2 \le C_{\text{geo}}^2 \left\{ \left( \kappa_A + \frac{\rho h_*}{\vartheta_*} + \left( \frac{\omega h_*}{\vartheta_*} \right)^2 \right) \|\xi\|_{\rho}^2 + \left( \frac{\rho h_*}{\vartheta_*} \right)^2 \left( \frac{\rho}{\omega} \right)^{2r} \text{osc}_{\rho,r+1}^2 \right\} $$

holds true for all $\omega > 0$ and $r \in \mathbb{N}$.

*Proof.* Recall the estimates in (2.1). By summing the local estimates in (6.3) over $a \in \mathcal{V}_h$, and using (6.2), we obtain

$$ \|\mathcal{R}(t)\|_{-1,\Omega}^2 \le (d+1) \max_{a \in \mathcal{V}_h} C_{\text{geo},a}^2 \\ \times \sum_{a \in \mathcal{V}_h} \left( \frac{h_a}{\vartheta_a} \|\ddot{\xi}_h(t)\|_{\mu,\omega_a} + \left( \frac{\rho h_a}{\vartheta_a} \right)^{1/2} \frac{1}{\rho^{1/2}} \|\dot{\xi}_h(t)\|_{\gamma,\omega_a \cap \Gamma_A} + \kappa_{A,a} \|\nabla \xi_h(t)\|_{A,\omega_a} \right)^2 \\ \le C_{\text{geo}}^2 \left( \left( \frac{h_*}{\vartheta_*} \right)^2 \|\ddot{\xi}_h(t)\|_{\mu,\Omega}^2 + \frac{\rho h_*}{\vartheta_*} \frac{1}{\rho} \|\dot{\xi}_h(t)\|_{\gamma,\Gamma_A}^2 + \kappa_{A}^2 \|\nabla \xi_h(t)\|_{A,\Omega}^2 \right), $$

and therefore

$$ \|\mathcal{R}e^{-\rho t}\|_{-1,\Omega}^2 \le C_{\text{geo}}^2 \left( \left( \frac{h_*}{\vartheta_*} \right)^2 \|\ddot{\xi}_h e^{-\rho t}\|_{\mu,\Omega}^2 + \frac{\rho h_*}{\vartheta_*} \frac{1}{\rho} \|\dot{\xi}_h e^{-\rho t}\|_{\gamma,\Gamma_A}^2 + \kappa_{A}^2 \|\nabla \xi_h e^{-\rho t}\|_{A,\Omega}^2 \right). $$
We then apply (6.4), showing that

$$ \left(\frac{h_{\star}}{\vartheta_{\star}}\right)^{2} \|\ddot{\xi}_{h} e^{-\rho t}\|_{\mu, \Omega}^{2} \leq \left(\frac{\omega h_{\star}}{\vartheta_{\star}}\right)^{2} \|\dot{\xi}_{h} e^{-\rho t}\|_{\mu, \Omega}^{2}+\left(\frac{\rho h_{\star}}{\vartheta_{\star}}\right)^{2}\left(\frac{\rho}{\omega}\right)^{2 r} \operatorname{osc}_{\rho, r+1}^{2}, $$

and leading to (6.5).

□

For the reader's convenience, we also provide a simplified version of our efficiency estimate. It is obtained from (6.5) by selecting $\omega > 0$ such that $(\omega h_{\star}/\vartheta_{\star})^2 = \rho h_{\star}/\vartheta_{\star}$ and $r = 2p$.

**Corollary 6.6 (Simplified efficiency estimate).** We have

$$ (6.6) \qquad \|\mathcal{R}e^{-\rho t}\|_{-1, \Omega}^2 \le C_{\text{geo}}^2 \left\{ \left( \kappa_A + 2 \frac{\rho h_{\star}}{\vartheta_{\star}} \right) \|\xi\|_{\rho}^2 + \left( \frac{\rho}{\omega} \right)^{2(p+1)} \operatorname{osc}_{\rho, 2p+1}^2 \right\}. $$

6.3. **Application to the equilibrated estimator.** Theorem 3.3 and Corollary 3.4 simply follow from Proposition 2.2, Theorem 6.5 and Corollary 5.2.

## 7. NUMERICAL EXAMPLES

This section presents a set of 2D numerical examples. We use elements of degree $p = 1$ or 2, coupled with an explicit leap-frog scheme for time integration [26]. The meshes are generated with the `mmg` software [18].

**7.1. Time discretization.** Considering a time step $\delta t > 0$ we set $t_n = n\delta t$, and build a sequence of approximation $(u_h^n)_{n \in \mathbb{N}}$, where $u_h^n$ is meant to approximate $u(t_n)$. We employ a leap-frog scheme to mimick the time-derivative [26]. Specifically, we introduce

$$ D_{\delta t}^2 u_h^n = \frac{u_h^{n+1} - 2u_h^n + u_h^{n-1}}{\delta t^2}, \qquad D_{\delta t} u_h^n = \frac{u_h^{n+1} - u_h^{n-1}}{2\delta t}. $$

The sequence $(u_h^n)_{n \in \mathbb{N}}$ is then defined by setting $u_h^0 = 0$, $u_h^1 = 0$, and then iteratively by requiring that

$$ (7.1) \qquad (\mu D_{\delta t}^2 u_h^n, v_h)_{\Omega} + (\gamma D_{\delta t} u_h^n, v_h)_{\Gamma_A} + (A \nabla u_h^n, \nabla v_h)_{\Omega} = (\mu f(t_n), v_h)_{\Omega} + (\gamma g(t_n), v_h)_{\Gamma_A} $$

for all $v_h \in V_h$. Classically, the choice of $\delta t$ is restricted by a so-called “CFL condition”, hence we set

$$ (7.2) \qquad \delta t = \alpha \min_{K \in \mathcal{T}_h} \frac{\rho_K}{\vartheta_K} $$

where $\alpha > 0$ is selected sufficiently small to obtain a stable discretization (recall that $\rho_K$ is inscribed diameter of the element $K$ and that $\vartheta_K$ is the wave speed inside $K$). For the meshes we employ, we empirically find that the values $\alpha = 1.5$ when $p = 1$ and $\alpha = 0.6$ when $p = 2$ are close to the CFL limit.

This leads to an explicit time-integration scheme where only the matrix representation of the form

$$ V_h \ni w_h, v_h \to \mathbb{R} \to \frac{1}{\delta t^2} (\mu w_h, v_h)_{\Omega} + \frac{1}{2\delta t} (\gamma w_h, v_h)_{\Gamma_A} \in \mathbb{R} $$

needs to be inverted. For the sake of simplicity, we employ the `mumps` package [5] to obtain the Choleski factorization of the above matrix in our computations. In practice however, it is possible to employ mass-lumping to avoid factorizing the matrix [16].
In practice, we can only compute a finite number of iterates, meaning that the computations stop after $N \in \mathbb{N}$ steps, with associated time $T := t_N$. Our estimates require an infinite simulation time, but this can be easily circumvented by monitoring

$$ \int_{0}^{T} \eta(t)^2 e^{-2\rho t} dt $$

as $T$ increases and stopping the simulation when it stagnates.

For each time step $n \in \{0, \dots, N\}$, and each vertex $\mathbf{a} \in \mathcal{V}_h$, we employ the time-discrete version of (3.3a) to define $\boldsymbol{\sigma}_h^{a,n}$

$$ (7.3) \quad \boldsymbol{\sigma}_h^{a,n} := \arg \min_{\substack{\tau_h \in \mathcal{R}T_{p+1}(\mathcal{T}_a) \cap \mathcal{H}_0(\mathrm{div}, \omega_a) \\ \nabla \cdot \boldsymbol{\tau}_h = \mu \psi_a (f(t_n) - D_{\delta t}^2 u_h^n) - \mathcal{A} \nabla \psi_a \cdot \nabla u_h^n \text{ in } \omega_a \\ \boldsymbol{\tau}_h \cdot \mathbf{n} = \gamma \psi_a (g(t_n) - D_{\delta t} u_h^n) \text{ on } \Gamma_A}} \| \mathcal{A}^{-1} \boldsymbol{\tau}_h + \nabla u_h^n \|_{\mathcal{A}, \omega_a} $$

and we set

$$ \eta_K^n := \| \mathcal{A}^{-1} \boldsymbol{\sigma}_h^n + \nabla u_h^n \|_{\mathcal{A}, K}, \quad \eta^n := \left( \sum_{K \in \mathcal{T}_h} (\eta_K^n)^2 \right)^{1/2}, $$

as well as

$$ (7.4) \quad \Lambda_{\rho}^{2} := \frac{1}{2\delta t^2} \sum_{n=0}^{N} ((\eta^{n+1})^2 e^{-2\rho t_{n+1}} + (\eta^n)^2 e^{-2\rho t_n}). $$

Notice that the compatibility condition in (7.3) is satisfied due to (7.1).

**7.2. Standing wave.** We consider the unit square $\Omega = (0,1)^2$ surrounded by a Dirichlet boundary condition with $\Gamma_D = \partial\Omega$. We consider the solution

$$ u(t, \boldsymbol{x}) = \chi(t) \sin(\sqrt{2\pi t}) \sin(\pi \boldsymbol{x}_1) \sin(\pi \boldsymbol{x}_2), $$

where $\chi$ is a cutoff function defined as follows. For $t \le 0$, $\chi(t) = 0$ and $\chi(t) = 1$ when $t \ge 1$. In the interval $[0, 1]$ $\chi$ is defined as the unique $\mathcal{P}_5$ polynomial that enables $C^2$ junctions. Observe that for $t \ge 1$, $u$ solves the wave equation with a vanishing right-hand side $f$. We perform the simulation on the interval $(0, T)$ with $T := 10$.

The goal of this example is to emphasize the role of the damping parameter $\rho > 0$. First, Figure 7.1 presents the evolution of the instantaneous error (recall that $\Gamma_A = \emptyset$ here)

$$ \mathcal{E}^2(t) := \|\dot{\xi}(t)\|_{\mu, \Omega}^2 + \|\nabla \xi(t)\|_{\mathcal{A}, \Omega}^2, $$

and the “instantaneous” estimator $\eta(t)$, as well as the cumulated error

$$ C_{\rho}^{2}(t) := \int_{0}^{t} \mathcal{E}^{2}(\tau) e^{-2\rho\tau} d\tau, $$

and the cumulated estimator $\Lambda_{\rho}(t)$ obtained by summing (7.4) up to the current time step. On Figure 7.1a, we can see that $\mathcal{E}(t)$ globally increases with $t$, whereas $\eta(t)$ remains globally constant. The reason is that, similar to what happens in the frequency domain [14], the estimator is insensitive to the dispersion error that accumulates over time. The effect is counter balanced by the parameter $\rho$ in the damped energy norm. As shown on Figure 7.1b, 7.1c and 7.1d, for a fixed mesh, the estimator tends to underestimate the error as $\rho$ is decreased.
<figure><img src="image_3.png" /><figcaption>FIGURE 7.1. Error evolution in the standing wave example for $h_{max} = 0.05$</figcaption></figure>

Figure 7.2 presents the behaviour of the estimator $\Lambda_{\rho}$ and the error as the mesh is refined. As claimed, our estimate is asymptotically constant-free: the two curves becomes indistinguishable as $h \to 0$. The effect is analyzed in more depth on Figure 7.3, where the effictivity index $\text{eff}_{\rho} := \Lambda_{\rho}/\|u - u_h\|_{\rho}$ is plotted against the mesh size for different values of $\rho$, $\alpha$ and $p$. We can see on Figures 7.3a and 7.3b that for $p = 1$ and $p = 2$ (with an over-refined time step), we observe exactly the behaviour predicted by our analysis: the estimate is asymptotically constant-free, and the asymptotic regime is achieved faster when $\rho$ and $p$ are larger. On Figure 7.3c, we display the same results for $p = 2$ and a “natural” time step close to the CFL stability limit. We can see in this case that the time-discretization error have a small contribution that is not capture by the estimator.

**7.3. Reflections by Dirichlet boundaries.** In this example, the domain is the unit square $\Omega := (0, 1)^2$. On the “bottom-left” boundary of $\Omega$, we impose a Dirichlet boundary condition, whereas on the “top-right” part, we place an absorbing boundary condition.
<figure><img src="image_4.png" /><figcaption>FIGURE 7.2. Convergence in the standing wave example for $\rho = 0.5$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 7.3. Effectivity index in the standing wave example</figcaption></figure>
Let us first consider the travelling wave

$$v_{\sigma,d}(t, \mathbf{x}) = p_{\sigma}((t - t_0) - \mathbf{d} \cdot \mathbf{x}),$$

where $t_0 = 4$, $\mathbf{d} = (\cos \theta, \sin \theta)$ with $\theta = 11\pi/8$, and $p_{\sigma}$ is the profile

$$p_{\sigma}(\tau) = \tau e^{-(\tau/\sigma)^2} \quad \tau \in \mathbb{R}.$$

We then construct our analytical solution by the principle of images, that is by “mirroring” $v_{\sigma,d}$ along the $\{\mathbf{x}_1 = 0\}$ and $\{\mathbf{x}_2 = 0\}$ hyperplanes. Specifically, we set

$$u(t, \mathbf{x}) = v_{\sigma,d}(t, \mathbf{x}_1, \mathbf{x}_2) - v_{\sigma,\mathbf{x}}(t, -\mathbf{x}_1, \mathbf{x}_2) - v_{\sigma,\mathbf{x}}(t, \mathbf{x}_1, -\mathbf{x}_2) + v_{\sigma,d}(t, -\mathbf{x}_1, -\mathbf{x}_2).$$

It is easily seen that $u$ satisfies the volume equation with $f = 0$ in $\Omega$ and the Dirichlet condition on $\Gamma_D$. We then set $g = \dot{u} + \nabla u \cdot \mathbf{n}$ on $\Gamma_A$. The initial condition $u(0, \mathbf{x})$ and $\dot{u}(0, \mathbf{x})$ are not exactly zero, but are sufficiently small that setting $u^0 = u^1 = 0$ in our computations generates a level of error less than numerical discretization. The simulation time is $T := 10$, and we set $\rho := 1/T$.

Here, our goal is to highlight the sensitivity of the estimator to the oscillations present in the right-hand side, that are here described by the parameter $\sigma$. Figure 7.4 represents the evaluation of the instantaneous and cumulated errors

$$\mathcal{E}_{\rho}^{2}(t) := \|\dot{\xi}_{h}(t)\|_{\mu,\Omega}^{2} + \frac{1}{\rho}\|\dot{u}(t)\|_{\gamma,\Gamma_{A}}^{2} + \|\nabla\xi_{h}(t)\|_{\mathbf{A},\Omega}^{2}, \quad \mathcal{C}_{\rho}^{2}(t) := \int_{0}^{t} \mathcal{E}_{\rho}^{2}(\tau)e^{-2\rho\tau}d\tau,$$

and estimators, whereas Figure 7.5 presents the behaviour of the estimator for $\sigma = 0.5$.

Figure 7.6 displays the effictivity index for different mesh sizes and values of $\sigma$ and $p, \alpha$. As predicted by our analysis, if the time-step is sufficiently fine (Figures 7.6a and 7.6b), the estimate becomes asymptotically constant-free as the mesh is refined for any value of $\sigma$, and the asymptotic regime is achieved faster if $p$ is large or if $\sigma$ is small (i.e. the data oscillate less). When considering $p = 2$ with a large time-step (Figure 7.6c), the contribution of the time-discretization, not included in our analysis, is clearly visible.

On Figure 7.7, we represent the true solution, the instantaneous error $\mathcal{E}_{\rho}(t)$, and $\eta(t)$ for different times $t$. We see that the estimators correctly locates the error at all times, although there is an underetimation close the abosrbing boundary. This is due to the spatial oscillation term that we have not included in the estimator (the right-hand side is not piecewise polynomial here). We also see that as the time increases, even if the error is correclty located, its magnitude becomes underestimated, which is in agreement with previous numerical observations and theoretical predictions. Again, this is because the estimator does not capture the (increasing) dispersion error.

### 7.4. Reflections by a penetrable obstacle.
We consider an example similar to the previous one, where an incident wave is reflected. As before, we set $\mathbf{d} = (\cos \theta, \sin \theta)$ with $\theta = 11\pi/8$, and consider the incident field $u_{\text{inc}} := v_{\sigma,d}$ for different values of $\sigma$. We then define $g := \dot{u}_{\text{inc}} + \nabla u_{\text{inc}} \cdot \mathbf{n}$ on $\Gamma_A$. The coefficients are define by

$$\mu := \begin{vmatrix} \mu_D & \text{in } D \\ 1 & \text{outside} \end{vmatrix} \quad \mathbf{A} := \begin{vmatrix} \mathbf{A}_D & \text{in } D \\ \mathbf{I} & \text{outside} \end{vmatrix},$$

with $\mu_D = 2$ and $\mathbf{A}_D = (1/2)\mathbf{I}$, where $D$ is the polygon defined by joining the vertices $(0, -0.5)$, $(0.5, 0.5)$, $(0, 0)$, $(0, 0.5)$ and $(-0.5, 0)$ (see Figure 7.9). As before, $T := 10$ and $\rho := 1/T$.

Here, the analytical solution is not available. As a result we limit our investigation to $p = 1$ with $\alpha = 1.5$. For a reference solution, we employ our numerical scheme with $p = 2$ on the
<figure><img src="image_6.png" /><figcaption>FIGURE 7.4. Error evolution in the reflection wave example for $\sigma = 0.5$, $p = 2$, $\alpha = 0.6$ and $h_{\max} = 0.05$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 7.5. Convergence in the reflection wave example for $\sigma = 0.5$</figcaption></figure>

same mesh with a time-step $\delta t$ divided by 3. As can be seen on Figures 7.8 and 7.9, the behaviour of the estimator is similar to the previous experiments, and in complete agreement with our analysis.

## 8. CONCLUSION

We present a construction of an equilibrated estimator for the scalar wave equation. Our construction avoids elliptic reconstructions, and is similar to [10]. The key novelty of our work is to employ a damped energy norm to measure the error together with a careful reliability and efficiency analysis providing a guaranteed and asymptotically constant-free upper bound as well as a polynomial-degree-robust lower bound. Numerical examples highlight the theory and suggest that it is sharp. This work is currently limited to the semi-discretization in space, and future work will be guided towards taking into account time-discretization, for instance
<figure><img src="image_8.png" /><figcaption>FIGURE 7.6. Effectivity index in the reflection example</figcaption></figure>

following [26, 28]. Besides, it would be interesting to provide equilibrated flux constructions that can operate in the presence of mass-lumping [16], or to consider discontinuous Galerkin schemes [21, 32].

REFERENCES

1. R. Adams and J. Fournier, *Sobolev spaces*, Academic Press, 2003.

2. S. Adjerid, *A posteriori finite element error estimation for second-order hyperbolic problems*, Comput. Meth. Appl. Mech. Engrg. **191** (2002), 4699–4719.

3. __________, *A posteriori error estimation for the method of lumped masses applied to second-order hyperbolic problems*, Comput. Meth. Appl. Mech. Engrg. **195** (2006), 4203–4219.

4. M. Ainsworth and J.T. Oden, *A posteriori error estimation in finite element analysis*, Wiley, 2000.

5. P.R. Amestoy, I.S. Duff, and J.Y. L'Excellent, *Multifrontal parallel distributed symmetric and unsymmetric solvers*, Comput. Methods Appl. Mech. Engrg. **184** (2000), 501–520.

6. R. Arcangeli and J.L. Gout, *Sur l'évaluation de l'erreur d'interpolation de lagrange dans un ouvert de Rⁿ*, R.A.I.R.O. Analyse numérique **10** (1976), no. 3, 5–27.

7. C. Baldassari, H. Barucq, H. Calandra, B. Denel, and J. Diaz, *Performance analysis of a high-order discontinuous Galerkin method application to the reverse time migration*, Commun. Comput. Phys. **11** (2012), no. 2, 660–673.
<figure><img src="image_9.png" /><figcaption>FIGURE 7.7. Reflection example with $\sigma = 2^{-3}$, $p = 2$, $\alpha = 0.6$, and $h = 0.05$</figcaption></figure>

8. W. Bangerth, M. Geiger, and R. Rannacher, *Adaptive Galerkin finite element methods for the wave equation*, Comput. Meth. Appl. Math. **10** (2010), no. 1, 3–4875–591.

9. W. Bangerth and R. Rannacher, *Adaptive finite element technique for the acoustic wave equation*, J. Comput. Acous. **9** (2001), no. 2, 575–591.

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

11. D. Boffi, F. Brezzi, and M. Fortin, *Mixed and hybrid finite element methods and applications*, Springer, 1991.

12. D. Braess, V. Pillwein, and J. Schöberl, *Equilibrated residual error estimates are p-robust*, Comput. Meth. Appl. Mech. Engrg. **198** (2009), 1189–1197.

13. E. Burman, O. Duran, A. Ern, and M. Steins, *Convergence analysis of hybrid high-order methods for the wave equation*, J. Sci. Comput. **87** (2021), no. 3.

14. T. Chaumont-Frelet, A. Ern, and M. Vohralík, *On the derivation of guaranteed and p-robust error estimates for the Helmholtz equation*, Numer. Math. **148** (2021), 525–573.
<figure><img src="image_10.png" /><figcaption>FIGURE 7.8. Errors and effectivity index in the obstacle example</figcaption></figure>

15. ________, Stable broken $H(\text{curl})$ polynomial extensions and p-robust a posteriori error estimates by broken patchwise equilibration for the curl-curl problem., Math. Comp. **91** (2022), 37–74.

16. G. Cohen, P. Joly, J.E. Roberts, and N. Tordjman, High order triangular finite element with mass lumping for the wave equation, SIAM J. Numer. Anal. **38** (2001), no. 6, 2047–2078.

17. P. Destuynder and B. Métivet, Explicit error bounds in a conforming finite element method, Math. Comp. **68** (1999), no. 228, 1379–1396.

18. C. Dobrzynski, MMG3D: User guide, Tech. Report 422, Inria, 2012.

19. W. Dörfler and S. Sauter, A posteriori error estimation for highly indefinite Helmholtz problems, Comput. Meth. Appl. Math. **13** (2013), 333–347.

20. A. Ern and M. Vohralík, A posteriori error estimation based on potential and flux reconstruction for the heat equation, SIAM J. Numer. Anal. **48** (2010), no. 1, 108–223.

21. ________, Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations, SIAM J. Numer. Anal. **53** (2015), no. 2, 1058–1081.

22. ________, Stable broken $H^1$ and $\mathbf{H}(\text{div})$ polynomial extensions for polynomial-degree-robust potential and flux reconstruction in three space dimensions, Math. Comp. **89** (2021), 551–594.

23. P. Fernandes and G. Gilardi, Magnetostatic and electrostatic problems in inhomogeneous anisotropic media with irregular boundary and mixed boundary conditions, Math. Meth. Appl. Sci. **47** (1997), no. 4, 2872–2896.

24. L. Fezoui, S. Lanteri, S. Lohrengel, and S. Piperno, Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes, ESAIM Math. Model. Numer. Anal. **39** (2005), no. 6, 1149–1176.

25. E.H. Georgoulis, O. Lakkis, and C. Makridakis, A posteriri $L^\infty(L^2)$-error bounds for finite element approximations to the wave equation, IMA J. Numer. Anal. **33** (2013), 1245–1264.

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

27. V. Girault and P.A. Raviart, Finite element methods for Navier-Stokes equations: theory and algorithms, Springer-Verlag, 1986.

28. 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.

29. R. Griesmaier and P. Monk, Discretization of the wave equation using continuous elements in time and a hybridizable discontinuous Galerkin method in space, J. Sci. Comput. **58** (2014), 472–498.

30. P. Grisvard, Elliptic problems in nonsmooth domains, Pitman, 1985.

31. ________, Singularities in boundary value problems, Springer-Verlag, 1992.

32. M.J. Grote, A. Schneebeli, and D. Schötzau, Discontinuous Galerkin finite element method for the wave equation, SIAM J. Numer. Anal. **44** (2006), no. 6, 2408–2431.
<figure><img src="image_11.png" /><figcaption>FIGURE 7.9. Obstacle example with $\sigma = 0.25$ and $h = 0.0125$</figcaption></figure>

33. J.S. Hesthaven and T. Warburton, *Nodal high-order methods on unstructured grids. Part I. Time-domain solution of Maxwell's equations*, J. Comput. Phys. **181** (2002), 1266–1288.

34. F. Jacobsen and P.M. Juhl, *Fundamentals of general linear acoustics*, Wiley, 2013.

35. X. Liu and F. Kikuchi, *Analysis and estimation of error constants for $P_0$ and $P_1$ interpolations over triangular finite elements*, J. Math. Sci. Univ. Tokyo **17** (2010), 27–78.

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

37. M. Picasso, *Numerical study of an anisotropic error estimator in the $L^2(H^1)$ norm for the finite element discretization of the wave equation*, SIAM J. Sci. Comput. **52** (2010), 2213–2234.

38. W. Prager and J.L. Synge, *Approximations in elasticity based on the concept of function space*, Quart. Appl. Math. **5** (1947), no. 3, 241–269.

39. F. Tantardini and A. Veeser, *The $L^2$-projection and quasi-optimality of Galerkin methods for parabolic equations*, SIAM J. Numer. Anal. **54** (2016), no. 1, 317–340.

40. R. Verfürth, *A posteriori error estimation techniques for finite element methods*, Oxford science publications, 2013.