# A paraxial approach for the inverse problem of vibro-acoustic imaging in
frequency domain

Teresa Rauscher
Department of Mathematics
University of Klagenfurt
Klagenfurt, Austria
<EMAIL_ADDRESS>

###### Abstract

This paper presents a paraxial modeling approach for vibro-acoustography, a
high-frequency ultrasound imaging technique that makes use of the excited low-
frequency field to achieve a higher resolution while avoiding speckles. We
start from a general second order wave equation, introduce a second order
perturbation and make use of a paraxial change of variables to pass over to
directed beams. Depending on the relation between the second order
perturbation and the directivity of the beams this leads to different systems
of PDEs that involve a spatially variable parameter governing the interaction
of the incoming beams to generate the propagating acoustic field. We then
consider the inverse problem of determining this parameter for one of these
models where we present a Landweber-Kaczmarz reconstruction method to obtain a
spatial image of the region of interest.

_Keywords_ vibro-acoustography $\cdot$ ultrasound imaging $\cdot$ paraxial
modeling $\cdot$ inverse problem $\cdot$ parameter identification

##  1 Introduction

In medical diagnostics the quest for non-invasive, safe and precise imaging
techniques has been relentless. While traditional ultrasound imaging is well-
established and widely used, there is a continuing interest in nonlinear
acoustics as the demand for higher image resolution and diagnostic accuracy
continues to grow. However, nonlinear models also come with drawbacks such as
scattering from small inclusions and stronger attenuation at higher
frequencies. In this paper we consider vibro-acoustic imaging by means of
ultrasound imaging [3, 4]. This technique avoids the drawbacks of nonlinear
models while achieving speckle-free images of high resolution and providing
valuable diagnostic information on tissue properties.

In contrast to conventional ultrasound imaging, where a beam of ultrasound
waves is sent into the tissue using a transducer, in vibro-acoustic imaging
two ultrasound beams at high but slightly different frequencies are sent into
the medium. At a focus they interact nonlinearly and this interaction excites
a wave field that basically propagates at the difference frequency and is
eventually measured by a hydrophone. For more details on the experiments see
[3, 4] and for a graphical illustration of the described experiment for image
acquisition see Figure 1.

Vibro-acoustography can be used in many biomedical and non-destructive
applications, especially, to image several types of tissues, such as breast,
liver, prostate, arteries, and thyroid. Typically, the ultrasound frequencies
used are in the range of $1$ to $10$ MHz and the difference frequencies range
form $10$ to $100$ kHz. For wide applicability in clinical settings one uses
linear-array transducers for forming vibro-acoustography beams. More
information on beamforming and implementation on clinical ultrasound systems
can be found in [18, 19, 20].

A modeling and simulation framework for vibro-acoustic imaging has been
presented in [11, 12]. The model consists of a system of partial differential
equations that involve the speed of sound $c=c(x)$ and the nonlinearity
parameter $\frac{B}{A}=\frac{B}{A}(x)$ which depend on space due to the
inhomogeneity of the medium. Therefore, their reconstruction yields a spatial
image of the region of interest. In case of reconstructing the individual
coefficients $c=c(x)$ or $\frac{B}{A}=\frac{B}{A}(x)$, this is related to
ultrasound tomography and nonlinearity parameter imaging, respectively, cf.
e.g. [1, 2, 7] and the citing literature. The corresponding inverse problem
for vibro-acoustic imaging has been put into a mathematical framework in [9].
In this model all waves are assumed to propagate omni-directional, although
the high frequency waves actually have a strongly preferred direction of
propagation, which can justify the use of a parabolic approximation, cf. e.g.
[17].

There are several methods of deriving parabolic approximations [8]. In [11,
12] a diagonalization procedure has been chosen for the parabolic approach and
for efficient numerical simulation a decomposition approach has been devised.
The forward problem is split into three components: (i) directed high
frequency propagation of the two ingoing beams $\phi_{1}$, $\phi_{2}$, (ii)
nonlinear interaction of these at the focus, and (iii) undirected low
frequency propagation of $\psi$ to the hydrophone.

In this paper the goal is to derive directed models making use of the paraxial
approximation to derive the Khokhlov-Zabolotskaya-Kuznetsov (KZK) equation
[21] from Kuznetsov’s equation. The KZK equation is a directed nonlinear wave
equation that is particularly important for modeling and understanding the
behavior of ultrasound waves in medical ultrasound applications as it takes
various effects including nonlinearity, diffraction and absorption into
account, see also, e.g. [15] for the analysis of the KZK equation. As we pass
from undirected beams to directed beams, the equations change significantly
and an additional small parameter for the order of the direction comes into
play.

Throughout this paper we assume the speed of sound to be constant and only the
nonlinearity parameter to be space dependent. Reconstructing both parameters
simultaneously the paraxial change of variables has to be adapted and
redefined, see equation (6), since this transformation of variables changes
the time variable to the retarded time, which means that the time gets space
dependent and therefore, also dependent on the speed of sound. Results on the
simultaneous reconstruction of sound speed and nonlinearity parameter can be
found in [10].

The plan of this paper is as follows. In Section 2, we first of all derive
several paraxial models starting from the second order wave equation,
introducing a second order perturbation for the velocity potential and then
performing the paraxial change of variables. Then we will transfer the time
domain formulations to frequency domain, equip the system with suitable
boundary conditions and give some remarks on adaptions and generalizations of
the model and its experimental setup. Focusing on one of these models in
Section 3, we will state the inverse problem of nonlinearity imaging and
present a Landweber-Kaczmarz method to solve the linear inverse problem of
reconstructing the nonlinearity parameter. Algorithms for the direct problem
and also the inverse problem can be found in Section 4.

##  2 The model

In vibro-acoustography two ultrasound beams with velocity potentials
$\phi_{1}$, $\phi_{2}$ are excited by a (linear) transducer. They are focused
at two slightly different high frequencies $\omega_{1}$, $\omega_{2}$
($\Delta\omega=\omega_{1}-\omega_{2}\ll\omega_{1},\,\omega_{2}$) with
$\omega_{1}>\omega_{2}>0$. These ingoing beams have a strong preferred
direction of propagation, here on the central axis $z$, away from the
transducer. At a focus they interact nonlinearly and this interaction excites
a wave field with velocity potential $\psi$ that basically propagates at the
difference frequency $\Delta\omega$, which is then quite low. Since the source
of the low frequency field is point-like the resulting field is omni-
directional. The acoustic response is eventually measured by a hydrophone
located far from the focus point and an image can be formed. The experimental
setup for image acquisition is illustrated in Figure 1.

Figure 1: experimental setup

We assume that the beam is formed by a linear transducer. Note that there are
different ways to arrange the two beams. We choose to split the transducer
array in left and right, since we want to keep the beams symmetrically to the
central axis and the focus. This is emphasized in Figure 1 by the frequency
$\omega_{1}$ on the left and $\omega_{2}$ on the right. Some further remarks
on adaptions and generalizations of the experimental setup can be found in
Section 2.4.

###  2.1 Derivation of paraxial models

We now want to look at the derivation of a directive model for the ingoing
beams in this methodology. We assume that the beams are naturally focused, so
by geometrical means and make use of the parabolic approximation for the
derivation of the KZK equation from Kuznetsov’s equation, cf. e.g., [15].

We assume that the direction of propagation is the $x_{1}$-axis. Then, the
acoustical properties of directive beam’s propagation are that the beams

  * •

are concentrated near the $x_{1}$-direction;

  * •

propagate along the $x_{1}$-direction;

  * •

are generated either by an initial condition or by a forcing term on the
boundary $x_{1}=0$.

Furthermore, it is assumed that the variation of beam’s propagation in the
direction $x^{\prime}=(x_{2},x_{3},\ldots,x_{n})$ perpendicular to the
$x_{1}$-axis is much larger than its variation along the $x_{1}$-axis. This
can mathematically be modeled by a paraxial change of variables. For the
initial variables $(t,x_{1},x^{\prime})$, where
$x^{\prime}=(x_{2},\ldots,x_{n})$, we perform the paraxial change of variables

$\displaystyle(\tau,z,y)=\left(t-\frac{x_{1}}{c},\tilde{\varepsilon}x_{1},\,\sqrt{\tilde{\varepsilon}}x^{\prime}\right),$
(1)

where $\tilde{\varepsilon}$ is a small parameter. This change of variables is
illustrated in Figure 2.

$x_{1}$$x^{\prime}$$t$

$\Rightarrow$
$z=\tilde{\varepsilon}x_{1}$$y=\sqrt{\tilde{\varepsilon}}x^{\prime}$$\tau=t-\frac{x_{1}}{c}$

Figure 2: paraxial change of variables

In the paraxial change of variables in equation (1) the speed of sound is
assumed to be constant. If considered it variable the definition (1), in
particular the retarded time variable, would have to be adapted such that

$\displaystyle(\tau,z,y)=\left(t-\check{\mathfrak{t}},\tilde{\varepsilon}x_{1},\sqrt{\tilde{\varepsilon}}x^{\prime}\right),$
(2)

with
$\check{\mathfrak{t}}=\check{\mathfrak{t}}(x_{1},x^{\prime}),\,\mathfrak{t}=\mathfrak{t}(z,y)$,
where

$\frac{1}{c^{2}}=|\nabla_{x}\check{\mathfrak{t}}|^{2}=|\partial_{x_{1}}\check{\mathfrak{t}}|^{2}+|\nabla_{x^{\prime}}\check{\mathfrak{t}}|^{2}=\tilde{\varepsilon}^{2}|\partial_{z}\mathfrak{t}|^{2}+\tilde{\varepsilon}|\nabla_{y}\mathfrak{t}|^{2},\quad\check{\mathfrak{t}}(0,x^{\prime})=0,\
x^{\prime}\in\tilde{\Omega}_{y}.$

In the following, we assume constant speed of sound and consider the
definition for the paraxial change of variables in (1).

To derive a model for vibro-acoustic imaging, we start with the general wave
equation in terms of the velocity potential $\phi$ given by

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

