Barbara Kaltenbacher$^*$
Department of Mathematics
Alpen-Adria-Universität Klagenfurt

9020 Klagenfurt, Austria

William Rundell
Department of Mathematics
Texas A&M University

Texas 77843, USA

(Communicated by Habib Ammari)

We consider an undetermined coefficient inverse problem for a non-
linear partial
differential equation occurring in high intensity ultrasound propagation
as used in acoustic tomography.
In particular, we investigate the recovery of the nonlinearity coefficient
commonly labeled as $B/A$ in the literature which is part of a space
dependent coefficient $\kappa$ in the  Westervelt equation governing nonlinear
Corresponding to the typical measurement setup, the overposed data consists
of time trace measurements on some zero or one dimensional set $\Sigma$
representing the receiving transducer array.
After an analysis of the map from $\kappa$ to the overposed data, we show
injectivity of its linearisation and use this as motivation for several
iterative schemes to recover $\kappa$.
Numerical simulations will also be shown to illustrate the efficiency of the

§ INTRODUCTION

The use of ultrasound is well established in the imaging of human tissue.
High intensity ultrasound is modeled by nonlinear wave equations, in which
a certain ratio of Taylor expansion coefficients $B/A$ governs the nonlinearity.
Recently, it has been shown that this $B/A$ parameter is sensitive to
differences in the tissue properties,
thus appropriate for characterization of biological tissues, see, e.g.,
[3, 7, 8, 22, 48, 52, 53].
Therefore, when viewed as a spatially varying coefficient, it can be used for
medical imaging purposes, known as  acoustic nonlinearity parameter tomography.
This parameter appears in the pdes describing high intensity ultrasound
propagation, thus the related imaging problem becomes a coefficient
identification for these pdes.

We will therefore give a brief introduction into the relevant
models and highlight the one  of our focus:
namely the Westervelt equation.
Then we will specify the possible physical measurements leading
to overposed data and state the resulting inverse problem.

§.§.§ The model

For a brief derivation of the fundamental acoustic equations, we refer to, e.g., the review [27].
More details can be found, for example, in [11, 15].

The main physical quantities involved in the description of
sound propagation are

* the acoustic particle velocity $\vec{v}$;

* the acoustic pressure $p$;

* the mass density $\varrho$;

that can be decomposed into their constant mean and a fluctuating part

\[
\vec{v}=\vec{v}_0+\vec{v}_\sim\,, \quad p=p_0+p_\sim\,, \quad \varrho=\varrho_0+\varrho_\sim,
\]

where $\vec{v}_0=0$ in the absence of a flow.

These quantities are interrelated by balance equations of momentum and mass, as well as the state equation relating the acoustic pressure and density fluctuations $p_\sim$ and $\varrho_\sim$. The latter contains a so-called parameter of nonlinearity $B/A$.
Combining these balance laws, inserting the state equation and and dropping certain terms according to a certain hierarchy, which in nonlinear acoustics is known as Blackstock's scheme [34, 4] one arrives at second order wave equations in terms of the space- and time dependent fluctuating quantities.

This first of all yields Kuznetsov's equation [33, 32]

\begin{equation}\label{Kuznetsov}
{p_\sim}_{tt}- c^2\triangle p_\sim - b \triangle {p_\sim}_t =
\left(\frac{1}{\varrho_0 c^2}\frac{B}{2A} p_\sim^2 + \varrho_0 |\vec{v}|^2\right)_{tt}
\end{equation}

where $b$ is the diffusivity of sound,
and we have related the velocity to the pressure via the linearised force balance

\begin{equation}\label{eq:forcebal}
\varrho_0\vec{v}_t=-\nabla p_\sim\,.
\end{equation}

rep C, 2.

If we ignore local nonlinear effects modeled by the quadratic velocity term, thus approximating
$\rho_0|\vec{v}|^2_{tt}\approx \frac{1}{\rho_0^2c^2}p^2_{tt}$,
we arrive at the Westervelt equation

\begin{equation}\label{Westervelt}
{p_\sim}_{tt}- c^2\triangle p_\sim - b \triangle {p_\sim}_t =
\frac{\beta_a}{\varrho_0 c^2} {p_\sim^2}_{tt}
\end{equation}

with $\beta_a = 1 + B/(2A)$, cf.,
Under the already made assumption $\nabla\times\vec{v}=0$ on a simply connected domain there exists an acoustic velocity potential $\psi$ with $\vec{v}=-\nabla\psi$, whose constant part by (<ref>) can be chosen such that

\begin{equation}\label{prhopsitp}
\varrho_0\psi_t= p\,.
\end{equation}

Hence both equations (<ref>) and (<ref>) can also
be written in terms of the  acoustic velocity potential $\psi$

\begin{equation}\label{WestKuz_psi}
\psi_{tt}-c^2\triangle\psi-b\triangle\psi_t=\frac{1}{c^2}\Bigl(\beta_a(\psi_t)^2
+s_{\mbox{\footnotesize WK}} \left[c^2|\nabla\psi|^2-(\psi_t)^2\right]\Bigr)_t
\end{equation}

with $s_{\mbox{\footnotesize WK}}=0$ for (<ref>) and $s_{\mbox{\footnotesize WK}}=1$ for (<ref>).

Further simplifications of the model lead to the Khokhlov-Zabolotskaya-Kuznet-sov (KZK) equation [51]
and the well-known Burgers' equation in one space dimension [6, 14, 42].

rep C, 2.

We mention in passing that there exist also more complex and higher order models. Since it takes into account most of the relevant physical effects, the Westervelt equation appears to be the best established model in the physics and engineering literature of nonlinear acoustics. We will therefore also adhere to this model here.

§.§.§ The inverse problem

The imaging task under consideration consists of identifying $\kappa=\kappa(x)$ in the Westervelt equation in pressure formulation (skipping the subscript $\sim$)

\begin{equation}\label{eqn:Westervelt_p}
p_{tt}-c^2\triangle p - b\triangle p_t = \kappa(x)(p^2)_{tt} + r(x,t)
\end{equation}

($\kappa(x)= \frac{1}{\rho_0c^2} \beta_a(x)$) or in velocity potential formulation

\begin{equation}\label{eqn:Westervelt_psi}
\psi_{tt}-c^2\triangle \psi - b\triangle \psi_t = \kappa(x)(\psi_t^2)_t + r(x,t)
\end{equation}

($\kappa(x)= \frac{1}{c^2} \beta_a(x)$)
where $r$ is a given excitation function and pressure and velocity potential are related by

The system will typically be at rest initially, leading to homogeneous initial conditions on $p$, $p_t$ or $\psi$, $\psi_t$, respectively.
The spatial domain $\Omega\subseteq\mathbb{R}^d$, $d\in\{1,2,3\}$ on which the pdes are supposed to hold will be assumed to be smooth and bounded and the Laplacian equipped with boundary conditions on $\partial\Omega$. For simplicity one might think of homogeneous Dirichlet ones here but note that also other boundary conditions – in-homogeneous Neumann for modeling excitation by a transducer array, absorbing or impedance conditions for modeling damping or reflections at the boundary – are relevant in this context. Excitation will here be modeled by an interior space and time dependent source term $r$, see also [28].

