# Quantitative inverse problem in visco-acoustic media under attenuation model uncertainty

Florian Faucher, Otmar Scherzer

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

## To cite this version:

Florian Faucher, Otmar Scherzer. Quantitative inverse problem in visco-acoustic media under attenuation model uncertainty. Journal of Computational Physics, 2023, 472, pp.111685. 10.1016/j.jcp.2022.111685. hal-03871831

HAL Id: hal-03871831

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

Submitted on 17 Sep 2024

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

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

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

ELSEVIER

Contents lists available at [ScienceDirect](https://www.sciencedirect.com)

Journal of Computational Physics

journal homepage: [www.elsevier.com/locate/jcp](https://www.elsevier.com/locate/jcp)

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

# Quantitative inverse problem in visco-acoustic media under attenuation model uncertainty

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

Check for
updates

Florian Faucher<sup>a,b,*</sup>, Otmar Scherzer<sup>a,c,d</sup>

## ARTICLE INFO

**Article history:**

Received 30 April 2022

Received in revised form 19 September 2022

Accepted 6 October 2022

Available online 19 October 2022

**Keywords:**

Inverse problem

Visco-acoustic models

Time-harmonic wave propagation

Full waveform inversion

## ABSTRACT

We consider the inverse problem of quantitative reconstruction of properties (e.g., bulk modulus, density) of visco-acoustic materials based on measurements of responding waves after stimulation of the medium. Numerical reconstruction is performed by an iterative minimization algorithm. Firstly, we investigate the robustness of the algorithm with respect to attenuation model uncertainty, that is, when different attenuation models are used to simulate synthetic observation data and for the inversion, respectively. Secondly, to handle data-sets with multiple reflections generated by wall boundaries around the domain, we perform inversion using complex frequencies, and show that it offers a robust framework that alleviates the difficulties of multiple reflections. To illustrate the efficiency of the algorithm, we perform numerical simulations of ultrasound imaging experiments to reconstruct a synthetic breast sample that contains an inclusion of high-contrast properties. We perform experiments in two and three dimensions, where the latter also serves to demonstrate the numerical feasibility in a large-scale configuration.

© 2022 The Author(s). Published by Elsevier Inc. This is an open access article under the

CC BY license ([http://creativecommons.org/licenses/by/4.0/](http://creativecommons.org/licenses/by/4.0/)).

## 1. Introduction

Attenuation refers to the loss of energy of waves propagating in a medium. The level of energy loss is quantified as a material property of a viscous medium. In the context of imaging with waves, visco-acoustic media represent additional challenges compared to an ideal (that is a non-attenuating) medium: For instance, the mathematical description of such materials involves more parameters, e.g., to account for the attenuation mechanisms. In this work we investigate the quantitative reconstruction of properties of a visco-acoustic medium, which is relevant, for instance, in ultrasound imaging, non-destructive testing, or seismic imaging, e.g., [15,27,65,34,5,51,3]. In such a setup, probing waves are sent to the area of interest, and the mechanical response of the medium is indirectly measured at the position of the receivers.

Attenuation is a phenomenon that intrinsically depends on the frequencies of the propagating waves, where each frequency component loses a different amount of energy [7,1,11]. It is common to categorize attenuation models into different families, such as those that dominantly attenuate low or high frequency components, see Section 2. Generalized models combine multiple mechanisms of attenuation and enable the consideration of signals with wide frequency ranges, [11,67], see

<footnote><sup>a</sup> Faculty of Mathematics, University of Vienna, Oskar-Morgenstern-Platz 1, A-1090 Vienna, Austria</footnote>

<footnote><sup>b</sup> Project-Team Makutu, Inria Bordeaux, Université de Pau et des Pays de l'Adour, France</footnote>

<footnote><sup>c</sup> Johann Radon Institute for Computational and Applied Mathematics (RICAM), Altenbergerstraße 69, A-4040 Linz, Austria</footnote>

<footnote><sup>d</sup> Christian Doppler Laboratory for Mathematical Modeling and Simulation of Next Generations of Ultrasound Devices (MaMSi), Oskar-Morgenstern-Platz 1, A-1090 Vienna, Austria</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: [florian.faucher@inria.fr](mailto:florian.faucher@inria.fr) (F. Faucher), [otmar.scherzer@univie.ac.at](mailto:otmar.scherzer@univie.ac.at) (O. Scherzer).</footnote>
Remark 2. Therefore, the wave equations that model the propagation of attenuated waves can have different forms, as introduced, for instance, in [62,11]. From a practical perspective, it is difficult to identify the appropriate mathematical model of attenuation that adequately describes a given material, cf. [11,57] for geophysical applications or [70] for medical ones. Consequently, in practical applications of inverse problems, the appropriate attenuation model to describe the material behavior represents an additional unknown. Note that despite the different nature of the media and of the acquisition setups, cf. [53], geophysical and medical imaging share the same models of visco-acoustic wave propagation.

In our work we assume that the excitation that generates the waves within the medium is a broad-band time-signal. Then, we apply a Fourier transform to the time measurements to work with frequency-domain data, as we illustrate in Subsection 2.1. To simulate the waves that will later be compared with the measurement data used to reconstruct the medium's properties, we use the frequency-domain wave propagation equations. The frequency-domain formulation has been recognized to be more convenient to handle attenuation, cf. [7,47,55], as it allows us to work with Partial Differential Equations (PDEs) regardless of the attenuation model. On the other hand, the time-domain formulation can lead to integro-differential equations instead. This is due to the fact that in the frequency-domain the attenuation is encoded via complex-valued parameters in the PDEs, e.g., [47,12,62]. However, causality principles must be carefully addressed in the frequency-domain as advocated by [31,11,18]. The level of attenuation of a medium is quantified by the frequency-dependent *quality factor*, see Definition 1. The higher it is, the less attenuating is the medium. In our experiments, we consider the reconstruction of breast tissues at ultrasound frequencies, where the quality factor is higher than 100 (see Table 3, Sections 4 and 5): this corresponds to weakly attenuating medium, [64].

For the quantitative reconstruction of medium properties, our experiments follow the principles of ultrasound tomography, which have been applied, for instance, in the context of breast imaging, e.g., [16,17,41], where simplifying assumptions regarding the models of wave propagation can be used. Reconstruction methods are further compared in [50,24]. We further use the full model of wave propagation for the simulations, and the reconstructions are carried out with an iterative method that minimizes a misfit function defined to evaluate the distance between the observed measurements and simulations. This approach is usually referred to as the *Full Waveform Inversion* (FWI) in seismic imaging, [40,61,54,65]. As its name indicates, FWI relies on the full measurement signals, e.g., contrary to approaches that only use the travel-times.

Imaging using the full waveform has been used in the context of ultrasound tomography, for instance with special emphasis on bone structures, [6], for the brain, [32] and for breast imaging considering the breast as a viscous medium, with implementations both in the frequency, [42], and time-domain, [53,51,3,44]. We can highlight two major differences between our work and these references:

1. The aforementioned references rely on a single attenuation model, which is fixed a-priori, prior to the reconstruction.

2. These papers assume wave propagation in free-space, that is, assuming that no reflection comes from the boundary of the acquisition setup. For practical experiments this requires particular experimental conditions, such as adding padding materials, or a non-trivial data post-processing step ([68]) to remove the multiple reflections coming from imperfect boundaries while preserving the ones that come from the sample.

More specifically, the contributions of our work are the following:

- We review and implement seven different attenuation models describing visco-acoustic wave propagation. Each model encompasses a different effect of wave dissipation and dispersion.

- We carry out inversion with *attenuation model uncertainty*, that is, we use a different attenuation model to generate the synthetic data (simulating measurement data) and to carry out the numerical reconstruction procedure. Despite the resulting changes in the forward PDE models, we show that FWI is a robust approach that does not suffer from inconsistency in the attenuation model.

- We start with experiments where we impose *absorbing boundary conditions* to constrain the numerical domain, hence mimicking a free-space wave propagation. We then investigate the consideration of wall boundary on the sides of the domain. The resulting multiple reflections are shown to strongly influence the accuracy of the reconstructions. To overcome the difficulty, we use complex frequencies. This approach is also referred to as the *Laplace-Fourier domain method*, [58,59] and it is shown to improve the convexity of the misfit function for inversion in [22,19]. In our work, we show that this transformation can alleviate the difficulty that occurs from the multiple reflections coming from wall boundaries, by enhancing wave first arrivals, [34].

- In the context of multi-parameter inversion, we investigate the choice of parametrization, that is, the choice of model parameters with respect to whom the gradient is computed, [8,37,19]. In particular, the density and attenuation properties are known to be hard to reconstruct, [65,33], and possibly require a specific misfit function, [35].

- We perform experiments in three dimensions to explore the feasibility of our methodology and detail the computational cost.

In Section 2, we provide the visco-acoustic equation for wave propagation in the frequency domain, and review seven different attenuation models from the literature, [62,11,18]. We highlight the benefit of working in the frequency domain, in which case all the attenuation models are encoded as PDEs with complex-valued parameters. Therefore the frequency domain allows for using PDE software infrastructure, which the time domain does not allow. The inverse procedure for
<figure><img src="image_6.png" /><figcaption>Fig. 1. Non-invasive measurement setup with Ricker source excitation: the data are obtained from a single point-source excitation or a transducer device, which corresponds to an array of fixed length composed of multiple point-sources, which are simultaneously excited. The source function is a Ricker-wavelet, which is a broad-band signal. The receivers are positioned around the sample (black dashed line). (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)</figcaption></figure>

the reconstruction of parameters using FWI is detailed in Section 3. In Section 4, we carry out reconstructions for a two-dimensional sample corresponding to a breast model, while a three-dimensional case is performed in Section 5, where we also provide the computational cost of the method.

## 2. Visco-acoustic forward wave problem

### 2.1. From time-domain to complex frequency-domain data-sets

We investigate the reconstruction of material parameters of a heterogeneous visco-acoustic medium. Waves are excited at the boundary of a sample and propagate through it; receivers positioned outside of the sample measure the response wavefield which is used to reconstruct the medium properties. In particular, this means that we investigate a non-invasive experiment. We assume that the exciting source has a compact support. Mathematically, it can be represented by a delta-Dirac in space and a time-domain wavelet: $\delta(x) S(t)$. We refer to such a function as a point-source, see Fig. 1a. We can also consider line-source excitations, see Fig. 1b. Mathematically it is written as a composition of point-sources that are excited at the same time-instant. In this work, the time-domain signal/wavelet $S(t)$ of the point-source excitation follows a Ricker wavelet function, which represents a broad-band signal containing multiple frequencies, as illustrated in Fig. 1d.

Time-domain signals are transformed with a Laplace-Fourier transform [59,22,25], introducing the complex frequency $\omega$ such that,

$$ \omega := \omega_R + i\omega_I \Rightarrow i\omega = i\omega_R - \omega_I \quad \text{with} \quad \omega_R > 0 \text{ and } \omega_I \ge 0. \qquad (2.1) $$

Namely, the discrete transform of a signal $S$ composed of $N_t$ time-steps such that $S = \{S_0, S_1, S_2, \dots, S_{N_t-1}\}$ is computed with

$$ \tilde{\mathcal{F}}(S)(\omega) = \sum_{k=0}^{N_t-1} S_k e^{i\omega k/N_t}, \qquad (2.2) $$

using the same convention as [1, Box 5.2]. Note that when the imaginary part is zero, $\omega_I = 0$, it corresponds to the ordinary frequency $\nu$ Hz, with $\omega_R = 2\pi\nu$. In the experiments, receivers record time signals (e.g., pressure fields), which are thus first transformed to complex-frequency data, see Fig. 2.

In the following of the paper, we only study wave propagation in the frequency-domain, or time-harmonic waves, omitting the preliminary step that consists in transforming the measured time-domain signals. As illustrated in Fig. 1d, the available frequency content depends on the Ricker source frequency peak.
<figure><img src="image_7.png" /></figure>

a) Data are obtained in the time-domain (left) from a point-source excitation (Figure 1a), and a complex Fourier transform is applied to produce the complex frequency data (right). Here, $y_t$ is the signal at time step $t$ with $N_t$ time steps in total, and $x_{rcv}$ is the position of a receiver.

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

b) Time-domain signals using a transducer source (Figure 1b) of length 3cm, composed of 128 Ricker point-sources simultaneously excited

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

c) Time-domain signals using a transducer source (Figure 1b) of length 9cm, composed of 128 Ricker-point sources simultaneously excited

**Fig. 2.** Illustration of the steps from time-domain measurements to the complex frequency data in a two-dimensional setup depending on the source excitation. The source is positioned at angle 0°, corresponding to the Cartesian coordinates $x = 17.5$ cm, $z = 9$ cm in Fig. 1. Inversion experiments following this setup are carried out in Section 4, 3D experiments are in Section 5.

## 2.2. Frequency-domain visco-acoustic wave equations

Following the complex Fourier transform of the time-domain data, we consider the propagation of waves in the frequency domain. Let us consider the domain $\Omega$ in dimension $d$ (2 or 3 in our experiments), with boundary $\Gamma$. The propagation of waves in a visco-acoustic medium is described by the particle velocity vector field $v: \Omega \to \mathbb{C}^d(\Omega)$ and the scalar pressure field $p: \Omega \to \mathbb{C}(\Omega)$, that satisfy the first-order system of equations, [25]:

$$ -i\omega \rho(\mathbf{x})\mathbf{v}(\mathbf{x}, \omega) + \nabla p(\mathbf{x}, \omega) = 0, \qquad (2.3a) $$

$$ -\frac{i\omega}{\kappa_{\dagger}(\mathbf{x}, \omega)} p(\mathbf{x}, \omega) + \nabla \cdot \mathbf{v}(\mathbf{x}, \omega) = g(\mathbf{x}, \omega), \qquad (2.3b) $$

where $\mathbf{x}$ is the space-coordinates and $g$ the interior source function. The derivation of the complex-frequency time-harmonic wave equation is further detailed in Appendix A.

The material properties describing the medium consist of the density $\rho: \Omega \to \mathbb{R}(\Omega)$ and of the bulk modulus $\kappa_{\dagger}: \Omega \to \mathbb{C}(\Omega)$, which is complex-valued and frequency dependent to acknowledge the effects of attenuation. The complex wave speed $c_{\dagger}: \Omega \to \mathbb{C}(\Omega)$ is defined from these two coefficients by

$$ c_{\dagger}(\mathbf{x}, \omega) = \sqrt{\frac{\kappa_{\dagger}(\mathbf{x}, \omega)}{\rho(\mathbf{x})}}, \qquad (2.4) $$
where $\sqrt{\gamma}$ uses the principal argument branch $[0, 2\pi)$, cf. [18, Section 3], that is

$$ \kappa_{\uparrow}^{\beta} = |\kappa_{\uparrow}|^{\beta} (\beta \arg(\kappa_{\uparrow})). \qquad (2.5) $$

**Free-space propagation** We refer to free-space propagation when waves are supposed to propagate up to infinity, without generating reflections. This is typically assumed in the ultrasound imaging setup, e.g., [53,51,3]. In this case, we impose artificial absorbing boundary conditions (corresponding to Sommerfeld's radiation condition) to constrain the numerical domain in our simulations. Denoting by $\mathbf{n}$ the normal direction, it corresponds to imposing the following relation on the numerical boundary $\Gamma$:

$$ \text{absorbing boundary condition:} \qquad -\frac{1}{\rho(\mathbf{x}) c_{\uparrow}(\mathbf{x}, \omega)} p(\mathbf{x}, \omega) + \mathbf{v}(\mathbf{x}, \omega) \cdot \mathbf{n} = 0, \quad \text{on } \Gamma. \qquad (2.6) $$

**Wall boundary conditions** We also consider a different experiment, where the sample is enclosed in a domain with reflecting boundaries. Mathematically, this means that Neumann boundary conditions for the pressure have to be imposed:

$$ \text{wall boundary condition:} \qquad \nabla p(\mathbf{x}, \omega) \cdot \mathbf{n} = 0, \quad \text{on } \Gamma. \qquad (2.7) $$

This case is more challenging as multiple wave reflections occur from the boundary and appear in the data, as we highlight in the numerical experiments.

## 2.3. Attenuation with complex coefficients

The attenuation of waves is a frequency-dependent phenomenon, conveniently encoded in the frequency-domain equations with complex-valued parameters, cf. [7,62,11], while it can be more challenging to implement with the time-domain formulation, cf. [7,47]. For propagation with attenuation, the complex wavenumber $k$ is decomposed into its real and imaginary parts as follows,

$$ k(\mathbf{x}, \omega) = \frac{\omega}{c_{\uparrow}(\mathbf{x}, \omega)} = \frac{\omega_0}{c_p(\mathbf{x})} + i\alpha(\mathbf{x}, \omega), \qquad (2.8) $$

where $c_p : \Omega \to \mathbb{R}_{>0}(\Omega)$ is the phase velocity, [62], the imaginary part $\alpha : \Omega \to \mathbb{R}_{>0}(\Omega)$ and $\omega_0 \in \mathbb{R}_{>0}$.

The positiveness of the coefficients leads to conditions on the real and imaginary parts of the wave speed and frequency. Considering $c_{\uparrow} = c_R + ic_I$, we have

$$ \begin{align} \frac{\omega}{c_{\uparrow}} &= \frac{\omega_R + i\omega_I}{c_R + ic_I} = \frac{\omega_0}{c_p} + i\alpha \nonumber \\ \Rightarrow \quad & \frac{\omega_R c_R - i\omega_R c_I + i\omega_I c_R + \omega_I c_I}{c_R^2 + c_I^2} = \frac{\omega_0}{c_p} + i\alpha \tag{2.9} \\ \Rightarrow \quad & \frac{\omega_R c_R + \omega_I c_I}{c_R^2 + c_I^2} = \frac{\omega_0}{c_p} \quad \text{and} \quad \frac{\omega_I c_R - \omega_R c_I}{c_R^2 + c_I^2} = \alpha. \nonumber \end{align} $$

Therefore, to ensure that $\omega_0$, $c_p$ and $\alpha$ are positive, the real and imaginary parts of the wave speed and frequency must verify

$$ \omega_R c_R + \omega_I c_I > 0, \quad \omega_I c_R - \omega_R c_I \ge 0. \qquad (2.10) $$

Firstly, to verify the first condition in the case where $\omega_I = 0$ (i.e., for ordinary frequencies), we must have that the real part of the wave speed is positive: $c_R > 0$. Now assuming $c_R > 0$, the second condition in (2.10) is always valid (in particular when $\omega_I = 0$) when

$$ c_I \le 0. \qquad (2.11) $$

Therefore we consider that the wave speed has a non positive imaginary part. This is a minimal condition for attenuation models, see [1]. We summarize the main assumptions regarding attenuation below.

**Assumption 1 (Validity of wave propagation with attenuation).** We assume the following conditions for the formulation of the wave problem (2.3):

$$ \text{1. The complex frequency verifies:} \qquad \omega_R > 0 \quad \text{and} \quad \omega_I \ge 0; \qquad (2.12a) $$

$$ \text{2. The real part of the wave speed verifies:} \qquad c_R > 0; \qquad (2.12b) $$

$$ \text{3. The imaginary part of the wave speed verifies:} \qquad c_I \le 0; \qquad (2.12c) $$
In addition, we introduce the *quality factor* Q, that quantifies the level of attenuation.

**Definition 1 (Quality factor).** The *quality factor* Q is defined by the ratio between the real and imaginary parts of the complex wave speed ([11]),

$$Q(\mathbf{x}, \omega) := \frac{\mathrm{Re}(\kappa_{\uparrow}(\mathbf{x}, \omega))}{-\mathrm{Im}(\kappa_{\uparrow}(\mathbf{x}, \omega))} \stackrel{(2.4)}{=} \frac{\mathrm{Re}(c_{\uparrow}(\mathbf{x}, \omega)^2)}{-\mathrm{Im}(c_{\uparrow}(\mathbf{x}, \omega)^2)} > 0. \qquad (2.13)$$

The inverse of the *quality factor* $Q^{-1}$ is referred to as the *dissipation factor* and is 0 in a non-attenuating media.

## 2.4. Models of attenuation

The choice of *attenuation model* indicates the interplay between the real and imaginary parts of the wavenumber *k* and its frequency dependency. Several models have been introduced in the literature, cf., e.g., [62,11,18]. In this work, we are considering seven different models, which are given in Table 1. In their formulations, we indicate explicitly the dependency in frequency, such that all of the variables ($\kappa_0$, $\eta$, $\eta$, $\tau$, $\tau_\epsilon$ and $\tau_\sigma$) are real and only depend on the space coordinates. We further indicate the condition on the parameters to ensure the positiveness and negativeness of the signs, respectively for the real and imaginary parts of the wave speed, see Assumption 1.

We note that in the case where $\beta = 1$, the Cole–Cole model is equivalent to the Zener model. The Kolsky–Futterman model (2.15) originates from [38,29] and in our work, we use a simplified version which is frequency independent. It is obtained under the assumptions of *weak attenuation* ($4Q^2 \gg 1$), as detailed in [19, Section 1.7]) This model (in its frequency-dependent version) is also used in [56,45]. In the definition of the modified Szabo model (2.21), $\tau$ is placed together with the frequency and both are at power ($\beta - 1$), contrary to [18] where $\tau$ is separated from the power of $\omega$; This is only motivated for unit consistency, such that $\tau$ is expressed in time unit in (2.21).

**Remark 1 (Causality).** The models summarized in Table 1 are causal as they are directly derived from their time-domain counterpart, cf. [31,62,11,18]. More generally, causality principles can be verified via the Kramers–Kronig relations, or ensuring that the parameters are analytic in the lower-half complex plane, and we refer to [31,11,18] for more details.

**Remark 2 (Generalized models).** The attenuation models given in Table 1 account for one frequency mechanism of attenuation, which is sufficient when the data contain a relatively narrow band of frequency, as highlighted in our experiments. For broadband signals, one has to consider multiple mechanisms each associated with a different frequency, introducing a generalized model, cf. [12, 48, 11, 46, 67]. For instance, considering *L* attenuation mechanisms each associated with frequency $\omega_L$, one can define, e.g., [67],

<caption>**Table 1**</caption>

Attenuation models used for the visco-acoustic time-harmonic propagation, extracted from [62,18] and [11, Section 2]. The parameters $\kappa_0$, $\tau$, $\tau_\epsilon$, $\tau_\sigma$, $\eta$ and $\eta$ are real-valued, only depend on the space variable, and are non-negative to validate the conditions given in (2.12). In terms of units, $\kappa_0$ as unit [Pa], $\tau$, $\tau_\epsilon$ and $\tau_\sigma$ are in time units (in [s]), $\eta$ is in [Pa s] and $\eta$ is unitless.

<table><thead><tr><th>model</th><th>coeff.</th><th>complex bulk modulus</th><th>condition</th></tr></thead><tbody><tr><td>no-attenuation</td><td>&kappa;<sub>0</sub></td><td>&kappa;<sub>&uparrow;</sub><sup>(no)</sup> = &kappa;<sub>0</sub></td><td>(2.14) n/a</td></tr><tr><td>simplified Kolsky-Futterman</td><td>&kappa;<sub>0</sub>, &eta;</td><td>&kappa;<sub>&uparrow;</sub><sup>(kf)</sup> = &kappa;<sub>0</sub> - i &kappa;<sub>0</sub>/&eta;</td><td>(2.15) &eta; &ge; 0</td></tr><tr><td>Cole-Cole, [14,62]</td><td>&kappa;<sub>0</sub>, &tau;<sub>&epsilon;</sub>, &tau;<sub>&sigma;</sub>, &beta;</td><td>&kappa;<sub>&uparrow;</sub><sup>(cc)</sup> = &kappa;<sub>0</sub> (1 + (-i&omega;<sub>R</sub> &tau;<sub>&epsilon;</sub>)<sup>&beta;</sup> / (1 + (-i&omega;<sub>R</sub> &tau;<sub>&sigma;</sub>)<sup>&beta;</sup>))</td><td>(2.16) &tau;<sub>&epsilon;</sub> &ge; &tau;<sub>&sigma;</sub> &ge; 0, 0 &le; &beta; &le; 1</td></tr><tr><td>Zener, [69], [11, Section 2.4.3]</td><td>&kappa;<sub>0</sub>, &tau;<sub>&epsilon;</sub>, &tau;<sub>&sigma;</sub></td><td>&kappa;<sub>&uparrow;</sub><sup>(z)</sup> = &kappa;<sub>0</sub> (1 - i&omega;<sub>R</sub> &tau;<sub>&epsilon;</sub>) / (1 - i&omega;<sub>R</sub> &tau;<sub>&sigma;</sub>)</td><td>(2.17) &tau;<sub>&epsilon;</sub> &ge; &tau;<sub>&sigma;</sub> &ge; 0</td></tr><tr><td>Kelvin-Voigt, [11, Section 2.4.2]</td><td>&kappa;<sub>0</sub>, &tau;<sub>&epsilon;</sub></td><td>&kappa;<sub>&uparrow;</sub><sup>(kv)</sup> = &kappa;<sub>0</sub> - i&omega;<sub>R</sub> &kappa;<sub>0</sub> &tau;<sub>&epsilon;</sub></td><td>(2.18) &tau;<sub>&epsilon;</sub> &ge; 0</td></tr><tr><td>Maxwell, [11, Section 2.4.1]</td><td>&kappa;<sub>0</sub>, &eta;</td><td>&kappa;<sub>&uparrow;</sub><sup>(m)</sup> = &kappa;<sub>0</sub> &eta; / (1 - i&omega;<sub>R</sub> &eta; / &kappa;<sub>0</sub>)</td><td>(2.19) &eta; &gt; 0</td></tr><tr><td>KSB (Kowar-Scherzer-Bonnefond), [18,39]</td><td>&kappa;<sub>0</sub>, &eta;, &tau;, &beta;</td><td>&kappa;<sub>&uparrow;</sub><sup>(ksb)</sup> = &kappa;<sub>0</sub> / (1 + &eta;/(&radic;1 + (-i&omega;<sub>R</sub> &tau;)<sup>&beta;</sup>))<sup>2</sup></td><td>(2.20) &eta; &gt; 0, &tau; &gt; 0, 0 &lt; &beta; &lt; 1</td></tr><tr><td>modified Szabo, [60,18]</td><td>&kappa;<sub>0</sub>, &tau;, &beta;</td><td>&kappa;<sub>&uparrow;</sub><sup>(sz)</sup> = &kappa;<sub>0</sub> / (1 + (-i&omega;<sub>R</sub> &tau;)<sup>&beta;-1</sup>)</td><td>(2.21) &tau; &gt; 0, 0 &lt; &beta; &lt; 1</td></tr></tbody></table>
<caption>**Table 2**</caption>

Choice of coefficients for the attenuation models to compute $Q$ in Fig. 3. These are chosen such that the quality factor at 300 kHz is $Q=118$. For models that involve multiple coefficients, there exist different combinations that would also give $Q = 118$ at 300 kHz and we only choose one of them.

<table><thead><tr><th>attenuation model</th><th>parameters</th></tr></thead><tbody><tr><td>simplified Kolsky-Futterman</td><td>η = 118</td></tr><tr><td>Kelvin-Voigt</td><td>τ<sub>ε</sub> = 4.5 ns</td></tr><tr><td>Maxwell</td><td>η = 1.4.10<sup>5</sup> Pa s</td></tr><tr><td>Zener</td><td>τ<sub>ε</sub> = 90 ns, τ<sub>σ</sub> = 85.4 ns</td></tr><tr><td>Cole-Cole</td><td>τ<sub>ε</sub> = 90.5 ns, τ<sub>σ</sub> = 85.5 ns, β = 0.8</td></tr><tr><td>KSB</td><td>τ = 2.10<sup>5</sup> s, η = 8.75, β = 0.5</td></tr><tr><td>modified Szabo</td><td>τ = 13.28 s, β = 0.6</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Fig. 3.** Comparisons of attenuation models with frequency using $c_0 = 1500 \text{ m s}^{-1}$ and $ρ = 1000 \text{ kg m}^{-3}$. The attenuation parameters for each model are chosen such that their quality factor $Q$ is equal to 118 at frequency 300 kHz, cf. Table 2.</figcaption></figure>

$$ \kappa_{\uparrow}^{(G)}(\mathbf{x}, \omega) = \kappa_0(\mathbf{x}) \left( 1 - \sum_{l=1}^{L} B_l(\mathbf{x}, \omega) \frac{\omega_l}{\omega_l + i\omega} \right), \qquad (2.22) $$

where $B_l$ represents the attenuation function of the $l^{th}$ mechanism.

## 2.5. Comparison of attenuation models

We illustrate the dependency of the quality factor $Q$ defined in (2.13) with the frequency for the different attenuation models in Fig. 3. We fix the values of the density to $ρ = 1000 \text{ kg m}^{-3}$, and of the bulk modulus $κ_0 = 2.25 \text{ GPa}$, such that the wave speed is $c_0 = \sqrt{κ_0/ρ} = 1500 \text{ m s}^{-1}$. To select and compare the parameters of each of the attenuation models, we impose the quality factor at the reference frequency $ω_R/(2π) = 300 \text{ kHz}$ to be $Q = 118$, leading to the values represented in Table 2.

One main difference between the attenuation models is that either the low or high frequencies are attenuated the most. For instance, the quality factor decreases with frequency for the Kelvin-Voigt, Zener and Cole-Cole models, while it increases with frequency for the Maxwell, KSB and modified Szabo ones. For the latter family, it means that low-frequency waves are more attenuated while it is the high-frequency waves that are attenuated the most with the first family of models. In the investigated frequency band, we observe that the Kelvin-Voigt, Zener and Cole-Cole models have the same pattern with high quality factors at low frequencies, a rapid decrease, and then a stabilization. The KSB and modified Szabo models instead have low quality factors at low frequencies, before it slowly increases with higher frequencies. One could also use the generalized models, as discussed in Remark 2, to allow more flexibility on the quality factor's evolution. Here, we see that, while the quality factors of the attenuation models all coincide at frequency 300 kHz, they become all different when stepping away from this reference frequency.

## 3. Quantitative inverse problem using iterative minimization

The inverse wave problem associated to the forward equation (2.3) corresponds to the reconstruction of the heterogeneous properties $κ_{\uparrow}(\mathbf{x})$ and $ρ(\mathbf{x})$ given some wave measurements. In this section, we review the steps for quantitative
imaging based upon iterative minimization. In the context of seismic imaging where one uses the phase and amplitude information recorded by the seismograms this is usually referred to as the Full Waveform Inversion (FWI), [40,61,65,19].

## 3.1. Forward problem

We assume measurements of the pressure field $p$ at the position of receivers. As described in Subsection 2.1, we consider that the time-domain measurements results from a source, which is mathematically represented as a Ricker wavelet. Then, a complex Fourier (or Laplace-Fourier) transform is applied to work with frequency-domain data, see Fig. 2. Therefore, we define the forward problem $\mathcal{F}$ directly in the frequency domain: For a source $g$ (which corresponds to the complex Fourier transform of the original Ricker source), it is defined such that,

$$ \mathcal{F}(\mathbf{m}, \omega, g) = \{ p(\mathbf{x}_1, \omega, g), p(\mathbf{x}_2, \omega, g), \dots, p(\mathbf{x}_{n_{\text{rcv}}}, \omega, g) \}, \qquad (3.1) $$

where $p$ solves (2.3) with right-hand side $g$. Here, $\mathbf{x}_1, \dots, \mathbf{x}_{n_{\text{rcv}}}$ denotes the position of the $n_{\text{rcv}}$ receivers. The model parameters are represented by $\mathbf{m}$ and correspond to the bulk modulus and density such that $\mathbf{m} := \{\kappa_\uparrow, \rho\}$. We have $n_{\text{src}}$ independent sources in the experiment, and define the data-set associated to the frequency $\omega$ by

$$ \mathcal{F}(\mathbf{m}, \omega) = \{ \mathcal{F}(\mathbf{m}, \omega, g_1), \mathcal{F}(\mathbf{m}, \omega, g_2), \dots, \mathcal{F}(\mathbf{m}, \omega, g_{n_{\text{src}}}) \}. \qquad (3.2) $$

Here, we keep the notation $\mathcal{F}$ when there is no ambiguity. Note that we have assumed that the receivers remain in the same position for each source (that is, they do not depend on $g$), but modifying the receiver positions together with the source would not modify the methodology.

## 3.2. Reconstruction with iterative minimization, FWI

The reconstruction is carried out following the minimization of a misfit functional $\mathcal{J}$, that evaluates a distance between measurements $\mathbf{d}$ and simulations:

$$ \mathcal{J}(\mathbf{m}, \omega) = \text{dist}(\mathcal{F}(\mathbf{m}, \omega), \mathbf{d}(\omega)), \qquad (3.3) $$

where the choice of misfit function, depending on the data-sets and inverted parameters, is the subject of numerous studies, e.g., [58,28,9,63,35,21,26,23]. Per simplicity, we rely on the $l^2$-distance of the difference, such that

$$ \mathcal{J}(\mathbf{m}, \omega) = \frac{1}{2} \sum_{k=1}^{n_{\text{src}}} \| \mathcal{F}(\mathbf{m}, \omega, g_k) - \mathbf{d}(\omega, g_k) \|_2^2. \qquad (3.4) $$

The reconstruction is performed by minimizing $\mathcal{J}$ with respect to the model parameters $\mathbf{m}$, with a truncated-Newton method, [49]. We further follow a progressive increase in the frequency-content to improve the convergence of the algorithm and mitigate ill-posedness, as advocated in [10,4,22].

The pressure wave propagating through the medium is broadband (see Fig. 1d). It is recorded at a finite number of receivers. Then the complex Fourier transform of the time signals is performed, and we get information for all frequencies $\omega = \omega_R + i\omega_I$ such that $\omega_R^{(\text{min})} \le \omega_R \le \omega_R^{(\text{max})}$ and $\omega_I^{(\text{min})} \le \omega_I \le \omega_I^{(\text{max})}$. Following the guidelines of [4,22], we further use a sequential frequency progression (i.e., one $\omega$ at a time) in the reconstruction algorithm (instead of a band of frequencies inverted at once). The progression in the selection of $\omega_R$ and $\omega_I$ is detailed and motivated below. The iterative minimization algorithm is summarized in Algorithm 1.

**Ordinary frequency progression** When the imaginary part of the frequency is zero, $\omega_I = 0$, (see Fig. 2), one works with ordinary frequencies only. In this case the progression in frequency content is chosen from low to high values, [10,4,22,19, 65]. That is, we start with $\omega^{(1)} = \omega_R^{(1)} \ge \omega_R^{(\text{min})}$, carry out minimization iterations, and then update to $\omega^{(2)} = \omega_R^{(2)} > \omega_R^{(1)}$, and repeat with the next frequency, see Algorithm 1.

**Complex frequency progression** In the case where the imaginary part of the frequency is non-zero (illustrated in Fig. 2), the progression in the frequency content follows the strategy of [22]: $\omega_I$ varies first, from high to low values with $\omega_R$ fixed. Then $\omega_R$ is increased and we repeat the sequence of $\omega_I$ (from high to low). This choice is motivated from the estimates of attraction basins size provided in [22]. Namely, we start with $\omega^{(1)} = \omega_R^{(1)} + i\omega_I^{(1)}$, such that $\omega_R^{(1)} \ge \omega_R^{(\text{min})}$ and $\omega_I^{(1)} \le \omega_I^{(\text{max})}$. The next frequency is $\omega^{(2)} = \omega_R^{(2)} + i\omega_I^{(2)}$, with $\omega_I^{(2)} \le \omega_I^{(1)}$, see Algorithm 1.
**Algorithm 1:** Iterative minimization procedure for quantitative reconstruction using complex frequencies. The real and imaginary parts of the frequency are ordered such that $\omega_l^{(1)} > \omega_l^{(2)} > \dots > \omega_l^{(N_{\omega_l})}$ and $\omega_R^{(1)} < \omega_R^{(2)} < \dots < \omega_R^{(N_{\omega_R})}$.

Initialization: starting model parameters $\mathbf{m}_1 = (\kappa_{\dagger,1}, \rho_1)$, and measurements $\mathbf{d}$.

for $i = 1, \dots, (N_{\omega_l})$ do

for $j = 1, \dots, (N_{\omega_R})$ do

Set $\omega := \omega_R^{(j)} + i\omega_l^{(i)}$;

Compute the complex Fourier transform of the time-domain measurement at $\omega$, see (2.2) and Fig. 2;

for $k = 1, \dots, n_{\text{iter}}$ do

Set $l := (i - 1) N_{\omega_R} n_{\text{iter}} + (j - 1)n_{\text{iter}} + k$;

Solve Problem (2.3) using current models $\mathbf{m}_l$ and frequency $\omega$ for all sources in the acquisition;

Compute the misfit function $\mathcal{J}(\mathbf{m}_l, \omega)$;

Compute the gradient of the misfit function with adjoint-state method, [25];

Compute the search direction (e.g., nonlinear conjugate gradient, [49]);

Compute the step length, $\alpha_l$, using line-search method, [49];

Update the model: $\mathbf{m}_{l+1} = \mathbf{m}_l - \alpha_l s_l$.

end

end

end

### 3.3. Numerical implementation

For the discretization of the forward wave problem (2.3), we use the Hybridizable Discontinuous Galerkin method (HDG), [13,25]. The HDG method used the mixed formulation and static condensation to solve Problem (2.3). It results in a linear system to be solved, which size is shown to be less than, e.g., for finite elements method, depending on the degree of the approximation polynomials, see, [36,66,25]. Our implementation for the forward and inverse problems follows the steps described in [25], where the only change is to consider a complex-valued bulk modulus $\kappa_{\dagger}$. The computation of the gradient is performed with the adjoint-state method ([54,52,4]) specifically derived for the HDG discretization in [25]. For the two and three dimensional computational experiments that are carried out in the following sections, we use the open-source software *hawen*,<fnref n="1" /> [20]. The precise computational cost is further discussed in Subsection 5.4.

## 4. Numerical experiments in two dimensions

In this section, we investigate a two-dimensional experiment, following the setup described in Fig. 2. The reconstruction is carried out following Algorithm 1, and we study the following:

1. We evaluate the robustness of the reconstruction procedure with attenuation model uncertainty, that is, we select an attenuation model to generate synthetic measurement data and use a different model to carry out the reconstruction.

2. While the primary focus of this paper is on the reconstruction of the wave speed variations, we also investigate the choice of parametrization to recover the medium's density, which is known to be harder to recover, see [65,33,19] and the references therein.

3. We investigate the effect of the domain boundary conditions, that is, whether we allow waves to freely escape the domain using absorbing boundary conditions (free-space) or assuming wall boundaries around the sample.

The following experiments have all been realized using software *hawen*, [20], see Footnote 1. We further detail the computational cost of the 2D and 3D reconstructions in Subsection 5.4.

### 4.1. Experimental setup

We consider a synthetic experiment of ultrasound imaging with a two-dimensional breast model, described by its wave speed $c_0$, density $\rho$, and quality factor $Q$, pictured in Fig. 4. The medium is of size $18 \times 18$ cm<fnref n="2" /> and is composed of tissues, blood and fat. These layers are obtained from a cross-section of the OA-Breast Phantom data-set,<fnref n="2" /> [43], where the values of the parameters depending on the layers (e.g., blood, skin) are taken following the IT'IS Database,<fnref n="3" /> see Table 3. In addition, we incorporate a high-contrast ellipsoid inclusion in the model parameters, which can be seen as a defect that has to be identified with ultrasound imaging. From the values given in Table 3, we see that the quality factor is the lowest in the skin layer, that is, waves are attenuated the most by skin. In the other layers, the quality factor is relatively high (at least 280), such that we have a weakly attenuating medium, [64].

To carry out the reconstruction, we assume measurements of the pressure field, following the setup described in Subsection 2.1: A Ricker source wavelet excites the medium and time-domain data are measured. The Ricker wavelet conveys

<footnote n="1">¹ https://ffaucher.gitlab.io/hawen-website/</footnote>

<footnote n="2">² https://anastasiolab.wustl.edu/downloadable-content/oa-breast-database/</footnote>

<footnote n="3">³ https://itis.swiss/virtual-population/tissue-properties/</footnote>
<caption>**Table 3**</caption>

Interval of the model values depending on the type of tissues. These follow the IT'S database, see Footnote 3.

<table><thead><tr><th>medium</th><th>wave speed c<sub>0</sub> (m s<sup>-1</sup>)</th><th>density ρ (kg m<sup>-3</sup>)</th><th>Q at 300 kHz</th></tr></thead><tbody><tr><td>background</td><td>1490</td><td>1000</td><td>800</td></tr><tr><td>skin</td><td>(1590, 1610)</td><td>(1100, 1120)</td><td>(100, 120)</td></tr><tr><td>blood</td><td>(1565, 1575)</td><td>(1090, 1110)</td><td>(290, 310)</td></tr><tr><td>fat</td><td>(1440, 1460)</td><td>(920, 940)</td><td>(410, 430)</td></tr><tr><td>glandular tissue</td><td>(1490, 1520)</td><td>(1030, 1050)</td><td>(280, 300)</td></tr><tr><td>inclusion</td><td>1550</td><td>1050</td><td>350</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Fig. 4.** Two-dimensional breast models of wave-speed, density and quality factor at frequency $\omega_R/(2\pi) = 300$ kHz including a high-contrast inclusion. The values of the parameters depending on the type of tissues are detailed in Table 3. The domain is of size $18 \times 18$ cm$^2$.</figcaption></figure>

<caption>**Table 4**</caption>

Interval of values for the coefficients of the attenuation models to represent the parameters given Fig. 4 and Table 3.

<table><thead><tr><th>model</th><th>parameters</th></tr></thead><tbody><tr><td>simplified Kolsky-Futterman</td><td>η ∈ (100, 800)</td></tr><tr><td>Kelvin-Voigt</td><td>τ<sub>ε</sub> ∈ (0.663, 5.305) ns</td></tr><tr><td>Maxwell</td><td>η ∈ (0.148, 0.942) MPa s</td></tr><tr><td>Zener</td><td>τ<sub>ε</sub> ∈ (4.729, 8.663) ns, τ<sub>σ</sub> ∈ (1, 8) ns</td></tr><tr><td>Cole-Cole</td><td>τ<sub>ε</sub> ∈ (3.150, 8.387) ns, τ<sub>σ</sub> ∈ (1, 8) ns, β = 0.8.</td></tr><tr><td>KSB</td><td>τ ∈ (1, 8) s, η ∈ (17.332, 3.036.10<sup>8</sup>), β = 0.5</td></tr><tr><td>modified Szabo</td><td>τ ∈ (5.905, 384.552) s, β = 0.5</td></tr></tbody></table>

multiple frequency contents (Fig. 1d) and one transforms the recorded time-domain data into complex frequency-domain ones, see Fig. 2. Nonetheless, as we consider a synthetic experiment, for simplicity we directly generate the data-sets in the frequency domain and incorporate white noise a-posteriori. The acquisition is composed of 36 independent point-sources positioned in a circular pattern, as illustrated in Fig. 1a. The receivers measuring the pressure field for each of the source are positioned onto the same disk, with a total of 360 receivers.

For the reconstruction procedure, the initial models consist in the constant values of the parameters in the background, that is, the inversion starts with $c_0 = 1500$ m s$^{-1}$, $\rho = 1000$ kg m$^{-3}$ and attenuation parameters chosen (depending on the attenuation model) such that $Q = 800$ at $\omega_R/(2\pi) = 300$ kHz.

## 4.2. Comparison of wave propagation with attenuation models and boundary conditions

We generate synthetic data for each of the seven attenuation models given in Table 1, leading to seven different data-sets. The choice of the parameters that describe the attenuation is such that the quality factors Q of each model are equal at 300 kHz, and can be seen in Fig. 4c. This amounts to the intervals of values given in Table 4.

For the models that have several degrees of freedom, we consider different magnitudes in the parameters. We illustrate the data depending on the medium attenuation model in Fig. 5, for a source positioned at angle 0°(see Fig. 2) and the medium properties pictured in Fig. 4. We further compare the simulations with absorbing or wall boundary conditions. For the sake of clarity, we only compare three ordinary (real) frequencies, that is, we keep $\omega_I = 0$, and picture the results for the Kolsky-Futterman, Kelvin-Voigt, and KSB models. The data used for inversion are obtained every 1°(i.e., 360 receivers, see Fig. 5e), but we also plot in Fig. 5 the signals with data-point every 0.1° for better visualization and comparisons.

We observe that the three selected attenuation models all give different signals. In Fig. 5a, we see that the difference exists for all angles $\theta$, and particularly at large offset, with drastic changes in terms of amplitude. In Figs. 5b and 5c, we picture the solution at 500 kHz for $\theta \in (90^\circ, 120^\circ)$ depending on the choice of the boundary condition: We observe that when using absorbing conditions (i.e., assuming free-space propagation, Fig. 5c), the difference between the signals is small
<figure><img src="image_12.png" /><figcaption>a) Signals at 100 kHz using wall boundary conditions on the sides, with one data point per 0.1°.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>b) Signals at 500 kHz using wall boundary conditions on the sides, with one data point per 0.1°.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>c) Signals at 500 kHz using absorbing boundary conditions on the sides, with one data point per 0.1°.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>d) Signals at 300 kHz using wall boundary conditions on the sides, with one data point per 0.1°.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>e) Signals at 300 kHz using wall boundary conditions on the sides, with one data point per 1°.</figcaption></figure>

