# A posteriori error estimates and space-adaptive mesh refinements for time-dependent scattering problems

Théophile Chaumont-Frelet, Heiko Gimperlein, Ignacio Labarca-Figueroa, Jörg Nick

► To cite this version:

Théophile Chaumont-Frelet, Heiko Gimperlein, Ignacio Labarca-Figueroa, Jörg Nick. A posteriori error estimates and space-adaptive mesh refinements for time-dependent scattering problems. 2025. [https://arxiv.org/abs/2010.05335](https://arxiv.org/abs/2010.05335).

HAL Id: hal-05335116

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

Preprint submitted on 28 Oct 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
A POSTERIORI ERROR ESTIMATES AND SPACE-ADAPTIVE
MESH REFINEMENTS FOR TIME-DEPENDENT SCATTERING
PROBLEMS

THÉOPHILE CHAUMONT-FRELET, HEIKO GIMPERLEIN,
IGNACIO LABARCA-FIGUEROA, AND JÖRG NICK

ABSTRACT. This work studies a posteriori error estimates and their use for time-dependent acoustic scattering problems, formulated as a time-dependent boundary integral equation based on a single-layer ansatz. The integral equation is discretized by the convolution quadrature method in time and by boundary elements in space. We prove the reliability of an error estimator of residual type and study the resulting space-adaptive mesh refinements. Moreover, we present a simple modification of the convolution quadrature method based on temporal shifts, which recovers, for the boundary densities, the full classical temporal convergence order $2m - 1$ of the temporal convolution quadrature method based on the $m$-stage convolution quadrature semi-discretization. We numerically observe that the adaptive scheme yields asymptotically optimal meshes for an acoustic scattering problem in two dimensions.

# 1. INTRODUCTION AND PROBLEM FORMULATION

The computational treatment of time-domain scattering problems by formula-
tions based on boundary integral equations is by now a classical approach. Provably
convergent numerical methods of those formulations mostly fall into two categories
and either discretize the time-domain fundamental solution via a Galerkin approach
(referred to as space-time Galerkin), or rely on underlying time stepping schemes
to approximate inverse Laplace transforms, which leads to convolution quadrature
methods. Galerkin schemes with adaptive local mesh refinements in space, in time
or in space-time have been investigated recently, e.g., in [2, 3, 22, 23, 26, 39, 40, 41].
Beyond boundary integral equations, similar spatially adaptive methods for the
wave equation based on a finite element discretization in space and time stepping
schemes or the time-continuous problem in time have been considered in [16, 17].
In this paper, we initiate the study of a posteriori error estimates and the resulting
adaptive mesh refinement procedures based on convolution quadrature methods.
For these methods, Laplace domain operators have to be approximated at a large
number of wavenumbers. The computational bottleneck of these schemes is the as-
sembly and storage of the discretized Laplace domain boundary integral operators,
most commonly in the form of boundary element matrices.

Different approaches have been considered to reduce the number of necessary
frequencies at which these matrices have to be approximated. Convolution quadra-
ture schemes that use variable time step sizes have been introduced in [28, 29, 38]
and were further used to derive a temporally adaptive scheme in [39]. Other ap-
proaches to reduce the computational cost in the context of wave-types problems
derive higher-order methods [10, 6, 30] or even or even p-versions [37], or alternatively reduce the number of necessary evaluations through specialized quadrature for the approximation of highly oscillatory integrals [8, 33]. These variants of convolution quadrature methods reduce the number of frequencies at which the Laplace domain operators have to be approximated, often at the price of some computational overhead.

Unlike for space-time Galerkin methods, improved spatial discretizations, however, have attracted less attention. Efforts on improving the efficiency include compression schemes [7, 15] or wavelet based discretizations [18]. Based on the analysis by Plamenevskii and collaborators [27], it is known that for scattering problems in polygonal domains quasi-optimal approximations of the solution are achieved by time-independent graded meshes in space [35]. These results were later extended to further problems with geometric singularities and to approximations by hp methods [4, 13, 20, 21, 32, 34].

**Our contribution and outline.** Reducing the number of necessary spatial degrees of freedom to compute approximations of a desired accuracy is therefore a key opportunity to decrease both the runtime and the memory requirements of fully discrete schemes based on the combination of the convolution quadrature method with boundary elements. In this paper, we provide an a posteriori error analysis of the error introduced by the spatial boundary element discretization applied to time-dependent and time-discrete boundary integral equations. We derive error indicators that are, up to a data oscillation term, equivalent to the error. These bounds are formulated in Theorems 1–2, the main theoretical results in this article. Based on these error indicators, we derive a spatially adaptive scheme that is numerically shown to produce grids that recover the optimal rate of convergence with respect to the spatial degrees of freedom. This numerical evidence indicates that the proposed procedure produces quasi-optimal time-independent meshes [20, 35], although a proof of quasi-optimality of the generated meshes is beyond the scope of this paper and left for future research.

1.1. **Problem formulation.** Let $\Omega \subset \mathbb{R}^d$, $d = 2, 3$, be the exterior of a bounded Lipschitz domain or of a Lipschitz screen. On this exterior domain, we impose the acoustic wave equation

$$ (1) \qquad \partial_t^2 u^{\text{tot}} - \Delta u^{\text{tot}} = 0, \quad \text{in } \Omega^+. $$

Homogeneous Dirichlet boundary conditions are imposed along the boundary of the scatterers $\Omega$,

$$ (2) \qquad u^{\text{tot}} = 0, \quad \text{on } \Gamma = \partial\Omega. $$

Initially, the total wave is given by an incoming wave $u^{\text{inc}}$ with support away from the boundaries of the scatterers. The objective of the scattering problem is then to compute the unknown scattered wave $u$, which initially vanishes, such that the total wave $u^{\text{tot}} = u + u^{\text{inc}}$ fulfills the above boundary value problem.

## 2. THE LAPLACE DOMAIN PROBLEM

2.1. **Boundary integral formulation in the Laplace domain.** This section describes the setting of the Laplace domain problem, for some fixed frequency $s \in \mathbb{C}$ with positive real part.
Applying the Laplace transform to the acoustic wave equation reveals the Helmholtz problem, which reads for $s \in \mathbb{C}$ with $\mathrm{Re}\,s > 0$

$$ (3) \qquad s^2 \hat{u} - \Delta \hat{u} = 0, \quad \text{in } \mathbb{R}^d \setminus \mathcal{O}, $$

$$ (4) \qquad \hat{u} = -\hat{g}^{\text{inc}}, \quad \text{on } \Gamma = \partial \mathcal{O}. $$

For $|x| \to \infty$, an asymptotic condition is further imposed, that ensures the square integrability of $\hat{u}$. The fundamental solutions of the Helmholtz problem in two and three dimensions are given by

$$ G(s, x) = \frac{i}{4} H_0^{(1)}(isr), \quad x \in \mathbb{R}^2, \quad \text{resp.} \quad G(s, x) = \frac{e^{-s|x|}}{4\pi|x|}, \quad x \in \mathbb{R}^3, $$

in terms of the Hankel function $H_0^{(1)}$. We reformulate this problem as a boundary integral equation, fully formulated on the boundary of the scattering obstacle $\partial \mathcal{O}$. The potential operator reads

$$ S(s)\varphi(x) = \int_{\Gamma} G(s, x - y)\varphi(y), \quad x \in \mathbb{R}^d \setminus \mathcal{O}. $$

Taking the trace towards the boundary $x \to \Gamma$, which gives the single-layer boundary operator, which reads

$$ V(s)\varphi(x) = \int_{\Gamma} G(s, x - y)\varphi(y), \quad x \in \Gamma. $$

The solution of the above scattering problem (3) is then determined by

$$ \hat{u} = -S(s)V(s)^{-1}\hat{g}^{\text{inc}}. $$

The boundary operators fulfill s-explicit bounds, which are described in the following Lemma.

**Lemma 1.** [5, Proposition 3] The single-layer boundary operator is bounded by

$$ \|V(s)\|_{H^{1/2}(\Gamma) \leftarrow \tilde{H}^{-1/2}(\Gamma)} \leq C_{\Gamma} \frac{|s|^2 + 1}{\mathrm{Re}\,s}, $$

$$ \|V(s)^{-1}\|_{\tilde{H}^{-1/2}(\Gamma) \leftarrow H^{1/2}(\Gamma)} \leq C_{\Gamma} \frac{|s|^2 + 1}{\mathrm{Re}\,s}, $$

where C only depends on the geometry.

An explicit estimate on the constant $C_\Gamma$ in terms of the norms of trace operators with a proof is found in [36, Lemmas 2.7 & 2.8].

**2.2. Boundary element discretization.** The weak formulation of the boundary integral equation determining $\hat{\varphi}$ with test functions in the appropriate test space $\tilde{H}^{-1/2}(\Gamma)$ yields the following weak formulation: Find $\hat{\varphi} \in \tilde{H}^{-1/2}(\Gamma)$, such that for all $\hat{v} \in \tilde{H}^{-1/2}(\Gamma)$,

$$ \langle \hat{v}, V(s) \hat{\varphi} \rangle_{\Gamma} = - \langle \hat{v}, \hat{u}^{\text{inc}} \rangle_{\Gamma}. $$

Consider the space of continuous piecewise constant elements $X_h \subset \tilde{H}^{-1/2}(\Gamma)$. On this finite dimensional subspace, we formulate the Galerkin approximation: find $\hat{\varphi}_h \in X_h \subset \tilde{H}^{-1/2}(\Gamma)$, such that for all $\hat{v}_h \in X_h \subset \tilde{H}^{-1/2}(\Gamma)$

$$ \langle \hat{v}_h, V(s) \hat{\varphi}_h \rangle_{\Gamma} = - \langle \hat{v}_h, \hat{u}^{\text{inc}} \rangle_{\Gamma}. $$
Here, $\langle \cdot, \cdot \rangle_{\Gamma}$ denotes the extension of the $L^2$-pairing onto the $\tilde{H}^{-1/2}(\Gamma) \times H^{1/2}(\Gamma)$. Due to the positive real part of $s$, the bilinear form of the right-hand side is coercive, which enables a classical a priori error analysis. Using the Lax-Milgram Lemma on the subspace $X_h$ shows that the spatially discrete resolvent $V_h(s)^{-1}$, defined by

$$ \hat{\varphi}_h = -V_h(s)^{-1} \hat{u}^{\text{inc}}. $$

inherits the Laplace domain bounds of the spatially continuous single-layer operator, in particular

$$ (5) \qquad \|V_h(s)^{-1}\|_{\tilde{H}^{-1/2}(\Gamma) \leftarrow H^{1/2}(\Gamma)} \le C_{\Gamma} \frac{|s|^2 + 1}{\text{Re } s}. $$

2.3. **A posteriori error estimate.** We then find the following time-harmonic a posteriori estimate.

$$ (6) \qquad \begin{aligned} \|\hat{\varphi} - \hat{\varphi}_h\|_{\tilde{H}^{-1/2}(\Gamma)} &= \|V(s)^{-1} (V(s)\hat{\varphi} - V(s)\hat{\varphi}_h)\|_{\tilde{H}^{-1/2}(\Gamma)} \\ &\le C_{\Gamma} \frac{|s|^2 + 1}{\text{Re } s} \|\hat{f} - V(s)\hat{\varphi}_h\|_{H^{1/2}(\Gamma)}. \end{aligned} $$

Let $\hat{R}_h = \hat{f} - V(s)\hat{\varphi}_h \in H^{1/2}(\Gamma)$ denote the time-harmonic residual. We note that by following the argument structure above, we further obtain

$$ (7) \qquad \begin{aligned} \|R_h(s)\|_{H^{1/2}(\Gamma)} &= \|V(s) (V(s)^{-1} \hat{f} - \hat{\varphi}_h)\|_{H^{1/2}(\Gamma)} \\ &\le C_{\Gamma} \frac{|s|^2 + 1}{\text{Re } s} \|\hat{\varphi} - \hat{\varphi}_h\|_{\tilde{H}^{-1/2}(\Gamma)}. \end{aligned} $$

Overall, we therefore arrive at the following equivalence of the norms of the error and of the residual, which holds up to a factor depending on the Laplace parameter $s$

$$ (8) \qquad \frac{\text{Re } s}{C_{\Gamma}(|s|^2 + 1)} \|R_h(s)\|_{H^{1/2}(\Gamma)} \le \|\hat{\varphi} - \hat{\varphi}_h\|_{\tilde{H}^{-1/2}(\Gamma)} \le C_{\Gamma} \frac{|s|^2 + 1}{\text{Re } s} \|R_h(s)\|_{H^{1/2}(\Gamma)}. $$

### 3. THE TIME-DEPENDENT SCATTERING PROBLEM

3.1. **Temporal convolution operators.** We describe the techniques to transfer the Laplace domain results to the time-domain. Let $V$ be an Hilbert space and consider an analytic family of bounded linear operators $K(s): V \to W$, which are defined for $\text{Re } s > 0$. We assume that $K$ is polynomially bounded with respect to the Laplace domain parameter $s$, in the following sense: There exist a real $\kappa$, $\nu \ge 0$ and, for every $\sigma > 0$, there exists $M_{\sigma} < \infty$, such that

$$ (9) \qquad \|K(s)\|_{W \leftarrow V} \le M_{\sigma} \frac{|s|^{\kappa}}{(\text{Re } s)\nu}, \quad \text{Re } s \ge \sigma > 0. $$

For a sufficiently regular function $g: [0, T] \to X$, which vanishes together with its first $m > \kappa$ derivatives, we write

$$ (10) \qquad K(\partial_t)g = (\mathcal{L}^{-1}K) * g $$

for the temporal convolution of the inverse Laplace transform of $K$ with $g$, where $g$ is extended by zero on the negative real axis. Due to the associativity of convolution
we obtain, for two such families of operators $K(s)$ and $L(s)$ mapping into compatible spaces, the composition rule

$$ (11) \qquad K(\partial_t)L(\partial_t)g = (KL)(\partial_t)g. $$

We denote the Sobolev space of real order $r$ of $V$-valued functions on $\mathbb{R}$ by $H^r(\mathbb{R}, V)$ and, on finite intervals $(0, T)$, we write

$$ H_0^r(0, T; V) = \{g|_{(0,T)} : g \in H^r(\mathbb{R}, V) \text{ with } g = 0 \text{ on } (-\infty, 0)\}. $$

The natural norm on $H_0^r(0, T; V)$ is, for $r \ge 0$, equivalent to the norm $\|\partial_t^r g\|_{L^2(0,T;V)}$. By the Plancherel formula, the temporal convolution operators $K(\partial_t)$ are bounded on the appropriate temporal Sobolev spaces. [31, Lemma 2.1]: If $K(s)$ is bounded by (9) in the half-plane $\Re s > 0$, then $K(\partial_t)$ fulfills

$$ (12) \qquad \|K(\partial_t)\|_{H_0^r(0,T;V') \to H_0^{r+\kappa}(0,T;V)} \le eM_{1/T}T^\nu $$

for arbitrary real $r$. (The bound on the right-hand side arises from the bound $e^{\sigma T}M_\sigma$ on choosing $\sigma = 1/T$.) We note that for any integer $k \ge 0$ and real $\alpha > \frac{1}{2}$, we have the continuous embedding $H_0^{k+\alpha}(0, T; V) \subset C^k([0, T]; V)$.

**3.2. Time-dependent boundary integral equations.** We turn towards the time-dependent scattering problem, whose solution is determined by the boundary integral equation

$$ V(\partial_t)\varphi = -u^{\text{inc}}. $$

The scattered wave can be evaluated, in the exterior domain $\mathbb{R}^3 \setminus \mathcal{O}$, by the representation formula

$$ u = S(\partial_t)\varphi. $$

By (12) and Laplace domain bounds, we obtain bounds on the (temporal) regularity of the scattered wave $u$, under assumptions on the incoming wave $u^{\text{inc}}$.

**3.3. Boundary element discretization.** The spatially discrete time-dependent boundary density is constructed by employing the Galerkin approach on the temporal convolution equation, which yields the following formulation: Find the spatially discrete boundary density $\varphi_h(t) \in X_h$ (for $t \in [0, T]$), such that

$$ \langle v_h, V(\partial_t)\varphi_h \rangle_{\Gamma} = - \langle v_h, \hat{u}^{\text{inc}} \rangle_{\Gamma} \quad \text{for all } v_h \in X_h. $$

An equivalent and compact characterization of this approximation can be formulated with the temporal convolution operator associated to $V_h(s)^{-1}$, namely we have

$$ \varphi_h = -V_h(\partial_t)^{-1}u^{\text{inc}}. $$

We then have the following a posteriori error estimate of the spatially discrete system.

**Proposition 1.** Let $\omega > 0$ be an arbitrary constant. Then, the error of the spatial-semidiscretization is bounded, for all $r \in \mathbb{N}$, by

$$ \int_0^\infty \|\varphi - \varphi_h\|_{L^2(0,T;\tilde{H}^{-1/2}(\Gamma))} \le C(1+\omega^2) \|f - V(\partial_t)\varphi_h\|_{L^2(0,T;H^{1/2}(\Gamma))} \\ + \frac{C}{\omega^r} \|f\|_{H_0^{r+6}(0,T;H^{1/2}(\Gamma))}, $$

if the right-hand side is finite.
Proof. We start the proof by showing the stated estimate on a weighted $L^2$-norm on the positive half line.

(i) Proving the estimate in the $L^2(\mathbb{R}_+, e^{-2\sigma t}dt)$-norm. In the following, we choose $\sigma = 1/T$ and use the Plancherel formula several times. We start from the standard estimate

$$ \int_0^\infty e^{-2\sigma t} \|\varphi - \varphi_h\|_{\dot{H}^{-1/2}(\Gamma)}^2 dt \le C_T T \left( \int_0^\infty e^{-2\sigma t} \|\partial_t^2 f - V(\partial_t)\partial_t^2 \varphi_h\|_{H^{1/2}(\Gamma)}^2 \right. \\ \left. \qquad + e^{-2\sigma t} \|f - V(\partial_t)\varphi_h\|_{H^{1/2}(\Gamma)}^2 dt \right). $$

We note that the residual is measured in the temporal $H^2$ norm, which we intend to avoid, at the cost of introducing a data oscillation term. Following [16], we split the corresponding contour integral in the frequency domain and obtain, for any $\omega > 0$, the estimate (for $\xi_h = f - V(\partial_t)\varphi_h$)

$$ \begin{aligned} \int_0^\infty e^{-2\sigma t} \|\partial_t^2 \xi_h\|_{H^{1/2}(\Gamma)}^2 \le \int_{\sigma+i\mathbb{R}} |s|^4 \|\hat{\xi}_h\|_{H^{1/2}(\Gamma)}^2 ds \\ \le \omega^4 \int_{\sigma+i\mathbb{R}} \|\hat{\xi}_h\|_{H^{1/2}(\Gamma)}^2 \chi_{|s|<\omega} ds + \frac{1}{\omega^{2r}} \int_{\sigma+i\mathbb{R}} |s|^{4+2r} \|\hat{\xi}_h\|_{H^{1/2}(\Gamma)}^2 \chi_{|s|>\omega} ds \\ \le \omega^4 \int_0^\infty e^{-2\sigma t} \|\xi_h\|_{H^{1/2}(\Gamma)}^2 dt + \frac{1}{\omega^{2r}} \int_0^\infty e^{-2\sigma t} \|\partial_t^{r+2} \xi_h\|_{H^{1/2}(\Gamma)}^2 dt. \end{aligned} $$

Under sufficient temporal regularity of the right-hand side, we can therefore split the error estimate into

$$ \int_0^\infty e^{-2\sigma t} \|\varphi - \varphi_h\|_{\dot{H}^{-1/2}(\Gamma)}^2 dt \le C(1+\omega^2) \int_0^\infty e^{-2\sigma t} \|f - V(\partial_t)\varphi_h\|_{H^{1/2}(\Gamma)}^2 dt \\ + \frac{C}{\omega^r} \int_0^\infty e^{-2\sigma t} \|\partial_t^{r+2} f - V(\partial_t)\partial_t^{r+2} \varphi_h\|_{H^{1/2}(\Gamma)}^2 dt. $$

The final term can be estimated by

$$ \int_0^\infty e^{-2\sigma t} \|\partial_t^{r+2} (I - V(\partial_t)V_h^{-1}(\partial_t))f\|_{H^{1/2}(\Gamma)}^2 dt \le C \int_0^\infty e^{-2\sigma t} \|\partial_t^{r+6} f\|_{H^{1/2}(\Gamma)}^2 dt, $$

which gives the stated result in the weighted norm on $\mathbb{R}_+$.

(ii) Transferring the result to $L^2(0, T)$. Stein's extension theorem [1, Theorem 5.24], shows the existence of an extension operator that is bounded by

$$ \|E\|_{H_0^k(\mathbb{R}_+; H^{1/2}(\Gamma)) \leftarrow H_0^k(0, T; H^{1/2}(\Gamma))} \le C, $$

for arbitrary $k$. We extend both $\xi_h$ and $f$ with this operator and further extend $\varphi_h$ and $\varphi$ by defining

$$ \tilde{\varphi} = V(\partial_t)^{-1}(E(f)) \quad \text{and,} \quad \tilde{\varphi}_h = V(\partial_t)^{-1}(E(f) - E(\xi_h)). $$

By causality, these functions coincide with $\varphi$ and $\varphi_h$ on $(0, T)$. Applying the bound of part (i) and the continuity of the extension operator now yields the stated result
by estimating

$$
\begin{align*}
\|\varphi - \varphi_h\|_{L^2(0,T;\tilde{H}^{-1/2}(\Gamma))}^2 &\le e \int_0^\infty e^{-2\sigma t} \|\tilde{\varphi} - \tilde{\varphi}_h\|_{\tilde{H}^{-1/2}(\Gamma)}^2 \, dt \\
&\le C \int_0^\infty e^{-2\sigma t} (1+\omega^2) \|E(\xi_h)\|_{H^{1/2}(\Gamma)}^2 + \frac{e^{-2\sigma t}}{\omega^{2r}} \|\partial_t^{r+6} E(f)\|_{H^{1/2}(\Gamma)}^2 \, dt \\
&\le C(1+\omega^2) \|f - V(\partial_t)\varphi_h\|_{L^2(0,T;H^{1/2}(\Gamma))}^2 + \frac{C}{\omega^r} \|f\|_{H_0^{r+6}(0,T;H^{1/2}(\Gamma))}^2.
\end{align*}
$$

□

4. FULL DISCRETIZATION

In the following, we extend our analysis to temporally discrete systems.

4.1. **Time discretization by Convolution Quadrature.** We recall the approximation of temporal convolutions $K(\partial_t)g$ by the convolution quadrature method based on underlying time stepping scheme and fix the associated notation. With the constant time step size $\tau > 0$, we construct approximations $K(\partial_t^r)g^n$ to $K(\partial_t)g(t_n)$ at equidistant time points $t_n = n\tau$. Most effective convolution quadrature schemes are based on Runge-Kutta methods, which are fully determined by their Butcher-tableau

$$
\mathcal{A} = (a_{ij})_{i,j=1}^{m}, \quad b = (b_1, \dots, b_m)^T, \quad \text{and} \quad c = (c_1, \dots, c_m)^T.
$$

The associated stability function of the Runge–Kutta method is $R(z) = 1 + zb^T(I - z\mathcal{A})^{-1}1$, where $1 = (1, 1, \dots, 1)^T \in \mathbb{R}^m$. Consider an temporal convolution operator in the sense of Section 3.1, namely $K(\partial_t) : H_0^{r+\kappa}(0, T; V) \to H_0^r(0, T; W)$ for arbitrary real $r$. Moreover, let $g : [0, T] \to V$ be a time-dependent function that is sufficiently regular and vanishes, together with sufficiently many derivatives at $t=0$ for the expression (10) to be well-defined. The convolution $(K(\partial_t)g)(t)$ is approximated at the stages of the Runge–Kutta methods, namely at the discrete times

$$
t_n = (t_n + c_l \tau)_{l=1}^{m}, \text{ where } t_n = n\tau.
$$

The Runge–Kutta differentiation symbol is defined by the expression

$$
\Delta(\zeta) := \left( \mathcal{A} + \frac{\zeta}{1-\zeta} \mathbb{1} b^T \right)^{-1} \in \mathbb{C}^{m \times m}, \quad \zeta \in \mathbb{C} \text{ with } |\zeta| < 1.
$$

As a consequence of the Sherman–Woodbury formula, this expression is well defined for $|\zeta| < 1$ if $R(\infty) = 1 - b^T \mathcal{A}^{-1} 1$ satisfies $|R(\infty)| \le 1$. For A-stable Runge–Kutta methods (e.g. the Radau IIA methods), the eigenvalues of the matrices $\Delta(\zeta)$ have positive real part for $|\zeta| < 1$ (see [10, Lem. 3]). The Sherman–Morrison formula then yields the expression

$$
\Delta(\zeta) = \mathcal{A}^{-1} - \frac{\zeta}{1 - R(\infty)\zeta} \mathcal{A}^{-1} 1 b^T \mathcal{A}^{-1}.
$$

We are now in a position to define the convolution quadrature weights $W_n(K) : V^m \to W^m$.

We replace the complex argument $s$ in $K(s)$ by the matrix-valued analytic function $\Delta(\zeta)/\tau$ and write down the power series expansion

$$
K\left(\frac{\Delta(\zeta)}{\tau}\right) = \sum_{n=0}^{\infty} W_n(K)\zeta^n.
$$
In the following, we use an upper index to denote a sequence element with $m$ components. Thus, for a sequence $g = (g^n)$ with $g^n = (g_l^n)_{l=1}^m \in V^m$ we arrive at the discrete convolution denoted by

$$ (13) \qquad (K(\partial_t^\tau)g)^n := \sum_{j=0}^{n} W_{n-j}(K)g^j \in W^m. $$

The notation $K(\partial_t^\tau)g$ in (13) indicates that the resulting vector contains approximations at the stages $t_n$. For functions $g: [0, T] \to V$, we use this notation for the vectors $g^n = g(t_n) = (g(t_n + c_i\tau))_{i=1}^m$ of values of $g$. The $\ell$-th component of the vector $(K(\partial_t^\tau)g)^n$, that we denote by $(K(\partial_t^\tau)g)^{n,\ell}$, is then an approximation to $(K(\partial_t)g)(t_n + c_\ell\tau)$, i.e. $(K(\partial_t^\tau)g)^{n,\ell} \approx (K(\partial_t)g)(t_n + c_\ell\tau)$ (see [9, Thm. 4.2]). In particular, if $c_m = 1$, as is the case with stiffly stable Runge-Kutta methods, which includes the Radau IIA methods, the continuous convolution at $t_n$ is approximated by the $m$-th, i.e. last component of the $m$-vector (13) for $n-1$:

$$ (K(\partial_t^\tau)g)_n := (K(\partial_t^\tau)g)^{n-1,m} \in W. $$

These components approximate the continuous convolution at the equidistant time points $t_n$, i.e. $(K(\partial_t^\tau)g)_n \approx (K(\partial_t)g)(t_n)$.

A property that is key to show stability in many settings is that the composition rule (11) is preserved under this discretization: For two compatible operator families $K(s)$ and $L(s)$, we have

$$ K(\partial_t^\tau)L(\partial_t^\tau)g = (KL)(\partial_t^\tau)g. $$

Such a property can only be formulated for the vector valued discrete convolution, which includes the approximations at the stages $t_n$, but can not be formulated for the approximation $K(\partial_t^\tau)g$ at the equidistant time points $t_n$.

**Remark 1** (Multistep based convolution quadrature). The setting of multistep-based convolution quadrature is the technically simplest setting. The convolution quadrature weights are then determined by the power series

$$ K\left(\frac{\delta(\zeta)}{\tau}\right) = \sum_{n=0}^{\infty} W_n(K)\zeta^n, $$

where $\delta$ is the scalar-valued characteristic function of the multistep method, which is in the case of the BDF-p method given by

$$ \delta(\zeta) = \sum_{\ell=1}^{p} \frac{1}{\ell} (1 - \zeta)^{\ell}, \quad \text{for } p = 1, 2. $$

For the method to be A-stable, we require $p \le 2$.

The following error bound for Runge-Kutta convolution quadrature from [10, Thm. 3], here directly stated for the Radau IIA methods [25, Sec. IV.5] and transferred to a Hilbert space setting, will be the basis for our error bounds of the time discretization.

**Lemma 2.** Let $K(s): V \to W$, $s \ge \sigma > 0$, be an analytic family of linear operators between Banach spaces $V$ and $W$ satisfying the bound (9) with exponents $\kappa$ and $\nu$. Consider the Runge-Kutta convolution quadrature based on the Radau IIA method with $m$ stages. Let $r > \max(2m + \kappa, 2m - 1, m + 1)$ and further let $g \in C^r([0, T], V)$
satisfy $g(0) = g'(0) = \dots = g^{(r-1)}(0) = 0$. Then, the following error bound holds at $t_n = n\tau \in [0, T]$:

$$ \begin{aligned} & \|(K(\partial_t^\tau)g)_i^n\}_{i=1,\dots,m} - (K(\partial_t)g)(\underline{t}_n)\|_{W^m} \\ & \le C M_{1/T} \tau^{\min(2m-1,m+1-\kappa+\nu)} \left( \|g^{(r)}(0)\|_V + \int_0^t \|g^{(r+1)}(t')\|_V dt' \right). \end{aligned} $$

The constant $C$ is independent of $\tau$ and $g$ and $M_\sigma$ of (9), but depends on the exponents $\kappa$ and $\nu$ in (9) and on the final time $T$.

### 4.2. Modified formulation.
In the context of transfer operators that are bounded only by a high power of $|s|$ (i.e. that fulfill (9) with a moderately large $\kappa$), the Runge-Kutta convolution quadrature approximation suffers from order reductions. We use the following modification to the convolution quadrature method to overcome this difficulty. Let $g \in H_0^1(\mathbb{R}_+; V)$ be some arbitrary density, for which its Laplace transform

$$ \mathcal{L}(g(\cdot - \eta))(s) = \int_0^\infty e^{-st} g(t - \eta) dt = e^{-s\eta} \int_{-\eta}^\infty e^{-st} g(t) dt = e^{-s\eta} \mathcal{L}(g)(s). $$

We observe that the $\eta$-shift in time-domain is a multiplication with $J_\eta(s) = e^{-s\eta}$ in the Laplace domain. The boundary integral equation can therefore be modified,

$$ \varphi = (J_{\eta}(\partial_t)V(\partial_t)^{-1})J_{-\eta}(\partial_t)f. $$

The modification

$$ (14) \qquad J_{-\eta}(\partial_t)f(t) = \begin{cases} 0, & \text{for } t \le \sigma, \\ f(t+\eta), & \text{for } t > \sigma. \end{cases} $$

can be realized exactly, by appending the function by zeros (and is valid since $f$ is assumed to vanish on the negative real-axis). In the Laplace domain, we now have that

$$ (15) \qquad \|J_{\eta}(s)V(s)^{-1}\|_{\bar{H}^{-1/2}(\Gamma)\leftarrow H^{1/2}(\Gamma)} \le C_{\Gamma} \frac{|s|^2+1}{\text{Re } s} e^{-\eta \text{Re } s}. $$

The shifted scheme then reads

$$ (16) \qquad \varphi_{\eta}^{\tau} = (V(\partial_t^{\tau})^{-1}) J_{\eta}(\partial_t) J_{-\eta}(\partial_t) \underline{f}. $$

Note that in order to construct a useful solution to the original convolution equation, we have to compute the solution to the convolution equation (16) until the final time $T + \eta$ and then disregard the numerical solution on the extended interval $[T, T + \eta]$. This modification has two desirable effects:

**Full classical convergence rate.** The exponential decay in the Laplace domain estimate (15) improves corresponding the error bound Lemma 2 and shows that the shifted density $\varphi_\eta^\tau$ converges at the full classical order $\tau^{2m-1}$. Without any shift, i.e. $\eta = 0$, this reduces to the stage order $\tau^m$.

**Fast evaluation due to the truncation of damped wave numbers.** The “all-at-once” implementation of the convolution quadrature method reduces the computational bottleneck to the evaluation of the Laplace domain operator $J_\eta(s)V(s)^{-1}$ at several different frequencies $s_j$ (see e.g. [12, Section 3.2]). Due to the exponential decay in the bound (15), the operator can be truncated for Laplace domain arguments $s$ with real part that is larger than some predetermined constant.
4.3. **Fully discrete a posteriori error analysis.** For arbitrary real-valued order $r > 0$, we have the following chain of equations

$$ (17) \qquad \begin{aligned} \|\varphi^\tau - \varphi_h^\tau\|_{H_0^r(0,T;\tilde{H}^{-1/2}(\Gamma))} &= \left\|V(\underline{\partial}_t^\tau)^{-1} \left(V(\underline{\partial}_t^\tau)\varphi^\tau - V(\underline{\partial}_t^\tau)\varphi_h^\tau\right)\right\|_{H_0^r(0,T;\tilde{H}^{-1/2}(\Gamma))} \\ &\le C \left\|\left(\underline{\partial}_t^\tau\right)^2 \left(J_{-\eta}(\underline{\partial}_t^\tau)J_\eta(\partial_t)\underline{f} - V(\underline{\partial}_t^\tau)\varphi_h^\tau\right)\right\|_{H_0^r(0,T;\tilde{H}^{-1/2}(\Gamma))}. \end{aligned} $$

Here, we used the operator valued bound

$$ \left\| V(\underline{\partial}_t^\tau)^{-1}(\underline{\partial}_t^\tau)^{-2} \right\|_{H_0^r(0,T;\tilde{H}^{-1/2}(\Gamma)) \leftarrow H_0^r(0,T;H^{1/2}(\Gamma))} \le C. $$

By combining the estimate with the temporal semi-discrete error bound yields the following error estimate.

**Theorem 1.** Under temporal regularity on the exact solution, we have the estimate

$$ \begin{aligned} \|\varphi - \varphi_h^\tau\|_{H_0^r(0,T;\tilde{H}^{-1/2}(\Gamma))} \le C \left\| \left(\underline{\partial}_t^\tau\right)^2 \left( J_{-\eta}(\underline{\partial}_t^\tau)J_\eta(\partial_t)\underline{f} - V(\underline{\partial}_t^\tau)\varphi_h^\tau \right) \right\|_{H_0^r(0,T;\tilde{H}^{-1/2}(\Gamma))} \right. \\ \left. + C\tau^{2m-1}, \right\} $$

where the constant $C$ depends on the geometry of $\Omega$, polynomially on the final time $T$ and on regularity of $f$.

For the stages, we further obtain

$$ \left\| \underline{\varphi} - \underline{\varphi}_h^\tau \right\|_{H_0^r(0,T;\tilde{H}^{-1/2}(\Gamma))} \le C \left\| \left( \underline{\partial}_t^\tau \right)^2 \left( \underline{f} - V(\underline{\partial}_t^\tau)\varphi_h^\tau \right) \right\|_{H_0^r(0,T;\tilde{H}^{-1/2}(\Gamma))} + C\tau^{m+1}. $$

4.4. **A posteriori error analysis of BDF methods.** In the case of multistep methods, we can improve the a posteriori error estimate by using the stability of the convolution quadrature method.

**Theorem 2.** Let $\omega > 0$ be an arbitrary constant. Then, the space discretization error of the temporal semi-discretization (with the BDF-1 or the BDF-2 method) is bounded, for all $r \in \mathbb{N}$, by

$$ \begin{aligned} \|\varphi^\tau - \varphi_h^\tau\|_{L^2(0,T;\tilde{H}^{-1/2}(\Gamma))} \le C(1+\omega^2) \|f - V(\underline{\partial}_t^\tau)\varphi_h^\tau\|_{L^2(0,T;H^{1/2}(\Gamma))} \\ + \frac{C}{\omega^r} \|f\|_{H_0^{r+6}(0,T;H^{1/2}(\Gamma))}, \end{aligned} $$

if the right-hand side is finite.
Proof. In the following, we write $R_h^\tau = f - V(\partial_t^\tau)\varphi_h^\tau$. We (formally) start with

$$
\begin{align*}
\int_0^\infty e^{-2\sigma t} \|\varphi^\tau - \varphi_h^\tau\|_{\tilde{H}^{-1/2}(\Gamma)}^2 \, dt &= \int_{\sigma+i\mathbb{R}} \|\hat{\varphi}^\tau - \hat{\varphi}_h^\tau\|_{\tilde{H}^{-1/2}(\Gamma)}^2 \, ds \\
&\le C \int_{\sigma+i\mathbb{R}} \left|\left(\frac{\delta(e^{-s\tau})}{\tau}\right)^2\right|^2 \left\|\hat{f} - V\left(\frac{\delta(e^{-s\tau})}{\tau}\right) \hat{\varphi}_h^\tau\right\|_{\tilde{H}^{-1/2}(\Gamma)}^2 \, ds \\
&\le C \sup_{s \in \sigma+i\mathbb{R}, |s| < \omega} \left|\left(\frac{\delta(e^{-s\tau})}{\tau}\right)^2\right|^2 \int_{\sigma+i\mathbb{R}} \left\|\hat{f} - V\left(\frac{\delta(e^{-s\tau})}{\tau}\right) \hat{\varphi}_h^\tau\right\|_{\tilde{H}^{-1/2}(\Gamma)}^2 \chi_{|s| < \omega} \, ds \\
&\quad + C \frac{1}{\omega^{2r}} \int_{\sigma+i\mathbb{R}} |s|^{2r} \left|\left(\frac{\delta(e^{-s\tau})}{\tau}\right)^2\right|^2 \left\|\hat{f} - V\left(\frac{\delta(e^{-s\tau})}{\tau}\right) \hat{\varphi}_h^\tau\right\|_{\tilde{H}^{-1/2}(\Gamma)}^2 \chi_{|s| > \omega} \, ds, \\
&\le C \omega^4 \int_0^\infty e^{-2\sigma t} \|f - V(\partial_t^\tau)\varphi_h^\tau\|_{\tilde{H}^{-1/2}(\Gamma)}^2 \, dt + \frac{C}{\omega^{2r}} \int_{\sigma+i\mathbb{R}} |s|^{2r+4+8} \left\|\hat{f}\right\|_{\tilde{H}^{-1/2}(\Gamma)}^2 \, ds.
\end{align*}
$$

The last estimate holds only in the case of $\delta(\zeta) = 1 - \zeta$ (implicit Euler method), or the BDF-2 method, with [11, Lemma 3.3]. The bound can then be carried over to the finite time interval $(0, T)$ by the extension argument of the proof of Proposition 1. $\square$

4.5. Localized error estimator. As in [24] (pp. 340-342 for 2d, pp. 349-352 for 3d), the right hand side is estimated by a local, computable quantity:

$$
\begin{align*}
& \left\| (\partial_t^\tau)^k f - V(\partial_t^\tau) (\partial_t^\tau)^k \varphi_h^\tau \right\|_{t^2(0,T;H^{1/2}(\Gamma))}^2 \\
& \le \tilde{C}_\Gamma \tau \sum_{t_i \in E_j} \sum_{E_j} h_{E_j} \int_{E_j} \left( (\partial_t^\tau)^k \nabla_\Gamma (f(t_i, x) - V(\partial_t^\tau) \varphi_h^\tau(t_i, x)) \right)^2 dx,
\end{align*}
$$

where $k \in \mathbb{N}_0$, $h_{E_j}$ is the diameter of the element $E_j$ and $t_i = i\tau$. If relevant, the proof also extends to $H^s(\Gamma)$ for $s \in [0, 1]$.

The space-adaptive algorithm relies on the error indicator

$$
\eta(E_j)^2 = \tau h_{E_j} \sum_{k \in \{0,2\}} \sum_{t_i} \int_{E_j} \left( (\partial_t^\tau)^k \nabla_\Gamma (f(t_i, x) - V(\partial_t^\tau) \varphi_h^\tau(t_i, x)) \right)^2 dx
$$

in the element $E_j$ of the spatial mesh. Here, $\nabla_\Gamma$ denotes the surface gradient on $\Gamma$.

5. ADAPTIVE ALGORITHM AND IMPLEMENTATION

We refer the reader to [22] for more details on adaptive schemes, there presented
in the context of a space-time Galerkin method. We follow the classical construc-
tion of adaptive finite element methods, which has the structure:

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

More precisely, we iterate over the following steps, until the error indicator falls
below a target tolerance:

(1) Solve the boundary integral equation

$$
\langle v_h, V(\partial_t^\tau) \varphi_h^\tau \rangle_\Gamma = \langle v_h, J_\eta(\partial_t^\tau) J_{-\eta}(\partial_t) f \rangle_\Gamma.
$$
(2) Compute error indicators $\eta(E_j)$,

(3) Use the Dörfler marking criterion [19] to mark those elements that have a relatively large error indicator, characterized that for some predetermined $\theta \in (0, 1)$ it holds that

$$\eta(E_j) > \theta\eta_{\max},$$

where $\eta_{\max}$ is the maximal error indicator of the mesh.

(4) Refine all marked elements $E_j$.

The practical computation and implementation of the error estimator involves some numerical challenges, for which we follow [14, Section 7.1].

To be specific, for a given discrete solution $\varphi_h^\tau(t_i) \in X_h$ at time step $t_i = i\tau$ we project $V(\partial_t^\tau)\varphi_h^\tau$ to a function $P_h(t_i)$ in the finite element space of piecewise linear functions $Y_h \subset H^1(\Gamma)$. We also determine the best approximation $f_h(t_i) \in Y_h$ for the right-hand side $f(t) \in H^{1/2}(\Gamma)$. Then $R_h^\tau(t_i) = P_h(t_i) - f_h(t_i) \in Y_h$. For every element $E_j$ of the spatial mesh we approximate

$$(18) \qquad \eta(E_j)^2 \approx \eta_h(E_j)^2 := \tau h_{E_j} \sum_{t_i} \int_{E_j} |\nabla_{\Gamma} R_h^\tau(t_i, x)|^2 dx,$$

where $\eta_h(E_j)$ can be computed exactly by numerical quadrature. In the numerical experiments below, 5 Gauss quadrature points are used per element.

## 6. NUMERICAL EXPERIMENTS

For the numerical experiments we consider a right hand side corresponding to an incoming plane wave. Define

$$f(t) = \sin(\omega(t - t_{\text{lag}}))H(t)H(L - t), \quad \omega = 2, L = 2,$$

where

$$H(t) = 1 - \frac{1}{1 + \exp(\beta t)}, \quad \beta = 5,$$

and $t_{\text{lag}} = 4$. For a given direction $d \in \mathbb{R}^2$ we define $u^{\text{inc}}$ as the plane wave

$$u^{\text{inc}}(x) := f(x \cdot d - t), \quad \text{for all } x \in \mathbb{R}^2.$$

We solve the Dirichlet scattering problem up to the final time $T = 10$ with a Galerkin boundary-element discretization in space with $\mathcal{P}_0$-basis functions and RKCQ with a two-step Radau IIA in time. The following schemes are compared: uniform refinement; adaptive refinement (described in Section 5); $\beta$-graded meshes with grading parameter $\beta \in \{2, 3\}$ [24]. Point evaluations are computed using the representation formula at the points $x_1 = (2, 2)$, $x_2 = (-2, 2)$, $x_3 = (-2, -2)$ and $x_4 = (2, -2)$. The corresponding error is computed as the square root of the sum of the squares of errors in each of these points.

We measure the errors in the energy space $L^2(0, T; \tilde{H}^{-1/2}(\Gamma))$ using the equivalent norm

$$(19) \qquad \|\varphi_h^\tau\|_{L^2(0,T; \tilde{H}^{-1/2}(\Gamma))}^2 = \int \|\varphi_h^\tau(t)\|_{\tilde{H}^{-1/2}(\Gamma)}^2 \, dt = \int \langle \varphi_h^\tau(t), V(1)\varphi_h^\tau(t) \rangle_{\Gamma} \, dt \\ \approx \tau \sum_{t_i} \langle \varphi_h^\tau(t_i), V(1)\varphi_h^\tau(t_i) \rangle_{\Gamma},$$
where $V(1)$ is the single-layer boundary integral operator with $s = 1$.

**6.1. Flat Screen.** We solve the problem at a flat screen $\Gamma := (-1, 1) \times \{0\}$. We show a sequence of meshes obtained by adaptive refinement in Figure 1. Results for convergence of point evaluation and convergence of densities in energy norm (19) are presented in Figure 2. The error estimator underestimates the error roughly by an order of magnitude. Nevertheless, we observe that the adaptive refinement displays the same behavior as the asymptotically optimal 3-graded mesh refinement. Figure 5 shows snapshots of the density $\varphi_h^T$ at times $t = 2, 4, 6$ and $8$. Near the endpoints, the solution behaves singularly during all the snapshot times, which is the key challenge that has to be overcome by the numerical approximation for an accurate solution. The sequence of refinements is visualized in Figure 1.

<figure><img src="image_3.png" /><figcaption>FIGURE 1. Scattering at a flat screen, $d = (-\frac{\sqrt{3}}{2}, \frac{1}{2})$: sequence of meshes obtained by adaptive refinement. The red lines denote the marked elements at each iteration.</figcaption></figure>

**6.2. Wedge.** We solve the problem at a wedge geometry $\Gamma := \Gamma_0 \cup \Gamma_1$, where $\Gamma_0 := (0, 1) \times \{0\}$ and $\Gamma_1 := \{0\} \times (0, 1)$. Results for convergence of point evaluation and convergence of densities in energy norm (19) are presented in Figure 4. As before, we observe that our adaptive refinement performs as the optimal 3-graded mesh refinement.

**6.3. Trapping Geometry.** We solve the problem at a screen with the trapping geometry

$$ \Gamma := \Gamma_0 \cup \Gamma_1 \cup \Gamma_2, $$

where each of the segments is defined by

$$ \Gamma_0 := \{0\} \times (-1, 1), $$

$$ \Gamma_1 := \{x \in \mathbb{R}^2 : x = \lambda x_{10} + (1 - \lambda)x_{11}, \lambda \in [0, 1]\}, $$

$$ \Gamma_2 := \{x \in \mathbb{R}^2 : x = \lambda x_{20} + (1 - \lambda)x_{21}, \lambda \in [0, 1]\}, $$

with $x_{10} = (0, 1)$, $x_{11} = (\frac{\sqrt{3}}{2}, \frac{1}{2})$, $x_{20} = (0, -1)$, and $x_{21} = (\frac{\sqrt{3}}{2}, -\frac{1}{2})$. The geometry is visualized in Figure 6.
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Scattering at a flat screen, $d = (-\frac{\sqrt{3}}{2}, \frac{1}{2})$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 3. Snapshots of the density $\varphi_h^T$: scattering at a flat screen.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIGURE 4. Scattering at a wedge, $d = (-\frac{\sqrt{3}}{2}, \frac{1}{2})$.</figcaption></figure>

Results for convergence of point evaluation and convergence of densities in energy norm (19) are presented in Figure 5. The additional complexity of the solution induced by the trapping does not seem to affect the performance of the scheme, as the adaptive method still converges with the optimal rate. Snapshots of the solution at times $t = 1, 4, 7$ and $10$ are visualized in Figure 6.

<figure><img src="image_7.png" /><figcaption>FIGURE 5. Scattering at a screen with trapping geometry, $d = (-\frac{\sqrt{3}}{2}, \frac{1}{2})$.</figcaption></figure>

### 6.4. Time discretization error.
In this section we study the convergence of the time-discretization. We consider the setting of a flat screen from Section 6.1 refining with the adaptive algorithm in space. The results are shown in Figure 7. As a benchmark we use the solution with $\tau = \frac{T}{N} = \frac{10}{6400} = 0.0015625$.
<figure><img src="image_8.png" /><figcaption>FIGURE 6. Snapshots of the solution $u_h^\tau$: scattering at a screen with trapping geometry.</figcaption></figure>

We observe that the standard formulation ($\eta \equiv 0$) leads to a reduced convergence rate of 2 for the density $\varphi_h^\tau$, due to the results from Lemma 2. The modified formulation from (16) with $\eta = \frac{1}{2T}$ shows the classical order of convergence $2m - 1 = 3$ for the two-stage Radau IIIA method, as Theorem 1 predicts. For both cases, point evaluations show the full classical order of convergence.

Finally, we investigate the efficiency of the a posteriori error estimate, i.e. the ratio error/estimator in Figure 8. The error (and estimator) for $M$ spatial degrees of freedom is written as $e_M(\eta_M)$. The computed ratios $e_M/(\eta_M + C\tau^2)$ remain bounded, both as $M \to \infty$ and the $\tau \to 0$. We conclude therefore that the efficiency constant remains bounded, even though this cannot be shown for multistage methods, with the current theory. For the constant in the ratio that estimates the time-discretization error, we used $C = 5$.

Similar results are obtained for the modified formulation with a shift of $\frac{1}{2T}$. The results are shown in Figure 9. Here, we empirically observe a convergence rate of 2.8, close to the classical convergence rate of 3. We plot $e_M/(\eta_M + C\tau^{2.8})$ with
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Convergence in time of the standard and modified formulations.</figcaption></figure>

a constant $C = 9$, and observe how most of the values remain around 1. It is important to mention that the mesh obtained by the adaptive algorithm contains elements of very small size towards the endpoints of the segment. This leads to ill-conditioned matrices during the computation of the reference solution. Further experiments require the use of preconditioning, and a higher number of quadrature points.

<figure><img src="image_10.png" /><figcaption>FIGURE 8. Error and efficiency for different time steps. Standard formulation.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIGURE 9. Error and efficiency for different time steps. Modified formulation.</figcaption></figure>

ACKNOWLEDGEMENTS

We thank C. Schwab for initiating early discussions, which led to the inception
of this project.

REFERENCES

[1] R. A. Adams and J. J. F. Fournier. Sobolev spaces, volume 140 of Pure and Applied Mathematics (Amsterdam). Elsevier/Academic Press, Amsterdam, second edition, 2003.

[2] A. Aimi, G. Di Credico, H. Gimperlein, and C. Guardasoni. Adaptive time-domain boundary element methods for the wave equation with Neumann boundary conditions. *Computers & Mathematics with Applications*, 198:196–213, 2025.

[3] A. Aimi, G. Di Credico, H. Gimperlein, and C. Guardasoni. A space-time adaptive boundary element method for the wave equation. 2025.

[4] A. Aimi, G. Di Credico, H. Gimperlein, and E. P. Stephan. Higher-order time domain boundary elements for elastodynamics: graded meshes and hp versions. *Numerische Mathematik*, 154(1):35–101, 2023.

[5] A. Bamberger and T. Ha Duong. Formulation variationnelle espace-temps pour le calcul par potentiel retardé de la diffraction d'une onde acoustique. I. *Math. Methods Appl. Sci.*, 8(3):405–435, 1986.

[6] L. Banjai and M. Ferrari. Runge-Kutta convolution quadrature based on Gauss methods. *Numer. Math.*, 156(5):1719–1750, 2024.

[7] L. Banjai and M. Kachanovska. Fast convolution quadrature for the wave equation in three dimensions. *J. Comput. Phys.*, 279:103–126, 2014.

[8] L. Banjai, M. López-Fernández, and A. Schädle. Fast and oblivious algorithms for dissipative and two-dimensional wave equations. *SIAM J. Numer. Anal.*, 55(2):621–639, 2017.

[9] L. Banjai and C. Lubich. Runge-Kutta convolution coercivity and its use for time-dependent boundary integral equations. *IMA J. Numer. Anal.*, 39(3):1134–1157, 2019.

[10] L. Banjai, C. Lubich, and J. M. Melenk. Runge-Kutta convolution quadrature for operators arising in wave propagation. *Numer. Math.*, 119(1):1–20, 2011.

[11] L. Banjai, C. Lubich, and J. Nick. Time-dependent acoustic scattering from generalized impedance boundary conditions via boundary elements and convolution quadrature. *IMA J. Numer. Anal.*, 150(1):1–26, 2022.
[12] L. Banjai and S. Sauter. Rapid solution of the wave equation in unbounded domains. *SIAM J. Numer. Anal.*, 47(1):227–249, 2008/09.

[13] P. Bansal, A. Moiola, I. Perugia, and C. Schwab. Space-time discontinuous galerkin approximation of acoustic waves with point singularities. *IMA Journal of Numerical Analysis*, 41(3):2056–2109, 2021.

[14] A. Bespalov, T. Betcke, A. Haberl, and D. Praetorius. Adaptive BEM with optimal convergence rates for the Helmholtz equation. *Computer Methods in Applied Mechanics and Engineering*, 346:260–287, 2019.

[15] S. Börm, M. Lopez-Fernandez, and S. A. Sauter. Variable order, directional $H^2$-matrices for Helmholtz problems with complex frequency. *IMA J. Numer. Anal.*, 41(4):2896–2935, 2021.

[16] T. Chaumont-Frelet. Asymptotically constant-free and polynomial-degree-robust a posteriori estimates for space discretizations of the wave equation. *SIAM Journal on Scientific Computing*, 45:A1591–A1620, 2023.

[17] T. Chaumont-Frelet and A. 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*, 59(4):1937–1972, 2025.

[18] L. Desiderio and S. Falletta. Efficient solution of two-dimensional wave propagation problems by CQ-wavelet BEM: algorithm and applications. *SIAM J. Sci. Comput.*, 42(4):B894–B920, 2020.

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

[20] H. Gimperlein, F. Meyer, C. Özdemir, D. Stark, and E. P. Stephan. Boundary elements with mesh refinements for the wave equation. *Numerische Mathematik*, 139:867–912, 2018.

[21] H. Gimperlein, C. Özdemir, D. Stark, and E. P. Stephan. hp-version time domain boundary elements for the wave equation on quasi-uniform meshes. *Computer Methods in Applied Mechanics and Engineering*, 356:145–174, 2019.

[22] H. Gimperlein, C. Özdemir, D. Stark, and E. P. Stephan. A residual a posteriori estimate for the time-domain boundary element method. *Numerische Mathematik*, 146:239–280, 2020.

[23] M. Glaefke. *Adaptive Methods for Time Domain Boundary Integral Equations*. PhD thesis, Brunel University, 2012.

[24] J. Gwinner and E. P. Stephan. Advanced boundary element methods – treatment of boundary value, transmission and contact problems. volume 52 of Springer Series in Computational Mathematics. 2018.

[25] E. Hairer and G. Wanner. *Solving ordinary differential equations. II*, volume 14 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, revised edition, 2010.

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

[27] D. Korikov, B. Plamenevskii, and O. Sarafanov. *Asymptotic theory of dynamic boundary value problems in irregular domains*, volume 284. Springer Nature, 2021.

[28] M. Lopez-Fernandez and S. Sauter. Generalized convolution quadrature with variable time stepping. *IMA J. Numer. Anal.*, 33(4):1156–1175, 2013.

[29] M. Lopez-Fernandez and S. Sauter. Generalized convolution quadrature with variable time stepping. Part II: Algorithm and numerical results. *Appl. Numer. Math.*, 94:88–105, 2015.

[30] M. Lopez-Fernandez and S. Sauter. Generalized convolution quadrature based on Runge-Kutta methods. *Numer. Math.*, 133(4):743–779, 2016.

[31] C. Lubich. On the multistep time discretization of linear initial-boundary value problems and their boundary integral equations. *Numer. Math.*, 67(3):365–389, 1994.

[32] V. T. Luong and N. T. Tung. The Dirichlet–Cauchy problem for nonlinear hyperbolic equations in a domain with edges. *Nonlinear Analysis*, 125:457–467, 2015.

[33] J. M. Melenk and J. Nick. Parsimonious convolution quadrature. *arXiv preprint arXiv:2410.15079*, 2024.

[34] F. Müller and C. Schwab. Finite elements with mesh refinement for elastic wave propagation in polygons. *Mathematical Methods in the Applied Sciences*, 39(17):5027–5042, 2016.

[35] F. L. Müller and C. Schwab. Finite elements with mesh refinement for wave equations in polygons. *Journal of Computational and Applied Mathematics*, 283:163–181, 2015.

[36] J. Nick. *Wave scattering from nontrivial boundary conditions*. PhD thesis, Universität Tübingen, 2023.
[37] A. Rieder. A p-version of convolution quadrature in wave propagation. *SIAM Journal on Numerical Analysis*, 63(4):1729–1756, 2025.

[38] S. Sauter and M. Schanz. Convolution quadrature for the wave equation with impedance boundary conditions. *J. Comp. Phys.*, 334:442–459, 2017.

[39] S. Sauter and A. Veit. Adaptive time discretization for retarded potentials. *Numer. Math.*, 132(3):569–595, 2016.

[40] M. Zank. *Inf-sup stable space-time methods for time-dependent partial differential equations*. Verlag der Technischen Universität Graz, 2020.

[41] M. Zank and O. Steinbach. Adaptive space-time boundary element methods for the wave equation. *PAMM*, 16:777–778, 2016.

INRIA UNIV. LILLE AND LABORATOIRE PAUL PAINLEVÉ, 59655 VILLENEUVE-D'ASCQ, FRANCE.

Email address: theophile.chaumont@inria.fr

ENGINEERING MATHEMATICS, UNIVERSITY OF INNSBRUCK, 6020 INNSBRUCK, AUSTRIA.

Email address: heiko.gimperlein@uibk.ac.at

INSTITUTE FOR MATHEMATICAL AND COMPUTATIONAL ENGINEERING, SCHOOL OF ENGINEERING AND FACULTY OF MATHEMATICS, PONTIFICIA UNIVERSIDAD CATÓLICA DE CHILE, SANTIAGO, CHILE, AND

ENGINEERING MATHEMATICS, UNIVERSITY OF INNSBRUCK, 6020 INNSBRUCK, AUSTRIA.

Email address: ignacio.labarca@uc.cl

ETH ZÜRICH, SEMINAR FOR APPLIED MATHEMATICS, RÄMISTRASSE 101, CH-8092 ZÜRICH, SWITZERLAND.

Email address: joerg.nick@math.ethz.ch