where $\gamma=\frac{B}{A}+1$ and $\phi=\phi(x,t)$ with
$x=(x_{1},\ldots,x_{n})\in\mathbb{R}^{n}$, $t\in\left[0,T\right]$ for $T>0$.
Equation (3) is exact for a perfect lossless gas and often serves as a
starting point for perturbation analyses in nonlinear acoustics.
Approximations lead to simplified wave equations, e.g. the Westervelt
equation, that are more amenable to analysis. For more details and the
derivation see [6].

Since the acoustic state variables are all $O(\varepsilon)$ or smaller, one
can expect that the velocity potential depends analytically on $\varepsilon$.
Therefore, we introduce the following second order approximation for the
velocity potential $\phi$ such that

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

where $\varepsilon$ is very small, such that $\phi_{1}$, $\phi_{2}$ represent
the high-frequency beams and $\psi$ the low-frequency field.

Now, plugging this power series expansion of $\phi$ relative to $\varepsilon$
into the general wave equation given by equation (3) yields

$\displaystyle\varepsilon\partial_{t}^{2}(\phi_{1}+\phi_{2})+\varepsilon^{2}\partial_{t}^{2}\psi$
$\displaystyle-c^{2}\varepsilon\nabla^{2}(\phi_{1}+\phi_{2})-c^{2}\varepsilon^{2}\nabla^{2}\psi$
$\displaystyle=\varepsilon^{2}\partial_{t}|\nabla(\phi_{1}+\phi_{2})|^{2}+\varepsilon^{2}(\gamma-1)\partial_{t}(\phi_{1}+\phi_{2})\nabla^{2}(\phi_{1}+\phi_{2})+O(\varepsilon^{3}).$

Considering only terms up to order $\varepsilon^{2}$ and rearranging the
equation above gives

$\displaystyle\varepsilon\left[\partial_{t}^{2}(\phi_{1}+\phi_{2})-c^{2}\nabla^{2}(\phi_{1}+\phi_{2})-\varepsilon\partial_{t}|\nabla(\phi_{1}+\phi_{2})|^{2}-\right.$
$\displaystyle\left.\varepsilon(\gamma-1)\partial_{t}(\phi_{1}+\phi_{2})\nabla^{2}(\phi_{1}+\phi_{2})\right]$
$\displaystyle=\varepsilon^{2}\left[-\partial_{t}^{2}\psi+c^{2}\nabla^{2}\psi\right].$
(5)

One can observe that we arrived at well known wave equations. The equation on
the left hand side in terms of the velocity potential $\phi_{1}+\phi_{2}$ is
known as the Kuznetsov equation (without viscosity) and the right hand side is
just the linear wave equation in terms of $\psi$.

The next step is to include the concept of directive beams in this model.
Therefore, we perform the paraxial change of variables (1), where
$\tilde{\varepsilon}\ll 1$, so
$\tilde{\varepsilon}\ll\sqrt{\tilde{\varepsilon}}$ and in general
$\tilde{\varepsilon}$ can be different from $\varepsilon$ that was introduced
for the second order approximation in (4) but for the model to make sense we
will see that the relation
$\varepsilon\leq\tilde{\varepsilon}\leq\sqrt{\varepsilon}$ should hold. For
the derivatives this transform of variables yields

$\partial_{t}=\partial_{\tau},\quad\nabla=\left(-\frac{1}{c}\partial_{\tau}+\tilde{\varepsilon}\partial_{z},\sqrt{\tilde{\varepsilon}}\nabla_{y}\right),\quad\nabla^{2}=\frac{1}{c^{2}}\partial_{\tau}^{2}-2\frac{\tilde{\varepsilon}}{c}\partial^{2}_{\tau
z}+\tilde{\varepsilon}^{2}\partial_{z}^{2}+\tilde{\varepsilon}\nabla_{y}^{2}$
(6)

and for the term $|\nabla\phi|^{2}$ one arrives at

$|\nabla\phi|^{2}=\nabla\phi\cdot\nabla\phi=\frac{1}{c^{2}}(\partial_{\tau}\phi)^{2}-2\frac{\tilde{\varepsilon}}{c}\partial_{\tau}\phi\partial_{z}\phi+\tilde{\varepsilon}^{2}(\partial_{z}\phi)^{2}+\tilde{\varepsilon}|\nabla_{y}\phi|^{2}.$

For the left hand side of equation (2.1) this change of variables then leads
to

$\displaystyle\varepsilon$
$\displaystyle\left[\partial_{t}^{2}(\phi_{1}+\phi_{2})-c^{2}\nabla^{2}(\phi_{1}+\phi_{2})-\varepsilon\partial_{t}|\nabla(\phi_{1}+\phi_{2})|^{2}-\varepsilon(\gamma-1)\partial_{t}(\phi_{1}+\phi_{2})\nabla^{2}(\phi_{1}+\phi_{2})\right]$
$\displaystyle=\varepsilon\left[2\tilde{\varepsilon}c\partial^{2}_{\tau
z}(\phi_{1}+\phi_{2})-\tilde{\varepsilon}^{2}c^{2}\partial_{z}^{2}(\phi_{1}+\phi_{2})-\tilde{\varepsilon}c^{2}\nabla_{y}^{2}(\phi_{1}+\phi_{2})\right.$
$\displaystyle\left.-\varepsilon\left(\frac{\gamma+1}{2c^{2}}\partial_{\tau}\left|\partial_{\tau}(\phi_{1}+\phi_{2})\right|^{2}-2\frac{\tilde{\varepsilon}}{c}\left(\partial_{\tau}^{2}(\phi_{1}+\phi_{2})\partial_{z}(\phi_{1}+\phi_{2})+\gamma\partial_{\tau}(\phi_{1}+\phi_{2})\partial_{\tau
z}^{2}(\phi_{1}+\phi_{2})\right)\right.\right.$
$\displaystyle\left.\left.+\tilde{\varepsilon}^{2}\left(\partial_{\tau}(\partial_{z}(\phi_{1}+\phi_{2}))^{2}+(\gamma-1)\partial_{\tau}(\phi_{1}+\phi_{2})\partial_{z}^{2}(\phi_{1}+\phi_{2})\right)\right.\right.$
$\displaystyle\left.\left.+\tilde{\varepsilon}\left(\partial_{\tau}|\nabla_{y}(\phi_{1}+\phi_{2})|^{2}+(\gamma-1)\partial_{\tau}(\phi_{1}+\phi_{2})\nabla_{y}^{2}(\phi_{1}+\phi_{2})\right)\right)\right].$

A paraxial change of variables for the velocity potential $\psi$ for the right
hand side of equation (2.1) yields

$\displaystyle\varepsilon^{2}\left[-\partial_{t}^{2}\psi+c^{2}\nabla^{2}\psi\right]$
$\displaystyle=\varepsilon^{2}\left[-2\tilde{\varepsilon}c\partial_{\tau
z}\psi+c^{2}\tilde{\varepsilon}^{2}\partial_{z}^{2}\psi+\tilde{\varepsilon}c^{2}\nabla_{y}^{2}\psi\right].$

Therefore, altogether we end up with the equation

$\displaystyle\varepsilon$
$\displaystyle\left[2\tilde{\varepsilon}c\partial^{2}_{\tau
z}(\phi_{1}+\phi_{2})-\tilde{\varepsilon}^{2}c^{2}\partial_{z}^{2}(\phi_{1}+\phi_{2})-\tilde{\varepsilon}c^{2}\nabla_{y}^{2}(\phi_{1}+\phi_{2})\right.$
$\displaystyle\left.-\varepsilon\left(\frac{\gamma+1}{2c^{2}}\partial_{\tau}\left|\partial_{\tau}(\phi_{1}+\phi_{2})\right|^{2}-2\frac{\tilde{\varepsilon}}{c}\left(\partial_{\tau}^{2}(\phi_{1}+\phi_{2})\partial_{z}(\phi_{1}+\phi_{2})+\gamma\partial_{\tau}(\phi_{1}+\phi_{2})\partial_{\tau
z}^{2}(\phi_{1}+\phi_{2})\right)\right.\right.$
$\displaystyle\left.\left.+\tilde{\varepsilon}^{2}\left(\partial_{\tau}(\partial_{z}(\phi_{1}+\phi_{2}))^{2}+(\gamma-1)\partial_{\tau}(\phi_{1}+\phi_{2})\partial_{z}^{2}(\phi_{1}+\phi_{2})\right)\right.\right.$
$\displaystyle\left.\left.+\tilde{\varepsilon}\left(\partial_{\tau}|\nabla_{y}(\phi_{1}+\phi_{2})|^{2}+(\gamma-1)\partial_{\tau}(\phi_{1}+\phi_{2})\nabla_{y}^{2}(\phi_{1}+\phi_{2})\right)\right)\right]$
$\displaystyle=$
$\displaystyle\,\varepsilon^{2}\left[-2\tilde{\varepsilon}c\partial_{\tau
z}\psi+c^{2}\tilde{\varepsilon}^{2}\partial_{z}^{2}\psi+\tilde{\varepsilon}c^{2}\nabla_{y}^{2}\psi\right].$

Reordering with respect to the different orders of $\varepsilon$ and
$\tilde{\varepsilon}$ yields

