On superconvergence of Runge-Kutta convolution quadrature for the wave equation

The semidiscretization of a sound soft scattering problem modelled by the wave equation is analyzed. The spatial treatment is done by integral equation methods. Two temporal discretizations based on Runge-Kutta convolution quadrature are compared: one relies on the incoming wave as input data and the other one is based on its temporal derivative. The convergence rate of the latter is shown to be higher than previously established in the literature. Numerical results indicate sharpness of the analysis. The proof hinges on a novel estimate on the Dirichlet-to-Impedance map for certain Helmholtz problems. Namely, the frequency dependence can be lowered by one power of $\abs{s}$(up to a logarithmic term for polygonal domains) compared to the Dirichlet-to-Neumann map.


Introduction
Boundary element methods have established themselves as one of the standard methods when dealing with scattering problems, especially if the domain of interest is unbounded.First introduced for stationary problems, beginning with the seminal works [BH86a,BH86b] these methods have steadily been extended to time-dependent problems; see [Say16] for an overview.The method of convolution quadrature (CQ), introduced by Lubich in [Lub88a,Lub88b], is a convenient way of extending the stationary results to a time-dependent setting.It is well-known that the convergence rate of a Runge-Kutta convolution quadrature (RK-CQ) as introduced in [LO93], is determined by bounds on the convolution symbol K in the Laplace domain.Namely, a bound of the form K(s) ≤ C |s| µ leads to convergence rate q + 1 − µ, as was proven in [BLM11], see also [BL11,LO93] for earlier results in this direction.Thus one might expect that changing the symbol to s −1 K(s) would increase the convergence order by one.When considering discretizations of the wave equation using boundary integral methods, this is not always the case.Instead, it has been observed that sometimes a "superconvergence phenomenon" appears, where the observed convergence rate surpasses those predicted, see [RSM20a,RSM20b,Rie17].In this paper, we give a first explanation why such a phenomenon occurs in the model problem of sound soft scattering, i.e., the discretization of the Dirichlet-to-Neumann map.We expect that similar phenomena can also explain the improved convergence rate for the Neumann problem or more complex scattering problems.The proof relies on the observation that the s −1 -weighted Dirichlet-to-Neumann map can be decomposed into a Dirichlet-to-Impedance map plus the identity operator.For the Dirichlet-to-Impedance operator, it was observed in [Ban14] that an improved bound holds compared to the Dirichlet-to-Neumann map as long as the geometry is given by the sphere or the half-space.It is then conjectured that a similar bound holds for smooth, convex geometries.In this paper, provided that we restrict the Laplace parameter s to a sector, we generalize this result to a much broader class of geometries, namely, smooth or polygonal geometries, without convexity assumption.This will then immediately give the stated improved bound for the convolution quadrature scattering problem.In the case of polygons, the result holds in a slightly weaker form in that it contains an additional logarithmic term.As a consequence of this observation, it may often be beneficial to select a problem formulation with an extra time derivative.In many situations, such formulations are even the natural choice, see, e.g., [BR18,BL18,BLS15], and especially when one works with the wave equation as a first order system as in [RSM20b].Another way of looking at this phenomenon is that when using a standard formulation (i.e., taking λ k as in Proposition 3.2 without using a time derivative on the data), then the discrete integral will exhibit a superconvergence effect.We point out that the present paper focuses on a semidiscretization of the problem with respect to the time variable.For practical purposes one would also have to take into account the discretization in space using boundary elements.We also mention that while popular, CQ is only one possibility to apply boundary integral techniques to wave propagation problems.Notably also space-time based methods have gained popularity [GNS17, GMO + 18, JR17] in recent years.

