# Inverse Boundary Value Problem For The Helmholtz Equation: Quantitative Conditional Lipschitz Stability Estimates

Elena Beretta, Maarten V. de Hoop, Florian Faucher, Otmar Scherzer

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

## To cite this version:

Elena Beretta, Maarten V. de Hoop, Florian Faucher, Otmar Scherzer. Inverse Boundary Value Problem For The Helmholtz Equation: Quantitative Conditional Lipschitz Stability Estimates. SIAM Journal on Mathematical Analysis, 2016, 48 (6), pp.3962 - 3983. 10.1137/15M1043856. hal-01407446

HAL Id: hal-01407446

https://hal.science/hal-01407446v1

Submitted on 6 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# INVERSE BOUNDARY VALUE PROBLEM FOR THE HELMHOLTZ EQUATION: QUANTITATIVE CONDITIONAL LIPSCHITZ STABILITY ESTIMATES *

ELENA BERETTA†, MAARTEN V. DE HOOP‡, FLORIAN FAUCHER§, AND OTMAR SCHERZER¶

**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 $L^2(\Omega)$.

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

<footnote>*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.</footnote>

<footnote>†Dipartimento di Matematica “Brioschi”, Politecnico di Milano, Italy (elena.beretta@polimi.it).</footnote>

<footnote>‡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.</footnote>

<footnote>§INRIA Bordeaux Sud-Ouest Research Center, Team Project Magique-3D, France (florian.faucher@inria.fr).</footnote>

<footnote>¶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.</footnote>
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 $F_\omega$ from $L^\infty(\Omega)$ to $\mathcal{L}(H^{1/2}(\partial\Omega), H^{-1/2}(\partial\Omega))$ defined as

$$F_\omega(c^{-2}) = \Lambda_{\omega^2 c^{-2}},$$

where $\Lambda_{\omega^2 c^{-2}}$ indicates the Dirichlet to Neumann operator. Indeed, at fixed frequency $\omega^2$, we consider the boundary value problem,

$$ (1) \qquad \begin{cases} (-\Delta - \omega^2 c^{-2}(x))u = 0, & \text{in } \Omega, \\ u = g & \text{on } \partial\Omega, \end{cases} $$

while $\Lambda_{\omega^2 c^{-2}}: g \to \frac{\partial u}{\partial \nu}|_{\partial \Omega}$, where $\nu$ denotes the outward unit normal vector to $\partial \Omega$. In this section, we will state some known results concerning the well-posedness of problem (1) (see, for example, [12]) and regularity properties of the nonlinear map $F_\omega$. 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 $B_1, B_2$ such that

$$ (2) \qquad B_1 \le c^{-2} \le B_2 \quad \text{in } \Omega. $$

In the sequel of Section 2 $C = C(a, b, c, \dots)$ indicates that $C$ depends only on the parameters $a, b, c, \dots$ and we will indicate different constants with the same letter $C$.

**PROPOSITION 2.** Let $\Omega$ be a bounded Lipschitz domain in $\mathbb{R}^3$, $f \in L^2(\Omega)$, $g \in H^{1/2}(\partial\Omega)$ and $c^{-2} \in L^\infty(\Omega)$ satisfying Assumption 1. Then, there exists a discrete set $\Sigma_{c^{-2}} := \{\tilde{\lambda}_n \mid \lambda_n > 0, \forall n \in \mathbb{N}\}$ such that, for every $\omega^2 \in \mathbb{C}\setminus\Sigma_{c^{-2}}$, there exists a unique solution $u \in H^1(\Omega)$ of

$$ (3) \qquad \begin{cases} (-\Delta - \omega^2 c^{-2}(x))u = f & \text{in } \Omega, \\ u = g & \text{on } \partial\Omega. \end{cases} $$
Furthermore, there exists a positive constant C such that

$$ (4) \qquad \|u\|_{H^1(\Omega)} \le C \left( 1 + \frac{\omega^2}{d(\omega^2, \Sigma_{c-2})} \right) \left( \|g\|_{H^{1/2}(\partial\Omega)} + \|f\|_{L^2(\Omega)} \right), $$

where $C = C(\Omega, B_2)$ and $d(\omega^2, \Sigma_{c-2})$ indicates the distance of $\omega^2$ from $\Sigma_{c-2}$.

*Proof.* We first prove the result for $g = 0$. Consider the linear operators $-\Delta : H_0^1(\Omega) \to H^{-1}(\Omega)$ and the multiplication operator

$$ (5) \qquad M_{c-2} : L^2(\Omega) \to L^2(\Omega), \\ u \to c^{-2}u $$

respectively. We can now consider the operator $K = \Delta^{-1}M_{c-2} : H_0^1(\Omega) \to H_0^1(\Omega)$. The equation

$$ (-\Delta - \omega^2 c^{-2}(x))u = f. $$

for $u \in H_0^1(\Omega)$ is equivalent to

$$ (6) \qquad (I - \omega^2 K)u = \Delta^{-1}f $$

Note that $K : H_0^1(\Omega) \to H_0^1(\Omega)$ is compact by Rellich-Kondrachov compactness theorem. Furthermore, by Assumption 1 and the properties of $\Delta^{-1}$ it follows that $K$ is self-adjoint and positive. Hence, $K$ has a discrete set of positive eigenvalues $\{\alpha_n\}_{n \in \mathbb{N}}$ such that $\alpha_n \to 0$ as $n \to \infty$. Let $\tilde{\lambda}_n := \frac{1}{\alpha_n}$, $n \in \mathbb{N}$ and define $\Sigma_{c-2} := \{\tilde{\lambda}_n : n \in \mathbb{N}\}$ and let $\omega^2 \in \mathbb{C} \setminus \Sigma_{c-2}$, and show that it satisfies the assumptions of this proposition. Then, by the Fredholm alternative, there exists a unique solution $u \in H_0^1(\Omega)$ of (6).

To prove estimate (4) we observe that

$$ u = \sum_{n=1}^{\infty} \langle u, e_n \rangle e_n, \quad Ku = \sum_{n=1}^{\infty} \alpha_n \langle u, e_n \rangle e_n $$

where $\{e_n\}_{n \in \mathbb{N}}$ is an orthonormal basis of $L^2(\Omega)$. Hence we can rewrite (6) in the form

$$ \sum_{n=1}^{\infty} (1 - \omega^2 \alpha_n) \langle u, e_n \rangle e_n = \sum_{n=1}^{\infty} \langle h, e_n \rangle e_n \quad \text{where } h = \Delta^{-1}f $$

Hence,

$$ \langle u, e_n \rangle = \frac{1}{1 - \frac{\omega^2}{\tilde{\lambda}_n}} \langle h, e_n \rangle, \quad \forall n \in \mathbb{N} $$

and

$$ u = \sum_{n=1}^{\infty} \frac{1}{1 - \frac{\omega^2}{\tilde{\lambda}_n}} \langle h, e_n \rangle e_n $$

so that

$$ (7) \qquad \|u\|_{L^2(\Omega)} \le \left( 1 + \frac{\omega^2}{d(\omega^2, \Sigma_{c-2})} \right) \|h\|_{L^2(\Omega)} \le C \left( 1 + \frac{\omega^2}{d(\omega^2, \Sigma_{c-2})} \right) \|f\|_{L^2(\Omega)} $$

where $C = C(\Omega, B_2)$.

Now, by multiplying equation (3) with $u$, integrating by parts, using Schwartz' inequality, Assumption 1 and (7) it follows in the case $g = 0$:

$$ (8) \qquad \|\nabla u\|_{L^2(\Omega)} \le C \left( 1 + \frac{\omega^2}{d(\omega^2, \Sigma_{c-2})} \right) \|f\|_{L^2(\Omega)} $$
Hence, by (7) and (8) we finally get

$$ \|u\|_{H^1(\Omega)} \le C \left( 1 + \frac{\omega^2}{d(\omega^2, \Sigma_{c-2})} \right) \|f\|_{L^2(\Omega)}. $$

If $g$ is not identically zero then we reduce the problem to the previous case by considering $v = u - \tilde{g}$ where $\tilde{g} \in H^1(\Omega)$ is such that $\tilde{g} = g$ on $\partial\Omega$ and $\|\tilde{g}\|_{H^1(\Omega)} \le \|g\|_{H^{1/2}(\partial\Omega)}$ and we derive easily the estimate

$$ \|u\|_{H^1(\Omega)} \le C \left( 1 + \frac{\omega^2}{d(\omega^2, \Sigma_{c-2})} \right) \left( \|f\|_{L^2(\Omega)} + \|g\|_{H^{1/2}(\partial\Omega)} \right) $$