$\displaystyle\varepsilon\tilde{\varepsilon}\left[2c\partial^{2}_{\tau
z}(\phi_{1}+\phi_{2})-c^{2}\nabla_{y}^{2}(\phi_{1}+\phi_{2})\right]$
$\displaystyle=0,$ (7)
$\displaystyle\varepsilon^{2}\left[-\frac{\gamma+1}{2c^{2}}\partial_{\tau}\left|\partial_{\tau}(\phi_{1}+\phi_{2})\right|^{2}\right]$
$\displaystyle=0,$ (8)
$\displaystyle\varepsilon\tilde{\varepsilon}^{2}\left[-c^{2}\partial_{z}^{2}(\phi_{1}+\phi_{2})\right]$
$\displaystyle=0,$ (9)
$\displaystyle\varepsilon^{2}\tilde{\varepsilon}\left[\frac{2}{c}\left(\partial_{\tau}^{2}(\phi_{1}+\phi_{2})\partial_{z}(\phi_{1}+\phi_{2})+\gamma\partial_{\tau}(\phi_{1}+\phi_{2})\partial_{\tau
z}^{2}(\phi_{1}+\phi_{2})\right)-\partial_{\tau}|\nabla_{y}(\phi_{1}+\phi_{2})|^{2}\right.$
$\displaystyle\left.-(\gamma-1)\partial_{\tau}(\phi_{1}+\phi_{2})\nabla_{y}^{2}(\phi_{1}+\phi_{2})+2c\partial^{2}_{\tau
z}\psi-c^{2}\nabla_{y}^{2}\psi\right]$ $\displaystyle=0,$ (10)
$\displaystyle\varepsilon^{2}\tilde{\varepsilon}^{2}\left[-\partial_{\tau}\left(\partial_{z}(\phi_{1}+\phi_{2})\right)^{2}-(\gamma-1)\partial_{\tau}(\phi_{1}+\phi_{2})\partial_{z}^{2}(\phi_{1}+\phi_{2})-c^{2}\partial_{z}^{2}\psi\right]$
$\displaystyle=0.$ (11)

The order of the system of equations should stay as they are stated above.
Therefore, the relation
$\varepsilon\leq\tilde{\varepsilon}\leq\sqrt{\varepsilon}$ should hold such
that we arrive at different cases. The following table gives an overview on
these.

Case |  | equ. (7) | equ. (8) | equ. (9) | equ. (10) | equ. (11)
---|---|---|---|---|---|---
Case 1: | $\tilde{\varepsilon}=\varepsilon^{1/2}$ | $\varepsilon^{3/2}$ | $\varepsilon^{2}$ | $\varepsilon^{2}$ | $\varepsilon^{5/2}$ | $\varepsilon^{3}$
Case 2: | $\varepsilon<\tilde{\varepsilon}<\varepsilon^{1/2}$ | $\varepsilon\tilde{\varepsilon}$ | $\varepsilon^{2}$ | $\varepsilon\tilde{\varepsilon}^{2}$ | $\varepsilon^{2}\tilde{\varepsilon}$ | $\varepsilon^{2}\tilde{\varepsilon}^{2}$
Case 3: | $\tilde{\varepsilon}=\varepsilon$ | $\varepsilon^{2}$ | $\varepsilon^{2}$ | $\varepsilon^{3}$ | $\varepsilon^{3}$ | $\varepsilon^{4}$
Table 1: Relation between $\varepsilon$ and $\tilde{\varepsilon}$ for
equations (7) - (11)

Note that the waves with velocity potentials $\phi_{1}$ and $\phi_{2}$ follow
a parabolic approximation whereas the resulting field is omni-directional and
the hydrophone is typically located far from the focal point, which means that
the parabolic approximation is not ideal for the resulting wave that is here
expressed in terms of $\psi$. So, since $\phi_{1}$, $\phi_{2}$ are paraxial
$\tilde{\varepsilon}$ is quite small from the point of view of $\phi_{1}$,
$\phi_{2}$ whereas for $\psi$ it holds that $\tilde{\varepsilon}$ is not that
small from its point of view.

The following table gives an overview of the different models that we end up
with. The subcases refer to a different accuracy or respectively to different
orders of $\varepsilon$ up to which terms are considered. The enumeration (i)
and (ii) refers to the first and second equation in the following case
distinctions, respectively. The colors emphasize that these cases are very
similar and just differ in some constants.

Case |  | equ. (7) | equ. (8) | equ. (9) | equ. (10) | equ. (11)
---|---|---|---|---|---|---
|  |  |  |  | $\psi$ | $\phi_{1}+\phi_{2}$ | $\psi$ | $\phi_{1}+\phi_{2}$
Case 1: A | $\tilde{\varepsilon}=\varepsilon^{1/2}$ | (i) | (ii) | (ii) | (ii) |  | (ii) |
Case 1: B | $\tilde{\varepsilon}=\varepsilon^{1/2}$ | (i) | (ii) | (ii) | (ii) | (ii) | (ii) |
Case 2: A | $\varepsilon<\tilde{\varepsilon}<\varepsilon^{1/2}$ | (i) | (ii) |  | (ii) |  | (ii) |
Case 2: B | $\varepsilon<\tilde{\varepsilon}<\varepsilon^{1/2}$ | (i) | (ii) | (ii) | (ii) |  | (ii) |
Case 2: C | $\varepsilon<\tilde{\varepsilon}<\varepsilon^{1/2}$ | (i) | (ii) | (ii) | (ii) | (ii) | (ii) |
Case 3: | $\tilde{\varepsilon}=\varepsilon$ | (i) | (i) | (ii) | (ii) | (ii) | (ii) |
Table 2: Overview of the following cases on the relation of $\varepsilon$ and
$\tilde{\varepsilon}$ for equations (7) - (11)

Case 1: A For the case $\tilde{\varepsilon}=\sqrt{\varepsilon}$ considering
only terms up to order $\varepsilon^{2}$ we end up with the model

$\displaystyle\frac{2}{c}\partial^{2}_{\tau z}\phi_{k}-\nabla_{y}^{2}\phi_{k}$
$\displaystyle=0,\qquad\text{ for }k\in\left\\{1,2\right\\}$ (12)
$\displaystyle\tilde{\varepsilon}\left[\frac{2}{c}\partial^{2}_{\tau
z}\psi-\nabla_{y}^{2}\psi\right]-\tilde{\varepsilon}^{2}\partial_{z}^{2}\psi$
$\displaystyle=\frac{\gamma+1}{2c^{4}}\partial_{\tau}\left|\partial_{\tau}(\phi_{1}+\phi_{2})\right|^{2}+\partial_{z}^{2}(\phi_{1}+\phi_{2}),$
(13)

where we split the first equation in $\phi_{1}+\phi_{2}$ into two separate
equations for $\phi_{1}$ and $\phi_{2}$ since they are linear and can be
considered independently. Observe that the way equation (13) is written
indicates the solution approach, in that equation (12) is solved for
$\phi_{1}$, $\phi_{2}$, which then build a source term driving the second
equation.

Case 1: B For the case $\tilde{\varepsilon}=\sqrt{\varepsilon}$ considering
terms up to order $\varepsilon^{5/2}$ we end up with the model

$\displaystyle\frac{2}{c}\partial^{2}_{\tau z}\phi_{k}-\nabla_{y}^{2}\phi_{k}$
$\displaystyle=0,\qquad\text{ for }k\in\left\\{1,2\right\\}$ (14)
$\displaystyle\tilde{\varepsilon}\left[\frac{2}{c}\partial^{2}_{\tau
z}\psi-\nabla_{y}^{2}\psi\right]-\tilde{\varepsilon}^{2}\partial_{z}^{2}\psi$
$\displaystyle=\frac{\gamma+1}{2c^{4}}\partial_{\tau}\left|\partial_{\tau}(\phi_{1}+\phi_{2})\right|^{2}+\partial_{z}^{2}(\phi_{1}+\phi_{2})$
$\displaystyle-\frac{\tilde{\varepsilon}}{c^{2}}\left[\frac{2}{c}\left(\partial_{\tau}^{2}(\phi_{1}+\phi_{2})\partial_{z}(\phi_{1}+\phi_{2})+\gamma\partial_{\tau}(\phi_{1}+\phi_{2})\partial_{\tau
z}^{2}(\phi_{1}+\phi_{2})\right)\right.$
$\displaystyle\left.-\partial_{\tau}|\nabla_{y}(\phi_{1}+\phi_{2})|^{2}-(\gamma-1)\partial_{\tau}(\phi_{1}+\phi_{2})\nabla_{y}^{2}(\phi_{1}+\phi_{2})\right].$
(15)

Using the first equation and plugging

$\nabla_{y}^{2}\phi_{k}=\frac{2}{c}\partial^{2}_{\tau z}\phi_{k}$

into the second equation, it reduces to

$\displaystyle\tilde{\varepsilon}\left[\frac{2}{c}\partial^{2}_{\tau
z}\psi-\nabla_{y}^{2}\psi\right]-\tilde{\varepsilon}^{2}\partial_{z}^{2}\psi$
$\displaystyle=\frac{\gamma+1}{2c^{4}}\partial_{\tau}\left|\partial_{\tau}(\phi_{1}+\phi_{2})\right|^{2}+\partial_{z}^{2}(\phi_{1}+\phi_{2})$
$\displaystyle-\frac{\tilde{\varepsilon}}{c^{2}}\left[\frac{2}{c}\left(\partial_{\tau}^{2}(\phi_{1}+\phi_{2})\partial_{z}(\phi_{1}+\phi_{2})+\partial_{\tau}(\phi_{1}+\phi_{2})\partial_{\tau
z}^{2}(\phi_{1}+\phi_{2})\right)\right.$
$\displaystyle\left.-\partial_{\tau}|\nabla_{y}(\phi_{1}+\phi_{2})|^{2}\right].$
(16)

Case 2: A In case $\varepsilon<\tilde{\varepsilon}<\varepsilon^{1/2}$ for
terms up to order $\varepsilon^{2}$ the system of equations is

$\displaystyle\frac{2}{c}\partial^{2}_{\tau z}\phi_{k}-\nabla_{y}^{2}\phi_{k}$
$\displaystyle=0,\qquad\text{ for }k\in\left\\{1,2\right\\}$ (17)
$\displaystyle\tilde{\varepsilon}\left[\frac{2}{c}\partial^{2}_{\tau
z}\psi-\nabla_{y}^{2}\psi\right]-\tilde{\varepsilon}^{2}\partial_{z}^{2}\psi$
$\displaystyle=\frac{\gamma+1}{2c^{4}}\partial_{\tau}\left|\partial_{\tau}(\phi_{1}+\phi_{2})\right|^{2}.$
(18)

Note, that the model contains the KZK equation

$\frac{2}{c}\partial^{2}_{\tau
z}\phi_{k}-\nabla_{y}^{2}\phi_{k}=\frac{\gamma+1}{2c^{4}}\partial_{\tau}\left|\partial_{\tau}(\phi_{1}+\phi_{2})\right|^{2},$

