Erik Burman

Guillaume Delay

Alexandre Ern

A hybridized high-order method for unique continuation subject to the Helmholtz equation

We design and analyze an arbitrary-order hybridized discontinuous Galerkin method to approximate the unique continuation problem subject to the Helmholtz equation. The method is analyzed using conditional stability estimates for the continuous problem, leading to error estimates in norms over interior subdomains of the computational domain. The convergence order reflects the Hölder continuity of the conditional stability estimates and the approximation properties of the finite element space for sufficiently smooth solutions. Under a certain convexity condition, the constant in the estimates is independent of the frequency. Moreover, certain weighted averages of the error are shown to converge independently of the stability properties of the continuous problem. Numerical examples illustrate the performances of the method with respect to the degree of ill-posedness of the problem, increasing polynomial order, and perturbations in the data.

Introduction

We are interested in the computational approximation of solutions to linear ill-posed problems based on the Helmholtz differential operator for acoustic wave propagation in the frequency domain. Ill-posed problems typically arise in the context of inverse problems and data assimilation, and have traditionally been approached through regularization at the continuous level, using either Tikhonov regularization [START_REF] Tikhonov | Solutions of ill-posed problems[END_REF] or quasi-reversibility [START_REF] Lattès | Méthode de quasi-réversibilité et applications[END_REF]. The regularized problem is well-posed and can be discretized by using any suitable numerical method. Since the regularization typically relies on a variational framework, this approach fits the "optimize then discretize" paradigm. The size of the regularization parameter is chosen depending on the noise in the data, and the mesh-size should then be chosen small with respect to the regularization parameter so that the discretization error is negligible. There is a very large literature on the topic of Tikhonov regularization and inverse problems, and we refer the reader to [START_REF] Ito | Inverse problems[END_REF] and the references therein for an overview on computational methods using this approach. Methods based on quasi-reversibility that are related to the present work can be found in [START_REF] Bourgeois | A mixed formulation of quasi-reversibility to solve the Cauchy problem for Laplace's equation[END_REF][START_REF] Bourgeois | Convergence rates for the quasi-reversibility method to solve the Cauchy problem for Laplace's equation[END_REF][START_REF] Dardé | An H div -based mixed quasi-reversibility method for solving elliptic Cauchy problems[END_REF].

An alternative route to the design of computational methods for ill-posed problems was proposed in [START_REF] Burman | Stabilized finite element methods for nonsymmetric, noncoercive, and ill-posed problems. Part I: Elliptic equations[END_REF] based on the "discretize then optimize" paradigm. The idea is to discretize the problem first, e.g., using the finite element method, and then to choose regularization terms known from the finite element approximation of well-posed problems with poor stability properties (as those encountered in fluid mechanics for instance). This leads to a stabilized finite element method cast in a primal-dual formulation, where weakly consistent stabilization terms are minimized under the constraint of the ill-posed weak formulation. The idea is introduced in [START_REF] Burman | Stabilized finite element methods for nonsymmetric, noncoercive, and ill-posed problems. Part I: Elliptic equations[END_REF] for standard H 1 -conforming finite element methods. The analysis of the ill-posed case is carried out in [START_REF] Burman | Error estimates for stabilized finite element methods applied to ill-posed problems[END_REF], and the approach is extended to nonconforming approximations in [START_REF] Burman | A stabilized nonconforming finite element method for the elliptic Cauchy problem[END_REF]. Notice that in both cases, the focus is on low-order approximation spaces. The error analysis relies on the existence of sharp conditional stability estimates for the continuous problem. The estimates are conditional in the sense that a certain a priori bound must be assumed to hold for the solution, and the continuity expressed in this bound is typically only Hölder [START_REF] John | Continuous dependence on data for solutions of partial differential equations with a prescribed bound[END_REF]. Such estimates are known in the literature as quantitative uniqueness results and use theoretical tools such as Carleman estimates or three-ball estimates [START_REF] Isakov | Inverse problems for partial differential equations[END_REF][START_REF] Alessandrini | The stability for the Cauchy problem for elliptic equations[END_REF]. Error bounds derived using conditional stability estimates can be optimal in the sense that they reflect the approximation order of the finite element space and the stability of the ill-posed problem. In particular, this means that when applied to a well-posed problem, the finite element method recovers optimal convergence.

The ill-posed problem that we consider here is the unique continuation problem. This is a model problem for data assimilation where the data on the boundary of the domain is unknown, and instead measurements are available in some subset of the bulk domain. It is known that such data can be extended in a unique way subject to an elliptic operator, and our objective is to compute an accurate approximation of this unique continuation by means of finite elements. In the case of ill-posed second-order elliptic problems, the constants in the conditional stability estimates depend in a nontrivial way on the physical parameters of the problem, on the geometry of the subdomain where data are available, and on the geometry of the target subdomain where accuracy of the computed solution is desired. In the case of the Helmholtz equation, it has been shown that the geometry of the data set relative to the geometry of the target domain is crucial for robustness of the estimates as the frequency grows [START_REF] Hrycak | Increased stability in the continuation of solutions to the Helmholtz equation[END_REF]. Bounds that are robust and also suitable to be used in the analysis of numerical methods were recently derived in [START_REF] Burman | Unique continuation for the Helmholtz equation using stabilized finite element methods[END_REF] for the Helmholtz equation and in [START_REF] Burman | A stabilized finite element method for inverse problems subject to the convection-diffusion equation. I: diffusion-dominated regime[END_REF] for convection-diffusion equations. In both of these references, conforming piecewise affine approximation spaces were considered. It is, however, well known that high-order methods are (much) more accurate for the approximation of the Helmholtz equation, in particular to control the pollution error [START_REF] Melenk | Wavenumber explicit convergence analysis for Galerkin discretizations of the Helmholtz equation[END_REF]. In the ill-posed regime, however, the source of the pollution is twofold and includes both the well-known dispersion error (appearing in well-posed Helmholtz problems as well) and the pollution due to unstable modes induced by the ill-posed character.

The objective of the present work is to explore if the use of high-order methods in the context of primal-dual stabilized Galerkin methods is equally advantageous in the illposed case as it is in the well-posed case. Inspired by the approach proposed in [START_REF] Burman | Unique continuation for the Helmholtz equation using stabilized finite element methods[END_REF] for the lowest-order finite element discretization of the unique continuation problem subject to the Helmholtz equation, we derive here a hybridized discontinuous Galerkin (dG) method of arbitrary order k ≥ 1 for the same problem. The stabilization consists of the standard penalty term on the solution jumps supplemented by a Tikhonov regularization term that is scaled with the mesh-size h in such a way that it does not perturb the accuracy of the method. We then prove error estimates following the two-step approach introduced in [START_REF] Burman | Error estimates for stabilized finite element methods applied to ill-posed problems[END_REF].

• Step 1: we derive a discrete stability estimate where we show that the method offers improved stability on certain residual quantities compared to the standard Galerkin method. Using this estimate we can prove that these residual quantities will converge with an optimal rate for a smooth solution (up to perturbations).

•

Step 2: we show that the conditional stability estimates of [START_REF] Burman | Unique continuation for the Helmholtz equation using stabilized finite element methods[END_REF] can be applied to the error equation whose right-hand side only depends on the above residual quantities.

The combination of the convergence rate derived in the previous step and an a priori bound on the discrete solution allows us to conclude.

Observe that Step 1 essentially makes use of fairly well-known finite element analysis techniques and gives no information on the size of the actual error. It only ensures that the residual of the computation goes to zero. To prove convergence in the second step, we need a stability estimate relating the residual to the computational error, and contrary to standard finite element methods, we also need here an a priori bound on the discrete solution. Moreover, deriving this a priori bound is more challenging in the high-order case than it is in the low-order case. This constitutes our main result; see Theorem 10. We notice that related approaches based on the weak Galerkin method were proposed in [START_REF] Wang | A new primal-dual weak Galerkin finite element method for ill-posed elliptic Cauchy problems[END_REF][START_REF] Wang | Primal-dual weak Galerkin finite element methods for elliptic Cauchy problems[END_REF] for the approximation of the elliptic Cauchy problem subject to the Poisson problem. In those references, however, the ill-posed character of the problem was not taken into account. That is, only Step 1 above was considered.

In our analysis, we also consider the estimation of the error in local averages of the solution and find that the averaging weight can be constructed in such a manner that the averaged error converges independently of the stability properties of the problem (see Proposition 11). To the best of our knowledge, this is the first result in computational ill-posed problems that shows that nonglobal average quantities can be computed in a stable way. In addition, we show that the averaging weight function may be constructed in such a way that the radius of its support decreases as the inverse of the frequency. In all of our estimates, we also track the effect of perturbations of data. This is a subtle exercise in the ill-posed case, since perturbations of data can make the computation diverge when the regularization parameter is much smaller than the perturbation level. Similarly to standard Tikhonov regularization, this allows one to connect the regularization parameter to the noise level, but here it also gives an estimate of what scales are possible to resolve with the given data, since the connection to the mesh parameter is immediate.

