plus1sp

[1]Barbara Kaltenbacher

# Simultaneous reconstruction of sound speed and nonlinearity parameter in a
paraxial model of vibro-acoustography in frequency domain

Teresa Rauscher  University of Klagenfurt, Department of Mathematics,
Klagenfurt, Austria, e-mail<EMAIL_ADDRESS>
###### Abstract

In this paper we consider the inverse problem of vibro-acoustography, a
technique for enhancing ultrasound imaging by making use of nonlinear effects.
It amounts to determining two spatially variable coefficients in a system of
PDEs describing propagation of two directed sound beams and the wave resulting
from their nonlinear interaction. To justify the use of Newton's method for
solving this inverse problem, on one hand we verify well-definedeness and
differentiability of the forward operator corresponding to two versions of the
PDE model; on the other hand we consider an all-at-once formulation of the
inverse problem and prove convergence of Newton's method for its solution.

##  1 Introduction

Recently, several approaches for enhancing ultrasound by means of nonlinear
effects have been proposed. In this paper we consider vibro-acoustography
which has originally been proposed in [4, 5] to achieve the enhanced
resolution by high frequency waves while avoiding the drawbacks of scattering
from small inclusions and of stronger attenuation at higher frequencies. The
experiment for image acquisition basically consists of three parts: (i) Two
ultrasound beams of high and slightly different frequencies $\omega_{1}$ and
$\omega_{2}$ are excited at two parts $\Sigma_{1}$, $\Sigma_{2}$ of an array
of piezoelectric transducers (emitters). (ii) They interact nonlinearly in a
focus region, thus exciting a wave that basically propagates at the difference
frequency $\omega_{1}-\omega_{2}$ (which allows to avoid strong scattering and
attenuation). (iii) The latter is eventually measured by a receiver array
$\Gamma$. We refer to, e.g., [10, 18] for a graphical illustration.
The fact that the value of the nonlinearity parameter $\frac{B}{A}$ governing
the interaction depends on the tissue properties and thus varies in space
$\frac{B}{A}=\frac{B}{A}(x)$ yields a means of imaging. Likewise, also the
sound speed typically exhibits spatial variation $c=c(x)$. A modeling and
simulation framework for this imaging technology has been devised in [16, 15].

The aim of this paper is to study the inverse problem of identifying both
$\frac{B}{A}(x)$ and $c(x)$ in an appropriate PDE model, in which these two
quantities appear as coefficients. Some preliminary results on this parameter
identification problem have been obtained in [10], however, without modeling
the excitation waves as sound beams. The latter (along with an analysis of
this model with variable coefficients) is one of the main novel contributions
of this paper. Another one is a convergence analysis of a Newton type method
for solving the inverse problem.

The model we will use here can be derived from a general wave equation for an
acoustic velocity potential $\phi$ given by

$\partial_{t}^{2}\phi-c^{2}\Delta\phi=\partial_{t}|\nabla\phi|^{2}-\frac{1}{2}(\nabla\phi)\cdot\nabla|\nabla\phi|^{2}+\frac{B}{A}\partial_{t}\phi\Delta\phi-\frac{B}{2A}|\nabla\phi|^{2}\Delta\phi,$
(1)

by means of the asymptotic ansatz

$\tilde{\phi}(t,x,\varepsilon)=\varepsilon(\tilde{\phi}_{1}(t,x)+\tilde{\phi}_{2}(t,x))+\varepsilon^{2}\tilde{\psi}(t,x),$
(2)

with small $\varepsilon$ and a paraxial approximation (to take into account
the fact that propagation is concentrated to a certain given axis) as well as
transformation to frequency domain. In (1), $c$ is the speed of sound and
$\frac{B}{A}$ the nonlinearity parameter.

### Paraxial approximation

We assume that the direction of propagation is the $x_{1}$-axis. Then, for the
original variables $(t,x)=(t,x_{1},x^{\prime})$, where
$x^{\prime}=(x_{2},\ldots,x_{d})$, (and marking functions of these variables
by a check) we perform the paraxial change of variables

$\displaystyle(\tau,z,y)=\left(t-\check{\mathfrak{t}},\tilde{\varepsilon}x_{1},\sqrt{\tilde{\varepsilon}}x^{\prime}\right),\quad\textup{
with
}\check{\mathfrak{t}}=\check{\mathfrak{t}}(x_{1},x^{\prime}),\quad\mathfrak{t}=\mathfrak{t}(z,y)$
(3) $\displaystyle\textup{ such that
}\frac{1}{c^{2}}=|\nabla_{x}\check{\mathfrak{t}}|^{2}=|\partial_{x_{1}}\check{\mathfrak{t}}|^{2}+|\nabla_{x^{\prime}}\check{\mathfrak{t}}|^{2}=\tilde{\varepsilon}^{2}|\partial_{z}\mathfrak{t}|^{2}+\tilde{\varepsilon}|\nabla_{y}\mathfrak{t}|^{2},\quad\check{\mathfrak{t}}(0,x^{\prime})=0,\
x^{\prime}\in\tilde{\Omega}_{y}$

where $\tilde{\varepsilon}\ll 1$, so
$\tilde{\varepsilon}\ll\sqrt{\tilde{\varepsilon}}$, and in general
$\tilde{\varepsilon}$ can be different from $\varepsilon$ that was introduced
for the second order approximation in (2) but for the model to make sense the
relation $\varepsilon\leq\tilde{\varepsilon}\leq\sqrt{\varepsilon}$ should
hold cf. [18].111Note that $\check{\mathfrak{t}}(x)$ is the distance of $x$
from the boundary $\\{0\\}\times\tilde{\Omega}_{y}$ in the Riemannian metric
determined by $\frac{1}{c^{2}}$, cf., e.g. [3, Section 7.2.3] For the
derivatives, this change of variables yields

$\displaystyle\partial_{t}=\partial_{\tau},\quad\nabla_{x}=\left(\tilde{\varepsilon}\bigl{(}-\partial_{z}\mathfrak{t}\,\partial_{\tau}+\partial_{z}\bigr{)},\sqrt{\tilde{\varepsilon}}\bigl{(}-\nabla_{y}\mathfrak{t}\,\partial_{\tau}+\nabla_{y})\right),$
(4)
$\displaystyle\Delta_{x}=-\tilde{\varepsilon}^{2}\partial_{z}^{2}\mathfrak{t}\,\partial_{\tau}+\tilde{\varepsilon}^{2}|\partial_{z}\mathfrak{t}|^{2}\partial_{\tau}^{2}-2\tilde{\varepsilon}^{2}\partial_{z}\mathfrak{t}\partial_{z}\partial_{\tau}+\tilde{\varepsilon}^{2}\partial_{z}^{2}-\tilde{\varepsilon}\Delta_{y}\mathfrak{t}\,\partial_{\tau}+\tilde{\varepsilon}|\nabla_{y}\mathfrak{t}|^{2}\partial_{\tau}^{2}-2\tilde{\varepsilon}\nabla_{y}\mathfrak{t}\cdot\nabla_{y}\partial_{\tau}+\tilde{\varepsilon}\Delta_{y}$
$\displaystyle\hskip
14.22636pt=\frac{1}{c^{2}}\partial_{\tau}^{2}+\tilde{\varepsilon}\Bigl{(}-(\tilde{\varepsilon}\partial_{z}^{2}\mathfrak{t}+\Delta_{y}\mathfrak{t})\,\partial_{\tau}-2(\tilde{\varepsilon}\partial_{z}\mathfrak{t}\,\partial_{z}+\nabla_{y}\mathfrak{t}\cdot\nabla_{y})\,\partial_{\tau}+\tilde{\varepsilon}\partial_{z}^{2}+\Delta_{y}\Bigr{)}$

Here and in (3) we have used the fact that $\partial_{t}\check{s}=0$,
$\partial_{t}\check{\mathfrak{t}}=0$, hence $\partial_{\tau}s=0$,
$\partial_{\tau}\mathfrak{t}=0$,
$\partial_{x_{1}}\check{s}=\tilde{\varepsilon}\partial_{z}s$,
$\partial_{x_{1}}\check{\mathfrak{t}}=\tilde{\varepsilon}\partial_{z}\mathfrak{t}$,
$\nabla_{x^{\prime}}\check{\mathfrak{t}}=\sqrt{\tilde{\varepsilon}}\nabla_{y}\mathfrak{t}$

### Frequency domain formulation

We make the time harmonic ansatz

$\displaystyle\phi_{k}(\tau,z,y)$
$\displaystyle=\hat{\phi}_{k}(z,y)e^{\imath\omega_{k}\tau}\qquad\text{ for
}k\in\left\\{1,2\right\\},$ $\displaystyle\psi(\tau,z,y)$
$\displaystyle=\Re{\left\\{\hat{\psi}(z,y)e^{\imath(\omega_{1}-\omega_{2})\tau}\right\\}}.$

Here the excitation frequencies are assumed to satisfy
$\omega_{1}>\omega_{2}>0$, with a comparably small difference frequency
$\omega_{d}:=\omega_{1}-\omega_{2}>0$, and we abbreviate
$\omega_{p}^{3}:=\omega_{1}\omega_{2}\omega_{d}$. (Typical values in
ultrasonics are $\omega_{k}\sim 2\,$MHz, $\omega_{d}\sim 50\,$KHz.)

Altogether, for
$\tilde{\psi}(t,x_{1},x^{\prime})=\psi(\tau,z,y)=\Re\left(\hat{\psi}(z,y)e^{\imath\omega_{d}\tau}\right)$
and
$\tilde{\psi}(t,x_{1},x^{\prime})=\Re\left({\check{\psi}}(x_{1},x^{\prime})e^{\imath\omega_{d}t}\right)$
with
$e^{\imath\omega_{d}t}=e^{\imath\omega_{d}x_{1}/c}e^{\imath\omega_{d}\tau}$ we
end up with the transformation in particular of the wave operator

$\displaystyle\Re\left(\left[-\tfrac{1}{c^{2}}\omega_{d}^{2}-\Delta_{x}\right]{\check{\psi}}(x_{1},x^{\prime})\,e^{\imath\omega_{d}\mathfrak{t}}\,e^{\imath\omega_{d}\tau}\right)=[\tfrac{1}{c^{2}}\partial_{t}^{2}-\Delta_{x}]\tilde{\psi}(t,x_{1},x^{\prime})$
(5)
$\displaystyle=\tilde{\varepsilon}\Re\left(\Bigl{[}\imath\omega_{d}(\tilde{\varepsilon}\partial_{z}^{2}\mathfrak{t}+\Delta_{y}\mathfrak{t})+2\imath\omega_{d}(\tilde{\varepsilon}\partial_{z}\mathfrak{t}\,\partial_{z}+\nabla_{y}\mathfrak{t}\cdot\nabla_{y})-\tilde{\varepsilon}\partial_{z}^{2}-\Delta_{y}\Bigr{]}\hat{\psi}(z,y)e^{\imath\omega_{d}\tau}\right)$

and likewise for $\tilde{\phi}_{k}$, $\hat{\phi}_{k}$.
Due to the boundary condition $\check{\mathfrak{t}}(0,x^{\prime})=0,\
x^{\prime}\in\tilde{\Omega}_{y}$, we also have
$\nabla_{y}\mathfrak{t}(0,x^{\prime})$ and
$\Delta_{y}\mathfrak{t}(0,x^{\prime})$, which allows us to neglect the terms
containing $\nabla_{y}\mathfrak{t}$ and $\Delta_{y}\mathfrak{t}$ to arrive at

$\displaystyle\left[-\tfrac{1}{c^{2}}\omega_{d}^{2}-\Delta_{x}\right]{\check{\psi}}(x_{1},x^{\prime})\,e^{\imath\omega_{d}\mathfrak{t}}\approx\tilde{\varepsilon}\Bigl{[}\imath\omega_{d}\partial_{z}s+2\imath\omega_{d}s\,\partial_{z}-\tilde{\varepsilon}\partial_{z}^{2}-\Delta_{y}\Bigr{]}\hat{\psi}(z,y)$
(6)

The outcome of the nonlinear term is more complicated and depends on the
interplay between the asymptotic ansatz (2) governed by $\varepsilon$ and the
paraxial approximation (3) governed by $\tilde{\varepsilon}$. Depending on the
relation between $\varepsilon$ and $\tilde{\varepsilon}$, this leads to
several different models; see [18] for more details.

### The considered models

The inverse problem of combined nonlinearity imaging and spound speed
reconstruction will be considered for one of these models, namely

$\displaystyle{\imath\omega_{k}\,\partial_{z}s\,\hat{\phi}_{k}+}2\imath\omega_{k}s\partial_{z}\hat{\phi}_{k}-\Delta_{y}\hat{\phi}_{k}=0\mbox{
in }\tilde{\Omega},\qquad k\in\\{1,2\\},$ (7)
$\displaystyle{\imath\omega_{d}\,\partial_{z}s\,\hat{\psi}+}2\imath\omega_{d}s\partial_{z}\hat{\psi}-\tilde{\varepsilon}\partial_{z}^{2}\hat{\psi}-\Delta_{y}\hat{\psi}=\imath\omega_{p}^{3}\tilde{\varepsilon}^{-1}\,\eta\,\hat{\phi}_{1}\,\overline{\hat{\phi}}_{2}\mbox{
in }\tilde{\Omega}$

with boundary conditions

$\displaystyle\hat{\phi}_{k}(0,y)=h_{k}(y),\quad y\in\tilde{\Omega}_{y}\,,$
(8)
$\displaystyle\partial_{\nu_{y}}\hat{\phi}_{k}(z,y)=-\imath\omega_{k}\sigma_{k}\hat{\phi}_{k}(z,y),\quad(z,y)\in(0,L)\times\partial\tilde{\Omega}_{y},$
$\displaystyle\quad k\in\\{1,2\\},$
$\displaystyle-\partial_{z}\hat{\psi}(0,y)=-\imath\omega_{d}\sigma_{0}\hat{\psi}(0,y),\quad\partial_{z}\hat{\psi}(L,y)=-\imath\omega_{d}\sigma_{L}\hat{\psi}(L,y),\quad
y\in\tilde{\Omega}_{y},$
$\displaystyle\partial_{\nu_{y}}\hat{\psi}(z,y)=-\imath\omega_{d}\sigma\hat{\psi}(z,y),\quad(z,y)\in(0,L)\times\partial\tilde{\Omega}_{y},$