but in our case the quadratic nonlinearity is decoupled and appears as a
source term for the equation in terms of $\psi$.

Case 2: B In case $\varepsilon<\tilde{\varepsilon}<\varepsilon^{1/2}$ for
terms up to order $\varepsilon\tilde{\varepsilon}^{2}$ the system of equations
is

$\displaystyle\frac{2}{c}\partial^{2}_{\tau z}\phi_{k}-\nabla_{y}^{2}\phi_{k}$
$\displaystyle=0,\qquad\text{ for }k\in\left\\{1,2\right\\}$ (19)
$\displaystyle\tilde{\varepsilon}\left[\frac{2}{c}\partial^{2}_{\tau
z}\psi-\nabla_{y}^{2}\psi\right]-\tilde{\varepsilon}^{2}\partial_{z}^{2}\psi$
$\displaystyle=\frac{\gamma+1}{2c^{4}}\partial_{\tau}\left|\partial_{\tau}(\phi_{1}+\phi_{2})\right|^{2}+\frac{\tilde{\varepsilon}^{2}}{\varepsilon}\partial_{z}^{2}(\phi_{1}+\phi_{2}).$
(20)

Note that since $\varepsilon<\tilde{\varepsilon}<\varepsilon^{1/2}$ it holds
that

$\varepsilon^{3}<\varepsilon^{2}\tilde{\varepsilon}<\varepsilon\tilde{\varepsilon}^{2}<\varepsilon^{2}$

and therefore,

$\varepsilon<\tilde{\varepsilon}<\frac{\tilde{\varepsilon}^{2}}{\varepsilon}<1.$

Case 2: C In case $\varepsilon<\tilde{\varepsilon}<\varepsilon^{1/2}$ for
terms up to order $\varepsilon^{2}\tilde{\varepsilon}$ the system of equations
is

$\displaystyle\frac{2}{c}\partial^{2}_{\tau z}\phi_{k}-\nabla_{y}^{2}\phi_{k}$
$\displaystyle=0,\qquad\text{ for }k\in\left\\{1,2\right\\}$ (21)
$\displaystyle\tilde{\varepsilon}\left[\frac{2}{c}\partial^{2}_{\tau
z}\psi-\nabla_{y}^{2}\psi\right]-\tilde{\varepsilon}^{2}\partial_{z}^{2}\psi$
$\displaystyle=\frac{\gamma+1}{2c^{4}}\partial_{\tau}\left|\partial_{\tau}(\phi_{1}+\phi_{2})\right|^{2}+\frac{\tilde{\varepsilon}^{2}}{\varepsilon}\partial_{z}^{2}(\phi_{1}+\phi_{2})$
$\displaystyle-\frac{\tilde{\varepsilon}}{c^{2}}\left[\frac{2}{c}\left(\partial_{\tau}^{2}(\phi_{1}+\phi_{2})\partial_{z}(\phi_{1}+\phi_{2})+\partial_{\tau}(\phi_{1}+\phi_{2})\partial_{\tau
z}^{2}(\phi_{1}+\phi_{2})\right)\right.$
$\displaystyle\left.-\partial_{\tau}|\nabla_{y}(\phi_{1}+\phi_{2})|^{2}\right],$
(22)

where we reduced the second equation by plugging the first one into the
second.

Case 3: For the case $\tilde{\varepsilon}=\varepsilon$ considering terms up to
order $\varepsilon^{3}$ we end up with the model

$\displaystyle\frac{2}{c}\partial^{2}_{\tau
z}(\phi_{1}+\phi_{2})-\nabla_{y}^{2}(\phi_{1}+\phi_{2})$
$\displaystyle=\frac{\gamma+1}{2c^{4}}\partial_{\tau}\left|\partial_{\tau}(\phi_{1}+\phi_{2})\right|^{2},$
(23) $\displaystyle\tilde{\varepsilon}\left[\frac{2}{c}\partial^{2}_{\tau
z}\psi-\nabla_{y}^{2}\psi\right]-\tilde{\varepsilon}^{2}\partial_{z}^{2}\psi$
$\displaystyle=\partial_{z}^{2}(\phi_{1}+\phi_{2})$
$\displaystyle-\frac{\tilde{\varepsilon}}{c^{2}}\left[\frac{2}{c}\left(\partial_{\tau}^{2}(\phi_{1}+\phi_{2})\partial_{z}(\phi_{1}+\phi_{2})+\partial_{\tau}(\phi_{1}+\phi_{2})\partial_{\tau
z}^{2}(\phi_{1}+\phi_{2})\right)\right.$
$\displaystyle\left.-\partial_{\tau}|\nabla_{y}(\phi_{1}+\phi_{2})|^{2}\right],$
(24)

where we reduced the second equation by plugging the first one into the
second. In this case the first equation in $\phi_{1}+\phi_{2}$ can not be
split since it is nonlinear in $\phi_{1}+\phi_{2}$.

###  2.2 Transformation into frequency domain

The next step is to transfer the time domain formulations into frequency
domain with the time harmonic ansatz

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

We will therefore look at the cases from above and deal with them separately.
Note that we assume $\omega_{1}>\omega_{2}>0$.

Case 1: A For equation (12) the time harmonic ansatz yields

$\displaystyle\frac{2}{c}\partial^{2}_{\tau
z}\phi_{k}-\nabla_{y}^{2}\phi_{k}=\frac{2}{c}\partial^{2}_{\tau
z}\hat{\phi}_{k}e^{\imath\omega_{k}\tau}-\nabla_{y}^{2}\hat{\phi}_{k}e^{\imath\omega_{k}\tau}=\left(\frac{2\imath\omega_{k}}{c}\partial_{z}\hat{\phi}_{k}-\nabla_{y}^{2}\hat{\phi}_{k}\right)e^{\imath\omega_{k}\tau}=0.$

Therefore,

$\frac{2\imath\omega_{k}}{c}\partial_{z}\hat{\phi}_{k}-\nabla_{y}^{2}\hat{\phi}_{k}=0\qquad\text{
for }k\in\left\\{1,2\right\\}.$ (25)

Now, going on to the second equation. Starting with the right hand side of
equation (13) the time harmonic ansatz yields

$\displaystyle\frac{\gamma+1}{2c^{4}}\partial_{\tau}\left|\partial_{\tau}(\phi_{1}+\phi_{2})\right|^{2}$
$\displaystyle+\partial_{z}^{2}(\phi_{1}+\phi_{2})$ $\displaystyle=$
$\displaystyle\frac{\gamma+1}{2c^{4}}\partial_{\tau}\left|\imath\omega_{1}\hat{\phi}_{1}e^{\imath\omega_{1}\tau}+\imath\omega_{2}\hat{\phi}_{2}e^{\imath\omega_{2}\tau}\right|^{2}+e^{\imath\omega_{1}\tau}\partial_{z}^{2}\hat{\phi}_{1}+e^{\imath\omega_{2}\tau}\partial_{z}^{2}\hat{\phi}_{2}$
$\displaystyle=$
$\displaystyle\frac{\gamma+1}{2c^{4}}\left[2(\omega_{1}-\omega_{2})\omega_{1}\omega_{2}\Re{\imath\hat{\phi}_{1}\overline{\hat{\phi}}_{2}e^{\imath(\omega_{1}-\omega_{2})\tau}}\right]+e^{\imath\omega_{1}\tau}\partial_{z}^{2}\hat{\phi}_{1}+e^{\imath\omega_{2}\tau}\partial_{z}^{2}\hat{\phi}_{2}.$

For the left hand side of equation (13) the time harmonic ansatz yields

$\displaystyle\tilde{\varepsilon}\left[\frac{2}{c}\partial^{2}_{\tau
z}\psi-\nabla_{y}^{2}\psi\right]-\tilde{\varepsilon}^{2}\partial_{z}^{2}\psi=\Re{\left[\tilde{\varepsilon}\left(\frac{2\imath(\omega_{1}-\omega_{2})}{c}\partial_{z}\hat{\psi}-\nabla_{y}^{2}\hat{\psi}-\tilde{\varepsilon}\partial_{z}^{2}\hat{\psi}\right)\right]e^{\imath(\omega_{1}-\omega_{2})\tau}}.$

So, putting the second equation together yields

$\displaystyle\real\left\\{\left(\frac{2\imath(\omega_{1}-\omega_{2})}{c}\partial_{z}\hat{\psi}-\nabla_{y}^{2}\hat{\psi}-\tilde{\varepsilon}\partial_{z}^{2}\hat{\psi}\right.\right.-$
$\displaystyle\left.\left.\frac{\gamma+1}{2\tilde{\varepsilon}c^{4}}\left[2(\omega_{1}-\omega_{2})\omega_{1}\omega_{2}\imath\hat{\phi}_{1}\overline{\hat{\phi}}_{2}\right]\right)e^{\imath(\omega_{1}-\omega_{2})\tau}\right\\}$
$\displaystyle=\frac{1}{\tilde{\varepsilon}}\left(e^{\imath\omega_{1}\tau}\partial_{z}^{2}\hat{\phi}_{1}+e^{\imath\omega_{2}\tau}\partial_{z}^{2}\hat{\phi}_{2}\right).$

Case 1: B The first equation in frequency domain is already given in (25).
Analogously to the case above, the time harmonic ansatz for the second
equation yields