**Fig. 5.** Comparison of the pressure field at the position of the receivers depending on the attenuation model of the medium, with medium parameters pictured in Fig. 4. The receiver devices are positioned around the sample, in the configuration of Fig. 2, and the source is at angle $θ=0°$. For visualization, we first consider one receiver every 0.1°. Nevertheless, the data used for inversion consider one receiver every 1° instead, i.e., 360 receivers.

both in terms of amplitude and phase. On the other hand, when assuming the wall boundaries on the sides of the domain, we see drastic changes in the wavefields, cf. Fig. 5c. In Fig. 5d, we plot the results for frequency 300 kHz, at which the quality factors of each attenuation model are the same, here we see that the signals corresponding to the Kolsky-Futterman and Kelvin-Voigt models are indeed very close. Nonetheless, the signal corresponding to the KSB model is different in amplitude and phase. Furthermore, to carry out the inversion, we see that we have a sub-sampled signal with one data-point every 1°, as plotted in Fig. 5e.

Note also that there is some flexibility with the position of the absorbing boundary conditions to constrain the computational domain, as waves are anyway supposed to freely propagate up until infinity. On the other hand, when considering wall boundaries, it is essential that the computational domain boundaries follow the actual setup, so that the reflections from the wall are adequately simulated.

### 4.3. Reconstruction with absorbing boundary conditions (free-space)