The paper is organized as follows. In Section 2, we introduce the model problem and its ill-posed weak formulation. We then recall some stability results for the continuous problem in Section 3. The numerical method is introduced in Section 4 and analyzed in Section 5. The error analysis includes error estimates hinging on conditional stability on the one hand and using local averages on the other hand. Finally, in Section 6, we study numerically the performance of the method on some test cases of varying difficulty concerning ill-posedness and data perturbation.

Model problem

Let Ω be a polygonal/polyhedral domain in R d , d ∈ {2, 3} (open, bounded, connected, Lipschitz subset of R d ) and let be some (open, bounded, connected, Lipschitz) subset of Ω. We consider the following continuation problem subjected to the Helmholtz equation:

-∆u -ω 2 u = f in Ω, (1a) u = g in , (1b)
with f ∈ L 2 (Ω) and g ∈ H 1 ( ). In this study, the frequency ω is a positive real number, and u is real-valued. To make problem (1) nontrivial, we assume that Ω. The function g is assumed to be the restriction to of a solution to (1a), so that (1) has a unique solution (ω can even be a resonance frequency since uniqueness is ensured by (1b)). However, problem (1) is ill-posed in the sense of Hadamard since there is no stability with respect to the data f and g (see [START_REF] Burman | Solving ill-posed control problems by stabilized finite element methods: an alternative to Tikhonov regularization[END_REF] for a detailed discussion). Notice in particular that no boundary condition is enforced to complement the PDE (1a), and that the use of the data g in (1b) compensates this lack of knowledge. The weak form of (1) can be written as follows: Find u ∈ H 1 (Ω) such that u| = g and

a(u, v) = (f, v) Ω , ∀v ∈ H 1 0 (Ω), (2)
where

a(u, v) := (∇u, ∇v) Ω -ω 2 (u, v) Ω , ∀u, v ∈ H 1 (Ω). (3)
Notice that in (2) we look for a function in H 1 (Ω) and use test functions in the smaller space H 1 0 (Ω). Here and in what follows, for a measurable subset S ⊆ Ω, we denote by (•, •) S the L 2 (S)-inner product with appropriate Lebesgue measure and by • S the induced norm. Moreover, we write S ⊂⊂ Ω whenever S ⊂ Ω.

Key stability estimates

We present here some analysis results from [START_REF] Alessandrini | The stability for the Cauchy problem for elliptic equations[END_REF][START_REF] Burman | Unique continuation for the Helmholtz equation using stabilized finite element methods[END_REF] on conditional stability estimates. As highlighted in the introduction, these results play an essential role in the error analysis. First, we give a stability estimate allowing one to bound some norm of any function in H 1 (Ω) in any subset B ⊂⊂ Ω (typically larger than ) by some norm defined in Ω and some norm defined in and involving the dual norm of the Helmholtz operator. For any subset S ⊆ Ω, we define

E S (v) := v S + Ω ∆v + ω 2 v H -1 (Ω) , ∀v ∈ H 1 (Ω), (4)
where the length scale Ω := diam(Ω) is introduced to make the expression of E S (v) dimensionally consistent. Notice that the dual residual norm is evaluated in Ω.

Lemma 1 (Conditional stability). Let Ω and B ⊂ Ω, with (B \ ) ⊂⊂ Ω. There are C(ω) and α ∈ (0, 1] such that

∇v B + ω v B ≤ C(ω)ωE Ω (v) 1-α E (v) α , ∀v ∈ H 1 (Ω). ( 5
)
Both C(ω) and α depend on the geometry of , B, and Ω, α is independent of the frequency ω, whereas C(ω) in general depends on ω. Left panel: specific configuration defined in [START_REF] Burman | A stabilized nonconforming finite element method for the elliptic Cauchy problem[END_REF]. Right panel: configuration also satisfying the bound [START_REF] Burman | Continuous interior penalty hp-finite element methods for advection and advection-diffusion equations[END_REF].

Proof. See [START_REF] Burman | Unique continuation for the Helmholtz equation using stabilized finite element methods[END_REF]. The reader can also refer to [1, Thm. 5.1] and [START_REF] Robbiano | Théorème d'unicité adapté au contrôle des solutions des problèmes hyperboliques[END_REF].

In the estimate (5), the constant C(ω) generally depends on ω (see [START_REF] Burman | Unique continuation for the Helmholtz equation using stabilized finite element methods[END_REF]Example 4]). In particular, it is shown in [START_REF] Burman | Unique continuation for the Helmholtz equation using stabilized finite element methods[END_REF] that for any N ∈ N, C(ω) ≤ ω N cannot hold uniformly w.r.t. ω, meaning that the bound (5) is useful only for low frequencies. Moreover, the size of the constant depends on the geometry of the subdomains and B in a nontrivial way. A remarkable result from [START_REF] Hrycak | Increased stability in the continuation of solutions to the Helmholtz equation[END_REF] states that under a certain convexity condition on the domains and B, essentially that any straight line intersecting B also intersects , the stability estimate (5) holds true with constant independent of ω. Similar results suitable in the context of finite element analysis were then derived in [START_REF] Burman | Unique continuation for the Helmholtz equation using stabilized finite element methods[END_REF] using the following special geometry (see the left panel of Figure 1). Consider the half-space

H := {(x 1 , ..., x d ) | x 1 < 0}, let r > 0, β > 0, R > r such that r 2 + β 2 < ρ < R 2 + β 2 , set y := (β, 0, ..., 0) ∈ R d and define Ω := H ∩ B(0, R), := Ω \ B(0, r), B := (Ω \ B(y, ρ)). (6)
For this geometry, the following bound holds true.

Lemma 2 (Robustness with respect to ω). Let Ω, Ω, and B ⊂ Ω with (B \ ) ⊂⊂ Ω be defined as in [START_REF] Burman | A stabilized nonconforming finite element method for the elliptic Cauchy problem[END_REF]. There are C > 0 and α ∈ (0, 1], both depending on the geometry of , B and Ω but independent of ω, such that

∇v B + ω v B ≤ CωE Ω (v) 1-α E (v) α , ∀v ∈ H 1 (Ω). (7)
Proof. See Remark 3.1. No proof is given in [START_REF] Burman | Unique continuation for the Helmholtz equation using stabilized finite element methods[END_REF] for general configurations satisfying the convexity condition from [START_REF] Hrycak | Increased stability in the continuation of solutions to the Helmholtz equation[END_REF], but it has been shown how to modify the proof for other domains such as, for instance, the geometry in the right panel of Figure 1 [START_REF] Nechita | Unique continuation problems and stabilised finite element methods[END_REF]Example 4.1]. This geometry will be used as an example in our numerical tests.

Discretization by a hybridized dG method

In this section, we describe the hybridized dG method used to discretize the unique continuation problem subject to the Helmholtz equation.

Discrete setting

Let (T h ) h>0 be a family of matching meshes of Ω. In principle, the meshes can have cells that are polyhedra with planar faces in R d , and hanging nodes are also possible. However, the analysis below requires the mesh to be such that the underlying discontinuous polynomial approximation space has a global H 1 -conforming subspace with optimal approximation properties. For simplicity, we will therefore restrict the discussion to meshes composed of simplices. One can also use meshes composed of cuboids by replacing the polynomial space P k by Q k in the scheme described below. The mesh cells are conventionally taken to be open subsets of R d , and n T denotes the unit outward normal to the generic mesh cell T ∈ T h . For a subset S ⊂ R d , h S denotes the diameter of S, and for a mesh T h , the index h refers to the maximal diameter of the mesh cells. The mesh faces are collected in the set F h which is split into the set of the mesh interfaces, F int h , and the set of the boundary faces, F ∂ h . Any mesh interface F ∈ F int h is oriented by a fixed unit normal vector n F , and for a piecewise smooth function φ defined on both cells sharing F , [[φ]] F denotes the jump of φ across F in the direction of n F . To avoid technicalities, we assume henceforth that the mesh family (T h ) h>0 is quasi-uniform, and that all the meshes are fitted to the subset .

Let k ≥ 1 be the polynomial degree of the hybridized dG method. We denote by P k (S) the set of polynomials of total degree at most k on the subset S ⊆ Ω. The discrete unknowns are piecewise polynomials of degree k attached to the mesh faces and of degree k attached to the mesh cells. We define the discrete spaces

Û k h := U k T × U k F , U k T := T ∈T h P k (T ), U k F := F ∈F h P k (F ). (8)
For a generic pair vh ∈ Û k h , we write

vh := (v T , v F ), v T := (v T ) T ∈T h ∈ U k T , v F := (v F ) F ∈F h ∈ U k F . (9)
We denote by Û k h0 the linear subspace of Û k h in which all the degrees of freedom attached to the boundary faces are null. Since the discrete bilinear forms are devised cellwise, it is convenient to adopt a local viewpoint. For a generic pair vh ∈ Û k h , its degrees of freedom associated with a generic mesh cell T ∈ T h are denoted by

vT := (v T , v ∂T := (v F ) F ∈F ∂T ) ∈ Û k T := P k (T ) × P k (F ∂T ), (10)
where P k (F ∂T ) := F ∈F ∂T P k (F ) and F ∂T := {F ∈ F h | F ⊂ ∂T } collects the mesh faces composing the boundary of T .