$\displaystyle\real\left\\{\left(\frac{2\imath(\omega_{1}-\omega_{2})}{c}\partial_{z}\hat{\psi}-\nabla_{y}^{2}\hat{\psi}-\tilde{\varepsilon}\partial_{z}^{2}\hat{\psi}-(\omega_{1}-\omega_{2})\left[\frac{\gamma+1}{\tilde{\varepsilon}c^{4}}\omega_{1}\omega_{2}\imath\hat{\phi}_{1}\overline{\hat{\phi}}_{2}+\frac{2}{c^{2}}\nabla_{y}\hat{\phi}_{1}\nabla_{y}\overline{\hat{\phi}}_{2}\right]\right)e^{\imath(\omega_{1}-\omega_{2})\tau}\right\\}$
$\displaystyle=\frac{1}{\tilde{\varepsilon}}\left(e^{\imath\omega_{1}\tau}\partial_{z}^{2}\hat{\phi}_{1}+e^{\imath\omega_{2}\tau}\partial_{z}^{2}\hat{\phi}_{2}\right)$
$\displaystyle-\frac{2}{c^{3}}\left[-\omega_{1}^{2}\partial_{z}\hat{\phi}_{1}^{2}e^{2\imath\omega_{1}\tau}-\left(\left(\omega_{1}^{2}+\omega_{1}\omega_{2}\right)\hat{\phi}_{1}\partial_{z}\hat{\phi}_{2}+\left(\omega_{2}^{2}+\omega_{1}\omega_{2}\right)\partial_{z}\hat{\phi}_{1}\hat{\phi}_{2}\right)e^{\imath\left(\omega_{1}+\omega_{2}\right)\tau}-\omega_{2}^{2}\partial_{z}\hat{\phi}_{2}^{2}e^{2\imath\omega_{2}\tau}\right].$

Case 2: A Using the derivations in case 1, the model is here given by

$\displaystyle\frac{2\imath\omega_{k}}{c}\partial_{z}\hat{\phi}_{k}-\nabla_{y}^{2}\hat{\phi}_{k}$
$\displaystyle=0\qquad\text{ for }k\in\left\\{1,2\right\\},$ (26)
$\displaystyle\frac{2\imath(\omega_{1}-\omega_{2})}{c}\partial_{z}\hat{\psi}-\nabla_{y}^{2}\hat{\psi}-\tilde{\varepsilon}\partial_{z}^{2}\hat{\psi}$
$\displaystyle=\frac{\gamma+1}{\tilde{\varepsilon}c^{4}}\left[(\omega_{1}-\omega_{2})\omega_{1}\omega_{2}\imath\hat{\phi}_{1}\overline{\hat{\phi}}_{2}\right].$
(27)

One can really observe here that the low frequency field propagates at the
difference frequency $\omega_{1}-\omega_{2}$.

Case 2: B The system of equations only differs in a constant factor from case
1: A, so, we refer to the time harmonic model there.

Case 2: C The system of equations only differs in a constant factor from case
1: B. Therefore, we refer to case 1: B for the time harmonic model.

Case 3: The time harmonic Ansatz yields for the first equation

$\displaystyle\left(\frac{2\imath\omega_{1}}{c}\partial_{z}\hat{\phi}_{1}-\nabla_{y}^{2}\hat{\phi}_{1}\right)e^{\imath\omega_{1}\tau}+\left(\frac{2\imath\omega_{2}}{c}\partial_{z}\hat{\phi}_{2}-\nabla_{y}^{2}\hat{\phi}_{2}\right)e^{\imath\omega_{2}\tau}=\frac{\gamma+1}{c^{4}}\left[(\omega_{1}-\omega_{2})\omega_{1}\omega_{2}\Re{\imath\hat{\phi}_{1}\overline{\hat{\phi}}_{2}e^{\imath(\omega_{1}-\omega_{2})\tau}}\right]$

and for the second

$\displaystyle\real\left\\{\left(\frac{2\imath(\omega_{1}-\omega_{2})}{c}\partial_{z}\hat{\psi}-\nabla_{y}^{2}\hat{\psi}-\tilde{\varepsilon}\partial_{z}^{2}\hat{\psi}-\frac{2(\omega_{1}-\omega_{2})}{c^{2}}\nabla_{y}\hat{\phi}_{1}\nabla_{y}\overline{\hat{\phi}}_{2}\right)e^{\imath(\omega_{1}-\omega_{2})\tau}\right\\}$
$\displaystyle=\frac{1}{\tilde{\varepsilon}}\left(e^{\imath\omega_{1}\tau}\partial_{z}^{2}\hat{\phi}_{1}+e^{\imath\omega_{2}\tau}\partial_{z}^{2}\hat{\phi}_{2}\right)$
$\displaystyle-\frac{2}{c^{3}}\left[-\omega_{1}^{2}\partial_{z}\hat{\phi}_{1}^{2}e^{2\imath\omega_{1}\tau}-\left(\left(\omega_{1}^{2}+\omega_{1}\omega_{2}\right)\hat{\phi}_{1}\partial_{z}\hat{\phi}_{2}+\left(\omega_{2}^{2}+\omega_{1}\omega_{2}\right)\partial_{z}\hat{\phi}_{1}\hat{\phi}_{2}\right)e^{\imath\left(\omega_{1}+\omega_{2}\right)\tau}-\omega_{2}^{2}\partial_{z}\hat{\phi}_{2}^{2}e^{2\imath\omega_{2}\tau}\right].$

###  2.3 Boundary conditions

We consider a bounded computational domain $\tilde{\Omega}$ and due to the
paraxial approach, where we assumed that the beams propagate in $z$-direction,
we write $\tilde{\Omega}$ as a product

$\tilde{\Omega}=\left(0,L\right)\times\tilde{\Omega}_{y}$

with $L>0$ and $\tilde{\Omega}_{y}$ a smooth subset of $\mathbb{R}^{d-1}$, so
that $\left(z,y\right)\in\left(0,L\right)\times\tilde{\Omega}_{y}$.

Since the ultrasound beams are excited by a transducer, with
$h_{k}(\tau,0,y)=\hat{h}_{k}(y)e^{i\omega_{k}\tau}$ being the given time
harmonic initial data at $z=0$ and absorbing boundary conditions on the
boundary parts parallel to the $z$-axis we arrive at

$\displaystyle\hat{\phi}_{k}(0,y)$
$\displaystyle=\hat{h}_{k}(y)\qquad\qquad\quad\text{ in }\tilde{\Omega}_{y},$
$\displaystyle\partial_{\nu_{y}}\hat{\phi}_{k}(z,y)$
$\displaystyle=-\imath\sigma_{k}\hat{\phi}_{k}(z,y)\qquad(z,y)\in\left(0,L\right)\times\partial\tilde{\Omega}_{y}.$
(28)

Also for the velocity potential $\psi$ we assume impedance boundary conditions
in order to avoid spurious reflections of waves on the boundary of the
computational domain, i.e.

$\displaystyle-\partial_{z}\hat{\psi}(0,y)$
$\displaystyle=-\imath\sigma_{0}\hat{\psi}(0,y)\qquad\,\,\,y\in\tilde{\Omega}$
$\displaystyle\partial_{z}\hat{\psi}(L,y)$
$\displaystyle=-\imath\sigma_{L}\hat{\psi}(L,y)\qquad\,\,y\in\tilde{\Omega},$
(29) $\displaystyle\partial_{\nu_{y}}\hat{\psi}(z,y)$
$\displaystyle=-\imath\sigma\hat{\psi}(z,y)\qquad\quad(z,y)\in\left(0,L\right)\times\partial\tilde{\Omega}_{y}.$

Here $\sigma_{0},\,\sigma_{L}>0$ and the coefficients $\sigma_{k}$ and
$\sigma$ are nonnegative $L^{\infty}$ impedance coefficients that are bounded
away from zero on an open subset of $\partial\tilde{\Omega}$.

###  2.4 Adaptions and generalizations of the model

#### Transducer setup.

There are several ways to arrange the two beams. In the considered model the
transducer is split symmetrically in the middle, so along the central axis
where the focus lies, such that the left and right side can be excited with
different frequencies. Another possible arrangement could be to split the
transducer symmetrically into three parts, so a broader part in the middle
where one beam is excited and smaller parts on the left and right side for the
other beams. It could also be interesting to use a random approach for the
excitement of the transducer elements. In general it is as well possible to
use two probes for the two different beams, but this involves further
challenges on the application side which is the reason why medical tests are
currently performed with only one probe.

#### Focusing.

A very important property of the beams in vibro-acoustography is their
focusing. In the different models derived above we assume that the beams are
naturally focused, so by geometrical means, which is assumed to be
incorporated in the initial data $h_{k}$. Another way to achieve focusing is
by adding time delays to achieve a sequence of delayed pulses which produce
the equivalent of a lens, cf. e.g. [16] . This is often referred to as
electronic focusing. We briefly discuss the questions how electronic focusing
can be accomplished and how it can be included in a paraxial model in vibro-
acoustography.

Figure 3 illustrates the geometry for the derivation of the time delay to
focus an element of the transducer. In this case $x_{1}$ is the axis of
propagation and the transducer is aligned to the $x_{2}$-direction.
Furthermore, $x_{2_{n}}$ is the distance from the origin to the center of the
$n$-th element of the transducer and at position $(x_{1_{r}},x_{2_{r}})$ is
the focus point.

Figure 3: Electronic focusing of transducer elements

The time delay to focus an element $x_{2_{n}}$ can be given by

$\tilde{\tau}_{n}=\frac{1}{c}\left[\sqrt{x_{2_{r}}^{2}+x_{1_{r}}^{2}}-\sqrt{(x_{2_{r}}-x_{2_{n}})^{2}+x_{1_{r}}^{2}}\right].$

Under the assumptions that lateral variations are smaller than the axial
distance, i.e. $x_{2_{r}}-x_{2_{n}}\ll x_{1_{r}}$ and the source is only
excited at points $x_{2_{r}}\ll x_{1_{r}}$ the expression for
$\tilde{\tau}_{n}$ can be approximated by

$\tilde{\tau}_{n}\approx\frac{1}{c}\left[\frac{x_{2_{r}}^{2}}{2x_{1_{r}}}-\frac{(x_{2_{r}}-x_{2_{n}})^{2}}{2x_{1_{r}}}\right]=\frac{1}{c}\left[-\frac{x_{2_{n}}^{2}}{2x_{1_{r}}}+x_{2_{n}}\frac{x_{2_{r}}}{x_{1_{r}}}\right]=\frac{1}{c}\left[-\frac{x_{2_{n}}^{2}}{2x_{1_{r}}}+x_{2_{n}}\tan(\theta)\right].$

