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
ω
1
and
ω
2
are excited at two parts
Σ
1
,
Σ
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
ω
1
−
ω
2
(which allows to avoid strong scattering and attenuation).
(iii) The latter is eventually measured by a receiver array Γ.
We refer to, e.g., [11, 20] for a graphical illustration.
The fact that the value of the nonlinearity parameter
B
A
governing the interaction depends on the tissue properties and thus varies in space
B
A
=
B
A
(
x
)
yields a means of imaging.
Likewise, also the speed of sound typically exhibits spatial variation
c
=
c
(
x
)
.
A modeling and simulation framework for this imaging technology has been devised in [18, 17].
The aim of this paper is to study the inverse problem of identifying both
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 [11], 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 𝜙 given by
(1.1)
∂
t
2
ϕ
−
c
2
Δ
ϕ
=
∂
t
|
∇
ϕ
|
2
−
1
2
(
∇
ϕ
)
⋅
∇
|
∇
ϕ
|
2
+
B
A
∂
t
ϕ
Δ
ϕ
−
B
2
A
|
∇
ϕ
|
2
Δ
ϕ
,
by means of the asymptotic ansatz
(1.2)
ϕ
̃
(
t
,
x
,
ε
)
=
ε
(
ϕ
̃
1
(
t
,
x
)
+
ϕ
̃
2
(
t
,
x
)
)
+
ε
2
ψ
̃
(
t
,
x
)
,
with small
ε
>
0
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.1), 𝑐 is the speed of sound and
B
A
the nonlinearity parameter.
First of all, plugging (1.2) into (1.1) yields
ε
∂
t
2
(
ϕ
1
+
ϕ
2
)
+
ε
2
∂
t
2
ψ
−
c
2
ε
Δ
(
ϕ
1
+
ϕ
2
)
−
c
2
ε
2
Δ
ψ
=
ε
2
∂
t
|
∇
(
ϕ
1
+
ϕ
2
)
|
2
+
ε
2
B
A
∂
t
(
ϕ
1
+
ϕ
2
)
Δ
(
ϕ
1
+
ϕ
2
)
+
O
(
ε
3
)
.
Considering only terms up to order
ε
2
and rearranging the equation above results in
(1.3)
ε
[
∂
t
2
(
ϕ
1
+
ϕ
2
)
−
c
2
Δ
(
ϕ
1
+
ϕ
2
)
−
ε
∂
t
|
∇
(
ϕ
1
+
ϕ
2
)
|
2
−
ε
B
A
∂
t
(
ϕ
1
+
ϕ
2
)
Δ
(
ϕ
1
+
ϕ
2
)
]
=
ε
2
[
−
∂
t
2
ψ
+
c
2
Δ
ψ
]
,
an expression containing the second order in space and time wave operator
∂
t
2
−
c
2
Δ
as well as further nonlinear terms; separating contributions that are asymptotic to different powers of 𝜀 (namely
ε
1
and
ε
2
in our case) allows to split (1.3) into separate equation for
ϕ
1
,
ϕ
2
and 𝜓.
However, these 𝜀 asymptotics interact with another scaling that results from the geometry of wave propagation as follows.
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
′
)
, where
x
′
=
(
x
2
,
…
,
x
d
)
, (and marking functions of these variables by a check), we perform the paraxial change of variables
(1.4)
(
τ
,
z
,
y
)
=
(
t
−
T
̌
,
ε
̃
x
1
,
ε
̃
x
′
)
,
with
T
̌
=
T
̌
(
x
1
,
x
′
)
,
T
=
T
(
z
,
y
)
chosen such that
1
c
2
=
|
∇
x
T
̌
|
2
=
|
∂
x
1
T
̌
|
2
+
|
∇
x
′
T
̌
|
2
=
ε
̃
2
|
∂
z
T
|
2
+
ε
̃
|
∇
y
T
|
2
,
T
̌
(
0
,
x
′
)
=
0
,
x
′
∈
Ω
̃
y
,
where
ε
̃
≪
1
, so
ε
̃
≪
ε
̃
, and in general,
ε
̃
can be different from 𝜀 that was introduced for the second order approximation in (1.2), but for the model to make sense, the relation
ε
≤
ε
̃
≤
ε
should hold; cf. [20].
Concerning well-definedness of
T
̌
, we note that it is actually the distance of 𝑥 from the boundary
{
0
}
×
Ω
̃
y
in the Riemannian metric determined by
1
c
2
; cf., e.g. [3, Section 7.2.3].
For the derivatives, this change of variables yields
(1.5)
∂
t
=
∂
τ
,
∇
x
=
(
ε
̃
(
−
∂
z
T
∂
τ
+
∂
z
)
,
ε
̃
(
−
∇
y
T
∂
τ
+
∇
y
)
)
,
Δ
x
=
−
ε
̃
2
∂
z
2
T
∂
τ
+
ε
̃
2
|
∂
z
T
|
2
∂
τ
2
−
2
ε
̃
2
∂
z
T
∂
z
∂
τ
+
ε
̃
2
∂
z
2
−
ε
̃
Δ
y
T
∂
τ
+
ε
̃
|
∇
y
T
|
2
∂
τ
2
−
2
ε
̃
∇
y
T
⋅
∇
y
∂
τ
+
ε
̃
Δ
y
=
1
c
2
∂
τ
2
+
ε
̃
(
−
(
ε
̃
∂
z
2
T
+
Δ
y
T
)
∂
τ
−
2
(
ε
̃
∂
z
T
∂
z
+
∇
y
T
⋅
∇
y
)
∂
τ
+
ε
̃
∂
z
2
+
Δ
y
)
.
Here
∇
x
=
(
∂
x
1
,
∇
x
′
)
and
(
∂
z
,
∇
y
)
are the nabla operators with respect to the original and paraxial variables, respectively.
Likewise for
Δ
x
, etc.
Moreover, in (1.4) and (1.5), we have used the fact that
∂
t
s
̌
=
0
,
∂
t
T
̌
=
0
; hence
∂
τ
s
=
0
,
∂
τ
T
=
0
,
∂
x
1
s
̌
=
ε
̃
∂
z
s
,
∂
x
1
T
̌
=
ε
̃
∂
z
T
,
∇
x
′
T
̌
=
ε
̃
∇
y
T
.
Frequency Domain Formulation
We make the time harmonic ansatz
ϕ
k
(
τ
,
z
,
y
)
=
ϕ
̂
k
(
z
,
y
)
e
ı
ω
k
τ
for
k
∈
{
1
,
2
}
,
ψ
(
τ
,
z
,
y
)
=
ℜ
{
ψ
̂
(
z
,
y
)
e
ı
(
ω
1
−
ω
2
)
τ
}
.
Here the excitation frequencies are assumed to satisfy
ω
1
>
ω
2
>
0
, with a comparably small difference frequency
ω
d
:=
ω
1
−
ω
2
>
0
, and we abbreviate
ω
p
3
:=
ω
1
ω
2
ω
d
.
(Typical values in ultrasonics are
ω
k
∼
2
MHz,
ω
d
∼
50
kHz.)
Here
ı
=
−
1
denotes the complex unit and
ϕ
̂
1
,
ϕ
̂
2
,
ψ
̂
are complex valued in general.
Altogether, for
ψ
̃
(
t
,
x
1
,
x
′
)
=
ψ
(
τ
,
z
,
y
)
=
ℜ
(
ψ
̂
(
z
,
y
)
e
ı
ω
d
τ
)
,
ψ
̃
(
t
,
x
1
,
x
′
)
=
ℜ
(
ψ
̌
(
x
1
,
x
′
)
e
ı
ω
d
t
)
with
e
ı
ω
d
t
=
e
ı
ω
d
x
1
/
c
e
ı
ω
d
τ
,
we end up with the transformation in particular of the wave operator
(1.6)
ℜ
(
[
−
1
c
2
ω
d
2
−
Δ
x
]
ψ
̌
(
x
1
,
x
′
)
e
ı
ω
d
T
e
ı
ω
d
τ
)
=
[
1
c
2
∂
t
2
−
Δ
x
]
ψ
̃
(
t
,
x
1
,
x
′
)
=
ε
̃
ℜ
(
[
ı
ω
d
(
ε
̃
∂
z
2
T
+
Δ
y
T
)
+
2
ı
ω
d
(
ε
̃
∂
z
T
∂
z
+
∇
y
T
⋅
∇
y
)
−
ε
̃
∂
z
2
−
Δ
y
]
ψ
̂
(
z
,
y
)
e
ı
ω
d
τ
)
and likewise for
ϕ
̃
k
,
ϕ
̂
k
.
Due to the boundary condition
T
̌
(
0
,
x
′
)
=
0
,
x
′
∈
Ω
̃
y
, we also have
∇
y
T
(
0
,
x
′
)
and
Δ
y
T
(
0
,
x
′
)
, which allows us to neglect the terms containing
∇
y
T
and
Δ
y
T
to arrive at
[
−
1
c
2
ω
d
2
−
Δ
x
]
ψ
̌
(
x
1
,
x
′
)
e
ı
ω
d
T
≈
ε
̃
[
ı
ω
d
∂
z
s
+
2
ı
ω
d
s
∂
z
−
ε
̃
∂
z
2
−
Δ
y
]
ψ
̂
(
z
,
y
)
.
The outcome of the nonlinear term is more complicated and depends on the interplay between the asymptotic ansatz (1.2) governed by 𝜀 and the paraxial approximation (1.4) governed by
ε
̃
.
Depending on the relation between 𝜀 and
ε
̃
, this leads to several different models; see [20] for more details.
The Considered Models
The inverse problem of combined nonlinearity imaging and speed of sound reconstruction will be considered for one of these models, namely
(1.7)
ı
ω
k
∂
z
s
ϕ
̂
k
+
2
ı
ω
k
s
∂
z
ϕ
̂
k
−
Δ
y
ϕ
̂
k
=
0
in
Ω
̃
,
k
∈
{
1
,
2
}
,
ı
ω
d
∂
z
s
ψ
̂
+
2
ı
ω
d
s
∂
z
ψ
̂
−
ε
̃
∂
z
2
ψ
̂
−
Δ
y
ψ
̂
=
ı
ω
p
3
ε
̃
−
1
η
ϕ
̂
1
ϕ
̂
̄
2
in
Ω
̃
,
with boundary conditions
(1.8)
ϕ
̂
k
(
0
,
y
)
=
h
k
(
y
)
,
y
∈
Ω
̃
y
,
∂
ν
y
ϕ
̂
k
(
z
,
y
)
=
−
ı
ω
k
σ
k
ϕ
̂
k
(
z
,
y
)
,
(
z
,
y
)
∈
(
0
,
L
)
×
∂
Ω
̃
y
,
k
∈
{
1
,
2
}
,
−
∂
z
ψ
̂
(
0
,
y
)
=
−
ı
ω
d
σ
0
ψ
̂
(
0
,
y
)
,
∂
z
ψ
̂
(
L
,
y
)
=
−
ı
ω
d
σ
L
ψ
̂
(
L
,
y
)
,
y
∈
Ω
̃
y
,
∂
ν
y
ψ
̂
(
z
,
y
)
=
−
ı
ω
d
σ
ψ
̂
(
z
,
y
)
,
(
z
,
y
)
∈
(
0
,
L
)
×
∂
Ω
̃
y
,
and
Ω
̃
=
(
0
,
L
)
×
Ω
̃
y
.
The coefficients we are interested in are related to the speed of sound and the nonlinearity parameter by
s
=
1
c
,
η
=
B
/
A
+
2
c
4
.
Due to space dependence of 𝑐 and
B
A
, they vary in space, thus on the paraxial variables
s
=
s
(
z
,
y
)
,
η
=
η
(
z
,
y
)
.
In (1.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
Ω
̃
.
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
−
∂
z
ϕ
̂
k
(
0
,
⋅
)
=
h
k
in
Ω
̃
y
.
In an analysis analogous to the one we provided here, this could be tackled by differentiating
ϕ
̂
k
with respect to 𝑧 and considering the corresponding initial value problems for
ϕ
̂
k
′
:=
∂
z
ϕ
̂
k
.
Dirichlet conditions as used here can be justified by a surface force – pressure continuity
F
=
[
σ
]
⋅
ν
=
p
ν
together with the relation
p
k
=
ρ
∂
t
ϕ
k
derived from momentum balance.
We refer to, e.g., [6, equations (12), (13)] and the references therein for the proper choice of interface conditions in structure-acoustics coupling.
Well-posedness of the forward problem (1.7)–(1.8) will be studied in Section 2.
From the point of view of the outgoing wave described by
ψ
̂
, the parameter
ε
̃
is not necessarily small since propagation of sound is non-directed for the 𝜓 field.
Therefore, alternatively to (1.7), (1.8), we will consider
(1.9)
ı
ω
k
∂
z
s
ϕ
̂
k
+
2
ı
ω
k
s
∂
z
ϕ
̂
k
−
Δ
y
ϕ
̂
k
=
0
in
Ω
̃
,
k
∈
{
1
,
2
}
,
−
ω
d
2
s
2
ψ
̌
−
ı
ω
d
b
Δ
x
ψ
̌
−
Δ
x
ψ
̌
=
ı
ω
p
3
1
Ω
(
η
̌
P
−
1
(
ϕ
̂
1
ϕ
̂
̄
2
)
)
in
Ω
,
for
ψ
̌
=
P
−
1
ψ
̂
,
η
̌
=
e
−
ı
ω
d
T
P
−
1
η
(where 𝑃 is defined by
(
P
ψ
̌
)
(
z
,
y
)
=
ψ
̌
(
x
1
,
x
′
)
, cf. (1.5)), with boundary conditions
(1.10)
ϕ
̂
k
(
0
,
⋅
)
=
h
k
in
Ω
̃
y
,
∂
ν
y
ϕ
̂
k
(
z
,
y
)
=
−
ı
ω
k
σ
k
ϕ
̂
k
in
(
0
,
L
)
×
∂
Ω
̃
y
,
k
∈
{
1
,
2
}
,
∂
ν
ψ
̌
=
−
ı
ω
d
σ
ψ
̌
−
β
ψ
̌
on
∂
Ω
;
see Section 2.2 for its well-posedness analysis.
Due to (1.5), the coefficients in the boundary conditions (1.8), (1.10) are related by
σ
0
(
y
)
=
σ
(
0
,
y
)
−
s
(
0
,
y
)
,
σ
L
(
y
)
=
σ
(
L
,
y
)
−
s
(
L
,
y
)
.
Thus, with the typical choice
σ
=
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
(1.11)
2
s
∂
τ
z
2
ψ
−
∇
y
2
ψ
=
η
2
∂
τ
|
∂
τ
ψ
|
2
(see [23] and also, e.g., [21] for its analysis).
Note, however, that, in our case, the quadratic nonlinearity is decoupled and appears as a source term for the 𝜓 equation, whereas (1.11) is a nonlinear (more precisely, quasilinear) equation, whose expansion in frequency domain would lead to an infinite system of space-dependent PDEs, similarly to [9].
The Inverse Problem
Our aim is to reconstruct 𝑠 and 𝜂 in the boundary value problem (1.7), (1.8) from measurements of the acoustic pressure
p
meas
=
ı
ω
d
tr
Γ
ψ
̂
in
Γ
,
where
Γ
=
P
(
Γ
̌
)
and
Γ
̌
is a manifold representing the receiver array immersed in the acoustic domain Ω.
This can be formulated as an operator equation
(1.12)
F
(
s
,
η
)
=
y
with
y
=
p
meas
and the forward operator
F
=
C
∘
S
being a concatenation of the (nonlinear) parameter-to-state map
(1.13)
S
:
(
s
,
η
)
↦
(
ϕ
̂
1
,
ϕ
̂
2
,
ψ
̂
)
solving (
1.7
), (
1.8
)
with the (linear) observation operator
(1.14)
C
:
(
ϕ
̂
1
,
ϕ
̂
2
,
ψ
̂
)
↦
ı
ω
d
tr
Γ
ψ
̂
.
We consider 𝐹 as an operator
F
:
X
s
×
X
η
→
Y
with the parameter space
X
:=
X
s
×
X
η
and the data space 𝑌 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 𝑆.
Identification of the nonlinearity coefficient 𝜂 in time domain models of nonlinear acoustics has been studied, e.g., in [1, 13, 12, 15, 8, 22]
and in frequency domain in [14]; in particular, [15] also considers simultaneous identification of 𝑠 and 𝜂.
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 [11].
In [20], models using a paraxial approximation are derived in time and frequency domain, and the inverse problem of reconstructing 𝜂 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 (1.13) for both options (1.7) and (1.9), 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 𝜂 and 𝑠 in (1.9).
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 𝑠 and 𝜂 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 [15].
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 (1.7), (1.8) and (1.9), (1.10), 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, 𝑠 and 𝜂 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 (1.12), we will also prove differentiability of the parameter-to-state map (1.13).
2.1 Well-Posedness of Paraxial Wave Propagation with Variable Coefficients
Consider the PDEs (1.7) with boundary conditions (1.8).
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., [16] and the references therein). Since we here require estimates that are explicit in terms of appropriate norms of 𝑠 and 𝜂, we will first of all provide some energy estimates for solutions of the linear variable coefficient problems
(2.1)
ı
ω
s
̃
u
+
2
ı
ω
s
∂
z
u
−
Δ
y
u
=
f
in
(
0
,
L
)
×
Ω
̃
y
,
u
(
0
,
y
)
=
h
(
y
)
,
y
∈
Ω
̃
y
,
∂
ν
y
u
(
z
,
y
)
=
−
ı
ω
σ
u
(
z
,
y
)
,
(
z
,
y
)
∈
(
0
,
L
)
×
∂
Ω
̃
y
,
and
(2.2)
ı
ω
s
̃
u
b
+
2
ı
ω
s
∂
z
u
b
−
b
∂
z
2
u
b
−
Δ
y
u
b
=
f
in
(
0
,
L
)
×
Ω
̃
y
,
−
∂
z
u
b
(
0
,
y
)
=
−
ı
ω
σ
0
u
b
(
0
,
y
)
,
∂
z
u
b
(
L
,
y
)
=
−
ı
ω
σ
L
u
b
(
L
,
y
)
,
y
∈
Ω
̃
y
,
∂
ν
y
u
b
(
z
,
y
)
=
−
ı
ω
σ
u
b
(
z
,
y
)
,
(
z
,
y
)
∈
(
0
,
L
)
×
∂
Ω
̃
y
.
Here we assume all coefficients 𝑠,
σ
0
,
σ
L
, 𝜎, 𝑏 to be real valued and positive with
s
,
1
s
∂
z
s
∈
L
∞
(
0
,
L
;
L
∞
(
Ω
̃
y
)
)
,
σ
,
1
σ
∈
L
∞
(
0
,
L
;
L
∞
(
∂
Ω
̃
y
)
)
,
σ
0
,
1
σ
0
,
σ
L
,
1
σ
L
∈
L
∞
(
Ω
̃
y
)
and 𝑏 constant.
The following two identities will be used repeatedly:
2
s
ℜ
(
∂
z
v
v
̄
)
=
s
d
d
z
|
v
|
2
=
d
d
z
(
s
|
v
|
2
)
−
∂
z
s
|
v
|
2
,
∫
Ω
̃
y
−
Δ
y
u
w
d
y
=
∫
Ω
̃
y
∇
y
u
⋅
∇
y
w
d
y
+
ı
ω
∫
∂
Ω
̃
y
σ
u
w
d
Γ
(
y
)
.
Moreover, for (2.2), we will assume the Poincaré–Friedrichs type estimate on the domain
Ω
̃
y
,
(2.3)
∥
v
∥
L
2
(
Ω
̃
y
)
2
≤
C
0
∥
σ
v
∥
L
2
(
∂
Ω
̃
y
)
2
+
C
1
∥
∇
y
v
∥
L
2
(
Ω
̃
y
)
2
for all
v
∈
H
1
(
Ω
̃
y
)
,
to hold for some
C
0
,
C
1
>
0
.
Proposition 1
If
s
∈
W
1
,
∞
(
0
,
L
;
L
∞
(
Ω
̃
y
)
)
,
∂
z
s
̃
∈
L
∞
(
Ω
̃
)
,
h
,
Δ
y
h
+
f
(
0
)
∈
L
2
(
Ω
̃
y
)
,
f
∈
L
∞
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
∩
W
1
,
1
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
,
then a solution of (2.1), exists, is unique and satisfies the estimate
∥
∇
y
u
(
z
)
∥
L
∞
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
2
+
ω
∥
s
∂
z
u
∥
L
∞
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
2
+
ω
∥
s
u
∥
L
∞
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
2
≤
C
e
L
(
s
∼
+
μ
)
(
ω
|
s
(
0
)
h
|
L
2
(
Ω
̃
y
)
2
+
|
1
2
s
(
0
)
(
Δ
y
h
+
f
(
0
)
)
|
L
2
(
Ω
̃
y
)
2
+
1
ω
∥
1
s
f
∥
L
∞
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
2
+
∥
1
s
f
∥
L
1
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
2
+
∥
1
s
∂
z
f
∥
L
1
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
2
)
with some constants
μ
>
0
,
C
>
0
independent of 𝜔.
If
s
∈
W
1
,
∞
(
0
,
L
;
L
∞
(
Ω
̃
y
)
)
,
f
∈
L
2
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
,
(2.3) holds
and the relative smallness conditions on 𝑠,
∂
z
s
, 𝜔,
(2.4)
C
0
|
s
̃
−
∂
z
s
|
L
∞
<
1
and
2
C
1
ω
2
|
s
|
L
∞
2
<
b
(
1
−
C
0
|
s
̃
−
∂
z
s
|
L
∞
)
,
are satisfied with
c
¯
, 𝜆 as in (2.3),
then a solution of (2.2) exists, is unique and satisfies the estimate
∥
u
b
∥
H
1
(
(
0
,
L
)
×
Ω
̃
y
)
2
≤
C
(
1
+
1
ω
2
)
∥
f
∥
L
2
(
0
,
L
;
L
2
(
Ω
)
)
2
with some constant
C
>
0
independent of 𝜔.
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
Ω
̃
y
, abbreviating by
|
⋅
|
L
2
(
Ω
̃
y
)
and
|
⋅
|
L
2
(
∂
Ω
̃
y
)
the
L
2
norm on
Ω
̃
y
and
∂
Ω
̃
y
, respectively.
Testing this way first of all (2.1) with
u
̄
and taking the real and imaginary parts, we obtain
(2.5)
|
∇
y
u
(
z
)
|
L
2
(
Ω
̃
y
)
2
=
2
ω
∫
Ω
̃
y
s
ℑ
(
∂
z
u
(
z
)
u
̄
(
z
)
)
)
d
y
+
∫
Ω
̃
y
ℜ
(
f
u
̄
)
d
y
,
(2.6)
ω
d
d
z
|
s
u
|
L
2
(
Ω
̃
y
)
2
(
z
)
+
ω
|
σ
u
(
z
)
|
L
2
(
∂
Ω
̃
y
)
2
=
−
ω
∫
Ω
̃
y
(
s
̃
−
∂
z
s
)
|
u
(
z
)
|
2
d
y
+
∫
Ω
̃
y
ℑ
(
f
u
̄
)
d
y
≤
ω
(
s
∼
+
μ
)
|
s
u
|
L
2
(
Ω
̃
y
)
2
(
z
)
+
1
4
μ
ω
|
1
s
f
|
L
2
(
Ω
̃
y
)
2
(
z
)
for all
z
∈
(
0
,
L
)
, with
s
∼
:=
∥
s
̃
−
∂
z
s
s
∥
L
∞
(
0
,
L
;
L
∞
(
Ω
̃
y
)
)
, where we have used Young’s inequality.
Doing the same with (2.2) and integrating over
(
0
,
L
)
, we obtain
(2.7)
∥
∇
y
u
b
∥
L
2
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
2
+
b
∥
∂
z
u
b
∥
L
2
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
2
=
2
ω
∫
0
L
∫
Ω
̃
y
s
ℑ
(
∂
z
u
b
u
̄
b
)
d
y
d
z
+
∫
0
L
∫
Ω
̃
y
ℜ
(
f
u
̄
b
)
d
y
d
z
,
(2.8)
ω
|
s
u
b
|
L
2
(
Ω
̃
y
)
2
(
L
)
+
ω
∥
σ
u
b
∥
L
2
(
0
,
L
;
L
2
(
∂
Ω
̃
y
)
)
2
+
b
ω
σ
0
|
u
b
(
0
)
|
L
2
(
Ω
̃
y
)
2
+
b
ω
σ
L
|
u
b
(
L
)
|
L
2
(
Ω
̃
y
)
2
=
ω
|
s
u
b
|
L
2
(
Ω
̃
y
)
2
(
0
)
−
ω
∫
0
L
∫
Ω
̃
y
(
s
̃
−
∂
z
s
)
|
u
b
|
2
d
y
d
z
+
∫
0
L
∫
Ω
̃
y
ℑ
(
f
u
̄
b
)
d
y
d
z
.
We also differentiate (2.1) with respect to 𝑧 and test with
∂
z
u
̄
.
Analogously to above, this yields
|
∇
y
∂
z
u
(
z
)
|
L
2
(
Ω
̃
y
)
2
=
ω
∫
Ω
̃
y
(
∂
z
s
̃
ℑ
(
∂
z
u
(
z
)
u
̄
(
z
)
)
+
2
s
ℑ
(
∂
z
2
u
(
z
)
∂
z
u
̄
(
z
)
)
)
d
y
+
∫
Ω
̃
y
ℜ
(
∂
z
f
∂
z
u
̄
)
d
y
,
(2.9)
ω
∫
Ω
̃
y
∂
z
s
̃
ℜ
(
∂
z
u
(
z
)
u
̄
(
z
)
)
d
y
+
ω
d
d
z
|
s
∂
z
u
|
L
2
(
Ω
̃
y
)
2
(
z
)
+
ω
|
σ
∂
z
u
(
z
)
|
L
2
(
∂
Ω
̃
y
)
2
≤
ω
(
s
+
+
μ
)
|
s
∂
z
u
|
L
2
(
Ω
̃
y
)
2
(
z
)
+
1
4
μ
ω
|
1
s
∂
z
f
|
L
2
(
Ω
̃
y
)
2
(
z
)
with
s
+
=
∥
s
̃
+
∂
z
s
s
∥
L
∞
(
0
,
L
;
L
∞
(
Ω
̃
y
)
)
.
For the initial value problem (with respect to 𝑧) for the Schrödinger equation (2.1), we combine
(
2.6
)
+
ρ
⋅
(
2.9
)
, which yields for
η
(
z
)
=
η
0
(
z
)
+
ρ
η
1
(
z
)
=
|
s
u
|
L
2
(
Ω
̃
y
)
2
(
z
)
+
ρ
|
s
∂
z
u
|
L
2
(
Ω
̃
y
)
2
(
z
)
that
d
d
z
η
0
≤
(
s
∼
+
μ
)
η
0
+
1
4
ω
2
μ
|
1
s
f
|
L
2
(
Ω
̃
y
)
2
,
d
d
z
η
1
≤
(
s
+
+
μ
)
η
1
+
1
4
ω
2
μ
|
1
s
∂
z
f
|
L
2
(
Ω
̃
y
)
2
+
∥
∂
z
s
̃
/
s
∥
L
∞
(
Ω
̃
)
η
0
η
1
,
thus, using Young’s inequality with factors
1
2
ρ
,
ρ
2
for the last term,
d
d
z
η
≤
β
η
+
α
,
where
β
=
max
{
s
∼
,
s
+
}
+
μ
+
ρ
/
2
∥
∂
z
s
̃
/
s
∥
L
∞
(
Ω
̃
)
,
α
(
z
)
=
1
4
ω
2
μ
(
|
1
s
f
(
z
)
|
L
2
(
Ω
̃
y
)
2
+
|
1
s
∂
z
f
(
z
)
|
L
2
(
Ω
̃
y
)
2
)
.
An application of Gronwall’s inequality yields
η
(
z
)
=
|
s
u
|
L
2
(
Ω
̃
y
)
2
(
z
)
+
ρ
|
s
∂
z
u
|
L
2
(
Ω
̃
y
)
2
(
z
)
≤
e
β
z
(
|
s
u
|
L
2
(
Ω
̃
y
)
2
(
0
)
+
ρ
|
s
∂
z
u
|
L
2
(
Ω
̃
y
)
2
(
0
)
)
∫
0
z
e
β
(
z
−
a
)
α
(
a
)
d
a
,
where we can insert the initial conditions and their differentiated version substituting from the PDE
s
(
0
)
u
(
0
)
=
s
(
0
)
h
,
s
(
0
)
∂
z
u
(
0
)
=
1
2
ı
ω
s
(
0
)
(
f
(
0
)
+
Δ
y
h
−
ı
ω
s
̃
(
0
)
h
)
.
Using this together with the Cauchy–Schwarz and Young’s inequalities in (2.5) yields the estimate
|
∇
y
u
(
z
)
|
L
2
(
Ω
̃
y
)
2
≤
2
ω
∫
Ω
̃
y
|
s
∂
z
u
|
(
z
)
|
s
u
|
(
z
)
d
y
+
∫
Ω
̃
y
|
1
s
f
|
|
s
u
|
d
y
≤
ω
λ
1
|
s
∂
z
u
|
L
2
(
Ω
̃
y
)
2
(
z
)
+
ω
λ
1
|
s
u
|
L
2
(
Ω
̃
y
)
2
(
z
)
+
λ
2
2
|
s
u
|
L
2
(
Ω
̃
y
)
2
(
z
)
+
1
2
λ
2
|
1
s
f
|
L
2
(
Ω
̃
y
)
2
(
z
)
≤
ω
λ
1
η
+
ρ
λ
1
4
ω
(
λ
1
2
−
ρ
)
|
1
s
f
|
L
2
(
Ω
̃
y
)
2
(
z
)
,
where we have chosen
λ
2
such that
ω
λ
1
=
ρ
(
ω
λ
1
+
λ
2
2
)
in order to make optimal use of 𝜂.
To obtain an estimate for the two point boundary value problem (with respect to 𝑧) for the perturbed Schrödinger equation (2.2), we combine
(
2.7
)
+
C
0
C
1
ω
⋅
(
2.8
)
and use (2.3), (2.4) to obtain
1
C
1
∥
u
b
∥
L
2
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
2
+
b
∥
∂
z
u
b
∥
L
2
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
2
≤
∥
∇
y
u
b
∥
L
2
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
2
+
b
∥
∂
z
u
b
∥
L
2
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
2
+
C
0
C
1
|
s
u
b
|
L
2
(
Ω
̃
y
)
2
(
L
)
+
C
0
C
1
∥
σ
u
b
∥
L
2
(
0
,
L
;
L
2
(
∂
Ω
̃
y
)
)
2
+
C
0
C
1
(
b
σ
0
−
|
s
|
L
∞
)
|
u
b
(
0
)
|
L
2
(
Ω
̃
y
)
2
+
C
0
C
1
b
σ
L
|
u
b
(
L
)
|
L
2
(
Ω
̃
y
)
2
≤
2
ω
∫
0
L
∫
Ω
̃
y
s
|
ℑ
(
∂
z
u
b
u
̄
b
)
|
d
y
d
z
+
∫
0
L
∫
Ω
̃
y
|
ℜ
(
f
u
̄
b
)
|
d
y
d
z
−
C
0
C
1
∫
0
L
∫
Ω
̃
y
(
s
̃
−
∂
z
s
)
|
u
b
|
2
d
y
d
z
+
C
0
C
1
ω
∫
0
L
∫
Ω
̃
y
|
ℑ
(
f
u
̄
b
)
|
d
y
d
z
≤
(
μ
+
C
0
C
1
|
s
̃
−
∂
z
s
|
L
∞
)
∥
u
b
∥
L
2
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
2
+
1
2
μ
(
4
ω
2
|
s
|
L
∞
2
∥
∂
z
u
b
∥
L
2
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
2
+
(
1
+
(
C
0
C
1
ω
)
2
)
∥
f
∥
L
2
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
2
)
,
where
C
0
,
C
1
are as in (2.3).
In order to be able to achieve
C
1
μ
+
C
0
|
s
̃
−
∂
z
s
|
L
∞
<
1
and
2
ω
2
|
s
|
L
∞
2
μ
<
b
,
we invoke (2.4) and choose
μ
∈
(
2
ω
2
|
s
|
L
∞
2
b
,
1
−
C
0
|
s
̃
−
∂
z
s
|
L
∞
C
1
)
.
∎
Remark 2
Since we will apply the part of Proposition 1 concerning (2.1) with large
ω
=
ω
k
,
k
=
1
,
2
, we track frequency dependence of the estimates explicitly there.
On the other hand, the energy estimate in Proposition 1 for (2.2) will only be used for the (small) difference frequency
ω
=
ω
d
, and therefore the restrictions on the size of 𝜔 made there do not compromise practical applicability in our context.
A wave number explicit estimate for (2.2) follows by application of the transformation (1.6) to the Helmholtz equation without having to impose a smallness condition (2.4); see Proposition 2 in the next subsection.
Using Proposition 1, we can conclude the following results on the parameter-to-state map:
(2.10)
S
:
(
s
,
η
)
↦
(
ϕ
̂
1
,
ϕ
̂
2
,
ψ
̂
)
solving (
1.7
), (
1.8
)
and its linearization
d
S
(
s
,
η
)
(
δ
s
¯
,
δ
η
¯
)
:=
(
δ
ϕ
¯
1
,
δ
ϕ
¯
2
,
δ
ψ
¯
)
defined by the system
(2.11)
ı
ω
k
∂
z
s
δ
ϕ
¯
k
+
2
ı
ω
k
s
∂
z
δ
ϕ
¯
k
−
Δ
y
δ
ϕ
¯
k
=
−
ı
ω
k
∂
z
δ
s
¯
ϕ
̂
k
−
2
ı
ω
k
δ
s
¯
∂
z
ϕ
̂
k
in
Ω
̃
,
k
∈
{
1
,
2
}
,
ı
ω
d
∂
z
s
δ
ψ
¯
+
2
ı
ω
d
s
∂
z
δ
ψ
¯
−
ε
̃
∂
z
2
δ
ψ
¯
−
Δ
y
δ
ψ
¯
=
−
ı
ω
d
∂
z
δ
s
¯
ψ
̂
−
2
ı
ω
d
δ
s
¯
∂
z
ψ
̂
+
ı
ω
p
3
ε
̃
−
1
(
δ
η
¯
ϕ
̂
1
ϕ
̂
̄
2
+
η
(
δ
ϕ
¯
1
ϕ
̂
̄
2
+
ϕ
̂
1
δ
ϕ
¯
̄
2
)
)
in
Ω
̃
,
with boundary conditions
(2.12)
δ
ϕ
¯
k
(
0
,
y
)
=
0
,
y
∈
Ω
̃
y
,
∂
ν
y
δ
ϕ
¯
k
(
z
,
y
)
=
−
ı
ω
k
σ
k
δ
ϕ
¯
k
(
z
,
y
)
,
(
z
,
y
)
∈
(
0
,
L
)
×
∂
Ω
̃
y
,
k
∈
{
1
,
2
}
,
−
∂
z
δ
ψ
¯
(
0
,
y
)
=
−
ı
ω
d
σ
0
δ
ψ
¯
(
0
,
y
)
,
∂
z
δ
ψ
¯
(
L
,
y
)
=
−
ı
ω
d
σ
L
δ
ψ
¯
(
L
,
y
)
,
y
∈
Ω
̃
y
,
∂
ν
y
δ
ψ
¯
(
z
,
y
)
=
−
ı
ω
d
σ
δ
ψ
¯
(
z
,
y
)
,
(
z
,
y
)
∈
(
0
,
L
)
×
∂
Ω
̃
y
.
Corollary 1
Assume
h
,
Δ
y
h
∈
L
2
(
Ω
̃
y
)
.
Then 𝑆 is well-defined by (2.10) as an operator
S
:
D
s
1
×
L
2
(
0
,
L
;
L
p
(
Ω
̃
y
)
)
→
V
and Gâteaux differentiable for all
(
s
,
η
)
∈
D
s
2
×
L
2
(
0
,
L
;
L
p
(
Ω
̃
y
)
)
as an operator
W
1
,
∞
(
0
,
L
;
L
∞
(
Ω
̃
y
)
)
×
L
2
(
0
,
L
;
L
p
(
Ω
̃
y
)
)
→
V
,
where
D
s
1
=
{
s
∈
W
1
,
∞
(
0
,
L
;
L
∞
(
Ω
̃
y
)
)
:
(
2.4
) holds for
ω
=
ω
d
=
ω
1
−
ω
2
}
,
D
s
2
=
D
s
1
∩
W
2
,
q
(
0
,
L
;
L
∞
(
Ω
̃
y
)
)
,
V
=
(
L
∞
(
0
,
L
;
H
1
(
Ω
̃
y
)
)
∩
W
1
,
∞
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
)
2
×
H
1
(
(
0
,
L
)
×
Ω
̃
y
)
for some
p
,
q
∈
{
(
1
,
∞
]
if
Ω
̃
y
⊆
R
2
,
[
1
,
∞
]
if
Ω
̃
y
⊆
R
1
,
with
d
S
(
s
,
η
)
(
δ
s
¯
,
δ
η
¯
)
:=
(
δ
ϕ
¯
1
,
δ
ϕ
¯
2
,
δ
ψ
¯
)
defined by (2.11), (2.12).
Proof
The first part of Proposition 1 with
ω
=
ω
k
implies that the solutions are
ϕ
̂
1
,
ϕ
̂
2
∈
L
∞
(
0
,
L
;
H
1
(
Ω
̃
y
)
)
∩
W
1
,
∞
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
.
Setting
f
:=
ı
ω
p
3
ε
̃
−
1
η
ϕ
̂
1
ϕ
̂
̄
2
∈
L
2
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
and since
H
1
(
Ω
̃
y
)
continuously embeds into
L
2
p
/
(
p
−
1
)
, we conclude
ψ
̂
∈
H
1
(
(
0
,
L
)
×
Ω
̃
y
)
from the second part of Proposition 1 with
ω
=
ω
d
.
It is readily checked that the first variation of 𝑆 at
(
s
,
η
)
in direction
(
δ
s
¯
,
δ
η
¯
)
is given by the solution of (2.11), (2.12) and that
d
S
(
s
,
η
)
is a linear operator.
To prove boundedness of
d
S
(
s
,
η
)
:
W
1
,
∞
(
0
,
L
;
L
∞
(
Ω
̃
y
)
)
→
V
for
(
s
,
η
)
∈
D
s
2
×
L
2
(
0
,
L
;
L
p
(
Ω
̃
y
)
)
, we apply the first part of Proposition 1 with
ω
=
ω
k
and
f
:=
−
2
ı
ω
k
δ
s
¯
∂
z
ϕ
̂
k
,
h
=
0
.
To prove that
f
∈
L
∞
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
∩
W
1
,
1
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
, we need somewhat more regularity of
ϕ
̂
k
.
The latter can be derived by considering the PDE for
ϕ
̂
k
:=
∂
z
ϕ
̂
k
, (2.1) with
f
=
−
2
ı
ω
k
∂
z
s
∂
z
ϕ
̂
k
∈
L
∞
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
∩
W
1
,
1
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
for
s
∈
D
s
2
and implies
ϕ
̂
k
∈
W
1
,
∞
(
0
,
L
;
H
1
(
Ω
̃
y
)
)
∩
W
1
,
∞
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
.
Finally, application of the second part of Proposition 1 with
ω
=
ω
d
and
f
:=
−
2
ı
ω
d
δ
s
¯
∂
z
ψ
̂
+
ı
ω
p
3
ε
̃
−
1
(
δ
η
¯
ϕ
̂
1
ϕ
̂
̄
2
+
η
(
δ
ϕ
¯
1
ϕ
̂
̄
2
+
ϕ
̂
1
δ
ϕ
¯
̄
2
)
)
∈
L
2
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
with
f
(
0
)
∈
L
2
(
Ω
̃
y
)
implies the required bound on
δ
ψ
¯
in
H
1
(
(
0
,
L
)
×
Ω
̃
y
)
.
∎
2.2 Well-Posedness, Using the Helmholtz Equation for the Propagating Wave
Alternatively to the paraxial form of the PDE for
ψ
̂
, we consider
(2.13)
ı
ω
k
∂
z
s
ϕ
̂
k
+
2
ı
ω
k
s
∂
z
ϕ
̂
k
−
Δ
y
ϕ
̂
k
=
0
in
Ω
̃
,
k
∈
{
1
,
2
}
,
−
ω
d
2
s
2
ψ
̌
−
Δ
x
ψ
̌
=
ı
ω
p
3
1
Ω
(
η
̌
P
−
1
(
ϕ
̂
1
ϕ
̂
̄
2
)
)
in
Ω
,
for
ψ
̌
=
P
−
1
ψ
̂
,
η
̌
=
e
−
ı
ω
d
T
P
−
1
η
(cf. (1.5)), with boundary conditions
(2.14)
ϕ
̂
k
(
0
,
⋅
)
=
h
k
in
Ω
̃
y
,
∂
ν
y
ϕ
̂
k
(
z
,
y
)
=
−
ı
ω
k
σ
k
ϕ
̂
k
in
(
0
,
L
)
×
∂
Ω
̃
y
,
k
∈
{
1
,
2
}
,
∂
ν
ψ
̌
=
−
ı
ω
d
σ
ψ
̌
on
∂
Ω
,
with
P
(
Ω
)
⊇
Ω
̃
=
(
0
,
L
)
×
Ω
̃
y
, where 𝑃 is the paraxial transform defined by (1.4) and
1
Ω
the extension by zero operator for a function defined on
P
−
1
(
Ω
̃
)
to all of Ω.
This allows to use a possibly larger propagation domain for the 𝜓 wave as compared to the one for the
ϕ
1
,
ϕ
2
waves, and to get rid of the smallness assumption (2.4).
Well-posedness and Gâteaux differentiability of the forward operator defined by system (2.13), (2.14) follows analogously to Corollary 1 by combining the first part of Proposition 1 with known results on the Helmholtz equation with impedance boundary conditions; see, e.g., [19, 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 [19, Chapter 8] (where the constant coefficient case with
σ
=
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 𝑠 visible for our purposes.
Testing the general Helmholtz equation on a smooth domain Ω with impedance boundary conditions
(2.15)
−
ω
2
s
2
u
−
Δ
x
u
=
f
in
Ω
,
∂
ν
u
+
ı
ω
σ
u
=
h
on
∂
Ω
,
with 𝑢 and taking real and imaginary parts yields
(2.16)
∥
∇
x
u
∥
L
2
(
Ω
)
2
−
ω
2
∥
s
u
∥
L
2
(
Ω
)
2
=
∫
Ω
ℜ
(
f
u
)
d
x
,
ω
∫
∂
Ω
σ
|
u
|
2
d
Γ
=
∫
Ω
ℑ
(
f
u
)
d
x
.
On a starlike domain Ω in two space dimensions
d
=
2
with
(2.17)
x
⋅
ν
(
x
)
≥
c
Ω
>
0
on
∂
Ω
, we can also use
v
(
x
)
:=
x
⋅
∇
u
(
x
)
as a test function (provided it is contained in
H
1
(
Ω
)
) and consider the real part to arrive at
(2.18)
ω
2
∫
Ω
(
1
+
x
⋅
∇
s
s
)
|
s
u
|
2
d
x
−
1
2
ω
2
∫
∂
Ω
|
s
u
|
2
x
⋅
ν
(
x
)
d
Γ
(
x
)
+
1
2
∫
∂
Ω
|
∇
u
|
2
x
⋅
ν
(
x
)
d
Γ
(
x
)
+
ω
ℜ
(
ı
σ
∫
∂
Ω
u
x
⋅
∇
u
̄
d
Γ
(
x
)
)
=
ℜ
(
∫
Ω
f
x
⋅
∇
u
̄
d
x
)
.
Here we have used the identities
2
ℜ
(
s
2
u
x
⋅
∇
u
)
=
x
⋅
(
∇
|
s
u
|
2
−
2
s
∇
s
|
u
|
2
)
,
∇
(
x
⋅
∇
u
̄
)
=
∇
u
̄
+
(
x
⋅
∇
)
∇
u
̄
,
2
ℜ
(
∇
u
⋅
(
x
⋅
∇
)
∇
u
̄
)
=
∇
⋅
(
|
∇
u
|
2
x
)
−
d
|
∇
u
|
2
as well as the divergence theorem and
d
=
2
.
This yields the following results in a low and higher regularity regime of 𝑠 and correspondingly estimates on 𝑢 with different frequency dependence of the constants.
Proposition 2
If Ω is a bounded
C
1
,
1
smooth or convex
domain,
s
∈
L
∞
(
Ω
)
,
f
∈
H
1
(
Ω
)
*
,
σ
≠
0
a.e on
∂
Ω
, then a solution of (2.15) exists, is unique and satisfies the estimate
∥
∇
x
u
∥
L
2
(
Ω
)
2
+
∥
u
∥
L
2
(
Ω
)
2
≤
C
(
ω
)
(
∥
f
∥
H
1
(
Ω
)
*
2
+
∥
h
∥
H
−
1
/
2
(
∂
Ω
)
2
)
with some constant
C
(
ω
)
>
0
depending on 𝜔.
If additionally
f
∈
L
2
(
Ω
)
,
h
∈
H
1
/
2
(
∂
Ω
)
, then
(2.19)
∥
D
x
2
u
∥
L
2
(
Ω
)
2
+
∥
∇
x
u
∥
L
2
(
Ω
)
2
+
∥
u
∥
L
2
(
Ω
)
2
≤
C
(
ω
)
(
∥
f
∥
L
2
(
Ω
)
2
+
∥
h
∥
H
1
/
2
(
∂
Ω
)
2
)
.
Here
D
x
2
u
denotes the (weak derivative) Hessian of 𝑢.
If additionally Ω is a star-shaped domain with (2.17),
s
≥
0
,
ln
(
s
)
∈
W
1
,
∞
(
Ω
)
,
f
∈
L
2
(
Ω
)
,
h
=
0
,
d
=
2
and the relative smallness/largeness conditions
(2.20)
c
s
=
sup
x
∈
Ω
(
−
x
⋅
∇
s
(
x
)
s
(
x
)
)
<
1
,
σ
≥
c
σ
>
0
are satisfied, then the solution of (2.15) satisfies the estimate
(2.21)
1
ω
2
∥
D
x
2
u
∥
L
2
(
Ω
)
2
+
∥
∇
x
u
∥
L
2
(
Ω
)
2
+
ω
2
∥
u
∥
L
2
(
Ω
)
2
≤
C
(
1
+
1
ω
2
)
∥
f
∥
L
2
(
Ω
)
2
with some 𝐶 independent of 𝜔.
Proof
The well-posedness proof in the low regularity regime follows analogously to the one in the constant coefficient case [19, 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 (2.15) as
(
−
Δ
+
id
)
u
=
F
+
G
u
with the bounded and boundedly invertible operator
(
−
Δ
+
id
)
:
H
1
(
Ω
)
→
H
1
(
Ω
)
*
, the bounded operator
G
:
H
1
(
Ω
)
→
H
1
(
Ω
)
*
and the element
F
∈
H
1
(
Ω
)
*
defined by
⟨
(
−
Δ
+
id
)
u
,
v
⟩
(
H
1
)
*
,
H
1
:=
∫
Ω
(
∇
u
⋅
∇
v
̄
+
u
v
̄
)
d
x
,
⟨
G
u
,
v
⟩
(
H
1
)
*
,
H
1
:=
∫
Ω
(
1
+
ω
2
s
2
)
u
v
̄
d
x
−
ı
ω
∫
∂
Ω
σ
u
v
̄
d
Γ
,
⟨
F
,
v
⟩
(
H
1
)
*
,
H
1
:=
∫
Ω
f
v
̄
d
x
+
∫
∂
Ω
h
v
̄
d
Γ
for any
v
∈
H
1
(
Ω
)
.
This is equivalent to
(
I
−
K
)
u
=
b
with
K
=
(
−
Δ
+
id
)
−
1
G
,
b
=
(
−
Δ
+
id
)
−
1
F
,
where
b
∈
H
1
(
Ω
)
and the operator 𝐾 is bounded from
L
2
(
Ω
)
to
H
1
(
Ω
)
and thus compact when considered as an operator from
L
2
(
Ω
)
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
(
Ω
)
and finally lift regularity of 𝑢 to
H
1
(
Ω
)
by means of the fixed point identity
u
=
K
u
+
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
∈
L
2
(
Ω
)
, which due to
w
=
K
w
is automatically contained in
H
1
(
Ω
)
.
The second equation of (2.16) implies that 𝑤 satisfies homogeneous Dirichlet, hence by the impedance condition also homogeneous Neumann boundary conditions and can therefore be extended by zero to an
H
2
(
R
3
)
function
w
̂
that then satisfies the homogeneous Helmholtz equation
−
Δ
x
w
̂
+
s
̂
w
̂
on all of
R
3
with
s
̂
:=
1
Ω
(
s
−
1
)
+
1
and Sommerfeld radiation conditions.
From the results in [2, Section 8.3], we conclude that
w
̂
≡
0
and thus
w
≡
0
.
The higher order regularity result (2.19) follows from elliptic regularity and the fact that 𝑢 satisfies
−
Δ
x
u
=
f
̃
,
∂
ν
u
=
h
̃
with
f
̃
=
f
+
ω
2
s
2
u
∈
L
2
(
Ω
)
,
h
̃
=
h
−
ı
ω
σ
tr
∂
Ω
u
∈
H
1
/
2
(
∂
Ω
)
.
From (2.18) with (2.17), (2.20), applying the Cauchy–Schwarz inequality, we conclude
ω
2
(
1
−
c
Ω
)
∥
s
u
∥
L
2
(
Ω
)
2
+
1
2
γ
∥
∇
u
∥
L
2
(
∂
Ω
)
2
≤
1
2
ω
2
∫
∂
Ω
|
s
u
|
2
x
⋅
ν
(
x
)
d
Γ
(
x
)
+
ω
2
(
∫
∂
Ω
σ
2
|
x
|
2
|
u
|
2
d
Γ
(
x
)
+
∫
∂
Ω
|
x
|
2
|
f
|
2
d
Γ
(
x
)
)
∥
∇
u
∥
L
2
(
∂
Ω
)
.
Hence, by Young’s inequality and the second equation in (2.16), that is,
ω
∥
σ
u
∥
L
2
(
∂
Ω
)
2
≤
∥
f
/
s
∥
L
2
(
Ω
)
∥
s
u
∥
L
2
(
Ω
)
,
ω
2
(
1
−
c
Ω
)
∥
s
u
∥
L
2
(
Ω
)
2
≤
ω
2
2
∫
∂
Ω
|
s
u
|
2
x
⋅
ν
(
x
)
d
Γ
(
x
)
+
ω
2
γ
∫
∂
Ω
σ
2
|
x
|
2
|
u
|
2
d
Γ
(
x
)
+
1
γ
∫
∂
Ω
|
x
|
2
|
f
|
2
d
Γ
(
x
)
≤
ω
(
d
(
Ω
)
2
∥
s
2
σ
∥
L
∞
(
Ω
)
+
d
(
Ω
)
2
γ
∥
σ
∥
L
∞
(
Ω
)
)
∥
f
/
s
∥
L
2
(
Ω
)
∥
s
u
∥
L
2
(
Ω
)
+
d
(
Ω
)
2
γ
∥
f
∥
L
2
(
Ω
)
2
with
d
(
Ω
)
=
sup
x
∈
Ω
|
x
|
; hence, by another application of Young’s inequality (with factors
1
−
c
Ω
2
and
1
2
(
1
−
c
Ω
)
),
ω
2
∥
s
u
∥
L
2
(
Ω
)
2
≤
C
∥
f
∥
L
2
(
Ω
)
2
with
C
=
1
(
1
−
c
Ω
)
2
(
d
(
Ω
)
2
∥
s
2
σ
∥
L
∞
(
Ω
)
+
d
(
Ω
)
2
γ
∥
σ
∥
L
∞
(
Ω
)
)
2
∥
1
s
∥
L
∞
(
Ω
)
2
+
2
(
1
−
c
Ω
)
d
(
Ω
)
2
γ
.
Together with the first identity in (2.16), this yields
∥
∇
x
u
∥
L
2
(
Ω
)
2
+
ω
2
∥
s
u
∥
L
2
(
Ω
)
2
≤
(
2
C
+
1
ω
2
∥
1
s
∥
L
∞
(
Ω
)
2
)
∥
f
∥
L
2
(
Ω
)
2
,
and again, elliptic regularity yields (2.21).
∎
Analogously to Corollary 1, we obtain well-definedness and differentiability of the parameter-to-state map.
Corollary 2
Under the assumptions of Corollary 1, but without (2.4) and with 𝑉 replaced by
V
=
(
L
∞
(
0
,
L
;
H
1
(
Ω
̃
y
)
)
∩
W
1
,
∞
(
0
,
L
;
L
2
(
Ω
̃
y
)
)
)
2
×
H
2
(
Ω
)
,
the parameter-to-state map
(2.22)
S
:
(
s
,
η
)
↦
(
ϕ
̂
1
,
ϕ
̂
2
,
ψ
̂
)
solving (
2.13
), (
2.14
)
is well-defined as an operator
S
:
W
1
,
∞
(
0
,
L
;
L
∞
(
Ω
̃
y
)
)
×
L
2
(
Ω
)
→
V
and Gâteaux differentiable for all
(
s
,
η
)
∈
(
W
1
,
∞
(
0
,
L
;
L
∞
(
Ω
̃
y
)
)
∩
W
2
,
q
(
0
,
L
;
L
∞
(
Ω
̃
y
)
)
)
×
L
2
(
Ω
)
as an operator
W
1
,
∞
(
0
,
L
;
L
∞
(
Ω
̃
y
)
)
×
L
2
(
Ω
)
→
V
.
3 Convergence of a Frozen Newton Method for the Inverse Problem
Based on the parameter-to-state map 𝑆 defined in (1.13) and analyzed in Section 2 as well as the observation operator (1.14), we can write the inverse problem as
(3.1)
F
(
s
,
η
)
=
y
with
y
=
p
meas
and the forward operator
F
=
C
∘
S
and apply Newton’s method to (3.1), using the fact that
F
′
(
s
,
η
)
=
C
∘
S
′
(
s
,
η
)
and relying on our verification of Gâteaux differentiability from Section 2.
Alternatively, we here consider an all-at-once formulation
(3.2)
F
(
s
,
η
,
ϕ
̂
1
,
ϕ
̂
2
,
ψ
)
=
y
with
y
=
(
0
,
0
,
0
,
0
,
0
,
p
meas
)
based on the variational form of the system of initial-boundary value/boundary value problems (1.7), (1.8),
(3.3)
⟨
F
k
(
q
⃗
)
,
w
⃗
⟩
X
:=
∫
0
L
{
∫
Ω
̃
y
(
ı
ω
k
∂
z
s
k
ϕ
̂
k
w
̄
k
+
2
ı
ω
k
s
k
∂
z
ϕ
̂
k
w
̄
k
+
∇
y
ϕ
̂
k
⋅
∇
y
w
̄
k
)
d
y
+
ı
ω
k
∫
∂
Ω
̃
y
σ
k
ϕ
̂
k
w
̄
k
d
Γ
(
y
)
}
d
z
,
⟨
F
3
(
q
⃗
)
,
w
⃗
⟩
X
:=
∫
0
L
∫
Ω
̃
y
{
ı
ω
d
∂
z
s
1
ψ
̂
v
̄
+
2
ı
ω
d
s
1
∂
z
ψ
̂
v
̄
+
∇
y
ψ
̂
⋅
∇
y
v
̄
+
ε
̃
∂
z
ψ
̂
∂
z
v
̄
−
(
ı
ω
p
3
ε
̃
−
1
η
ϕ
̂
1
ϕ
̂
̄
2
)
v
̄
}
d
y
d
z
+
ı
ω
d
∫
0
L
∫
∂
Ω
̃
y
σ
ψ
̂
v
̄
d
Γ
(
y
)
d
z
+
ı
ω
d
ε
̃
(
σ
0
∫
Ω
̃
y
ψ
̂
(
0
)
v
̄
(
0
)
d
y
+
σ
L
∫
Ω
̃
y
ψ
̂
(
L
)
v
̄
(
L
)
d
y
)
,
⟨
F
k
+
3
(
q
⃗
)
,
w
⃗
⟩
X
:=
∫
Ω
̃
y
(
ϕ
̂
k
(
0
)
−
h
k
)
u
̄
k
d
y
,
⟨
F
6
(
q
⃗
)
,
w
⃗
⟩
X
:=
∫
Γ
ı
ω
d
ψ
̂
u
̄
d
Γ
for
k
∈
{
1
,
2
}
, where we have used the abbreviation
q
⃗
=
(
s
1
,
s
2
,
η
,
ϕ
̂
1
,
ϕ
̂
2
,
ψ
̂
)
and where
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 𝑠 in order to be able to verify the range invariance condition from [10] that allows to prove convergence of a frozen Newton method.
The additional constraint
0
=
P
(
q
⃗
)
:=
s
1
−
s
2
that will be imposed via a penalty term during the Newton type iteration (3.7) takes care of merging these two 𝑠 versions into one as the iteration proceeds.
The formal linearization of
F
=
(
F
1
,
…
,
F
6
)
at some
q
⃗
0
is given by
⟨
F
k
′
(
q
⃗
)
δ
q
⃗
¯
,
w
⃗
⟩
X
=
∫
0
L
{
∫
Ω
̃
y
(
ı
ω
k
(
∂
z
δ
s
¯
k
ϕ
̂
k
+
∂
z
s
k
δ
ϕ
¯
k
+
2
δ
s
¯
k
∂
z
ϕ
̂
k
+
2
s
k
∂
z
δ
ϕ
¯
k
)
w
̄
k
+
∇
y
δ
ϕ
¯
k
⋅
∇
y
w
̄
k
)
d
y
+
ı
ω
k
∫
∂
Ω
̃
y
σ
k
δ
ϕ
¯
k
w
̄
k
d
Γ
(
y
)
}
d
z
,
⟨
F
3
′
(
q
⃗
)
δ
q
⃗
¯
,
w
⃗
⟩
X
=
∫
0
L
∫
Ω
̃
y
{
ı
ω
d
(
∂
z
δ
s
¯
1
ψ
̂
+
∂
z
s
1
δ
ψ
¯
+
2
δ
s
¯
1
∂
z
ψ
̂
+
2
s
1
∂
z
δ
ψ
¯
)
v
̄
+
∇
y
δ
ψ
¯
⋅
∇
y
v
̄
+
ε
̃
∂
z
δ
ψ
¯
∂
z
v
̄
−
ı
ω
p
3
ε
̃
−
1
(
δ
η
¯
ϕ
̂
1
ϕ
̂
̄
2
+
η
δ
ϕ
¯
1
ϕ
̂
̄
2
+
η
ϕ
̂
1
δ
ϕ
¯
̄
2
)
v
̄
}
d
y
d
z
+
ı
ω
d
∫
0
L
∫
∂
Ω
̃
y
σ
δ
ψ
¯
v
̄
d
Γ
(
y
)
d
z
+
ı
ω
d
ε
̃
(
σ
0
∫
Ω
̃
y
δ
ψ
¯
(
0
)
v
̄
(
0
)
d
y
+
σ
L
∫
Ω
̃
y
δ
ψ
¯
(
L
)
v
̄
(
L
)
d
y
)
,
⟨
F
k
+
3
′
(
q
⃗
)
δ
q
⃗
¯
,
w
⃗
⟩
X
=
∫
Ω
̃
y
δ
ϕ
¯
k
(
0
)
u
̄
k
d
y
,
⟨
F
6
′
(
q
⃗
)
δ
q
⃗
¯
,
w
⃗
⟩
X
=
ı
ω
d
∫
Γ
δ
ψ
¯
u
̄
d
Γ
.
It is a bounded linear operator, and thus it is the Gâteaux derivative of 𝔽 when considered, e.g., as a mapping
F
:
X
→
Y
with
X
:=
W
1
,
p
(
0
,
L
;
L
p
(
Ω
̃
y
)
)
×
(
L
p
(
0
,
L
;
L
p
(
Ω
̃
y
)
)
)
2
×
(
W
1
,
∞
(
0
,
L
;
L
∞
(
Ω
̃
y
)
)
∩
L
2
(
0
,
L
;
H
1
(
Ω
̃
y
)
)
)
3
,
Y
:=
(
L
2
(
0
,
L
;
H
1
(
Ω
̃
y
)
*
)
)
2
×
(
H
1
(
(
0
,
L
)
×
Ω
̃
)
)
*
×
(
L
2
(
Ω
̃
)
)
2
×
L
2
(
Γ
)
for some
p
∈
[
1
,
∞
]
.
Note that this slightly differs from the function space setting suggested for the reduced setting by Corollary 1, but these differences are essential for the verification of convergence conditions for the frozen Newton method; see (3.5) 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 (3.2) as compared to a reduced one (3.1).
We can achieve the range invariance relation
(3.4)
F
(
q
⃗
)
−
F
(
q
⃗
0
)
=
F
′
(
q
⃗
0
)
r
(
q
⃗
)
by setting
r
(
q
⃗
)
:=
δ
q
⃗
¯
=
(
δ
s
¯
1
,
δ
s
¯
2
,
δ
η
¯
,
δ
ϕ
¯
1
,
δ
ϕ
¯
2
,
δ
ψ
¯
)
with
δ
s
¯
k
=
ϕ
̂
k
ϕ
̂
k
0
⋅
(
s
k
−
s
k
0
)
+
1
(
ϕ
̂
k
0
)
2
∫
0
z
(
s
k
−
s
k
0
)
(
z
′
)
(
ϕ
̂
k
0
∂
z
ϕ
̂
k
−
ϕ
̂
k
∂
z
ϕ
̂
k
0
)
(
z
′
)
d
z
′
,
δ
η
¯
=
ϕ
̂
1
ϕ
̂
2
̄
ϕ
̂
1
0
ϕ
̂
2
0
̄
⋅
(
η
−
η
0
)
+
(
ϕ
̂
1
−
ϕ
̂
1
0
)
(
ϕ
̂
2
̄
−
ϕ
̂
2
0
̄
)
ϕ
̂
1
0
ϕ
̂
2
0
̄
⋅
η
0
−
ω
d
ε
̃
ω
p
3
ϕ
̂
1
0
ϕ
̂
2
0
̄
(
ψ
̂
∂
z
(
s
1
−
s
1
0
)
+
2
(
s
1
−
s
1
0
)
∂
z
ψ
̂
−
ψ
̂
0
∂
z
δ
s
¯
1
−
2
δ
s
¯
1
∂
z
ψ
̂
0
)
,
δ
ϕ
¯
k
=
ϕ
̂
k
−
ϕ
̂
k
0
,
δ
ψ
¯
=
ψ
̂
−
ψ
̂
0
.
To this end, we have to choose
ϕ
̂
1
0
,
ϕ
̂
2
0
,
ψ
̂
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
δ
s
¯
1
0
,
δ
s
¯
2
0
,
δ
η
¯
0
.
Relying on this, we can also establish closeness of 𝑟 to the identity in the sense that
r
(
q
⃗
)
−
(
q
⃗
−
q
⃗
0
)
=
(
δ
s
¯
1
−
(
s
1
−
s
1
0
)
,
δ
s
¯
2
−
(
s
2
−
s
2
0
)
,
δ
η
¯
−
(
η
−
η
0
)
,
0
,
0
,
0
)
,
where
∥
δ
s
¯
k
−
(
s
k
−
s
k
0
)
∥
L
p
(
L
p
)
≤
∥
1
ϕ
̂
k
0
∥
L
∞
(
L
∞
)
∥
ϕ
̂
k
−
ϕ
̂
k
0
∥
L
∞
(
L
∞
)
∥
s
k
−
s
k
0
∥
L
p
(
L
p
)
+
∥
1
ϕ
̂
k
0
∥
L
∞
(
L
∞
)
2
L
1
/
p
∥
s
k
−
s
k
0
∥
L
p
(
L
p
)
∥
ϕ
̂
k
0
∂
z
(
ϕ
̂
k
−
ϕ
̂
k
0
)
−
(
ϕ
̂
k
−
ϕ
̂
k
0
)
∂
z
ϕ
̂
k
0
∥
L
p
*
(
L
∞
)
and similarly for
∥
∂
z
(
δ
s
¯
1
−
(
s
1
−
s
1
0
)
)
∥
L
p
(
L
p
)
, using the identity
∂
z
(
δ
s
¯
1
−
(
s
1
−
s
1
0
)
)
=
2
(
s
1
−
s
1
0
)
∂
z
(
ϕ
̂
k
ϕ
̂
k
0
)
+
ϕ
̂
k
−
ϕ
̂
k
0
ϕ
̂
k
0
∂
z
(
s
1
−
s
1
0
)
+
∂
z
(
1
(
ϕ
̂
k
0
)
2
)
∫
0
z
(
s
k
−
s
k
0
)
(
z
′
)
(
ϕ
̂
k
0
∂
z
ϕ
̂
k
−
ϕ
̂
k
∂
z
ϕ
̂
k
0
)
(
z
′
)
d
z
′
,
as well as
∥
δ
η
¯
−
(
η
−
η
0
)
∥
L
p
(
L
p
)
≤
∥
1
ϕ
̂
1
0
ϕ
̂
2
0
̄
∥
L
∞
(
L
∞
)
⋅
{
∥
(
ϕ
̂
1
−
ϕ
̂
1
0
)
ϕ
̂
2
0
̄
+
ϕ
̂
1
(
ϕ
̂
2
̄
−
ϕ
̂
2
0
̄
)
∥
L
∞
(
L
∞
)
∥
η
−
η
0
∥
L
p
(
L
p
)
+
ω
d
ε
̃
ω
p
3
(
∥
(
ψ
̂
−
ψ
̂
0
)
∂
z
(
s
1
−
s
1
0
)
∥
L
p
(
L
p
)
+
2
∥
(
s
1
−
s
1
0
)
(
∂
z
ψ
̂
−
∂
z
ψ
̂
0
)
∥
L
p
(
L
p
)
+
∥
ψ
̂
0
∂
z
(
δ
s
¯
1
−
(
s
1
−
s
1
0
)
)
∥
L
p
(
L
p
)
+
2
∥
(
δ
s
¯
1
−
(
s
1
−
s
1
0
)
)
∂
z
ψ
̂
0
∥
L
p
(
L
p
)
)
}
;
hence
(3.5)
∥
r
(
q
⃗
)
−
(
q
⃗
−
q
⃗
0
)
∥
X
≤
C
∥
q
⃗
−
q
⃗
0
∥
X
2
.
Based on the range invariance condition (3.4), we can rewrite the inverse problem of reconstructing
(
s
,
η
)
as a combination of an ill-posed linear and a well-posed nonlinear problem
(3.6)
F
′
(
q
⃗
0
)
r
̂
=
y
−
F
(
q
⃗
0
)
,
r
(
q
⃗
)
=
r
̂
,
P
q
⃗
=
0
for the unknowns
(
r
̂
,
q
⃗
)
.
Here
q
⃗
0
is fixed and
U
⊆
X
is a neighborhood of
q
⃗
0
.
For its regularized iterative solution, we consider the frozen Newton method
(3.7)
q
⃗
n
+
1
∈
arg
min
q
⃗
∈
U
∥
F
′
(
q
⃗
0
)
(
q
⃗
−
q
⃗
n
)
+
F
(
q
⃗
n
)
−
y
δ
∥
Y
2
+
α
n
∥
q
⃗
−
q
⃗
0
∥
X
2
+
∥
P
q
⃗
∥
Z
2
,
where
y
δ
≈
y
is the noisy data,
α
n
→
0
as
n
→
∞
and
Z
=
(
L
p
(
0
,
L
;
L
p
(
Ω
̃
y
)
)
)
.
To work in Hilbert spaces, we set
p
=
2
.
In order to prove convergence of (3.7), we additionally need a condition on compatibility of the linear operators
F
′
(
q
⃗
0
)
and 𝒫; see [10, Lemma B.1].
Sufficient for this condition is linearized uniqueness with
s
1
=
s
2
, that is, triviality of the intersection of the nullspaces
(3.8)
N
(
F
′
(
q
⃗
0
)
)
∩
N
(
P
)
=
{
0
}
.
For the model (1.7), (1.8) described by (3.3), this would likely require more than one excitation as well as observations at several frequencies, along with the corresponding extension of dependency of 𝑠 in order to allow for range invariance (3.4).
We here rigorously establish (3.8) for the alternative model (2.13), (2.14) using the Helmholtz equation for the outgoing wave
ψ
̂
and making the simplifying but still very realistic assumption that propagation of the excitation waves described by
ϕ
̂
1
,
ϕ
̂
2
takes place in a homogeneous domain with known and constant speed of sound.
The unknown parameters are then the squared slowness
s
̌
and the nonlinearity coefficient
η
̌
and the model reads as
−
ω
d
2
M
[
s
̌
ψ
̌
]
+
A
ψ
̌
+
ı
ω
d
D
ψ
̌
=
ı
ω
p
3
η
̌
f
,
where 𝑓 is given by the excitation waves
ϕ
̂
1
,
ϕ
̂
2
that can be precomputed under the assumption of
s
=
s
0
=
1
c
being known in the first equations of (2.13).
Moreover, on the boundary
∂
Ω
, we assume to know
s
̌
=
s
0
2
=
1
c
2
, and define the operators 𝒜, 𝒟, ℳ by
A
u
=
(
v
↦
∫
Ω
∇
u
⋅
∇
v
d
x
+
∫
∂
Ω
β
u
v
d
s
)
,
D
u
=
(
v
↦
b
∫
Ω
∇
u
⋅
∇
v
d
x
+
∫
∂
Ω
(
σ
+
b
β
)
u
v
d
s
)
,
M
u
=
(
v
↦
∫
Ω
u
v
d
x
+
∫
∂
Ω
σ
b
s
0
2
u
v
d
s
)
.
In our linearized uniqueness proof, we will assume that 𝒜, 𝒟 and ℳ are simultaneously diagonalizable.
This holds true, e.g., in the case
σ
=
0
, where
M
=
id
and
D
=
b
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 𝐼 of difference frequencies
ω
d
and with at least two different sets of excitation frequencies; more precisely,
ω
2
=
ω
1
ω
̃
ω
1
+
ω
̃
ω
̃
,
ω
̃
∈
{
ω
̃
1
,
ω
̃
2
}
⊆
R
+
,
ω
d
=
ω
1
−
ω
2
∈
I
⊆
R
+
with 𝐼 countable containing at least one accumulation point, which implies
ω
p
3
=
ω
d
2
ω
̃
.
Correspondingly, we consider
ψ
̌
=
ψ
̌
(
ω
d
,
ω
̃
)
as a function of
ω
d
, while the original coefficients
s
̌
,
η
̌
of course stay independent of
(
ω
d
,
ω
̃
)
.
In order to satisfy range invariance, we will introduce an artificial dependence
s
̌
=
s
̌
(
ω
d
,
ω
̃
)
, while keeping
η
̌
independent of
ω
d
and
ω
̃
.
Therewith, the inverse problem reads as
F
(
s
̌
,
η
̌
,
ψ
̌
)
=
y
with
F
1
(
q
⃗
)
=
−
ω
d
2
(
M
(
s
̌
ψ
̌
)
+
ı
ω
̃
η
̌
f
)
+
A
ψ
̌
+
ı
ω
d
D
ψ
̌
,
F
2
(
q
⃗
)
=
ı
ω
d
tr
Γ
̌
ψ
̌
,
where
q
⃗
=
(
s
̌
,
η
̌
,
ψ
̌
)
,
and
y
=
(
0
,
p
meas
)
.
It is readily checked that range invariance (3.4) holds with
r
(
s
̌
,
η
̌
,
ψ
̌
)
=
(
(
s
̌
−
s
̌
0
)
ψ
̌
ψ
̌
0
,
η
̌
−
η
̌
0
,
ψ
̌
−
ψ
̌
0
)
and
(3.9)
∥
r
(
q
⃗
)
−
(
q
⃗
−
q
⃗
0
)
∥
X
≤
∥
1
ψ
̌
0
∥
L
∞
(
Ω
)
⋅
∥
(
s
̌
−
s
̌
0
)
(
ψ
̌
−
ψ
̌
0
)
∥
L
2
(
Ω
)
≤
C
∥
q
⃗
−
q
⃗
0
∥
X
2
≤
c
∥
q
⃗
−
q
⃗
0
∥
X
with
X
=
L
2
(
Ω
)
×
L
2
(
Ω
)
×
H
2
(
Ω
)
, provided
∥
1
/
ψ
̌
0
∥
L
∞
(
Ω
)
<
∞
, that is,
ψ
̌
0
is bounded away from zero.
Here 𝑐 can be made small by restricting
q
⃗
to a sufficiently small neighborhood of
q
⃗
0
.
Defining 𝒫, 𝕏, 𝕐, ℤ by
P
(
q
⃗
)
(
ω
d
)
:=
(
s
̌
(
ω
d
,
ω
̃
1
)
−
s
̌
(
ω
d
,
ω
̃
2
)
,
s
̌
(
ω
d
,
ω
̃
1
)
−
1
μ
(
I
)
∫
I
s
̌
(
⋅
,
ω
̃
1
)
d
μ
)
,
X
=
L
μ
2
(
I
;
L
2
(
Ω
)
)
2
×
L
2
(
Ω
)
×
H
2
(
Ω
)
,
Y
=
L
μ
2
(
I
;
L
2
(
Γ
)
)
2
,
Z
=
L
μ
2
(
I
;
L
2
(
Ω
)
)
2
for some finite measure 𝜇 on 𝐼, we can write the inverse problem as (3.6) and use (3.7) for its regularized numerical solution.
To prove convergence of (3.7) we require linearized uniqueness (3.8), which we verify as follows.
For
(
δ
s
̌
¯
,
δ
η
̌
¯
,
δ
ψ
̌
¯
)
∈
N
(
P
)
, that is,
s
̌
independent of
ω
d
and
ω
̃
, the assumption
F
′
(
s
̌
0
,
η
̌
0
,
ψ
̌
0
)
(
δ
s
̌
¯
,
δ
η
̌
¯
,
δ
ψ
̌
¯
)
=
0
with
s
̌
0
:
≡
1
c
2
and
η
̌
0
arbitrary reads as
(
−
ω
d
2
1
c
2
M
+
A
+
ı
ω
d
D
)
δ
ψ
̌
¯
(
ω
d
,
ω
̃
)
=
ω
d
2
(
M
[
δ
s
̌
¯
ψ
̌
0
(
ω
d
,
ω
̃
)
]
+
ı
ω
̃
δ
η
̌
¯
f
(
ω
d
,
ω
̃
)
)
,
ı
ω
d
tr
Γ
̌
δ
ψ
̌
¯
(
ω
d
,
ω
̃
)
=
0
for all
ω
d
∈
I
,
ω
̃
∈
{
ω
̃
1
,
ω
̃
2
}
.
We now assume that
ψ
̌
0
and 𝑓 have been chosen independent of
ω
d
and
ω
̃
,
(3.10)
ψ
̌
0
(
ω
d
,
ω
̃
)
≡
ψ
̌
0
,
f
(
ω
d
,
ω
̃
)
≡
f
with
ψ
̌
≠
0
,
f
≠
0
a.e.
,
and diagonalize the operators 𝒜, 𝒟, ℳ, by means of their eigensystems
(
λ
j
,
φ
j
k
)
j
∈
N
,
k
∈
K
j
,
(
ρ
j
,
φ
j
k
)
j
∈
N
,
k
∈
K
j
,
(
μ
j
,
φ
j
k
)
j
∈
N
,
k
∈
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
,
σ
,
β
,
b
≥
0
, we have
μ
j
≥
1
,
λ
j
,
ρ
j
≥
0
,
λ
j
→
∞
as
j
→
∞
.
We will additionally assume
(3.11)
(
ρ
j
μ
j
=
ρ
ℓ
μ
ℓ
and
λ
j
μ
j
2
=
λ
ℓ
μ
ℓ
2
)
⟹
j
=
ℓ
,
which is the case, e.g., if
σ
=
0
since then
μ
j
≡
1
.
This yields
(3.12)
0
=
δ
ψ
̌
¯
(
x
0
;
ω
d
,
ω
̃
)
=
c
2
ω
d
2
∑
j
=
1
∞
∑
k
∈
K
j
1
−
ω
d
2
μ
j
+
c
2
λ
j
+
ı
ω
d
c
2
ρ
j
(
a
j
k
+
ı
ω
̃
b
j
k
)
φ
j
k
(
x
0
)
for all
x
0
∈
Γ
,
ω
d
∈
I
,
ω
̃
∈
{
ω
̃
1
,
ω
̃
2
}
,
for
a
j
k
=
μ
j
⟨
δ
s
̌
¯
ψ
̌
0
,
φ
j
k
⟩
L
2
(
Ω
)
,
b
j
k
=
⟨
δ
η
̌
¯
f
,
φ
j
k
⟩
L
2
(
Ω
)
.
It is straightforward to show that the functions
ω
d
↦
1
w
j
(
ı
ω
d
)
with
w
j
(
ı
ω
d
)
:=
−
ω
d
2
μ
j
+
c
2
λ
j
+
ı
ω
d
c
2
ρ
j
are linearly independent on 𝐼.
Indeed, assuming
∑
j
=
1
∞
1
w
j
(
ı
ω
d
)
c
j
=
0
for all
ω
d
∈
I
and a sequence
c
⃗
=
(
c
j
)
j
∈
N
, after multiplication with
∏
ℓ
∈
N
w
ℓ
(
ı
ω
d
)
, we conclude
0
=
W
c
⃗
(
z
)
:=
∑
j
=
1
∞
∏
ℓ
≠
j
w
ℓ
(
z
)
c
j
for all
z
∈
ı
I
,
and thus, since the function
W
c
⃗
is analytic, it vanishes on all of ℂ.
Inserting the roots
z
±
k
=
−
c
ρ
j
2
±
c
2
ρ
j
2
4
−
λ
j
of
w
k
and using the fact that, under condition (3.11), they are distinct (in the sense that
z
+
j
=
z
+
ℓ
and
z
−
j
=
z
−
ℓ
implies
j
=
ℓ
), we obtain
c
⃗
=
0
.
By the linear independence of
(
1
w
j
)
j
∈
N
, we conclude from (3.12) that
∑
k
∈
K
j
(
a
j
k
+
ı
ω
̃
b
j
k
)
φ
j
k
(
x
0
)
=
0
for all
x
0
∈
Γ
,
j
∈
N
,
ω
̃
∈
{
ω
̃
1
,
ω
̃
2
}
.
Under a linear independence assumption on the individual eigenspaces (cf., e.g., [12]),
for all
j
∈
N
,
(
∑
k
∈
K
j
γ
k
φ
j
k
(
x
)
=
0
for all
x
∈
Γ
)
⟹
(
γ
k
=
0
for all
k
∈
K
j
)
,
we conclude
{
a
j
k
+
ı
ω
̃
1
b
j
k
=
0
,
a
j
k
+
ı
ω
̃
2
b
j
k
=
0
,
for all
j
∈
N
,
k
∈
K
j
,
which implies
a
j
k
=
0
,
b
j
k
=
0
for all
j
∈
N
,
k
∈
K
j
; thus, under assumption (3.10), we get
δ
s
̌
¯
=
0
,
δ
η
̌
¯
=
0
.
According to [10, Theorem 2.6], we obtain the following.
Theorem 1
Let
q
⃗
†
=
(
s
̌
,
η
̌
,
ψ
̌
)
be a solution to (3.6) and let for the noise level
δ
≥
∥
y
δ
−
y
∥
Y
the stopping index
n
*
=
n
*
(
δ
)
be chosen such that
n
*
(
δ
)
→
0
,
δ
∑
j
=
0
n
*
(
δ
)
−
1
c
j
α
n
*
(
δ
)
−
j
−
1
−
1
/
2
→
0
as
δ
→
0
,
with 𝑐 as in (3.9).
Moreover, let (3.10) with
1
/
ψ
̌
0
∈
L
∞
(
Ω
)
hold and let the operators 𝒜, 𝒟, ℳ be jointly diagonalizable with eigenvalues satisfying (3.11), and
N
(
F
′
(
q
⃗
0
)
)
+
⊆
N
(
P
)
.
Then there exists
ϱ
>
0
sufficiently small such that, for
q
⃗
0
∈
B
ϱ
(
q
⃗
†
)
, the iterates
(
q
⃗
n
δ
)
n
∈
{
1
,
…
,
n
*
(
δ
)
}
are well-defined by (3.7), remain in
B
ϱ
(
q
⃗
†
)
and converge in 𝕏,
∥
q
⃗
n
*
(
δ
)
δ
−
q
⃗
†
∥
X
→
0
as the noise level
δ
→
0
.
In the noise free case
δ
=
0
,
n
*
(
δ
)
=
∞
, we have
∥
q
⃗
n
−
q
⃗
†
∥
X
→
0
as
n
→
∞
.