Since the actuating piezoelectric transducers are typically arranged in an array, that is, a surface $\Gamma$ lying in the interior of the computational domain, modeling excitation by means of interior sources appears to be justified as follows.
We can consider $r$ as an approximation of a source $q\cdot\delta_\Gamma$  concentrated on $\Gamma$, with $q=(c^2\tilde{r}+b\tilde{r}_t)$, (that is, $\tilde{r}(t)= \frac{1}{b}\int_0^t\exp(-\frac{c^2}{b}(t-s))\tilde{r}(s)\, ds$, in view of the fact that formally

\begin{eqnarray}
&&\left\{\begin{array}{r} p_{tt}-c^2\triangle p - b\triangle p_t - \kappa(x)(p^2)_{tt} =0\mbox{ in }\Omega\\ p=0\mbox{ on }\partial\Omega\\
{}[[\partial_\nu p]]=\tilde{r}\mbox{ on }\Gamma\end{array}\right.
\label{eqn:interiorsource}\\
&&\Leftrightarrow \
\int_\Omega
\Bigl((p_{tt}- \kappa(x)(p^2)_{tt})v + (c^2\nabla p+b\nabla p_t)\cdot\nabla v\Bigr)\, dx = \int_{\Gamma} (c^2\tilde{r}+b\tilde{r}_t) v\, ds
\nonumber\\
&& \hspace*{7cm}\mbox{ for all } v\in H_0^1(\Omega)
\nonumber\\
&&\Leftrightarrow \
\left\{\begin{array}{r} p_{tt}-c^2\triangle p - b\triangle p_t - \kappa(x)(p^2)_{tt} =(c^2\tilde{r}+b\tilde{r}_t) \delta_{\Gamma}\mbox{ in }\Omega\\ p =0\mbox{ on }\partial\Omega\end{array}\right.
\nonumber
\end{eqnarray}

where $[[\partial_\nu p]]$ denotes the jump of the normal derivative of $p$ over the interface $\Gamma$.

With $\tilde{r}/\rho_0$ being the acceleration of the transducer in normal direction, via the force balance (<ref>), the jump condition in (<ref>) after time integration corresponds to a matching of velocities at the solid-fluid interface.

rep C, 1.

As a counterpart to this imposed excitation, measurements of the acoustic pressure at an array of transducers or hydrophones are made.
Thus the overposed data consists of time trace measurements at some point $x_0$ or on some surface $\Sigma$ contained in $\overline{\Omega}$, see also Figure <ref> below
rep M, 1.

\begin{equation}\label{eqn:obs}
g(t)=p(x_0,t) \quad \mbox{ or } \quad g(x,t)=p(x,t) \,,\ x\in\Sigma\,, t\in(0,T).
\end{equation}

Thus the pressure formulation (<ref>) appears to be the more direct one in the sense that the observations are just point evaluations of the state, whereas (<ref>) would require time integration of the data in order to relate state and observations via (<ref>). Moreover the quadratic nonlinearity comes in terms of a higher derivative in (<ref>) than in (<ref>), which makes analysis and numerics slightly more involved.
Henceforth we will focus on the pressure formulation of the Westervelt equation (<ref>).

Thus the inverse problem under consideration here is the following.

Identify $\kappa=\kappa(x)$ in

\[
\begin{array}{rcll} p_{tt}-c^2\triangle p - b\triangle p_t &=& \kappa(x)(p^2)_{tt} + r&\mbox{ in }\Omega\times(0,T)\\
p &=&0&\mbox{ on }\partial\Omega\times(0,T)\\
\end{array}
\]

from observations (<ref>). Here, $c^2$, $b$, and $r=r(x,t)$ are known constants and functions, respectively.

Since this paper is a first step into the mathematics of this inverse problem
we provide some background information.
Thus in Section sec:analysis_forward, we give some regularity results
on the forward operator.
Section sec:injectivity deals with uniqueness.
Note that much of the results from recovery of space dependent coefficients
from time trace data in
hyperbolic or parabolic pdes, see, e.g., [24] and the
references therein, do not apply here.
These are the main differences:
first, the wave equations appearing as models are strongly damped and therefore
behave to some extent like a parabolic pde, see Remark rem:parabolic below;
second, the reformulation as parabolic pde contains a nonlocal in time term due to the integration operator;
third is the fact that the base pde is nonlinear and furthermore the
coefficient  $\kappa$ to be determined is intrinsically coupled to this
Thus, as we will see, the statement
“but linear inverse equations behave even better” certainly holds true.
We also show injectivity of the linearised inverse problem at $\kappa=0$
in case the excitation is chosen as an appropriate combination
of products of a space and time dependent functions.
rep C, 4.
In this section we will also comment on ill-posedness of the problem,
which is expected to be exponential.
In Section sec:regmeth we discuss some of the classical regularisation
paradigms in the context of the nonlinearity imaging task.
Section sec:reconstructions is devoted to numerical experiments
demonstrating the above iterative reconstruction methods.

We mention that an inverse source problem related to the linearisation of this
inverse problem with a higher order model of nonlinear acoustics has recently
been considered in [50].

Moreover, we point to [23], where in case of different observations, namely data $p(x, t_0)$ with some $t_0 > 0$, global Lipschitz stability is shown even in the more complicated setting of the Lamé operator replacing the Laplacian for the corresponding linearised inverse problem (but probably extendable to the nonlinear setting).
rep M, 2.

§ ANALYSIS OF THE FORWARD PROBLEM

We consider the operator $F:\mathcal{D}(F)\to Y$, $F=\mbox{tr}_\Sigma\circ G$,
$G(\kappa)=p$, where $p$ solves (<ref>) with homogeneous
and linear boundary conditions
– which we will assume to hold for all equations appearing in this section
without explicitly mentioning it, thus we will write $\mathcal{A}$ for
$-\triangle$ (or more generally a second order elliptic differential operator)
equipped with these boundary conditions on a sufficiently smooth boundary
$\partial\Omega$ so that for $s\in[0,2]$,
$\dot{H}^s(\Omega):=\mathcal{D}(\mathcal{A}^{s/2})\subseteq H^s(\Omega)$
with continuous embedding and equivalent norms.
This includes, e.g., the Laplacian with Dirichlet
($\dot{H}^1(\Omega)=H_0^1(\Omega)$) or impedance
($\dot{H}^1(\Omega)=H^1(\Omega)$) boundary conditions,
but excludes pure Neumann conditions.

§.§ Well-definedness and Fréchet differentiability

For taking differences between $\tilde{p}=G(\tilde{\kappa})$ and $p=G(\kappa)$ the following identity on the quadratic nonlinearity will be useful.

\[
\begin{aligned}
\tilde{\kappa}(\tilde{p}^2)_{tt}-\kappa(p^2)_{tt}
=& 2\tilde{\kappa}\tilde{p}\tilde{p}_{tt}-2\kappa p p_{tt} + 2 \tilde{\kappa}\tilde{p}_t^2-2\kappa p_t^2 \\
=& 2\underline{d\kappa}\tilde{p}\tilde{p}_{tt} + 2\kappa v\tilde{p}_{tt} + 2\kappa p v_{tt}
+ 2\underline{d\kappa}\,\tilde{p}_t^2 + 2 \kappa (\tilde{p}_t+ p_t) v_t
\end{aligned}
\]

where $\underline{d\kappa}=\tilde{\kappa}-\kappa$, $v = \tilde{p}-p$.
This implies that $v=G(\tilde{\kappa})-G(\kappa)$ solves

\begin{equation}\label{eqn:v}
(1-2\kappa p)v_{tt}+c^2\mathcal{A} v +b\mathcal{A} v_t - 2 \kappa (\tilde{p}_t+ p_t)\, v_t - 2\kappa \tilde{p}_{tt} \, v
= 2\underline{d\kappa}(\tilde{p}\tilde{p}_{tt} + \tilde{p}_t^2)\,,
\end{equation}

and the (so far formal) linearisation $z=G'(\kappa)\underline{d\kappa}$ solves

\begin{equation}\label{eqn:z}
(1-2\kappa p)z_{tt}+c^2\mathcal{A} z +b\mathcal{A} z_t - 4 \kappa p_t\, z_t - 2\kappa p_{tt} \, z
= 2\underline{d\kappa}(p p_{tt} + p_t^2)\,,
\end{equation}

so that the first order Taylor remainder $w=G(\tilde{\kappa})-G(\kappa)-G'(\kappa)(\tilde{\kappa}-\kappa)$ satisfies

\begin{equation}\label{eqn:w}
\begin{aligned}
(1-2\kappa p)w_{tt}&+c^2\mathcal{A} w +b\mathcal{A} w_t - 4 \kappa p_t\, w_t - 2\kappa p_{tt} \, w\\
&= 2\underline{d\kappa}(v \tilde{p}_{tt} + p v_{tt} + (\tilde{p}_t+p_t) v_t ) + 2\kappa (v v_{tt} + v_t^2) \,,
\end{aligned}
\end{equation}

in all three cases with homogeneous initial and boundary conditions.
Therefore the following lemma will be useful for estimating these differences
and establishing Fréchet differentiability of the forward map.
Note that the existing energy estimates for the linearised version of the Westervelt
equation provide higher regularity and exponential decay, but their transfer to
space dependent parameters would also involve derivatives of the coefficients and
in particular of $\kappa$. However, we wish to allow $\kappa$ to be discontinuous
(e.g., piecewise constant of piecewise continuous) as relevant for the underlying
imaging task in order to reproduce sharp interfaces between tissue with different
properties. Thus the use of some $L^p$ space for $\kappa$ is essential.
On the other hand, as the estimates below show, actually $L^\infty$ appears to be
the minimal requirement on $\kappa$.

In order to handle both well-definedness and differentiability of $G$, we first of all establish some energy estimates for a linear version of the Westervelt equation.

In the estimates below we will make use of continuity of the embeddings of $\Htwodot\to L^\infty(\Omega)$ and $\Honedot\to L^6(\Omega)$, more precisely

\begin{equation}\label{eqn:C1C2}
\begin{aligned}
&\|\phi\|_{L^6(\Omega)}\leq C_1 \|\mathcal{A}^{1/2}\phi\|_{L^2(\Omega)} \,, \quad \phi \in \mathcal{D}(\mathcal{A}^{1/2})=\Honedot  \\
&\|\phi\|_{L^\infty(\Omega)}\leq C_2 \|\mathcal{A}\phi\|_{L^2(\Omega)} \,, \quad \phi \in \mathcal{D}(\mathcal{A})=\Htwodot  \,.
\end{aligned}
\end{equation}

Moreover, we will impose smallness of a certain combination of the coefficients $\alpha,\gamma$

\begin{equation}\label{eqn:smallness}
\bar{\gamma}:= \|\tfrac{d}{dt}\ln(\tfrac{\gamma}{\alpha})\|_{L^1(0,T;L^\infty(\Omega))} <\tfrac{1}{16}
\end{equation}

Note that since these may contain common physical coefficients,  the quotient may enable some cancellations.

rep C, 3.

For $\alpha,\beta,\gamma\in L^\infty(0,T;L^\infty(\Omega))$, $\alpha,\beta,\gamma>0$, $\delta\in L^\infty(0,T;L^3(\Omega))$,  $\mu\in L^\infty(0,T;L^2(\Omega))$,
$\ln(\tfrac{\gamma}{\alpha}) \in W^{1,1}(0,T;L^\infty(\Omega))$
with (<ref>)
any solution $u$ to

\begin{equation}\label{eqn:linWest}
\alpha u_{tt}+\beta\mathcal{A} u_t + \gamma \mathcal{A} u + \delta u_t + \mu u = f
\end{equation}

satisfies the estimates

\begin{equation}
\|\mathcal{A}^{1/2}u_t\|_{L^\infty(0,t;L^2(\Omega))}^2
+ \|\sqrt{\tfrac{\beta}{\alpha}}\mathcal{A}u_t\|_{L^2(0,t;L^2(\Omega))}^2
\end{equation}

\begin{equation*}
\begin{aligned}
&\leq 8\Bigl(
\|\mathcal{A}^{1/2}u_t(0)\|_{L^2(\Omega)}^2 + \|\sqrt{\tfrac{\gamma(0)}{\alpha(0)}}\mathcal{A}u(0)\|_{L^2(\Omega)}^2
+\|\tfrac{1}{\alpha\beta}\|_{L^\infty(0,t;L^\infty(\Omega))} \|f\|_{L^2(0,t;L^2(\Omega))}^2\\
&+ C_1^2\|\tfrac{\delta}{\sqrt{\alpha\beta}}\|_{L^\infty(0,t;L^{3}(\Omega))}^2  \|\mathcal{A}^{1/2}u_t\|_{L^2(0,t;L^2(\Omega))}^2\\
\end{aligned}
\end{equation*}

\begin{equation}
\|u_{tt}\|_{L^2(0,t;L^2(\Omega))} \
\leq
\|\tfrac{1}{\alpha}\|_{L^\infty(0,t;L^\infty(\Omega))} \|f\|_{L^2(0,t;L^2(\Omega))}\label{eqn:enest_utt} \qquad \qquad \qquad \qquad
\end{equation}

\begin{equation*}
\begin{aligned}
&+ \|\tfrac{\beta}{\alpha}\|_{L^\infty(0,t;L^\infty(\Omega))}^{1/2}
\|\sqrt{\tfrac{\beta}{\alpha}}\mathcal{A}u_t\|_{L^2(0,t;L^2(\Omega))}\\
&+ \|\tfrac{\gamma}{\alpha}\|_{L^\infty(0,t;L^\infty(\Omega))}^{1/2} \sqrt{T}
\|\sqrt{\tfrac{\gamma}{\alpha}}\mathcal{A}u\|_{L^\infty(0,t;L^2(\Omega))}\\
&+ C_1 \|\tfrac{\delta}{\alpha}\|_{L^2(0,t;L^3(\Omega))} \|\mathcal{A}^{1/2}u_t\|_{L^\infty(0,t;L^2(\Omega))} + C_2 \|\tfrac{\mu}{\alpha}\|_{L^\infty(0,t;L^2(\Omega))}^2  \|\mathcal{A}u\|_{L^2(0,t;L^2(\Omega))}\,,
\end{aligned}
\end{equation*}

with $C_1$, $C_2$ as in (<ref>).

Dividing (<ref>) by $\alpha$, multiplying with $\mathcal{A} u_t$ and integrating over $(0,t)\times \Omega$, using the identity

\[
\tfrac{\gamma}{\alpha} \mathcal{A}u \mathcal{A}u_t = \tfrac12 \tfrac{d}{dt} \Bigl(\tfrac{\gamma}{\alpha} (\mathcal{A}u)^2\Bigr) - \tfrac12\Bigl[\tfrac{\gamma}{\alpha}\Bigr]_t (\mathcal{A}u)^2
\]

as well as Young's inequality we obtain

\[
\begin{aligned}
+ \int_0^t\|\sqrt{\tfrac{\beta}{\alpha}}\mathcal{A}u_t(\tau)\|_{L^2(\Omega)}^2\, d\tau
+ \tfrac12\|\sqrt{\tfrac{\gamma}{\alpha}}\mathcal{A}u(t)\|_{L^2(\Omega)}^2
\\
&=\tfrac12\|\mathcal{A}^{1/2}u_t(0)\|_{L^2(\Omega)}^2 + \tfrac12\|\sqrt{\tfrac{\gamma(0)}{\alpha(0)}}\mathcal{A}u(0)\|_{L^2(\Omega)}^2\\
&\qquad+ \int_0^t \Bigl( (\tfrac{1}{\alpha}f(\tau) - \tfrac{\delta}{\alpha} u_t(\tau) - \tfrac{\mu}{\alpha} u(\tau) ) \mathcal{A}u_t(\tau)
+  \tfrac12\Bigl[\tfrac{\gamma}{\alpha}\Bigr]_t (\mathcal{A}u(\tau))^2\Bigr)\, d\tau\\
\tfrac12\|\mathcal{A}^{1/2}u_t(0)\|_{L^2(\Omega)}^2 + \tfrac12\|\sqrt{\tfrac{\gamma(0)}{\alpha(0)}}\mathcal{A}u(0)\|_{L^2(\Omega)}^2\\
&\qquad+\tfrac12\int_0^t\|\sqrt{\tfrac{\alpha}{\beta}}(\tfrac{1}{\alpha}f - \tfrac{\delta}{\alpha} u_t - \tfrac{\mu}{\alpha} u )\|_{L^2(\Omega)}^2\, d\tau
+ \tfrac12 \int_0^t\|\sqrt{\tfrac{\beta}{\alpha}}\mathcal{A}u_t(\tau)\|_{L^2(\Omega)}^2\, d\tau\\
&\qquad + \tfrac12\int_0^t\|\tfrac{\alpha}{\gamma}\, \Bigl[\tfrac{\gamma}{\alpha}\Bigr]_t \|_{L^\infty(\Omega)}\, d\tau
\sup_{\tau\in(0,t)}\|\sqrt{\tfrac{\gamma}{\alpha}}(\mathcal{A}u(\tau))\|_{L^2(\Omega)}^2
\end{aligned}
\]

for any $t\in[0,T]$.
After multiplication by two and moving some terms we get, with $\tfrac{\alpha}{\gamma}\, \Bigl[\tfrac{\gamma}{\alpha}\Bigr]_t = \tfrac{\gamma_t}{\gamma}-\tfrac{\alpha_t}{\alpha}$$=\tfrac{d}{dt}\ln(\tfrac{\gamma}{\alpha})$,

\[
\begin{aligned}
+ \|\sqrt{\tfrac{\beta}{\alpha}}\mathcal{A}u_t\|_{L^2(0,t;L^2(\Omega))}^2
+ \|\sqrt{\tfrac{\gamma}{\alpha}}\mathcal{A}u(t)\|_{L^2(\Omega)}^2 \\
&\qquad\qquad- 2\|\revision{\tfrac{d}{dt}\ln(\tfrac{\gamma}{\alpha})}\|_{L^1(0,T;L^\infty(\Omega))} \|\sqrt{\tfrac{\gamma}{\alpha}}\mathcal{A}u\|_{L^\infty(0,t;L^2(\Omega))}^2
\\
\|\mathcal{A}^{1/2}u_t(0)\|_{L^2(\Omega)}^2 + \|\sqrt{\tfrac{\gamma(0)}{\alpha(0)}}\mathcal{A}u(0)\|_{L^2(\Omega)}^2
+\|\tfrac{1}{\alpha\beta}\|_{L^\infty(0,t;L^\infty(\Omega))} \|f\|_{L^2(0,t;L^2(\Omega))}^2\\
&\qquad+ \|\tfrac{\delta}{\sqrt{\alpha\beta}}\|_{L^\infty(0,t;L^3(\Omega))}^2  \|u_t\|_{L^2(0,t;L^6(\Omega))}^2
+ \|\tfrac{\mu}{\sqrt{\alpha\beta}}\|_{L^\infty(0,t;L^2(\Omega))}^2  \|u\|_{L^2(0,t;L^\infty(\Omega))}^2\,.
\end{aligned}
\]

Thus, for any $t'\in[0,T]$ taking the supremum over $t\in[0,t']$, using the fact that
$\sup_{t\in[0,t']}\sum_{i=1}^m a_i(t)\geq 2^{-m}\sum_{i=1}^m\sup_{t_i\in[0,t']} a_i(t_i)$, and assuming $\|\revision{\tfrac{d}{dt}\ln(\tfrac{\gamma}{\alpha})}\|_{L^1(0,T;L^\infty(\Omega))}$ to be smaller than $\frac{1}{16}$, we obtain (<ref>).

To obtain equation (<ref>), we just insert
(<ref>) after applying the $L^2(0,T;L^2(\Omega))$ norm to both sides
of the pde identity $u_{tt}= \tfrac{1}{\alpha} f
-\frac{\beta}{\alpha}\mathcal{A} u_t - \frac{\gamma}{\alpha} \mathcal{A} u - \frac{\delta}{\alpha} u_t - \frac{\mu}{\alpha} u$.

We first of all apply this lemma together with a fixed point argument to conclude well-posedness of the nonlinear equation (<ref>), i.e., of

\[
(1-2\kappa p)p_{tt}+c^2\mathcal{A} p +b\mathcal{A} p_t = r+2\kappa p_t^2
\]

by means of Banach's Fixed Point Theorem applied to the operator $\mathcal{T}$
mapping $p$ to a solution $\mathcal{T}(p)=p^+$ of the linear problem

\begin{equation}\label{eqn:calT}
(1-2\kappa p)p^+_{tt}+c^2\mathcal{A} p^+ +b\mathcal{A} p^+_t = r+2\kappa p_t^2\,,
\end{equation}

The difference $\hat{p}^+=\tilde{p}^+-p^+$ between values $\tilde{p}^+=\mathcal{T}(\tilde{p})$ and $p^+=\mathcal{T}(p)$ can be characterized by the pde

\begin{equation}\label{eqn:calT_diff}
(1-2\kappa p)\hat{p}^+_{tt}+c^2\mathcal{A} \hat{p}^+ +b\mathcal{A} \hat{p}^+_t
= 2\kappa (\tilde{p}^+_{tt}\, \hat{p} + (p_t+\tilde{p}_t)\, \hat{p}_t)\,,
\end{equation}

where $\hat{p}=\tilde{p}-p$.
Therefore we will make use of Lemma lem:enest with $\alpha=1-2\kappa p$, $\beta\equiv b$, $\gamma\equiv c^2$, $\delta\equiv\mu\equiv0$, both for proving that $\mathcal{T}$ is a self-mapping and for establishing its contractivity.
For this purpose it will be convenient to note down the estimates (<ref>), (<ref>) in the following somewhat compressed form for this particular setting, assuming additionally that $-\tfrac12\leq 2\kappa p\leq \tfrac12$, which we will actually guarantee by a proper choice of the domain of $\mathcal{T}$, and which implies that $\tfrac12\leq\alpha\leq\tfrac32$.

\begin{eqnarray}
+ \tfrac23 b\|\mathcal{A}u_t\|_{L^2(0,t;L^2(\Omega))}^2
+ \tfrac{2(1-16\bar{\gamma})}{3} c^2\|\mathcal{A}u\|_{L^\infty(0,t;L^2(\Omega))}^2
\nonumber\\
&&\quad\leq 8\Bigl(\|\mathcal{A}^{1/2}u_t(0)\|_{L^2(\Omega)}^2
+\tfrac{2}{b} \|f\|_{L^2(0,t;L^2(\Omega))}^2\Bigr)\,,
\label{eqn:enest_calT}\\
\nonumber\\
&&\quad \leq 2 \|f\|_{L^2(0,t;L^2(\Omega))}
+ \sqrt{2b}\|\sqrt{\tfrac{\beta}{\alpha}}\mathcal{A}u_t\|_{L^2(0,t;L^2(\Omega))}^2
\nonumber\\
\nonumber\\
&&\quad\leq 2\Bigl( 2 +\sqrt{\tfrac{c^2 T}{(1-16\bar{\gamma})b}}\Bigr)
\|f\|_{L^2(0,t;L^2(\Omega))}\,.
\label{eqn:enest_utt_calT}
\end{eqnarray}

We aim at keeping track of the constants $b$, $c$, $T$, $|\Omega|$ in our
estimates (at least those for establishing $\mathcal{T}$ as a self-mapping),
since in reality they can be of very different orders of magnitude.

In view of the estimates (<ref>), (<ref>), we will work on the function space

\[
V = H^2(0,T;L^2(\Omega))\cap W^{1,\infty}(0,T;\Honedot  )\cap H^1(0,T;\Htwodot  )
\]

and consider the fixed point operator defined in (<ref>) $\mathcal{T}:M\to M$ (with the self-mapping property yet to be established) on a bounded subset of $V$

\begin{equation}\label{eqn:M}
\begin{aligned}
M:=\{&u\in V\,:\, \|u\|_{L^\infty(0,T;L^\infty(\Omega))}\leq R_0,\\
+ \tfrac23 b\|\mathcal{A}u_t\|_{L^2(0,T;L^2(\Omega))}^2\\
&\qquad\qquad\qquad+ \tfrac{2(1-16\bar{\bar{\gamma}})}{3} c^2\|\mathcal{A}u\|_{L^\infty(0,T;L^2(\Omega))}^2
\leq\frac{R_1^2}{b},\\
&\|u_{tt}\|_{L^2(0,t;L^2(\Omega))}\leq R_2\}
\end{aligned}
\end{equation}

\begin{equation}\label{eqn:R0}
R_0\leq\frac{1}{4\|\kappa\|_{L^\infty}}\,, \quad
\bar{\bar{\gamma}}<\frac{1}{16}\,,
\end{equation}

and further conditions will be imposed on $R_1,R_2$, cf. (<ref>), (<ref>), (<ref>), (<ref>) below.

Now fix $p\in M$. We are going to prove that then $p^+ = \mathcal{T}p\in M$, that is, $\mathcal{T}$ is a self-mapping on $M$, by applying the estimates (<ref>), (<ref>) to (<ref>).
To this end, observe that (<ref>) immediately implies

\[
\tfrac12\leq\alpha\leq\tfrac32 \mbox{ for }\alpha=1-2\kappa p\,.
\]

\[
\bar{\gamma} = \|\tfrac{\alpha_t}{\alpha}\|_{L^1(0,T;L^\infty(\Omega))}
= \|\tfrac{2\kappa p_t}{\alpha}\|_{L^1(0,T;L^\infty(\Omega))}
\leq 4\|\kappa\|_{L^\infty(\Omega)}\|p_t\|_{L^1(0,T;L^\infty(\Omega))}\,,
\]

where by interpolation and with the constant $C_{(1+s)/2}$ of the embedding
$\mathcal{D}(\mathcal{A}^{(1+s)/2})$ $(=\dot{H}^{1+s}(\Omega))\to L^\infty(\Omega)$
for $s>\frac12$ (note that we are trying to be minimal with respect to negative powers of the typically small constant $b$ here)

\begin{equation}\label{eqn:L1}
\begin{aligned}
\leq C_{(1+s)/2} \sqrt{T}\|\mathcal{A}^{s/2}\mathcal{A}^{1/2}p_t\|_{L^2(0,T;L^2(\Omega))}\\
&\qquad\leq C_{(1+s)/2} \sqrt{T}
\|\mathcal{A}p_t\|_{L^2(0,T;L^2(\Omega))}^s
\|\mathcal{A}^{1/2}p_t\|_{L^2(0,T;L^2(\Omega))}^{1-s}\\
&\qquad\leq C_{(1+s)/2} \sqrt{T}
\|(b\mathcal{A})p_t\|_{L^2(0,T;L^2(\Omega))}^s
\|(b\mathcal{A})^{1/2}p_t\|_{L^2(0,T;L^2(\Omega))}^{1-s}\, b^{-(1+s)/2}\\
&\qquad\leq C_{(1+s)/2} \sqrt{T}
(\sqrt{T}R_1)^{1-s}\, b^{-(1+s)/2}
\end{aligned}
\end{equation}

(where we have used the first two terms in the $R_1$ estimate of (<ref>), after multiplying this estimate with $b$)
so that we get $\bar{\gamma}\leq\bar{\bar{\gamma}}$, provided

\begin{equation}\label{eqn:R1}
4\|\kappa\|_{L^\infty(\Omega)} C_{(1+s)/2}
(\tfrac32)^{s/2} T^{1-s/2}\, R_1\, b^{-(1+s)/2}
\leq\bar{\bar{\gamma}}\,.
\end{equation}

Thus, the assumptions of Lemma lem:enest are satisfied and we can make use of estimates (<ref>), (<ref>) with $f=r+2\kappa p_t^2$ cf. (<ref>).
For this purpose we estimate

\[
\|f\|_{L^2(0,T;L^2(\Omega))} \leq
\|r\|_{L^2(0,T;L^2(\Omega))} + 2\|\kappa\|_{L^\infty(\Omega)} \|p_t\|_{L^4(0,T;L^4(\Omega))}^2
\]

where by Hölder's inequality with exponent $\frac32$

\begin{equation}\label{eqn:L4}
\begin{aligned}
\|p_t\|_{L^4(0,T;L^4(\Omega))}^2 =& \Bigl(\int_0^T\int_\Omega p_t^4\, dx\, dt\Bigr)^{1/2}
\leq |\Omega|^{1/6} \Bigl(\int_0^T\|p_t\|_{L^6(\Omega)}^4\, dt\Bigr)^{1/2}\\
\leq& |\Omega|^{1/6} T^{1/2} C_1^2\|\mathcal{A}^{1/2}p_t\|_{L^\infty(0,T;L^2(\Omega))}^2
\leq |\Omega|^{1/6} T^{1/2} C_1^2 \frac{R_1^2}{b}
\end{aligned}
\end{equation}

Thus (<ref>), (<ref>) yield the $R_1$ and $R_2$ estimates on $p^+$ in (<ref>) provided that for

\begin{equation}\label{eqn:fbar}
\bar{f}:= \|r\|_{L^2(0,T;L^2(\Omega))} + 2\|\kappa\|_{L^\infty(\Omega)} |\Omega|^{1/6} T^{1/2} C_1^2 \frac{R_1^2}{b}
\end{equation}

we can guarantee

\begin{equation}\label{eqn:R2a}
\|\mathcal{A}^{1/2}u_1\|_{L^2(\Omega)}^2 + 2c^2\|\mathcal{A}u_0\|_{L^2(\Omega)}^2
+ \frac{2}{b}\bar{f}^2\Bigr) \leq \frac{R_1^2}{b}
\end{equation}

\begin{equation}\label{eqn:R2b}
2\Bigl( 2 +\sqrt{\tfrac{c^2 T}{(1-16\bar{\gamma})b}}\Bigr)\bar{f} \leq R_2\,.
\end{equation}

Finally, the $R_0$ estimate on $p^+$ in (<ref>) follows from the $R_1$ estimate on $p^+$ in (<ref>) (which we have just established) provided that

\begin{equation}\label{eqn:R3}
C_2 \Bigl({\frac{3}{2(1-16\bar{\bar{\gamma}})bc^2}} R_1\Bigr)^{1/2} \leq R_0\,,
\end{equation}

since $\|p^+\|_{L^\infty(0,T;L^\infty(\Omega))}\leq C_2 \|\mathcal{A}p^+\|_{L^\infty(0,T;L^2(\Omega))}$.
Thus we have shown that $\mathcal{T}$ is a self-mapping, provided (<ref>), (<ref>), (<ref>), (<ref>), (<ref>) hold. These can be achieved by making $\|\kappa\|_{L^\infty(\Omega)}$ small:
Choose $R_1 > 4 \|r\|_{L^2(0,T;L^2(\Omega))}$, set $R_0:=C_2 \sqrt{\frac{3}{2(1-16\bar{\bar{\gamma}})bc^2}}$
$R_2:=2\Bigl( 2 +\sqrt{\tfrac{c^2 T}{(1-16\bar{\gamma})b}}\Bigr)\bar{f}$, so that we have (<ref>), (<ref>);
then possibly decrease $\|\kappa\|_{L^\infty(\Omega)}$ to achieve (<ref>), (<ref>), (<ref>).
These requirements also show that smallness of $b$ can to some extent be compensated by smallness of $T$, cf. (<ref>), (<ref>), (<ref>).

We now proceed to derive contractivity of $\mathcal{T}$ by applying the
estimates (<ref>), (<ref>) to the pde
In the above, we have already shown that for any $p\in M$,
under conditions (<ref>), (<ref>), (<ref>),
(<ref>), (<ref>), the coefficients $\alpha=1-2\kappa p_t$,
$\beta\equiv b$, $\gamma\equiv c^2$, $\delta\equiv0$, $\mu\equiv0$ satisfy the
assumptions of Lemma lem:enest with $\frac12\leq\alpha\leq\frac32$,
$\bar{\gamma}\leq\bar{\bar{\gamma}}<\frac{1}{16}$, so that the estimates
(<ref>), (<ref>) with $u(0)=0$, $u_t(0)=0$
apply to (<ref>). Using them together with Gronwall's inequality
yields existence of  constants $C,\;\tilde{C}>0$
(depending only on the constants $R_0,R_1,R_2,\bar{\bar{\gamma}}$
in the definition of $M$ as well as $b,c^2$ and $T$), such that

\[
\begin{aligned}
\leq C \|2\kappa (\tilde{p}^+_{tt}\, \hat{p} + (p_t+\tilde{p}_t)\, \hat{p}_t)\|_{L^2(0,T;L^2(\Omega))}\\
&\qquad\leq 2C \|\kappa\|_{L^\infty(\Omega)}
\Bigl(\|\tilde{p}^+_{tt}\|_{L^2(0,T;L^2(\Omega))}\|\hat{p}\|_{L^\infty(0,T;L^\infty(\Omega))}
\\&\qquad\qquad\qquad\qquad+ \|p_t+\tilde{p}_t\|_{L^4(0,T;L^4(\Omega))} \|\hat{p}_t\|_{L^4(0,T;L^4(\Omega))}\Bigr)\\
&\qquad\leq 2C \|\kappa\|_{L^\infty(\Omega)}
\Bigl(R_2C_2\|\mathcal{A}\hat{p}\|_{L^\infty(0,T;L^2(\Omega))}
\\&\qquad\qquad\qquad\qquad+ 2|\Omega|^{1/6} T^{1/2} C_1^2 \frac{R_1}{\sqrt{b}} \|\mathcal{A}^{1/2}\hat{p}_t\|_{L^\infty(0,T;L^2(\Omega))}\Bigr)\\
&\qquad\leq \tilde{C} \|\kappa\|_{L^\infty(\Omega)} \|\hat{p}\|_V
\end{aligned}
\]

by (<ref>), (<ref>) and an estimate analogous to (<ref>).
For $\|\kappa\|_{L^\infty(\Omega)}$ small enough this implies contractivity.

For any $u_0\in \Htwodot $, $u_1\in \Honedot  $, $r\in L^2(0,T;L^2(\Omega))$, $T>0$ and any $C^{1,1}$ domain $\Omega$ there exists $R>0$ such that for all $\kappa\in L^\infty(\Omega)$, $\|\kappa\|_{L^\infty(\Omega)}<R$ the quasilinear pde (<ref>) with homogeneous Dirichlet boundary conditions and initial conditions $u(0)=u_0$, $u_t(0)=u_1$ is uniquely solvable on the space $M$ defined in (<ref>).

Hence, for any $C^{1,1}$ manifold $\Sigma\subseteq\overline{\Omega}$,
the forward operator $F:\mathcal{D}(F)\to Y$, $F(\kappa)=\mbox{tr}_\Sigma p$
is well-defined and bounded on $\mathcal{D}(F)=B_R(0)\subseteq L^\infty(\Omega)$
and maps into any space $Y$ in which the space
$W^{1,\infty}(0,T;H^{1/2}(\Sigma))\cap H^1(0,T;H^{3/2}(\Sigma))$ is continuously embedded.

Since measurements are done on values and not on derivatives of $p$, a natural image space $Y$ of $F$ is actually $L^Q(0,T;L^P(\Omega))$ for some $P,Q\in[1,\infty]$.

To prove Fréchet differentiability of $F$, we apply the estimates from Lemma lem:enest together with Gronwall's inequality to  (<ref>), (<ref>), (<ref>) with $\alpha=1-\kappa p$, $\beta\equiv b$, $\gamma\equiv c^2$,

\[
\delta = \begin{cases}
2\kappa(\tilde{p}_t+p_t) \mbox{ for } \eqref{eqn:v}\\
4\kappa p_t \mbox{ for } \eqref{eqn:z}, \eqref{eqn:w}
\end{cases}\,, \quad
\mu = \begin{cases}
2\kappa\tilde{p}_{tt} \mbox{ for } \eqref{eqn:v}\\
2\kappa p_{tt} \mbox{ for } \eqref{eqn:z}, \eqref{eqn:w}
\end{cases}\,,
\]

\[
2\underline{d\kappa}(\tilde{p}\tilde{p}_{tt} + \tilde{p}_t^2) \mbox{ for } \eqref{eqn:v}\\
2\underline{d\kappa}(p p_{tt} + p_t^2)  \mbox{ for } \eqref{eqn:z}\\
2\underline{d\kappa}(v \tilde{p}_{tt} + p v_{tt} + (\tilde{p}_t+p_t) v_t ) + 2\kappa (v v_{tt} + v_t^2)
\mbox{ for } \eqref{eqn:w}\end{cases}\,.
\]

This together with estimates similar to those in the proof of
Proposition prop:forward_welldef yields existence of a constant
$C>0$ (depending only on the constants $R_0,R_1,R_2,\bar{\bar{\gamma}}$
in the definition of $M$ as well as $b,\;c^2$, the time horizon $T$, and the
radius $R$ according to Proposition prop:forward_welldef), such that
for any $\kappa\in B_R(0)\subseteq L^\infty(\Omega)$

\[
\begin{aligned}
&\|v\|_V\leq C\|\underline{d\kappa}\|_{L^\infty(\Omega)}\,, \quad
\|z\|_V\leq C\|\underline{d\kappa}\|_{L^\infty(\Omega)}\,, \\
&\|w\|_V\leq C\Bigl(\|\underline{d\kappa}\|_{L^\infty(\Omega)}\|v\|_V+ R\|v\|_V^2)\leq C^2(1+CR)\|\underline{d\kappa}\|_{L^\infty(\Omega)}^2 \,.
\end{aligned}
\]

For any $u_0\in \Htwodot $, $u_1\in \Honedot  $, $r\in L^2(0,T;L^2(\Omega))$, $T>0$ any $C^{1,1}$ domain $\Omega$ and any $C^{1,1}$ manifold $\Sigma\subseteq\overline{\Omega}$, with $R$ as in Proposition prop:forward_welldef, the forward operator $F:\mathcal{D}(F)\to Y$, $F(\kappa)=\mbox{tr}_\Sigma p$ is Fréchet differentiable on $\mathcal{D}(F)=B_R(0)\subseteq L^\infty(\Omega)$ with derivative given by $F'(\kappa)\underline{d\kappa} = \mbox{tr}_\Sigma z$ for $z$ solving (<ref>) with homogeneous initial and boundary conditions.

Applying a similar reasoning to the difference $\bar{w}=\tilde{z}-z$ between $G'(\tilde{\kappa})=\tilde{z}$ and $G'(\kappa)=z$, which then with $\underline{d\kappa}=\tilde{\kappa}-\kappa$ and $v$ as in (<ref>) satisfies (similarly to $w$ in (<ref>))

\begin{equation}\label{eqn:wbar}
\begin{aligned}
(1-2\kappa p)\bar{w}_{tt}&+c^2\mathcal{A} \bar{w} +b\mathcal{A} \bar{w}_t - 4 \kappa p_t\, \bar{w}_t - 2\kappa p_{tt} \, \bar{w}\\
&= 2\underline{d\kappa}(v \tilde{p}_{tt} + p v_{tt} + (\tilde{p}_t+p_t) v_t ) + 2\kappa (\tilde{z} v_{tt} + 2v_t\tilde{z}_t)
\end{aligned}
\end{equation}

we obtain Lipschitz continuity of $F'$ on $B_R(0)\subseteq L^\infty(\Omega)$.

§.§ Adjoint of $F'(\kappa)$

The goal of this section is to compute the adjoint of $F'(\kappa)$, which is
needed, for example, for formulating Landweber iteration,
a Gauss-Newton method, or the necessary optimality conditions for Tikhonov
regularisation, as well as for assessing source conditions.
From the implementation point of view (see Section sec:regmeth),
we prefer a Hilbert space setting to a general Banach space one and therefore

\[
X=\dot{H}^s(\Omega)\,, \quad Y=L^2(0,T;L^2(\Sigma))\,,
\]

where $\dot{H}^s(\Omega)=\mathcal{D}(\mathcal{A}^{s/2})$ with $s\geq0$ appropriately chosen, see Section sec:regmeth.

The adjoint must satisfy the identity

\begin{equation}\label{eqn:adj_def}
\langle \underline{d\kappa}, F'(\kappa)^*y\rangle_X \stackrel{!}{=} \langle F'(\kappa)\underline{d\kappa}, y\rangle_Y =
\int_0^T\int_\Sigma z\, y \, dS \, dt
\end{equation}

for any $\underline{d\kappa}\in X$, $y\in Y$.
Taking into account (<ref>) and the fact that $(1-2\kappa p)z_{tt} - 4 \kappa p_t\, z_t - 2\kappa p_{tt} \, z = \Bigl((1-2\kappa p)z\Bigr)_{tt}$, this leads us to defining an adjoint state as the solution $a$ to

\begin{equation*}
a(T)=0, \ a'(T)=0 \mbox{ and for all } t\in(0,T)\,: \  a(t)\in \Honedot   \mbox{ and for all } \phi\in \Honedot  :
\end{equation*}

\begin{equation}\label{eqn:a_weak}
\int_\Omega\Bigl[ (1-2\kappa p) a_{tt} \, \phi - (b \nabla a_t - c^2 \nabla a)\cdot \nabla \phi\Bigr]\, dx = \int_\Sigma y\, \phi\, dS\,,
\end{equation}

which is the variational form of

\begin{equation}\label{eqn:a}
\begin{aligned}
(1-2\kappa p) a_{tt} - b\mathcal{A} a_t + c^2\mathcal{A} a = 0 \mbox{ in }\Omega,\
[[\partial_{\nu_\Sigma} (b\,a_t - c^2a)]] = y \mbox{ on }\Sigma
\end{aligned}
\end{equation}

with homogeneous end and boundary conditions.

To see that (<ref>) is the weak form of (<ref>), we consider a partition of $\Omega$ into two subdomains, whose interface contains $\Sigma$,
see Figure fig:Sigma,
with $\Omega_{in}$, $\Omega_{out}$ open, $\Omega_{in}\cap\Omega_{out}=\emptyset$,
$\Sigma\subseteq \overline{\Omega}_{in}\cup\overline{\Omega}_{out}$,
so that the outer normal vector $\nu$ coincides with $\nu_\Sigma$ on $\Sigma\cup\overline{\Omega}_{in}$ and with $-\nu_\Sigma$ on $\Sigma\cup\overline{\Omega}_{out}$ and integrate by parts to obtain

The surface $\Sigma$

\[
\begin{aligned}
&\int_\Omega \nabla a\cdot\nabla\phi \, dx =
\int_{\Omega_{in}} \nabla a\cdot\nabla\phi \, dx + \int_{\Omega_{out}} \nabla a\cdot\nabla\phi \, dx\\
&= \int_{\Omega_{in}} \mathcal{A} a \, \phi \, dx + \int_{\Omega_{out}} \mathcal{A} a \, \phi \, dx
- \int_{\partial\Omega_{in}}\partial_\nu a \, \phi\, dS - \int_{\partial\Omega_{out}}\partial_\nu a \, \phi\, dS\\
&= \int_\Omega \mathcal{A} a\,\phi \, dx - \int_\Sigma [[\partial_{\nu_\Sigma} a]] \, \phi\, dS\,,
\end{aligned}
\]

where ${\displaystyle [[\partial_{\nu_\Sigma} a]](x)
=\Bigl(\lim_{\tilde{x}\to x,\, \tilde{x}\in\Omega_{in}}\nabla a - \lim_{\tilde{x}\to x,\, \tilde{x}\in\Omega_{out}}\nabla a\Bigr)\cdot \nu_\Sigma}$.

With $a$ according to (<ref>), we can indeed establish the identity (<ref>): Inserting (<ref>) with $\phi=z(t)$ into (<ref>) and integrating over $(0,T)$, using integration by parts and the initial/end conditions on $z$ and $a$, respectively, we obtain

\[
\begin{aligned}
\langle F'(\kappa)\underline{d\kappa}, y\rangle_Y =&
\int_0^T\int_\Omega\Bigl[ (1-2\kappa p) a_{tt} \, z - (b \nabla a_t - c^2 a)\cdot \nabla z\Bigr]\, dx\, dt\\
\int_0^T\int_\Omega\Bigl[ ((1-2\kappa p) z)_{tt} \, a + (b \nabla z_t + c^2 z)\cdot \nabla a\Bigr]\, dx\, dt\\
=& \int_0^T\int_\Omega\underline{d\kappa}(p^2)_{tt}\, a \, dx\, dt
=\int_\Omega\underline{d\kappa} \int_0^T(p^2)_{tt}\, a \, dt\, dx\\
=&\int_\Omega\mathcal{A}^{s/2} \underline{d\kappa} \, \int_\Omega\mathcal{A}^{s/2}\Bigl[\mathcal{A}^{-s}\Bigl(\int_0^T(p^2)_{tt}\, a \, dt\Bigr)\Bigr]\, dx\,,
\end{aligned}
\]

where we have also used the weak form of (<ref>) with $a$ as a test function.
Thus with $\langle \kappa, \mu\rangle_X=: \int_\Omega \mathcal{A}^{s/2} \kappa\, \mathcal{A}^{s/2} \mu\,dx$ (where we can replace $\mathcal{A}^{s/2}$ by $\nabla$ in case $s=1$ and $\mathcal{A}$ is equipped with homogeneous Dirichlet conditions), we get the following.

Under the assumptions of Proposition prop:forward_diff, the $\dot{H}^s(\Omega) - L^2(0,T;L^2(\Sigma))$ adjoint of $F'(\kappa)$ is given by $F'(\kappa)^*y = \mathcal{A}^{-s}\Bigl(\int_0^T(p^2)_{tt}\, a \, dt\Bigr)$, where $a$ solves (<ref>) with $a(T)=a_t(T)=0$.

§.§ Weak sequential closedness

Application of convergence results for variational regularisation methods like Tikhonov, Ivanov, or Morozov regularisation [47, 13, 19, 25, 30, 35, 37] requires $F$ to be weak* sequentially closed, i.e.,

\begin{equation}\label{eqn:Fseqclos}
\Bigl(\kappa_j\stackrel{*}{\rightharpoonup} \kappa \mbox{ in }X\ \mbox{ and }
F(\kappa_j)\stackrel{*}{\rightarrow} y \mbox{ in }Y\Bigr) \ \Longrightarrow \
\Bigl(\kappa \in \mathcal{D}(F) \mbox{ and } F(\kappa)=y\Bigr)\,.
\end{equation}

Here we return to $X=L^\infty(\Omega)$, which
also implies (<ref>) with  $X=\dot{H}^s(\Omega)$ for $s>d/2$.
We verify (<ref>) by showing

For any $u_0\in\Htwodot $, $u_1\in \Honedot  $, $r\in L^2(0,T;L^2(\Omega))$, $T>0$ any $C^{1,1}$ domain $\Omega$ and $R$, $\mathcal{D}(F)$ as in Proposition prop:forward_welldef, the operator $G:\mathcal{D}(F)\to Y$, $G(\kappa)=p$  solving (<ref>) satisfies

\begin{equation}\label{eqn:Gseqclos}
\Bigl(\kappa_j\stackrel{*}{\rightharpoonup} \kappa \mbox{ in }X\ \mbox{ and }
G(\kappa_j)\stackrel{*}{\rightharpoonup} g \mbox{ in }V\Bigr) \ \Longrightarrow \
\Bigl(\kappa \in \mathcal{D}(F) \mbox{ and } G(\kappa)=g\Bigr)\,.
\end{equation}

Thus for any $C^{1,1}$ manifold $\Sigma\subseteq\overline{\Omega}$, $Y$ as in Proposition prop:forward_welldef, the forward operator $F:\mathcal{D}(F)\to Y$, $F(\kappa)=\mbox{tr}_\Sigma p$ satisfies (<ref>).

Abbreviating $\tilde{p}_j=G(\kappa_j)$, $p=G(\kappa)$, $v_j=\tilde{p}_j-p$, $\underline{d\kappa}_j=\kappa_j-\kappa$, we get, analogously to (<ref>) that for any $\phi\in C_0^\infty(\Omega)$, $\psi\in C_0^\infty(0,T)$

\begin{equation}\label{eqn:vj}
\begin{aligned}
(1-2\kappa p)v_{j,tt}+c^2\mathcal{A} v +b\mathcal{A} v_{j,t} \\
&\hspace*{3cm}- 2 \kappa (\tilde{p}_{j,t}+ p_t)\, v_{j,t} - 2\kappa \tilde{p}_{j,tt} \, v_j
- \underline{d\kappa}_j(\tilde{p}_j^2)_{tt}\Bigr)\, dx\, dt=0\,.
\end{aligned}
\end{equation}

For the linear (with respect to the sequences) terms we clearly get convergence from $v_j\stackrel{*}{\rightharpoonup} g-p=:\underline{dp}$  in $V$, that is,

\[
\begin{aligned}
(1-2\kappa p)v_{j,tt}+c^2\mathcal{A} v +b\mathcal{A} v_{j,t} - 2 \kappa p_t\, v_{j,t} - 2\kappa p_{tt} \, v_j\Bigr)\, dx\, dt\\
&\ \to \
\int_0^T\int_\Omega\psi(t)\phi(x)\Bigl(
(1-2\kappa p)\underline{dp}_{tt}+c^2\mathcal{A} \underline{dp} +b\mathcal{A} \underline{dp}_t - 2 \kappa p_t\, \underline{dp}_t - 2\kappa p_{tt} \, \underline{dp}\Bigr)\, dx\, dt
\,.
\end{aligned}
\]