The discrete problem

Let us first define the bilinear forms to discretize the Helmholtz equation. The local hybridized dG bilinear form on a generic mesh cell T ∈ T h is such that for all vT , ŵT ∈ Û k T ,

a T (v T , ŵT ) := (∇v T , ∇w T ) T -ω 2 (v T , w T ) T -(∇v T •n T , w T -w ∂T ) ∂T -(v T -v ∂T , ∇w T •n T ) ∂T , (11)
and the local stabilization bilinear form is defined such that for all vT , ŵT ∈ Û k T ,

s T (v T , ŵT ) := h -1 T (v T -v ∂T , w T -w ∂T ) ∂T . (12)
The global bilinear forms are then defined such that a h (v h , ŵh ) := T ∈T h a T (v T , ŵT ) and

s h (v h , ŵh ) := T ∈T h s T (v T , ŵT ) for all vh , ŵh ∈ Û k h .
In what follows, we consider that the only available measurements are perturbed by some noise. We denote by

g δ := g + δ (13)
the available measurements, where g is the value of the exact solution in (see (1b)) and δ is some random noise. We will see that the a priori error bound depends on the amplitude of δ (see Theorem 10). Noise can also be considered on the right-hand side f of the Helmholtz equation. The adaptation of the analysis is straightforward, but has been omitted to enhance the readability of the proofs.

We introduce the discrete Lagrangian

L h : Û k h × Û k h0 → R such that for all (v h , ζh ) ∈ Û k h × Û k h0 , L h (v h , ζh ) := 1 2 -2 Ω v T -g δ 2 + 1 2 s γ h (v h , vh ) - 1 2 σ h ( ζh , ζh ) + a h (v h , ζh ) -(f, ζ T ) Ω , (14)
where g δ is defined in [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems[END_REF], and the length scale Ω makes the expression of L h dimensionally consistent. We refer to the first argument of L h as the primal variable and to the second argument as the dual variable. Moreover in [START_REF] Dardé | An H div -based mixed quasi-reversibility method for solving elliptic Cauchy problems[END_REF], the bilinear form s γ h : Û k h × Û k h → R is built using the above-defined stabilization bilinear form s h and a discrete Tikhonov regularization in such a way that for all vh , ŵh ∈ Û k h ,

s γ h (v h , ŵh ) := s h (v h , ŵh ) + γ -2(k+1) Ω h 2k (v T , w T ) Ω , (15)
with γ ≥ 0 (we will need to assume γ > 0 to derive the error estimates), and the bilinear form σ h : Û k h0 × Û k h0 → R acts as a stabilization term on the dual variable and is defined such that for all ζh , ηh ∈ Û k h0 ,

σ h ( ζh , ηh ) := T ∈T h (∇ζ T , ∇η T ) T + s h ( ζh , ηh ). ( 16
)
Notice that ( 12) is a natural stabilization term in the context of hybridized dG methods, and its role is to weakly enforce the matching of the face unknowns and the trace of the cell unknowns in each mesh cell. Moreover, the Tikhonov regularization γ 15) provides a weak control on the L 2 -norm of the solution that is needed to bound the error term E (θ c T ); see step (iii) in the proof of Theorem 10.

