Inverse boundary value problem for the Helmholtz equation: quantitative conditional Lipschitz stability estimates This research was supported in part by the members, BGP, ExxonMobil, PGS, Statoil and Total, of the Geo-Mathematical Imaging Group now at Rice University.

# Inverse boundary value problem for the Helmholtz equation: quantitative conditional Lipschitz stability estimates ††thanks: This research was supported in part by the members, BGP, ExxonMobil, PGS, Statoil and Total, of the Geo-Mathematical Imaging Group now at Rice University.

Elena Beretta Dipartimento di Matematica “Brioschi”, Politecnico di Milano, Italy (elena.beretta@polimi.it).    Maarten V. de Hoop Department of Computational and Applied Mathematics and Department of Earth Science, Rice University, 6100 Main Street, Houston TX 77005, USA (mdehoop@rice.edu). The research of this author is supported in part by the Simons Foundation.    Florian Faucher INRIA Bordeaux Sud-Ouest Research Center, Team Project Magique-3D, France (florian.faucher@inria.fr).    Otmar Scherzer Computational Science Center, University of Vienna, Oskar-Morgenstern Platz 1, A-1090 Vienna, Austria (otmar.scherzer@univie.ac.at). The research of this author is supported by the Austrian Science Fund (FWF), Project P26687-N25 Interdisciplinary Coupled Physics Imaging.
###### Abstract

We study the inverse boundary value problem for the Helmholtz equation using the Dirichlet-to-Neumann map at selected frequencies as the data. A conditional Lipschitz stability estimate for the inverse problem holds in the case of wavespeeds that are a linear combination of piecewise constant functions (following a domain partition) and gives a framework in which the scheme converges. The stability constant grows exponentially as the number of subdomains in the domain partition increases. We establish an order optimal upper bound for the stability constant. We eventually realize computational experiments to demonstrate the stability constant evolution for three dimensional wavespeed reconstruction.

Key words. Inverse problems, Helmholtz equation, stability and convergence of numerical methods.

AMS subject classifications. 35R30, 86A22, 65N12, 35J25

## 1 Introduction

In this paper we study the inverse boundary value problem for the Helmholtz equation using the Dirichlet-to-Neumann map at selected frequencies as the data. This inverse problem arises, for example, in reflection seismology and inverse obstacle scattering problems for electromagnetic waves [3, 22, 4]. We consider wavespeeds containing discontinuities.

Uniqueness of the mentioned inverse boundary value problem was established by Sylvester & Uhlmann [21] assuming that the wavespeed is a bounded measurable function. This inverse problem has been extensively studied from an optimization point of view. We mention, in particular, the work of [5].

It is well known that the logarithmic character of stability of the inverse boundary value problem for the Helmholtz equation [1, 19] cannot be avoided, see also [14, 15]. In fact, in [17] Mandache proved that despite of regularity a priori assumptions of any order on the unknown wavespeed, logarithmic stability is the best possible. However, conditional Lipschitz stability estimates can be obtained: accounting for discontinuities, such an estimate holds if the unknown wavespeed is a finite linear combination of piecewise constant functions with an underlying known domain partitioning [6]. It was obtained following an approach introduced by Alessandrini and Vessella [2] and further developed by Beretta and Francini [7] for Electrical Impedance Tomography (EIT) based on the use of singular solutions. If, on one hand, this method allows to use partial data, on the other hand it does not allow to find an optimal bound of the stability constant. Here, we revisit the Lipschitz stability estimate for the full Dirichlet-to-Neumann map using complex geometrical optics (CGO) solutions which give rise to a sharp upper bound of the Lipschitz constant in terms of the number of subdomains in the domain partitioning. We develop the estimate in .

Unfortunately, the use of CGO’s solutions leads naturally to a dependence of the stability constant on frequency of exponential type. This is clearly far from being optimal as it is also pointed out in the paper of Nagayasu, Uhlmann and Wang [18]. There the authors prove a stability estimate, in terms of Cauchy data instead of the Dirichlet-to-Neumann map using CGO solutions. They derive a stability estimate consisting of two parts: a Lipschitz stability estimate and a Logarithmic stability estimate. When the frequency increases the logarithmic part decreases while the Lipschitz part becomes dominant but with a stability constant which blows up exponentially in frequency.