and $\tilde{\Omega}=(0,L)\times\tilde{\Omega}_{y}$.

The coefficients we are interested in are related to the sound speed and the
nonlinearity parameter by

$s=\frac{1}{c},\qquad\eta=\frac{B/A+2}{c^{4}},$ (9)

In (8), $h_{k}$, $k=1,2$ models excitation on part of the boundary and the
remaining boundary conditions are supposed to prevent spurious reflections on
the boundary of the computational domain $\tilde{\Omega}$.

######  Remark 1.

To model excitation by an array of piezoelectric transducers, continuity of
the normal velocity over the transducer-fluid interface would induce Neumann
boundary conditions on the velocity potential. In our setting these would read
as $-\partial_{z}\hat{\phi}_{k}(0,\cdot)=h_{k}\mbox{ in }\tilde{\Omega}_{y}$.
In an analysis analogous to the one we provided here, this could be tackled by
differentiating $\hat{\phi}_{k}$ with respect to $z$ and considering the
corresponding initial value problems for
$\hat{\phi}_{k}^{\prime}:=\partial_{z}\hat{\phi}_{k}$. Dirichlet conditions as
used here can be justified a surface force – pressure continuity
$F=[\sigma]\cdot\nu=p\nu$ together with the relation
$p_{k}=\rho\partial_{t}\phi_{k}$ derived from momentum balance. We refer to,
e.g., [6, eqs (12),(13)] and the references therein for the proper choice of
interface conditions in structure-acoustics coupling.

Well-posedness of the forward problem (7), (8) will be studied in Section 2.

From the point of view of the outgoing wave described by $\hat{\psi}$, the
parameter $\tilde{\varepsilon}$ is not necessarily small, since propagation of
sound is non-directed for the $\psi$ field. Therefore, alternatively to (7),
(8), we will consider

$\displaystyle{\imath\omega_{k}\partial_{z}s\,\hat{\phi}_{k}+}2\imath\omega_{k}s\partial_{z}\hat{\phi}_{k}-\Delta_{y}\hat{\phi}_{k}=0\mbox{
in }\tilde{\Omega},\qquad k\in\\{1,2\\},$ (10)
$\displaystyle-\omega_{d}^{2}s^{2}\check{\psi}-\imath\omega_{d}b\Delta_{x}\check{\psi}-\Delta_{x}\check{\psi}=\imath\omega_{p}^{3}\,1_{\Omega}\Bigl{(}\check{\eta}\,P^{-1}\bigl{(}\hat{\phi}_{1}\,\overline{\hat{\phi}}_{2}\bigr{)}\Bigr{)}\mbox{
in }\Omega$

for $\check{\psi}=P^{-1}\hat{\psi}$,
$\check{\eta}={e^{-\imath\omega_{d}\mathfrak{t}}}P^{-1}\eta$ (where $P$ is
defined by $(P\check{\psi})(z,y)=\check{\psi}(x_{1},x^{\prime})$, cf. (4))
with boundary conditions

$\displaystyle\hat{\phi}_{k}(0,\cdot)=h_{k}\mbox{ in
}\tilde{\Omega}_{y}\,,\quad\partial_{\nu_{y}}\hat{\phi}_{k}(z,y)=-\imath\omega_{k}\sigma_{k}\hat{\phi}_{k}$
$\displaystyle\mbox{ in }(0,L)\times\partial\tilde{\Omega}_{y}\quad
k\in\\{1,2\\}$ (11)
$\displaystyle\partial_{\nu}\check{\psi}=-\imath\omega_{d}\sigma\check{\psi}{-\beta\check{\psi}}\mbox{
on }\partial\Omega$

see Section 2.2 for its well-posedness analysis.

Due to (4), the coefficients in the boundary conditions (8), (11) are related
by

$\sigma_{0}(y)=\sigma(0,y)-s(0,y),\quad\sigma_{L}(y)=\sigma(L,y)-s(L,y).$

Thus with the typical choice $\sigma=s$, in the propagation direction
$z||x_{1}$, the absorbing boundary condition in the original coordinates
becomes a Neumann condition in paraxial coordinates.

We wish to mention that the paraxial approximation is also made use of in the
derivation of the Khokhlov-Zabolotskaya-Kuznetsov (KZK) equation

$2s\partial^{2}_{\tau
z}\psi-\nabla_{y}^{2}\psi=\frac{\eta}{2}\partial_{\tau}\left|\partial_{\tau}\psi\right|^{2}$
(12)

[21], see also, e.g., [19] for its analysis. Note however, that in our case,
the quadratic nonlinearity is decoupled and appears as a source term for the
$\psi$ equation, whereas (12) is a nonlinear (more precisely, quasilinear)
equation, whose expansion in frequency domain would lead to an inifinite
system of space-dependent PDEs, similarly to [8].

### The inverse problem

Our aim is to reconstruct $s$ and $\eta$ in the boundary value problem (17),
(18) from measurements of the acoustic pressure

$p^{\textup{meas}}=\imath\omega_{d}\mbox{tr}_{\Gamma}\hat{\psi}\quad\mbox{ in
}\Gamma$ (13)

where $\Gamma=P(\check{\Gamma})$ and $\check{\Gamma}$ is a manifold
representing the receiver array immersed in the acoustic domain $\Omega$.

This can be formulated as an operator equation

$F(s,\eta)=y$ (14)

with the forward operator $F=C\circ S$ being a concatenation of the
(nonlinear) parameter-to-state map

$S:(s,\eta)\mapsto(\hat{\phi}_{1},\hat{\phi}_{2},\hat{\psi})\mbox{ solving
\eqref{PDEs_intro}, \eqref{BCs_intro}}$ (15)

with the (linear) observation operator

$C:(\hat{\phi}_{1},\hat{\phi}_{2},\hat{\psi})\mapsto\imath\omega_{d}\mbox{tr}_{\Gamma}\hat{\psi}.$
(16)

We consider $F$ as an operator $F:X_{s}\times X_{\eta}\to Y$ with the
parameter space $X:=X_{s}\times X_{\eta}$ and the data space $Y$ yet to be
specified.

In Section 3, we will alternatively consider an all-at-once formulation of
this problem, that keeps the parameters and the state as simultaneous
unknowns, thus avoiding the use of a parameter-to-state map $S$.

Identification of the nonlinearity coefficient $\eta$ in time domain models of
nonlinear acoustics has been studied, e.g., in [1, 11, 12, 13, 20]; in
particular [13] also considers simultaneous identification of $s$ and $\eta$.
However, the physical background and therefore also the model differs from the
ones we consider here.

A preliminary analysis of the inverse problem of ultrasound vibro-
acoustography with models similar to those considered here, but still without
the paraxial approximation, can be found in [10]. In [18], models using a
paraxial approximation are derived in time and frequency domain, and the
inverse problem of reconstructing $\eta$ is studied.

The plan of this paper is as follows. Section 2 is devoted to the forward
problem of solving the PDE model for given coefficients. We prove well-
definedness and differentiability of the parameter-to-state map (15) for both
options (7) and (10), in order to justify the application of Newton type
methods for the inverse problems. These are discussed in Section 3, where we
establish convergence of a frozen Newton method for reconstructing $\eta$ and
$s$ in (10). We point to the fact that proving convergence of iterative
regularization methods is notoriously difficult in inverse problems with
boundary observations, as typical for tomographic imaging. This is due to the
fact that convergence criteria, such as the so-called tangential cone
condition, can usually not be verified in the situation of restricted
observations. We therefore here work with a range invariance condition
instead, that indeed can be rigorously verified here and allows to conclude
convergence.

An implementation and numerical experiments with the methods analyzed here is
subject of future work. Some numerical results on the simultaneous
reconstruction of $s$ and $\eta$ in the Westervelt equation (thus related to
this work, but using a model in time domain with a single PDE rather than a
system) based on a Newton type iteration as well, can be found in [13].

##  2 Well-posedness of the forward problem

In this section we will prove well-definedness and differentiability of the
parameter-to-state map for the systems (7), (8) and (10), (11), respectively.

In doing so, we put a particular emphasis on monitoring the smoothness
assumptions on the coefficients, which we aim at keeping minimal in view of
the fact that in practice, $s$ and $\eta$ tend to be only piecewise smooth and
also the inverse problem becomes more ill-posed the higher order the norm in
preimage space needs to be chosen.

To justify the use of Newton's method for solving the inverse problem (14), we
will also prove differentiability of the parameter-to-state map (15).

###  2.1 Well-posedness of paraxial wave propagation with variable
coefficients

Consider

$\displaystyle{\imath\omega_{k}\partial_{z}s\,\hat{\phi}_{k}+}2\imath\omega_{k}s\partial_{z}\hat{\phi}_{k}-\Delta_{y}\hat{\phi}_{k}=0\mbox{
in }\tilde{\Omega},\qquad k\in\\{1,2\\},$ (17)
$\displaystyle{\imath\omega_{d}\partial_{z}s\,\hat{\psi}+}2\imath\omega_{d}s\partial_{z}\hat{\psi}-\tilde{\varepsilon}\partial_{z}^{2}\hat{\psi}-\Delta_{y}\hat{\psi}=\imath\omega_{p}^{3}\tilde{\varepsilon}^{-1}\,\eta\,\hat{\phi}_{1}\,\overline{\hat{\phi}}_{2}\mbox{
in }\tilde{\Omega}$

with boundary conditions

$\displaystyle\hat{\phi}_{k}(0,y)=h_{k}(y),\quad y\in\tilde{\Omega}_{y}\,,$
(18)
$\displaystyle\partial_{\nu_{y}}\hat{\phi}_{k}(z,y)=-\imath\omega_{k}\sigma_{k}\hat{\phi}_{k}(z,y),\quad(z,y)\in(0,L)\times\partial\tilde{\Omega}_{y},$
$\displaystyle\quad k\in\\{1,2\\}$
$\displaystyle-\partial_{z}\hat{\psi}(0,y)=-\imath\omega_{d}\sigma_{0}\hat{\psi}(0,y),\quad\partial_{z}\hat{\psi}(L,y)=-\imath\omega_{d}\sigma_{L}\hat{\psi}(L,y),\quad
y\in\tilde{\Omega}_{y},$
$\displaystyle\partial_{\nu_{y}}\hat{\psi}(z,y)=-\imath\omega_{d}\sigma\hat{\psi}(z,y),\quad(z,y)\in(0,L)\times\partial\tilde{\Omega}_{y}$

and $\tilde{\Omega}=(0,L)\times\tilde{\Omega}_{y}$.

The equations take the form of a (perturbed) Schrödinger equation, hence well-
known techniques for that equation can be adopted here (see, e.g., [14] and
the references therein). Since we here require estimates that are explicit in
terms of appropriate norms of $s$ and $\eta$, we will first of all provide
some energy estimates for solutions of the linear variable coefficient
problems

$\displaystyle\imath\omega\,\widetilde{s}\,u+2\imath\omega
s\partial_{z}u-\Delta_{y}u=f\quad\textup{ in }(0,L)\times\tilde{\Omega}_{y},$
(19) $\displaystyle u(0,y)=h(y),\quad y\in\tilde{\Omega}_{y}\,,$
$\displaystyle\partial_{\nu_{y}}u(z,y)=-\imath\omega\sigma
u(z,y),\quad(z,y)\in(0,L)\times\partial\tilde{\Omega}_{y}$

and

$\displaystyle{\imath\omega\,\widetilde{s}\,u^{b}+}2\imath\omega
s\partial_{z}u^{b}-b\partial_{z}^{2}u^{b}-\Delta_{y}u^{b}=f\quad\textup{ in
}(0,L)\times\tilde{\Omega}_{y},$ (20)
$\displaystyle-\partial_{z}u^{b}(0,y)=-\imath\omega\sigma_{0}u^{b}(0,y),\quad\partial_{z}u^{b}(L,y)=-\imath\omega\sigma_{L}u^{b}(L,y),\quad
y\in\tilde{\Omega}_{y},$
$\displaystyle\partial_{\nu_{y}}u^{b}(z,y)=-\imath\omega\sigma
u^{b}(z,y),\quad(z,y)\in(0,L)\times\partial\tilde{\Omega}_{y}.$

Here we assume all coefficients $s$, $\sigma_{0}$, $\sigma_{L}$, $\sigma$, $b$
to be real valued and positive with $s,\,\frac{1}{s}\,\partial_{z}s\,\in
L^{\infty}(0,L;L^{\infty}(\tilde{\Omega}_{y}))$, $\sigma,\,\frac{1}{\sigma}\in
L^{\infty}(0,L;L^{\infty}(\partial\tilde{\Omega}_{y}))$,
$\sigma_{0},\,\frac{1}{\sigma_{0}},\,\sigma_{L},\,\frac{1}{\sigma_{L}}\,\in
L^{\infty}(\tilde{\Omega}_{y})$ and $b$ constant. The following two identities
will be used repeatedly

$2s\Re(\partial_{z}v\,\overline{v})=s\frac{d}{dz}|v|^{2}=\frac{d}{dz}\Bigl{(}s|v|^{2}\Bigr{)}-\partial_{z}s\,|v|^{2},$
$\int_{\tilde{\Omega}_{y}}-\Delta_{y}u\,w\,dy=\int_{\tilde{\Omega}_{y}}\nabla_{y}u\cdot\nabla_{y}w\,dy+\imath\omega\int_{\partial\tilde{\Omega}_{y}}\sigma
u\,w\,d\Gamma(y)\,.$

Moreover, for (20), we will assume the Poincaré-Friedrichs type estimate on
the domain $\tilde{\Omega}_{y}$

$\|v\|_{L^{2}(\tilde{\Omega}_{y})}^{2}\leq
C_{0}\|\sqrt{\sigma}v\|_{L^{2}(\partial\tilde{\Omega}_{y})}^{2}+C_{1}\|\nabla_{y}v\|_{L^{2}(\tilde{\Omega}_{y})}^{2}\qquad\forall
v\in H^{1}(\tilde{\Omega}_{y})$ (21)