which concludes the proof.

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

The constants appearing in the estimate of Proposition 2 depends on $c^{-2}$ and $\Sigma_{c-2}$ which are unknown. To our purposes it would be convenient to have constants depending only on a priori parameters $B_1$, $B_2$ and other known parameters. Let us denote by $\Sigma_0$ the spectrum of $-\Delta$. Then, we have the following

**PROPOSITION 3.** Suppose that the assumptions of Proposition 2 are satisfied. Let $\{\lambda_n\}_{n \in \mathbb{N}}$ denote the Dirichlet eigenvalues of $-\Delta$. Then, for any $n \in \mathbb{N}$,

$$ (9) \qquad \frac{\lambda_n}{B_2} \le \tilde{\lambda}_n \le \frac{\lambda_n}{B_1}. $$

If $\omega^2$ is such that,

$$ (10) \qquad 0 < \omega^2 < \frac{\lambda_1}{B_2}, $$

or, for some $n \ge 1$,

$$ (11) \qquad \frac{\lambda_n}{B_1} < \omega^2 < \frac{\lambda_{n+1}}{B_2}, $$

then there exists a unique solution $u \in H^1(\Omega)$ of Problem (1) and the following estimate holds

$$ \|u\|_{H^1(\Omega)} \le C \left( \|g\|_{H^{1/2}(\partial\Omega)} + \|f\|_{L^2(\Omega)} \right), $$

where $C = C(B_1, B_2, \omega^2, \Sigma_0)$.

*Proof.* To derive estimate (9) we consider the Rayleigh quotient related to equation (1)

$$ \frac{\int_{\Omega} |\nabla v|^2}{\int_{\Omega} c^{-2} v^2}. $$

By Assumption 1, for any non trivial $v \in H_0^1(\Omega)$ we have

$$ \frac{1}{B_2} \frac{\int_{\Omega} |\nabla v|^2}{\int_{\Omega} v^2} \le \frac{\int_{\Omega} |\nabla v|^2}{\int_{\Omega} c^{-2} v^2} \le \frac{1}{B_1} \frac{\int_{\Omega} |\nabla v|^2}{\int_{\Omega} v^2}. $$

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 Assumption 1 the Hilbert space

$$ L_c^2(\Omega) = \{v : \int_{\Omega} c^{-2} v^2 < \infty\}, $$
with norm $||v||_{L^2_c} = \int_{\Omega} v^2 c^{-2}$ is equivalent to $L^2(\Omega)$.

$$ \tilde{\lambda}_n := \inf_{\{\tilde{u}_1, \dots, \tilde{u}_n \in H_0^1(\Omega)\}} \sup_{v \in \operatorname{span}\{\tilde{u}_1, \dots, \tilde{u}_n\}: ||v||_{L^2_c} \le 1} \frac{\int_{\Omega} |\nabla v|^2}{\int_{\Omega} c^{-2} v^2}, $$

$$ \lambda_n := \inf_{\{u_1, \dots, u_n \in H_0^1(\Omega)\}} \sup_{v \in \operatorname{span}\{u_1, \dots, u_n\}: ||v||_{L^2} \le 1} \frac{\int_{\Omega} |\nabla v|^2}{\int_{\Omega} v^2}. $$

Note that $||v||_{L^2} \le 1$ implies that $||v||_{L^2_c}^2 \le B_2$ and that $L_c^2(\Omega) = L^2(\Omega)$. Therefore

$$ \lambda_n \le \inf_{\{u_1, \dots, u_n \in H_0^1(\Omega)\}} \sup_{v \in \operatorname{span}\{u_1, \dots, u_n\}: ||v||_{L^2_c}^2 \le B_2} \frac{\int_{\Omega} |\nabla v|^2}{\int_{\Omega} v^2}. $$

Now, using the scale invariance of $\frac{\int_{\Omega} |\nabla v|^2}{\int_{\Omega} v^2}$ and that $c^{-2} \le B_2$, we get

$$ \lambda_n \le B_2 \inf_{\{u_1, \dots, u_n \in H_0^1(\Omega)\}} \sup_{v \in \operatorname{span}\{u_1, \dots, u_n\}: ||v||_{L^2_c} \le 1} \frac{\int_{\Omega} |\nabla v|^2}{\int_{\Omega} c^{-2} v^2} = B_2 \tilde{\lambda}_n. $$

To get lower bound estimate for $\tilde{\lambda}_n$ observe that if $||v||_{L^2_c} \le 1$ then $||v||_{L^2}^2 \le \frac{1}{B_1}$. Hence

$$ \tilde{\lambda}_n \le \inf_{\{\tilde{u}_1, \dots, \tilde{u}_n \in H_0^1(\Omega)\}} \sup_{v \in \operatorname{span}\{\tilde{u}_1, \dots, \tilde{u}_n\}: ||v||_{L^2}^2 \le \frac{1}{B_1}} \frac{\int_{\Omega} |\nabla v|^2}{\int_{\Omega} c^{-2} v^2}. $$

Now, using the scale invariance of $\frac{\int_{\Omega} |\nabla v|^2}{\int_{\Omega} v^2}$ and that $c^{-2} \ge B_1$, we get

$$ \tilde{\lambda}_n \le \inf_{\{\tilde{u}_1, \dots, \tilde{u}_n \in H_0^1(\Omega)\}} \sup_{v \in \operatorname{span}\{\tilde{u}_1, \dots, \tilde{u}_n\}: ||v||_{L^2} \le 1} \frac{\int_{\Omega} |\nabla v|^2}{\int_{\Omega} c^{-2} v^2} = \frac{1}{B_1} \lambda_n. $$

Thus we have shown that

$$ \frac{\lambda_n}{B_2} \le \tilde{\lambda}_n \le \frac{\lambda_n}{B_1}, \quad \forall n \in \mathbb{N}. $$

Hence, we have well-posedness of problem (1) if we select an $\omega^2$ satisfying (10) or (11) and the claim follows. □

We observe that in order to derive the uniform estimates of Proposition 3 we need to assume that either the frequency is small (10) or that the oscillation of $c^{-2}$ is sufficiently small (11). 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 $c_0^{-2}$. The following local result holds

**PROPOSITION 4.** Let $\Omega$ and $c_0^{-2}$ satisfy the assumptions of Proposition 2 and let $\omega^2 \in \mathbb{C} \setminus \Sigma_{c_0^{-2}}$ where $\Sigma_{c_0^{-2}}$ is the Dirichlet spectrum of equation (1) corresponding to $c_0^{-2}$. Then, there exists $\delta = \delta(\Omega, \omega^2, B_2, \Sigma_{c_0^{-2}}) > 0$ such that, if

$$ \|c^{-2} - c_0^{-2}\|_{L^{\infty}(\Omega)} \le \delta, $$

then $\omega^2 \in \mathbb{C} \setminus \Sigma_{c^{-2}}$ and the solution $u$ of Problem (3) corresponding to $c^{-2}$ satisfies

$$ \|u\|_{H^1(\Omega)} \le C \left( 1 + \frac{\omega^2}{d(\omega^2, \Sigma_{c_0^{-2}})} \right) \left( \|f\|_{L^2(\Omega)} + \|g\|_{H^{1/2}(\partial\Omega)} \right), $$

$$ C = C(\Omega, B_2). $$
Proof. Let $\delta_c := c^{-2} - c_0^{-2}$ and consider $u_0 \in H^1(\Omega)$ the unique solution of (3) for $c_0^{-2}$ and consider the problem

$$ (12) \qquad \begin{cases} -\Delta v - \omega^2 c_0^{-2} v - \omega^2 \delta_c v = \omega^2 u_0 \delta_c & \text{in } \Omega, \\ v = 0, & \text{on } \partial\Omega. \end{cases} $$

Let now

$$ L := -\Delta - \omega^2 c_0^{-2} $$

then, by assumption, it is invertible from $H_0^1(\Omega)$ to $L^2(\Omega)$ and we can rewrite problem (12) in the form

$$ (13) \qquad (I - K)v = h, $$

where $K = \omega^2 L^{-1} M_{\delta_c}$ and $M_{\delta_c}$ is the multiplication operator defined in (5) and $h = L^{-1}(\omega^2 u_0 \delta_c)$. Observe now that from (4) $\|L^{-1}\| \le C(1 + \frac{\omega^2}{d_0})$ with $C = C(\Omega, B_2)$ and where $d_0 = \text{dist}(\omega^2, \Sigma_{c_0^{-2}})$. Hence, we derive