We consider the medium in free-space, such that acoustic waves propagate according to (2.3) in the domain $Ω$ with absorbing boundary conditions (2.6) on $Γ$. In this configuration, one assumes that either the propagation domain is sufficiently large so that reflections from the boundaries are negligible, or that these reflections can somehow be removed in a pre-processing stage of the data. Note that the later option is *not* a trivial task in general (cf. the suppression of the so-called multiples in seismic acquisition) and can lead to inconsistency in the data. We further investigate the case of reflecting boundaries, which is more challenging, in Subsection 4.4.

#### 4.3.1. Wave speed reconstruction with attenuation model uncertainty

In this test-case, we investigate the robustness of the iterative reconstruction procedure with respect to the attenuation model. We generate a synthetic data-set for each of the attenuation model (see Subsection 4.2) and carry out the recon-
struction using a different one. All of the data-sets have the same acquisition with 36 independent point-sources and 360 receivers measuring the pressure field, as illustrated in Figs. 2 and 5. We use synthetic data, but the numerical setup is changed between the forward and inverse problems. Namely, the data are generated using a discretization mesh of more than 200000 cells to ensure the accuracy of the different layers, and the solution is approximated with polynomials of order 5. On the other hand, the inversion is carried out on a mesh of less than 10000 cells, with variable orders using *p*-adaptivity for efficiency, cf. [25]. In addition, white noise is added to the synthetic data, using a signal-to-noise ratio of 20 dB. During inversion, the model parameter is represented as a piecewise-polynomial function per cell, using Lagrange basis function of linear order. Therefore, the model parameters are represented by 3 unknowns on each cell of the discretized mesh.