-2(k+1) Ω h 2k (v T , w T ) Ω in (
The discrete unique continuation problem subject to the Helmholtz equation is defined by looking for a critical point of the discrete Lagrangian. In other words, we seek (û

h , ξh ) ∈ Û k h × Û k h0 such that m (u T , w T ) + s γ h (û h , ŵh ) + a h ( ŵh , ξh ) = m (g δ , w T ), ∀ ŵh ∈ Û k h , (17a)
a h (û h , ηh ) -σ h ( ξh , ηh ) = (f, η T ) Ω , ∀η h ∈ Û k h0 , (17b)
with m (φ, ψ) := -2 Ω (φ, ψ) for all φ, ψ ∈ L 2 (Ω). We will prove in Section 5 that the square linear problem ( 17) has a unique solution, and that (û h , ηh ) → (u, 0) as h → 0 when δ = 0. Since the discrete Lagrangian is convex w.r.t. the primal variable and concave w.r.t. the dual variable, the pair (û h , ξh ) ∈ Û k h × Û k h0 solves [START_REF] Ern | Finite elements. I Approximation and Interpolation[END_REF] iff it is a saddlepoint of L h . Notice that the convexity of L h w.r.t. vh involves only the L 2 ( )-norm and the stabilization bilinear form s h , whereas the concavity of L h w.r.t. ηh involves a stronger H 1 -like norm. Notice also that if σ h := 0, (û h , ξh ) solves [START_REF] Ern | Finite elements. I Approximation and Interpolation[END_REF]

iff ûh minimizes over Û k h the functional vh → 1 2 -2 Ω v T -g δ 2 + 1 2 s γ h (v h , vh ) under the constraint that ûh satisfies the discrete Helmholtz equation a h (û h , ηh ) = (f, η T ) Ω for all ηh ∈ Û k h0 .
Remark 4.1. (Mesh assumptions) Since the mesh T h is assumed to be fitted to the subset , we have

(v T , w T ) = T ∈T h (v T , w T ) T for all v T , w T ∈ U k T , where T h := {T ∈ T h | T ⊂ }.
Observe that even if the mesh is not fitted to , the subset of composed of the cells in T h can be used in the analysis, possibly assuming h small enough if the geometry condition from Section 3 must hold. Furthermore, since the mesh family is assumed to be quasi-uniform, we use the global mesh-size h in the discrete Tikhonov regularization (15).

Static condensation

The discrete problem [START_REF] Ern | Finite elements. I Approximation and Interpolation[END_REF] can be solved efficiently by eliminating locally all the cell unknowns using static condensation. This local elimination leads to a global transmission problem on the mesh skeleton involving only the face unknowns with a stencil that couples unknowns attached to neighboring faces (in the sense of cells). Once this global transmission problem is solved, the cell unknowns are recovered by local solves. We refer the reader to, e.g., [START_REF] Cicuttin | Implementation of Discontinuous Skeletal methods on arbitrary-dimensional, polytopal meshes using generic programming[END_REF].

Let us verify that the cell unknowns can indeed be eliminated locally by expressing them in terms of the data and the face unknowns on F ∂T . Fix T ∈ T h and (u ∂T , ξ ∂T ) ∈ [P k (F ∂T )] 2 . To eliminate the cell unknowns, we want to solve the following problem: Find (u T , ξ T ) ∈ [P k (T )] 2 such that, setting ûT := (u T , u ∂T ) and ξT := (ξ T , ξ ∂T ), the pair (û T , ξT ) solves

       -2 Ω (u T , w T ) T ∩ + s T (û T , ŵT ) + γ -2(k+1) Ω h 2k (u T , w T ) T + a T ( ŵT , ξT ) = -2 Ω (g δ , w T ) T ∩ , a T (û T , ηT ) -(∇ξ T , ∇η T ) T -s T ( ξT , ηT ) = (f, η T ) T , (18)
for all ŵT := (w T , 0) with w T ∈ P k (T ) and all ηT := (η T , 0) with η T ∈ P k (T ).

Lemma 3 (Local problem). Problem (18) admits a unique solution for every γ ≥ 0.

Proof. Since [START_REF] Hadamard | Lectures on Cauchy's Problem in Linear Partial Differential Equations[END_REF] amounts to a square linear system, we only need to show uniqueness of the solution for zero data and u ∂T = ξ ∂T = 0. First, if γ > 0, we take w T := u T and η T = -ξ T and get

-2 Ω u T 2 T ∩ + s T (û T , ûT ) + γ -2(k+1) Ω h 2k u T 2 T + ∇z T 2 T + s T ( ξT , ξT ) = 0.
This implies that u T = ξ T = 0. In the case where γ = 0, we can prove again that u T = ξ T = 0 by choosing v T := h 2 T (-∆u T -ω 2 u T ). Proceeding as in step (iii) of the proof of Lemma 5 below, one shows that ∆u T + ω 2 u T T = 0, which proves that u T = 0 since ω > 0 by assumption. (Notice that if ω = 0, we get ∆u T T = 0, and since s T (û T , ûT ) = 0, we have u T | ∂T = 0, and we conclude that u T = 0.)

Analysis

In this section, we analyze the convergence of the discretization method introduced above. In what follows, we often abbreviate A B the inequality A ≤ CB for positive real numbers A and B, where the constant C is independent of h and ω. At this stage, the constant C only depends on the mesh shape-regularity and the polynomial degree k. For the reader's convenience, we collect in the following lemma several classical inequalities used in the analysis; see, e.g., [START_REF] Ern | Finite elements. I Approximation and Interpolation[END_REF]Chap. 12] for proofs.

Lemma 4 (Analysis tools). Let k ∈ N. There is C > 0, depending on k, the mesh regularity, and d, such that for all q ∈ P k (T ) and all T ∈ T h , the following discrete trace and inverse inequalities hold:

q ∂T ≤ Ch -1 2 T q T , ∇q T ≤ Ch -1 T q T . (19)
Moreover, there is C, depending on the mesh regularity and k, such that for all v ∈ H 1 (T ) and all T ∈ T h , the trace inequality

v ∂T ≤ C h -1 2 T v T + h 1 2
T ∇v T [START_REF] Isakov | Inverse problems for partial differential equations[END_REF] holds, as does the Poincaré inequality v -v T ≤ Ch T ∇v T , where v denotes the mean value of v over T .

Residual stability

In this section, we will prove the inf-sup stability (and hence the well-posedness) of the discrete problem (17) using a norm based on some residual quantities. At this stage, we can just assume that the regularization parameter satisfies γ ≥ 0. We first define for all

v T ∈ U k T + H 1 (Ω), v T 2 R := T ∈T h h 2 T ∆v T + ω 2 v T 2 T + F ∈F int h h F ∇v T F •n F 2 F . (21)
The notation U k T + H 1 (Ω) stands for the sum of the two linear spaces. Notice that the map v T → v T R defines a norm on U k T since a polynomial satisfying ∆v T + ω 2 v T = 0 is necessarily null (since ω > 0 by assumption). Consequently, the map vh

→ v T R + s h (v h , vh )
1 2 defines a norm on Û k h . We equip the product space Û k h × Û k h0 with the following norm:

|||v h , ζh ||| 2 := -2 Ω v T 2 + s γ h (v h , vh ) + v T 2 R + σ h ( ζh , ζh ). ( 22
)
Since σ h (•, •)

1 2 defines a norm on Û k h0 , |||•, •||| indeed defines a norm on Û k h × Û k h0 (notice that s h ≤ s γ h since γ ≥ 0).
To allow for a more compact notation, it is convenient to introduce the global bilinear form A h such that for all (v h , ζh ), ( ŵh , ηh )

∈ Û k h × Û k h0 , A h ((v h , ζh ), ( ŵh , ηh )) := m (v T , w T ) + s γ h (v h , ŵh ) + a h ( ŵh , ζh ) + a h (v h , ηh ) -σ h ( ζh , ηh ), (23)
so that ( 17) can be rewritten as follows:

Find (û h , ξh ) ∈ Û k h × Û k h0 such that A h ((û h , ξh ), ( ŵh , ηh )) = m (g δ , w T ) + (f, η T ) Ω ∀( ŵh , ηh ) ∈ Û k h × Û k h0 . (24)
Let us now prove an inf-sup condition on the bilinear form A h w.r.t. the triple norm defined in [START_REF] John | Continuous dependence on data for solutions of partial differential equations with a prescribed bound[END_REF].

Lemma 5 (Inf-sup stability). Assume that γ ≥ 0. The following holds true for all

(v h , ζh ) ∈ Û k h × Û k h0 : |||v h , ζh ||| sup ( ŵh ,η h )∈ Û k h × Û k h0 \{(0,0)} |A h ((v h , ζh ), ( ŵh , ηh ))| ||| ŵh , ηh ||| . (25)
Proof. Let (v h , ζh ) ∈ Û k h × Û k h0 and let S denote the supremum on the right-hand side of [START_REF] Lehrenfeld | High order exactly divergence-free hybrid discontinuous Galerkin methods for unsteady incompressible flows[END_REF]. (i) Owing to the definition (23) of the bilinear form A h , we infer that -2

Ω v T 2 + s γ h (v h , vh ) + σ h ( ζh , ζh ) = A h ((v h , ζh ), (v h , -ζh )) S|||v h , ζh |||, (26)
so that it only remains to bound v T 2 R . (ii) We consider ηh := (0, (η

F ) F ∈F h ) with η F := h F ∇v T F •n F for all F ∈ F int h and η F := 0 for all F ∈ F ∂ h (notice that ηh ∈ Û k h0 ). Since |||0, ηh ||| = σ h (η h , ηh ) 1 2
, we have

F ∈F int h h F ∇v T F •n F 2 F = a h (v h , ηh ) = A h ((v h , ζh ), (0, ηh )) + σ h ( ζh , ηh ) ≤ Sσ h (η h , ηh ) 1 2 + σ h ( ζh , ηh ) ≤ S + σ h ( ζh , ζh ) 1 2 σ h (η h , ηh ) 1 2
,

where the last bound follows from the Cauchy-Schwarz inequality applied to σ h . Moreover, we have

σ h (η h , ηh ) = T ∈T h h -1 T η ∂T 2 ∂T F ∈F int h h F ∇v T F •n F 2
F . This implies that

F ∈F int h h F ∇v T F •n F 2 F S 2 + σ h ( ζh , ζh ). ( 27
) (iii) We now consider ηh := ((η T ) T ∈T h , 0) with η T := -h 2 T (∆v T + ω 2 v T ) for all T ∈ T h (notice that ηh ∈ Û k h0
). Recalling the definition [START_REF] Burman | Space time stabilized finite element methods for a unique continuation problem subject to the wave equation[END_REF] of the local bilinear form a T and integrating by parts in each mesh cell, we have

a h (v h , ηh ) = T ∈T h (∇v T , ∇η T ) T -ω 2 (v T , η T ) T -(∇v T •n T , η T ) ∂T -(v T -v ∂T , ∇η T •n T ) ∂T = T ∈T h (-∆v T -ω 2 v T , η T ) T -(v T -v ∂T , ∇η T •n T ) ∂T .

This implies that

T ∈T h h 2 T ∆v T + ω 2 v T 2 T = A h ((v h , ζh ), (0, ηh )) + σ h ( ζh , ηh ) + T ∈T h (v T -v ∂T , ∇η T •n T ) ∂T Sσ h (η h , ηh ) 1 2 + σ h ( ζh , ζh ) + s h (v h , vh ) 1 2 σ h (η h , ηh ) 1 2
,

where we used the discrete trace inequality h

1 2
T ∇η T ∂T ∇η T T (see Lemma 4) and the definitions of s h and σ h . Moreover, invoking the inverse inequality ∇η T T h -1 T η T T and the discrete trace inequality η T ∂T h

-1 2 T η T T , we infer that σ h (η h , ηh ) T ∈T h h -2 T η T 2
T . This implies that

T ∈T h h 2 T ∆v T + ω 2 v T 2 T S 2 + s h (v h , vh ) + σ h ( ζh , ζh ). ( 28
)
(iv) Combining the estimates ( 26), [START_REF] Nechita | Unique continuation problems and stabilised finite element methods[END_REF], and [START_REF] Robbiano | Théorème d'unicité adapté au contrôle des solutions des problèmes hyperboliques[END_REF], and recalling that s h ≤ s γ h , we infer that |||v h , ζh ||| 2 S|||v h , ζh ||| + S 2 , which readily proves the expected inf-sup condition.

Lemma 6 (Unique solvability). The discrete problem (17) admits a unique solution for all γ ≥ 0.

Proof. The proof is a direct consequence of Lemma 5 since ( 17) is equivalent to a square linear system.

Consistency and a priori residual bound

We denote by u the solution to (1) and by (û h , ξh ) the solution to [START_REF] Ern | Finite elements. I Approximation and Interpolation[END_REF]. Let Îk

h (u) := (I k T (u), I k F (u)) ∈ Û k h be defined in each mesh cell T ∈ T h by (I k T (u), I k ∂T (u)) with I k T (u) the L 2 -projection of u onto P k (T ) and I k ∂T (u)| F the L 2 -projection onto P k (F ) for all F ∈ F ∂T .
We define the discrete approximation error and the interpolation error as follows:

êh := ûh -Îk h (u) ∈ Û k h , θ T := u -I k T (u) ∈ H 1 (Ω) + U k T . ( 29
)
Lemma 7 (Consistency). Assume that there exists ε > 0 such that u ∈ H 3/2+ε (Ω).

Assume that γ ≥ 0. For all ηh ∈ U k h0 , we have the estimate

|a h (ê h , ηh ) -σ h ( ξh , ηh )| θ T # σ h (η h , ηh ) 1 2 , ( 30
)
with θ T 2 # := T ∈T h ∇θ T 2 T + (ωh T ) 4 h -2 T θ T 2 T + h T ∇θ T 2 ∂T + h -1 T θ T 2 ∂T . (31)
Proof. (i) Using (17b), (1a), and integrations by parts, we get

a h (ê h , ηh ) -σ h ( ξh , ηh ) = (f, η T ) Ω -a h ( Îk h (u), ηh ) = T ∈T h -(∆u + ω 2 u, η T ) T -(∇I k T (u), ∇η T ) T + ω 2 (I k T (u), η T ) T + (∇I k T (u)•n T , η T -η ∂T ) ∂T + (I k T (u) -I k ∂T (u), ∇η T •n T ) ∂T = T ∈T h (∇θ T , ∇η T ) T -ω 2 (θ T , η T ) T -(∇θ T •n T , η T -η ∂T ) ∂T -(θ T , ∇η T •n T ) ∂T , (32)
where we used that T ∈T h (∇u•n T , η ∂T ) ∂T = 0 since ∇u and η ∂T are single-valued on all F ∈ F int h and η ∂T | F = 0 for all F ∈ F ∂ h , and that (u

-I k ∂T (u), ∇η T •n T ) ∂T = 0 since ∇η T •n T ∈ P k (F ∂T ).
(ii) We need to bound the four terms on the right-hand side of (32). For the first and the third terms, we just invoke the Cauchy-Schwarz inequality. For the second term, we observe that (θ T , η T ) T = (θ T , η T -η T ) T , where η T is the mean value of η T over T . Owing to the Cauchy-Schwarz inequality and the Poincaré inéquality (see Lemma 4), we infer that