$$ \|K\| \le \omega^2 \|L^{-1}\| \|M_{\delta_c}\| \le \omega^2 \|L^{-1}\| \delta \le C \omega^2 \left(1 + \frac{\omega^2}{d_0}\right) \delta. $$

Hence, choosing $\delta = \frac{1}{2}(C\omega^2(1 + \frac{\omega^2}{d_0}))^{-1}$ the bounded operator $K$ has norm smaller than one. Hence, $I - K$ is invertible and there exists a unique solution $v$ of (13) in $H_0^1$ satisfying (4) with $C = C(B_2, \omega^2, \Omega, d_0)$ and since $u = u_0 + v$ the statement follows. $\square$

Let $\omega^2$ be such that either

$$ 0 < \omega^2 < \frac{\lambda_1}{B_2}, $$

or for some $n \ge 1$

$$ \frac{\lambda_n}{B_1} < \omega^2 < \frac{\lambda_{n+1}}{B_2}, $$

and let

$$ W := \{c^{-2} \in L^{\infty}(\Omega) : B_1 \le c^{-2} \le B_2\}. $$

Then the direct operator

$$ F_{\omega}: W \to \mathcal{L}(H^{1/2}(\partial\Omega), H^{-1/2}(\partial\Omega)), \\ c^{-2} \mapsto \Lambda_{\omega^2 c^{-2}}, $$

is well defined.

We will examine regularity properties of $F_{\omega}$ in the following lemmas. We will show the Fréchet differentiability of it.

**LEMMA 5 (Fréchet differentiability).** Let $c^{-2} \in L^{\infty}(\Omega)$ satisfy Assumption 1. Assume that $\omega^2 \in \mathbb{C} \setminus \Sigma_{c^{-2}}$. Then, the direct operator $F_{\omega}$ is Fréchet differentiable at $c^{-2}$ and its Fréchet derivative $DF_{\omega}(c^{-2})$ satisfies

$$ (14) \qquad \|DF_{\omega}[c^{-2}]\|_{L(L^{\infty}(\Omega), L(H^{1/2}(\partial\Omega), H^{-1/2}(\partial\Omega)))} \le C\omega^2 \left(1 + \frac{\omega^2}{d(\omega^2, \Sigma_{c^{-2}})}\right)^2 $$

where $C = C(\Omega, B_2)$.
Proof. Consider $c^{-2} + \delta c^{-2}$. Then, from Proposition 4, if $\|\delta c^{-2}\|_{L^{\infty}(\Omega)}$ is small enough, $\omega^2 \notin \Sigma_{c^{-2}+\delta c^{-2}}$. An application of Alessandrini's identity then gives

$$ (15) \qquad \langle (\Lambda_{\omega^2(c^{-2}+\delta c^{-2})} - \Lambda_{\omega^2 c^{-2}})g, h \rangle = \omega^2 \int_{\Omega} \delta c^{-2} uv \, dx, $$

where $\langle \cdot, \cdot \rangle$ is the dual pairing with respect to $H^{-1/2}(\partial\Omega)$ and $H^{1/2}(\partial\Omega)$ and $u$ and $v$ solve the boundary value problems,