In this experiment we only use ordinary (real) frequencies (i.e., regular Fourier transform of the time-domain data) such that $\omega_l = 0$. We consider the increasing frequency progression as follows: $\omega_R/(2\pi) = \{200, 300, 400, 600\}$ kHz. The evolution in the frequency content of the data is chosen sequentially, [22], and 30 iterations are performed per frequency for the minimization, following Algorithm 1. Therefore, there is a total of 120 iterations. The reconstruction is carried out with respect to the bulk modulus $\kappa_0$, that is, the density $\rho$ and quality factor $Q$ remain at their constant initial values. For visualization, we instead picture the resulting wave speed $c_0 = \sqrt{\kappa_0/\rho}$, assembled from the reconstructed real part of the bulk modulus $\kappa_0$ and the (homogeneous) initial density $\rho$. In Fig. 6 are shown the reconstructed wave speed models, comparing for all data-sets the use of a different model for inversion where, for visualization, the raw results are slightly smoothed a-posteriori, see Remark 3.

While each attenuation model gives a quite different data-set, as highlighted in Fig. 5, we see that all of the configurations provide the appropriate layers in the sample, and have similar accuracy. The main tissue features are reconstructed with the appropriate values and the skin of the sample is correctly obtained. Namely, using an attenuation model for the reconstruction that is different from the one of the original medium *does not* alter the accuracy of the features reconstructed, demonstrating the robustness of the iterative minimization procedure for this two-dimensional test-case. In addition, we see that while the density $\rho$ and attenuation $Q$ are not inverted and remain to their initial values, the reconstructed profiles of wave speed $c_0$ are accurate.