It remains to consider the nonlinear terms.

For $\int_0^T\int_\Omega\psi(t)\phi(x) \kappa \tilde{p}_{j,t}\, v_{j,t} \, dx\, dt$,
we can make use of compactness of the embedding $V\to W^{1,4}(0,T;L^4(\Omega))$
so that there exists a subsequence $j_\ell$ along which the convergence
$\tilde{p}_{j,t}\to g_t$ and $v_j\to \underline{dp}$ is actually strong in
$L^4(0,T;L^4(\Omega))$ and we get

\[
\begin{aligned}
&\Big|\int_0^T\int_\Omega\psi(t)\phi(x) \kappa \Bigl(\tilde{p}_{j_\ell,t}\, v_{j_\ell,t} - g_t\underline{dp}_t\Bigr)\, dx\, dt\Big|\\
&\leq \|\kappa\|_{L^\infty(\Omega)} \|\phi\|_{L^2(\Omega)} \|\psi\|_{L^2(0,T)}
\Bigl(
\|\tilde{p}_{j_\ell,t}\|_{L^4(0,T;L^4(\Omega))} \|v_{j_\ell,t} - \underline{dp}_t\|_{L^4(0,T;L^4(\Omega))}\\
+ \|\tilde{p}_{j_\ell,t} - g_t\|_{L^4(0,T;L^4(\Omega))} \|\underline{dp}_t\|_{L^4(0,T;L^4(\Omega))}\Bigr)
\to0
\end{aligned}
\]