Model Problem and Notation
We consider a sound soft scattering problem for acoustic waves.For a bounded Lipschitz domain Ω − ⊆ R d with Ω + := R d \ Ω − , the problem reads: Find u tot such that ütot = ∆u tot in Ω + , u tot (t)| Γ = 0 for t > 0, u tot (t) = u inc (t) for t ≤ 0.
Here u inc is a given incoming wave, i.e., u inc also solves the wave equation, and we assume that for t ≤ 0 it has not reached the scatterer yet.The problem can be recast by decomposing the total wave into the incoming and outgoing wave, u tot = u inc + u, where u solves: ü = ∆u in Ω + , u(t)| Γ = −u inc (t)| Γ for t > 0, u(t) = 0 for t ≤ 0. (2.1) This will be the problem we are discretizing.
For simplicity, we consider two possible cases.The bounded Lipschitz domain Ω − ⊆ R d has either a smooth boundary or Ω − ⊆ R 2 is a polygon.While we expect that the results and techniques can be generalized to the case of piecewise smooth geometries, such an extension would lead to a much higher level of technicality in the present paper.Although we focus on the exterior scattering problem as our motivating model problem all of the main results also hold for the interior Dirichlet problem.We end the section by fixing some notation.We write H m (Ω ± ) for the usual (complex valued) Sobolev spaces on Ω + or Ω − .On the interface Γ := ∂Ω we also need fractional spaces H s (Γ) for s ∈ [−1, 1], see, e.g., [McL00,AF03] for precise definitions.We also set for the exterior and interior trace operator, and for the normal derivative.We note that in both cases, we take the normal to point out of the bounded domain Ω − .We write γu := γ + u − γ − u and { {γu} } := 1 2 (γ + u + γ − u) for the trace jump and mean, and

Boundary Integral Methods and Convolution Quadrature
It is well-known that scattering problems of the form presented in Section 2 can be solved by employing boundary integral methods, see [Say16] for a detailed time domain treatment.For the frequency domain, results can be found in most textbooks on the subject, see [SS11,Ste08,McL00,GS18,HW08].The use of boundary integral methods for discretizing the time domain scattering problem dates back to the works [BH86a,BH86b], where also important Laplace domain estimates of the form (3.2) were first shown.For s ∈ C + := z ∈ C : Re(z) > 0 , we introduce the single and double layer potentials where Φ is the fundamental solution for the operator −∆ + s 2 : Here H (1) 0 denotes the first kind Hankel function of order zero, see [McL00,Chap. 9].Finally, we introduce the boundary integral operators induced by the potentials: In practice, these operators can be realized explicitly as integrals over the boundary Γ since for sufficiently smooth functions ψ, ϕ the following equations hold: The operator we consider for discretizing (2.1) is the Dirichlet-to-Neumann map.
We then define the operators (2.6) In practice, the following well-known proposition gives an explicit way to calculate DtN.
Proposition 2.2 (see, e.g., [LS09, Appendix 2]).The Dirichlet-to-Neumann map can be written as (2.7) RK-CQ was introduced by Lubich and Ostermann in [LO93].It provides a simple and general way of approximating convolution integrals by a high order method and has the great advantage that only the Laplace transform of the convolution symbol is required.We only very briefly introduce the method and notation.Let K be a holomorphic function in the half plane Re(s) > σ 0 > 0. Let L denote the Laplace transform and L −1 its inverse.We (formally) introduce the operational calculus by defining where g ∈ dom (K(∂ t )) is such that the inverse Laplace transform exists and the expression above is well defined.For a Runge-Kutta method given by the Butcher tableau A, b T , c, the convolution quadrature approximation of K(∂ t ) with time step size k > 0 is given, for any function g : R → R with g(t) = 0 for t ≤ 0, by the expression where the matrix valued function ∆ is given by The extension to operator valued functions K is straight forward.In practice, we only consider evaluating K(∂ k t )g at the discrete time steps t j := j k.We make the following assumptions on the Runge-Kutta method, slightly stronger than [BLM11].
Assumption 2.3.(i) The Runge-Kutta method is A-stable with (classical) order p ≥ 1 and stage order q ≤ p.
(ii) The stability function R(z (iii) The Runge-Kutta coefficient matrix A is invertible.
Remark 2.4.Assumption 2.3 is satisfied by the Radau IIA and Lobatto IIIC methods, see [HW10].Also note that the order conditions imply that c m = 1 for such methods.

Main results
To simplify the notation, we introduce a symbol for the sectors in Proposition 2.5.Throughout this work we fix σ 0 > 0 and δ > 0 and set Remark 3.1.The choice of σ 0 > 0 and δ > 0 in the definition of S is arbitrary, and all our estimates will hold for any choice, although all the constants will depend on σ 0 and δ.
We are now able to state the main result of the paper.We start by stating the standard convergence result for discretizing the Dirichlet-to-Neumann map.
Let λ := DtN ± (∂ t )g be the exact normal derivative and λ k := DtN ± (∂ k t )g denote the standard CQ-approximation.Then there exist constants C(T ) > 0 and k > 0 such that the following estimate holds for 0 ≤ t ≤ T and 0 < k < k: with a constant C(T ) depending on the terminal time T , the Runge-Kutta method, Γ, and k.
Proof.Follows from the well-known bound (see for example [LS09]) and Proposition 2.5.
We will observe numerically in Section 5 that Proposition 3.2 is essentially sharp.Thus, when considering the differentiated equation, one expects an increase in order by one, which follows directly from Proposition 2.5.But for the Dirichlet-to-Neumann map the increase of order is even greater, as long as one assumes slightly higher regularity of the data.
] ġ denote the CQ-approximation using ġ as input data.Then, for all ε > 0, there exist constants C(T, ε) > 0 and k > 0 such that the following estimate holds for 0 ≤ t ≤ T and 0 < k < k: (3.3) The constant C(T, ε) depends on ε, the terminal time T , the Runge-Kutta method, Γ, and k.If Γ is smooth, one can take ε = 0.
Proof.We apply Proposition 2.5.By linearity, we can write the Dirichlet-to-Neummann operator as s −1 DtN(s) = s −1 DtI(s) + I or, in the time domain, The second operator (in frequency domain) is independent of s.It is a simple calculation that in such cases, i.e., if K(s) = T for all s, the convolution weights satisfy W j = δ j,0 T .Thus, we have Since stiff accuracy implies b T A −1 = (0, . . ., 0, 1) and c m = 1, the operator I is reproduced exactly by the CQ.A similar decomposition was already invoked in [BLM11] to explain a superconvergence phenomenon for a scalar problem.Combining standard estimates, e.g., [LS09, Table 1], with Theorem 3.4 shows that the Dirichlet-to-Impedance map satisfies By Proposition 2.5 and by estimating the logarithmic term by C|s| ε for arbitrary ε > 0, we obtain (3.3).
While Theorem 3.3 is the main motivation for this paper, its proof is based on another result, which may be of independent interest.
Theorem 3.4.Let s ∈ S .Let Ω ⊆ R d be a bounded smooth Lipschitz domain.The following estimate holds for the Dirichlet-to-Neumann map: If Ω ⊆ R 2 is a bounded Lipschitz polygon, then one has The constant C depends only on Ω and the parameters σ 0 , δ defining the sector S .
Proof.Due to its lengthy and technical nature, we defer the proof to Section 4. For smooth geometries it is shown as Corollary 4.11.Polygonal domains are handled in Corollary 4.19.
Remark 3.5.The regularity requirement g ∈ H 1 (Γ) is stronger than the expected requirement g ∈ H 1/2 (Γ).This is due to the construction of the boundary layer function (see (4.15)).
Since all our results hold for both the interior and exterior problem, we can also easily treat the case of an indirect BEM formulation.
Corollary 3.6 (Indirect formulation).Let s ∈ S and assume that Ω ⊆ R d is smooth.Then, the operator V −1 (s) − 2s satisfies the bound ] ġ be its CQ-approximation.Then there exist constants C(T ) > 0 and k > 0 such that the following estimate holds for 0 ≤ t ≤ T and 0 < k < k: (3.7) The constant C(T ) depends on T , the Runge-Kutta method, Γ, and k.
If Ω ⊆ R 2 is a polygon, then and where ε > 0 is arbitrary, and C(T, ε) depends additionally on ε.