In all of the cases, the ellipse-shaped included defect is found at the right position and with its precise shape. We note that the speed in the skin layer in slightly less than expected, while the one in the contrasting object is slightly higher than expected. The robustness of the reconstruction algorithm with respect to changes in attenuation models could be explained by two reasons:

1. Data of relatively limited frequency bandwidth are sufficient for the reconstruction (from 200 to 600 kHz) and the difference between attenuation models is less significant is such a narrow band, see Fig. 3.

2. The reconstruction with FWI is known to be more sensitive to phase shifts rather than amplitude, hence robust with respect to incorrect attenuation which mostly affects the signal amplitude. Indeed, the dispersion resulting from the change of attenuation models remains small compared to the one that comes from changes in the bulk modulus $\kappa_0$.

In addition, we are in a case of tissues which have weak attenuation properties, which may also justify the lack of sensitivity to change in attenuation model, even though the wave fields of Fig. 5 show strong differences. In the following experiments, for the sake of clarity, we only consider one of the cases, with the data-set generated with the Kolsky-Futterman model, while FWI is conducted using the Kelvin-Voigt model.

**Remark 3 (Post-processing visualization).** The visualization of the reconstruction may suffer from the coarse discretized mesh (about 10000 cells) employed for the numerical discretization, which is chosen to reduce the numerical cost. To enhance the visualization, we use the function `imgaussfilt` of MATLAB, which corresponds to a Gaussian smoothing filter, that we illustrate in Fig. 7. This procedure is done a posteriori, independently of the reconstruction algorithm, and is therefore costless. One could instead rely on an extra regularization term in the optimization, such that Tikhonov or Total Variation regularization, but it needs an additional threshold in the formulation of the minimization problem, which can be difficult to select, e.g., [26]. This post-processing is only performed to enhance the visual aspect of the reconstructions by smoothing some numerical artifacts.

### 4.3.2. Multi-parameter reconstruction

We have inverted with respect to model parameter $\kappa_0$ in the previous experiments, and now investigate the reconstruction of two model parameters simultaneously. We refer to the choices of the parameters to invert as the *parametrization*, which is shown to have a strong influence on the reconstruction, cf. [8,37,19]. Namely, the physical properties characterizing the medium can be expressed with different physical parameters, for instance, omitting the attenuation property, we have the bulk modulus $\kappa_0$, the density $\rho$, the wave speed $c_0 = \sqrt{\kappa_0/\rho}$, or the impedance $I_0 = \sqrt{\kappa_0/\rho}$. For inversion, here we select two parameters among the aforementioned four, that will be reconstructed simultaneously, following Algorithm 1. It amounts to twelve possibilities but, for the sake of conciseness, we consider four cases:

- Inversion with respect to $(\kappa_0, \rho)$, - Inversion with respect to $(I_0, \rho)$,

- Inversion with respect to $(c_0, \rho)$, - Inversion with respect to $(I_0, c_0)$.

In particular, the last choice is motivated by [65]. To evaluate the gradient of the misfit function with respect to the selected model parameter, we start by computing the derivatives with respect to $\kappa_0$ and $\rho$ (which are the main unknowns
<figure><img src="image_17.png" /><figcaption>**Fig. 6.** Reconstruction of the wave speed $c_0 = \sqrt{\kappa_0/\rho}$ using the iterative minimization Algorithm 1 with sequential frequency progression $\omega_R/(2\pi) = (200, 300, 400, 600)$ kHz and fixed $\omega_l = 0$, starting from homogeneous backgrounds (Table 3). The visualization corresponds to $x$ from 2 to 16 cm and $z$ from 4 to 14.5 cm, with color scale similar to Fig. 4. Each line corresponds to a different data-set, each column corresponds to a different attenuation model used for FWI: (column 1) Kolsky–Futterman model; (column 2) Kelvin–Voigt model; (column 3) Maxwell model; (column 4) Zener model; (column 5) Cole–Cole model; (column 6) KSB model; (column 7) Szabo model. Diagonal elements are reconstructions with the “right” model.</figcaption></figure>

appearing in the wave equation (2.3)), and then use the chain rule formula, cf. [19, Sections 5.5 and 5.6]. In Fig. 8, we compare the reconstructions: For all cases, we picture the wave speed $c_0$ and density $\rho$ that are reassembled from the reconstructed model parameters. In addition, we maintain the attenuation model uncertainty, with data generated with the Kolsky–Futterman attenuation model, while inversion is carried out with the Kelvin–Voigt model.

The assembled wave speed is similarly accurate with all choices of parametrizations, with the different tissue features correctly retrieved and the ellipse-shaped defect resulting in a contrast of wave speed is visually clearly identified. However, we observe major differences in the reconstruction of the density: The parametrization $(\kappa_0, \rho)$, $(l_0, c_0)$ and $(l_0, \rho)$ are able to identify the main layers and $(l_0, \rho)$ provides the best resolution, with smoother contrasts. Nonetheless, the values of the density are incorrect in all cases, and remain near to the initial value. Here, the parametrization $(c_0, \rho)$ is the worst regarding the density reconstruction, with a noisy background where only the skin delimitation barely appears.
<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 7.** Illustration of the Gaussian filtering of the reconstruction to a-posteriori improve the visualization. We use a standard deviation for the Gaussian $\sigma = 2$. This reconstruction corresponds to the upper left image of Fig. 6.</figcaption></figure>

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

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

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

<figure><img src="image_23.png" /><figcaption>**Fig. 8.** Comparison of the reconstruction depending on the parametrization of the inverse problem. The reconstruction follows Algorithm 1 with sequential frequency progression $\omega_R/(2\pi) = (200, 300, 400, 600)$ kHz and fixed $\omega_I = 0$ with 30 iterations per frequency, starting from homogeneous backgrounds (Table 3). The domain units are given in cm.</figcaption></figure>

**Remark 4.** We note that the difficulty of recovering the density is similar in seismic imaging, [65,19]. This is explained as density perturbations mostly modify the amplitude of the signals and, as indicated above, the FWI focuses on the differences in the signals phase-shift. As an alternative, one could carry the FWI reconstruction of one main parameter, such as the bulk modulus, wave speed or impedance, and then infer the density via analytic formula inherited from physics, such as Gardner's relation, [30]. Another alternative is to select a different misfit function.

### 4.4. Reconstruction with wall boundary conditions and complex frequencies