We can exploit the quantitative stability estimate, via a Fourier transform, in the corresponding time-domain inverse boundary value problem with bounded frequency data. Datchev and De Hoop [9] showed how to choose classes of non-smooth coefficient functions, one of which is consistent with the class considered here, so that optimization formulations of inverse wave problems satisfy the prerequisites for application of steepest descent and Newton-type iterative reconstruction methods. The proof is based on resolvent estimates for the Helmholtz equation. Thus, one can allow approximate localization of the data in selected time windows, with size inversely proportional to the maximum allowed frequency. This is of importance to applications in the context of reducing the complexity of field data. We note that no information is lost by cutting out a (short) time window, since the boundary source functions (and wave solutions), being compactly supported in frequency, are analytic with respect to time. We cannot allow arbitrarily high frequencies in the data. This restriction is reflected, also, in the observation by Blazek, Stolk & Symes [8] that the adjoint equation, which appears in the mentioned iterative methods, does not admit solutions.

As a part of the analysis, we study the Fréchet differentiability of the direct problem and obtain the frequency and domain partitioning dependencies of the relevant constants away from the Dirichlet spectrum. Our results hold for finite fixed frequency data including frequencies arbitrarily close to zero while avoiding Dirichlet eigenfrequencies; in view of the estimates, inherently, there is a finest scale which can be reached. Finally we estimate the stability numerically and demonstrate the validity of the bounds, in particular in the context of reflection seismology.

## 2 Inverse boundary value problem with the Dirichlet-to-Neumann map as the data

### 2.1 Direct problem and forward operator

We describe the direct problem and some properties of the data, that is, the Dirichlet-to-Neumann map. We will formulate the direct problem as a nonlinear operator mapping from to defined as

 Fω(c−2)=Λω2c−2,