Since $\cos(\theta)=\frac{x_{2_{r}}}{x_{1_{r}}}\approx 1$ and
$\tan(\theta)\approx\sin(\theta)$, it follows that

$\tilde{\tau}_{n}\approx\frac{1}{c}\left[-\frac{x_{2_{n}}^{2}}{2x_{1_{r}}}+x_{2_{n}}\sin(\theta)\right].$
(30)

For a focused beam at least three points are necessary. Figure 4 gives an
overview of a model where we consider three velocity potentials $\phi_{1}$,
$\phi_{2}$ and $\phi_{3}$. Here, $\phi_{1}$ starts at position $x_{2_{r}}$ and
to keep it more general we assume that the focus point and $x_{2_{r}}$ are not
on the central axis $x_{1}$ and we also do not assume symmetry. The other two
beams start at position $x_{2_{n}}$ and $x_{2_{m}}$ respectively.

Figure 4: Focusing of three beams with time delays

The next step is to include these time delays. Since $\phi_{1}$ is at the
center in our model no additional time delay is needed. Therefore, we perform
the paraxial change of variables

$\phi_{1}(\tau,z,x)=\phi_{1}\left(t-\frac{x_{1}}{c},\tilde{\varepsilon}x_{1},\sqrt{\tilde{\varepsilon}}x_{2}\right).$
(31)

The velocity potential $\phi_{2}$ starts at position $x_{2_{n}}$, which gives
a time delay of $\tilde{\tau}_{n}$ and the direction of propagation is
$\sqrt{(x_{2}-x_{2_{n}})^{2}+x_{1}^{2}}$. Similar to what we have seen above
when we derived the approximated time delay $\tilde{\tau}_{n}$, also the
direction of propagation can be approximated, which yields

$\sqrt{(x_{2}-x_{2_{n}})^{2}+x_{1}^{2}}\approx
x_{1}+\frac{(x_{2}-x_{2_{n}})^{2}}{2x_{1}}.$

Therefore, we perform the following paraxial approach

$\phi_{2}(\tau_{n},z,x)=\phi_{2}\left(t-\frac{x_{1}}{c}-\frac{(x_{2}-x_{2_{n}})^{2}}{2cx_{1}}-\tilde{\tau}_{n},\tilde{\varepsilon}x_{1},\sqrt{\tilde{\varepsilon}}x_{2}\right)=\phi_{2}\left(\tau-\tau_{r_{n}}-\tilde{\tau}_{n},\tilde{\varepsilon}x_{1},\sqrt{\tilde{\varepsilon}}x_{2}\right),$

where $\tau_{r_{n}}=\frac{(x_{2}-x_{2_{n}})^{2}}{2cx_{1}}$ and
$\tilde{\tau}_{n}$ from (30). Similarly for the velocity potential $\phi_{3}$
we get

$\phi_{3}(\tau_{m},z,x)=\phi_{3}\left(t-\frac{x_{1}}{c}-\frac{(x_{2}-x_{2_{m}})^{2}}{2cx_{1}}-\tilde{\tau}_{m},\tilde{\varepsilon}x_{1},\sqrt{\tilde{\varepsilon}}x_{2}\right)=\phi_{3}\left(\tau-\tau_{r_{m}}-\tilde{\tau}_{m},\tilde{\varepsilon}x_{1},\sqrt{\tilde{\varepsilon}}x_{2}\right).$

Now, analogously to above one could derive a model for vibro-acoustic imaging
with electronic focusing.

#### Multiple beams.

Consider for example the model in Case 2: A given by equation (17), (18). This
model can simply be extended to an arbitrary number of beams that are excited
with two different frequencies $\omega_{1},\,\omega_{2}$. Suppose that $n$
beams are excited with frequency $\omega_{1}$ and $m$ beams with $\omega_{2}$.
The model in time domain is

$\displaystyle 2\partial^{2}_{\tau z}\phi_{j}-c\Delta_{y}\phi_{j}$
$\displaystyle=0,$ (32)
$\displaystyle\tilde{\varepsilon}\left[\frac{2}{c}\partial^{2}_{\tau
z}\psi-\nabla_{y}^{2}\psi\right]-\tilde{\varepsilon}^{2}\partial_{z}^{2}\psi$
$\displaystyle=\frac{\gamma+1}{2c^{4}}\partial_{\tau}\left|\partial_{\tau}\sum_{i=1}^{n+m}\phi_{i}\right|^{2},$
(33)

for $j\in\left\\{1,\ldots,n+m\right\\}$. In frequency domain with a time
harmonic ansatz this yields

$\displaystyle\imath\omega_{1}\partial_{z}\hat{\phi}_{k}-c\Delta_{y}\hat{\phi}_{k}$
$\displaystyle=0,$ (34)
$\displaystyle\imath\omega_{2}\partial_{z}\hat{\phi}_{l}-c\Delta_{y}\hat{\phi}_{l}$
$\displaystyle=0,$ (35)
$\displaystyle\frac{2\imath(\omega_{1}-\omega_{2})}{c}\partial_{z}\hat{\psi}-\nabla_{y}^{2}\hat{\psi}-\tilde{\varepsilon}\partial_{z}^{2}\hat{\psi}$
$\displaystyle=\frac{\gamma+1}{c^{4}}(\omega_{1}-\omega_{2})\omega_{1}\omega_{2}\imath\left(\sum_{k=1}^{n}\hat{\phi}_{k}\right)\left(\sum_{l=n+1}^{n+m}\overline{\hat{\phi}}_{l}\right)$
(36)

for $k\in\left\\{1,\ldots n\right\\}$, $l\in\left\\{n+1,\ldots n+m\right\\}$
and arbitrary $m,n\in\mathbb{N}$. Taking multiple beams into account can be of
interest to reduce the ill-posedness of the inverse problem.

##  3 The inverse problem

The inverse problem of nonlinearity parameter imaging will be considered for
one of these models, namely Case 2A

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

on $\tilde{\Omega}=(0,L)\times\tilde{\Omega}_{y}$, where $s=\frac{1}{c}$ and
$\eta=\frac{\gamma+1}{c^{4}}$ with boundary conditions

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

The parameter we want to identify is $\eta$, which is directly connected to
the space dependent nonlinearity parameter $B/A$, since

$\eta=\frac{\gamma+1}{c^{4}}=\frac{2\beta}{c^{4}}=\frac{2+B/A}{c^{4}}.$

In well-known nonlinear wave equations the nonlinearity parameter is often
denoted by $\beta$ in the context of liquids or $\gamma$ for gases.

We assume the speed of sound $c$ to be constant. Identifying the nonlinearity
parameter and the speed of sound simultaneously leads to a nonlinear inverse
problem that is dealt with in [10].

Since the excited low-frequency field in terms of $\psi$ propagates omni-
directionally, alternatively to the paraxial form of the partial differential
equation in terms of $\hat{\psi}$ one can also consider the Helmholtz equation

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

for $\check{\psi}=P^{-1}\hat{\psi}$, $\check{\eta}=P^{-1}\eta$ (cf. (6)) with
boundary conditions

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

with $P(\Omega)\supseteq\tilde{\Omega}=(0,L)\times\tilde{\Omega}_{y}$, where
$P$ is the paraxial transform defined by (1) and $1_{\Omega}$ the extension by
zero operator for a function defined on $P^{-1}(\tilde{\Omega})$ to all of
$\Omega$. This would allow to use a possibly larger propagation domain for the
$\psi$ wave as compared to the one for the $\phi_{1}$, $\phi_{2}$ waves, but
come at the price of a more involved implementation.

In the following we deal with the inverse problem of nonlinearity imaging with
the alternative formulation in (39). Since $\check{\eta}$ only appears in the
second equation, only this equation is relevant for reconstruction, so that
the system of equations reduces to an inhomogeneous Helmholtz equation and
$\hat{\phi}_{1},\,\hat{\phi}_{2}$ can be computed independently beforehand.

###  3.1 Formulation of the inverse problem

For the case of $s$ being known and constant, consider the reconstruction of
$\check{\eta}=\check{\eta}(x)$ in

$\displaystyle-\omega_{d}^{2}s^{2}\check{\psi}-\Delta_{x}\check{\psi}=\check{\eta}\,f$
$\displaystyle\mbox{ in }\Omega$ (41)
$\displaystyle\partial_{\nu}\check{\psi}=-\imath\sigma\check{\psi}$
$\displaystyle\mbox{ on }\partial\Omega$

with
$f=\imath\omega_{p}\,P^{-1}\bigl{(}\hat{\phi}_{1}\,\overline{\hat{\phi}}_{2}\bigr{)}$
cf. (39) from observations

$\check{p}^{\textup{meas}}=\imath\omega_{d}\mbox{tr}_{\check{\Gamma}}\check{\psi}\,,$
(42)

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

This inverse problem can be written as a linear operator equation

$T\check{\eta}=y$ (43)

with $y=\check{p}^{\textup{meas}}$ and with the forward operator $T=C\circ S$
being a concatenation of the parameter-to-state map

$S:\check{\eta}\mapsto\check{\psi}\text{ solving \eqref{bvp_psicheck}}$

with the (linear) observation operator

$C:\check{\psi}\mapsto\imath\omega_{d}\mbox{tr}_{\check{\Gamma}}\check{\psi}.$

Note that since $(\hat{\phi}_{1},\hat{\phi}_{2})$ can be computed
independently of $\eta$, the forward operator $S$ is linear. So, this means
that the linear operator $T$ is given by

$T:L^{2}(\Omega)\to
L^{2}(\check{\Gamma}),\quad\check{\eta}\mapsto\imath\omega_{d}\mbox{tr}_{\check{\Gamma}}\check{\psi},$
(44)

where $\check{\psi}$ solves (41). Well-definedness and continuity of the
linear operator $T$ for $f\in L^{p}(\Omega)$ follows from known results on
well-posedness of the Helmholtz equation with impedance boundary conditions;
cf. e.g. [13, 14] and the references therein.

###  3.2 Solving the inverse problem

For solving the inverse problem we want to use the Landweber-Kaczmarz method.
A Landweber step for solving (43) is defined by a gradient descent step for
the least squares functional