In the previous experiments, we have seen that the wave speed is accurately reconstructed, even in the case of an attenuation model error. We have assumed that waves can escape the medium without incurring reflections, by using in the numerical simulations absorbing boundary conditions on all sides. We now instead consider wall boundary conditions on the domain boundary $\Gamma$ with (2.7), hence leading to multiple wave reflections from each side. The background coefficients also have low attenuation (see Table 3), further supporting these multiple reflections. We shall see that the modification of the boundary condition drastically alters the performance of the reconstruction, and that complex frequencies can be used to overcome the difficulties. For the sake of conciseness, we only consider the data-set using the Kolsky–Futterman attenuation model, while the reconstruction procedure is carried out with the Kelvin–Voigt model.

#### 4.4.1. Reconstruction with ordinary frequencies ($\omega_R > 0, \omega_I = 0$)

We keep the setup of the previous experiment where we only use Fourier frequencies with the progression $\omega_R/(2\pi) = (200, 300, 400, 600)$ kHz, with a fixed $\omega_I = 0$. We invert with respect to the bulk modulus only, and picture the reconstruction in Fig. 9, where we compare different choices of fixed quality factor for the reconstruction, that is, different choices of initial values for the attenuation parameters.
<figure><img src="image_24.png" /><figcaption>**Fig. 9.** Reconstruction of the wave speed $c_0 = \sqrt{\kappa_0/\rho}$ in the case of wall boundary conditions with different choices of (fixed) attenuation parameters. The reconstruction follows Algorithm 1 with sequential frequency progression $\omega_R/(2\pi) = \{200, 300, 400, 600\}$ kHz and fixed $\omega_I = 0$ with 30 iterations per frequency. The reconstruction starts from homogeneous wave speed and density (Table 3), and we compare different choices of constant attenuation parameters to observe the effect of the initial guess for the attenuation model onto the reconstruction.</figcaption></figure>

We see that the consideration of wall boundary on the sides of the medium completely changes the behavior and prevents any accurate reconstruction when using the constant value of background attenuation as initial guess, Fig. 9a, while it is the “natural” choice of background value. It appears that the multiple reflections from the boundaries totally annihilate the possibility of reconstruction, with information from the sample merged within the wave reflections. Interestingly, when the initial guess for the quality factor is much lower than expected (i.e., the initial guess is too attenuating compared to the real medium), the FWI reconstruction is slightly better, with the tissue features appearing, see Figs. 9b and 9c. However, the reconstruction remains with a poor resolution compared to the previous experiment assuming free-space propagation.

This experiment highlights the difficulty of having strong multiple reflections from the medium boundary, drastically reducing the performance of the reconstruction procedure. It would be necessary to treat those multiples in the measurement data, [68], to clear the reflections coming from the domain boundaries, but this task is not trivial in practice. As an alternative, one could perform the experiments in a medium with strong attenuation, to reduce or prevent the multiple reflections. In the next section, we show how to resolve the issue of the multiple boundary reflections, by encoding an artificial damping via the use of complex frequencies.

### 4.4.2. Reconstruction using complex frequencies ($\omega_R > 0, \omega_I \ge 0$)

We now investigate the use of complex frequencies for the case with wall boundary conditions, where we keep the same set of Fourier frequency, enriched with a Laplace component. In this case, it corresponds to applying a complex Fourier transform to the original time-domain observed data, as illustrated in Fig. 2.

For the reconstruction, we use the following sets:

- Fourier frequency $\omega_R/(2\pi) : \{200, 300, 400, 600\}$ kHz,

- Laplace damping $\omega_I : \{20.10^3, 15.10^3, 10.10^3\}$ s$^{-1}$.

Following the analysis of the inverse problem for complex frequencies carried out in [22], the progression of content first varies the Laplace damping from high to low for a fixed Fourier frequency, which then varies from low to high. Namely, the inverted complex frequencies follow the order: (200 kHz, $20.10^3$), (200 kHz, $15.10^3$), (200 kHz, $10.10^3$), (300 kHz, $20.10^3$), (300 kHz, $15.10^3$), etc. The reconstructed wave speed $c_0$ using the complex frequency set is pictured in Fig. 10.

The reconstruction using complex frequencies is accurate, with all the layers and the contrasting tumor appearing with the same high resolution as in the case of absorbing boundary conditions. Therefore the artificial damping has allowed to overcome the multiple boundary reflections. Note that this procedure can be applied for any data measured in the time-domain, and only necessitates to apply the complex Fourier transform as illustrated in Fig. 2. The performance of the use of the Laplace-Fourier transform with several complex frequencies can further be explained as the transform introduces artificial damping with the complex frequency, hence enhancing the first arrivals ([34]) by reducing the multiple reflections.

## 5. Numerical experiments in three dimensions

### 5.1. Experimental setup

We set up a three-dimensional experiment with a breast model extracted from the OA-Breast Phantom data-set, see Footnote 2, and further include a contrasting 3D object embedded in the medium. In Fig. 11 are displayed the wave speed, density and quality factor, which are encompassed in a rectangular domain of size $10 \times 10 \times 16$ cm$^3$. The values of the parameters in the different layers of the sample follow the values of the previous experiment, and are prescribed in Table 3.

The data are generated from the sides of this bounding box, except for the plane where $x = 10$ cm, which correspond to the patient's chest in Fig. 11. That is, the data are generated from the five remaining sides. We consider 224 point-sources
<figure><img src="image_25.png" /><figcaption>**Fig. 10.** Reconstruction of medium with wall boundary condition using complex frequencies. The starting models correspond to the constant background values of wave speed, density and quality factor of Table 3. (i.e., using fixed attenuation parameters chosen such that $Q = 800$ at 300 kHz). The corresponding case without using the complex frequencies is pictured Fig. 9a.</figcaption></figure>

and 1533 receivers equally partitioned between the different sides, with the receivers measuring the pressure field for each of the sources. Note that the choice of using a rectangular acquisition setup is only motivated per simplicity, and we could also employ a spherical acquisition, only requiring to adapt the discretization mesh.

For the reconstruction, we start with constant background parameters with $c_0 = 1490 \text{ m s}^{-1}$, $\rho = 1000 \text{ kg m}^{-3}$ and constant attenuation parameters chosen such that $Q (300 \text{ kHz}) = 800$. Furthermore, we only invert for the bulk modulus and keep the density and quality factor at their initial values during the entire iterations of FWI. We also incorporate attenuation model error in the experiment: the data are generated using the simplified Kolsky-Futterman model while the inversion is carried out with wave propagation using the Kelvin-Voigt attenuation model.

## 5.2. Reconstruction with absorbing boundary conditions (free-space)

In the free-space problem, absorbing boundary conditions (2.6) are implemented on all sides of surrounding box. In this case, following the observations of Section 4, we expect that only the Fourier frequencies would be needed (i.e., taking $\omega_l = 0$ to transform time-domain signals in Fig. 2), and we select the following three frequencies $\omega_R/(2\pi)$ for the reconstruction: $\{100, 200, 300\}$ kHz. We perform 20 minimization iterations per frequency, hence a total of 60 iterations. The wave-speed (assembled from the reconstructed bulk modulus) is shown in Fig. 12.

We see that the main features of the sample appear correctly on the reconstruction, with the appropriate values of speed recovered. The contrasting shape appears appropriately, in its correct location, as seen on the different cross-sections of Fig. 12. We also observe some background oscillatory artifacts in the surroundings of the medium but the skin contour of the breast is well identified. Here, smaller details and tissue features are missing and may require higher frequency-contents in the data (if available), to refine the accuracy of the reconstructed model.

## 5.3. Reconstruction with wall boundary conditions and complex frequencies

We consider wall boundary conditions (2.7) on the sides of the acquisition box. Here, we maintain an absorbing condition on the side of the chest (the plane at fixed $x = 10 \text{ cm}$ in Fig. 11) and use wall conditions on all of the other sides. To alleviate the difficulties of handling the multiple reflections generated from the different sides, we use complex frequencies. The Fourier content remains the same as the one used above, with $\omega_R/(2\pi) = \{100, 200, 300\}$ kHz, and we further incorporate two damping coefficients: $\omega_l = \{10^4, 5.10^3\} \text{ s}^{-1}$. The sequential progression of frequency follows a low-to-high Fourier content and high-to-low Laplace one, [22], such that we use $(\omega_R/(2\pi), \omega_l) = \{(100 \text{ kHz}, 10^4), (100 \text{ kHz}, 5.10^3), (200 \text{ kHz}, 10^4), \dots\}$.

The reconstructed wave speed is pictured in Fig. 13, where we see that we obtain an accuracy relatively similar to the case assuming the free-space propagation (Fig. 12). The main tissue features are visible in the cross-sections and the skin layer is recovered. In addition, the position, size and amplitude of the contrast are obtained and clearly stand out on the reconstruction. Similar to the free-space medium reconstruction, we observe some background perturbations, which however does not prevent us from clearly identifying the main components of the sample.

## 5.4. Computational cost

The reconstruction relies on the iterative algorithm depicted in Algorithm 1, where the computational cost comes from the resolution of the forward problems which has to be repeated with iterations. For the resolution of the time-harmonic problem, the discretization leads to relatively large matrix, and one further needs to solve the resulting linear system for each of the source in the acquisition. In our implementation, we made the following choices regarding the computational implementation:
<figure><img src="image_26.png" /><figcaption>**Fig. 11.** 3D breast sample in a bounding domain of size 10 × 10 × 16 cm³. The values of the parameters depending on the type of tissues are given in Table 3.</figcaption></figure>

- We use HDG discretization method, see Subsection 3.3, which can help reduce the computational cost compared to other approaches, see [25].

- Using HDG, the order of the approximation polynomial on each cell of the discretization mesh is chosen independently depending on the local (to the cell) wavelength, to ensure the linear system is as small as possible, cf. [25].

- For the resolution of the linear systems, we use the direct solver MUMPS, [2], that is particularly efficient for solving systems with multiple right-hand sides, contrary to iterative solvers. Namely, once the global matrix is factorized, the resolution for all of the sources (i.e., 224 point-sources in our 3D experiments) is fast. On the other hand, the matrix factorization is a computationally intensive step, especially regarding the memory consumption.

We use software *hawen* (see [20] and Footnote 1) in all of the experiments. They are carried out on the cluster PlaFRIM.<fnref n="4" />
The computational costs are the following:

<footnote n="4">⁴ https://www.plafrim.fr/.</footnote>
<figure><img src="image_27.png" /><figcaption>**Fig. 12.** Three-dimensional reconstruction of the breast model Fig. 11 assuming free-space propagation. The reconstruction starts from constant background parameters. The iterative minimization uses three frequencies from 100 to 300 kHz, with a total of 60 iterations.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Fig. 13.** Three-dimensional reconstruction of the breast model Fig. 11 assuming wall boundary on the sides of the acquisition box. The reconstruction starts from constant background parameters and only the bulk modulus is inverted while the density and quality factor remain constant. The iterative minimization uses six (complex) frequencies with $\omega_R/(2\pi) = (100, 200, 300)$ kHz and $\omega_l = (10^4, 5.10^3)$ s$^{-1}$, for a total of 120 iterations.</figcaption></figure>

- The 2D experiments of Section 4 are performed on 36 physical cores, with 18 processors and 2 threads per processors. The computational time to obtain a reconstruction as in Fig. 6 is about 10 min and corresponds to 120 inversion iterations. The memory required to factorize the matrix is about 1 GiB for all frequencies.