to hold for some $C_{0}$, $C_{1}>0$.

######  Proposition 1.

If $s\in W^{1,\infty}(0,L;L^{\infty}(\tilde{\Omega}_{y}))$,
$\partial_{z}\widetilde{s}\in L^{\infty}(\tilde{\Omega})$,
$h,\,\Delta_{y}h+f(0)\in L^{2}(\tilde{\Omega}_{y})$, $f\in
L^{\infty}(0,L;L^{2}(\tilde{\Omega}_{y}))\cap
W^{1,1}(0,L;L^{2}(\tilde{\Omega}_{y}))$ then a solution of (19), exists, is
unique and satisfies the estimate

$\displaystyle\|\nabla_{y}u(z)\|_{L^{\infty}(0,L;L^{2}(\tilde{\Omega}_{y}))}^{2}+\omega\|\sqrt{s}\partial_{z}u\|_{L^{\infty}(0,L;L^{2}(\tilde{\Omega}_{y}))}^{2}+\omega\|\sqrt{s}u\|_{L^{\infty}(0,L;L^{2}(\tilde{\Omega}_{y}))}^{2}$
$\displaystyle\leq
Ce^{L(s_{\sim}+\mu)}\Bigl{(}\omega|\sqrt{s(0)}h|_{L^{2}(\tilde{\Omega}_{y})}^{2}+|\frac{1}{2\sqrt{s(0)}}(\Delta_{y}h+f(0))|_{L^{2}(\tilde{\Omega}_{y})}^{2}$
$\displaystyle\hskip
85.35826pt+\tfrac{1}{\omega}\left\|\tfrac{1}{\sqrt{s}}f\right\|_{L^{\infty}(0,L;L^{2}(\tilde{\Omega}_{y}))}^{2}+\left\|\tfrac{1}{\sqrt{s}}f\right\|_{L^{1}(0,L;L^{2}(\tilde{\Omega}_{y}))}^{2}+\left\|\tfrac{1}{\sqrt{s}}\partial_{z}f\right\|_{L^{1}(0,L;L^{2}(\tilde{\Omega}_{y}))}^{2}\Bigr{)}$

with some constants $\mu>0$, $C>0$ independent of $\omega$.

If $s\in W^{1,\infty}(0,L;L^{\infty}(\tilde{\Omega}_{y}))$, $f\in
L^{2}(0,L;L^{2}(\tilde{\Omega}_{y}))$, (21) holds and the relative smallness
conditions on $s$, $\partial_{z}s$, $\omega$

$C_{0}|{\widetilde{s}-}\partial_{z}s|_{L^{\infty}}<1\textup{ and
}2C_{1}\omega^{2}|s|_{L^{\infty}}^{2}<b(1-C_{0}|{\widetilde{s}-}\partial_{z}s|_{L^{\infty}})$
(22)

are satisfied with $\underline{c}$, $\lambda$ as in (21), then a solution of
(20), exists, is unique and satisfies the estimate

$\|u^{b}\|_{H^{1}((0,L)\times\tilde{\Omega}_{y})}^{2}\leq
C\left(1+\frac{1}{\omega^{2}}\right)\|f\|_{L^{2}(0,L;L^{2}(\Omega)}^{2}$

with some constant $C>0$ independent of $\omega$.

  * Proof.

As in standard evolutionary PDEs, the proof is based on a Galerkin
approximation by eigenfunctions of the negative Laplacian, energy estimates
and taking weak limits, cf. e.g., [3]. We will here focus on the energy
estimates, since the other steps of the proof are relatively straightforward
for the linear problems under consideration.

We will multiply the PDEs with appropriate test functions and integrate over
$\tilde{\Omega}_{y}$, abbreviating by $|\cdot|_{L^{2}(\tilde{\Omega}_{y})}$,
$|\cdot|_{L^{2}(\partial\tilde{\Omega}_{y})}$ the $L^{2}$ norm on
$\tilde{\Omega}_{y}$ and $\partial\tilde{\Omega}_{y}$, respectively.

Testing this way first of all (19) with $\overline{u}$ and taking the real and
imaginary parts we obtain

$|\nabla_{y}u(z)|_{L^{2}(\tilde{\Omega}_{y})}^{2}=2\omega\int_{\tilde{\Omega}_{y}}s\Im(\partial_{z}u(z)\,\overline{u}(z)))\,dy+\int_{\tilde{\Omega}_{y}}\Re(f\,\overline{u})\,dy,$
(23)
$\displaystyle\omega\frac{d}{dz}|\sqrt{s}u|_{L^{2}(\tilde{\Omega}_{y})}^{2}(z)+\omega|\sqrt{\sigma}u(z)|_{L^{2}(\partial\tilde{\Omega}_{y})}^{2}$
$\displaystyle=-\omega\int_{\tilde{\Omega}_{y}}({\widetilde{s}-}\partial_{z}s)\,|u(z)|^{2}\,dy+\int_{\tilde{\Omega}_{y}}\Im(f\,\overline{u})\,dy$
(24)
$\displaystyle\leq\omega(s_{\sim}+\mu)|\sqrt{s}u|_{L^{2}(\tilde{\Omega}_{y})}^{2}(z)+\frac{1}{4\mu\omega}|\frac{1}{\sqrt{s}}f|_{L^{2}(\tilde{\Omega}_{y})}^{2}(z)$

for all $z\in(0,L)$, with
$s_{\sim}:=\|\frac{\widetilde{s}-\partial_{z}s}{s}\|_{L^{\infty}(0,L;L^{\infty}(\tilde{\Omega}_{y}))}$,
where we have used Young's inequality.

Doing the same with (20) and integrating over $(0,L)$ we obtain

$\displaystyle\|\nabla_{y}u^{b}\|_{L^{2}(0,L;L^{2}(\tilde{\Omega}_{y}))}^{2}+b\|\partial_{z}u^{b}\|_{L^{2}(0,L;L^{2}(\tilde{\Omega}_{y}))}^{2}=2\omega\int_{0}^{L}\int_{\tilde{\Omega}_{y}}s\Im(\partial_{z}u^{b}\,\overline{u}^{b})\,dy\,dz+\int_{0}^{L}\int_{\tilde{\Omega}_{y}}\Re(f\,\overline{u}^{b})\,dy\,dz,$
(25)
$\displaystyle\omega|\sqrt{s}u^{b}|_{L^{2}(\tilde{\Omega}_{y})}^{2}(L)+\omega\|\sqrt{\sigma}u^{b}\|_{L^{2}(0,L;L^{2}(\partial\tilde{\Omega}_{y}))}^{2}+b\omega\sigma_{0}|u^{b}(0)|_{L^{2}(\tilde{\Omega}_{y})}^{2}+b\omega\sigma_{L}|u^{b}(L)|_{L^{2}(\tilde{\Omega}_{y})}^{2}$
(26)
$\displaystyle=\omega|\sqrt{s}u^{b}|_{L^{2}(\tilde{\Omega}_{y})}^{2}(0)-\omega\int_{0}^{L}\int_{\tilde{\Omega}_{y}}({\widetilde{s}-}\partial_{z}s)\,|u^{b}|^{2}\,dy\,dz+\int_{0}^{L}\int_{\tilde{\Omega}_{y}}\Im(f\,\overline{u}^{b})\,dy\,dz.$

We also differentiate (19) with respect to $z$ and test with
$\partial_{z}\overline{u}$. Analogously to above, this yields

$|\nabla_{y}\partial_{z}u(z)|_{L^{2}(\tilde{\Omega}_{y})}^{2}=\omega\int_{\tilde{\Omega}_{y}}\Bigl{(}{\partial_{z}\widetilde{s}\Im(\partial_{z}u(z)\,\overline{u}(z)))+}2s\Im(\partial_{z}^{2}u(z)\,\partial_{z}\overline{u}(z)))\Bigr{)}\,dy+\int_{\tilde{\Omega}_{y}}\Re(\partial_{z}f\,\partial_{z}\overline{u})\,dy,$
(27)
$\displaystyle{\omega\int_{\tilde{\Omega}_{y}}\partial_{z}\widetilde{s}\Re(\partial_{z}u(z)\,\overline{u}(z)))\,dy+}\omega\frac{d}{dz}|\sqrt{s}\partial_{z}u|_{L^{2}(\tilde{\Omega}_{y})}^{2}(z)+\omega|\sqrt{\sigma}\partial_{z}u(z)|_{L^{2}(\partial\tilde{\Omega}_{y})}^{2}$
(28)
$\displaystyle\leq\omega(s_{+}+\mu)|\sqrt{s}\partial_{z}u|_{L^{2}(\tilde{\Omega}_{y})}^{2}(z)+\frac{1}{4\mu\omega}|\frac{1}{\sqrt{s}}\partial_{z}f|_{L^{2}(\tilde{\Omega}_{y})}^{2}(z).$

with
$s_{+}=\|\frac{\widetilde{s}+\partial_{z}s}{s}\|_{L^{\infty}(0,L;L^{\infty}(\tilde{\Omega}_{y}))}$.

For the initial value problem (with respect to $z$) for the Schrödinger
equation (19), we combine (24)+$\rho\cdot$(28) yields, for
$\eta(z)=\eta_{0}(z)+\rho\eta_{1}(z)=|\sqrt{s}u|_{L^{2}(\tilde{\Omega}_{y})}^{2}(z)+\rho|\sqrt{s}\partial_{z}u|_{L^{2}(\tilde{\Omega}_{y})}^{2}(z)$
that

$\displaystyle\frac{d}{dz}\eta_{0}$
$\displaystyle\leq(s_{\sim}+\mu)\eta_{0}+\frac{1}{4\omega^{2}\mu}|\frac{1}{\sqrt{s}}f|_{L^{2}(\tilde{\Omega}_{y})}^{2}$
$\displaystyle\frac{d}{dz}\eta_{1}$
$\displaystyle\leq(s_{+}+\mu)\eta_{1}+\frac{1}{4\omega^{2}\mu}|\frac{1}{\sqrt{s}}\partial_{z}f|_{L^{2}(\tilde{\Omega}_{y})}^{2}+\|\partial_{z}\widetilde{s}/s\|_{L^{\infty}(\tilde{\Omega})}\sqrt{\eta_{0}\eta_{1}}$

thus, using Young's inequality with factors $\frac{1}{2\sqrt{\rho}}$,
$\frac{\sqrt{\rho}}{2}$ for the last term

$\displaystyle\frac{d}{dz}\eta\leq\beta\eta+\alpha\textup{ where }$
$\displaystyle\beta=\max\\{s_{\sim},s_{+}\\}+\mu+\sqrt{\rho}/2\,\|\partial_{z}\widetilde{s}/s\|_{L^{\infty}(\tilde{\Omega})},$
$\displaystyle\alpha(z)=\frac{1}{4\omega^{2}\mu}\Bigl{(}|\frac{1}{\sqrt{s}}f(z)|_{L^{2}(\tilde{\Omega}_{y})}^{2}+|\frac{1}{\sqrt{s}}\partial_{z}f(z)|_{L^{2}(\tilde{\Omega}_{y})}^{2}\Bigr{)}.$

An application of Gronwall's inequality yields

$\eta(z)=|\sqrt{s}u|_{L^{2}(\tilde{\Omega}_{y})}^{2}(z)+\rho|\sqrt{s}\partial_{z}u|_{L^{2}(\tilde{\Omega}_{y})}^{2}(z)\leq
e^{\beta
z}\Bigl{(}|\sqrt{s}u|_{L^{2}(\tilde{\Omega}_{y})}^{2}(0)+\rho|\sqrt{s}\partial_{z}u|_{L^{2}(\tilde{\Omega}_{y})}^{2}(0)\Bigr{)}\int_{0}^{z}e^{\beta(z-a)}\alpha(a)\,da,$

where we can insert the initial conditions and their differentiated version
substituting from the PDE

$\sqrt{s(0)}u(0)=\sqrt{s(0)}h,\qquad\sqrt{s(0)}\partial_{z}u(0)=\frac{1}{2\imath\omega\sqrt{s(0)}}\Bigl{(}f(0)+\Delta_{y}h-\imath\omega\,\widetilde{s}(0)\,h\Bigr{)}.$

Using this together with the Cauchy-Schwarz and Young's inequalities in (23)
yields the estimate

$\displaystyle|\nabla_{y}u(z)|_{L^{2}(\tilde{\Omega}_{y})}^{2}\leq
2\omega\int_{\tilde{\Omega}_{y}}|\sqrt{s}\partial_{z}u|(z)\,|\sqrt{s}u|(z)\,dy+\int_{\tilde{\Omega}_{y}}|\tfrac{1}{\sqrt{s}}f|\,|\sqrt{s}u|\,dy$
(29)
$\displaystyle\leq\omega\lambda_{1}|\sqrt{s}\partial_{z}u|_{L^{2}(\tilde{\Omega}_{y})}^{2}(z)+\tfrac{\omega}{\lambda_{1}}|\sqrt{s}u|_{L^{2}(\tilde{\Omega}_{y})}^{2}(z)+\tfrac{\lambda_{2}}{2}|\sqrt{s}u|_{L^{2}(\tilde{\Omega}_{y})}^{2}(z)+\tfrac{1}{2\lambda_{2}}|\tfrac{1}{\sqrt{s}}f|_{L^{2}(\tilde{\Omega}_{y})}^{2}(z)$
$\displaystyle{\leq\omega\lambda_{1}\eta+\tfrac{\rho\lambda_{1}}{4\omega(\lambda_{1}^{2}-\rho)}|\tfrac{1}{\sqrt{s}}f|_{L^{2}(\tilde{\Omega}_{y})}^{2}(z)}$

where we have chosen $\lambda_{2}$ such that
$\omega\lambda_{1}=\rho(\tfrac{\omega}{\lambda_{1}}+\tfrac{\lambda_{2}}{2})$
in order to make optimal use of $\eta$.