$\left\|T(\check{\eta}^{(n)})-y\right\|^{2}_{L^{2}}$

i.e., by

$\check{\eta}^{(n+1)}=\check{\eta}^{(n)}+\mu T^{*}(y-T(\check{\eta}^{(n)})).$
(45)

with an appropriately chosen step size $\mu$. Therefore, we want to derive the
adjoint operator of $T(\check{\eta})$.

Observe that for all $\check{\eta}\in L^{2}(\Omega),\,q\in
L^{2}(\check{\Gamma})$ it holds that

$\left\langle
T\check{\eta},q\right\rangle_{L^{2}(\check{\Gamma})}=\int_{\check{\Gamma}}T\check{\eta}\cdot\overline{q}\differential{\check{\Gamma}}=\int_{\check{\Gamma}}\imath\omega_{d}\check{\psi}\cdot\overline{q}\differential{\check{\Gamma}}=\int_{\Omega}\check{\eta}\cdot
T^{*}q\differential{x}=\left\langle\check{\eta},T^{*}q\right\rangle_{L^{2}(\Omega)}$

and where $\check{\psi}\in H^{1}(\Omega;\mathbb{C})$ satisfies the weak direct
equation

$\displaystyle\int_{\partial\Omega}\imath\sigma\check{\psi}\overline{w}\differential{S}+\int_{\Omega}\nabla\check{\psi}\nabla\overline{w}\differential{x}-\int_{\Omega}\omega_{d}^{2}s^{2}\check{\psi}\overline{w}\differential{x}=\int_{\Omega}\check{\eta}\,f\overline{w}\differential{x}\qquad\text{for
all }w\in H^{1}(\Omega;\mathbb{C}).$

Thus $T^{*}q=\overline{f}r$ where $r$ solves the weak form of the adjoint
equation

$\int_{\partial\Omega\setminus\check{\Gamma}}\imath\sigma
v\overline{r}\differential{S}+\int_{\Omega}\nabla
v\nabla\overline{r}\differential{x}-\int_{\Omega}\omega_{d}^{2}s^{2}v\overline{r}\differential{x}=\int_{\check{\Gamma}}\imath\omega_{d}v\overline{q}\differential{\check{\Gamma}}\qquad\text{for
all }v\in H^{1}(\Omega;\mathbb{C})$

and the adjoint operator is given by

$T^{*}:L^{2}(\check{\Gamma})\to L^{2}(\Omega),\quad q\mapsto f\overline{r},$
(46)

where $r$ solves

$\displaystyle-\omega_{d}^{2}s^{2}r-\Delta r$ $\displaystyle=0\quad\mbox{ in
}\Omega$ $\displaystyle\partial_{\nu}r$ $\displaystyle=\imath\sigma
r\quad\mbox{ on }\partial\Omega\setminus\check{\Gamma}$ (47)
$\displaystyle\left[\partial_{\nu}r\right]_{\check{\Gamma}}=-\imath\omega_{d}q\quad\text{
on }\check{\Gamma}\cap\Omega,$
$\displaystyle\qquad\partial_{\nu}r=-\imath\omega_{d}q\quad\text{ on
}\check{\Gamma}\cap\partial\Omega,$

where $\left[\partial_{\nu}r\right]$ denotes the jump of the normal derivative
over the interface $\check{\Gamma}\cap\Omega$ (in case $\check{\Gamma}$ or
part of it is immersed in the computational domain $\Omega$).

Finally, to recover $\check{\eta}$ we take into account measurements at
multiple frequencies $\omega_{1,j}$, $\omega_{2,j}$,
$j\in\left\\{1,\ldots,J\right\\}$ and also several receiver array locations
$\check{\Gamma}^{m}$, $m\in\left\\{1,\ldots,M\right\\}$ might be used. Note
that for elementary dimensionality reasons a general function
$\check{\eta}=\check{\eta}(x)$ cannot be expected to be uniquely determined
from a single boundary observation.

Thus, we actually deal with a set of several model and observation operators,
which gives the forward operator $T_{p}=C^{m}\circ S^{j}$ and data $y_{p}$ for
$p=(m-1)J+j$. So, we can write the inverse problem of reconstructing
$\check{\eta}$ as a system of operator equations

$T_{p}(\check{\eta})=y_{p},\qquad p\in\left\\{1,\ldots,P=J\cdot M\right\\}$

and apply Kaczmarz type method by sequentially performing one Landweber step
in a cyclically repeated manner

$\check{\eta}^{(n+1)}=\check{\eta}^{(n)}+\mu
T_{p}^{*}(y_{p}-T_{p}(\check{\eta}^{(n)})),$ (48)

where $p=\mod(n-1,P)+1$.

##  4 Algorithms

We present algorithms for solving the direct and also the inverse problem for
the model in Case 2A given by (17), (18). Therefore, the following modules are
used as subroutines to the algorithms:

  * •

$\hat{\phi}_{k}=$solveParaxialPDE$(\omega_{k},s,\sigma_{k},h_{k})$: solve

$\displaystyle
2\imath\omega_{k}s\partial_{z}\hat{\phi}_{k}-\Delta_{y}\hat{\phi}_{k}=0\qquad$
$\displaystyle\mbox{ in }\tilde{\Omega},$
$\displaystyle\hat{\phi}_{k}(0,\cdot)=\hat{h}_{k}\qquad$ $\displaystyle\mbox{
in }\tilde{\Omega}_{y}\,,$ (49)
$\displaystyle\partial_{\nu_{y}}\hat{\phi}_{k}(z,y)=-\imath\sigma_{k}\hat{\phi}_{k}\qquad$
$\displaystyle\mbox{ in }(0,L)\times\partial\tilde{\Omega}_{y};$

  * •

$\check{\psi}=$solveHelmholtz$(\omega_{k},s,\sigma,\check{\eta},f)$: solve

$\displaystyle-\omega_{d}^{2}s^{2}\check{\psi}-\Delta\check{\psi}=\check{\eta}\,f\qquad$
$\displaystyle\mbox{ in }\Omega,$
$\displaystyle\partial_{\nu}\check{\psi}=-\imath\sigma\check{\psi}\qquad$
$\displaystyle\mbox{ on }\partial\Omega,$

with
$f=\imath\omega_{p}\,P^{-1}\bigl{(}\hat{\phi}_{1}\,\overline{\hat{\phi}}_{2}\bigr{)}$;

  * •

$r=$solveHelmholtz$\Gamma(\omega_{k},s,\sigma,q)$: solve

$\displaystyle-\omega_{d}^{2}s^{2}r-\Delta r$ $\displaystyle=0\quad\mbox{ in
}\Omega$ $\displaystyle\partial_{\nu}r$ $\displaystyle=\imath\sigma
r\quad\mbox{ on }\partial\Omega\setminus\check{\Gamma}$
$\displaystyle\left[\partial_{\nu}r\right]_{\check{\Gamma}}=-\imath\omega_{d}q\quad\text{
on }\check{\Gamma}\cap\Omega,$
$\displaystyle\qquad\partial_{\nu}r=-\imath\omega_{d}q\quad\text{ on
}\check{\Gamma}\cap\partial\Omega,$

The submodules
$\check{\psi}=$solveHelmholtz$(\omega_{k},s,\sigma,\check{\eta},f)$ and
$r=$solveHelmholtz$\Gamma(\omega_{k},s,\sigma,q)$ can be implemented using
standard solvers for Helmholtz equation.

For the submodule $\hat{\phi}_{k}=$solveParaxialPDE$(...)$ an individual
algorithm has been designed. We assume that the beams are excited by a linear
transducer, which means that in one dimension the transducer can be
represented as $\tilde{\Omega}_{k}=[-l_{k},l_{k}]\subseteq\tilde{\Omega}_{y}$
with $l_{k}>0$. We can make use of this geometry and consider the
corresponding eigenvalue problem for the Laplace operator $\Delta_{y}$ with
Neumann boundary conditions, i.e.

$\displaystyle-\Delta_{y}v_{k}^{j}(y)=\lambda_{k}^{j}v_{k}^{j}(y)$
$\displaystyle\qquad y\in\tilde{\Omega}_{k},$
$\displaystyle\partial_{\nu}v_{k}^{j}(y)=0$ $\displaystyle\qquad
y\in\partial\tilde{\Omega}_{k}.$

for $k=1,2$, $j=1,2,\ldots$. Then, the explicit solution is given by

$v_{k}^{j}(y)=\sqrt{\frac{1}{l_{k}}}\cos(j\pi\frac{y}{2\,l_{k}}),\qquad\lambda_{k}^{j}=\frac{\pi^{2}j^{2}}{4\,l_{k}^{2}}.$

Generalization to higher dimensions can be achieved easily, since on a
rectangular domain separation of variables is possible, e.g. for two
dimensions

$v^{mn}(y_{1},y_{2})=v^{m}(y_{1})v^{n}(y_{2}),\qquad\lambda^{mn}=\lambda^{m}+\lambda^{n}$

where $v^{m}(y_{1})$ and $\lambda^{m}$ are the eigenfunctions and eigenvalues
corresponding to the one-dimensional problem on $[-b_{1},b_{1}\,]$,
respectively, $v^{n}(y_{2})$ and $\lambda^{n}$ are those corresponding to the
one-dimensional problem on$[-b_{2},b_{2}\,]$ with $b_{1},\,b_{2}\,>0$.
Furthermore, this can simply be adapted to other boundary conditions and
geometries such as disks or spheres, where the explicit eigensystem of the
eigenvalue problem of the Laplace operator is known, cf. [5].

Overall, the velocity potential $\hat{\phi}_{k}$ at any position $(z,y)$ is
given by

$\hat{\phi}_{k}(z,y)=\sum_{j}a_{k}^{j}(z)v_{k}^{j}(y).$ (50)

Now, plugging this Ansatz into equation (• ‣ 4) yields

$\sum_{j}\left(2\imath\omega_{k}sa_{z}^{j}(z)+\lambda^{j}a^{j}(z)\right)v^{j}_{k}(y)=0,$

which means that we have to solve the ODE