|(θ T , η T ) T | θ T T η T -η T T h T θ T T ∇η T T .
Finally, the fourth term is bounded by means of the Cauchy-Schwarz inequality and the discrete trace inequality from Lemma 4.

We can now derive an a priori residual bound, that is, a bound on the discrete approximation error (ê h , ξh ) in the triple norm defined in [START_REF] John | Continuous dependence on data for solutions of partial differential equations with a prescribed bound[END_REF]. First, we recall an approximation estimate for the interpolation operator Îk h .

Lemma 8 (Approximation). Let s > 3 2 and set r := min(s, k + 1). Assume that ωh ≤ 1. The following holds true for all v ∈ H s (Ω):

v -I k T (v) # + v -I k T (v) R + s h ( Îk h (v), Îk h (v)) 1 2 + v -I k T (v) h r-1 |v| H r (Ω) . ( 33
)
Proof. The claim follows using the trace inequality [START_REF] Isakov | Inverse problems for partial differential equations[END_REF] and standard approximation results for the L 2 -projection.

Lemma 9 (A priori residual bound). Assume that γ ≥ 0 and ωh ≤ 1. Assume that the exact solution is in H k+1 (Ω), and set C(u)

:= |u| H k+1 (Ω) + -(k+1) Ω
u Ω . The following residual estimate holds true:

|||ê h , ξh ||| + s h (û h , ûh ) 1 2 h k max(1, γ) 1 2 C(u) + -1 Ω δ . (34)
Proof. (i) Using inf-sup stability (Lemma 5), we have

|||ê h , ξh ||| sup ( ŵh ,η h )∈ Û k h × Û k h0 \{(0,0)} |A h ((ê h , ξh ), ( ŵh , ηh ))| ||| ŵh , ηh ||| .
(ii) Using the definition of A h and consistency (Lemma 7; notice that by assumption u ∈ H 1+ (Ω) with > 1 2 since k ≥ 1), we infer that for all ηh ∈ Û k h0 ,

|A h ((ê h , ξh ), (0, ηh ))| = |a h (ê h , ηh ) -σ h ( ξh , ηh )| θ T # σ h (η h , ηh ) 1 2 ,
where θ T is defined in ( 29) and θ T # in [START_REF] Wang | Primal-dual weak Galerkin finite element methods for elliptic Cauchy problems[END_REF]. Moreover, using again the definition of A h , the discrete equation (17a), and the fact that u| = g, we have for all ŵh ∈ Û k h ,

|A h ((ê h , ξh ), ( ŵh , 0))| = |a h ( ŵh , ξh ) + s γ h (ê h , ŵh ) + m (e T , w T )| = | -s γ h ( Îk h (u), ŵh ) + m (θ T , w T ) + m (δ, w T )| ≤ s γ h ( Îk h (u), Îk h (u)) 1 2 + -1 Ω θ T + -1 Ω δ ||| ŵh , 0|||, since s γ h ( ŵh , ŵh ) + -2 Ω w T 2 ≤ ||| ŵh , 0||| 2 . (iii) It follows from the three above bounds that |||ê h , ξh ||| θ T # + s γ h ( Îk h (u), Îk h (u)) 1 2 + -1 Ω θ T + -1 Ω δ .
Moreover, the above right-hand side also bounds s h (û h , ûh )

1 2 since s h (û h , ûh ) 1 2 ≤ s h (ê h , êh ) 1 2 + s h ( Îk h (u), Îk h (u)) 1 2 ≤ |||ê h , ξh ||| + s γ h ( Îk h (u), Îk h (u)) 1 2 .
Recalling that s γ h ( Îk

h (u), Îk h (u)) = s h ( Îk h (u), Îk h (u)) + γ -2(k+1) Ω h 2k I k T (u)
2 Ω , the claim now follows since on the one hand

θ T # + s h ( Îk h (u), Îk h (u)) 1 2 + -1 Ω θ T h k |u| H k+1 (Ω) ,
where we used Lemma 8 and h -1 Ω ≤ 1 (since h T ≤ Ω = diam(Ω) for all T ∈ T h ), and on the other hand γ

-2(k+1) Ω h 2k I k T (u) 2 Ω ≤ γ -2(k+1) Ω h 2k u 2 Ω , owing to the L 2 -stability of I k T .

Error estimate using conditional stability

The objective of this section is to combine the conditional stability estimates on the continuous problem from Section 3 with the a priori residual bound from Lemma 9 to prove an error estimate in some subset B ⊂⊂ Ω. Since we are manipulating piecewise smooth functions (e.g., those in U k T ), we introduce the broken gradient operator ∇ T :

H 1 (Ω) + U k T → L 2 (Ω; R d ) acting cellwise. For a function v T ∈ H 1 (Ω) + U k T , we consider the norm v T 1,B := ∇ T v T B + ω v T B . (35)
Since the conditional stability bounds from Section 3 hold for functions in H 1 (Ω), we need to introduce a continuous interpolation of the discrete solution. Since the mesh is assumed to be simplicial, there exists an H 1 -conforming subspace of U k T with optimal approximation properties; this space can be built using Lagrange finite elements. Let π c

T denote the L 2 -orthogonal projection onto this continuous subspace. We will use the following well-known approximation result:

h -1 π c T (v T ) -v T Ω + ∇ T (π c T (v T ) -v T ) Ω ≤ s h (v h , vh ) 1 2 , ∀v h ∈ Û k h . (36)
The estimate (36) is proved using the construction of a particular interpolant into the space of continuous functions, together with the stability of the L 2 -projection (recall that we assume that the mesh sequence is quasi-uniform). For details, we refer the reader to, e.g., [7, Lemma 3.2 and Lemma 5.3].