as $\ell\to\infty$. A subsequence-subsequence argument yields convergence of
$(x) \kappa \tilde{p}_{j,t}\, v_{j,t} \, dx\, dt$
to $\int_0^T\int_\Omega\psi(t)\phi(x) \kappa g_t\, \underline{dp} \, dx\, dt$.

We cannot make immediate use of such a compactness argument for
$(x) \kappa \tilde{p}_{j,tt} \, v_j \, dx\, dt$,
since $\tilde{p}_{j,tt}$ is just in $L^2(0,T;L^2(\Omega))$ and no more.
Instead we integrate by parts with respect to time to obtain,
for any $W^{1,4}(0,T;L^4(\Omega))$ convergent subsequence

\[
\begin{aligned}
&\int_0^T\int_\Omega\psi(t)\phi(x) \kappa \Bigl(\tilde{p}_{j_\ell,tt}\, v_{j_\ell} - g_{tt}\underline{dp}\Bigr)\, dx\, dt\\
\psi(t)\phi(x) \kappa \Bigl(\tilde{p}_{j_\ell,t}\, v_{j_\ell,t} - g_t\underline{dp}_t\Bigr)
+\psi'(t)\phi(x) \kappa \Bigl(\tilde{p}_{j_\ell,t}\, v_j - g_t\underline{dp}\Bigr)
\Bigr\}\, dx\, dt\,,
\end{aligned}
\]