$$ \left\{ \begin{array}{ll} (-\Delta - \omega^2(c^{-2} + \delta c^{-2}))u = & 0, \quad x \in \Omega, \\ u = & g, \quad x \in \partial\Omega, \end{array} \right. $$

and

$$ \left\{ \begin{array}{ll} (-\Delta - \omega^2 c^{-2})v = & 0, \quad x \in \Omega, \\ v = & h, \quad x \in \partial\Omega, \end{array} \right. $$

respectively. We first show that the map $F_{\omega}$ is Fréchet differentiable and that the Fréchet derivative is given by

$$ (16) \qquad \langle DF_{\omega}[c^{-2}](\delta c^{-2})g, h \rangle = \omega^2 \int_{\Omega} \delta c^{-2} \tilde{u} v \, dx, $$

where $\tilde{u}$ solves the equation

$$ \left\{ \begin{array}{ll} (-\Delta - \omega^2 c^{-2})\tilde{u} = & 0, \quad x \in \Omega, \\ \tilde{u} = & g, \quad x \in \partial\Omega. \end{array} \right. $$

In fact, by (15), we have that

$$ (17) \qquad \langle (\Lambda_{\omega^2(c^{-2}+\delta c^{-2})} - \Lambda_{\omega^2 c^{-2}})g, h \rangle - \omega^2 \int_{\Omega} \delta c^{-2} \tilde{u} v \, dx = \omega^2 \int_{\Omega} \delta c^{-2} (u - \tilde{u}) v \, dx. $$

We note that $u - \tilde{u}$ solves the equations

$$ \left\{ \begin{array}{ll} (-\Delta - \omega^2 c^{-2})(u - \tilde{u}) = & -\omega^2 \delta c^{-2} u, \quad x \in \Omega, \\ u - \tilde{u} = & 0, \quad x \in \partial\Omega. \end{array} \right. $$

Using the fact that $u - \tilde{u}$ and $v$ are in $H^1(\Omega)$ and that $\delta c^{-2} \in L^\infty(\Omega)$ and applying Cauchy-Schwarz inequality, we get

$$ (18) \qquad \left| \omega^2 \int_{\Omega} \delta c^{-2} (u - \tilde{u}) v \, dx \right| \le \omega^2 \|\delta c^{-2}\|_{L^{\infty}(\Omega)} \|u - \tilde{u}\|_{L^2(\Omega)} \|v\|_{L^2(\Omega)}. $$

Finally, using the stability estimates of Proposition 2 applied to $u - \tilde{u}$ and to $v$ and the stability estimates of Proposition 4 applied to $u$ we derive

$$ (19) \qquad \left| \omega^2 \int_{\Omega} \delta c^{-2} (u - \tilde{u}) v \, dx \right| \le C \omega^4 \left( 1 + \frac{\omega^2}{d(\omega^2, \Sigma_{c^{-2}})} \right)^3 \|\delta c^{-2}\|_{L^{\infty}(\Omega)}^2 \|g\|_{H^{1/2}(\partial\Omega)} \|h\|_{H^{1/2}(\partial\Omega)}. $$

Hence

$$ \left| \langle (\Lambda_{\omega^2(c^{-2}+\delta c^{-2})} - \Lambda_{\omega^2 c^{-2}})g, h \rangle - \omega^2 \int_{\Omega} \delta c^{-2} \tilde{u} v \, dx \right| \le C \omega^4 \left( 1 + \frac{\omega^2}{d(\omega^2, \Sigma_{c^{-2}})} \right)^3 \|\delta c^{-2}\|_{L^{\infty}(\Omega)}^2 \|g\|_{H^{1/2}(\partial\Omega)} \|h\|_{H^{1/2}(\partial\Omega)}, $$

which proves differentiability.

Finally by

$$ \langle DF_{\omega}[c^{-2}](\delta c^{-2})g, h \rangle = \omega^2 \int_{\Omega} \delta c^{-2} \tilde{u} v \, dx, $$
and we get

$$ \begin{aligned} |\langle DF_{\omega}[c^{-2}](\delta c^{-2})g, h \rangle| & \le \omega^2 \|\delta c^{-2}\|_{L^{\infty}(\Omega)} \|\tilde{u}\|_{L^2(\Omega)} \|v\|_{L^2(\Omega)} \\ & \le \omega^2 \left(1 + \frac{\omega^2}{d(\omega^2, \Sigma_{c-2})}\right)^2 \|\delta c^{-2}\|_{L^{\infty}(\Omega)} \|g\|_{H^{1/2}(\partial\Omega)} \|h\|_{H^{1/2}(\partial\Omega)}. \end{aligned} $$

from which (14) follows. □

**2.2. Conditional quantitative Lipschitz stability estimate.** Let $B_2, r_0, r_1, A, L, N$ be positive with $N \in \mathbb{N}, N \ge 2, r_0 < 1$. 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. $\Omega \subset \mathbb{R}^3$ is a bounded domain such that

$$ |x| \le Ar_1, \quad \forall x \in \Omega. $$

Moreover,

$\partial\Omega$ of Lipschitz class with constants $r_1$ and $L$.

Let $\mathcal{D}_N$ be a partition of $\Omega$ given by

$$ (20) \qquad \mathcal{D}_N \triangleq \left\{ \{D_1, D_2, \dots, D_N\} \mid \bigcup_{j=1}^N \bar{D}_j = \Omega, (D_j \cap D_{j'})^\circ = \emptyset, j \ne j' \right\} $$

such that

$\{\partial D_j\}_{j=1}^N$ is of Lipschitz class with constants $r_0$ and $L$.

ASSUMPTION 7. The function $c^{-2} \in W_N$, that is, it satisfies

$$ B_1 \le c^{-2} \le B_2, \quad \text{in } \Omega $$

and is of the form

$$ c^{-2}(x) = \sum_{j=1}^{N} c_j \chi_{D_j}(x), $$

where $c_j, j = 1, \dots, N$ are unknown numbers and $\{D_1, \dots, D_N\} \in \mathcal{D}_N$.

ASSUMPTION 8. Assume

$$ 0 < \omega^2 < \frac{\lambda_1}{B_2}, $$

or, for some $n \ge 1$,

$$ \frac{\lambda_n}{B_1} < \omega^2 < \frac{\lambda_{n+1}}{B_2}. $$

Under the above assumptions we can state the following preliminary result

LEMMA 9. Let $\Omega$ and $\mathcal{D}_N$ satisfy Assumption 6 and let $c^{-2} \in W_N$. Then, for every $s' \in (0, 1/2)$, there exists a positive constant $C$ with $C = C(L, s')$ such that

$$ (21) \qquad \|c^{-2}\|_{H^{s'}(\Omega)} \le C(L, s') \frac{1}{r_0^{s'}} \|c^{-2}\|_{L^2(\Omega)}. $$
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

$$ (22) \qquad \|\chi_{D_j}\|_{H^{s'}(\Omega)}^2 \le \frac{16\pi}{(1-2s')(2s')^{1+2s'}} |D_j|^{1-2s'} |\partial D_j|^{2s'}. $$

We now use the fact that {$D_j$}$_{j=1}^N$ is a partition of disjoint sets of $\Omega$ to show the following inequality

$$ (23) \qquad \|c^{-2}\|_{H^{s'}(\Omega)}^2 \le 2 \sum_{j=1}^{N} c_j^2 \|\chi_{D_j}\|_{H^{s'}(\Omega)}^2 $$

In fact, in order to prove (23) recall that

$$ \|c^{-2}\|_{H^{s'}(\Omega)}^2 = \int_{\Omega} \int_{\Omega} \frac{|\sum_{j=1}^{N} c_j(\chi_{D_j}(x) - \chi_{D_j}(y))|^2}{|x-y|^{3+2s'}} dx dy $$

and observe that, since the {$D_j$}$_{j=1}^N$ is a partition of disjoint sets of $\Omega$, we get

$$ \left| \sum_{j=1}^{N} c_j (\chi_{D_j}(x) - \chi_{D_j}(y)) \right|^2 = \sum_{j=1}^{N} c_j^2 (\chi_{D_j}(x) - \chi_{D_j}(y))^2 - \sum_{i \neq j} c_i c_j \chi_{D_i}(x) \chi_{D_j}(y) $$

Again, by the fact that the {$D_j$}$_{j=1}^N$ are disjoint sets, we have

$$ \begin{align*} \sum_{i \neq j} |c_i c_j |\chi_{D_i}(x)\chi_{D_j}(y)| &\le \sum_{i \neq j} \frac{c_i^2 + c_j^2}{2} |\chi_{D_i}(x)\chi_{D_j}(y)| \\ &= \sum_{i \neq j} \frac{c_i^2}{2} (\chi_{D_i}(x) - \chi_{D_i}(y))^2 \chi_{D_i}(x)\chi_{D_j}(y) + \sum_{i \neq j} \frac{c_j^2}{2} (\chi_{D_j}(x) - \chi_{D_j}(y))^2 \chi_{D_i}(x)\chi_{D_j}(y) \\ &\le \sum_{i \neq j} \frac{c_i^2}{2} (\chi_{D_i}(x) - \chi_{D_i}(y))^2 \chi_{D_j}(y) + \sum_{i \neq j} \frac{c_j^2}{2} (\chi_{D_j}(x) - \chi_{D_j}(y))^2 \chi_{D_i}(x) \\ &\le \sum_{i=1}^{N} \frac{c_i^2}{2} (\chi_{D_i}(x) - \chi_{D_i}(y))^2 \sum_{j=1}^{N} \chi_{D_j}(y) + \sum_{j=1}^{N} \frac{c_j^2}{2} (\chi_{D_j}(x) - \chi_{D_j}(y))^2 \sum_{i=1}^{N} \chi_{D_i}(y) \\ &\le \sum_{i=1}^{N} \frac{c_i^2}{2} (\chi_{D_i}(x) - \chi_{D_i}(y))^2 + \sum_{j=1}^{N} \frac{c_j^2}{2} (\chi_{D_j}(x) - \chi_{D_j}(y))^2 \\ &= \sum_{i=1}^{N} c_i^2 (\chi_{D_i}(x) - \chi_{D_i}(y))^2 \end{align*} $$

where we have used the fact that $\sum_{i=1}^N \chi_{D_i} \le 1$. So, we have derived that

$$ \left| \sum_{j=1}^{N} c_j (\chi_{D_j}(x) - \chi_{D_j}(y)) \right|^2 \le 2 \sum_{j=1}^{N} c_j^2 (\chi_{D_j}(x) - \chi_{D_j}(y))^2 $$

from which it follows that

$$ \begin{align*} \|c^{-2}\|_{H^{s'}(\Omega)}^2 &= \int_{\Omega} \int_{\Omega} \frac{|\sum_{j=1}^{N} c_j(\chi_{D_j}(x) - \chi_{D_j}(y))|^2}{|x-y|^{3+2s'}} dx dy \\ &\le 2 \int_{\Omega} \int_{\Omega} \frac{\sum_{j=1}^{N} c_j^2 (\chi_{D_j}(x) - \chi_{D_j}(y))^2}{|x-y|^{3+2s'}} dx dy \\ &\le 2 \sum_{j=1}^{N} c_j^2 \int_{\Omega} \int_{\Omega} \frac{(\chi_{D_j}(x) - \chi_{D_j}(y))^2}{|x-y|^{3+2s'}} dx dy = 2 \sum_{j=1}^{N} c_j^2 \|\chi_{D_j}\|_{H^{s'}(\Omega)}^2 \end{align*} $$
which proves (23). so that finally from (22) and (23) and Assumption 6 we get

$$ \|c^{-2}\|_{H^{s'}(\Omega)}^2 \le 2 \sum_{j=1}^{N} c_j^2 \|\chi_{D_j}\|_{H^{s'}(\Omega)}^2 \le C(s') \sum_{j=1}^{N} c_j^2 |D_j| \left( \frac{|\partial D_j|}{|D_j|} \right)^{2s'} \le \frac{C(L, s')}{r_0^{2s'}} \|c^{-2}\|_{L^2(\Omega)}^2. \quad \square $$

We are now ready to state and prove our main stability result

**PROPOSITION 10.** Assume Assumption 6 and let $c_1^{-1}, c_2^{-1} \in W_N$ and let $\omega^2$ satisfy Assumption 8. Then, there exists a positive constant $K$, depending on $A, r_1, L$, such that,

$$ (24) \qquad \|c_1^{-2} - c_2^{-2}\|_{L^2(\Omega)} \le \frac{1}{\omega^2} e^{K(1+\omega^2 B_2)(|\Omega|/r_0^3)^{\frac{4}{3}}} \|\Lambda_{\omega^2 c_1^{-2}} - \Lambda_{\omega^2 c_2^{-2}}\|_{L(H^{1/2}(\partial\Omega), H^{-1/2}(\partial\Omega))}. $$

*Proof.* To prove our stability estimate we follow the idea of Alessandrini of using CGO solutions but we use slightly different ones than those introduced in [21] and in [1] to obtain better constants in the stability estimates as proposed by [11]. We also use the estimates proposed in [11] (see Theorem 4.4) and due to [13] concerning the case of bounded potentials.

In fact, by Theorem 4.3 of [11], since $c^{-2} \in L^\infty(\Omega)$, $\|c^{-2}\|_{L^\infty(\Omega)} \le B_2$, there exists a positive constant $C = C(\omega^2, B_2, A, r_1)$ such that for every $\zeta \in \mathbb{C}^3$ satisfying $\zeta \cdot \zeta = 0$ and $|\zeta| \ge C$ the equation

$$ -\Delta u - \omega^2 c^{-2} u = 0 $$

has a solution of the form

$$ u(x) = e^{ix \cdot \zeta} (1 + R(x)) $$

where $R \in H^1(\Omega)$ satisfies

$$ \|R\|_{L^2(\Omega)} \le \frac{C}{|\zeta|}, \quad \|\nabla R\|_{L^2(\Omega)} \le C. $$

Let $\xi \in \mathbb{R}^3$ and let $\tilde{\omega}_1$ and $\tilde{\omega}_2$ be unit vectors of $\mathbb{R}^3$ such that $\{\tilde{\omega}_1, \tilde{\omega}_2, \xi\}$ is an orthogonal set of vectors of $\mathbb{R}^3$. Let $s$ be a positive parameter to be chosen later and set for $k = 1, 2$,

$$ (25) \qquad \zeta_k = \begin{cases} (-1)^{k-1} \frac{s}{\sqrt{2}} \left( \sqrt{1 - \frac{|\xi|^2}{2s^2}} \tilde{\omega}_1 + (-1)^{k-1} \frac{1}{\sqrt{2}s} \xi + i \tilde{\omega}_2 \right) & \text{for } \frac{|\xi|}{\sqrt{2}s} < 1, \\ (-1)^{k-1} \frac{s}{\sqrt{2}} \left( (-1)^{k-1} \frac{1}{\sqrt{2}s} \xi + i \left( \sqrt{\frac{|\xi|^2}{2s^2} - 1} \tilde{\omega}_1 + \tilde{\omega}_2 \right) \right) & \text{for } \frac{|\xi|}{\sqrt{2}s} \ge 1. \end{cases} $$

Then a straightforward computation gives

$$ \zeta_k \cdot \zeta_k = 0 $$

for $k = 1, 2$ and

$$ \zeta_1 + \zeta_2 = \xi. $$

Furthermore, for $k = 1, 2$,

$$ (26) \qquad |\zeta_k| = \begin{cases} s & \text{for } \frac{|\xi|}{\sqrt{2s}} < 1, \\ \frac{|\xi|}{\sqrt{2}} & \text{for } \frac{|\xi|}{\sqrt{2s}} \ge 1. \end{cases} $$

Hence,

$$ (27) \qquad |\zeta_k| = \max\{s, \frac{|\xi|}{\sqrt{2}}\}. $$

Then, by Theorem 4.3 of [11], for $|\zeta_1|, |\zeta_2| \ge C_1 = \max\{C_0\omega^2 B_2, 1\}$, with $C_0 = C_0(A, r_1)$, there exist $u_1, u_2$, solutions to $-\Delta u_k - \omega^2 c_k^{-2} u_k = 0$ for $k = 1, 2$ respectively, of the form

$$ (28) \qquad u_1(x) = e^{ix \cdot \zeta_1} (1 + R_1(x)), \quad u_2(x) = e^{ix \cdot \zeta_2} (1 + R_2(x)) $$
with

$$ (29) \qquad \|R_k\|_{L^2(\Omega)} \le \frac{C_0 \sqrt{|\Omega|}}{s} \omega^2 B_2 $$

and

$$ (30) \qquad \|\nabla R_k\|_{L^2(\Omega)} \le C_0 \sqrt{|\Omega|} \omega^2 B_2 $$

for $k = 1, 2$. It is common in the literature to use estimates which contain $\sqrt{|\Omega|}$; Different estimates in terms of $|\Omega|$ are possible and just change the leading constant $C_0$.

Consider again Alessandrini's identity

$$ \int_{\Omega} \omega^2(c_1^{-2} - c_2^{-2})u_1 u_2 dx = \langle (\Lambda_1 - \Lambda_2)u_1|_{\partial\Omega}, u_2|_{\partial\Omega} \rangle, $$

where $u_k \in H^1(\Omega)$ is any solution of $-\Delta u_k - \omega^2 c_k^{-2} u_k = 0$ and $\Lambda_k = \Lambda_{\omega^2 c_k^{-2}}$ for $k = 1, 2$. Inserting the solutions (28) in Alessandrini's identity we derive

$$ (31) \qquad \left| \int_{\Omega} \omega^2 (c_1^{-2} - c_2^{-2}) e^{i\xi \cdot x} dx \right| \le \|\Lambda_1 - \Lambda_2\| \|u_1\|_{H^{1/2}(\partial\Omega)} \|u_2\|_{H^{1/2}(\partial\Omega)} + \left| \int_{\Omega} \omega^2 (c_1^{-2} - c_2^{-2}) e^{i\xi \cdot x} (R_1 + R_2 + R_1 R_2) dx \right| \le \|\Lambda_1 - \Lambda_2\| \|u_1\|_{H^1(\Omega)} \|u_2\|_{H^1(\Omega)} + E(\|R_1\|_{L^2(\Omega)} + \|R_2\|_{L^2(\Omega)} + \|R_1\|_{L^4(\Omega)} \|R_2\|_{L^4(\Omega)}). $$

where $E := \|\omega^2(c_1^{-2} - c_2^{-2})\|_{L^2(\Omega)}$. By (27), (29), and (30) and since $\Omega \subset B_{2R}(0)$ we have

$$ \|u_k\|_{H^1(\Omega)} \le C \sqrt{|\Omega|} (s + |\xi|) e^{Ar_1(s+|\xi|)}, \quad k = 1, 2. $$

Let $s \ge C_2$ so that $s + |\xi| \le e^{Ar_1(s+|\xi|)}$. Then, for $s \ge C_3 = \max(C_1, C_2)$, using (29) and (30) and the standard interpolation inequality $(\|u\|_{L^4(\Omega)} \le \|u\|_{L^6(\Omega)}^{3/4} \|u\|_{L^2(\Omega)}^{1/4})$ we get

$$ (32) \qquad |\omega^2(c_1^{-2} - c_2^{-2})^{\wedge}(\xi)| \le C \sqrt{|\Omega|} \left( e^{4Ar_1(s+|\xi|)} \|\Lambda_1 - \Lambda_2\| + \frac{\omega^2 B_2 E}{s} \right) $$

where the $\omega^2 c_k^{-2}$'s have been extended to all $\mathbb{R}^3$ by zero and $^{\wedge}$ denotes the Fourier transform. Hence, from (32), we get

$$ \int_{|\xi| \le \rho} |\omega^2(c_1^{-2} - c_2^{-2})^{\wedge}(\xi)|^2 d\xi \le C|\Omega|\rho^3 \left( e^{8Ar_1(s+\rho)} \|\Lambda_1 - \Lambda_2\|^2 + \frac{\omega^4 B_2^2 E^2}{s^2} \right) $$

and hence

$$ (33) \qquad \begin{aligned} \|\omega^2(c_1^{-2} - c_2^{-2})^{\wedge}\|_{L^2(\mathbb{R}^3)}^2 \le C|\Omega|\rho^3 & \left( e^{8Ar_1(s+\rho)} \|\Lambda_1 - \Lambda_2\|^2 + \frac{\omega^4 B_2^2 E^2}{s^2} \right) \\ & + \int_{|\xi| \ge \rho} |\omega^2(c_1^{-2} - c_2^{-2})^{\wedge}(\xi)|^2 d\xi \end{aligned} $$

where $C = C(A, r_1)$. By (21) and (23) we have that

$$ \|\omega^2(c_1^{-2} - c_2^{-2})\|_{H^{s'}(\Omega)}^2 \le \frac{C}{r_0^{2s'}} E^2, $$
where $C$ depends on $L, s'$ and hence

$$ \begin{aligned} \rho^{2s'} \int_{|\xi| \ge \rho} |\omega^2(c_1^{-2} - c_2^{-2})^{\wedge}(\xi)|^2 d\xi &\le \int_{|\xi| \ge \rho} |\xi|^{2s'} |\omega^2(c_1^{-2} - c_2^{-2})^{\wedge}(\xi)|^2 d\xi \\ &\le \int_{\mathbb{R}^3} (1 + |\xi|^2)^{s'} |\omega^2(c_1^{-2} - c_2^{-2})^{\wedge}(\xi)|^2 d\xi \le \frac{C}{r_0^{2s'}} E^2. \end{aligned} $$

Hence, we get

$$ \int_{|\xi| \ge \rho} |\omega^2(c_1^{-2} - c_2^{-2})^{\wedge}(\xi)|^2 d\xi \le \frac{CE^2}{r_0^{2s'} \rho^{2s'}} $$

for every $s' \in (0, 1/2)$. Inserting last bound in (33) we derive

$$ \|\omega^2(c_1^{-2} - c_2^{-2})^{\wedge}\|_{L^2(\mathbb{R}^3)}^2 \le C \left( \rho^3 |\Omega| e^{8Ar_1(s+\rho)} \|\Lambda_1 - \Lambda_2\|^2 + \rho^3 |\Omega| \frac{\omega^4 B_2^2 E^2}{s^2} + \frac{E^2}{r_0^{2s'} \rho^{2s'}} \right). $$

where $C = C(L, s')$. To make the last two terms in the right-hand side of the inequality of equal size we pick up

$$ \sqrt[3]{|\Omega|} \rho = \left( \frac{|\Omega|}{r_0^3} \right)^{\frac{2s'}{3(3+2s')}} \left( \frac{1}{\alpha} \right)^{\frac{1}{3+2s'}} s^{\frac{2}{3+2s'}} $$

with $\alpha = \max\{1, \omega^4 B_2^2\}$. Then, by Assumption 6 and observing that we might assume without loss of generality that $|\frac{|\Omega|}{r_0^3}| > 1$. In fact, if this is not the case we can choose a smaller value of $r_0$ so that the condition is satisfied.

$$ \|\omega^2(c_1^{-2} - c_2^{-2})\|_{L^2(\Omega)}^2 \le CE^2 \left(\frac{|\Omega|}{r_0^3}\right)^{\frac{2s'}{3+2s'}} \left( e^{C_4\left(\frac{|\Omega|}{r_0^3}\right)^{\frac{2s'}{3(3+2s')}} s} \left(\frac{\|\Lambda_1 - \Lambda_2\|}{E}\right)^2 + \left(\frac{\alpha}{s^2}\right)^{\frac{2s'}{3+2s'}} \right) $$

for $s \ge C_3$ and where $C$ depends on $s', L, A, r_1$ and $C_4$ depends on $L, A, r_1$. We now make the substitution

$$ s = \frac{1}{C_4 \left( \frac{|\Omega|}{r_0^3} \right)^{\frac{2s'}{3(3+2s')}}} \left| \log \frac{\| \Lambda_1 - \Lambda_2 \|}{E} \right| $$

where we assume that

$$ \frac{\|\Lambda_1 - \Lambda_2\|}{E} < c := e^{-\bar{C} \max\{1, \omega^2 B_2\} \left(\frac{|\Omega|}{r_0^3}\right)^{\frac{2s'}{3(3+2s')}}} $$

with $\bar{C} = \bar{C}(R)$ in order that the constraint $s \ge C_3$ is satisfied. Under this assumption,

$$ (34) \quad \|\omega^2(c_1^{-2} - c_2^{-2})\|_{L^2(\Omega)} \le C(\sqrt{\alpha})^{\frac{2s'}{3+2s'}} \left(\frac{|\Omega|}{r_0^3}\right)^{\frac{2s'}{3+2s'}} \frac{9+10s'}{6(3+2s')} E \left( \left| \log \frac{\|\Lambda_1 - \Lambda_2\|}{E} \right| \right)^{-\frac{2s'}{3+2s'}} $$

where $C = C(L, s', A, r_1)$ and we can rewrite last inequality in the form

$$ (35) \quad E \le C(1 + \omega^2 B_2)^{\frac{2s'}{3+2s'}} \left( \frac{|\Omega|}{r_0^3} \right)^{\frac{2s'}{3+2s'}} \frac{9+10s'}{6(3+2s')} E \left( \left| \log \frac{\|\Lambda_1 - \Lambda_2\|}{E} \right| \right)^{-\frac{2s'}{3+2s'}} $$

which gives

$$ (36) \quad E \le e^{C(1+\omega^2 B_2)\left(\frac{|\Omega|}{r_0^3}\right)^{\frac{9+10s'}{6(3+2s')}}} \|\Lambda_1 - \Lambda_2\| $$
where $C = C(L, s', A, r_1)$. On the other hand if

$$ \frac{\|\Lambda_1 - \Lambda_2\|}{E} \geq c, $$

then

$$ (37) \qquad \|\omega^2(c_1^{-2} - c_2^{-2})\|_{L^2(\Omega)} \le c^{-1}\|\Lambda_1 - \Lambda_2\| \le e^{\bar{C}(1+\omega^2 B_2)\left(\frac{|\Omega|}{r_0^3}\right)^{\frac{1}{3(3+2s')}}}\|\Lambda_1 - \Lambda_2\| $$

Hence, from (36) and (37) and recalling that $s' \in (0, \frac{1}{2})$, we have that

$$ (38) \qquad E \le e^{C(1+\omega^2 B_2)\left(\frac{|\Omega|}{r_0^3}\right)^{\frac{9+10s'}{6(3+2s')}}} \|\Lambda_1 - \Lambda_2\| $$

Choosing $s' = \frac{1}{4}$, we derive

$$ \|c_1^{-2} - c_2^{-2}\|_{L^2(\Omega)} \le \frac{1}{\omega^2} e^{K(1+\omega^2 B_2)(|\Omega|/r_0^3)^{\frac{4}{7}}} \|\Lambda_1 - \Lambda_2\| $$

where $K = K(L, A, r_1, s')$ and the claim follows. $\square$

REMARK 11. Here we state an $L^\infty$-stability estimate, in contrast to the $L^2$-stability estimate in Proposition 10.

Observing that

$$ \frac{1}{\sqrt{|\Omega|}} \|c_1^{-2} - c_2^{-2}\|_{L^2(\Omega)} \le \|c_1^{-2} - c_2^{-2}\|_{L^\infty(\Omega)} \le \frac{C}{r_0^{3/2}} \|c_1^{-2} - c_2^{-2}\|_{L^2(\Omega)}, $$

where $C = C(L)$, and we immediately get the following stability estimate in the $L^\infty$ norm

$$ \|c_1^{-2} - c_2^{-2}\|_{L^\infty(\Omega)} \le \frac{C}{\omega^2} e^{K(1+\omega^2 B_2)(|\Omega|/r_0^3)^{\frac{4}{7}}} \|\Lambda_1 - \Lambda_2\| $$

with $C = C(L)$.

REMARK 12. In [6] the following lower bound of the stability constant has been obtained in the case of a uniform polyhedral partition $D_N$

$$ (39) \qquad C_N \ge \frac{1}{4\omega^2} e^{K_1 N^{\frac{1}{3}}} $$

Choose a uniform cubical partition $D_N$ of $\Omega$ of mesh size $r_0$. Then,

$$ (40) \qquad |\Omega| = Nr_0^3 $$

and estimate (24) of Proposition 10 gives

$$ (41) \qquad C_N = \frac{1}{\omega^2} e^{K(1+\omega^2 B_2)N^{\frac{4}{7}}}, $$

which proves a sharp bound on the Lipschitz constant with respect to $N$ when the global DtN map is known. In [6] a Lipschitz stability estimate has been derived in terms of the local DtN map using singular solutions. This type of solutions allows to recover the unknown piecewise constant wavespeeds by determining it on the outer boundary of the domain and then, by propagating the singularity inside the domain, to recover step by step the wavespeed on the interface of all subdomains of the partition. This iterative procedure does not lead to sharp bounds of the Lipschitz constant appearing in the stability estimate. It would be interesting if one can get a better bound of the Lipschitz constant using oscillating solutions.
REMARK 13. In Lemma 5 we have seen that $F_\omega$ is Fréchet differentiable with Lipschitz derivative $DF_\omega$ for which we have derived an upper bound in terms of the apriori data. From the stability estimates we can easily derive the following lower bound

$$ (42) \qquad \min_{c^{-2} \in W_N; h \in \mathbb{R}^N, \|h\|_{L^\infty(\Omega)} = 1} \|DF_\omega[c^{-2}]h\|_* \ge \omega^2 e^{-K(1+\omega^2 B_2)(\frac{|\Omega|}{r_0^3})^{4/7}}. $$

where $K = K(L, A, r_1)$ and $\| \cdot \|_*$ indicates the norm in $\mathcal{L}(H^{1/2}(\partial\Omega), H^{-1/2}(\partial\Omega))$ i.e.

$$ \|T\|_* = \sup\{\langle Tg, f \rangle | : g, f \in H^{1/2}(\partial\Omega), \|g\|_{H^{1/2}(\partial\Omega)} = \|f\|_{H^{1/2}(\partial\Omega)} = 1\} $$

In fact, by the injectivity of $DF_\omega$

$$ \min_{c^{-2} \in W_N; h \in \mathbb{R}^N, \|h\|_{L^\infty(\Omega)} = 1} \|DF_\omega[c^{-2}]h\|_* = m_0/2 > 0 $$

Then, there exists $h_0$ satisfying $\|h_0\|_{L^\infty(\Omega)} = 1$ and $c_0^{-2} \in W_N$ such that

$$ \|DF_\omega[c_0^{-2}]h_0\|_* \le m_0. $$

Hence, by the definition of $\| \cdot \|_*$ it follows that

$$ |\langle DF_\omega[c_0^{-2}](h_0)g, f \rangle| = \left|\int_{\Omega} h_0 \tilde{u}_0 v_0 \right| \le m_0 \|\tilde{u}_0\|_{H^{1/2}(\partial\Omega)} \|v_0\|_{H^{1/2}(\partial\Omega)} $$

where $\tilde{u}_0$ and $v_0$ are solutions to the equation $(-\Delta - \omega^2 c_0^{-2})u = 0$ in $\Omega$ with boundary data $g$ and $f$, respectively. Proceeding like in the proof of the stability result Proposition 10 and Remark 11 we derive that

$$ 1 = \|h_0^{-2}\|_{L^\infty(\Omega)} \le \frac{1}{\omega^2} e^{K(1+\omega^2 B_2)(\frac{|\Omega|}{r_0^3})^{4/7}} m_0 $$

which gives the lower bound (42).

**3. Computational experiments.** In this section, we numerically compute the stability constant for the inverse problem associated with the Dirichlet-to-Neumann map. We illustrate the stability behaviour and compare it with the analytical bounds derived in section 2. The estimates we provide here are obtained from the definition of the stability constant,

$$ (43) \qquad \|c_1^{-2} - c_2^{-2}\|^2 < \mathcal{C} \|F_\omega(c_1^{-2}) - F_\omega(c_2^{-2})\|^2, $$

where $\|c_1^{-2} - c_2^{-2}\|$ denotes the $L^2$-norm of the functions from the finite dimensional Ansatz space. In particular we consider here a geophysical example of reconstruction where normal data are collected on the boundary. In this situation $c_1$ and $c_2$ are assimilated to two different wavespeeds. Hence the boundary value problem (1) corresponds to the propagation of acoustic wave in the media for a boundary source $g$ using the wavespeeds $c_1$ and $c_2$ respectively. In our experiments, Gaussian shaped (spatial) source functions (see Figure 1) are applied. Then the normal data (measurements of the normal derivative of the field) are acquired on the boundary in order to generate the forward operator. The numerical stability estimates are finally obtained by the knowledge of all quantities of equation (43).

In the Remark 12, we have formulated the stability constant depending on the number of cubical partitions $N$ in the model representation equation (40). This situation is well adapted for numerical applications where the domain is commonly discretized. Hence we want to verify the (exponential) dependence of the stability constant with $N$.

The model (assimilated to a wavespeed here) is defined on a cubical (structured) domain partition of a rectangular block. With increasing $N$, the size of the cubes decreases, possibly non uniformly. We use piecewise constant functions on the cubes to define the wavespeeds following
<figure><img src="image_4.png" /><figcaption>FIG. 1. Illustration of the source shape for a localized boundary source.</figcaption></figure>

the main assumption for the Lipschitz stability to hold. Such a partition can be related to Haar wavelets, where N determines the scale. These naturally introduce approximate representations, that is, when the scale of the approximation is coarser than the finest scale contained in the model.

In order to solve the forward problem, the numerical discretization of the operator is realized using discontinuous Galerkin method, where Dirichlet boundary conditions are invoked. The Dirichlet sources at the top boundary introduce Identity block in the discretized Helmholtz operator and give the following linear problem

$$ (44) \qquad \begin{pmatrix} A_{ii} & A_{i\partial} \\ A_{\partial i} & A_{\partial\partial} \end{pmatrix} \begin{pmatrix} u_i \\ u_\partial \end{pmatrix} = \begin{pmatrix} A_{ii} & A_{i\partial} \\ 0 & I_d \end{pmatrix} \begin{pmatrix} u_i \\ u_\partial \end{pmatrix} = \begin{pmatrix} 0 \\ g \end{pmatrix}, $$

where A represents the discretized operator, i labels interior points and $\partial$ labels boundary points, g has values at the source location and is zero elsewhere. This system verifies $u_\partial = u_{\partial\Omega} = g$ (i.e. Dirichlet boundary condition) and $A_{ii}u_i + A_{i\partial}u_{\partial\Omega} = 0$. The normal derivative data are generated by taking the normal derivative of the solution wavefield u on the surface.

Our experiments use a three dimensional model of size 2.55×1.45×1.22km. The wavespeed $c_1$ is viewed as a reference model (which is known in this test case) and is represented Figure 2 (courtesy Statoil). We also illustrate the different partitions of a model and the notion of approximation. Obviously the larger the number of subdomains is, the more precise will the representation be.

For the computation of the stability estimates we consider $c_2$ as the model shown in Figure 3. This setup can be associated with the ‘true’ subsurface Figure 2 and starting model Figure 3. In this context we have chosen the initial guess with no knowledge of any structures by simply considering a one dimensional variation in depth.

**3.1. Estimates using the full Dirichlet-to-Neumann map.** We consider the full data case where the Gaussian sources (see Figure 1) are positioned on each surface following a regular map. For each source, the data are acquired all over the boundary. We introduce a total of 630 sources and 76538 data points for each.

At a selected partition (number of domains) and frequency, we simulate the data for the two media $c_1$ and $c_2$ and compute the difference, from which we deduce the stability constant following equation (43). The main difference with the standard seismic setup is that we consider data on all the boundary and not only at the top. This last case will be mentioned in the subsection 3.2.

The numerical estimates for the stability constant $C$ should depend on the number of domains $N$ following the expression of the lower and upper bounds defined in the Remark 12, equations (39) and (41). Thus we fix the frequency and estimate the stability for different partitions. The evolution of the estimates and underlying bounds are presented in Figure 4 at two selected frequencies, 5 and 10Hz. We plot on a log log scale the function $\log(C\omega^2)$ to focus on the power of $N$ in the estimates,
<figure><img src="image_5.png" /><figcaption>FIG. 2. Three dimensional representations and horizontal sections at 800m depth of the reference wavespeed ($c_1$) using different partition, i.e. scales. Every scale has a structured (rectangular) decomposition using piecewise constant. The size of the rectangular boxes defines the scale of the wavespeed.</figcaption></figure>

which is the slope of the lines (4/7 for the upper bound and 1/5 for the lower bound).

Regarding the different coefficients in the analytical bounds, $K$ and $K_1$ remain undecided and are numerically approximated so that the bounds match the estimates at best. For instance the numerical value for $K_1$ is obtained from equation ((39)) by computing the average value based on
<figure><img src="image_6.png" /><figcaption>FIG. 3. Three dimensional wavespeed used for the successive estimation of the stability constant ($c_2$), 3D representation (left) and horizontal sections at 800m depth (right).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 4. The black squares represent the computational estimates of the stability constant (■) depending on the number of domains $N$ at selected frequency. The dashed line (—) represents the analytical lower bound and the dotted line (·) the upper bound, estimated with equation (45).</figcaption></figure>

the numerical stability estimates and $K$ is approximated following the same principle:

$$ (45) \qquad K_1 = \frac{1}{n_{st}} \sum_{i=1}^{n_{st}} \frac{\log(4\omega^2 C_i)}{N_i^{1/5}}, \qquad K = \frac{1}{n_{st}} \sum_{i=1}^{n_{st}} \frac{\log(\omega^2 C_i)}{(1 + \omega^2 B_2) N_i^{4/7}} $$

Here, $n_{st}$ is the number of numerical stability constant estimates and $C_i$ the corresponding estimate for partitioning $N_i$. We actually limit the computation of $K$ to use the first scales as it grows too rapidly. The numerical values obtained are given Table 1. We also note that the term $\omega^2 B_2$ of the upper bound equation (41) is relatively small in the geophysical context as we have here $B_2 = 5.10^{-7}$.

We can see that the stability constant increases with the number of subdomains, as expected. There are clearly two states in the evolution of the estimates at the highest frequency (10Hz, Figure 4(b)). For a low number of partitions $N$ the numerical estimates match particularly well the upper bound while at finer scale it follows accurately the lower bound. This is illustrated in Figure 5 where we decompose the two parts of the estimates between the low and high number of domains.
<figure><img src="image_8.png" /><figcaption>FIG. 5. The black squares represent the computational estimates of the stability constant (??) depending on the number of domains N at 10Hz. The left part shows the coarsest scales which match accurately the upper bound (dotted line, ??). On the right the finer scale estimates are accurately anticipated by the lower bound (dashed line, ??). The constants K and K₁ for the computation of the lower and upper bounds are numerically approximated with values given Table 1, following equation (45).</figcaption></figure>

<table><caption>TABLE 1
Numerical estimation of the constant in the analytical bounds formulation for the numerical estimates of the stability (Figure 4), with $B_2 = (1/1400)^2$.</caption><thead><tr><th></th><th>5Hz</th><th>10Hz</th></tr></thead><tbody><tr><th>K<sub>1</sub></th><td>1</td><td>0.7</td></tr><tr><th>K</th><td>0.15</td><td>0.05</td></tr></tbody></table>

Alternatively for a lower frequency, i.e. 5Hz on Figure 4(a), the upper bound appears to increase too rapidly while the lower bound matches accurately the evolution of the stability constant estimates. Hence the upper bound we have obtained here is particularly appropriate for coarse scale and high frequency: when the variation of model is much coarser compared to the wavelength.

**3.2. Seismic inverse problem using partial data.** In realistic geophysical experiments for the reconstruction of subsurface area (seismic tomography), it is more appropriate not to consider the full data but partial data only located on the upper surface. The data obtain from $c_1$ can be seen as field observation (sensor measurement of a seismic event at the surface). The data using $c_2$ are simulation using an ‘initial guess’. For the reconstruction, we mention the full waveform inversion method, where the recovery follows an iterative minimization of the difference between the measurements and simulations, to successively update the initial guess (see [23, 20]). There is also the difference in the boundary conditions where perfectly matched layers (PMLs) or absorbing boundary conditions are invoked instead of the Dirichlet boundary condition for the lateral and bottom boundaries. However the top boundary is a free surface and remains a Dirichlet boundary condition.

For this test case we reproduce the same experiments but limiting the set of sources and the collected data to be at the top boundary only. We define a set of sources at the surface, separated by 160m along the x-axis and 150m along y-axis to generate a regular map of 16 × 10 points. The receivers (data location) are positioned in the same fashion every 60m along the x-axis and 45m along y-axis and generate a regular map of 43 × 32 points, see Figures 6(a) and 6(b). The partial boundary data computed are illustrated for a single centered boundary shot at 5Hz frequency Fig-
ure 6(c).

<figure><img src="image_9.png" /><figcaption>(a) the crosses represent the boundary sources locations.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) The lattice represents the discretization of the data, i.e. the receivers location.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c) Data recovered from a boundary centered shot, i.e. wavefield measured at the receivers location.</figcaption></figure>

<caption>FIG. 6. Illustration of the seismic acquisition set.</caption>

In Figure 7 we compare the stability constant estimates using partial data with the stability constant estimates obtained when considering the full Dirichlet-to-Neumann map as the data. We incorporate the analytical lower bound that was computed in the previous test case.

The numerical estimates of the stability constants for the full and partial data in a log log-scale differ by a constant. This leads us to our conjecture that the log log of the stability constants (as a function of N) of the full and partial data case in the continuous setting differ by a constant.

## REFERENCES

[1] G. ALESSANDRINI, *Stable determination of conductivity by boundary measurements*, Appl. Anal., 27 (1988), pp. 153–172, doi:10.1080/00036818808839730, http://dx.doi.org/10.1080/00036818808839730.

[2] G. ALESSANDRINI AND S. VESSELLA, *Lipschitz stability for the inverse conductivity problem*, Adv. in Appl. Math., 35 (2005), pp. 207–241, doi:10.1016/j.aam.2004.12.002, http://dx.doi.org/10.1016/j.aam.2004.12.002.

[3] G. BAO AND P. LI, *Inverse medium scattering problems for electromagnetic waves*, SIAM J. Appl. Math., 65 (2005), pp. 2049–2066 (electronic), doi:10.1137/040607435, http://dx.doi.org/10.1137/040607435.

[4] G. BAO AND F. TRIKI, *Error estimates for the recursive linearization of inverse medium problems*, J. Comput. Math., 28 (2010), pp. 725–744, doi:10.4208/jcm.1003-m0004, http://dx.doi.org/10.4208/jcm.1003-m0004.

[5] H. BEN-HADJ-ALI, S. OPERTO, AND J. VIRIEUX, *Velocity model-building by 3d frequency-domain, full-waveform inversion of wide-aperture seismic data*, Geophysics, 73 (2008), doi:10.1190/1.2957948.

[6] E. BERETTA, M. V. DE HOOP, AND L. QIU, *Lipschitz stability of an inverse boundary value problem for a Schrödinger-type equation*, SIAM J. Math. Anal., 45 (2013), pp. 679–699, doi:10.1137/120869201, http://dx.doi.org/10.1137/120869201.

[7] E. BERETTA AND E. FRANCINI, *Lipschitz stability for the electrical impedance tomography problem: the complex case*, Comm. Partial Differential Equations, 36 (2011), pp. 1723–1749, doi:10.1080/03605302.2011.552930, http://dx.doi.org/10.1080/03605302.2011.552930.
<figure><img src="image_12.png" /><figcaption>FIG. 7. Comparison of the computational stability estimates using partial data only located on the top boundary (×) and using the full boundary data (■). The dashed line (—) represents the analytical lower bound as found in Figure 4.</figcaption></figure>

[8] K. D. BLAZEK, C. STOLK, AND W. W. SYMES, *A mathematical framework for inverse wave problems in heterogeneous media*, Inverse Problems, 29 (2013), p. 065001.

[9] K. DATCHEV AND M. V. DE HOOP, *Iterative reconstruction of the wavespeed for the wave equation with bounded frequency boundary data*, arXiv preprint arXiv:1506.09014, (2015).

[10] E. B. DAVIES, *Spectral theory and differential operators*, vol. 42 of Cambridge Studies in Advanced Mathematics, Cambridge University Press, Cambridge, 1995, doi:10.1017/CBO9780511623721, http://dx.doi.org/10.1017/CBO9780511623721.

[11] J. FELDMAN, M. SALO, AND G. UHLMANN, *The Calderón problem — An Introduction to Inverse Problems*, unpublished ed., 2015, http://www.math.ubc.ca/~feldman/ibook/.

[12] D. GILBARG AND N. S. TRUDINGER, *Elliptic partial differential equations of second order*, vol. 224 of Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences], Springer-Verlag, Berlin, second ed., 1983.

[13] P. HÄHNER, *A periodic Faddeev-type solution operator*, J. Differential Equations, 128 (1996), pp. 300–308, doi:10.1006/jdeq.1996.0096, http://dx.doi.org/10.1006/jdeq.1996.0096.

[14] P. HÄHNER AND T. HOHAGE, *New stability estimates for the inverse acoustic inhomogeneous medium problem and applications*, SIAM Journal on Mathematical Analysis, 33 (2001), pp. 670–685, doi:10.1137/S0036141001383564, http://dx.doi.org/10.1137/S0036141001383564, arXiv: http://dx.doi.org/10.1137/S0036141001383564.

[15] T. HOHAGE, *Logarithmic convergence rates of the iteratively regularized gauss- newton method for an inverse potential and an inverse scattering problem*, Inverse Problems, 13 (1997), p. 1279, http://stacks.iop.org/0266-5611/13/i=5/a=012.

[16] R. MAGNANINI AND G. PAPI, *An inverse problem for the Helmholtz equation*, Inverse Problems, 1 (1985), pp. 357–370, doi:10.1088/0266-5611/1/4/007, http://dx.doi.org/10.1088/0266-5611/1/4/007.

[17] N. MANDACHE, *Exponential instability in an inverse problem for the Schrödinger equation*, Inverse Problems, 17 (2001), pp. 1435–1444, doi:10.1088/0266-5611/17/5/313, http://dx.doi.org/10.1088/0266-5611/17/5/313.

[18] S. NAGAYASU, G. UHLMANN, AND J.-N. WANG, *Increasing stability in an inverse problem for the acoustic equation*, Inverse Problems, 29 (2013), p. 025012, doi:10.1088/0266-5611/29/2/025012, http://dx.doi.org/10.1088/0266-5611/29/2/025012.

[19] R. G. NOVIKOV, *New global stability estimates for the Gel'fand-Calderon inverse problem*, Inverse Problems, 27 (2011), p. 015001, doi:10.1088/0266-5611/27/1/015001, http://dx.doi.org/10.1088/0266-5611/27/1/015001.

[20] R. G. PRATT, C. SHIN, AND G. J. HICKS, *Gauss-newton and full newton methods in frequency-space seismic waveform inversion*, Geophysical Journal International, 133 (1998), pp. 341–362.

[21] J. SYLVESTER AND G. UHLMANN, *A global uniqueness theorem for an inverse boundary value problem*, Ann. of Math. (2), 125 (1987), pp. 153–169, doi:10.2307/1971291, http://dx.doi.org/10.2307/1971291.

[22] W. W. SYMES, *The seismic reflection inverse problem*, Inverse Problems, 25 (2009), pp. 123008, 39, doi:10.1088/0266-5611/25/12/123008, http://dx.doi.org/10.1088/0266-5611/25/12/123008.

[23] A. TARANTOLA, *Inversion of seismic reflection data in the acoustic approximation*, Geophysics, 49 (1984),
pp. 1259–1266.