To obtain an estimate for the two point boundary value problem (with respect
to $z$) for the perturbed Schrödinger equation (20) we combine (25)
$+\frac{C_{0}}{C_{1}\omega}\cdot$ (26) and use (21), (22) to obtain

$\displaystyle\frac{1}{C_{1}}\|u^{b}\|_{L^{2}(0,L;L^{2}(\tilde{\Omega}_{y}))}^{2}+b\|\partial_{z}u^{b}\|_{L^{2}(0,L;L^{2}(\tilde{\Omega}_{y}))}^{2}$
$\displaystyle\leq\|\nabla_{y}u^{b}\|_{L^{2}(0,L;L^{2}(\tilde{\Omega}_{y}))}^{2}+b\|\partial_{z}u^{b}\|_{L^{2}(0,L;L^{2}(\tilde{\Omega}_{y}))}^{2}+\frac{C_{0}}{C_{1}}|\sqrt{s}u^{b}|_{L^{2}(\tilde{\Omega}_{y})}^{2}(L)$
$\displaystyle+\frac{C_{0}}{C_{1}}\|\sqrt{\sigma}u^{b}\|_{L^{2}(0,L;L^{2}(\partial\tilde{\Omega}_{y}))}^{2}+\frac{C_{0}}{C_{1}}(b\sigma_{0}-|s|_{L^{\infty}})|u^{b}(0)|_{L^{2}(\tilde{\Omega}_{y})}^{2}+\frac{C_{0}}{C_{1}}b\sigma_{L}|u^{b}(L)|_{L^{2}(\tilde{\Omega}_{y})}^{2}$
$\displaystyle\leq
2\omega\int_{0}^{L}\int_{\tilde{\Omega}_{y}}s|\Im(\partial_{z}u^{b}\,\overline{u}^{b})|\,dy\,dz+\int_{0}^{L}\int_{\tilde{\Omega}_{y}}|\Re(f\,\overline{u}^{b})|\,dy\,dz-\frac{C_{0}}{C_{1}}\int_{0}^{L}\int_{\tilde{\Omega}_{y}}({\widetilde{s}-}\partial_{z}s)\,|u^{b}|^{2}\,dy\,dz+\frac{C_{0}}{C_{1}\omega}\int_{0}^{L}\int_{\tilde{\Omega}_{y}}|\Im(f\,\overline{u}^{b})|\,dy\,dz$
$\displaystyle\leq(\mu+\tfrac{C_{0}}{C_{1}}|{\widetilde{s}-}\partial_{z}s|_{L^{\infty}})\|u^{b}\|_{L^{2}(0,L;L^{2}(\tilde{\Omega}_{y}))}^{2}$
$\displaystyle\qquad+\frac{1}{2\mu}\Bigl{(}4\omega^{2}|s|_{L^{\infty}}^{2}\|\partial_{z}u^{b}\|_{L^{2}(0,L;L^{2}(\tilde{\Omega}_{y}))}^{2}+(1+(\tfrac{C_{0}}{C_{1}\omega})^{2})\|f\|_{L^{2}(0,L;L^{2}(\tilde{\Omega}_{y}))}^{2}\Bigr{)},$

where $C_{0}$, $C_{1}$ are as in (21). In order to be able to achieve

$C_{1}\mu+C_{0}|{\widetilde{s}-}\partial_{z}s|_{L^{\infty}}<1\quad\textup{ and
}\quad\frac{2\omega^{2}|s|_{L^{\infty}}^{2}}{\mu}<b,$

we invoke (22) and choose
$\mu\in(\frac{2\omega^{2}|s|_{L^{\infty}}^{2}}{b},\frac{1-C_{0}|{\widetilde{s}-}\partial_{z}s|_{L^{\infty}}}{C_{1}})$.
∎

######  Remark 2.

Since we will apply the part of Proposition 1 concerning (19) with large
$\omega=\omega_{k}$, $k=1,2$, we track frequency dependence of the estimates
explicitly there. On the other hand, the energy estimate in Proposition 1 for
(20) will only be used for the (small) difference frequency
$\omega=\omega_{d}$ and therefore the restrictions on the size of $\omega$
made there do not compromise practical applicability in our context.

A wavenumber explicit estimate for (20) follows by application of the
transformation (5) to the Helmholtz equation without having to impose a
smallness condition (22), see Proposition 2 in the next subsection.

Using Proposition 1, we can conclude the following results on the parameter-
to-state map

$S:(s,\eta)\mapsto(\hat{\phi}_{1},\hat{\phi}_{2},\hat{\psi})\mbox{ solving
\eqref{PDEs}, \eqref{BCs}}$ (30)

and its linearization $dS(s,\eta)(\underline{\delta
s},\underline{\delta\eta}):=(\underline{\delta\phi}_{1},\underline{\delta\phi}_{2},\underline{\delta\psi})$
defined by the system

$\displaystyle\imath\omega_{k}\partial_{z}s\,\underline{\delta\phi}_{k}+2\imath\omega_{k}s\partial_{z}\underline{\delta\phi}_{k}-\Delta_{y}\underline{\delta\phi}_{k}={-\imath\omega_{k}\partial_{z}\underline{\delta
s}\,\hat{\phi}_{k}}-2\imath\omega_{k}\underline{\delta
s}\partial_{z}\hat{\phi}_{k}\mbox{ in }\tilde{\Omega},\qquad k\in\\{1,2\\},$
(31)
$\displaystyle\imath\omega_{d}\partial_{z}s\,\underline{\delta\psi}+2\imath\omega_{d}s\partial_{z}\underline{\delta\psi}-\tilde{\varepsilon}\partial_{z}^{2}\underline{\delta\psi}-\Delta_{y}\underline{\delta\psi}={-\imath\omega_{d}\partial_{z}\underline{\delta
s}\,\hat{\psi}}-2\imath\omega_{d}\underline{\delta s}\partial_{z}\hat{\psi}$
$\displaystyle\hskip
99.58464pt+\imath\omega_{p}^{3}\tilde{\varepsilon}^{-1}\,\Bigl{(}\underline{\delta\eta}\,\hat{\phi}_{1}\,\overline{\hat{\phi}}_{2}+\eta(\underline{\delta\phi}_{1}\,\overline{\hat{\phi}}_{2}+\hat{\phi}_{1}\,\overline{\underline{\delta\phi}}_{2})\Bigr{)}\mbox{
in }\tilde{\Omega}$

with boundary conditions

$\displaystyle\underline{\delta\phi}_{k}(0,y)=0,\quad
y\in\tilde{\Omega}_{y}\,,$ (32)
$\displaystyle\partial_{\nu_{y}}\underline{\delta\phi}_{k}(z,y)=-\imath\omega_{k}\sigma_{k}\underline{\delta\phi}_{k}(z,y),\quad(z,y)\in(0,L)\times\partial\tilde{\Omega}_{y},$
$\displaystyle\quad k\in\\{1,2\\},$
$\displaystyle-\partial_{z}\underline{\delta\psi}(0,y)=-\imath\omega_{d}\sigma_{0}\underline{\delta\psi}(0,y),\quad\partial_{z}\underline{\delta\psi}(L,y)=-\imath\omega_{d}\sigma_{L}\underline{\delta\psi}(L,y),\quad
y\in\tilde{\Omega}_{y},$
$\displaystyle\partial_{\nu_{y}}\underline{\delta\psi}(z,y)=-\imath\omega_{d}\sigma\underline{\delta\psi}(z,y),\quad(z,y)\in(0,L)\times\partial\tilde{\Omega}_{y}.$

######  Corollary 1.

Assume $h,\,\Delta_{y}h\in L^{2}(\tilde{\Omega}_{y})$. Then $S$ is well-
defined by (44) as an operator

$S:\mathcal{D}_{s}^{1}\times L^{2}(0,L;L^{p}(\tilde{\Omega}_{y}))\to V$

and Gâteaux differentiable for all $(s,\eta)\in\mathcal{D}_{s}^{2}\times
L^{2}(0,L;L^{p}(\tilde{\Omega}_{y}))$ as an operator
$W^{1,\infty}(0,L;L^{\infty}(\tilde{\Omega}_{y}))\times
L^{2}(0,L;L^{p}(\tilde{\Omega}_{y}))\to V$, where