- The 3D experiments are performed on 540 physical cores, with 90 processors and 6 threads per processors. The computational time to obtain a reconstruction as in Fig. 12 is of 5 h and corresponds to 60 inversion iterations. The memory required to factorize the matrix depends on the frequency (high frequency needs higher polynomial order hence generates larger systems), and varies from 200 GiB at 100 kHz to 800 GiB at 300 kHz.
## 6. Conclusion

We have performed imaging of visco-acoustic media using an iterative reconstruction procedure. Attenuation can be encoded in the propagation via different models, and we have the following results.

1. We have implemented seven different models of attenuation and compared the wave propagation. It highlights that, while they can coincide at a reference frequency, each model leads to different wave patterns.

2. We have shown that our algorithm is robust with respect to attenuation model uncertainty, where we have used a relatively narrow band of frequency (here between 100 and 600 kHz). Namely, if one does not know the attenuation model that corresponds to the sample, it does not prevent from reconstructing its main properties.

3. While the wave sped and bulk modulus are accurately reconstructed and allow to describe the key-features of the sample, the recovery of the density and quality factor is more difficult, cf. [65,35] and the references therein. To recover them, one would need a parametrization that is more suitable to focus on amplitude variation, with an appropriate cost function (i.e., instead of the L2 norm), [35].

4. We have shown that the reconstruction assuming a propagation in free-space performs well but considering wall boundary around the samples leads to multiple reflections that drastically alter the reconstruction. As an alternative, we have proposed the use of complex frequencies which, by introducing an artificial damping, alleviate the difficulties and enable the accurate discoveries of the sample features. It provides us with a candidate to further consider media containing objects of high contrast that generate strong reflections, such as bones and skull, which are known to be harder to image. This is the subject of ongoing studies.

Our experiments use physical properties of tissues which, in this range of ultrasonic frequencies, are weakly attenuating (quality factor higher than 100). In the opposite situation where the medium is highly attenuating and prevent waves to propagate in the tissues with sufficient energy, the strategy we propose may have to be modified, this is part of ongoing works. In addition, the consideration of visco-elasticity is part of our ongoing research, where the difficulties come on the one hand from the incorporation of additional unknowns for inversion and on the other hand by the increased computational cost of solving a vector-wave problem.

### CRediT authorship contribution statement

**Florian Faucher:** Conceptualization, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. **Otmar Scherzer:** Conceptualization, Formal analysis, Supervision, Writing – review & editing.

### Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

Florian Faucher reports financial support was provided by Austrian Science Fund. Otmar Scherzer reports financial support was provided by Austrian Science Fund. Otmar Scherzer reports financial support was provided by Austrian Federal Government Ministry of Digital and Economic Affairs.

### Data availability

Data will be made available on request.

### Acknowledgements

The authors thank the reviewers and associate editor for helpful comments that have improved the quality of the manuscript. FF is funded by the Austrian Science Fund (FWF) under the Lise Meitner fellowship M 2791-N. For the experiments, we acknowledge the use of the cluster PlaFRIM,<fnref n="5" /> and of the Vienna Scientific Cluster VSC4.<fnref n="6" /> OS is supported by the Austrian Science Fund (FWF) with SFB F68 (Tomography Across the Scales), project F6807-N36 (Tomography with Uncertainties). The financial support by the Austrian Federal Ministry for Digital and Economic Affairs, the National Foundation for Research, Technology and Development and the Christian Doppler Research Association is gratefully acknowledged.

<footnote n="5">⁵ supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine, see https://www.plafrim.fr/</footnote>

<footnote n="6">⁶ https://vsc.ac.at/</footnote>
## Appendix A. Time-harmonic formulation with complex frequencies

To illustrate the derivation of the frequency-domain equation (2.3) for complex frequencies, we start with the time-domain wave equation. For simplicity, we consider the case without attenuation to avoid possibly integro-differential equation. In this case, time-dependent acoustic waves are given by the velocity and pressure fields, respectively $\hat{v}(x, t)$ and $\hat{p}(x, t)$, solutions to

$$ \rho(\mathbf{x}) \partial_t \hat{v}(\mathbf{x}, t) + \nabla \hat{p}(\mathbf{x}, t) = 0, \qquad (A.1a) $$

$$ \frac{1}{\kappa_0(\mathbf{x})} \partial_t \hat{p}(\mathbf{x}, t) + \nabla \cdot \hat{v}(\mathbf{x}, t) = \hat{g}(\mathbf{x}, t), \qquad (A.1b) $$

where $\hat{g}$ is the time-dependent source, that is, the Ricker wavelet in the context of Fig. 1.

We consider time-harmonic solutions of the following form:

$$ \hat{p}(\mathbf{x}, t) = p(\mathbf{x}, \omega) e^{-i(\omega_R + i\omega_l)t}, \qquad (A.2a) $$

$$ \hat{v}(\mathbf{x}, t) = v(\mathbf{x}, \omega) e^{-i(\omega_R + i\omega_l)t}, \qquad (A.2b) $$

$$ \hat{g}(\mathbf{x}, t) = g(\mathbf{x}, \omega) e^{-i(\omega_R + i\omega_l)t}, \qquad (A.2c) $$

where we remind that by definition, $\omega := \omega_R + i\omega_l$. Deriving with respect to time gives,

$$ \partial_t \hat{p}(\mathbf{x}, t) = -i\omega p(\mathbf{x}, \omega) e^{-i\omega t}, \qquad (A.3a) $$

$$ \partial_t \hat{v}(\mathbf{x}, t) = -i\omega v(\mathbf{x}, \omega) e^{-i\omega t}. \qquad (A.3b) $$

Replacing in (A.1), we obtain,

$$ -i\omega \rho(\mathbf{x}) v(\mathbf{x}, \omega) e^{-i\omega t} + \nabla(p(\mathbf{x}, \omega) e^{-i\omega t}) = 0, \qquad (A.4a) $$

$$ \frac{-i\omega}{\kappa_0(\mathbf{x})} p(\mathbf{x}, \omega) e^{-i\omega t} + \nabla \cdot (\mathbf{v}(\mathbf{x}, \omega) e^{-i\omega t}) = g(\mathbf{x}, \omega) e^{-i\omega t}. \qquad (A.4b) $$

Simplifying the terms $e^{-i\omega t}$ (that do not depend on space), we obtain the time-harmonic equation (2.3) in the case without attenuation. We note that this transformation introduces two parameters ($\omega_R$ and $\omega_l$), similarly to the continuous Gabor or the Wavelet transforms.

## References

[1] K. Aki, P.G. Richards, Quantitative Seismology, second ed., University Science Books, 2002.

[2] P.R. Amestoy, A. Buttari, J.-Y. L'excellent, T. Mary, Performance and scalability of the block low-rank multifrontal factorization on multicore architectures, ACM Trans. Math. Softw. 45 (2019) 1–26, https://doi.org/10.1145/3242094.

[3] E. Bachmann, J. Tromp, Source encoding for viscoacoustic ultrasound computed tomography, J. Acoust. Soc. Am. 147 (2020) 3221–3235, https://doi.org/10.1121/10.0001191.

[4] H. Barucq, G. Chavent, F. Faucher, A priori estimates of attraction basins for nonlinear least squares, with application to Helmholtz seismic inverse problem, Inverse Probl. 35 (2019) 115004, https://doi.org/10.1088/1361-6420.

[5] S. Beniwal, A. Ganguli, Defect detection around rebars in concrete using focused ultrasound and reverse time migration, Ultrasonics 62 (2015) 112–125, https://doi.org/10.1016/j.ultras.2015.05.008.

[6] S. Bernard, V. Monteiller, D. Komatitsch, P. Lasaygues, Ultrasonic computed tomography based on full-waveform inversion for bone quantitative imaging, Phys. Med. Biol. 62 (2017) 7011, https://doi.org/10.1088/1361-6560/aa7e5a.

[7] D.R. Bland, The Theory of Linear Viscoelasticity, Courier Dover Publications, 1960.

[8] R. Brossier, Two-dimensional frequency-domain visco-elastic full waveform inversion: parallel algorithms, optimization and performance, Comput. Geosci. 37 (2011) 444–455, https://doi.org/10.1016/j.cageo.2010.09.013.

[9] R. Brossier, S. Operto, J. Virieux, Which data residual norm for robust elastic frequency-domain full waveform inversion?, Geophysics 75 (2010) R37–R46, https://doi.org/10.1190/1.3379323.

[10] C. Bunks, F.M. Saleck, S. Zaleski, G. Chavent, Multiscale seismic waveform inversion, Geophysics 60 (1995) 1457–1473, https://doi.org/10.1190/1.1443880.

[11] J.M. Carcione, Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media, third ed., Elsevier, 2015.

[12] J.M. Carcione, D. Kosloff, R. Kosloff, Wave propagation simulation in a linear viscoelastic medium, Geophys. J. Int. 95 (1988) 597–611, https://doi.org/10.1111/j.1365-246X.1988.tb06706.x.

[13] B. Cockburn, J. Gopalakrishnan, R. Lazarov, Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems, SIAM J. Numer. Anal. 47 (2009) 1319–1365, https://doi.org/10.1137/070706616.

[14] K.S. Cole, R.H. Cole, Dispersion and absorption in dielectrics I. Alternating current characteristics, J. Chem. Phys. 9 (1941) 341–351.

[15] I. Donald, J. Macvicar, T. Brown, Investigation of abdominal masses by pulsed ultrasound, Lancet 271 (1958) 1188–1195, https://doi.org/10.1016/s0140-6736(58)91905-6.

[16] N. Duric, P. Littrup, A. Babkin, D. Chambers, S. Azevedo, E. Kalinin, R. Pevzner, M. Tokarev, E. Holsapple, O. Rama, et al., Development of ultrasound tomography for breast imaging: technical assessment, Med. Phys. 32 (2005) 1375–1386, https://doi.org/10.1118/1.1897463.
[17] N. Duric, P. Littrup, L. Poulo, A. Babkin, R. Pevzner, E. Holsapple, O. Rama, C. Glide, Detection of breast cancer with ultrasound tomography: first results with the computed ultrasound risk evaluation (CURE) prototype, Med. Phys. 34 (2007) 773-785, https://doi.org/10.1118/1.2432161.

[18] P. Elbau, O. Scherzer, C. Shi, Singular values of the attenuated photoacoustic imaging operator, J. Differ. Equ. 263 (2017) 5330-5376, https://doi.org/10.1016/j.jde.2017.06.018.

[19] F. Faucher, Contributions to seismic full waveform inversion for time harmonic wave equations: Stability estimates, convergence analysis, numerical experiments involving large scale optimization algorithms, PhD thesis, Université de Pau et Pays de l'Ardour, 2017, https://hal.archives-ouvertes.fr/tel-01807861.

[20] F. Faucher, *hawen*: time-harmonic wave modeling and inversion using hybridizable discontinuous Galerkin discretization, J. Open Sour. Softw. 6 (2021), https://doi.org/10.21105/joss.02699.

[21] F. Faucher, G. Alessandrini, H. Barucq, M. de Hoop, R. Gaburro, E. Sincich, Full reciprocity-gap waveform inversion, enabling sparse-source acquisition, Geophysics 85 (2020) R461-R476, https://doi.org/10.1190/geo2019-05271.

[22] F. Faucher, G. Chavent, H. Barucq, H. Calandra, A priori estimates of attraction basins for velocity model reconstruction by time-harmonic Full Waveform Inversion and Data-Space Reflectivity formulation, Geophysics 85 (2020) R223-R241, https://doi.org/10.1190/geo2019-02511.