where indicates the Dirichlet to Neumann operator. Indeed, at fixed frequency , we consider the boundary value problem,

 {(−Δ−ω2c−2(x))u=0, in Ω,u=g on ∂Ω, (1)

while , where denotes the outward unit normal vector to . In this section, we will state some known results concerning the well-posedness of problem LABEL:Helmholtz (see, for example, [12]) and regularity properties of the nonlinear map . We will sketch the proofs of these results because we will need to keep track of the dependencies of the constants involved on frequency. We invoke

###### Assumption 1.

There exist two positive constants such that

 B1≤c−2≤B2in Ω. (2)

In the sequel of Section 2 indicates that depends only on the parameters and we will indicate different constants with the same letter .

###### Proposition 2.

Let be a bounded Lipschitz domain in , , and satisfying LABEL:Aprioribound1. Then, there exists a discrete set such that, for every , there exists a unique solution of

 {(−Δ−ω2c−2(x))u=f   in % Ω,u=g on ∂Ω. (3)

Furthermore, there exists a positive constant such that

 ∥u∥H1(Ω)≤C(1+ω2d(ω2,Σc−2))(∥g∥H1/2(∂Ω)+∥f∥L2(Ω)), (4)

where and indicates the distance of from .

###### Proof.

We first prove the result for . Consider the linear operators and the multiplication operator

 Mc−2:L2(Ω) →L2(Ω), (5) u →c−2u

respectively. We can now consider the operator . The equation

 (−Δ−ω2c−2(x))u=f.

for is equivalent to

 (I−ω2K)u=Δ−1f (6)

Note that is compact by Rellich–Kondrachov compactness theorem. Furthermore, by LABEL:Aprioribound1 and the properties of it follows that is self-adjoint and positive. Hence, has a discrete set of positive eigenvalues such that as . Let and define and let , and show that it satisfies the assumptions of this proposition. Then, by the Fredholm alternative, there exists a unique solution of LABEL:opeq.

To prove estimate LABEL:energy_est1 we observe that

 u=∞∑n=1⟨u,en⟩en,Ku=∞∑n=1αn⟨u,en⟩en

where is an orthonormal basis of . Hence we can rewrite LABEL:opeq in the form

 ∞∑n=1(1−ω2αn)⟨u,en⟩en=∞∑n=1⟨h,en⟩en where h=Δ−1f

Hence,

 ⟨u,en⟩=11−ω2~λn⟨h,en⟩,∀n∈N

and

 u=∞∑n=111−ω2~λn⟨h,en⟩en

so that

 ∥u∥L2(Ω)≤(1+ω2d(ω2,Σc−2))∥h∥L2(Ω)≤C(1+ω2d(ω2,Σc−2))∥f∥L2(Ω) (7)

where .

Now, by multiplying equation LABEL:Helmholtz1 with , integrating by parts, using Schwartz’ inequality, LABEL:est_u and LABEL:Aprioribound1 it follows in the case :

 ∥∇u∥L2(Ω)≤C(1+ω2d(ω2,Σc−2))∥f∥L2(Ω) (8)

Hence, by LABEL:est_grad and LABEL:est_u we finally get

 ∥u∥H1(Ω)≤C(1+ω2d(ω2,Σc−2))∥f∥L2(Ω).

If is not identically zero then we reduce the problem to the previous case by considering where is such that on and and we derive easily the estimate

 ∥u∥H1(Ω)≤C(1+ω2d(ω2,Σc−2))(∥f∥L2(Ω)+∥g∥H1/2(∂Ω))

which concludes the proof.

The constants appearing in the estimate of LABEL:2-energy depends on and which are unknown. To our purposes it would be convenient to have constants depending only on a priori parameters , and other known parameters. Let us denote by the spectrum of . Then, we have the following

###### Proposition 3.

Suppose that the assumptions of LABEL:2-energy are satisfied. Let denote the Dirichlet eigenvalues of . Then, for any ,

 λnB2≤~λn≤λnB1. (9)

If is such that,

 0<ω2<λ1B2, (10)

or, for some ,

 λnB1<ω2<λn+1B2, (11)

then there exists a unique solution of Problem LABEL:Helmholtz and the following estimate holds

 ∥u∥H1(Ω)≤C(∥g∥H1/2(∂Ω)+∥f∥L2(Ω)),

where .

###### Proof.

To derive estimate LABEL:eigenvaluebound we consider the Rayleigh quotient related to equation LABEL:Helmholtz

 ∫Ω|∇v|2∫Ωc−2v2.

By LABEL:Aprioribound1, for any non trivial we have

 1B2∫Ω|∇v|2∫Ωv2≤∫Ω|∇v|2∫Ωc−2v2≤1B1∫Ω|∇v|2∫Ωv2.

Now, we apply Courant-Rayleigh minimax principle (see for instance [10, Theorem 4.5.1], where the infinite dimensional Courant-Rayleigh minimax principle has been considered): The following arguments are similar as in the simple one-dimensional Example of Davies’ book [10, Example 4.6.1]. Due to LABEL:Aprioribound1 the Hilbert space

 L2c(Ω)={v:∫Ωc−2v2<∞},

with norm is equivalent to .

 ~λn λn :=inf{u1,⋯,un∈H10(Ω)}supv∈span{u1,⋯,un}:∥v∥L2≤1∫Ω|∇v|2∫Ωv2.

Note that implies that and that . Therefore

 λn≤inf{u1,⋯,un∈H10(Ω)}supv∈span{u1,⋯,un}:∥v∥2L2c≤B2∫Ω|∇v|2∫Ωv2.

Now, using the scale invariance of and that , we get

 λn≤B2inf{u1,⋯,un∈H10(Ω)}supv∈span{u1,⋯,un}:∥v∥L2c≤1∫Ω|∇v|2∫Ωc−2v2=B2~λn.

To get lower bound estimate for observe that if then . Hence

 ~λn≤inf{~u1,⋯,~un∈H10(Ω)}supv∈span{~u1,⋯,~un}:∥v∥2L2≤1B1∫Ω|∇v|2∫Ωc−2v2.

Now, using the scale invariance of and that , we get

 ~λn≤inf{~u1,⋯,~un∈H10(Ω)}supv∈span{~u1,⋯,~un}:∥v∥L2≤1∫Ω|∇v|2∫Ωc−2v2=1B1λn.

Thus we have shown that

 λnB2≤~λn≤λnB1,∀n∈N.

Hence, we have well-posedness of problem LABEL:Helmholtz if we select an satisfying LABEL:smallfr or LABEL:higherfr and the claim follows.

We observe that in order to derive the uniform estimates of LABEL:uniformconstants we need to assume that either the frequency is small LABEL:smallfr or that the oscillation of is sufficiently small LABEL:higherfr. This observation can also been found in Davies’ book [10].

In the seismic application we have in mind we might know the spectrum of some reference wavespeed . The following local result holds

###### Proposition 4.

Let and satisfy the assumptions of LABEL:2-energy and let where is the Dirichlet spectrum of equation LABEL:Helmholtz corresponding to . Then, there exists such that, if

 ∥c−2−c−20∥L∞(Ω)≤δ,

then and the solution of Problem LABEL:Helmholtz1 corresponding to satisfies

 ∥u∥H1(Ω)≤C⎛⎝1+ω2d(ω2,Σc−20)⎞⎠(∥f∥L2(Ω)+∥g∥H1/2(∂Ω)),

.

###### Proof.

Let and consider the unique solution of LABEL:Helmholtz1 for and consider the problem

 {−Δv−ω2c−20v−ω2δcv=ω2u0δcin Ω,v=0,on ∂Ω. (12)

Let now

 L:=−Δ−ω2c−20

then, by assumption, it is invertible from to and we can rewrite problem LABEL:auxiliarypr in the form

 (I−K)v=h, (13)

where and is the multiplication operator defined in LABEL:eq:mult and . Observe now that from LABEL:energy_est1 with and where . Hence, we derive

 ∥K∥≤ω2∥L−1∥∥Mδc∥≤ω2∥L−1∥δ≤Cω2(1+ω2d0)δ.

Hence, choosing the bounded operator has norm smaller than one. Hence, is invertible and there exists a unique solution of LABEL:auxiliaryeq in satisfying LABEL:energy_est1 with and since the statement follows.

Let be such that either

 0<ω2<λ1B2,

or for some

 λnB1<ω2<λn+1B2,

and let

 W:={c−2∈L∞(Ω):B1≤c−2≤B2}.

Then the direct operator

 Fω:W →L(H1/2(∂Ω),H−1/2(∂Ω)), c−2 ↦Λω2c−2,

is well defined.

We will examine regularity properties of in the following lemmas. We will show the Fréchet differentiability of it.

###### Lemma 5 (Fréchet differentiability).

Let satisfy LABEL:Aprioribound1. Assume that . Then, the direct operator is Fréchet differentiable at and its Fréchet derivative satisfies

 ∥DFω[c−2]∥L(L∞(Ω),L(H1/2(∂Ω),H−1/2(∂Ω)))≤Cω2(1+ω2d(ω2,Σc−2))2 (14)

where .

###### Proof.

Consider . Then, from LABEL:continuity, if is small enough, . An application of Alessandrini’s identity then gives

 ⟨(Λω2(c−2+δc−2)−Λω2c−2)g,h⟩=ω2∫Ωδc−2uvdx, (15)

where where is the dual pairing with respect to and and and solve the boundary value problems,

 {(−Δ−ω2(c−2+δc−2))u=0,x∈Ω,u=g,x∈∂Ω,

and

 {(−Δ−ω2c−2)v=0,x∈Ω,v=h,x∈∂Ω,

respectively. We first show that the map is Fréchet differentiable and that the Fréchet derivative is given by

 ⟨DFω[c−2](δc−2)g,h⟩=ω2∫Ωδc−2~uvdx, (16)

where solves the equation

 {(−Δ−ω2c−2)~u=0,x∈Ω,~u=g,x∈∂Ω.

In fact, by LABEL:ident-1, we have that

 ⟨(Λω2(c−2+δc−2)−Λω2c−2)g,h⟩−ω2∫Ωδc−2~uvdx=ω2∫Ωδc−2(u−~u)vdx. (17)

We note that solves the equations

 {(−Δ−ω2c−2)(u−~u)=−ω2δc−2u,x∈Ω,u−~u=0,x∈∂Ω.

Using the fact that and are in and that and applying Cauchy-Schwarz inequality, we get

 ∣∣∣ω2∫Ωδc−2(u−~u)vdx∣∣∣≤ω2∥δc−2∥L∞(Ω)∥u−~u∥L2(Ω)∥v∥L2(Ω). (18)

Finally, using the stability estimates of LABEL:2-energy applied to and to and the stability estimates of LABEL:continuity applied to we derive

 ∣∣∣ω2∫Ωδc−2(u−~u)vdx∣∣∣≤Cω4(1+ω2d(ω2,Σc−2))3∥δc−2∥2L∞(Ω)∥g∥H1/2(∂Ω)∥h∥H1/2(∂Ω). (19)

Hence

 ∣∣∣⟨(Λω2(c−2+δc−2)−Λω2c−2)g,h⟩−ω2∫Ωδc−2~uvdx∣∣∣ ≤ Cω4(1+ω2d(ω2,Σc−2))3∥δc−2∥2L∞(Ω)∥g∥H1/2(∂Ω)∥h∥H1/2(∂Ω),

which proves differentiability.
Finally by

 ⟨DFω[c−2](δc−2)g,h⟩=ω2∫Ωδc−2~uvdx,

and we get

 ∣∣⟨DFω[c−2](δc−2)g,h⟩∣∣≤ ω2∥δc−2∥L∞(Ω)∥~u∥L2(Ω)∥v∥L2(Ω) ≤ ω2(1+ω2d(ω2,Σc−2))2∥δc−2∥L∞(Ω)∥g∥H1/2(∂Ω)∥h∥H1/2(∂Ω).

from which LABEL:L-bd follows.

### 2.2 Conditional quantitative Lipschitz stability estimate

Let be positive with , . In the sequel we will refer to these numbers as to the a priori data. To prove the results of this section we invoke the following common assumptions

###### Assumption 6.

is a bounded domain such that

 |x|≤Ar1,∀x∈Ω.

Moreover,

 ∂Ω of Lipschitz class with constants r1 and L.

Let be a partition of given by

 DN≜{{D1,D2,…,DN}∣N⋃j=1¯¯¯¯¯Dj=Ω,(Dj∩Dj′)∘=∅,j≠j′} (20)

such that

 {∂Dj}Nj=1 is of Lipschitz class with constants r0 and L.

###### Assumption 7.

The function , that is, it satisfies

 B1≤c−2≤B2,in Ω

and is of the form

 c−2(x)=N∑j=1cjχDj(x),

where are unknown numbers and .

###### Assumption 8.

Assume

 0<ω2<λ1B2,

or, for some ,

 λnB1<ω2<λn+1B2.

Under the above assumptions we can state the following preliminary result

###### Lemma 9.

Let and satisfy LABEL:assumption_domain and let . Then, for every , there exists a positive constant with such that

 ∥c−2∥Hs′(Ω)≤C(L,s′)1rs′0∥c−2∥L2(Ω). (21)

###### Proof.

The proof is based on the extension of a result of Magnanini and Papi in [16] to the three dimensional setting. In fact, following the argument in [16], one has that

 ∥χDj∥2Hs′(Ω)≤16π(1−2s′)(2s′)1+2s′|Dj|1−2s′|∂Dj|2s′. (22)

We now use the fact that is a partition of disjoint sets of to show the following inequality

 ∥c−2∥2Hs′(Ω)≤2N∑j=1c2j∥χDj∥2Hs′(Ω) (23)

In fact, in order to prove LABEL:Hsbound2 recall that

 ∥c−2∥2Hs′(Ω)=∫Ω∫Ω|∑Nj=1cj(χDj(x)−χDj(y))|2|x−y|3+2s′dxdy

and observe that, since the is a partition of disjoint sets of , we get

 |N∑j=1cj(χDj(x)−χDj(y))|2=N∑j=1c2j(χDj(x)−χDj(y))2−∑i≠jcicjχDi(x)χDj(y)

Again, by the fact that the are disjoint sets, we have

 ∑i≠j|cicj|χDi(x)χDj(y)≤∑i≠jc2i+c2j2χDi(x)χDj(y) = ∑i≠jc2i2(χDi(x)−χDi(y))2χDi(x)χDj(y)+∑i≠jc2j2(χDj(x)−χDj(y))2χDi(x)χDj(y) ≤ ∑i≠jc2i2(χDi(x)−χDi(y))2χDj(y)+∑i≠jc2j2(χDj(x)−χDj(y))2χDi(x) ≤ N∑i=1c2i2(χDi(x)−χDi(y))2N∑j=1χDj(y)+N∑j=1c2j2(χDj(x)−χDj(y))2N∑i=1χDi(y) ≤ N∑i=1c2i2(χDi(x)−χDi(y))2+N∑j=1c2j2(χDj(x)−χDj(y))2 = N∑i=1c2i(χDi(x)−χDi(y))2

where we have used the fact that . So, we have derived that

 |N∑j=1cj(χDj(x)−χDj(y))|2≤2N∑j=1c2j(χDj(x)−χDj(y))2

from which it follows that

 ∥c−2∥2Hs′(Ω) = ∫Ω∫Ω|∑Nj=1cj(χDj(x)−χDj(y))|2|x−y|3+2s′dxdy ≤ 2∫Ω∫Ω∑Nj=1c2j(χDj(x)−χDj(y))2|x−y|3+2s′dxdy ≤ 2N∑j=1c2j∫Ω∫Ω(χDj(x)−χDj(y))2|x−y|3+2s′dxdy=2N∑j=1c2j∥χDj∥2Hs′(Ω)

which proves LABEL:Hsbound2. so that finally from LABEL:assumption_domain, LABEL:bound1 and LABEL:Hsbound2 we get

 ∥c−2∥2Hs′(Ω)≤2N∑j=1c2j∥χDj∥2Hs′(Ω)≤C(s′)N∑j=1c2j|Dj|(|∂Dj||Dj|)2s′