Proofs
The proof of Theorem 3.4 hinges on three main observations, which require some technical work to be made rigorous: 1.In 1d on R + , the interior Dirichlet-to-Neumann map is given by g → sg.
2. The existing DtN-estimate's poor s dependence is mainly caused by boundary layers.
3. Boundary layers are essentially a 1d phenomenon, so observation 1 applies.

Preliminaries
There are many ways of defining fractional order Sobolev spaces.A convenient way of working with them is by introducing them via the real method of interpolation.Given Banach spaces X 1 ⊆ X 0 with continuous embedding and parameters θ ∈ (0, 1), q ∈ [1, ∞), we define the interpolation norm and space as follows: For θ ∈ (0, 1), the corresponding norms are defined via interpolation.If we want to include boundary conditions, we write • |s|,θ,∼,O for the interpolation space between L 2 (O) and H 1 0 (O), but equipped with the weighted norms (4.2).The dual norms are defined by θ,∼,O .
For the most part we will be working with the closed surface ∂Ω.There, the norms • |s|,θ,∂Ω and • |s|,θ,∼,∂Ω coincide.By Lemma A.1 we also have for θ ∈ (0, 1) and bounded domains O the norm equivalence with implied constants independent of |s|.
We start with some well-known s-explicit estimates for the (modified) Helmholtz equation.
Lemma 4.2 (Well posedness).Let s ∈ S .The sesquilinear form Proof.We calculate: Since Re(s) ∼ |s| in the sector S this concludes the proof.
Lemma 4.3 (Trace estimates).For Re(s) ≥ 0 and |s| > σ 0 , let u ∈ H 1 (Ω ± ) satisfy Then the following estimates hold for the traces of u: Proof.We start with the normal derivative.For any ξ ∈ H 1/2 (Γ) and any v with γ ± v = ξ we calculate: where in the last step we select v as the minimal energy function, satisfying , and (4.4) follows.
For the Dirichlet trace, we get using the multiplicative trace estimate and the same lifting v: The estimate for the impedance trace then follows trivially.
The previous lemma shows that for a priori estimates in terms of standard Sobolev norms the constants involved have some s dependence.The next lemmas show that the use of the weighted norms introduced in Definition 4.1 avoids such dependencies: Lemma 4.4.The operators γ ± : H 1 (Ω ± ) → H 1/2 (Γ) satisfy the bounds: Proof.The multiplicative trace estimate and Young's inequality give Combining this with the standard trace estimate concludes the proof in view of (4.3).
Lemma 4.5 (Dirichlet problem).Let g ∈ H 1/2 (Γ), f ∈ L 2 (Ω ± ).For any s ∈ S there exists a unique solution to the problem The function satisfies the a priori bound The implied constant depends only on Ω ± and the constants σ 0 , δ characterizing S .
Proof.Existence follows using the usual theory of elliptic problems.For the a priori bound, we first note that by [MS11, Lemma 4.22], there exists a lifting u D satisfying Thus the remainder u := u − u D solves: As the sesquilinear form a s from Lemma 4.2 is elliptic, we get with ζ defined there . Then for every s ∈ S there exists a unique solution to the problem u satisfies the a priori bound (4.10) The implied constant depends only on Ω ± and on σ 0 , δ characterizing S .
Proof.Follows easily from the weak formulation and (4.8).
We also have the following trace inequality in a weighted H −1/2 -norm: Lemma 4.7.If −∆u + s 2 u = 0 we can estimate: Proof.Follows easily from the weak definition of ∂ − n u, the Cauchy-Schwarz inequality, and (4.8).