where the first term can be tackled with exactly the same
$L^4(0,T;L^4(\Omega))$ convergence argument as above, and the second term can
be estimated even slightly easier in an analogous way:

\[
\begin{aligned}
&\Big|\int_0^T\int_\Omega\psi'(t)\phi(x) \kappa \Bigl(\tilde{p}_{j_\ell,t}\, v_{j_\ell} - g_t\underline{dp}\Bigr)\, dx\, dt\Big|\\
&\quad\leq \|\kappa\|_{L^\infty(\Omega)} \|\phi\|_{L^2(\Omega)} \|\psi'\|_{L^2(0,T)}
\Bigl(
\|\tilde{p}_{j_\ell,t}\|_{L^4(0,T;L^4(\Omega))} \|v_{j_\ell} - \underline{dp}\|_{L^4(0,T;L^4(\Omega))}\\
+ \|\tilde{p}_{j_\ell,t} - g_t\|_{L^4(0,T;L^4(\Omega))} \|\underline{dp}\|_{L^4(0,T;L^4(\Omega))}\Bigr)
\to0\qquad \mbox{as\ }
\ell\to\infty.
\end{aligned}
\]

Similarly, for
$\int_0^T\int_\Omega\psi(t)\phi(x) \underline{d\kappa}_j(\tilde{p}_j^2)_{tt} \, dx\, dt$ we get

\[
\begin{aligned}
&\Big|\int_0^T\int_\Omega\psi(t)\phi(x) \underline{d\kappa}_j (\tilde{p}_j^2)_{tt}\, dx\, dt\Big|
=\Big|\int_0^T\int_\Omega\psi''(t)\phi(x) \underline{d\kappa}_j \tilde{p}_j^2\, dx\, dt\Big|\\
&\leq\Big|\int_0^T\int_\Omega\psi''(t)\phi(x) \underline{d\kappa}_j g^2\, dx\, dt\Big|\\
&\qquad+\|\psi''\|_{L^\infty(0,T)} \|\phi\|_{L^\infty(\Omega)} \|\underline{d\kappa}_j\|_{L^\infty(\Omega)}
\int_0^T\int_\Omega |\tilde{p}_j^2-g^2|\, dx\, dt\\
\end{aligned}
\]

where the first term tends to zero by weak convergence of $\underline{d\kappa}_j$ to zero and the second one by $L^4(0,T;L^4(\Omega))$ convergence of $\tilde{p}_j$ to $g$ along subsequences.

Taking the limit $j\to\infty$ in (<ref>) and making use of the fact that $\phi\in C_0^\infty(\Omega)$, $\psi\in C_0^\infty(0,T)$ are arbitrary, we therefore get

\[
(1-2\kappa p)\underline{dp}_{tt}+c^2\mathcal{A} \underline{dp} +b\mathcal{A} \underline{dp}_t - 2 \kappa (p_t+g_t)\, \underline{dp}_t - 2\kappa g_{tt} \, \underline{dp}
= 0\,.
\]