Theorem 10 (Error estimate). Let u solve (1) and let (û h , ξh ) solve [START_REF] Ern | Finite elements. I Approximation and Interpolation[END_REF]. Recall the notation ûh := (u T , u F ) (see [START_REF] Burman | Unique continuation for the Helmholtz equation using stabilized finite element methods[END_REF]) and that δ denotes the noise in the measurements (see [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems[END_REF]).

Assume that u ∈ H k+1 (Ω) and recall that C(u)

:= |u| H k+1 (Ω) + -(k+1) Ω
u Ω . Assume that γ > 0 and that ωh ≤ 1. Let C(ω) > 0 and α ∈ (0, 1] result from Lemma 1 (recall that both real numbers depend on the geometry of and B, that α is independent of the frequency ω, and that the same holds true for C(ω) for the geometry defined in (6)). The following holds true:

u -u T 1,B C(ω)ω(h/ Ω ) kα k+1 Ω C(u) + (h/ Ω ) -k δ , (37)
where the hidden constant depends linearly on max(γ, 1, γ -1 ) 1 2 (and as above on the mesh shape-regularity and the polynomial degree k).

Proof. (i) Adding and subtracting π c

T (u T ), followed by the triangle inequality, and using (36), we infer that

u -u T 1,B ≤ u -π c T (u T ) 1,B + (1 + ωh)s h (û h , ûh ) 1 2 .
Owing to the assumption ωh ≤ 1 and invoking Lemma 9, we obtain

(1 + ωh)s h (û h , ûh ) 1 2 h k max(1, γ) 1 2 C(u) + -1 Ω δ . (ii) It remains to bound θ c T 1,B with θ c T := u -π c T (u T ). Since θ c T ∈ H 1 (Ω), we have θ c T 1,B = ∇θ c T B + ω θ c T B
, and Lemma 1 gives

θ c T 1,B ≤ C(ω)ωE Ω (θ c T ) 1-α E (θ c T ) α . ( 38
)
Thus, we need bounds on θ c T , θ c T Ω , and ∆θ c T + ω 2 θ c T H -1 (Ω) . (iii) Using the triangle inequality, we obtain

θ c T ≤ u -I k T (u) + I k T (u) -u T + u T -π c T (u T ) h k+1 |u| H k+1 ( ) + e T + hs h (û h , ûh ) 1 2 Ω h k max(1, γ) 1 2 C(u) + δ , (39)
where we used Lemma 8 and (36) on the second line, and the residual a priori estimate from Lemma 9 on the third line (notice that e T ≤ Ω |||ê h , ξh ||| and h ≤ Ω ). Proceeding in a similar way leads to

θ c T Ω h k+1 |u| H k+1 (Ω) + e T Ω + hs h (û h , ûh ) 1 2 Ω h k max(1, γ) 1 2 C(u) + δ + e T Ω .
This time we invoke the discrete Tikhonov regularization to control e T Ω . Owing to Lemma 9, we infer that

e T Ω ≤ γ -1 2 k+1 Ω h -k s γ h (ê h , êh ) 1 2
k+1

Ω max(γ -1 , 1)

1 2 C(u) + γ -1 2 k Ω h -k δ .
Combining the last two bounds and since h ≤ Ω , we obtain

θ c T Ω k+1 Ω max(γ -1 , 1, γ) 1 2 C(u) + max(γ -1 , 1) 1 2 k Ω h -k δ . ( 40
)
(iv) For the dual residual norm, we observe that

∆θ c T + ω 2 θ c T H -1 (Ω) = sup η∈H 1 0 (Ω), ∇η Ω =1 a(θ c T , η),
with the bilinear form a defined in (3). To accommodate the use of a piecewise smooth function as the first argument in the bilinear form, let us set a T (v T , η)

:= (∇ T v T , ∇η) Ω - ω 2 (v T , η) Ω for all v T ∈ H 1 (Ω) + U k T and all η ∈ H 1 0 (Ω). We have a(θ c T , η) = a T (θ c T , η) = a T (u T -π c T (u T ), η) + a T (u -u T , η),
and we bound the two terms separately. For the first term on the right-hand side, we use the Cauchy-Schwarz inequality, (36), that π c T is the L 2 -orthogonal projection onto the continuous polynomial subspace, and that η -π c T (η) Ω h ∇η Ω . This gives

a T (u T -π c T (u T ), η) = (∇ T (u T -π c T (u T )), ∇η) Ω + ω 2 (u T -π c T (u T ), η) Ω = (∇ T (u T -π c T (u T )), ∇η) Ω + ω 2 (u T -π c T (u T ), η -π c T (η)) Ω (1 + ω 2 h 2 )s h (û h , ûh ) 1 2 ∇η Ω s h (û h , ûh ) 1 2 ∇η Ω ,
where we used the assumption that ωh ≤ 1. Let us now consider the second term. Let ηh := Î1 h (η) ∈ Û k h0 , i.e., we have η T := I 1 T (η) for all T ∈ T h and η F := I 1 F (η) for all F ∈ F h (recall that k ≥ 1). Invoking ( 2) and (17b), we infer that

a T (u -u T , η) = a(u, η) -a T (u T , η) = (f, η) Ω -a T (u T , η) = (f, η -η T ) Ω + a h (û h , ηh ) -σ h ( ξh , ηh ) -a T (u T , η) = -(∆u + ω 2 u, η -η T ) Ω + a h (û h , ηh ) -a T (u T , η) -σ h ( ξh , ηh ).
Recalling the definitions of the bilinear forms a h and a T and integrating by parts in each mesh cell, we infer that

a h (û h , ηh ) -a T (u T , η) = T ∈T h (∇u T , ∇(η T -η)) T -ω 2 (u T , η T -η) T -(∇u T •n T , η T -η ∂T ) ∂T -(u T -u ∂T , ∇η T •n T ) ∂T = T ∈T h (∆u T + ω 2 u T , η -η T ) T -(u T -u ∂T , ∇η T •n T ) ∂T -(∇u T •n T , η -η ∂T ) ∂T .
Since η -η ∂T is single-valued at any interface F ∈ F int h and vanishes at any boundary face F ∈ F ∂ h , the last term on the above right-hand side can be rearranged as follows:

T ∈T h (∇u T •n T , η -η ∂T ) ∂T = F ∈F int h ([[∇u T ]] F •n F , η -η F ) F .
Putting the above identities together shows that

a T (u -u T , η) = T ∈T h (∆(u -u T ) + ω 2 (u -u T ), η T -η) T + F ∈F int h ([[∇u T ]] F •n F , η F -η) F - T ∈T h (u T -u ∂T , ∇η T •n T ) ∂T -σ h ( ξh , ηh ).
Let us denote by L 1 , L 2 , L 3 the terms on the three lines of the above right-hand side.

Owing to the Cauchy-Schwarz inequality and the approximation properties of Î1 h , we infer that

|L 1 + L 2 | u -u T R ∇η Ω , where • R is defined in (21). Moreover, since σ h (η h , ηh ) + T ∈T h h T ∇η T 2 ∂T ∇η 2 Ω , we infer that |L 3 | s h (û h , ûh ) 1 2 + σ h ( ξh , ξh ) 1 2
∇η Ω .

In conclusion, we can estimate the dual residual norm as follows:

∆θ c T + ω 2 θ c T H -1 (Ω) u -u T R + s h (û h , ûh ) 1 2 + σ h ( ξh , ξh ) 1 2
, and using the triangle inequality and recalling the definition of the triple norm, we obtain

∆θ c T + ω 2 θ c T H -1 (Ω) u -I k T (u) R + s h (û h , ûh ) 1 2 + |||ê h , ξh |||.
Invoking the residual a priori estimate from Lemma 9 and the approximation result from Lemma 8 which imply that u -I k T (u) R h k |u| H k+1 (Ω) (recall that ωh ≤ 1 by assumption), we infer that

Ω ∆θ c T + ω 2 θ c T H -1 (Ω) Ω h k max(1, γ) 1 2 C(u) + δ . (41)
(v) Collecting the bounds (39), (40), and (41), we see that

E (θ c T ) α E Ω (θ c T ) 1-α max(γ, 1, γ -1 ) 1 2 (h k Ω C(u) + δ ) α ( k+1 Ω C(u) + (h/ Ω ) -k δ ) 1-α max(γ, 1, γ -1 ) 1 2 (h/ Ω ) αk ( k+1 Ω C(u) + (h/ Ω ) -k δ ).
Inserting this bound into (38) concludes the proof.

Remark 5.1. (Estimate (37)) The error estimate (37) indicates two situations regarding mesh refinement (h → 0) with a fixed amount of data perturbation (quantified by δ ). In the well-posed case (α = 1), the error upper bound reaches a plateau when the approximation error lies below the data perturbation. Instead, in the ill-posed case (α < 1), mesh refinement will eventually lead to poorer accuracy. This degeneration of the bound is a consequence of the fact that, in general, a perturbation of the data results in nonexistence of a solution to the unique continuation problem, since the perturbed g δ is not the restriction of a solution of (1a) to . Stagnation can be achieved also in the ill-posed case by imposing a lower bound depending on C(u), Ω , and δ on the coefficient in the Tikhonov regularization term of [START_REF] Di Pietro | A Hybrid High-Order locking-free method for linear elasticity on general meshes[END_REF]. If, for instance, one wishes to ensure that no more than a factor two can be lost due to data perturbations, the lower bound is obtained by equating the two contributions k+1 Ω C(u) and (h/ Ω ) -k δ . The second term on the right-hand side of ( 15) then becomes γ

-2(k+1) Ω max(h, h min ) 2k (v T , w T ) Ω with h min := δ Ω C(u) 1 k
, and stagnation occurs when h ≤ h min (in practice, C(u) must of course be estimated).

Stable approximation of averaged quantities

The ill-posed problem can be seen as a problem that is close to resonance for all values of ω. To see this, consider the following Helmholtz Cauchy problem [9, Example 3] with

n ∈ N \ {0}:      ∆z + ω 2 z = 0 in Ω := (0, π) × (0, 1), z(x, 0) = 0 for x ∈ [0, π], z y (x, 0) = sin(nx) for x ∈ [0, π].
For n > ω, the solution is given by

z(x, y) = 1 √ n 2 -ω 2 sin(nx) sinh( n 2 -ω 2 y), (42)
for n = ω by z(x, y) = sin(nx)y, and for n < ω by

z(x, y) = 1 √ ω 2 -n 2 sin(nx) sin( ω 2 -n 2 y). (43)
This is the Helmholtz equivalent to the classical Hadamard solution for the elliptic Cauchy problem [START_REF] Hadamard | Lectures on Cauchy's Problem in Linear Partial Differential Equations[END_REF], which is recovered for ω = 0. We observe that, contrary to the Laplace case, for the Helmholtz equation perturbations with n < ω do not grow exponentially, whereas for n > ω the exponential growth is actually moderated by the frequency ω. In particular, exponentially growing modes must oscillate with a higher frequency than the solution.

For the unique continuation problem, this means that regardless of ω, there are always functions in the kernel of the Helmholtz operator that are small in the vicinity of but grow exponentially away from the data domain. However, since these unstable modes are more oscillating than the solution and are characterized by the kernel of the Helmholtz operator, a natural question is whether some averaged quantities can be approximated with better stability. That is, can we design averages for which functions such as z are filtered out? Whether this is possible hinges upon the existence of certain functions that we will now characterize.

We assume that there exists a function φ fulfilling the following assumptions:

φ ∈ C 1 (Ω) ∩ H 2 (Ω), φ| ∂Ω = ∇ φ • n| ∂Ω = 0, (44a)
∆ φ + ω 2 φ ≤ 0 in Ω, ∆ φ + ω 2 φ L 1 (Ω) > 0. ( 44b
)
An explicit construction is proposed in Lemma 12 below. We define the normalized function ϕ := φ/ ∆ φ + ω 2 φ L 1 (Ω) . Let ψ := -∆ϕ -ω 2 ϕ. Then ψ ≥ 0 and Ω ψ = 1, so that it makes sense to consider the local average functional

j ψ (v) := (v, ψ) Ω . (45)
Let us now prove that the average error |j ψ (u -u T )| can be bounded independently of the stability properties of the unique continuation problem.

Proposition 11 (Average error). Let u solve (1) and let (û h , ξh ) solve [START_REF] Ern | Finite elements. I Approximation and Interpolation[END_REF]. Recall the notation ûh := (u T , u F ) (see [START_REF] Burman | Unique continuation for the Helmholtz equation using stabilized finite element methods[END_REF]) and that δ denotes the noise in the measurements (see [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems[END_REF]). Assume that u ∈ H k+1 (Ω) and recall that C(u)

:= |u| H k+1 (Ω) + -(k+1) Ω
u Ω . The following holds true:

|j ψ (u -u T )| (1 + ω 2 h 2 )h k C(u) + -1 Ω δ , (46)
where the hidden constant depends on φ (and as above on the mesh shape-regularity and the polynomial degree k).

Proof. Using the definition of ψ and integrating by parts cellwise, we infer that

j ψ (u -u T ) = T ∈T h (u -u T , -∆ϕ -ω 2 ϕ) T = T ∈T h (∇(u -u T ), ∇ϕ) T -ω 2 (u -u T , ϕ) T + (∇ϕ•n T , u T -u ∂T ) ∂T .
Since ûh satisfies (17b), using ηh := Îk h (ϕ) ∈ Û k h0 , we get

T ∈T h -(∇(u -u T ), ∇I k T (ϕ)) T + ω 2 (u -u T , I k T (ϕ)) T -(u T -u ∂T , ∇I k T (ϕ)•n T ) ∂T +(∇(u -u T )•n T , I k T (ϕ) -I k ∂T (ϕ)) ∂T -σ h ( ξh , Îk h (ϕ)) = 0.
We then have

j ψ (u -u T ) = T ∈T h (∇(u -u T ), ∇(ϕ -I k T (ϕ))) T -ω 2 (u -u T , ϕ -I k T (ϕ)) T + (∇(ϕ -I k T (ϕ))•n T , u T -u ∂T ) ∂T + (∇(u -u T )•n T , I k T (ϕ) -I k ∂T (ϕ)) ∂T -σ h ( ξh , Îk h (ϕ)).
Integrating by parts cellwise once again, we obtain the error representation

j ψ (u -u T ) = -σ h ( ξh , Îk h (ϕ)) + T ∈T h (∆(u T -u) + ω 2 (u T -u), ϕ -I k T (ϕ)) T + (∇(u -u T )•n T , ϕ -I k ∂T (ϕ)) ∂T + (∇(ϕ -I k T (ϕ))•n T , u T -u ∂T ) ∂T .
We observe that

T ∈T h (∇(u -u T )•n T , ϕ -I k ∂T (ϕ)) ∂T = F ∈F int h ( ∇(u -u T ) F •n F , ϕ -I k F (ϕ)) F .
Invoking the Cauchy-Schwarz inequality and the approximation properties of Îk h (ϕ), we see that

|j ψ (u -u T )| u -u T R + s h (û h , ûh ) + σ h ( ξh , ξh ) 1 2
,

where the hidden constant has the dependencies made in the statement. Invoking the a priori residual estimate from Lemma 9 proves the claim.

We close this section by giving an elementary construction of a suitable function φ to be used in the definition of the local average functional under the assumption that ω is large enough. Recall that d denotes the space dimension and that we assume d ≥ 2.

Lemma 12 (Suitable function φ). Assume (without loss of generality) that 0 ∈ Ω and let r 3 > 0 be such that B(0, r 3 ) ⊂ Ω. Assume that ωr 3 ≥ 2d. Then there exists φ ∈ C 1 (Ω) supported in B(0, r 3 ) fulfilling (44).

Proof. We build a suitable function φ with radial symmetry, i.e., φ ≡ φ(r). The assumption ωr 3 ≥ 2d implies that there exist δr > 0 and ζ ≥ d -2 such that

r 3 = (2 + ζ)δr, δr ≥ 1 ω 2 d + ζ 1 + ζ . ( 47
)
Indeed, these conditions can be realized if ωr 3 ≥ κ(ζ) := 2 d+ζ 1+ζ (2 + ζ), and since the function κ is increasing, a sufficient condition for the existence of ζ ≥ d -2 is ωr 3 ≥ κ(d -2) = 2d. Let us then set r 0 := 0, r 1 := ζδr, r 2 := (1 + ζ)δr, and recall that r 3 = (2 + ζ)δr. We build a piecewise quadratic function such that φ ∈ C 1 (R -), φ (r) = 0 for all r ∈ [r 0 , r 1 ], φ(r) = φ (r) = 0 for all r > r 3 , and φ (r

)| [r 1 ,r 2 ] = -φ (r)| [r 2 ,r 3 ] = χ > 0.
The parameter χ > 0 is irrelevant owing to the normalization step, and it is kept only to be dimensionally consistent. Using integration and the above design conditions leads to

φ(r) =          -χδr 2 r ∈ [r 0 , r 1 ], χ (r-r 2 +δr) 2 2 -χδr 2 r ∈ [r 1 , r 2 ], -χ (r-r 3 ) 2 2 r ∈ [r 2 , r 3 ], 0 r > r 3 .
Recall that the Helmholtz operator for functions with radial symmetry in d space dimensions is given by H r := d 2 dr 2 + (d-1) r d dr + ω 2 I, where I is the identity operator. Applying H r to φ leads to

H r ( φ) =            -ω 2 χδr 2 r ∈ [r 0 , r 1 ], χ d + (d -1) (δr-r 2 ) r + ω 2 (r-r 2 +δr) 2 2 -δr 2 r ∈ [r 1 , r 2 ], χ -d + (d -1) r 3 r -ω 2 (r-r 3 ) 2 2 r ∈ [r 2 , r 3 ], 0 r ≥ r 3 .
We have for all r ∈ [r 2 , r 3 ],

χ -1 H r ( φ) = -d + (d -1) r 3 r -ω 2 (r -r 3 ) 2 2 ≤ -d + (d -1) r 3 r ≤ -d + (d -1) (2 + ζ) (1 + ζ) ≤ 0, since r 3 = (2 + ζ)δr and 1 r ≤ 1 (1+ζ)δr in [r 2 , r 3 ]. This implies that H r ( φ)| [r 2 ,r 3 ] ≤ 0 under the condition ζ ≥ d -2.
Moreover, using that δr -r 2 = -ζδr, -1 r ≤ -1 (1+ζ)δr , and (r -r 2 + δr) 2 ≤ δr 2 for all r ∈ [r 1 , r 2 ], we have

χ -1 H r ( φ)| [r 1 ,r 2 ] ≤ d + ζ 1 + ζ -ω 2 δr 2 2 .
A sufficient condition for H r ( φ)| [r 1 ,r 2 ] ≤ 0 is then given by the second condition in (47), and altogether this proves that H r ( φ) ≤ 0. Finally, H r ( φ) L 1 (Ω) > 0 since (for instance) H r ( φ)| [r 0 ,r 1 ] < 0. An example of functions -φ and -H r ( φ) is given in Figure 2.

Remark 5.2. (Localization, normalization) The assumption ωr 3 ≥ 2d shows that the above construction breaks down in the limit ω → 0; i.e., it is not clear if there are any stable averages for the ill-posed Poisson problem. On the other hand, if ω is large enough for the above condition to be fulfilled, it is typically interesting to choose ζ as large as possible, i.e., ζ = χ -1 (ωr 3 ) so that the transition layer in the graph of φ is as thin as possible. Finally, we notice that straightforward computations show that the normalization factor ∆ φ + ω 2 φ L 1 (Ω) is of the order O(δr d (1 + (ωδr) 2 )), where δr d is proportional to the volume of the support of φ.

Numerical tests

In this section, we present numerical simulations to verify the results of Section 5. We consider the domain Ω := (0, 1) 2 so that Ω = √ 2. We use four meshes that are triangulations built from successive refinements of Cartesian meshes (see Figure 3). The corresponding maximal diameters are h = 0.05689, h = 0.02845, h = 0.01422 and h = 0.007111.

We consider two geometries for and B. The first geometry, with := Ω\ (0, 0.875)× (0.125, 0.875) and B := Ω \ (0, 0.125) × (0.125, 0.875) (see Figure 4), satisfies the bound (7) which is robust with respect to ω. The second geometry, with := (0.25, 0.75)× (0, 0.5) and B := (0.125, 0.875) × (0, 0.875) (see Figure 5), satisfies the more general bound (5) that is sensitive to the value of ω. Notice that the second geometry violates the condition (B \ ) ⊂⊂ Ω on a small part of the lower boundary. This allows us to verify numerically that such a small violation of the condition does not necessarily destroy the convergence behavior, even in the case of high-order approximation since a convergence of order h kα is nonetheless observed.

The numerical simulations are run by using the exact solution defined in (42) with n := 5 for ω = 1 and n := 11 for ω = 10, so that both solutions exhibit a similar exponential blowup (we have √ n 2 -ω 2 = √ 24 in the first case and √ n 2 -ω 2 = √ 21 in the second case). Note that f := 0 in both cases. We use the value γ := 1 for the Tikhonov regularization. All the errors are measured in the

H 1 -seminorm |u -u T | 2 H 1 (T h ) := T ∈T h ∇(u -u T ) 2 T ,
and to facilitate comparisons, we report relative errors, i.e., rescaled by |u| H 1 (Ω) . All the numerical tests are run using the static condensation procedure described in Section 4.3; this is a key step to ensuring the computational efficiency of the scheme.

Moderate perturbations

Perturbations are used to pollute the data in the following way. At every Gauss point, we add to the exact value of g a perturbation equal to rand() * ( -1 Ω h T ) k , where h T denotes the diameter of the cell T ∈ T h containing the Gauss point and rand() is a C++ function that returns a random number in (-1, 1).

The relative errors in B in the H 1 -seminorm for the first geometry are reported in Figure 6. We observe that we obtain optimal convergence rates for ω = 1 and ω = 10 (i.e., convergence of order k). The dual variable also converges at optimal rate k. Note that for ω = 10, the method with k = 1 does not immediately reach optimal convergence, contrary to the orders k ∈ {2, 3}. This corroborates the fact that high-order polynomials are more efficient when ω grows. Overall, as expected, the results are better for ω = 1.

The relative errors in B in the H 1 -seminorm for the second geometry are reported in Figure 7. Here, we observe that the ill-posedness of the problem is more severe than for the first geometry, so that the parameter α in the conditional stability estimate is expected to be smaller than one. This is confirmed by the numerical results since the convergence rates are reduced compared to Figure 6. We observe a convergence rate of about 0.5 for k = 2 and 1 for k = 3, and the convergence is very slow for k = 1. Notice that here the dual variable converges with a slightly reduced optimal rate (about k -1

2 ). This may indicate that the considered meshes are not fine enough to reach the asymptotic convergence rates.

to consider hybrid high-order (HHO) methods, as originally devised in [START_REF] Di Pietro | An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators[END_REF][START_REF] Di Pietro | A Hybrid High-Order locking-free method for linear elasticity on general meshes[END_REF]. The two main differences with the above hybridized dG method is that the HHO method uses a local gradient reconstruction in each mesh cell from the cell and the face unknowns and that it employs a stabilization bilinear form such that a higher-order form of consistency is achieved. One possibility in the HHO method is to consider a mixed-order setting where the cell unknowns are polynomials of degree (k+1) and the face unknowns are polynomials of degree k ≥ 0. The advantage of this choice is that the stabilization bilinear form is simpler to formulate and amounts to the Lehrenfeld-Schöberl (LS) stabilization [START_REF] Lehrenfeld | Hybrid Discontinuous Galerkin methods for solving incompressible flow problems[END_REF][START_REF] Lehrenfeld | High order exactly divergence-free hybrid discontinuous Galerkin methods for unsteady incompressible flows[END_REF] devised in the context of hybridizable discontinuous Galerkin (HDG) methods [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems[END_REF].

The HHO method with the LS stabilization is defined as follows. For every T ∈ T h , the local unknowns belong to Û k T := P k+1 (T ) × P k (F ∂T ).

(48)

We define a gradient reconstruction G k T : Û k T → P k (T ; R 2 ) such that

(G k T (v T ), q) T := (∇v T , q) T -(v T -v ∂T , q•n T ) ∂T , ∀q ∈ P k (T ; R 2 ), (49)
for all vT ∈ Û k T , and the local LS stabilization bilinear form is such that

s ls T (v T , ŵT ) := h -1 T (Π k ∂T (v T ) -v ∂T , Π k ∂T (w T ) -w ∂T ) ∂T , ∀v T , ŵT ∈ Û k T , (50)
where Π k ∂T stands for the L 2 -projection onto P k (F ∂T ). We then set

where the stabilization bilinear form sh results from the elementwise summation of the local LS stabilization defined in (50). At the theoretical level, the LS stabilization, which is crucial to delivering convergence rates of order (k + 1) in the H 1 -seminorm in the well-posed case, is problematic in the ill-posed case since the proof of Theorem 10 fails, especially the estimate (36), because we only control the projection Π k ∂T of the jump and not the full jump. Our numerical results indicate though that this theoretical bottleneck may be pessimistic, at least in some situations. For purposes of comparison, we also consider the above modifications in (52) but with the full stabilization bilinear form s h defined in [START_REF] Cicuttin | Implementation of Discontinuous Skeletal methods on arbitrary-dimensional, polytopal meshes using generic programming[END_REF]. Then the proof of Theorem 10 goes through, but the price to pay is that the H 1 -convergence rate is expected to be of order k only.

We consider perturbations of amplitude (h/ Ω ) k , defined as in the previous numerical tests (with θ = 0). We study only the first geometry. The results are reported in Figure 9 for ω = 1 and in Figure 10 for ω = 10. When using the LS stabilization, we observe better convergence rates for k ∈ {1, 2}, but convergence stops on the last refinement for k = 3. The results obtained for the HHO method with full stabilization are close to those of the hybridized dG scheme studied in the present work. These numerical tests let us think that it may be possible to obtain higher-order convergence on the present continuation problem with an HHO method; a further study could include a Tikhonov regularization and a noise level of order (k + 1) instead of k. However, a theoretical gap still remains, and we postpone further analysis to future work.

Conclusion

In this paper, we have proposed and analyzed a high-order hybridized discontinuous Galerkin method for the approximation of an ill-posed wave propagation problem in the frequency domain. Using conditional stability estimates, we have proven error estimates that reflect the approximation order of the polynomial space, the degree of ill-posedness as expressed by the available conditional stability estimates and the effect of perturbations in data. We have also shown that certain averaged quantities can be approximated with Lipschitz stability. The results were illustrated by a series of numerical examples. The advantages of the high-order approximation are clearly seen, in particular for higher wave numbers. However, care must be taken to stop the mesh refinement when the perturbations are stronger than the threshold allowed by the estimates, since otherwise the solution quality degenerates. Alternatively, the coefficient of the Tikhonov regularization term can be made lower bounded, on the scale of the perturbations, leading to stability on finer scales (without convergence). Several future research directions are possible-for instance, unique continuation subject to the wave equation (see also [START_REF] Burman | Space time stabilized finite element methods for a unique continuation problem subject to the wave equation[END_REF] for a related space-time approach) in the time domain or Maxwell's equations in the frequency domain.

Figure 1 :

 1 Figure 1: Configurations from [9] satisfying the convexity condition necessary for Lemma 2.Left panel: specific configuration defined in[START_REF] Burman | A stabilized nonconforming finite element method for the elliptic Cauchy problem[END_REF]. Right panel: configuration also satisfying the bound[START_REF] Burman | Continuous interior penalty hp-finite element methods for advection and advection-diffusion equations[END_REF].

[ 9 ,

 9 Cor. 3 & Lem. 2].

Figure 2 :

 2 Figure 2: Example of radial cross section of -φ (left) and -H r ( φ) (right). The parameter values are r 3 = 1.0, ζ = 3, δr = 0.2, ω = δr √ 2 (2 + ζ)/(1 + ζ), χ = 5.

Figure 3 :

 3 Figure 3: First two meshes used. Left: h = 0.05689. Right: h = 0.02845.

Figure 4 :

 4 Figure 4: First geometry. Left:(highlighted in dark gray). Right: B (highlighted in light gray).

Figure 5 :

 5 Figure 5: Second geometry. Left:(highlighted in dark gray). Right: B (highlighted in light gray).

Figure 6 :

 6 Figure 6: Relative H 1 -error in B for the first geometry and moderate data perturbation. Left: ω = 1. Right: ω 10.

a

  T (v T , ŵT ) := (G k T (v T ), G k T ( ŵT )) T -ω 2 (v T , w T ) T ,(51)and we consider the problem[START_REF] Ern | Finite elements. I Approximation and Interpolation[END_REF] with the following modifications:a h (v h , ŵh ) := T ∈T h a T (v T , ŵT ),(52a)s γ h (v h , ŵh ) := sh (v h , ŵh ) + γ -2(k+1) Ω h 2k (v T , w T ) Ω ,(52b)σ h ( ζh , ηh ) :=T ∈T h (∇ζ T , ∇η T ) T + sh ( ζh , ηh ),

Acknowledgment

EB was partly supported by EPSRC grants EP/P01576X/1 and EP/P012434/1.

Larger perturbations

We now pollute the data with larger perturbations. At every Gauss point, we add to the exact value of g a perturbation that is equal to rand() * ( -1 Ω h T ) k-θ with θ ∈ {1, 2}. Notice that the results for θ = 0 were presented in Section 6.1. We report the results for the first geometry and ω = 1 in Figure 8. We observe similar results for θ = 0 and θ = 1. However, for θ = 2 the convergence is slower. Very similar results are obtained for ω = 10. The second geometry is, however, less sensitive to the perturbations, and the results are similar for θ = 0, θ = 1, and θ = 2. This can be a consequence of the larger stability constants (the second geometry fulfills (5) instead of ( 7)).

On the use of a higher-order method

We want to evaluate the possibility of using a higher-order method. Since we are working with discrete unknowns attached to the mesh cells and the mesh faces, one possibility is