[23] F. Faucher, M.V. De Hoop, O. Scherzer, Reciprocity-gap misfit functional for distributed acoustic sensing, combining data from passive and active sources, Geophysics 86 (2020) 1-46, https://doi.org/10.1190/geo2020-0305.1.

[24] F. Faucher, C. Kirisits, M. Quellmalz, O. Scherzer, E. Setterqvist, Diffraction tomography, Fourier reconstruction, and full waveform inversion, arXiv preprint, https://arxiv.org/abs/2110.07921, 2021.

[25] F. Faucher, O. Scherzer, Adjoint-state method for Hybridizable Discontinuous Galerkin discretization, application to the inverse acoustic wave problem, Comput. Methods Appl. Mech. Eng. 372 (2020) 113406, https://doi.org/10.1016/j.cma.2020.113406.

[26] F. Faucher, O. Scherzer, H. Barucq, Eigenvector models for solving the seismic inverse problem for the Helmholtz equation, Geophys. J. Int. 221 (2020) 394-414, https://doi.org/10.1093/gji/gga009.

[27] A. Fenster, S. Tong, H.N. Cardinal, C. Blake, D.B. Downey, Three-dimensional ultrasound imaging system for prostate cancer diagnosis and treatment, IEEE Trans. Instrum. Meas. 47 (1998) 1439-1447, https://doi.org/10.1109/19.746709.

[28] A. Fichtner, B.L. Kennett, H. Igel, H.-P. Bunge, Theoretical background for continental- and global-scale full-waveform inversion in the time-frequency domain, Geophys. J. Int. 175 (2008) 665-685, https://doi.org/10.1111/j.1365-246X.2008.03923.x.

[29] W. Futterman, Dispersive body waves, J. Geophys. Res. 67 (1962) 5279-5291, https://doi.org/10.1029/JZ067i013p05279.

[30] G.H.F. Gardner, L.W. Gardner, A.R. Gregory, Formation velocity and density – the diagnostic basics for stratigraphic traps, Geophysics 39 (1974) 770-780, https://doi.org/10.1190/1.1440465.

[31] J.M. Golden, G.A. Graham, Boundary Value Problems in Linear Viscoelasticity, first ed., Springer-Verlag, Berlin Heidelberg, 1988.

[32] L. Guasch, O.C. Agudo, M.-X. Tang, P. Nachev, M. Warner, Full-waveform inversion imaging of the human brain, npj Digit. Med. 3 (2020) 1-12, https://doi.org/10.1038/s41746-020-0240-8.

[33] W. Jeong, H.-Y. Lee, D.-J. Min, Full waveform inversion strategy for density in the frequency domain, Geophys. J. Int. 188 (2012) 1221-1241, https://doi.org/10.1111/j.1365-246X.2011.05314.x.

[34] R. Kamei, R.G. Pratt, Inversion strategies for visco-acoustic waveform inversion, Geophys. J. Int. 194 (2013) 859-884, https://doi.org/10.1093/gji/ggt109.

[35] H. Karaoğlu, B. Romanowicz, Global seismic attenuation imaging using full-waveform inversion: a comparative assessment of different choices of misfit functionals, Geophys. J. Int. 212 (2017) 807-826, https://doi.org/10.1093/gji/ggx442.

[36] R.M. Kirby, S.J. Sherwin, B. Cockburn, To CG or to HDG: a comparative study, J. Sci. Comput. 51 (2012) 183-212, https://doi.org/10.1007/s10915-011-9501-7.

[37] D. Köhn, D. De Nil, A. Kurzmann, A. Przebindowska, T. Bohlen, On the influence of model parametrization in elastic full waveform tomography, Geophys. J. Int. 191 (2012) 325-345, https://doi.org/10.1111/j.1365-246X.2012.05633.x.

[38] H. Kolsky, The propagation of stress pulses in viscoelastic solids, Philos. Mag. 1 (1956) 693-710, https://doi.org/10.1080/14786435608238144.

[39] R. Kowar, O. Scherzer, X. Bonnefond, Causality analysis of frequency-dependent wave attenuation, Math. Methods Appl. Sci. 34 (2011) 108-124, https://doi.org/10.1002/mma.1344.

[40] P. Lailly, The seismic inverse problem as a sequence of before stack migrations, in: J.B. Bednar (Ed.), Conference on Inverse Scattering: Theory and Application, Society for Industrial and Applied Mathematics, 1983, pp. 206-220.

[41] C. Li, N. Duric, P. Littrup, L. Huang, In vivo breast sound-speed imaging with ultrasound tomography, Ultrasound Med. Biol. 35 (2009) 1615-1628, https://doi.org/10.1016/j.ultrasmedbio.2009.05.011.

[42] C. Li, G.S. Sandhu, O. Roy, N. Duric, V. Allada, S. Schmidt, Toward a Practical Ultrasound Waveform Tomography Algorithm for Improving Breast Imaging, Medical Imaging 2014: Ultrasonic Imaging and Tomography, vol. 9040, International Society for Optics and Photonics, 2014, p. 90401P.

[43] Y. Lou, W. Zhou, T.P. Matthews, C.M. Appleton, M.A. Anastasio, Generation of anatomically realistic numerical phantoms for photoacoustic and ultrasonic breast imaging, J. Biomed. Opt. 22 (2017) 041015, https://doi.org/10.1117/1.JBO.22.4.041015.

[44] F. Lucka, M. Pérez-Liva, B.E. Treeby, B.T. Cox, High resolution 3D ultrasonic breast imaging by time-domain full waveform inversion, arXiv preprint, https://arxiv.org/abs/2012.00755, 2021.

[45] M. Malinowski, S. Operto, A. Ribodetti, High-resolution seismic attenuation imaging from wide-aperture onshore data by visco-acoustic frequency-domain full-waveform inversion, Geophys. J. Int. 186 (2011) 1179-1204, https://doi.org/10.1111/j.1365-246X.2011.05098.x.

[46] J. McLaughlin, A. Thomas, J.-R. Yoon, Basic theory for generalized linear solid viscoelastic models, Tomogr. Inverse Transp. Theory 559 (2011) 101-134, https://doi.org/10.1090/conm/559/11075.

[47] G. Müller, Rheological properties and velocity dispersion of a medium with power-law dependence of Q on frequency, J. Geophys. 54 (1983) 20-29, https://geophysicsjournal.com/article/82.

[48] A.I. Nachman, J.F. Smith III, R.C. Waag, An equation for acoustic propagation in inhomogeneous media with relaxation losses, J. Acoust. Soc. Am. 88 (1990) 1584-1595, https://doi.org/10.1121/1.400317.

[49] J. Nocedal, S.J. Wright, Numerical Optimization, second ed., Springer Series in Operations Research, 2006.

[50] N. Ozmen, R. Dapp, M. Zapf, H. Gemmeke, N.V. Ruiter, K.W. van Dongen, Comparing different ultrasound imaging methods for breast cancer detection, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 62 (2015) 637-646, https://doi.org/10.1109/TUFFC.2014.006707.

[51] M. Pérez-Liva, J. Herraiz, J. Udías, E. Miller, B. Cox, B. Treeby, Time domain reconstruction of sound speed and attenuation in ultrasound computed tomography using full wave inversion, J. Acoust. Soc. Am. 141 (2017) 1595-1604, https://doi.org/10.1121/1.4976688.

[52] R.-E. Plessix, A review of the adjoint-state method for computing the gradient of a functional with geophysical applications, Geophys. J. Int. 167 (2006) 495-503, https://doi.org/10.1111/j.1365-246X.2006.02978.x.

[53] R.G. Pratt, L. Huang, N. Duric, P. Littrup, Sound-Speed and Attenuation Imaging of Breast Tissue Using Waveform Tomography of Transmission Ultrasound Data, Medical Imaging 2007: Physics of Medical Imaging, vol. 6510, International Society for Optics and Photonics, 2007, p. 65104S.

[54] R.G. Pratt, C. Shin, G.J. Hick, Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion, Geophys. J. Int. 133 (1998) 341-362, https://doi.org/10.1046/j.1365-246X.1998.00498.x.

[55] R.G. Pratt, R.M. Shipp, Seismic waveform inversion in the frequency domain, part 2: fault delineation in sediments using crosshole data, Geophysics 64 (1999) 902-914, https://doi.org/10.1190/1.1444598.
[56] A. Ribodetti, J. Virieux, Asymptotic theory for imaging the attenuation factor Q, Geophysics 63 (1998) 1767-1778, https://doi.org/10.1190/1.1444471.

[57] T.J. Royston, Z. Dai, R. Chaunsali, Y. Liu, Y. Peng, R.L. Magin, Estimating material viscoelastic properties based on surface wave measurements: a comparison of techniques and modeling assumptions, J. Acoust. Soc. Am. 130 (2011) 4126-4138, https://doi.org/10.1121/1.3655883.

[58] C. Shin, Y.H. Cha, Waveform inversion in the Laplace domain, Geophys. J. Int. 173 (2008) 922-931, https://doi.org/10.1111/j.1365-246X.2008.03768.x.

[59] C. Shin, Y.H. Cha, Waveform inversion in the Laplace Fourier domain, Geophys. J. Int. 177 (2009) 1067-1079, https://doi.org/10.1111/j.1365-246X.2009.04102.x.

[60] T.L. Szabo, Time domain wave equations for lossy media obeying a frequency power law, J. Acoust. Soc. Am. 96 (1994) 491-500, https://doi.org/10.1121/1.410434.

[61] A. Tarantola, Inversion of seismic reflection data in the acoustic approximation, Geophysics 49 (1984) 1259-1266, https://doi.org/10.1190/1.1441754.

[62] B. Ursin, T. Toverud, Comparison of seismic dispersion and attenuation models, Stud. Geophys. Geod. 46 (2002) 293-320, https://doi.org/10.1023/A:1019810305074.

[63] T. Van Leeuwen, W. Mulder, A correlation-based misfit criterion for wave-equation traveltime tomography, Geophys. J. Int. 182 (2010) 1383-1394, https://doi.org/10.1111/j.1365-246X.2010.04681.x.

[64] V. Vavryčuk, Weak anisotropy-attenuation parameters, Geophysics 74 (2009) WB203-WB213, https://doi.org/10.1190/1.3173154.

[65] J. Virieux, S. Operto, An overview of full-waveform inversion in exploration geophysics, Geophysics 74 (2009) WCC1-WCC26, https://doi.org/10.1190/1.3238367.

[66] S. Yakovlev, D. Moxey, R.M. Kirby, S.J. Sherwin, To CG or to HDG: a comparative study in 3D, J. Sci. Comput. 67 (2016) 192-220, https://doi.org/10.1007/s10915-015-0076-6.

[67] P. Yang, R. Brossier, L. Métivier, J. Virieux, A review on the systematic formulation of 3-D multiparameter full waveform inversion in viscoelastic medium, Geophys. J. Int. 207 (2016) 129-149, https://doi.org/10.1093/gji/ggw262.

[68] Ö. Yilmaz, Seismic Data Analysis, vol. 1, Society of exploration geophysicists, Tulsa, OK, 2001, https://doi.org/10.1190/1.9781560801580.

[69] C.M. Zener, Elasticity and Anelasticity of Metals, University of Chicago Press, 1948.

[70] B. Zhou, X. Zhang, Comparison of five viscoelastic models for estimating viscoelastic parameters using ultrasound shear wave elastography, J. Mech. Behav. Biomed. Mater. 85 (2018) 109-116, https://doi.org/10.1016/j.jmbbm.2018.05.041.