(This identity actually holds in $L^2(0,T;L^2(\Omega)$.)
Moreover, like $\tilde{p}_j$ and $p$, also $g$ satisfies the initial conditions $g(0)=p_0$, $g_t(0)=p_1$, so $\underline{dp}=g-p$ satisfies homogeneous initial conditions.
Lemma lem:enest – which is applicable due to $v,g\in V$, $p\in M$ – together with Gronwall's inequality therefore yields $\underline{dp}=0$, thus $g=p$.

The conclusion (<ref>) on $F$ follows by boundedness of the trace operator $\mbox{tr}_\Sigma $, which is linear and therefore weak* continuous.

We mention in passing that due to the strong damping present in the equation, one might, alternatively to the second order wave equation, use a first order in time heat equation type reformulation.
Integrating with respect to time we get, in place of (<ref>),
the parabolic pde with memory

\begin{equation}\label{eqn:Westervelt_p_int}
p_t - b\triangle p -c^2\triangle \int_0^t p(\tau)\, d\tau = \kappa(x)(p^2)_t + R(x,t)
\end{equation}

where $R(x,t)=\int_0^t r(x,\tau)\, d\tau$.
Note that the right hand side also contains the first time derivative of the state, so we rewrite the equation as

\[
(1-2\kappa(x)p) p_t - b\triangle p -c^2\int_0^t\triangle p(\tau)\, d\tau =  R(x,t)
\]

Although the underlying pde is clearly supposed to model wave propagation,
this parabolic reformulation is justified mathematically by the fact that due to the strong damping, the
linearisation of (<ref>) is known to give rise to an analytic
semigroup and maximal parabolic regularity in the appropriate spaces, cf.
[29, 36].
We make use of this reformulation mainly for numerical computations where it allows
us to adapt a highly efficient Crank-Nicolson solver for the forward problem.

§ UNIQUENESS AND ILL-POSEDNESS

In either of the formulations the nonlocal parabolic (<ref>)
or the hyperbolic (<ref>) a time trace data can be expected to give valuable
information on unknown spatially-dependent coefficients occurring in the
operator.  For linear equations this idea goes back to at least the work of
Pierce [43], where the unique determination of a potential
$q$ in $\,u_t -u_{xx} + q(x) u = 0$ could be obtained from measurements of the
overposed value of $h(t) = u(1,t)$ give Neumann base conditions at that point.
The technique was to convert the problem to one of Sturm-Liouville type and
thus the method worked for all the physically important coefficients in a
parabolic equation in one space dimension by means of the Liouville transform.
As we will see later in this section the recovery of the spectral information
relies on an analytic continuation argument.
In its practical aspect it amounts to recovering the coefficients
and exponents in a Dirichlet series representing an analytic function
(which is merely sampled at a discrete set of points).
Even after this has been accomplished what must be recovered is the variation
of the actual sequence of eigenvalues of the operator from their expected
asymptotic form given by the Weyl formula.
However,  this means the information content is a decreasing sequence
(at least in $\ell^2$ but with additional smoothness on the
coefficients can have a much faster decay) that has to be
recovered when masked by an increasing sequence that is quadratically growing.
See, [45] for details.
Thus even in the case of a linear equation the recovery of a spatial
coefficient from time trace data is exponentially ill-conditioned and
any constants that one might hope work in ones favour turn out to be in

However the method breaks down completely in the case of nonlinear equations.
In particular, our unknown coefficient $\kappa$ appears intrinsically coupled
to the nonlinearity.
Thus instead we shall make do with a local injectivity argument for the
recovery of $\kappa$ leaving the broader issue of actual uniqueness
of the inverse problem for future work.
Although this is certainly a lesser result,
it does provide hope and insight on possible uniqueness
for the original nonlinear inverse problems and is also crucial for
linearisation methods such as  Newton and Halley.
In particular, we shall show that it leads to the frozen Newton and
frozen Halley schemes being well-defined and justifies our reconstruction
methods to be presented in Section sec:reconstructions.

From (<ref>), (<ref>) we have $F'(0)\underline{d\kappa} = \mbox{tr}_\Sigma z_0$ where $p_0$ and $z_0$ solve

\begin{equation}\label{eqn:p0}
p_{0,tt}+c^2\mathcal{A} p_0 + b\mathcal{A} p_{0,t} = r
\end{equation}

\begin{equation}\label{eqn:z0}
z_{0,tt}+c^2\mathcal{A} z_0 +b\mathcal{A} z_{0,t}
= \underline{d\kappa}(p_0^2)_{tt}\,,
\end{equation}

both with homogeneous initial conditions.

In setting up an experiment designed to recover $\kappa$,
since the initial conditions are zero, we can use non-homogeneous boundary
conditions to drive the system or a nonzero forcing function $r(x,t)$.
We assume here that the latter has been taken and specifically that $r$
has the form

\begin{equation}\label{eqn:choicer}
r(x,t) = f(x) \beta''(t) + \mathcal{A}f(x)(c^2\beta(t) + b\beta'(t))
\end{equation}

with some function $f$ in the domain of $\mathcal{A}$ vanishing only on a set of
measure zero and some twice differentiable function $\beta$ of time such that $(\beta^2)''(t_0)\not=0$ for some $t_0>0$.
Considering, for example the simple case $\beta(t)=t$, one can see that due to
the properties of solutions to Poisson's equation (that is, $\triangle f\leq0$ in $int(\Gamma)\subset\Omega$ and  $f > 0$ on $\Gamma$ implies $f>0$ in $int(\Gamma)$ by the maximum principle), this is possible even with
excitations concentrated on a surface $\Gamma$ like those pointed to in the introduction.
With (<ref>), the solution $p_0$ of equation (<ref>) is clearly given by
$p_0(x,t) = f(x) \beta(t)$.

With this, now $\underline{d\kappa}(p_0^2)_{tt}$ can be written in the form

\begin{equation}\label{eqn:p0_fac2}
\bigl(\underline{d\kappa}(x)p_0^2(x,t)\bigr)_{tt} = \sum_{j=1}^\infty a_j\phi_j(x)\psi_j(t)
\end{equation}

where $a_j$ are the coefficients of $\underline{d\kappa} f$ with respect to the basis $(\phi_j)_{\j\in\mathbb{N}}$, $\psi_j=(\beta^2)''$,

and the eigensystem $(\lambda_j,\phi_j)_{j\in\mathbb{N}}$ of the elliptic operator
$\mathcal{A}$ and the functions $\psi_j$ are known and we normalise $\phi_j$
to have $L^2$ norm unity.
Note that $r$ is the chosen excitation of the system. So the above condition (<ref>) – although it might appear as a restriction at first glance – actually gives a hint on how to select the input $r$ to enable sufficient sensitivity of the observations with respect to the searched for coefficient.
Also note that condition (<ref>) might be achieved even in more general settings than (<ref>).

rep C, 4.

We can rewrite equation (<ref>) as

\begin{equation}\label{eqn:z0j2}
z_j''(t)+c^2\lambda_j z_j(t) +b\lambda_j z_j'(t)
= a_j \psi_j(t)\,, \quad t>0\,, \qquad z_j(0)=0\,, \ z_j'(0)=0
\end{equation}

for all $j\in\mathbb{N}$, where

\[
z_0(x,t)=\sum_{j=1}^\infty z_j(t)\phi_j(x)\, .
\qquad
\]

Applying the Laplace transform to both sides of (<ref>) yields

\begin{equation}\label{eqn:z0jhat2}
\hat{z}_j(s)= \hat w_j(s) a_j \widehat{\psi}_j(s), \qquad
\mbox{where }\ \hat w_j(s) = \frac{1}{s^2+c^2\lambda_j  +b\lambda_j s}
\,, \qquad s\in\mathbb{C}\,,
\end{equation}

where we have used homogeneity of the initial conditions.
The poles of the function $\hat{w}_j(s)=\frac{1}{s^2+c^2\lambda_j  +b\lambda_j s}$ can be explicitly computed

\[
p_{j,\pm} = \tfrac12\left(-b\pm\sqrt{b^2-\tfrac{4c^2}{\lambda_j}}\right)\lambda_j
= \tfrac{2c^2}{-b\mp\sqrt{b^2-\tfrac{4c^2}{\lambda_j}}}
\]

are single, real up to finitely many complex conjugate pairs and they lie strictly in the left half of the complex plane, accumulating only at $-\tfrac{c^2}{b}$ and $-\infty$, more precisely,

\[
p_{j,+} \to -\tfrac{c^2}{b}\,, \quad p_{j,-} \to -\infty \mbox{ with } \frac{p_{j,-}}{\lambda_j}\to-b
\mbox{ as }j\to\infty.
\]

Moreover, they are different for different $\lambda_j$, i.e.,

\begin{equation}
p_{j,+}=p_{k,+} \ \Rightarrow \ \lambda_j=\lambda_k\,.
\end{equation}

Indeed, if $p_{j,+}=p_{k,+}$ then $-\lambda_k(c^2+bp_{k,+})=p_{k,+}^2=p_{j,+}^2=-\lambda_j(c^2+bp_{j,+})=-\lambda_j(c^2+bp_{k,+})$, thus $(\lambda_j-\lambda_k)(c^2+bp_{k,+})=0$ where $(c^2+bp_{k,+})$ does not vanish for otherwise $p_{k,+}^2$ would vanish.

Thus, assuming that $F'(0)\underline{d\kappa} = \mbox{tr}_\Sigma z_0=0$ implies that

\[
0 = \hat{z}_0(x_0,s) = \sum_{j=1}^\infty a_j \phi_j(x_0) \hat{w}_j(s) \widehat{\psi}_j(s)\,, \qquad
\mbox{ for all }s\in\mathbb{C}\,, \ x_0\in  \Sigma\,.
\]

Considering the residues at some pole $p_{m,+}$ yields

\[
\begin{aligned}
0 &= \mbox{Res}(\hat{z}_0(x_0;p_{m,+})
= \sum_{j=1}^\infty a_j \phi_j(x_0) \lim_{s\to p_{m,+}}(s-p_{m,+}) \hat{w}_j(s) \widehat{\psi}_j(s)\\
&=  \sum_{k\in K_m} a_k \phi_k(x_0) \mbox{Res}(\hat{w}_k;p_{m,+}) \widehat{\psi}_k(p_{m,+})
=  \frac{\widehat{\psi}_k(p_{m,+})}{\sqrt{b^2-\tfrac{4c^2}{\lambda_m}}\lambda_m} \sum_{k\in K_m} a_k \phi_k(x_0) \,,
\end{aligned}
\]

where $K_m\subseteq\mathbb{N}$ is an enumeration of the eigenspace basis $(\phi_k)_{k\in K_m}$ corresponding to the eigenvalue $\lambda_m$.
Assuming now that

\begin{equation}\label{eqn:ass_inj_psi}
\widehat{\psi}_k(p_{m,+})\not=0
\end{equation}

and there exists points $x_{0,m,1}, \ldots x_{0,m,N_m}\in \Sigma$, $N_m\geq \# K_m$ such that

\begin{equation}\label{eqn:ass_inj_Sigma}
\mbox{the matrix }\phi_k(x_{0,m,i})_{k\in K_m,i\in\{1,\ldots,N_m\}} \mbox{ has full rank }\# K_m
\end{equation}

we can conclude that $a_k=0$ for all $k\in K_m$.
The same argument goes through with $p_{m,-}$ in place of $p_{m,+}$.

Now since $(p_0^2)_{tt}(t_0)=f (\beta^2)''(t_0)$ only vanishes on a set of measure zero
we can conclude that $\underline{d\kappa}=0$ almost everywhere.

Under the above assumptions (<ref>), (<ref>), (<ref>) for all $m\in\mathbb{N}$, $k\in K_m$, the linearised derivative at $\kappa=0$, $F'(0)$ is injective.

In particular, (<ref>) is satisfied in the spatially 1-dimensional case $\Sigma=\{x_0\}$, where all eigenvalues of $\mathcal{A}$ are single, i.e., $\# K_m=1$ for all $m$, provided none of the eigenfunctions vanish at $x_0$; this can be achieved by taking $x_0$ on the boundary and where $\phi_j$ is subject to non-Dirichlet conditions.

Note that the very reasonable assumption that the eigenvalues of the
operator $\mathcal{A}$ are known is in fact partially redundant:
these can be obtained from the above argument in the case of one spatial
dimension and up to multiplicity in the case of higher space dimensions.

Although the ability to obtain the pole locations and hence $\{\lambda_j\}$
given the values of $b$ and $c^2$
from the functions $\hat w_j(s)$ appearing in equation (<ref>)
is redundant, we can nevertheless use the fact that $\{\lambda_j\}$ is
already assumed known to obtain the values of both $b$ and $c$ from
the Laplace transform of the overposed data $h(t)$.
Once again, this is an argument based on analytic continuation of
an analytic function and hence severely ill-conditioned.

§ REGULARIZATION METHODS

In this section we revisit some of the commonly used regularisation paradigms and discuss their application to the inverse problems of nonlinearity imaging.
To this end, we will write the inverse problem as an operator equation

\[
\]

where $F:X\to L^2(0,T;L^2(\Sigma))$ is the forward operator analyzed in Section sec:analysis_forward and we consider the choices $X=\dot{H}^s(\Omega)$ or $X=L^\infty(\Omega)$. In place of the exact time trace data $h$, we are given a noisy version $h^\delta$ or actually usually only its sample at a finite number of time instances, which is the setting we use in our numerical experiments in Section sec:reconstructions. The superscript $\delta$ indicates the noise level in

\begin{equation}\label{eqn:delta}
\|F(\kappa_{act})-h^\delta\|\leq\delta
\end{equation}

where $\kappa_{act}$ denotes an exact solution of the inverse problem, i.e., such that $F(\kappa_{act})=h$.

§.§ Tikhonov regularisation

rep C, 7.
Certainly the most well-known, most widely used, and – via the choice of data misfit and regularisation functionals – most versatile approach is Tikhonov-Philips regularisation.

Using norms for both functionals we define a regulariser $\kappa_\alpha$ by minimizing $J_\alpha$ given by

\[
\]

where $\kappa_0$ is an a priori guess and we either choose
$X=\dot{H}^s(\Omega)$ with $s>\frac{d}{2}$ to be a Hilbert space that is
embedded in $L^\infty(\Omega)$, or directly use the Banach space
$X=L^\infty(\Omega)$, although this makes the minimization of $J_\alpha$
more challenging, due to additional factors of nonlinearity and nonsmoothness.

The actual computation of the Tikhonov regulariser will usually be based on
descent methods using the gradient of the cost function $J_\alpha$.
The use of the Hessian in principle gives better search directions but comes
at a high computational cost, as we also point out in
Section sec:reconstructions in the context of Halley's method.

The gradient of $J_\alpha$ in the case of $X=\dot{H}^s(\Omega)$ is given by

\[
\]

rep C, 5.
with the Hilbert space adjoint as in Proposition prop:adj,
see also the lines following (<ref>) for its implementation in the context of the inverse problem under consideration here.

For the use of $X=L^\infty(\Omega)$, considering the special case $s=0$
in Proposition prop:adj, we immediately obtain the
$L^\infty(\Omega)$-$L^2(0,T;L^2(\Omega))$ Banach space adjoint of $F'(\kappa)$,
which by  applying a duality mapping on $L^\infty(\Omega)$, in principle allows
one to construct the gradient of $J_\alpha$.

However, convergence of gradient type methods in nonreflexive Banach spaces is in general a highly nontrivial question.
If $X$ is the dual of a separable space, as is the case for $X=L^\infty(\Omega)$, the Tikhonov minimization problem can be tackled by means of a duality approach. This also naturally leads to an appropriate discretization that lends itself to iterative minimization by, e.g., a semismooth Newton method, see [10, 9].

rep C, 6.

Propositions prop:forward_welldef, prop:Fseqclos together with <cit.>, <cit.> or <cit.> respectively, yield well-definedness and convergence of Tikhonov regularisation.

Tikhonov regularisation with $X=\dot{H}^s(\Omega)$, $s>\frac{d}{2}$ or $X=L^\infty(\Omega)$ is well-defined and converges in the sense that for data $h^\delta$ satisfying

\begin{equation*} \|h^\delta-\mbox{tr}_\Sigma p_{act}\|_{L^2(0,T;L^2(\Sigma))}\leq\delta
\end{equation*}

\begin{equation*}
\alpha(\delta)\to0, \frac{\delta}{\alpha(\delta)}\to0
\end{equation*}

we have subsequential[every sequence has a convergent subsequence and the limit of every subsequence is a solution to the inverse problem; in case of uniqueness of $\kappa_{act}$, the whole sequence converges] convergence of $\kappa_{\alpha(\delta)}$ to $\kappa_{act}$ as $\delta\to0$ where convergence takes place strongly in case of $X=\dot{H}^s(\Omega)$, and weakly* in case of $X=L^\infty(\Omega)$.

An alternative option is to use $X=L^2(\Omega)$ and constrain minimization to a subset $\{\kappa\in L^2(\Omega)\, : \, \underline{\kappa}\leq\kappa\leq\overline{\kappa}\}$ of $L^\infty(\Omega)$, see [40].

§.§ Iterative methods

As already mentioned in the previous section, any variational regularisation method such as Tikhonov, Ivanov or Morozov regularisation and versions thereof, that are based on minimization of some cost functional, require employment of some iterative descent algorithm for their numerical implementation.
Alternatively, one can directly use iterative solution methods, equipped with some regularisation, as reconstruction methods.

Due to the above mentioned difficulties with constructing duality mappings in $L^\infty(\Omega)$, we here consider iterative methods in a Hilbert space setting $X=\dot{H}^s(\Omega)$ only.

§.§.§ Landweber iteration

The most simple iterative approach consists of performing gradient descent for the least squares functional $\frac12\|F(\kappa)-h^\delta\|^2$, which results in Landweber iteration

\begin{equation}\label{eqn:LW}
\begin{aligned}
&\kappa_{n+1}=\kappa_n+\mu_n F'(\kappa_n)^*\bigl(h^\delta-F(\kappa_n)\bigr).
\end{aligned}
\end{equation}

Therefore, one step of Landweber iteration (see, e.g., [17]) reads as follows

* given $\kappa=\kappa_n$, solve (<ref>) to obtain $p$;

* set $y=h^\delta-p$;

* solve (<ref>) to obtain $a$ and set $\delta\kappa=\mathcal{A}^{-s}\int_0^Ta(t)(p^2)_{tt}(t)\, dt$;

* choose a stepsize $\mu_n$ and set $\kappa_{n+1}=\kappa_n+\mu_n\delta\kappa$.

Freezing the derivative at a fixed argument, e.g., the starting value, still yields a fixed point iteration, which one might call frozen Landweber:

\begin{equation}\label{eqn:frozenLW}
\begin{aligned}
&\kappa_{n+1}=\kappa_n+\mu_n F'(\kappa_0)^* (h^\delta-F(\kappa_n))
\,.
\end{aligned}
\end{equation}

By Propositions prop:forward_welldef, prop:forward_diff, prop:adj, Landweber iteration is well-defined.
As they stand, the Landweber steps do not require further regularisation (although there do exist modified versions that also contain regularisation terms, see, e.g., [46]).
Still, in order to avoid explosion of the noise propagated through the iterations, one has to stop after an appropriately chosen number of steps.
A well-established criterion for this purpose is the discrepancy principle, which terminates the iteration as soon as the residual is of the order of the noise level.

\begin{equation}\label{eqn:discrprinc}
k_*=\min\{k\in \mathbb{N}\, : \, \|F(\kappa)-h^\delta\|\leq\tau\delta\}\,,
\end{equation}

for some fixed safety factor $\tau>1$. (In our computations we used $\tau=2$ which is known to be the minimal choice for nonlinear problems [17].)
The availability of the noise level $\delta$ in practical applications might be
debatable; yet its necessity for establishing convergence guarantees is known
(as Bakushinski's veto [2], see also <cit.>).
The stepsize $\mu_k$ might be chosen as a sufficiently small constant cf. [17], or adapted in each step, cf., e.g., [41].

Landweber iteration is known to be notoriously slow. However, with the ingredients above, also state-of the art accelerated versions of Landweber iteration, such as the steepest descent or the minimal error method [41] or Nesterov iteration and versions thereof [21, 39] can be implemented in a straightforward manner.

Also note that alternatively to choosing $s$ is sufficiently large to
enforce continuity of the embedding $H^s(\Omega)\to L^\infty(\Omega)$,
we can just take $s=0$, that is, $X=L^2(\Omega)$, and project the regularised
iterates into an $L^\infty(\Omega)$ ball in order to guarantee well-definedness
of $F$ and $F'$ according to Propositions prop:forward_welldef,

§.§.§ Newton's method

Significantly faster methods result from using the first order Taylor expansion of the forward operator which results in Newton's method

\begin{equation}\label{eqn:Newton}
\begin{aligned}
&\kappa_{n+1}=\kappa_n+F'(\kappa_n)^{-1} (h-F(\kappa_n))
\end{aligned}
\end{equation}

or its frozen version

\begin{equation}\label{eqn:frozenNewton}
\begin{aligned}
&\kappa_{n+1}=\kappa_n+F'(\kappa_0)^{-1} (h-F(\kappa_n))
\,.
\end{aligned}
\end{equation}

In the context of nonlinearity imaging, $F(\kappa_n)=\mbox{tr}_\Sigma p$ with $p$ solving (<ref>) and $F'(\kappa_n)[\underline{d\kappa}]=\mbox{tr}_\Sigma z$ with $z$ solving (<ref>).
In particular, in the frozen Newton case at $\kappa=0$, both equations needed to evaluate $F'(0)[\underline{d\kappa}]=\mbox{tr}_\Sigma z_0$, namely (<ref>) and (<ref>) are linear, while of course we still have to solve the nonlinear pde to obtain $F(\kappa_n)$.

By Propositions prop:forward_welldef, prop:forward_diff, and Theorem th:lininj, Newton's method is well-defined; however, in case of noisy data $h^\delta\approx h$ it needs regularisation.
Regularized versions of Newton usually rely on Tikhonov's method applied to the linearised problem, which in a Hilbert space setting leads to the Levenberg-Marquardt method

\begin{equation}\label{eqn:regNewton}
\begin{aligned}
&\kappa_{n+1}=\kappa_n+(F'(\kappa_n)^*F'(\kappa_n)+\alpha_nI)^{-1} F'(\kappa_n)^*(h^\delta-F(\kappa_n))
\end{aligned}
\end{equation}

or its frozen version

\begin{equation}\label{eqn:regfrozenNewton}
\begin{aligned}
&\kappa_{n+1}=\kappa_n+(F'(\kappa_0)^*F'(\kappa_0)+\alpha_nI)^{-1} F'(\kappa_0)^*(h^\delta-F(\kappa_n))
\,.
\end{aligned}
\end{equation}

cf., e.g., [16, 44] or the iteratively regularised Gauss-Newton method, cf., e.g., [1, 5, 20].

Here $\alpha_k$ is a sequence of regularisation parameters that can be chosen a posteriori according to some inexact Newton strategy (see [16]) or simply as a geometrically decaying sequence $\alpha_n=\theta^n\alpha_0$.
The overall iteration again needs to be appropriately stopped, e.g. by the discrepancy principle (<ref>).

§.§.§ Halley's method

An even faster iterative method can be achieved by including Hessian information on the forward operator, which leads to Halley's method:

\begin{equation}\label{eqn:Halley}
\begin{aligned}
&\kappa_{n+\frac12}=\kappa_n+F'(\kappa_n)^{-1} (h-F(\kappa_n))\\
&\kappa_{n+1}=\kappa_n+(F'(\kappa_n)+\tfrac12 F''(\kappa_n)[\kappa_{n+\frac12}-\kappa_n,\cdot])^{-1} (h-F(\kappa_n))\,.
\end{aligned}
\end{equation}

Here it can also make sense to freeze evaluation of $F'$ and $F''$ to some
point $\kappa_0$:

\begin{equation}\label{eqn:frozenHalley}
\begin{aligned}
&\kappa_{n+\frac12}=\kappa_n+F'(\kappa_0)^{-1} (h-F(\kappa_n))\\
&\kappa_{n+1}=\kappa_n+(F'(\kappa_0)+\tfrac12 F''(\kappa_0)[\kappa_{n+\frac12}-\kappa_n,\cdot])^{-1} (h-F(\kappa_n))\,.
\end{aligned}
\end{equation}

In this case, the function values and derivatives of $F$ can be computed as

* $F(\kappa_n)=\mbox{tr}_\Sigma p$ with $p$ solving (<ref>);

* $F'(\kappa_0)[\underline{d\kappa}]=\mbox{tr}_\Sigma z$ with $z$ solving (<ref>) with $\kappa=\kappa_0$, $p=p_0$;

* $F''(\kappa_0)[\underline{d\kappa}^{(1)},\underline{d\kappa}^{(2)}]=\mbox{tr}_\Sigma w$ with $w$ solving

\begin{equation}\label{eqn:Fpp}
\Bigl((1-2\kappa_0 p_0) w\Bigr)_{tt} + b\mathcal{A} w_t +c^2\mathcal{A} w =
2\Bigl( \kappa_0 \, z^{(1)} \, z^{(2)} + p_0 ( \underline{d\kappa}^{(1)} z^{(2)} + \underline{d\kappa}^{(2)} z^{(1)})\Bigr)_{tt}\,,
\end{equation}

and $z^{(i)}$ satisfying (<ref>) with $\kappa=\kappa_0$, $p=p_0$, $\underline{d\kappa}=\underline{d\kappa}^{(i)}$.

In particular note that the evaluation of $F'$ and $F''$ at $\kappa_n$ lead to linear initial boundary value problems for the same pde just with different right hand sides.

Halley's method also needs to be regularised when used with noisy data,
cf. [18, 26].

Some remarks on convergence of iterative regularisation methods for the inverse
problem under consideration are in order.
Restrictions on the nonlinearity of $F$ such as the tangential cone condition
are not likely to hold here, in view of the fact that only boundary
observations are available.
Thus, convergence of Newton or gradient type (Landweber) regularisation methods cannot be proven without assuming further regularity.
According to, e.g., <cit.>, for the regularised Landweber iteration and the iteratively regularized Gauss-Newton method, respectively,
rep C, 8.
there is still the option of obtaining convergence (with rates) under a source condition
$\kappa_{act}-\kappa_0=F'(\kappa)^*w$ for some $w\in Y$ assuming (instead of
the tangential cone condition) only Lipschitz continuity of $F'$, cf. Remark rem:contFprime.
However due to the severe ill-posedness, that is, the infinite smoothing
properties of $F'(\kappa)^*$, this would only yield convergence (with rates)
in case of infinitely smooth initial error

§ RECONSTRUCTIONS

In this section we show reconstructions of $\kappa$ from time trace data.
The spatial set will be the interval $[0,1]$ and we will take the measurement
point of $h(t)=p(1,t)$ to be the right-hand endpoint $x=1$.

Our numerical implementation uses (<ref>)
in the transformed version as in (<ref>)
and so treat it as a parabolic equation with nonlocal memory
term $c^2\int_0^t \triangle p(\tau)\,d\tau$.
A Crank-Nicolson integrator was used with an inner iteration loop
to handle the nonlinear term $\,-2\kappa\,p\,p_t$.
A Neumann boundary condition was imposed at the right hand endpoint;
the left hand condition could be Dirichlet, Neumann or impedance type.
Typically, in the physical model one would have zero initial conditions but
this isn't necessary for the mathematical formulation.

Data consisted of the time trace measurement $h(t) = p(1,t)$.
As a practical matter we used the above mentioned solver to obtain this data
and collected a sample at $50$ equally spaced points on the interval $[0,T]$.
Uniformly distributed random noise was then added to these values to obtain
$h_{\rm meas}(t)$.
This was then pre-filtered by smoothing and up-resolving to the working
resolution of the number of time points taken ($\sim$400 for the interval
$t\in [0,1]$) for the direct solver used in the inversion routine.
Our reconstructions are mainly
based on two noise levels: $0.1\%$ and $1\%$, which in view of the exponential
ill-posedness appear to be in the most reasonable range.

The unknown $\kappa$ was represented in terms of given basis functions.
Since we wish to make no constraints on the form of
$\kappa$ other than sufficient regularity and positivity,
we do not choose a basis with in-built restrictions as would be
obtained from an eigenfunction expansion.
Instead we used a radial basis set consisting of either
shifted Gaussian functions $b_j(x) := e^{-(x-x_j)^2/\sigma}$
centered at nodal points  $\{x_j\}$
and with width specified by the parameter $\sigma$ or,
since we also wish to reconstruct non-smooth $\kappa$,
we also used chapeau piecewise linear functions.
Many of the figures shown below use the latter type and the interval
$x\in [0,1]$ contained $41$ such basis functions.
We also show a reconstruction
of a piecewise constant $\kappa$ since this situation is physically
meaningful and for this case we also adopted a Haar basis.
In all cases the starting approximation for the iterative methods used
was the value $\kappa=0$.

The assumption of $1-2\kappa p$ being positive and bounded away from zero is easily fulfilled in practice, see, e.g., <cit.>. Also in our numerical experiments we chose coefficient functions $\kappa$ and excitations $r$ such that this condition holds.

From a physical standpoint $\kappa(x)$ should be nonnegative and in some cases
(for example in Figure fig:Newton_iteration_kappa6) we forced this
situation by truncating all negative values to zero at the end of
each iteration of the scheme.
Typically, the need for this decreased as the iteration scheme progressed.

The Newton scheme as described in Section sec:Newton, in particular
the frozen version about $\kappa=0$, performed reliably and
converged rapidly for a wide range of trial $\kappa$ functions.
Regularization was by the Tikhonov method with regularisation parameter
$\alpha_k=0.5^k\alpha_0$ and as a stopping criterion we used the
Discrepancy Principle (<ref>).

The Landweber scheme from Section sec:LW was also effective and
able to give reconstructions under higher noise levels than Newton.
We do not dwell on this situation as the convergence rate was indeed
notoriously slow and we did not implement any
of the possible acceleration schemes mentioned earlier.

We also show reconstructions using Halley's method, again frozen about a fixed
$\kappa$ value (here $\kappa=0$) taking the predictor-corrector approach
of [18] and as shown in equation (<ref>).

§.§ Newton iteration

Figures fig:Newton_iteration_kappa4 and fig:Newton_iteration_kappa6
below show reconstructions of  two actual $\kappa$ functions:
one $C^\infty$ and the other just piecewise linear.
In each case the frozen Newton scheme was used with the leftmost graphs
for the case of a noise level of $0.1\%$ and the rightmost for $1\%$ noise.
Note the very rapid convergence of the scheme.
In the case of $1\%$ noise the scheme was stopped by the Discrepancy Principle
at the second iteration, whereas with the lower noise level a third iteration
was reached.

Reconstructions of a smooth $\kappa(x)$ from time trace data at $\,x=1\,$
under 0.1%  (left) and 1% (right) noise using Newton's method.

Reconstructions of piecewise linear $\kappa(x)$ from time trace data at $\,x=1$
under 0.1%  (left) and 1% (right) noise using Newton's method.

As a final example for Newton's method we show a reconstruction of a piecewise
constant function in Figure fig:kappa9_reconstruction.
This uses both piecewise linear or chapeau basis functions and also a
Haar piecewise constant basis.
We were careful to ensure that the basis breakpoints did not align with the
discontinuities of $\kappa$.
Clearly, the reconstruction here is much poorer but still able to follow the
significant features of the actual $\kappa(x)$ function.

Reconstructions of a piecewise constant $\kappa(x)$ from time trace data at
$x=1$ under $0.1\%$ noise using Newton iteration.

§.§ Halley iteration

Theoretically, the frozen Newton scheme is only first order accurate and
thus it is natural to seek a higher order of convergence method.
Thus comparisons with Halley should be made here.
The frozen Halley with the added corrector step is second order and might
seem the obvious choice.
Implementation requires the computation of the Hessian and this if implemented
in the most straightforward way takes approximately $M$ times longer than
just computing the Jacobian where $M$ is the number of basis elements.
Thus, given the rapid convergence of the Newton scheme with smooth $\kappa$
the use of Halley to improve the convergence rate seems counterproductive
given what can be a considerable computational expense.
However, there is another factor to consider.
In many cases, (see for example, [18]) the use of the corrector can
sometimes lead to a slightly improved reconstruction and this is
the case here as we show in Figure fig:Newton_Halley_compare.

Comparison of Newton (in red) and Halley (in blue) final reconstructions under $0.1\%$ noise.

These two reconstructions (from the same data subject to $0.1\%$ noise)
are hardly distinguishable except for the extreme left hand endpoint
where Halley outperforms Newton.
The final errors:
$\;\|\kappa - \kappa_{\rm act}\|_{\infty}$,
$\;\|\kappa - \kappa_{\rm act}\|_{L^2}\,$
for Newton were $0.0084$ and $0.0059$ respectively.
The corresponding errors for Halley were
$0.0066$ and $0.0045$.

For this numerical run
there is a homogeneous Dirichlet condition at the left hand endpoint
and the small values of $p(x,t)$ near $x=0$ are a multiplier for the $\kappa$
to be reconstructed. With these boundary conditions reconstructions were
always worse at this endpoint.

Figure fig:Newton_Halley_norm_compare below shows another
Newton-Halley comparison this time using a piecewise linear function for
Here we took Neumann boundary conditions at both endpoints and so the above
issue of the small multiplier is lessened.
Again, there is a slight improvement using Halley's method.

Comparison of Newton (in red) and Halley (in blue) final reconstructions and
norm differences of the $n^{\rm th}$ iterate $\kappa_n$ and the actual $\kappa$.
Noise level was $1\%$.

The rightmost figure shows the $L^\infty$ norm difference between the
computed and the actual $\kappa$ at each iteration. Several features are worthy of note.

First, both Newton and Halley schemes wander on the first  few iterations
and this is much more noticeable in the Halley.
There is no stepsize control implemented and the actual $\kappa$ is quite
far from the initial approximation $\kappa=0$.
This effect would disappear with proper stepsize control.

Second, after this initial phase the Newton progresses with approximately
linear convergence whereas the Halley decreases roughly quadratically.
In this simulation the use of the Discrepancy Principle to terminate the
iteration process was turned off and so after a certain point the noise in
the data plays a role and the norm differences start to increase again.
This point is after about 15 iterations with Halley but much later with Newton;
somewhere around 35 iterations.

Third, note that after terminating the schemes by some effective mechanism
the overall error in the Halley method is less than that of the Newton by a
factor of approximately the same amount as noted for
Figure fig:Newton_Halley_compare.

Thus in summary here, the computational cost of a Halley implementation
and the need to compute the Hessian is not repaid in time to reach convergence
as opposed to Newton, but there is a relatively small but consistent advantage
of a superior reconstruction.

§.§ Landweber iteration

As a point of reference, Figure fig:Landweber_iteration_kappa4
shows the piecewise linear actual $\kappa$ reconstructed using Landweber
After 5,000 steps the scheme was still converging so computational efficiency
is orders of magnitude below either Newton or Halley methods.

The rate was better for a smooth $\kappa$ as shown in
Figure fig:Landweber_iteration_kappa6 but typical of this scheme
it lags orders of magnitude behind Newton-schemes if these are applicable to
the problem.

Reconstructions of a piecewise linear $\kappa(x)$ from time trace data at $x=1$
under $1\%$ noise using Landweber iteration.

The leftmost figure shows reconstructions of $\kappa(x)$
under $0.1\%$ noise using Landweber iteration.
The rightmost figure shows the decay of the norm
$\kappa_n(x)-\kappa_{\rm act}(x)$

The main advantage to the Landweber scheme is its greater robustness against
noise than the Newton scheme and certainly so against the Halley which requires
two regularisation constants, one for each of the predictor and the corrector

As noted in Section sec:LW, there are acceleration methods available
for this scheme, but given the performance of the Newton and Halley schemes
shown here, these seem the methods of choice.

§ ACKNOWLEDGMENTS

The work of the first author was supported by the Austrian Science Fund fwf
under the grants P30054 and DOC78. The work of the second author was supported
in part by the National Science Foundation through award dms-1620138. We thank both reviewers for their careful reading of the manuscript and their detailed reports with valuable comments and suggestions that have led to an improved version of the paper.

[1] (MR1185952)
A. B. Bakushinskiĭ,

On a convergence problem of the iterative-regularised Gauss-Newton method,
Comput. Math. Math. Phys., 32 (1992), 1353–1359.

[2] [10.1016/0041-5553(84)90253-2]
A. B. Bakushinskii,
Remarks on choosing a regularisation parameter using the quasi-optimality and ratio criterion,
USSR Comput. Math. Math. Phys., 24 (1984), 181–182.
[3] [10.1016/0041-624x(86)90102-2]
L. Bjørnø,
Characterization of biological media by means of their non-linearity,
Ultrasonics, 24 (1986), 254–259.

[4]
D. T. Blackstock,
Approximate equations governing finite-amplitude sound in thermoviscous fluids,
Tech Report, GD/E Report, GD-1463-52, General Dynamics Corp., Rochester, NY, 1963.

[5] (MR1459331) [10.1093/imanum/17.3.421]
B. Blaschke, A. Neubauer and O. Scherzer,
On convergence rates for the iteratively regularised Gauss-Newton method,
IMA J. Numer. Anal., 17 (1997), 421–436.

[6] [10.1007/978-94-010-1745-9]
J. M. Burgers,
The Nonlinear Diffusion Equation,
Springer, Netherlands, 1974.

[7]
V. Burov, I. Gurinovich, O. Rudenko and E. Tagunov,
Reconstruction of the spatial distribution of the nonlinearity parameter and sound velocity in acoustic nonlinear tomography,
Acoustical Physics, 40 (1994), 816–823.

[8] [10.1109/ULTSYM.1985.198640]
C. A. Cain,
Ultrasonic reflection mode imaging of the nonlinear parameter B/A: A theoretical basis,
IEEE 1985 Ultrasonics Symposium, San Francisco, CA, USA, 1985.

[9] (MR3859297) [10.1515/jiip-2018-0026]
C. Clason and A. Klassen,
Quasi-solution of linear inverse problems in non-reflexive Banach spaces,
J. Inverse Ill-Posed Probl., 26 (2018), 689–702.

[10] (MR2775195) [10.1051/cocv/2010003]
C. Clason and K. Kunisch,
A duality-based approach to elliptic control problems in non-reflexive Banach spaces,
ESAIM Control Optim. Calc. Var., 17 (2011), 243–266.

[11] [10.1146/annurev.fl.11.010179.000303]
D. G. Crighton,
Model equations of nonlinear acoustics,
Ann. Rev. Fluid Mech., 11 (1979), 11–33.

[12] (MR1408680)
H. W. Engl, M. Hanke and A. Neubauer,
Regularization of Inverse Problems,
Mathematics and its Applications, 375, Kluwer Academic Publishers Group, Dordrecht, 1996.

[13] (MR1009037) [10.1088/0266-5611/5/4/007]
H. W. Engl, K. Kunisch and A. Neubauer,
Convergence rates for Tikhonov regularisation of non-linear ill-posed problems,
Inverse Problems, 5 (1989), 523–540.

[14] (MR1625845)
L. C. Evans,
Partial Differential Equations,
Graduate Studies in Mathematics, 19,  American Mathematical Society, Providence, RI, 1998.

[15]
M. F. Hamilton and D. T. Blackstock,
Nonlinear Acoustics, Vol. 1,
Academic Press, San Diego, 1998.

[16] (MR1435869) [10.1088/0266-5611/13/1/007]
M. Hanke,
A regularizing Levenberg–Marquardt scheme, with applications to inverse groundwater filtration problems,
Inverse Problems, 13 (1997), 79–95.

[17] (MR1359706) [10.1007/s002110050158]
M. Hanke, A. Neubauer and O. Scherzer,
A convergence analysis of the Landweber iteration for nonlinear ill-posed problems,
Numer. Math., 72 (1995), 21–37.

[18] (MR1740766) [10.1137/S0036142998341246]
F. Hettlich and W. Rundell,
A second degree method for nonlinear inverse problems,
SIAM J. Numer. Anal., 37 (2000), 587–620.

[19] (MR2329928) [10.1088/0266-5611/23/3/009]
B. Hofmann, B. Kaltenbacher, C. Pöschl and O. Scherzer,

A convergence rates result for Tikhonov regularisation in Banach spaces with non-smooth operators,
Inverse Problems, 23 (2007), 987–1010.

[20] (MR1474369) [10.1088/0266-5611/13/5/012]
T. Hohage,
Logarithmic convergence rates of the iteratively regularised Gauß-Newton method for an inverse potential and an inverse scattering problem,
Inverse Problems, 13 (1997), 1279–1299.

[21] (MR3830146) [10.1088/1361-6420/aacebe]
S. Hubmer and R. Ramlau,
Nesterov's accelerated gradient method for nonlinear ill-posed problems with a locally convex residual functional,
Inverse Problems, 34 (2018), 30pp.

[22] [10.1177/016173468300500401]
N. Ichida, T. Sato and M. Linzer,
Imaging the nonlinear ultrasonic parameter of a medium,
Ultrasonic Imaging, 5 (1983), 295–299.

[23] (MR4041949) [10.1088/1361-6420/ab323e]
O. Y. Imanuvilov and M. Yamamoto,
Carleman estimate and an inverse source problem for the Kelvin-Voigt model for viscoelasticity,
Inverse Problems, 35 (2019), 45pp.

[24] (MR2193218) [10.1007/0-387-32183-7]
V. Isakov,
Inverse Problems for Partial Differential Equations,
Applied Mathematical Sciences, 127, Springer, New York, 2006.

[25] (MR0140944)
V. K. Ivanov,
On linear problems which are not well-posed,
Dokl. Akad. Nauk SSSR, 145 (1962), 270–272.

[26] (MR3383327) [10.1007/s00211-014-0682-5]
B. Kaltenbacher,
An iteratively regularized Gauss-Newton-Halley method for solving nonlinear ill-posed problems,
Numer. Math., 131 (2015), 33–57.

[27] (MR3461695) [10.3934/eect.2015.4.447]
B. Kaltenbacher,
Mathematics of nonlinear acoustics,
Evol. Equ. Control Theory, 4 (2015), 447–491.

[28] [10.3934/eect.2020063]
B. Kaltenbacher,
Periodic solutions and multiharmonic expansions for the Westervelt equation,
to appear, Evol. Equ. Control Theory.

[29] (MR2525765) [10.3934/dcdss.2009.2.503]
B. Kaltenbacher and I. Lasiecka,
Global existence and exponential decay rates for the Westervelt equation,
Discrete Contin. Dyn. Syst. Ser. S, 2 (2009), 503–523.

[30] (MR3788158) [10.1088/1361-6420/aab739]
B. Kaltenbacher and A. Klassen,
On convergence and convergence rates for Ivanov and Morozov regularisation and application to some parameter identification problems in elliptic PDEs,
Inverse Problems, 34 (2018), 24pp.

[31] (MR2459012) [10.1515/9783110208276]
B. Kaltenbacher, A. Neubauer and O.Scherzer,
Iterative Regularization Methods for Nonlinear Ill-Posed Problems,
Radon Series on Computational and Applied Mathematics, 6, Walter de Gruyter GmbH & Co. KG, Berlin, 2008.

[32]
V. Kuznetsov,
Equations of nonlinear acoustics,
Soviet Physics - Acoustics, 16 (1971), 467–470.

[33] [10.1017/S0022112068000303]
M. B. Lesser and R. Seebass,
The structure of a weak shock wave undergoing reflexion from a wall,
J. Fluid Mech., 31 (1968), 501–528.

[34] (MR0077346)
M. J. Lighthill,
Viscosity effects in sound waves of finite amplitude,
in Surveys in Mechanics, Cambridge, at the University Press, 1956, 250–351.

[35] (MR3080477) [10.1088/0266-5611/29/7/075016]
D. Lorenz and N. Worliczek,
Necessary conditions for variational regularisation schemes,
Inverse Problems, 29 (2013), 19pp.

[36] (MR2822410) [10.1007/s00245-011-9138-9]
S. Meyer and M. Wilke,
Optimal regularity and long-time behavior of solutions for the Westervelt equation,
Appl. Math. Optim., 64 (2011), 257–271.

[37] (MR0208819)
V. A. Morozov,
On the solution of functional equations by the method of regularisation,
Soviet Math. Dokl., 7 (1966), 414–417.

[38] (MR3923817) [10.3934/eect.2019010]
M. Muhr, V. Nikolić, B. Wohlmuth and L. Wunderlich,
Isogeometric shape optimization for nonlinear ultrasound focusing,
Evol. Equ. Control Theory, 8 (2019), 163–202.

[39] (MR3652264) [10.1515/jiip-2016-0060]
A. Neubauer,
On Nesterov acceleration for Landweber iteration of linear ill-posed problems,
J. Inverse Ill-Posed Probl., 25 (2017), 381–390.

[40] (MR947434) [10.1016/0021-9045(88)90025-1]
A. Neubauer,
Tikhonov-regularisation of ill-posed linear operator equations on closed convex sets,
J. Approx. Theory, 53 (1988), 304–320.

[41] (MR1337265) [10.4171/ZAA/679]
A. Neubauer and O. Scherzer,
A convergent rate result for a steepest descent method and a minimal error method for the solution of nonlinear ill-posed problems,
Z. Anal. Anwendungen, 14 (1995), 369–377.

[42] (MR2027679) [10.1007/b97537]
H. Ockendon and J. R. Ockendon,
Waves and Compressible Flow,
Texts in Applied Mathematics, 47, Springer-Verlag, New York, 2004.

[43] (MR534419) [10.1137/0317035]
A. Pierce,
Unique identification of eigenvalues and coefficients in a parabolic problem,
SIAM J. Control Optim., 17 (1979), 494–499.

[44] (MR1826857) [10.1007/PL00005448]
A. Rieder,
On convergence rates of inexact Newton regularizations,
Numer. Math., 88 (2001), 347–365.

[45] (MR1106979) [10.1090/S0025-5718-1992-1106979-0]
W. Rundell and P. E. Sacks,
Reconstruction techniques for classical inverse Sturm-Liouville problems,
Math. Comp., 58 (1992), 161–183.

[46] (MR1620795) [10.1007/s002459900081]
O. Scherzer,
A modified Landweber iteration for solving parameter estimation problems,
Appl. Math. Optim., 38 (1998), 45–68.

[47] (MR991919) [10.1088/0266-5611/5/2/008]
T. I. Seidman and C. R. Vogel,
Well-posedness and convergence of some regularisation methods for non-linear ill posed problems,
Inverse Problems, 5 (1989), 227–238.

[48] [10.1109/TUFFC.2011.1933]
F. Varray, O. Basset, P. Tortoli and C. Cachard,
Extensions of nonlinear B/A parameter imaging methods for echo mode,
IEEE Trans. Ultrasonics, Ferroelectrics, and Frequency Control, 58 (2011), 1232–1244.

[49] [10.1121/1.1918525]
P. J. Westervelt,
Parametric acoustic array,
J. Acoustical Soc. Amer., 35 (1963), 535–537.

[50]
M. Yamamoto and B. Kaltenbacher,
An inverse source problem related to acoustic nonlinearity parameter imaging,
to appear, Time-Dependent Problems in Imaging and Parameter Identification, Springer, 2021.

[51]
E. A. Zabolotskaya and R. V. Khokhlov,
Quasi-plane waves in the non-linear acoustics of confined beams,
Soviet Physics - Acoustics, 15 (1969), 35–40.

[52] [10.1121/1.1344160]
D. Zhang, X. Chen and X.-F. Gong,
Acoustic nonlinearity parameter tomography for biological tissues via parametric array from a circular piston source - Theoretical analysis and computer simulations,
J. Acoustical Soc. Amer., 109 (2001), 1219–1225.

[53] [10.1121/1.415427]
D. Zhang, X. Gong and S. Ye,
Acoustic nonlinearity parameter tomography for biological specimens via measurements of the second harmonics,
J. Acoustical Soc. Amer., 99 (1996), 2397–2402.

Received August 2020; 1st revision August 2020; 2nd revision November 2020.

E-mail address<EMAIL_ADDRESS>
E-mail address<EMAIL_ADDRESS>