Smooth geometries
In order to prove a first version of Theorem 3.4, we consider a simplified setting of smooth geometry and Dirichlet trace.Closely following the ideas from [MS99, Mel02], we construct a lowest order boundary layer function that will be the basis for all further estimates.
where n( x) is the outer normal vector to Ω − at the point T ( x).
For ε > 0 sufficiently small, F is a smooth diffeomorphism onto where T and R are smooth and T ( x) is invertible.
Here T 1 := (e 1 , . . ., e d−1 ) T D x T ( x), and thus T 1 2 ≤ D x T ( x) 2 .We further compute: where R 1 collects the remaining terms.For sufficiently small ρ > 0, depending only on D x T 2 and D x n 2 , we can linearize the inverse in (4.13) to get (4.12) with T := T 1 T T 1 −1 (the latter inverse exists since D x T and thus also T 1 has full rank).
(iv) For ε > 0 define the set Then, the following estimates hold for all ∈ R with constants independent of s: (v) The analogous statement also holds for the exterior problem Ω + , replacing −s by s in (ii).
Proof.We only show the case of the interior problem and abbreviate g := γ − u.We work in boundary fitted coordinates ( x, ρ) as described in Lemma 4.8.First assume, that supp(g) ⊂ T (O), i.e., g is supported by the part of the boundary parametrized by T .The change of variables formula shows that if u solves −∆u + s 2 u = f , then u := u • F solves: with J := det(DF ) and f = f • F (see, e.g., [Eva98, Step 7 of proof of Thm. 4, Sec.6.3.2]).On the other hand, if u BL satisfies We set A := DF −1 DF −T and define with g := g • T the function where, in the last step, we exploited the definition of u BL .From its definition, one can easily see that u BL satisfies the estimates Re(s) , Re(s) .
Transforming back gives (4.14) for the part of Ω − parametrized by F .Assertion (iv) follows easily from the definition, as the exponential decay dominates all powers of |s|.This allows us to smoothly cut off u BL for large ρ and extend it by 0 to the whole domain.For general g, we use a smooth partition of unity to decompose g into functions with local support.
As the next step, we lower the regularity requirement on γ − u.
Corollary 4.10.Let Ω − have a smooth boundary Γ.For any s ∈ S and for every u ∈ H 1 (Ω − ) with γ − u ∈ H 1 (Γ) solving −∆u + s 2 u = 0 there exists a function u BL ∈ H 2 (Ω − ) with the following properties: The implied constant depends only on Ω − and the constants σ 0 , δ characterizing S . (iv Then, the following estimates hold for all ∈ R with constants independent of s: (v) The analogous statement also holds in the case of the exterior problem upon replacing s by −s in (i) and (ii).
Proof.In order to apply Lemma 4.9, we need H 2 -regularity of g := γ − u.We fix a function g ∈ H 2 (Γ) with the following properties: This can be either seen by realizing H 1 (Γ) as the interpolation space between L 2 (Γ) and H 2 (Γ) and using [BS78,Lemma] or constructed directly via the usual mollifiers as done in [AF03, Thm.2.29]: The approximation estimate follows from [AF03,Eqn (20)] and an interpolation argument.See also [AF03,Sec. 7.48] for how to trade Sobolev regularity for approximation properties of the mollified function.
Let u denote the solution to Since g ∈ H 2 (Γ), we can apply Lemma 4.9 to construct u BL .Assertion (i) then follows by construction.For (iii): we note that by Lemmas 4.5 and 4.9: For (ii), we use Lemma 4.3 and (4.9) to get that Similarly, we have Assertion (iv) follows directly from Lemma 4.9 (iv) and (4.16).
Corollary 4.11.Let Ω − ⊂ R d be smooth and s ∈ S .Let g ∈ H 1 (Γ) and u solve Then (4.17) The analogous statement holds for the exterior problem upon replacing s by −s in (4.17).
Proof.Follows by writing u = u BL + (u − u BL ).The impedance trace of u BL vanishes by Corollary 4.10 (i).The impedance trace of the remainder is uniformly bounded with respect to s via Corollary 4.10 (ii).

Polygons
In this section, we consider a polygonal domain Ω ⊂ R 2 as an example of a non-smooth domain.
In order to match the boundary layer solutions from Lemma 4.9 at corners, we solve an appropriate transmission problem, similarly to what was done in [Mel02].We refer to Figure 4.1b for the geometric situation.We first need one additional Sobolev space.For a smooth curve Γ and θ ∈ [0, 1], we introduce where d ∂Γ denotes the distance to the endpoints of Γ .

A transmission problem in a cone
In this section, we investigate certain transmission problems.These will allow us to match different boundary layer functions in the vicinity of a corner of the domain.We start by investigating the special case of a transmission problem on a sector or an infinite cone.Due to its special structure, we can derive sharper estimates for the normal derivative than what can be obtained from the energy methods used in Lemma 4.16 below.We introduce some notation.Given ω ∈ (0, π), we define the infinite cone with opening angle 2ω and C by removing from C its bisector: Next, we define the sector S ω := {(r cos ϕ, r sin ϕ), r ∈ (0, 1), |ϕ| ∈ (0, ω)}, which is just the truncated cone C ∩ B 1 (0).For its boundary, we write Γ ±ω := (r cos(±ω), r sin(±ω)), r ∈ (0, 1) for the two parts of the boundary of the sector that are adjacent to the origin and set Γ S := Γ ω ∪ Γ −ω .Finally, we need to define the normal jump across interfaces.If Γ denotes a smooth interface separating domains O 1 and O 2 we define the normal jump across Γ via where the normal vectors n j are taken to point out of O j respectively.
Lemma 4.12.Consider the solution u ∈ H 1 (C) to the following problem on the infinite cone for µ > 0 and s ∈ S with | s| = 1: Then, the following statements hold for u: (i) For each ∈ N there exist constants C , α > 0 such that for all r ≥ 1 u W ,∞ (C \Br(0)) ≤ C e −α r .
(ii) There exists a constant C > 0 such that ∂ n u satisfies the estimates The constants depend only on the opening angle 2ω, the parameter µ, and the choices of σ 0 and δ in the definition of S .
Proof.We first show (ii) in Steps 1-3 and then (i) in Step 4.
Step 1: We start with an energy estimate in exponentially weighted spaces, namely, for any 0 < α < µ Re( s) the following estimate holds: with a constant C only depending on α, µ, ω, and s.We fix some notation.We write u 2 1,α := e αr ∇u 2 L 2 (C) + e αr u 2 L 2 (C) , and analogously for u 1,−α .Also we set h(x 1 , x 2 ) := e − sµx 1 for the transmission data.The proof follows [Mel02, Prop.6.4.6] verbatim.The sesquilinear form B(u, v) := (∇u, ∇u) L 2 (C) + ŝ2 (u, v) L 2 (C) satisfies an inf-sup condition: There is c > 0 depending only on α ∈ [0, 1) such that inf This can be seen by taking, for given u in the infimum, the function v := se 2r u in the supremum and performing elementary calculations.Next, we show that |(h, γv)  We conclude that the solution u satisfies u 1,α ≤ C for some constant C > 0 depending only on the choice of α < µ Re( s).
Step 5: Inspection of the proof reveals that all constants (if at all) depend continuously on ŝ.Since we are only interested in ŝ in a compact set determined by the constants from the definition of S we can make all the constants independent of ŝ.
Having studied the transmission problem in a dimensionless form in Lemma 4.12, we can transfer the results to the setting we actually require using a scaling argument.
Lemma 4.13.Fix ω ∈ (0, π).For s ∈ S and µ > 0, let u ∈ H 1 (S ω ) solve the transmission problem on the sector S ω : Recall that Γ S denotes the parts of ∂S ω adjacent to the origin.Then The constants C > 0 depend only on ω, the parameter µ, and the choices of σ 0 and δ in the definition of S .
We need the following modification of [MPW17, Lemma 3.13].
Proposition 4.14.Let O be a Lipschitz domain.Define the Besov space (cf.(4.1)) For ε > 0 and every w ∈ H 1/2 (O), there exists a function The constant depends only on the domain O.
Proof.This is essentially [MPW17, Lemma 3.13].The only modification needed is that we consider the H 1/2 (O)-norm on the right-hand side instead of the B 1/2 2,∞ -norm, which is the reason for getting penalized by a factor |log(ε)| instead of a factor |log(ε)| as in [MPW17, Lemma 3.13].The result follows from the same proof, only noting that one can bound using the Cauchy-Schwarz inequality , and the last factor produces a factor 1 + | log(ε)|.
Lemma 4.15.Let u solve (4.25).Then, there exists a constant C > 0 depending only on ω, µ and the parameters in the definition of S such that Proof.For w ∈ H 1/2 (Γ), select w ε as in Proposition 4.14 with ε > 0 to be chosen later.We calculate for Γ ±ω , i.e., the two parts of ∂Ω adjacent to the origin: Since Γ ±ω is a one-dimensional line segment, we can use the Sobolev embedding [Tar07, Sec.32] to estimate . Overall, we get using the properties of w ε from Proposition 4.14 and the estimates on ∂ n u from (4.24): Choosing ε = |s| −1 completes the proof.
We are now in a position to study a more general transmission problem, namely, allowing for Dirichlet jumps and more general Neumann transmission data.Given g ∈ H 1/2 (Γ ), h ∈ H −1/2 (Γ ), there exists a unique solution u ∈ H 1 (O 1 ∪ O 2 ) to the following problem: Additionally, the following estimate holds: If O 1 and O 2 are polygons, Γ is a straight line, and h can be decomposed as for some h 1 ∈ R, µ > 0, and h 2 ∈ H −1/2 (Γ ), then ∂ n u exists pointwise almost everywhere and Proof.Proof of (4.27):Since g is assumed in H 1/2 (Γ ), we can extend it by 0 to a function g ∈ H 1/2 (∂O 1 ) such that (see for example [McL00,Thm. 3.33]) . We solve a Dirichlet problem on O 1 with data g to obtain u 1 and extend it by 0 to O 2 .Then we solve the following problem on O: , where γ Γ denotes the trace operator on Γ .The function u := u 1 + u 2 then solves the transmission problem.The estimate (4.27) follows from Lemmas 4.5 to bound u 1 and, in order to bound u 2 , the ellipticity of the sesquilinear form (see Lemma 4.2) together with the trace estimate (Lemma 4.4) to estimate the contribution h, γ Γ v Γ .Proof of (4.28) for h 1 = 0: Introduce V := ∂Γ ∪ {vertices of Ω}.Since O 1 and O 2 are piecewise smooth and u = 0 on ∂O and the right-hand side is homogeneous, the solution u is smooth up to the boundary with the exception of the vertices of Ω and near the interface Γ .Hence, ∂ n u exists pointwise everywhere on ∂O \ V. To show the estimate (4.28), we consider test functions v ∈ V := {v ∈ C ∞ (O) | v vanishes in a neighborhood of V}.Since ∂ n u exists pointwise on ∂O \ V and v ∈ V vanishes in a neighborhood of V, the duality pairing ∂ n u, v ∂O is well-defined and an integration by parts gives Since V is dense in H 1 (Ω) (because V consists of finitely many points), the equation (4.29) actually holds for all v ∈ H 1 (Ω).Given ξ ∈ H 1/2 (∂O) we select v ξ ∈ H 1 (Ω) with v| ∂O = ξ as the lifting given by (4.7) in Lemma 4.3, which satisfies By the trace theorem we have v ξ |s|,1/2,Γ v ξ |s|,1,O .Taking the supremum over all ξ ∈ H 1/2 (∂O) yields (4.28) for the case h 1 = 0. Proof of (4.28) for h 1 = 0: If h 1 = 0, we lift this contribution separately by solving the corresponding transmission problem (4.25).The estimate follows by applying Lemma 4.15 and a suitable localization.

Corner layers
Before we can bound the functions used to match boundary layers, we must control the jump between two boundary layer solutions.We start with a very simple geometric situation.
The first term is handled analogously to the L 2 -term.For the second term we use the crude estimate |e −sµrc 2 | 1 and get: Interpolating (4.32) and (4.33) then gives (4.31).

Decomposing the DtN-operator
In Section 4.2 we discussed the DtN-operator for smooth geometries.Here, we study the case of polygonal domains.We will do so by introducing corner layers, similarly to what was done in [Mel02, Sec.7.4.3].The following Theorem 4.18 presents a decomposition of the DtN-Operator into several contributions.To describe them, we need some notation as illustrated in Fig. 4.1b.The polygon Ω has vertices A 1 , . . ., A J and edge Γ j connects A j with A j+1 (we set A J+1 := A 1 and Γ J+1 := Γ 1 and, for simplicity of notation, we assume that ∂Ω consists of a single component of connectedness).Γ j is the bisector of the angle at vertex A j .The subdomains Ω j are confined by four curves: Γ j , the bisectors at A j and A j+1 (dashed black in Figure 4.1a), and a fourth curve completely contained in Ω and sufficiently close to Γ j (dashed blue in Figure 4.1a) .We set Ω 0 := Ω\∪ J i=1 Ω i and, for convenience Ω J+1 := Ω 1 .We fix χ BL ∈ C ∞ (R 2 ) with supp χ BL ⊂ ∪ i Ω i and χ BL ≡ 1 near Γ.Finally, for each vertex j we let χ CL,j ∈ C ∞ (R 2 ) be a cut-off function with supp χ CL,j ∩ {A j } = ∅ for j = j such that χ CL,j ≡ 1 on Γ Then u can be decomposed as u = χ BL u BL + J j=1 χ CL,j u CL,j + r such that for a C > 0 depending only on Ω: (ii) For each j ∈ {1, . . ., J} the function u CL,j is in H 1 (Ω i ) for each i = 0, . . ., J and u CL,j | Γ = 0. Furthermore, −∆u CL,j + s 2 u CL,j = 0 on Ω j ∪ Ω j+1 and ∂ n χ CL,j u CL,j exists on each edge Γ i , i = 1, . . ., J, and (iii) The remainder r satisfies The analogous statement holds for the exterior problem upon replacing s by −s in (i)-(iii).
Proof. 1. step (construction of u BL ): For each Ω i , let (θ i , ρ i ) be the boundary fitted coordinates obtained by an affine parametrization of the line that contains Γ i by θ i and denoting by ρ i the (signed) distance from that line.Write g(θ i ) for the function g on Γ i in the coordinates (θ i , ρ i ) and extend it H 1 and H 2 -stable to the line.We define, in boundary fitted coordinates (θ i , ρ i ), the function u BL (θ i , ρ i ) := g(θ i )e −sρ i .That is, the function u BL is given by applying the construction from Lemma 4.9.We have by construction ∂ n u BL − su BL = 0 on Γ and u BL |s|,1,Ω i |s| g H 1 (Γ i ) . 2. step (construction of u CL ): The function u BL is discontinuous across the bisectors Γ j .The corner layers u CL,j corrects this.Focussing on the bisector Γ j , let Γ j and Γ j+1 be the edges meeting at A j .Fix κ such that χ BL u BL ≡ 0 on Γ ∩ R 2 \ B κ (A j ).On the sector S ω = B κ (A j ) ∩ Ω − define u CL,j as the solution of the following transmission problem: Up to translation and rotation, we are essentially in the setting of Lemmas 4.16 and 4.17.That is, on Ω j and Ω j+1 the function u BL has the form given in Lemma 4.17 so that (taking additionally the effect of χ BL into account) we arrive at where r j = dist(•, A j ).For the normal derivative, Lemma 4.17 provides the representation By Lemma 4.16, we therefore get Furthermore, Lemma 4.16 provides for ∂ n (χ CL,j u CL,j ) on (Γ j ∪ Γ j+1 ) ∩ B κ (A j ) the bound Noting that u CL,j = 0 on Γ, the asserting of (ii) follows.

3.
Step (Construction of r): The function r ∈ H 1 0 (Ω − ) is defined as r := u−u BL − J j=1 χ CL,j u CL,j .It satisfies the equation The bounds of Lemma 4.3 and 4.5 then conclude the proof.

Step (exterior domains):
The result for the exterior problem follows along the same lines.
The analogous statement holds for the exterior problem upon replacing s by −s.
Proof.We employ the smoothing technique as in Corollary 4.10: By smoothing, on a length scale |s| ≥ 1 or by interpolation, one can construct a function g ∈ H 1 (Γ) such that g ∈ H 2 (Γ i ) for each edge Γ i and such that This can be done in two steps: first, one defines the approximation edgewise and in a second step ensure continuity at the vertices of Ω by introducing an appropriate correction, e.g., by a piecewise linear function.The remainder of the proof is then as in Corollary 4.10.

Numerical Examples
In this section, we compare the performance of the numerical schemes of Theorem 3.3 with the more standard method of Proposition 3.2 for an interior scattering problem.That is, we compare the Runge-Kutta convolution quadrature approximation by the following two methods: • [DtN − (∂ k t )]u inc , which is denoted "standard method", and ] uinc , which is denoted "differentiated method".
We use two different Runge-Kutta methods of the Radau IIA family, one with 3 and one with 5 stages.For the 3-stage version, we have q = 3 and p = 5.We therefore expect a convergence rate of order 3 for the standard method and full classical order 5 up to logarithmic terms for the differentiated scheme.
In order to show that our theoretical estimates are sharp, we also look at the 5-stage method.There, the stage order is q = 5 and the classical order p = 9.The expected rates are therefore 5 and 7 respectively for the two numerical schemes up to logarithmic terms.
As the space discretization, we employ a Galerkin boundary element method of order 5, based on a code developed by F.-J. Sayas and his group at the University of Delaware.A sufficiently refined grid is employed to be able to focus on the temporal error.Instead of evaluating the H −1/2 -error, we compute the quantity max j=0,...n V (1)e j , e j Γ with e j := Π L 2 λ(t j ) − λ k (t j ).
Here Π L 2 denotes the L 2 -orthogonal projection onto the BEM space.Since the grid is sufficiently fine and fixed, this should not impact the observed convergence rates.The operator V (1) was taken because it gives an (s-independent) equivalent norm on H −1/2 (Γ).In Figure 5.1, we observe that the rates from Proposition 3.2 and Theorem 3.3 are obtained as predicted.We conclude that while the fact that the rate jumps by order 2, even though the modification of the scheme is of order one, is at first surprising, this can be rigorously explained by Theorem 3.3.Observations of this type provided the main motivation for the investigations in this work.
the Poincaré inequality: Conversely, for the jump of the normal derivative.The notation A B abbreviates A ≤ CB with implied constants independent of critical parameters, in particular the parameter s that appears throughout this work.For (relatively) open sets O we introduce the L 2 (O) scalar product (u, v) L 2 (O) := O uv.
.1b) When working with the Helmholtz equation, it is convenient to work with |s|-weighted norms: Definition 4.1.For an open (or relatively open) set O, parameters s ∈ C + and θ ∈ {0, 1}, we define the weighted Sobolev norms

Lemma 4. 16 (
Transmission problem).Let O ⊂ R 2 be an open Lipschitz domain.Let Γ ⊂ O be a smooth interface that splits O into two disjoint Lipschitz domains O 1 and O 2 .
Figure 5.1.:Comparison of the standard and differentiated method for different RK schemes
we have by the mapping properties of the single layer potential (see [CWGLS12, Thm.2.17]) that j| ∂(C∩B 2 Subdomains for defining the corner layers u CL and boundary layers u BL : Solid blue regions indicate supp(u CL ).The regions Ω j on which u BL is defined are confined by the dashed lines.Situation at a corner A j .Marked in color is the support of the cut-off functions χ BL ; solid green: Ω j , dotted red: Ω j+1 .