$a_{z}^{j}=\imath\frac{\lambda_{k}^{j}}{2\omega_{k}s}a^{j}.$

The explicit solution is then given by

$a^{j}=a^{j}(0)e^{\imath\frac{\lambda_{k}^{j}}{2\omega_{k}s}z}.$ (51)

with

$a^{j}(0)=\frac{1}{2l_{k}}\int_{-l_{k}}^{l_{k}}\hat{\phi}_{k}(0,y)v^{j}_{k}(y)\differential{y}.$

Therefore, $\hat{\phi}_{k}$ is given by

$\displaystyle\hat{\phi}_{k}(z,y)=\sum_{j}\frac{1}{2l_{k}}\int_{-l_{k}}^{l_{k}}\hat{\phi}_{k}(0,y)v^{j}_{k}(y)\differential{y}e^{\imath\frac{\lambda_{k}^{j}}{2\omega_{k}s}z}v_{k}^{j}(y),\qquad
k\in\left\\{1,2\right\\}.$

In the discrete case calculating the initial data corresponds to applying a
discrete cosine transform (DCT) and using the Ansatz means applying the
inverse discrete cosine transform (IDCT).

Input: $\hat{h}_{k}(0,y)$ | % transducer excitation
---|---
Output: $\hat{\phi}_{k}(z_{\text{out}},y)$ |
$a^{j}_{k}(0)=\text{DCT}(\hat{h}_{k}(0,y))$ | % Ansatz for linear transducer
for $z=0:\Delta z:z_{\text{out}}$ |
$a^{j}_{k}(z)=a^{j}_{k}(0)e^{\imath\frac{\lambda_{k}^{j}}{2\omega_{k}s}z}$ |
$\hat{\phi}_{k}(z_{\text{out}},y)=\text{IDCT}(a^{j}_{k}(z_{\text{out}}))$ |
Table 3: Pseudo-Code for
$\hat{\phi}_{k}=$solveParaxialPDE$(\omega_{k},s,\sigma_{k},\hat{h}_{k})$

###  4.1 Direct problem

The direct problem consists of solving (39), (40) for $\check{\psi}$, where
$\check{\eta}$ is known. The initial transducer signal is given by
$\hat{h}_{k}$ and designed such that the beam focuses, which can be achieved
by time delayed pulses, see section 2.4 for more details.

Input: $\hat{h}_{k}(0,y)$ | % transducer excitation
---|---
Output: $\check{\psi}(x)$ |
for $z=0:\Delta z:z_{\text{out}}$ |
$\hat{\phi}_{k}(z,\tau)=$solveParaxialPDE$(\omega_{k},s,\sigma_{k},h_{k})$ | % paraxial variables
$\check{\phi}_{12}=P^{-1}\bigl{(}\hat{\phi}_{1}\,\overline{\hat{\phi}}_{2}\bigr{)}$ | % change of variables
for $x=x_{\text{in}}:\Delta x:x_{\text{out}}$ |
$\check{\psi}(x)=$solveHelmholtz$(\omega_{k},s,\sigma,\check{\eta},\check{\phi}_{12})$ | % original variables
Table 4: Pseudo-Code for the forward problem

###  4.2 Inverse problem

For the inverse problem we want to reconstruct $\check{\eta}$ in (41) from
given observations, cf. (42). Since $\check{\eta}$ only appears in the
equation for $\check{\psi}$, we can compute $\hat{\phi}_{2},\,\hat{\phi}_{2}$
separately as presented above, apply the inverse operator for the paraxial
change of variables and use this as an input for the inhomogeneous Helmholtz
equation in the following system.

Input: $\check{\eta}^{(n)},\,\check{\phi}_{12}$ |
---|---
Output: $\check{\eta}^{(n+1)}$ |
$\check{\psi}^{(n)}=$solveHelmholtz$(\omega_{k},s,\sigma,\check{\eta},\check{\phi}_{12})$ | % computation of $\check{\psi}=S(\check{\eta}^{(n)})$
$q=y-\mbox{tr}_{\check{\Gamma}}\check{\psi^{(n)}}$ | % computation of adjoint states
$r^{(n)}=$solveHelmholtz$\Gamma(\omega_{k},s,\sigma,q)$ |
$\xi^{(n)}=\imath\omega_{p}\check{\phi}_{12}r^{(n)}$ |
set $\check{\eta}^{(n+1)}=\check{\eta}^{(n)}+\mu\xi^{(n)}$ |
Table 5: Pseudo-Code for the inverse problem

##  5 Outlook

In this paper we have introduced a new approach for directive models in vibro-
acoustic imaging, where we used a paraxial change of variables. We have
transferred the systems of equations to frequency domain, studied the
identification of the nonlinearity parameter for one of these models, and
derived a Landweber-Kaczmarz method for reconstruction.

Concerning forward simulation we have presented a numerical method that takes
the geometry of the transducer into account and explicitly solves the paraxial
partial differential equation. Next steps are to implement the devised
algorithms for the direct problem and the Landweber-Kaczmarz method for the
reconstruction of the nonlinearity parameter.

Future work would also consist of adapting the numerical method for the direct
problem and also the inverse problem for variable speed of sound and variable
nonlinearity parameter. While the reconstruction of the nonlinearity parameter
alone yields a linear inverse problem where the first equation of the system
in $\phi_{1}$, $\phi_{2}$ can be computed independently and the inverse
problem then simply consists of the second equation, for the simultaneous
reconstruction of both parameters the whole system has to be taken into
account, which results in a nonlinear inverse problem.

## Acknowledgment

The author wishes to thank her supervisor Barbara Kaltenbacher, University of
Klagenfurt, for fruitful and inspiring discussions. This work was supported by
the Austrian Science Fund (FWF) under the grant DOC78.

## References

  * [1] L. Bjørnø.  Characterization of biological media by means of their non-linearity.  Ultrasonics, 24(5):254–259, 1986.
  * [2] Charles A. Cain.  Ultrasonic reflection mode imaging of the nonlinear parameter b/a: I. a theoretical basis.  The Journal of the Acoustical Society of America, 80(1):28–32, 1986\.
  * [3] Mostafa Fatemi and James F. Greenleaf.  Ultrasound-stimulated vibro-acoustic spectrography.  Science, 280(5360):82–85, 1998.
  * [4] Mostafa Fatemi and James F. Greenleaf.  Vibro-acoustography: An imaging modality based on ultrasound-stimulated acoustic emission.  Proceedings of the National Academy of Sciences, 96(12):6603–6608, 1999.
  * [5] D. S. Grebenkov and B.-T. Nguyen.  Geometrical structure of laplacian eigenfunctions.  SIAM Review, 55(4):601–667, jan 2013.
  * [6] M.F. Hamilton and D.T. Blackstock.  Nonlinear Acoustics.  Academic Press, 1998.
  * [7] Nobuyuki Ichida, Takuso Sato, and Melvin Linzer.  Imaging the nonlinear ultrasonic parameter of a medium.  Ultrasonic Imaging, 5(4):295–299, 1983.  PMID: 6686896.
  * [8] F.B. Jensen, W.A. Kuperman, M.B. Porter, and H. Schmidt.  Computational Ocean Acoustics.  Modern Acoustics and Signal Processing. Springer New York, 2011.
  * [9] Barbara Kaltenbacher.  On the inverse problem of vibro-acoustography.  Meccanica, pages 1–12, 2022.  arXiv:2109.01907 [math.AP].
  * [10] Barbara Kaltenbacher and Teresa Rauscher.  Simultaneous reconstruction of sound speed and nonlinearity parameter in a paraxial model of vibro-acoustography in frequency domain, submitted 2023\.  arXiv:2303.17236 [math.NA].
  * [11] A.E. Malcolm, F. Reitich, J. Yang, J.F. Greenleaf, and M. Fatemi.  A combined parabolic-integral equation approach to the acoustic simulation of vibro-acoustic imaging.  Ultrasonics, 48(6):553–558, 2008.  Selected Papers from ICU 2007.
  * [12] A.E. Malcolm, F. Reitich, J. Yang, J.F. Greenleaf, and M. Fatemi.  Numerical Modeling for Assessment and Design of Ultrasound Vibroacoustography Systems.  2008\.
  * [13] Jens Markus Melenk.  On generalized finite-element methods.  ProQuest LLC, Ann Arbor, MI, 1995.  Thesis (Ph.D.)–University of Maryland, College Park.
  * [14] Andrea Moiola and Euan A. Spence.  Is the Helmholtz equation really sign-indefinite?  SIAM Rev., 56(2):274–312, 2014.
  * [15] Anna Rozanova.  The Khokhlov-Zabolotskaya-Kuznetsov equation.  Comptes Rendus Mathematique, 344(5):337–342, 2007.
  * [16] T.L. Szabo.  Diagnostic Ultrasound Imaging: Inside Out.  Biomedical Engineering. Elsevier Science, 2013.
  * [17] Fred D. Tappert.  The parabolic approximation method, pages 224–287.  Springer Berlin Heidelberg, Berlin, Heidelberg, 1977.
  * [18] Matthew Urban, Carl Chalek, Bruno Haider, Kai Thomenius, Mostafa Fatemi, and Azra Alizad.  A beamforming study for implementation of vibro-acoustography with a 1.75-d array transducer.  IEEE transactions on ultrasonics, ferroelectrics, and frequency control, 60:535–51, 03 2013.
  * [19] Matthew W. Urban, Carl Chalek, Randall R. Kinnick, Thomas M. Kinter, Bruno Haider, James F. Greenleaf, Kai E. Thomenius, and Mostafa Fatemi.  Implementation of vibro-acoustography on a clinical ultrasound system.  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 58(6):1169–1181, 2011.
  * [20] M.W. Urban, G.T. Silva, M. Fatemi, and J.F. Greenleaf.  Multifrequency vibro-acoustography.  IEEE Transactions on Medical Imaging, 25(10):1284–1295, 2006.
  * [21] Evgenia A. Zabolotskaya.  Quasi-plane waves in the nonlinear acoustics of confined beams.  volume 15, pages 35–40, 1969.