$\displaystyle\mathcal{D}_{s}^{1}$ $\displaystyle=\\{s\in
W^{1,\infty}(0,L;L^{\infty}(\tilde{\Omega}_{y}))\,:\,\textup{
\eqref{smallness_perturbedSchr} holds for
}\omega=\omega_{d}=\omega_{1}-\omega_{2}\\},$ (33)
$\displaystyle\mathcal{D}_{s}^{2}$ $\displaystyle=\mathcal{D}_{s}^{1}\cap
W^{2,q}(0,L;L^{\infty}(\tilde{\Omega}_{y})),$ $\displaystyle V$
$\displaystyle=\Bigl{(}L^{\infty}(0,L;H^{1}(\tilde{\Omega}_{y})\cap
W^{1,\infty}(0,L;L^{2}(\tilde{\Omega}_{y})\Bigr{)}^{2}\times
H^{1}((0,L)\times\tilde{\Omega}_{y})$

for some $p,q\in\begin{cases}&(1,\infty]\textup{ if
$\tilde{\Omega}_{y}\subseteq\mathbb{R}^{2}$}\\\ &[1,\infty]\textup{ if
$\tilde{\Omega}_{y}\subseteq\mathbb{R}^{1}$}\end{cases}$, with
$dS(s,\eta)(\underline{\delta
s},\underline{\delta\eta}):=(\underline{\delta\phi}_{1},\underline{\delta\phi}_{2},\underline{\delta\psi})$
defined by (31), (32).

  * Proof.

The first part of Proposition 1 with $\omega=\omega_{k}$ implies that the
solutions $\hat{\phi}_{1},\,\hat{\phi}_{2}\,\in
L^{\infty}(0,L;H^{1}(\tilde{\Omega}_{y}))\cap
W^{1,\infty}(0,L;L^{2}(\tilde{\Omega}_{y}))$.
Setting
$f:=\imath\omega_{p}^{3}\tilde{\varepsilon}^{-1}\,\eta\,\hat{\phi}_{1}\,\overline{\hat{\phi}}_{2}\,\in
L^{2}(0,L;L^{2}(\tilde{\Omega}_{y}))$ and since $H^{1}(\tilde{\Omega}_{y})$
continuously embeds into $L^{2p/(p-1)}$, we conclude $\hat{\psi}\in
H^{1}((0,L)\times\tilde{\Omega}_{y})$ from the second part of Proposition 1
with $\omega=\omega_{d}$.

It is readily checked that the first variation of $S$ at $(s,\eta)$ in
direction $(\underline{\delta s},\underline{\delta\eta})$ is given by the
solution of (31), (32) and that $dS(s,\eta)$ is a linear operator. To prove
boundedness of $dS(s,\eta):W^{1,\infty}(0,L;L^{\infty}(\tilde{\Omega}_{y}))\to
V$ for $(s,\eta)\in\mathcal{D}_{s}^{2}\times
L^{2}(0,L;L^{p}(\tilde{\Omega}_{y}))$ we apply the first part of Proposition 1
with $\omega=\omega_{k}$ and $f:=-2\imath\omega_{k}\underline{\delta
s}\partial_{z}\hat{\phi}_{k}$, $h=0$. To prove that $f\in
L^{\infty}(0,L;L^{2}(\tilde{\Omega}_{y}))\cap
W^{1,1}(0,L;L^{2}(\tilde{\Omega}_{y}))$ we need somewhat more regularity of
$\hat{\phi}_{k}$. The latter can be derived by considering the PDE for
$\hat{\phi}_{k}:=\partial_{z}\hat{\phi}_{k}$, (19) with
$f=-2\imath\omega_{k}\partial_{z}s\partial_{z}\hat{\phi}_{k}\in
L^{\infty}(0,L;L^{2}(\tilde{\Omega}_{y}))\cap
W^{1,1}(0,L;L^{2}(\tilde{\Omega}_{y}))$ for $s\in\mathcal{D}_{s}^{2}$ and
implies $\hat{\phi}_{k}\in W^{1,\infty}(0,L;H^{1}(\tilde{\Omega}_{y}))\cap
W^{1,\infty}(0,L;L^{2}(\tilde{\Omega}_{y}))$.
Finally, application of the second part of Proposition 1 with
$\omega=\omega_{d}$ and $f:=-2\imath\omega_{d}\underline{\delta
s}\partial_{z}\hat{\psi}+\imath\omega_{p}^{3}\tilde{\varepsilon}^{-1}\,\Bigl{(}\underline{\delta\eta}\,\hat{\phi}_{1}\,\overline{\hat{\phi}}_{2}+\eta(\underline{\delta\phi}_{1}\,\overline{\hat{\phi}}_{2}+\hat{\phi}_{1}\,\overline{\underline{\delta\phi}}_{2})\Bigr{)}\in
L^{2}(0,L;L^{2}(\tilde{\Omega}_{y}))$ with $f(0)\in L^{2}(\tilde{\Omega}_{y})$
implies the required bound on $\underline{\delta\psi}$ in
$H^{1}((0,L)\times\tilde{\Omega}_{y})$. ∎

###  2.2 Well-posedness, using the Helmholtz equation for the propagating wave

Alternatively to the paraxial form of the PDE for $\hat{\psi}$ we consider

$\displaystyle{\imath\omega_{k}\partial_{z}s\,\hat{\phi}_{k}+}2\imath\omega_{k}s\partial_{z}\hat{\phi}_{k}-\Delta_{y}\hat{\phi}_{k}=0\mbox{
in }\tilde{\Omega},\qquad k\in\\{1,2\\},$ (34)
$\displaystyle-\omega_{d}^{2}s^{2}\check{\psi}-\Delta_{x}\check{\psi}=\imath\omega_{p}^{3}\,1_{\Omega}\Bigl{(}\check{\eta}\,P^{-1}\bigl{(}\hat{\phi}_{1}\,\overline{\hat{\phi}}_{2}\bigr{)}\Bigr{)}\mbox{
in }\Omega$

for $\check{\psi}=P^{-1}\hat{\psi}$,
$\check{\eta}={e^{-\imath\omega_{d}\mathfrak{t}}}P^{-1}\eta$ (cf (4)) with
boundary conditions

$\displaystyle\hat{\phi}_{k}(0,\cdot)=h_{k}\mbox{ in
}\tilde{\Omega}_{y}\,,\quad\partial_{\nu_{y}}\hat{\phi}_{k}(z,y)=-\imath\omega_{k}\sigma_{k}\hat{\phi}_{k}$
$\displaystyle\mbox{ in }(0,L)\times\partial\tilde{\Omega}_{y},\quad
k\in\\{1,2\\},$ (35)
$\displaystyle\partial_{\nu}\check{\psi}=-\imath\omega_{d}\sigma\check{\psi}\mbox{
on }\partial\Omega$

with $P(\Omega)\supseteq\tilde{\Omega}=(0,L)\times\tilde{\Omega}_{y}$, where
$P$ is the paraxial transform defined by (3) and $1_{\Omega}$ the extension by
zero operator for a function defined on $P^{-1}(\tilde{\Omega})$ to all of
$\Omega$. This allows to use a possibly larger propagation domain for the
$\psi$ wave as compared to the one for the $\phi_{1}$, $\phi_{2}$ waves, and
to get rid of the smallness assumption (22).

Well-posedness and Gâteaux differentiability of the forward operator defined
by the system (34), (35) follows analogously to Corollary 2 by combining the
first part of Proposition 1 with known results on the Helmholtz equation with
impedance boundary conditions, see, e.g., [17, Chapter 8]. For completeness
and to track the required coefficient regularity we here provide the essential
arguments.

We start with some energy estimates for the Helmholtz equation that can be
obtained by applying the testing strategy from [17, Chapter 8] (where the
constant coefficient case with $\sigma=s$ is considered) to our a variable
coefficient setting. Again, these results are in principle available in the
literature, see, e.g., [7], but we aim at making the required regularity of
$s$ visible for our purposes.

Testing the general Helmholtz equation on a smooth domain $\Omega$ with
impedance boundary conditions

$\displaystyle-\omega^{2}s^{2}u-\Delta_{x}u$ $\displaystyle=f\mbox{ in
}\Omega$ (36) $\displaystyle\partial_{\nu}u+\imath\omega\sigma u$
$\displaystyle=h\mbox{ on }\partial\Omega$

with $u$ and taking real and imaginary parts yields

$\displaystyle\|\nabla_{x}u\|_{L^{2}(\Omega)}^{2}-\omega^{2}\|su\|_{L^{2}(\Omega)}^{2}=\int_{\Omega}\Re(f\,u)\,dx,$
(37)
$\displaystyle\omega\int_{\partial\Omega}\sigma\,|u|^{2}\,d\Gamma=\int_{\Omega}\Im(f\,u)\,dx.$

On a starlike domain $\Omega$ in two space dimensions $d=2$ with

$x\cdot\nu(x)\geq c_{\Omega}>0$ (38)

on $\partial\Omega$ we can also use $v(x):=x\cdot\nabla u(x)$ as a test
function (provided it is contained in $H^{1}(\Omega)$) and consider the real
part to arrive at

$\displaystyle\omega^{2}\int_{\Omega}(1+x\cdot\tfrac{\nabla
s}{s})|s\,u|^{2}\,dx-\frac{1}{2}\omega^{2}\int_{\partial\Omega}|s\,u|^{2}\,x\cdot\nu(x)\,d\Gamma(x)+\frac{1}{2}\int_{\partial\Omega}|\nabla
u|^{2}\,x\cdot\nu(x)\,d\Gamma(x)$ (39)
$\displaystyle+\omega\Re\left(\imath\sigma\int_{\partial\Omega}u\
x\cdot\overline{\nabla
u}\,d\Gamma(x)\right)=\Re\left(\int_{\Omega}f\,x\cdot\overline{\nabla
u}\,dx\right).$

Here we have used the identities

$\displaystyle 2\Re\left(s^{2}u\,x\cdot\nabla
u\right)=x\cdot\Bigl{(}\nabla|su|^{2}-2s\nabla s|u|^{2}\Bigr{)},$
$\displaystyle\nabla(x\cdot\overline{\nabla u})=\overline{\nabla
u}+(x\cdot\nabla)\overline{\nabla u},\quad 2\Re\left(\nabla
u\cdot(x\cdot\nabla)\overline{\nabla u}\right)=\nabla\cdot(|\nabla
u|^{2}x)-d|\nabla u|^{2}$

as well as the divergence Theorem and $d=2$.

This yields the following results in a low and higher regularity regime of $s$
and correspondingly estimates on $u$ with different frequency dependence of
the constants.

######  Proposition 2.

If $\Omega$ is a bounded $C^{1,1}$ smooth or convex domain, $s\in
L^{\infty}(\Omega)$, $f\in H^{1}(\Omega)^{*}$, $\sigma\not=0$ a.e on
$\partial\Omega$, then a solution of (36) exists, is unique and satisfies the
estimate

$\|\nabla_{x}u\|_{L^{2}(\Omega)}^{2}+\|u\|_{L^{2}(\Omega)}^{2}\leq
C(\omega)\bigl{(}\|f\|_{H^{1}(\Omega)^{*}}^{2}+\|h\|_{H^{-1/2}(\partial\Omega)}^{2}\bigr{)}$

with some constant $C(\omega)>0$ depending on $\omega$.

If additionally $f\in L^{2}(\Omega)$, $h\in H^{1/2}(\partial\Omega)$ then

$\|D_{x}^{2}u\|_{L^{2}(\Omega)}^{2}+\|\nabla_{x}u\|_{L^{2}(\Omega)}^{2}+\|u\|_{L^{2}(\Omega)}^{2}\leq
C(\omega)\bigl{(}\|f\|_{L^{2}(\Omega)}^{2}+\|h\|_{H^{1/2}(\partial\Omega)}^{2}\bigr{)}.$
(40)

Here $D_{x}^{2}u$ denotes the (weak derivative) Hessian of $u$.

If additionally $\Omega$ is a star-shaped domain with (38), $s\geq 0$,
$\ln(s)\in W^{1,\infty}(\Omega)$, $f\in L^{2}(\Omega)$, $h=0$, $d=2$ and the
relative smallness/largeness conditions

$c_{s}=\sup_{x\in\Omega}(-x\cdot\tfrac{\nabla s(x)}{s(x)})<1,\quad\sigma\geq
c_{\sigma}>0$ (41)

are satisfied, then the solution of (36) satisfies the estimate

$\frac{1}{\omega^{2}}\|D_{x}^{2}u\|_{L^{2}(\Omega)}^{2}+\|\nabla_{x}u\|_{L^{2}(\Omega)}^{2}+\omega^{2}\|u\|_{L^{2}(\Omega)}^{2}\leq
C\left(1+\frac{1}{\omega^{2}}\right)\|f\|_{L^{2}(\Omega)}^{2}$ (42)

with some $C$ independent of $\omega$.

  * Proof.

The well-posedness proof in the low regularity regime follows analogously to
the one in the constant coefficient case [17, Chapter 8]. Since it can hardly
be found in the literature for our setting (variably slowness, impedance
boundary conditions) we provide it here. We rewrite (36) as

$(-\Delta+\text{id})u=F+Gu$

with the bounded and boundedly invertible operator
$(-\Delta+\text{id}):H^{1}(\Omega)\to H^{1}(\Omega)^{*}$, the bounded operator
$G:H^{1}(\Omega)\to H^{1}(\Omega)^{*}$ and the element $F\in
H^{1}(\Omega)^{*}$ defined by

$\displaystyle\langle(-\Delta+\text{id})u,v\rangle_{(H^{1})^{*},H^{1}}:=\int_{\Omega}(\nabla
u\cdot\nabla\overline{v}+u\,\overline{v})\,dx,$ $\displaystyle\langle
Gu,v\rangle_{(H^{1})^{*},H^{1}}:=\int_{\Omega}(1+\omega^{2}s^{2})\,u\,\overline{v}\,dx-\imath\omega\int_{\partial\Omega}\sigma
u\overline{v}\,d\Gamma,$ $\displaystyle\langle
F,v\rangle_{(H^{1})^{*},H^{1}}:=\int_{\Omega}f\,\overline{v}\,dx+\int_{\partial\Omega}h\overline{v}\,d\Gamma$

for any $v\in H^{1}(\Omega)$. This is equivalent to

$(I-K)u=b\quad\text{ with }K=(-\Delta+\text{id})^{-1}G,\
b=(-\Delta+\text{id})^{-1}F,$

where $b\in H^{1}(\Omega)$ and the operator $K$ is bounded from
$L^{2}(\Omega)$ to $H^{1}(\Omega)$ and thus compact when considered as an
operator from $L^{2}(\Omega)$ into itself. We can therefore apply Fredholm's
alternative for showing bijectivity of $I-K$, then the Bounded Inverse Theorem
on the Banach space $L^{2}(\Omega)$ and finally lift regularity of $u$ to
$H^{1}(\Omega)$ by means of the fixed point identity $u=Ku+b$. In order to use
Fredholm's alternative, we have to prove injectivity of $I-K$. To this end,
assume that $(I-K)w=0$ for some $w\in L^{2}(\Omega)$, which due to $w=Kw$ is
automatically contained in $H^{1}(\Omega)$. The second line of (37) implies
that $w$ satisfies homogeneous Dirichlet, hence by the impedance condition
also homogeneous Neumann boundary conditions and can therefore be extended by
zero to an $H^{2}(\mathbb{R}^{3})$ function222$w\in H^{2}(\Omega)$ follows
from elliptic regularity and $-\Delta_{x}w=\omega^{2}s^{2}w\in L^{2}(\Omega)$
provided $\Omega$ is $C^{1,1}$ or convex; This $H^{2}$ requirement is only
needed because of the used uniqueness results from [2, Section 8.3];
$\hat{w}$ that then satisfies the homogeneous Helmholtz equation
$-\Delta_{x}\hat{w}+\hat{s}\hat{w}$ on all of $\mathbb{R}^{3}$ with
$\hat{s}:=1_{\Omega}(s-1)+1$ and Sommerfeld radiation conditions. From the
results in [2, Section 8.3] we conclude that $\hat{w}\equiv 0$ and thus
$w\equiv 0$.

The higher order regularity result (40) follows from elliptic reguarity and
the fact that $u$ satisfies $-\Delta_{x}u=\tilde{f}$,
$\partial_{\nu}u=\tilde{h}$ with $\tilde{f}=f+\omega^{2}s^{2}u\in
L^{2}(\Omega)$,
$\tilde{h}=h-\imath\omega\sigma\textup{tr}_{\partial\Omega}u\in
H^{1/2}(\partial\Omega).$

From (39) with (38), (41), applying Cauchy-Schwarz' inequality we conclude

$\displaystyle\omega^{2}(1-c_{\Omega})\|su\|_{L^{2}(\Omega)}^{2}+\frac{1}{2}\gamma\|\nabla
u\|_{L^{2}(\partial\Omega)}^{2}\leq\frac{1}{2}\omega^{2}\int_{\partial\Omega}|s\,u|^{2}\,x\cdot\nu(x)\,d\Gamma(x)$
(43)
$\displaystyle\qquad+\omega^{2}\Bigl{(}\sqrt{\int_{\partial\Omega}\sigma^{2}|x|^{2}|u|^{2}\,d\Gamma(x)}\,+\sqrt{\int_{\partial\Omega}|x|^{2}|f|^{2}\,d\Gamma(x)}\Bigr{)}\|\nabla
u\|_{L^{2}(\partial\Omega)}.$

Hence by Young's inequality and the second line in (37), that is,
$\omega\|\sqrt{\sigma}u\|_{L^{2}(\partial\Omega)}^{2}\leq\|f/s\|_{L^{2}(\Omega)}\|su\|_{L^{2}(\Omega)}$,

$\displaystyle\omega^{2}(1-c_{\Omega})\|su\|_{L^{2}(\Omega)}^{2}$
$\displaystyle\leq\frac{\omega^{2}}{2}\int_{\partial\Omega}|s\,u|^{2}\,x\cdot\nu(x)\,d\Gamma(x)+\frac{\omega^{2}}{\gamma}\int_{\partial\Omega}\sigma^{2}|x|^{2}|u|^{2}\,d\Gamma(x)\,+\frac{1}{\gamma}\int_{\partial\Omega}|x|^{2}|f|^{2}\,d\Gamma(x)$
$\displaystyle\leq\omega\Bigl{(}\frac{\textup{d}(\Omega)}{2}\,\left\|\frac{s^{2}}{\sigma}\right\|_{L^{\infty}(\Omega)}+\frac{\textup{d}(\Omega)^{2}}{\gamma}\|\sigma\|_{L^{\infty}(\Omega)}\Bigr{)}\|f/s\|_{L^{2}(\Omega)}\|su\|_{L^{2}(\Omega)}\
+\frac{\textup{d}(\Omega)^{2}}{\gamma}\|f\|_{L^{2}(\Omega)}^{2}$

with $d(\Omega)=\sup_{x\in\Omega}|x|$, hence by another application of Young's
inequality (with factors $\frac{1-c_{\Omega}}{2}$ and
$\frac{1}{2(1-c_{\Omega})}$)

$\omega^{2}\|su\|_{L^{2}(\Omega)}^{2}\leq C\|f\|_{L^{2}(\Omega)}^{2}$

with
$C=\frac{1}{(1-c_{\Omega})^{2}}\Bigl{(}\frac{\textup{d}(\Omega)}{2}\,\left\|\frac{s^{2}}{\sigma}\right\|_{L^{\infty}(\Omega)}+\frac{\textup{d}(\Omega)^{2}}{\gamma}\|\sigma\|_{L^{\infty}(\Omega)}\Bigr{)}^{2}\|\frac{1}{s}\|_{L^{\infty}(\Omega)}^{2}+\frac{2}{(1-c_{\Omega})}\frac{\textup{d}(\Omega)^{2}}{\gamma}.$
Together with the first identity in (37), this yields

$\|\nabla_{x}u\|_{L^{2}(\Omega)}^{2}+\omega^{2}\|su\|_{L^{2}(\Omega)}^{2}\leq\Bigl{(}2C+\frac{1}{\omega^{2}}\left\|\frac{1}{s}\right\|_{L^{\infty}(\Omega)}^{2}\Bigr{)}\|f\|_{L^{2}(\Omega)}^{2}$

and again, elliptic regularity yields (42). ∎

Analogously to Corollary 2 we obtain well-definedness and differentiability of
the parameter-to-state map.

######  Corollary 2.

Under the assumptions of Corollary 2, but without (22) and with $V$ replaced
by $V=\Bigl{(}L^{\infty}(0,L;H^{1}(\tilde{\Omega}_{y})\cap
W^{1,\infty}(0,L;L^{2}(\tilde{\Omega}_{y})\Bigr{)}^{2}\times H^{2}(\Omega)$,
the parameter-to-state map

$S:(s,\eta)\mapsto(\hat{\phi}_{1},\hat{\phi}_{2},\hat{\psi})\mbox{ solving
\eqref{PDEsTransf}, \eqref{BCsTransf}}$ (44)

is well-defined as an operator
$S:W^{1,\infty}(0,L;L^{\infty}(\tilde{\Omega}_{y}))\times L^{2}(\Omega)\to V$
and Gâteaux differentiable for all
$(s,\eta)\in(W^{1,\infty}(0,L;L^{\infty}(\tilde{\Omega}_{y}))\cap
W^{2,q}(0,L;L^{\infty}(\tilde{\Omega}_{y}))\times L^{2}(\Omega)$ as an
operator $W^{1,\infty}(0,L;L^{\infty}(\tilde{\Omega}_{y}))\times
L^{2}(\Omega)\to V$.

##  3 Convergence of a frozen Newton method for the inverse problem

Based on the parameter-to-state map $S$ defined in (15) and analyzed in
Section 2 as well as the observation operator (16), we can write the inverse
problem as

$F(s,\eta)=y$ (45)

with $y=p^{\textup{meas}}$ and the forward operator $F=C\circ S$ and apply
Newton's method to (45), using the fact that $F^{\prime}(s,\eta)=C\circ
S^{\prime}(s,\eta)$ and relying on our verification of Gâteaux
differentiability from Section 2.

Alternatively we here consider an all-at-once formulation

$\mathbb{F}(s,\eta,\hat{\phi}_{1},\hat{\phi}_{2},\psi)=y$ (46)

with $y=(0,0,0,0,0,p^{\textup{meas}})$ based on the variational form of the
system of initial-boundary value / boundary value problems (17), (18)

$\displaystyle\langle\mathbb{F}_{k}(\vec{q}),\vec{w}\rangle_{\mathbb{X}}:=$
$\displaystyle\int_{0}^{L}\Bigl{\\{}\,\int_{\tilde{\Omega}_{y}}\Bigl{(}{\imath\omega_{k}\partial_{z}s_{k}\,\hat{\phi}_{k}\overline{w}_{k}+}2\imath\omega_{k}s_{k}\partial_{z}\hat{\phi}_{k}\,\overline{w}_{k}+\nabla_{y}\hat{\phi}_{k}\cdot\nabla_{y}\overline{w}_{k}\,\Bigr{)}dy\,+\imath\omega_{k}\int_{\partial\tilde{\Omega}_{y}}\sigma_{k}\hat{\phi}_{k}\,\overline{w}_{k}\,d\Gamma(y)\,\Bigr{\\}}\,dz,$
(47)
$\displaystyle\langle\mathbb{F}_{3}(\vec{q}),\vec{w}\rangle_{\mathbb{X}}:=$
$\displaystyle\int_{0}^{L}\int_{\tilde{\Omega}_{y}}\Bigl{\\{}{\imath\omega_{d}\partial_{z}s_{1}\,\hat{\psi}\overline{v}+}2\imath\omega_{d}s_{1}\partial_{z}\hat{\psi}\,\overline{v}+\nabla_{y}\hat{\psi}\cdot\nabla_{y}\overline{v}+\tilde{\varepsilon}\,\partial_{z}\hat{\psi}\,\partial_{z}\overline{v}-(\imath\omega_{p}^{3}\tilde{\varepsilon}^{-1}\,\eta\,\hat{\phi}_{1}\,\overline{\hat{\phi}}_{2})\,\overline{v}\Bigr{\\}}\,dy\,dz$
$\displaystyle+\imath\omega_{d}\int_{0}^{L}\int_{\partial\tilde{\Omega}_{y}}\sigma\hat{\psi}\,\overline{v}\,d\Gamma(y)\,dz+\imath\omega_{d}\,\tilde{\varepsilon}\,\Bigl{(}\sigma_{0}\int_{\tilde{\Omega}_{y}}\hat{\psi}(0)\,\overline{v}(0)\,dy+\sigma_{L}\int_{\tilde{\Omega}_{y}}\hat{\psi}(L)\,\overline{v}(L)\,dy\Bigr{)},$
$\displaystyle\langle\mathbb{F}_{k+3}(\vec{q}),\vec{w}\rangle_{\mathbb{X}}:=$
$\displaystyle\int_{\tilde{\Omega}_{y}}(\hat{\phi}_{k}(0)-h_{k})\,\overline{u}_{k}\,dy,$
$\displaystyle\langle\mathbb{F}_{6}(\vec{q}),\vec{w}\rangle_{\mathbb{X}}:=$
$\displaystyle\int_{\Gamma}\imath\omega_{d}\hat{\psi}\,\overline{u}\,d\Gamma$

for $k\in\\{1,2\\}$, where we have used the abbreviation
$\vec{q}=(s_{1},s_{2},\eta,\hat{\phi}_{1},\hat{\phi}_{2},\hat{\psi})$ and
where $\vec{w}=(w_{1},w_{2},v,u_{1},u_{2},u)$ varies over some space of test
functions. Note that we have introduced a second copy of the variable $s$ in
order to be able to verify the range invariance condition from [9] that allows
to prove convergence of a frozen Newton method. The additional constraint

$0=\mathcal{P}(\vec{q}):=s_{1}-s_{2},$

that will be imposed via a penalty term during the Newton type iteration (53)
takes care of merging these two $s$ versions into one as the iteration
proceeds.

The formal linearization of
$\mathbb{F}=(\mathbb{F}_{1},\ldots,\mathbb{F}_{6})$ at some $\vec{q}^{0}$ is
given by

$\displaystyle\langle\mathbb{F}_{k}^{\prime}(\vec{q})\underline{\delta\vec{q}},\vec{w}\rangle_{\mathbb{X}}=$
$\displaystyle\int_{0}^{L}\Bigl{\\{}\int_{\tilde{\Omega}_{y}}\Big{(}\imath\omega_{k}({\partial_{z}\underline{\delta
s}_{k}\hat{\phi}_{k}+\partial_{z}s_{k}\underline{\delta\phi}_{k}}+2\underline{\delta
s}_{k}\partial_{z}\hat{\phi}_{k}+2s_{k}\partial_{z}\underline{\delta\phi}_{k})\,\overline{w}_{k}+\nabla_{y}\underline{\delta\phi}_{k}\cdot\nabla_{y}\overline{w}_{k}\,\Big{)}dy\,$
(48) $\displaystyle\hskip
28.45274pt+\imath\omega_{k}\int_{\partial\tilde{\Omega}_{y}}\sigma_{k}\underline{\delta\phi}_{k}\,\overline{w}_{k}\,d\Gamma(y)\Bigr{\\}}dz,$
$\displaystyle\langle\mathbb{F}_{3}^{\prime}(\vec{q})\underline{\delta\vec{q}},\vec{w}\rangle_{\mathbb{X}}=$
$\displaystyle\int_{0}^{L}\int_{\tilde{\Omega}_{y}}\Bigl{\\{}\imath\omega_{d}({\partial_{z}\underline{\delta
s}_{1}\hat{\psi}+\partial_{z}s_{1}\underline{\delta\psi}}+2\underline{\delta
s}_{1}\partial_{z}\hat{\psi}+2s_{1}\partial_{z}\underline{\delta\psi})\,\overline{v}+\nabla_{y}\underline{\delta\psi}\cdot\nabla_{y}\overline{v}+\tilde{\varepsilon}\,\partial_{z}\underline{\delta\psi}\,\partial_{z}\overline{v}$
$\displaystyle\qquad-\imath\omega_{p}^{3}\tilde{\varepsilon}^{-1}\,(\underline{\delta\eta}\,\hat{\phi}_{1}\,\overline{\hat{\phi}}_{2}+\eta\,\underline{\delta\phi}_{1}\,\overline{\hat{\phi}}_{2}+\eta\,\hat{\phi}_{1}\,\overline{\underline{\delta\phi}}_{2})\,\overline{v}\Bigr{\\}}\,dy\,dz$
$\displaystyle+\imath\omega_{d}\int_{0}^{L}\int_{\partial\tilde{\Omega}_{y}}\sigma\underline{\delta\psi}\,\overline{v}\,d\Gamma(y)\,dz+\imath\omega_{d}\,\tilde{\varepsilon}\,\Bigl{(}\sigma_{0}\int_{\tilde{\Omega}_{y}}\underline{\delta\psi}(0)\,\overline{v}(0)\,dy+\sigma_{L}\int_{\tilde{\Omega}_{y}}\underline{\delta\psi}(L)\,\overline{v}(L)\,dy\Bigr{)},$
$\displaystyle\langle\mathbb{F}_{k+3}^{\prime}(\vec{q})\underline{\delta\vec{q}},\vec{w}\rangle_{\mathbb{X}}=$
$\displaystyle\int_{\tilde{\Omega}_{y}}\underline{\delta\phi}_{k}(0)\,\overline{u}_{k}\,dy,$
$\displaystyle\langle\mathbb{F}_{6}^{\prime}(\vec{q})\underline{\delta\vec{q}},\vec{w}\rangle_{\mathbb{X}}=$
$\displaystyle\imath\omega_{d}\int_{\Gamma}\underline{\delta\psi}\,\overline{u}\,d\Gamma.$

It is a bounded linear operator and thus it is the Gâteaux derivative of
$\mathbb{F}$ when considered, e.g., as a mapping

$\displaystyle\mathbb{F}:\mathbb{X}\to\mathbb{Y}\quad\textup{ with }$ (49)
$\displaystyle\mathbb{X}:={W^{1,p}(0,L;L^{p}(\tilde{\Omega}_{y}))\times(L^{p}(0,L;L^{p}(\tilde{\Omega}_{y})))^{2}}\times(W^{1,\infty}(0,L;L^{\infty}(\tilde{\Omega}_{y}))\cap
L^{2}(0,L;H^{1}(\tilde{\Omega}_{y})))^{3},$
$\displaystyle\mathbb{Y}:=(L^{2}(0,L;H^{1}(\tilde{\Omega}_{y})^{*}))^{2}\times(H^{1}((0,L)\times\tilde{\Omega}))^{*}\times(L^{2}(\tilde{\Omega}))^{2}\times
L^{2}(\Gamma)$

for some $p\in[1,\infty]$. Note that this slightly differs from the function
space setting suggested for the reduced setting by Corollary 2, but these
differences are essential for the verification of convergence conditions for
the frozen Newton method, see (51) below. Thus, we are here making use of the
additional freedom in choosing the function spaces that we gain by using an
all-at-once formulation (46) as compared to a reduced one (45).

We can achieve the range invariance relation

$\mathbb{F}(\vec{q})-\mathbb{F}(\vec{q}^{0})=\mathbb{F}^{\prime}(\vec{q}^{0})r(\vec{q})$
(50)

by setting $r(\vec{q}):=\underline{\delta\vec{q}}=(\underline{\delta
s}_{1},\underline{\delta
s}_{2},\underline{\delta\eta},\underline{\delta\phi}_{1},\underline{\delta\phi}_{2},\underline{\delta\psi})$
with

$\displaystyle\underline{\delta s}_{k}=$
$\displaystyle{\frac{\hat{\phi}_{k}}{\hat{\phi}_{k}^{0}}\cdot(s_{k}-s_{k}^{0})+\frac{1}{(\hat{\phi}_{k}^{0})^{2}}\int_{0}^{z}(s_{k}-s_{k}^{0})(z^{\prime})\Bigl{(}\hat{\phi}_{k}^{0}\,\partial_{z}\hat{\phi}_{k}-\hat{\phi}_{k}\,\partial_{z}\hat{\phi}_{k}^{0}\Bigr{)}(z^{\prime})\,dz^{\prime}}$
$\displaystyle\underline{\delta\eta}=$
$\displaystyle\frac{\hat{\phi}_{1}\,\overline{\hat{\phi}_{2}}}{\hat{\phi}_{1}^{0}\,\overline{\hat{\phi}_{2}^{0}}}\cdot(\eta-\eta^{0})+\frac{(\hat{\phi}_{1}-\hat{\phi}_{1}^{0})(\overline{\hat{\phi}_{2}}-\overline{\hat{\phi}_{2}^{0}})}{\hat{\phi}_{1}^{0}\,\overline{\hat{\phi}_{2}^{0}}}\cdot\eta^{0}$
$\displaystyle{-\frac{\omega_{d}\,\tilde{\varepsilon}}{\omega_{p}^{3}\,\hat{\phi}_{1}^{0}\,\overline{\hat{\phi}_{2}^{0}}}\left(\hat{\psi}\,\partial_{z}(s_{1}-s_{1}^{0})+2(s_{1}-s_{1}^{0})\,\partial_{z}\hat{\psi}-\hat{\psi}^{0}\,\partial_{z}\underline{\delta
s}_{1}-2\underline{\delta s}_{1}\,\partial_{z}\hat{\psi}^{0}\right),}$
$\displaystyle\underline{\delta\phi}_{k}=$
$\displaystyle\hat{\phi}_{k}-\hat{\phi}_{k}^{0},\qquad\qquad\underline{\delta\psi}=\hat{\psi}-\hat{\psi}^{0}.$

To this end, we have to choose $\hat{\phi}_{1}^{0}$, $\hat{\phi}_{2}^{0}$,
$\hat{\psi}^{0}$ such that the denominators in the above expression are
bounded away from zero. This is possible, due to the fact that in the all-at-
once formulation they need not be PDE solutions corresponding to the
coefficients $\underline{\delta s}_{1}^{0}$, $\underline{\delta s}_{2}^{0}$,
$\underline{\delta\eta}^{0}$. Relying on this, we can also establish closeness
of $r$ to the identity in the sense that

$r(\vec{q})-(\vec{q}-\vec{q}^{0})=\bigl{(}\underline{\delta
s}_{1}-(s_{1}-s_{1}^{0}),\underline{\delta
s}_{2}-(s_{2}-s_{2}^{0}),\underline{\delta\eta}-(\eta-\eta^{0}),0,0,0\bigr{)}$

where

$\displaystyle\|\underline{\delta
s}_{k}-(s_{k}-s_{k}^{0})\|_{L^{p}(L^{p})}\leq\|\tfrac{1}{\hat{\phi}_{k}^{0}}\|_{L^{\infty}(L^{\infty})}\|\hat{\phi}_{k}-\hat{\phi}_{k}^{0}\|_{L^{\infty}(L^{\infty})}\|s_{k}-s_{k}^{0}\|_{L^{p}(L^{p})}$
$\displaystyle\qquad+\|\tfrac{1}{\hat{\phi}_{k}^{0}}\|_{L^{\infty}(L^{\infty})}^{2}L^{1/p}\|s_{k}-s_{k}^{0}\|_{L^{p}(L^{p})}\|\hat{\phi}_{k}^{0}\,\partial_{z}(\hat{\phi}_{k}-\hat{\phi}_{k}^{0})-(\hat{\phi}_{k}-\hat{\phi}_{k}^{0})\,\partial_{z}\hat{\phi}_{k}^{0}\|_{L^{p*}(L^{\infty})}$

and similarly for $\|\partial_{z}(\underline{\delta
s}_{1}-(s_{1}-s_{1}^{0})\|_{L^{p}(L^{p})}$, using the identity
$\partial_{z}(\underline{\delta
s}_{1}-(s_{1}-s_{1}^{0})=2(s_{1}-s_{1}^{0})\,\partial_{z}(\frac{\hat{\phi}_{k}}{\hat{\phi}^{0}_{k}})+\frac{\hat{\phi}_{k}-\hat{\phi}^{0}_{k}}{\hat{\phi}_{k}^{0}}\,\partial_{z}(s_{1}-s_{1}^{0})+\partial_{z}(\frac{1}{(\hat{\phi}_{k}^{0})^{2}})\,\int_{0}^{z}(s_{k}-s_{k}^{0})(z^{\prime})\Bigl{(}\hat{\phi}_{k}^{0}\,\partial_{z}\hat{\phi}_{k}-\hat{\phi}_{k}\,\partial_{z}\hat{\phi}_{k}^{0}\Bigr{)}(z^{\prime})\,dz^{\prime}$,
as well as

$\displaystyle\|\underline{\delta\eta}-(\eta-\eta^{0})\|_{L^{p}(L^{p})}\leq\|\tfrac{1}{\hat{\phi}_{1}^{0}\,\overline{\hat{\phi}_{2}^{0}}}\|_{L^{\infty}(L^{\infty})}\
\cdot\Bigl{\\{}\|(\hat{\phi}_{1}-\hat{\phi}_{1}^{0})\overline{\hat{\phi}_{2}^{0}}+\hat{\phi}_{1}(\overline{\hat{\phi}_{2}}-\overline{\hat{\phi}_{2}^{0}})\|_{L^{\infty}(L^{\infty})}\|\eta-\eta^{0}\|_{L^{p}(L^{p})}$
$\displaystyle\hskip
113.81102pt+\frac{\omega_{d}\,\tilde{\varepsilon}}{\omega_{p}^{3}}\Bigl{(}\|(\hat{\psi}-\hat{\psi}^{0})\,\partial_{z}(s_{1}-s_{1}^{0})\|_{L^{p}(L^{p})}+2\|(s_{1}-s_{1}^{0})\,(\partial_{z}\hat{\psi}-\partial_{z}\hat{\psi}^{0})\|_{L^{p}(L^{p})}$
$\displaystyle\hskip
142.26378pt+\|\hat{\psi}^{0}\,\partial_{z}(\underline{\delta
s}_{1}-(s_{1}-s_{1}^{0})\|_{L^{p}(L^{p})}+2\|(\underline{\delta
s}_{1}-(s_{1}-s_{1}^{0}))\,\partial_{z}\hat{\psi}^{0}\|_{L^{p}(L^{p})}\Bigr{)}\Bigr{\\}}$

hence

$\displaystyle\|r(\vec{q})-(\vec{q}-\vec{q}^{0})\|_{\mathbb{X}}\leq
C\|\vec{q}-\vec{q}^{0}\|_{\mathbb{X}}^{2}.$ (51)

Based on the range invariance condition (50) we can rewrite the inverse
problem of reconstructing $(s,\eta)$ as a combination of an ill-posed linear
and a well-posed nonlinear problem

$\displaystyle\mathbb{F}^{\prime}(\vec{q}^{0})\hat{r}=y-\mathbb{F}(\vec{q}^{0}),$
(52) $\displaystyle r(\vec{q})=\hat{r},$ $\displaystyle\mathcal{P}\vec{q}=0$

for the unknowns $(\hat{r},\vec{q})$. Here $\vec{q}^{0}$ is fixed and
$U\subseteq\mathbb{X}$ is a neighborhood of $\vec{q}^{0}$.

For its regularized iterative solution we consider the frozen Newton method

$\vec{q}^{n+1}\in\mbox{argmin}_{\vec{q}\in
U}\|\mathbb{F}^{\prime}(\vec{q}^{0})(\vec{q}-\vec{q}^{n})+\mathbb{F}(\vec{q}^{n})-y^{\delta}\|_{\mathbb{Y}}^{2}+\alpha_{n}\|\vec{q}-\vec{q}^{0}\|_{\mathbb{X}}^{2}+\|\mathcal{P}\vec{q}\|_{\mathbb{Z}}^{2},$
(53)

where $y^{\delta}\approx y$ is the noisy data, $\alpha_{n}\to 0$ as
$n\to\infty$ and $\mathbb{Z}=(L^{p}(0,L;L^{p}(\tilde{\Omega}_{y})))$. To work
in Hilbert spaces, we set $p=2$.

In order to prove convergence of (53), we additionally need a condition on
compatibility of the linear operators $\mathbb{F}^{\prime}(\vec{q}^{0})$ and
$\mathcal{P}$, see [9, Lemma B.1]. Sufficient for this condition is linearized
uniqueness with $s_{1}=s_{2}$, that is, triviality of the intersection of the
nullspaces

$\mathcal{N}(\mathbb{F}^{\prime}(\vec{q}^{0}))\,\cap\,\mathcal{N}(\mathcal{P})\
=\\{0\\}.$ (54)

For the model (17), (18) described by (47) this would likely require more than
one excitation as well as observations at several frequencies, along with the
corresponding extension of dependency of $s$ in order to allow for range
invariance (50).

We here rigorously establish (54) for the alternative model (34), (35) using
the Helmholtz equation for the outgoing wave $\hat{\psi}$ and making the
simplifying but still very realistic assumption that propagation of the
excitation waves described by $\hat{\phi}_{1}$, $\hat{\phi}_{2}$ takes place
in a homogeneous domain with known and constant sound speed.

The unknown parameters are then the squared slowness $\check{s}$ and the
nonlinearity coefficient $\check{\eta}$ and the model reads as

$-\omega_{d}^{2}\mathcal{M}[\check{s}\check{\psi}]+\mathcal{A}\check{\psi}+\imath\omega_{d}\mathcal{D}\check{\psi}=\imath\omega_{p}^{3}\,\check{\eta}\,f,$

where $f$ is given by the excitation waves $\hat{\phi}_{1}$, $\hat{\phi}_{2}$
that can be precomputed under the assumption of $s=s_{0}=\frac{1}{c}$ being
known in the first equations of (34). Moroever, on the boundary
$\partial\Omega$ we assume to know $\check{s}=s_{0}^{2}=\frac{1}{c^{2}}$, and
define the operators $\mathcal{A}$, $\mathcal{D}$, $\mathcal{M}$ by

$\displaystyle\mathcal{A}u=\Bigl{(}v\mapsto\int_{\Omega}\nabla u\cdot\nabla
v\,dx+\int_{\partial\Omega}\beta\,u\,v\,ds\Bigr{)},$ (55)
$\displaystyle\mathcal{D}u=\Bigl{(}v\mapsto b\int_{\Omega}\nabla u\cdot\nabla
v\,dx+\int_{\partial\Omega}(\sigma+b\beta)u\,v\,ds\Bigr{)}\,,$
$\displaystyle\mathcal{M}u=\Bigl{(}v\mapsto\int_{\Omega}u\,v\,dx+\int_{\partial\Omega}\frac{\sigma
b}{s_{0}^{2}}\,u\,v\,ds\Bigr{)}.$

In our linearized uniqueness proof we will assume that $\mathcal{A}$,
$\mathcal{D}$ and $\mathcal{M}$ are simultaneously diagonalizable. This holds
true, e.g., in the case $\sigma=0$, where $\mathcal{M}=\textup{id}$ and
$\mathcal{D}=b\mathcal{A}$. Note that in the all-at-once setting considered
here, we are not bound to well-posedness theory of the underlying PDE problems
and therefore have more freedom in choosing the coefficients in the boundary
conditions.

To obtain uniqueness, we will need observations on an interval $I$ of
difference frequencies $\omega_{d}$ and with at least two different sets of
excitation frequencies, more precisely

$\displaystyle\omega_{2}=\frac{\omega_{1}\tilde{\omega}}{\omega_{1}+\tilde{\omega}}\,\tilde{\omega},\quad\tilde{\omega}\in\\{\tilde{\omega}_{1},\tilde{\omega}_{2}\\}\subseteq\mathbb{R}^{+},\quad\omega_{d}=\omega_{1}-\omega_{2}\in
I\subseteq\mathbb{R}^{+}$ (56) with $I$ countable containing at least one
accumumation point

which implies $\omega_{p}^{3}=\omega_{d}^{2}\tilde{\omega}$. Correspondingly,
we consider $\check{\psi}=\check{\psi}(\omega_{d},\tilde{\omega})$ as a
function of $\omega_{d}$, while the original coefficients $\check{s}$,
$\check{\eta}$ of course stay independent of $(\omega_{d},\tilde{\omega})$. In
order to satisfy range invariance, we will introduce an artificial dependence
$\check{s}=\check{s}(\omega_{d},\tilde{\omega})$, while keeping $\check{\eta}$
independent of $\omega_{d}$ and $\tilde{\omega}$.

Therewith, the inverse problem reads as

$\mathbb{F}(\check{s},\check{\eta},\check{\psi})=y$

with

$\displaystyle\mathbb{F}_{1}(\vec{q})=-\omega_{d}^{2}\,\bigl{(}\mathcal{M}(\check{s}\,\check{\psi})+\imath\tilde{\omega}\,\check{\eta}\,f\bigr{)}+\mathcal{A}\check{\psi}+\imath\omega_{d}\mathcal{D}\check{\psi},$
$\displaystyle\mathbb{F}_{2}(\vec{q})=\imath\omega_{d}\textup{tr}_{\check{\Gamma}}\check{\psi},$

where $\vec{q}=(\check{s},\check{\eta},\check{\psi})$, and
$y=(0,p^{\textup{meas}})$.

It is readily checked that range invariance (50) holds with

$r(\check{s},\check{\eta},\check{\psi})=\Bigl{(}(\check{s}-\check{s}_{0})\frac{\check{\psi}}{\check{\psi}_{0}},\check{\eta}-\check{\eta}_{0},\check{\psi}-\check{\psi}_{0}\Bigr{)}$

and

$\displaystyle\|r(\vec{q})-(\vec{q}-\vec{q}^{0})\|_{\mathbb{X}}\leq\left\|\frac{1}{\check{\psi}_{0}}\right\|_{L^{\infty}(\Omega)}\cdot\|(\check{s}-\check{s}_{0})(\check{\psi}-\check{\psi}_{0})\|_{L^{2}(\Omega)}\leq
C\|\vec{q}-\vec{q}^{0}\|_{\mathbb{X}}^{2}\leq
c\|\vec{q}-\vec{q}^{0}\|_{\mathbb{X}}$ (57)

with $\mathbb{X}=L^{2}(\Omega)\times L^{2}(\Omega)\times H^{2}(\Omega)$,
provided $\|\frac{1}{\check{\psi}_{0}}\|_{L^{\infty}(\Omega)}<\infty$, that
is, $\check{\psi}_{0}$ is bounded away from zero. Here $c$ can be made small
by restricting $\vec{q}$ to a sufficiently small neighborhood of
$\vec{q}^{0}$.

Defining $\mathcal{P}$, $\mathbb{X}$, $\mathbb{Y}$, $\mathbb{Z}$ by

$\displaystyle\mathcal{P}(\vec{q})(\omega_{d}):=\Bigl{(}\check{s}(\omega_{d},\tilde{\omega}_{1})-\check{s}(\omega_{d},\tilde{\omega}_{2}),\
\check{s}(\omega_{d},\tilde{\omega}_{1})-\frac{1}{\mu(I)}\int_{I}\check{s}(\cdot,\tilde{\omega}_{1})\,d\mu\Bigr{)},$
(58) $\displaystyle\mathbb{X}=L^{2}_{\mu}(I;L^{2}(\Omega))^{2}\times
L^{2}(\Omega)\times
H^{2}(\Omega),\quad\mathbb{Y}=L^{2}_{\mu}(I;L^{2}(\Gamma))^{2},\quad\mathbb{Z}=L^{2}_{\mu}(I;L^{2}(\Omega))^{2}$

for some finite measure $\mu$ on $I$, we can write the inverse problem as (52)
and use (53) for its regularized numerical solution.

To prove convergence of (53) we require linearized uniqueness (54), which we
verify as follows. For
$(\underline{\delta\check{s}},\underline{\delta\check{\eta}},\underline{\delta\check{\psi}})\in\mathcal{N}(\mathcal{P})$,
that is, $\check{s}$ independent of $\omega_{d}$ and $\tilde{\omega}$, the
assumption
$\mathbb{F}^{\prime}(\check{s}_{0},\check{\eta}_{0},\check{\psi}_{0})(\underline{\delta\check{s}},\underline{\delta\check{\eta}},\underline{\delta\check{\psi}})=0$
with $\check{s}_{0}:\equiv\frac{1}{c^{2}}$ and $\check{\eta}_{0}$ arbitrary
reads as

$\displaystyle(-\omega_{d}^{2}\,\tfrac{1}{c^{2}}\mathcal{M}+\mathcal{A}+\imath\omega_{d}\mathcal{D})\,\underline{\delta\check{\psi}}(\omega_{d},\tilde{\omega})$
$\displaystyle=\omega_{d}^{2}\bigl{(}\mathcal{M}[\underline{\delta\check{s}}\,\check{\psi}_{0}(\omega_{d},\tilde{\omega})]+\imath\tilde{\omega}\,\underline{\delta\check{\eta}}\,f(\omega_{d},\tilde{\omega})\bigr{)},$
$\displaystyle\imath\omega_{d}\textup{tr}_{\check{\Gamma}}\underline{\delta\check{\psi}}(\omega_{d},\tilde{\omega})$
$\displaystyle=0$ $\displaystyle\hskip 71.13188pt\textup{for all
}\omega_{d}\in I,\
\tilde{\omega}\in\\{\tilde{\omega}_{1},\tilde{\omega}_{2}\\}.$

We now assume that $\check{\psi}_{0}$ and $f$ have been chosen independent of
$\omega_{d}$ and $\tilde{\omega}$

$\check{\psi}_{0}(\omega_{d},\tilde{\omega})\equiv\check{\psi}_{0},\quad
f(\omega_{d},\tilde{\omega})\equiv f\textup{ with }\check{\psi}\not=0,\
f\not=0\textup{ a.e.}$ (59)

and diagonalize the operators $\mathcal{A}$, $\mathcal{D}$, $\mathcal{M}$, by
means of their eigensystems
$(\lambda_{j},\varphi_{j}^{k})_{j\in\mathbb{N},k\in K^{j}}$,
$(\rho_{j},\varphi_{j}^{k})_{j\in\mathbb{N},k\in K^{j}}$
$(\mu_{j},\varphi_{j}^{k})_{j\in\mathbb{N},k\in K^{j}}$ — recall that we have
assumed them to be jointly diagonalizable and that they are symmetric and
positive semidefinite. Note that with $s_{0}>0$, $\sigma,\,\beta,\,b\geq 0$,
we have

$\mu_{j}\geq 1,\quad\lambda_{j},\,\rho_{j}\geq
0,\quad\lambda_{j}\to\infty\textup{ as }j\to\infty.$

We will additionally assume

$\Bigl{(}\frac{\rho_{j}}{\mu_{j}}=\frac{\rho_{\ell}}{\mu_{\ell}}\textup{ and
}\frac{\lambda_{j}}{\mu_{j}^{2}}=\frac{\lambda_{\ell}}{\mu_{\ell}^{2}}\Bigr{)}\
\Rightarrow\ j=\ell,$ (60)

which is the case, e.g., if $\sigma=0$, since then $\mu_{j}\equiv 1$.

This yields

$\displaystyle
0=\underline{\delta\check{\psi}}(x_{0};\omega_{d},\tilde{\omega})=c^{2}\omega_{d}^{2}\sum_{j=1}^{\infty}\sum_{k\in
K^{j}}\frac{1}{-\omega_{d}^{2}\,\mu_{j}+c^{2}\,\lambda_{j}+\imath\omega_{d}c^{2}\,\rho_{j}}(a_{j}^{k}+\imath\tilde{\omega}b_{j}^{k})\,\varphi_{j}^{k}(x_{0})$
(61) $\displaystyle\textup{ for all }x_{0}\in\Gamma,\ \omega_{d}\in I,\
\tilde{\omega}\in\\{\tilde{\omega}_{1},\tilde{\omega}_{2}\\}$

for

$a_{j}^{k}=\mu_{j}\langle\underline{\delta\check{s}}\,\check{\psi}_{0},\varphi_{j}^{k}\rangle_{L^{2}(\Omega)},\quad
b_{j}^{k}=\langle\underline{\delta\check{\eta}}\,f,\varphi_{j}^{k}\rangle_{L^{2}(\Omega)}.$

It is straightforward to show that the functions
$\omega_{d}\mapsto\frac{1}{w_{j}(\imath\omega_{d})}$ with
$w_{j}(\imath\omega_{d}):=-\omega_{d}^{2}\,\mu_{j}+c^{2}\,\lambda_{j}+\imath\omega_{d}c^{2}\,\rho_{j}$
are linearly independent on $I$. Indeed, assuming
$\sum_{j=1}^{\infty}\frac{1}{w_{j}(\imath\omega_{d})}\,c_{j}=0$ for all
$\omega_{d}\in I$ and a sequence $\vec{c}=(c_{j})_{j\in\mathbb{N}}$, after
multiplication with $\prod_{\ell\in\mathbb{N}}w_{\ell}(\imath\omega_{d})$, we
conclude
$0=W^{\vec{c}}(z):=\sum_{j=1}^{\infty}\prod_{\ell\not=j}w_{\ell}(z)\,c_{j}$
for all $z\in\imath I$ and thus, since the function $W^{\vec{c}}$ is analytic,
it vanishes on all of $\mathbb{C}$. Inserting the roots
$z^{k}_{\pm}=-\frac{c\rho_{j}}{2}\pm\sqrt{\frac{c^{2}\rho_{j}^{2}}{4}-\lambda_{j}}$
of $w_{k}$ and using the fact that under condition (60) they are distinct (in
the sense that $z^{j}_{+}=z^{\ell}_{+}$ and $z^{j}_{-}=z^{\ell}_{-}$ implies
$j=\ell$), we obtain $\vec{c}=0$.
By the linear independence of $(\frac{1}{w_{j}})_{j\in\mathbb{N}}$ we conclude
from (61)

$\sum_{k\in
K^{j}}(a_{j}^{k}+\imath\tilde{\omega}b_{j}^{k})\,\varphi_{j}^{k}(x_{0})=0\quad\textup{
for all }x_{0}\in\Gamma,\ j\in\mathbb{N},\
\tilde{\omega}\in\\{\tilde{\omega}_{1},\tilde{\omega}_{2}\\}.$

Under a linear independence assumption on the individual eigenspaces (cf.,
e.g., [12])

$\textup{ for all }j\in\mathbb{N}\ :\quad\left(\sum_{k\in
K^{j}}\gamma_{k}\varphi_{j}^{k}(x)=0\ \mbox{ for all }x\in\Gamma\right)\
\Longrightarrow\ \left(\gamma_{k}=0\mbox{ for all }k\in K^{j}\right)$ (62)

we conclude

$\begin{cases}a_{j}^{k}+\imath\tilde{\omega}_{1}b_{j}^{k}=0\\\
a_{j}^{k}+\imath\tilde{\omega}_{2}b_{j}^{k}=0\end{cases}\quad\textup{ for all
}j\in\mathbb{N},\ k\in K^{j},$

which implies $a_{j}^{k}=0$, $b_{j}^{k}=0$, for all $j\in\mathbb{N}$, $k\in
K^{j}$, thus under assumption (59), we get $\underline{\delta\check{s}}=0$,
$\underline{\delta\check{\eta}}=0$.

According to [9, Theorem 2], we obtain the following.

######  Theorem 1.

Let $\vec{q}^{\dagger}=(\check{s},\check{\eta},\check{\psi})$ be a solution to
(52) and let for the noise level $\delta\geq\|y^{\delta}-y\|_{\mathbb{Y}}$ the
stopping index $n_{*}=n_{*}(\delta)$ be chosen such that

$n_{*}(\delta)\to
0,\quad\delta\sum_{j=0}^{n_{*}(\delta)-1}c^{j}\alpha_{n_{*}(\delta)-j-1}^{-1/2}\to
0\qquad\textup{ as }\delta\to 0$ (63)

with $c$ as in (57). Moreover, let (59) with $\frac{1}{\check{\psi}_{0}}\in
L^{\infty}(\Omega)$ hold and let the operators $\mathcal{A}$, $\mathcal{D}$,
$\mathcal{M}$ be jointly diagonalizable with eigenvalues satisfying (60).

Then there exists $\varrho>0$ sufficiently small such that for
$\vec{q}_{0}\in\mathcal{B}_{\varrho}(\vec{q}^{\dagger})$ the iterates
$(\vec{q}_{n}^{\delta})_{n\in\\{1,\ldots,n_{*}(\delta)\\}}$ are well-defined
by (53), remain in $\mathcal{B}_{\varrho}(\vec{q}^{\dagger})$ and converge in
$\mathbb{X}$,
$\|\vec{q}_{n_{*}(\delta)}^{\delta}-\vec{q}^{\dagger}\|_{\mathbb{X}}\to 0$ as
the noise level $\delta\to 0$. In the noise free case $\delta=0$,
$n_{*}(\delta)=\infty$ we have
$\|\vec{q}_{n}-\vec{q}^{\dagger}\|_{\mathbb{X}}\to 0$ as $n\to\infty$.

This work was supported by the Austrian Science Fund fwf
http://dx.doi.org/10.13039/501100002428 under the grant doc78.

## Bibliography

  * [1] Sebastian Acosta, Gunther Uhlmann, and Jian Zhai.  Nonlinear ultrasound imaging modeled by a Westervelt equation.  SIAM Journal on Applied Mathematics, 82(2):408–426, 2022.
  * [2] David Colton and Rainer Kress.  Inverse acoustic and electromagnetic scattering theory, volume 93 of Applied Mathematical Sciences.  Springer, New York, third edition, 2013.
  * [3] Lawrence C. Evans.  Partial Differential Equations.  Graduate studies in mathematics. American Mathematical Society, Providence (R.I.), 1998.
  * [4] Mostafa Fatemi and James F Greenleaf.  Ultrasound-stimulated vibro-acoustic spectrography.  Science, 280:82–85, 1998.
  * [5] Mostafa Fatemi and James F. Greenleaf.  Vibro-acoustography: An imaging modality based on ultrasound-stimulated acoustic emission.  Proceedings of the National Academy of Sciences, 96(12):6603–6608, 1999.
  * [6] Bernd Flemisch, Manfred Kaltenbacher, and Barbara I. Wohlmuth.  Elasto-Acoustic and Acoustic-Acoustic Coupling on Nonmatching Grids.  Int. J. Numer. Meth. Engng, 67(13):1791–1810, 2006.
  * [7] Ivan G. Graham and Stefan A. Sauter.  Stability and finite element error analysis for the Helmholtz equation with variable coefficients.  Math. Comp., 89(321):105–138, 2020.
  * [8] Barbara Kaltenbacher.  Periodic solutions and multiharmonic expansions for the Westervelt equation.  Evolution Equations and Control Theory EECT, 2020.  19 pages.
  * [9] Barbara Kaltenbacher.  Convergence guarantees for coefficient reconstruction in PDEs from boundary measurements by variational and Newton type methods via range invariance.  2022\.  submitted; see also arXiv:2209.12596 [math.NA].
  * [10] Barbara Kaltenbacher.  On the inverse problem of vibro-acoustography.  Meccanica, 2022.  see also arXiv:2109.01907 [math.AP].
  * [11] Barbara Kaltenbacher and William Rundell.  On the identification of the nonlinearity parameter in the westervelt equation from boundary measurements.  Inverse Problems & Imaging, 15:865–891, 2021.
  * [12] Barbara Kaltenbacher and William Rundell.  On an inverse problem of nonlinear imaging with fractional damping.  Mathematics of Computation, 91:245–276, 2022.  see also arXiv:2103.08965 [math.AP].
  * [13] Barbara Kaltenbacher and William Rundell.  On the simultanenous identification of two space dependent coefficients in a quasilinear wave equation.  2022\.  submitted; see also arXiv:2210.08063 [math.NA].
  * [14] Carlos E. Kenig.  Lecture notes : Global Well-posedness, scattering and blow up for the energy-critical, focusing, non-linear Schrödinger and wave equations.  Journées équations aux dérivées partielles, 2007.
  * [15] Alison E. Malcolm, Fernando Reitich, J. Yang, J.F. Greenleaf, and Mostafa Fatemi.  A combined parabolic-integral equation approach to the acoustic simulation of vibro-acoustic imaging.  Ultrasonics, 48:553–558, 2008.
  * [16] Alison E. Malcolm, Fernando Reitich, Jiaqi Yang, James F. Greenleaf, and Mostafa Fatemi.  Numerical modeling for assessment and design of ultrasound vibro-acoustography systems.  In Biomedical Applications of Vibration and Acoustics for Imaging and Characterizations. ASME Press, New York, 2007.
  * [17] Jens Markus Melenk.  On generalized finite-element methods.  ProQuest LLC, Ann Arbor, MI, 1995.  Thesis (Ph.D.)–University of Maryland, College Park.
  * [18] Teresa Rauscher.  A paraxial approach for the inverse problem of vibroacoustic imaging in frequency domain.  in preparation.
  * [19] Anna Rozanova.  The Khokhlov-Zabolotskaya-Kuznetsov equation.  C. R. Math. Acad. Sci. Paris, 344(5):337–342, 2007.
  * [20] Masahiro Yamamoto and Barbara Kaltenbacher.  An inverse source problem related to acoustic nonlinearity parameter imaging.  In Barbara Kaltenbacher, Anne Wald, and Thomas Schuster, editors, Time-dependent Problems in Imaging and Parameter Identification. Springer, New York, 2021.
  * [21] E. A. Zabolotskaya and R. V. Khokhlov.  Quasi-plane waves in the non-linear acoustics of confined beams.  Soviet Physics-Acoustics, 15:35–40, 1969.