# Probing the speckle to estimate the effective speed of sound, a first step towards quantitative ultrasound imaging

Josselin Garnier, Laure Giovangigli, Quentin Goepfert, Pierre Millien

► To cite this version:

Josselin Garnier, Laure Giovangigli, Quentin Goepfert, Pierre Millien. Probing the speckle to estimate the effective speed of sound, a first step towards quantitative ultrasound imaging. 2025. hal-05064390

HAL Id: hal-05064390

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

Preprint submitted on 12 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PROBING THE SPECKLE TO ESTIMATE THE EFFECTIVE SPEED OF SOUND, A FIRST STEP
TOWARDS QUANTITATIVE ULTRASOUND IMAGING.

JOSSELIN GARNIER<fnref n="1" />, LAURE GIOVANGIGLI<fnref n="2" />, QUENTIN GOEPFERT<fnref n="2" />, AND PIERRE MILLIEN<fnref n="3" />

**ABSTRACT.** In this paper, we present a mathematical model and analysis for a new experimental method [Bureau and al., arXiv:2409.13901, 2024] for effective sound velocity estimation in medical ultrasound imaging. We perform a detailed analysis of the point spread function of a medical ultrasound imaging system when there is a mismatch between the effective sound speed in the medium and the one used in the backpropagation imaging functional. Based on this analysis, an estimator for the speed of sound error is introduced. Using recent results on stochastic homogenization of the Helmholtz equation, we provide a representation formula for the field scattered by a random multi-scale medium (whose acoustic behavior is similar to a biological tissue) in the time-harmonic regime. We then prove that statistical moments of the imaging function can be accessed from data collected with only one realization of the medium. We show that it is possible to locally extract the point spread function from an image constituted only of speckle and build an estimator for the effective sound velocity in the micro-structured medium. Some numerical illustrations are presented at the end of the paper.

**Mathematics Subject Classification (MSC2020).** 35R30, 92C55,35R60,35J05

**Keywords.** Helmholtz equation, quantitative ultrasound imaging, stochastic homogenization.

# 1. INTRODUCTION

**Scientific context.** Ultrasound imaging is a cheap, safe portable way of imaging soft tissues in real time. Medical ultrasound exams account for about 20% of the imaging exams [10], and can be used for numerous different applications, from obstetric to cardiology or dermatology [30]. The principle is the following: an ultrasonic wave is emitted by a probe to insonify a specific organ (womb, liver, brain...). The inhomogeneities in the tissues produce echoes (*backscattered wave*) that are measured on the probe. An image showing the localization of the sources of the echoes is computed using a *backpropagation* type algorithm (Sum and delay, Kirchoff migration...). These types of *time reversal* based algorithms are all relying on the principle of a *space-time* correspondence (travel time) and require some *a priori* hypothesis on the speed of sound in the medium. The speed of sound contrast between different types of soft tissues rarely exceeds 10% and in practice most commercial devices use a constant value for the speed of sound, usually close to the one observed in water. Nevertheless this approximation can deteriorate the quality of the images significantly in certain situations by introducing artefacts, distorting geometric features or reducing the contrast. Beyond the improvement of the images, the knowledge of the speed of sound in certain organs can be used as an indicator of a pathological state (such as hepatic steatosis [8]). We refer to the review [9] for a complete overview of the current state of the art on quantitative ultrasound imaging technics. The stakes of measuring *in situ* the propagation speed of ultrasounds in soft tissues is therefore twofold: improve the reliability of the ultrasound images and provide new sensitive diagnostic tools.

**The reflexion matrix approach.** An important evolution in medical ultrasound imaging was the development of new insonification processes. By using a spherical wave or a plane wave, it is possible to illuminate the whole medium with one wave and to get confocal images in post-processing : instead of focusing physically in the medium like in conventional ultrasound imaging modalities, the focusing is done numerically by exploiting the linearity of the wave equation and the superposition principle.This method has allowed for a

<footnote n="1">¹CMAP, CNRS, ECOLE POLYTECHNIQUE, INSTITUT POLYTECHNIQUE DE PARIS, 91120 PALAISEAU, FRANCE.</footnote>

<footnote n="2">²POEMS, CNRS, INRIA, ENSTA PARIS, INSTITUT POLYTECHNIQUE DE PARIS, 91120 PALAISEAU, FRANCE</footnote>

<footnote n="3">³INSTITUT LANGEVIN, ESPCI PARIS, PSL UNIVERSITY, CNRS, 1 RUE JUSSIEU, 75005 PARIS, FRANCE.</footnote>

<footnote>E-mail addresses: pierre.millien@espci.fr.</footnote>
drastic improvement in the frame rate (from 24fps up to 1000 fps) with no loss in contrast or resolution [32]. Another advantage is that since the image processing is done numerically, any extra information about the speed of sound can be incorporated *a posteriori* to improve the images. Speed of sound estimation is an active area of research and several groups have proposed innovative methods based on adjusting the speed of sound by maximizing a *quality indicator*, such as echo amplitude [34], a coherence factor [18, 17], or the brightness of a strong reflector [28]. The *computed ultrasound tomography in echo mode* (CUTE) method [19, 20, 5] uses a localized detection of aberration phase shifts under different angles of illumination and inverts a linear model that relates these measurements to the speed of sound in tissues.

Over the last few years, the team of A. Aubry and M. Fink have obtained spectacular experimental results in the domain of aberration correction and speed of sound estimation using a *reflection matrix* approach [3, 24, 23] (see section 2.2). In [7] a new method to estimate the speed of sound in soft tissues has been introduced. The idea, inspired from adaptive optics, is also to assess the speed of sound in the medium by maximizing a *quality indicator* (the amplitude at the center of the focal spot). The main strength of this method is that it does not require the presence of a strong isolated reflector (guide star). The idea is to probe the point spread function directly in the speckle part of the image by using a clever focusing sequence that mimics the presence of an isolated reflector. The goal of the present work is to present a robust mathematical model and analysis to support their experimental proof of concept of the method and quantify precisely some of the effects observed experimentally.

**Outline of the article and main contributions.** In this paper, we start by presenting a general model for medical ultrasound imaging set ups. Particular attention is given to incorporating the specific features of these experiments, whether it is in the description of the acoustic properties of soft tissues with a stochastic model or the use of an asymptotic regime that describes accurately the different scales inherent to the physics of the problem. Building on our recent results on the high order stochastic homogenization of Helmholtz equation [13] we give an approximation of the measurements in the time-harmonic regime for a random multi-scale medium. In section 3, we perform an asymptotic analysis of the point spread function of a medical ultrasound imaging system in the case where the speed of sound used in the backpropagation imaging functional is not the effective speed of sound in the medium. We define a general notion of *focal spot* and quantify the effects of the misfit in speed of sound on the focal spot (displacement, deformation...). Based on this analysis, an estimator for the speed of sound error is introduced (proposition 3.1). In section 4, we perform a mathematical analysis of the method developed by A. Aubry and his collaborators [7]. We show that in a micro-structured medium, it is possible to access the average of the point spread function over the focal spot via the variance of the value of the imaging function at well chosen pixels (lemma 4.1). We exhibit an estimator for the effective speed of sound (proposition 4.2). In section 4.4, we show how to approximate statistical moments of the imaging function with data collected from a single realization of the medium using stationarity properties of the medium. Proposition 4.4 gives a quantification of the approximation error for the variance of the imaging function. The last section includes some numerical simulations of the problem to illustrate the claims of the paper.

**Position with respect to the state of the art.** In the mathematics community, the problem of velocity estimation from boundary measurements for various types of waves has been an active area of research for decades. The usefulness of any method depends strongly on the practical experimental constraints. For example, methods based on the inversion of the Eikonal equation [31] are of little use in situations where the sensors cannot surround the medium to image. Similarly, powerful *full wave inversion* methods [11] that can perform remarkably to reconstruct the velocity of seismic waves inside earth's crust will fall short if the imaging process is required to be in real time, due to their algorithmic complexity. The rule of thumb is that the more *a priori* information a method can incorporate, the better it will perform in specific situations. For example, in the case of medical ultrasound imaging, the data that are measured (back-scattered field) are generated by the biological tissue's micro-structure, while the properties of interest (effective speed of sound and travel times) appear at a macroscopic scale. To the best of our knowledge, there is no recent mathematical theory studying the effective velocity estimation problem from boundary measurements in a reflection geometry for a random micro-structured medium, particularly using reconstruction methods with a low algorithmic complexity.
## 2. MODEL FOR MEDICAL ULTRASOUND EXPERIMENTS

**2.1. Geometry and governing equations.** The general model for the governing equation is a divergence form wave equation in an unbounded domain.

**The coefficients.** We consider a $d$-dimensional framework. Take a simply connected open bounded domain $D$ with at least Lipschitz boundary regularity, modeling a medium embedded in an infinite homogeneous isotropic medium. Introduce the two coefficients:

$$ \begin{cases} a(x) := a_m \mathbb{1}_{\mathbb{R}^d \setminus \bar{D}}(x) + a_D(x) \mathbb{1}_D(x), \\ n(x) := n_m \mathbb{1}_{\mathbb{R}^d \setminus \bar{D}}(x) + n_D(x) \mathbb{1}_D(x), \end{cases} \quad x \in \mathbb{R}^d, \quad (1) $$

where $a_D, n_D \in L^\infty(D)$. Moreover, we assume that the functions $a$ and $n$ are uniformly positive. In acoustics (which is the context of medical ultrasound imaging), the coefficients $a$ and $n$ are related to the inverse of the density (usually denoted $\rho$) and the inverse of the bulk modulus (usually denoted $\kappa$). See for instance [27, section 3.3.3]. For simplicity, in this paper we make the restrictive assumption that

$$ a_D \equiv a_m = 1. $$

**Remark 2.1.** This assumption seems restrictive at first as the model would only be accurate for a medium with a constant density. Moreover, mathematically, the representation formulas for the scattered field are very different when the contrast is in the divergence part of the equation. However, this is a common hypothesis in the medical ultrasound community, where most models do not have contrast in the divergence. This is motivated by the fact that when the scatterers are dilute (at distances larger than the wavelength) and do not exhibit a particular orientation or elongation, it is impossible to distinguish a dipole from a monopole in the reflection geometry. Note that outside the dilute case the question under which conditions are the two contrast models equivalent ? remains open.

**Time domain versus frequency domain.** All the experiments are done in the time domain, and backscattered echoes are measured as time-dependent functions. Nevertheless, besides the simple sum and delay algorithms [26], most of the advanced post-processing is usually done in the frequency domain. The introduction of the notion of wavelength is critical to quantify the resolution of the methods so our mathematical analysis will be done in the time harmonic domain. A Laplace-Fourier transform can be used to switch from one domain to another. In the rest of the chapter the frequency will be denoted by $\omega$. Since the transducers have a finite bandwith $BW := [\omega_-, \omega_+]$ we will always assume $\Re\omega \in BW$. We will only write the dependency of the wavefields with respect to $\omega$ when it is needed. Usually, the time profile of the incident wave at a given point has a central frequency $\omega_c$:

$$ BW = [\omega_c - B, \omega_c + B], $$

where $B$ is the bandwidth. If the pulse consists of a few oscillations at the central frequency (see an example fig. 5) then we are in a broadband situation $B = O(\omega_c)$. If the pulse consists of oscillations at the central frequency with an enveloppe decaying slowly compared to the central period then we are in a narrowband case $B \ll \omega_c$.

**The incident wave.** The choice of the incident wave $u^i$ has little impact on the mathemati-cal model for the direct problem but is of crucial importance in the applications and for the inverse problem. The simplest model for an incident wave is to consider an eigenfunction for the free space laplacian, i.e. a plane wave of amplitude $U^i$ propagating in direction $\theta \in S^{d-1}$ with wavenumber $k_m = \sqrt{n_m}\omega$:

$$ u^i(x) = U^i e^{ik_m \theta \cdot x}, \quad x \in \mathbb{R}^d. $$

This is a good model when the source is far enough from the domain $D$ to be considered at infinity.

It is possible to use more sophisticated models for the incident wave. For example a point source outside of the domain $\bar{D}$ or a finite sum of point sources. In this case the incident wave is expressed as a linear combination of Green functions:

$$ u^i(x) = \sum_{j=1}^{N} \alpha_j \Gamma^{k_m}(x, x_j), \quad x \in \mathbb{R}^d \setminus \cup_j \{x_j\}, $$
with $\alpha_j \in \mathbb{C}$. This is a good model for the medical ultrasound experiments, where probes, which can be close to the medium, are made of a finite number of transducers. When the incident wave is emitted via a point source located at emission point $x_e$ we will denote it $u^i(x_e, \cdot)$ to keep track of the dependency on the location of the source when it is needed, especially for the inverse problem. Sometimes though a continuous model for the probe can be more suited for a mathematical analysis. In this case, one can model the probe by a compact manifold $\mathcal{P}$ outside of $\bar{D}$ and the incident wave can be modeled by a single or a double layer potential on that manifold:

$$u^i(x) = S_{\mathcal{P}}^{k_m}(\varphi)(x) \quad \text{or} \quad u^i(x) = D_{\mathcal{P}}^{k_m}(\varphi)(x) \quad x \in \mathbb{R}^d \setminus \mathcal{P},$$

where the layer potentials are given by [2]:

$$S_{\mathcal{P}}^{k_m}(\varphi)(x) = \int_{\mathcal{P}} \varphi(y) \Gamma^{k_m}(x, y) \mathrm{d}\sigma(y) \quad \text{and} \quad D_{\mathcal{P}}^{k_m}(\varphi)(x) = \int_{\mathcal{P}} \varphi(y) \frac{\partial \Gamma^{k_m}}{\partial v(x)}(x, y) \mathrm{d}\sigma(y).$$

The profile of the incident wave in the time domain $f(t)$ can be easily added to this model by adding its Fourier-Laplace transform $\hat{f}(\omega)$ as a prefactor. The frequency content does not affect much the direct problem but has a major impact on the reconstructed images' resolution as will be detailed in section 3.3

**The governing equation.** The scattering problem by the medium $D$ is modeled by the following Helmholtz equation in $\mathbb{R}^d$:

$$\left\{
\begin{array}{ll}
\Delta u + \omega^2 n u = 0 & \text{in } \mathbb{R}^d \\
\lim_{|x| \to \infty} |x|^{\frac{d-1}{2}} \left( \partial_{|x|} (u - u^i) - ik_m (u - u^i) \right) = 0. &
\end{array}
\right.
\quad (2)$$

This equation has a unique solution in $H_{\text{loc}}^2(\mathbb{R}^d)$. The well-posedness is well known and can be found in textbooks [21]. The difficulty of equation (2) lies in getting either sharp regularity results for the solution or some explicit control with respect to the frequency. Usually some extra assumptions on the coefficients are needed. These questions have been a very active area of research over the last decades. We refer to the excellent paper [25] and references within for a recent overview of these questions.

## 2.2. Measurement model: the reflection matrix.
As mentioned previously, in practice the probe has a finite number $N$ of transducers that can be modeled as a finite collection of point sources $(x_i)_{i=1...N}$ each on the manifold $\mathcal{P}$ describing the shape of the probe. For medical imaging modelization it is reasonable to assume the probe is linear and positioned just above the medium. We set

$$\mathcal{P} = [-\ell, \ell]^{d-1} \times \{0\}.$$

Each transducer can be controlled independently. The measurements are obtained in the following way:

* one transducer at $x_e$ is activated to emit an incident wave $u^i(x_e, \cdot, \omega)$,

* the scattered wave $u^s(x_e, \cdot, \omega) = u(x_e, \cdot, \omega) - u^i(x_e, \cdot, \omega)$ is recorded by all the transducers $x_r$ on the probe for all frequencies in the bandwidth,

* the operation is repeated with another transducer $x_e'$ and so on.

In the end one has access to the measurement matrix at each frequency

$$M_{e,r}(\omega) = u^s(x_e, x_r, \omega) \in \mathbb{C}^{N \times N}.$$

The continuous model to describe this type of measurements is to assume we are given the measurement map:

$$M \in L^2(\mathcal{P} \times \mathcal{P} \times \mathcal{B}W)$$

given by

$$M(x_e, x_r, \omega) = u^s(x_e, x_r, \omega). \quad (3)$$

**Remark 2.2.** In order to improve the signal-to-noise ratio in the presence of additive noise, multiplexing measurement strategies can be used [15], in which the transducers transmit orthogonal waveforms and the measurement matrix can then be extracted numerically.
**2.3. Model for the biological tissues.** We now want to give a more precise model for the coefficient $n$ inside the medium $D$. The signal measured by the probe comes from the local discontinuities in the acoustic impedance. To correctly capture the nature of the experiment, the model for the coefficient has to take into account both types of local discontinuities: the ones at the interface between different types of soft tissues (fat, liver, muscle etc...) and the ones occurring inside each tissue type due to the presence of small unresolved scatterers. The model was introduced in [13], we present a simple version that we will use in this work.

**Simple first model.** In this context, the propagating medium could be described as a composite medium, with a few components, each describing a different tissue type. Then the presence of small unresolved scatterers could be modeled by randomly sprinkling each component with a large but finite number of small or point-like scatterers. This looks like a perfectly sound model and has been successfully used to describe blood flow imaging experiments [1]. If there are enough small scatterers the numerical simulation of the propagation in such a medium will look like an ultrasound imaging experiment. Most estimates and tools available for the study of composite mediums are valid for a fixed number of scatterers and we do not know how the constants in the different estimates will behave when the number $N$ of scatterers grows very large. Moreover, the state of the art models do not account for the fact that a change in the properties of the scatterers (their density in the medium, their contrast...) induces a change in the acoustic properties of the medium at the macroscopic scale. A good way to convince oneself of this is to look at how the experimental phantoms are made. They usually consist of an incompressible gel with acoustic properties similar to those of water, with microscopic acoustic grains embedded inside the gel everywhere. Variations at the macro scale in acoustic properties are obtained by changing locally the repartition of the unresolved acoustic grains inside. The other strong hint towards the fact that macro-scale acoustic properties of tissues can be derived by an homogeneization process is the experimental observation of anisotropy. Since the compressibility and density are both scalar parameters, an homogeneization process seems like a necessary way to understand the acoustic anisotropy that can sometimes be observed in certain biological tissues like muscles.

**Advanced model.** In light of the previous paragraph, we aim at providing a model that can link the properties of the scatterers to the effective acoustic parameters of the medium. The good way to do that is to use a model that preserves the proportion of volume occupied by the scatterers even when the size of the scatterers decreases.

We start by considering a random distribution of scatterers of characteristic size $r_0 > 0$ in $\mathbb{R}^d$ and of average density of order one. We then rescale this distribution using a scaling parameter $\epsilon > 0$ to place the scatterers inside the medium $D$.

Let $(x_i)_{i \in \mathbb{N}}$ be the point process in $\mathbb{R}^d$ corresponding to the centers of the scatterers. Consider an open connected domain $Q$ of measure $r_0^d$ centered at 0. Each scatterer is then represented by $S_i = \{x_i + Q\}$. We make the following assumptions on $(x_i)_{i \in \mathbb{N}}$:

- $(x_i)_{i \in \mathbb{N}}$ is stationary, i.e. its probability distribution is invariant by translation and ergodic;

- the scatterers lie at a distance at least $\delta_{min} > 0$ from one another, i.e. there exists $\delta_{min} > 0$ such that

$$ \forall i \neq j, \text{dist}(s_i, s_j) > \delta_{min} \quad \text{a.s.} $$

We denote by $S = \cup_i \{x_i + Q\}$. We introduce the parameter $n : \mathbb{R}^d \times \Omega \to [n_-, n_+]$

$$ n := n_0 \mathbb{1}_{\mathbb{R}^d \setminus \bar{S}} + \sum_i n_{S_i} \mathbb{1}_{S_i}, \quad (4) $$

where $(n_{S_i})_i$ are independent and identically distributed. We then fix $\epsilon > 0$ such that $\epsilon r_0$ is now a length representing the characteristic size of the scatterers in the medium. The natural definition for the parameter $n_D$ in the domain $D$ in the case where scatterers are of size of order $\epsilon$ is:

$$ n_D(x) = n_{\epsilon}(x) := n\left(\frac{x}{\epsilon}\right). $$

**Remark 2.3.** Even though the characteristic size of the scatterers is going to zero as the scaling parameter $\epsilon$ is going to zero, the volume fraction they occupy in $D$ stays constant as the number of scatterers inside $D$ grows as $|D|r_0^{-d}\epsilon^{-d}$. This is the main difference with the other small scatterer models like the Foldy-Lax model.

**Remark 2.4.** There is no need to take the same shape for all the scatterers. The shape $Q$ can be replaced by $Q_i$ where the shapes $Q_i$ are realizations of an i.i.d. process.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Illustration of the geometry for an ultrasound imaging experiment.</figcaption></figure>

**Remark 2.5.** Our results apply to any random medium described by an index of the form $f_{\epsilon}(x) = f(\frac{x}{\epsilon})$ where $f: \mathbb{R}^d \times \Omega \to [n_{-}, n_{+}]$ is stationary and verifies (6).

We assume moreover that the distribution of scatterers possess some decorrelation properties. Namely, let us denote by $C$ the covariance function of $n$:

$$C(x) := \mathrm{Cov}(n(x), n(0)), \quad \forall x \in \mathbb{R}^d. \qquad (5)$$

We assume that $C$ can be controlled by a bounded sufficiently decaying decreasing function $\Phi: \mathbb{R}^{+} \to \mathbb{R}^{+}$ i.e such that for some $p > 2(d-1)$ (see remark A.1):

$$\int_{\mathbb{R}^d} \Phi(|x|)|x|^p dx < \infty. \qquad (6)$$

In particular, $\Phi$ and $\Phi^{\frac{1}{2}}$ are in $L^1(\mathbb{R}^+)$.

**2.4. Homogenized representation formula.** The corner stone of our analysis is the simplified asymptotic expression for the measurements obtained with quantitative stochastic homogenization techniques in our previous paper [13]. More precisely, we use the integral representation formula of [13, theorem 10]. We introduce $G^*$ the Green function associated to the homogenized problem i.e. $G^*$ is the outgoing solution in $D' (\mathbb{R}^d)$ of all

$$-\Delta G^*(\cdot, y) - \tilde{n}^* \omega^2 G^*(\cdot, y) = \delta(\cdot - y) \quad \text{in } \mathbb{R}^d, \qquad (7)$$

with $\tilde{n}^*$ being the homogenized coefficient in the object $D$ and the coefficient of the surrounding medium outside of $D$:

$$\tilde{n}^*(x) := n_m \mathbb{1}_{\mathbb{R}^d \setminus \bar{D}}(x) + n^* \mathbb{1}_D(x), \quad x \in \mathbb{R}^d. \qquad (8)$$
For the simplicity of the computations we can make the following assumption that $n_0 = n_m$ and that $\mathbb{E}[n_S] = 0$ so that $n^* = n_0$, $\tilde{n}^* = n_0$, and $G^* = \Gamma_{c^*}^{\omega}$ with $c^* = \frac{1}{\sqrt{n^*}}$. In that case the coherent wave is the incident field $\Gamma_{c^*}^{\omega}(\cdot, x_e)$. This can be viewed as an effective matched boundary condition.

**Remark 2.6.** This condition is restrictive in theory but mildly restrictive in practice. In case of a mismatch of the effective speed of sound, some reflections will occur at the boundary of the domain. The first one from the top of the domain D can be filtered in the time domain. The ones occurring at the side of the domain will not pollute too much the measurements in the geometry considered, due to the directivity of the sensors.

We recall that in that case, the scattered field can be approximated by [13, theorem 10]

$$ \mathcal{U}^s(x_e, y, \omega) := \omega^2 \int_D (n_{\varepsilon}(x) - n^*) \Gamma^{k^*}(x, x_e) \Gamma^{k^*}(y, x) dx, \quad y \in \mathbb{R}^d \setminus \bar{D}, \qquad (9) $$

and the pointwise error control is given in proposition B.1.

## 2.5. Asymptotic regime.

**Typical values for the parameters.** In practical situations, the typical central frequency of the transducers is of the order of a few Mhz, the speed of sound varies from $1400 - 1600 \text{ ms}^{-1}$, and thus the wavelength is of order $\lambda \sim 1 \text{mm}$ while the typical diameter of D is of order $10 \text{cm}$. Finally the size of the transducer array is usually a few centimeters $2\ell \sim 4 \text{cm}$ (see [6, Table 2.2]). The scatterers are usually of typical size $\varepsilon r_0$ of a few micron (think for example of a red blood cell), around two order of magnitude below the central wavelength. Mathematically this can be expressed by

$$ \varepsilon r_0 \ll \frac{c^*}{\omega_c} \ll \ell \ll \text{diam}(D). $$

**The paraxial regime.** In order to get analytic expressions, we will consider a classical asymptotic regime called the *paraxial regime* [14, chapter 6]. This will be done by introducing a small scaling parameter denoted $\eta \ll 1$ and expressing the geometric quantities of the model as functions of $\eta$. We introduce

$$ \omega_c := \frac{\omega_0}{\eta}, \quad \mathcal{B} := \frac{\mathcal{B}_0}{\eta}, \quad \mathcal{P} := \eta^{\frac{1}{2}} \mathcal{P}_0, \quad \ell := \eta^{\frac{1}{2}} \ell_0, \quad \varepsilon := \eta^{\alpha}, \quad \alpha > 0. $$

Regarding the points of the domain, we will say $y \in D$ satisfies the paraxial approximation if its transverse and axial coordinates can be written as:

$$ y := \left( \eta^{\frac{1}{2}} y^{\perp}, y^{\parallel} \right) \quad y^{\perp} \in \mathbb{R}^{d-1}, y^{\parallel} \in \mathbb{R}^{+}. $$

For the reminder of the manuscript the notation $y^{\perp}$ will designate the projection on the transverse component rescaled with respect to the asymptotic parameter $\eta$:

$$ \mathbb{R}^d \rightarrow \mathbb{R}^{d-1} \\ y \mapsto y^{\perp} := \eta^{-\frac{1}{2}} \mathbf{P}^{\perp}(y) $$

where $\mathbf{P}^{\perp}$ is the projection on the transverse plane.

**Remark 2.7.** For the approximation of eq. (9) to be accurate the error term of equation (31) in appendix A needs to be small, i.e. goes to 0 when $\eta$ goes to zero. This imposes a lower bound on $\alpha$. The bound depends on the dimension $d$ and the exponent $p'$ introduced in the appendix. In dimension $d = 3$, $\alpha > \frac{5}{2}$ is a sufficient condition.

## 3. THE IMAGING FUNCTIONAL

### 3.1. The classical confocal imaging function.
In the case where the effective speed of sound $c^*$ in the medium D is known the classical confocal imaging function can be written:

$$ I^{c^*}(z) := \int_{(\mathcal{P} \times \mathcal{P} \times B)} \overline{M(x_e, x_r, \omega)} G^*(z, x_e) G^*(z, x_r) dx_e dx_r d\omega. \quad z \in D \qquad (10) $$

This function, also known as the *Kirchhoff migration function* has been extensively studied in various cases (see [12] and references within).
In this context it is straightforward from the representation formula (9) that the measurement matrix M can be written as

$$M(x_e, x_r, \omega) = \omega^2 \int_D (n(y) - n^*) G^*(x_r, y) G^*(x_e, y) dy. \quad (11)$$

**Remark 3.1.** Even though $I^c$ is defined for $z \in D$ it is important to keep in mind that the imaging domain is a virtual domain while the integral in equation (11) happens physically in the real object D. This will be of importance in the rest of the work and justifies the following definition to avoid the confusion between the imaging points and the physical points in the object.

**Definition 3.1 (Virtual Domain).** The definition domain for $I^c$ will be denoted $D'$.

The imaging function can then be expressed as

$$I^c(z) = \int_D (n_e(y) - n^*) F^c(z, y) dy, \quad z \in D'$$

which is an integral operator acting on $n_e - n^*$ with kernel:

$$F^c(z, y) = \int_B \omega^2 \left( \int_P \Gamma^{\frac{\omega}{c^*}}(z, x_r) \overline{G^*(y, x_r)} dx_r \right)^2 d\omega, \quad (z, y) \in D' \times D.$$

which acts as the so-called *point spread function* for the imaging system.

**3.2. Imaging with an unknown speed of sound.** In our case the effective speed of sound in the medium D is unknown and we can parametrize the imaging function by the speed of sound c used in the backpropagation.

**Definition 3.2.**

$$I^c(z) := \int_{(P \times P \times B)} \overline{M(x_e, x_r, \omega)} \Gamma^{\frac{\omega}{c}}(z, x_e) \Gamma^{\frac{\omega}{c}}(z, x_r) dx_e dx_r d\omega, \quad z \in D'. \quad (12)$$

**Remark 3.2.** The measurements are done in the time domain on a finite interval $[0, T]$. The imaging domain $D'$'s axial size is determined by the longest distance a wave can travel in $\frac{T}{2}$. Note that now that the speed of sound in the virtual domain is modified the domain $D'$ is now different from D. The explicit expression for the mapping $\varphi_c$ between D and $D'$ is given in lemma 3.1 (eq. (15)).

We can still express the imaging functional as an integral operator

$$I^c(z) = \int_D (n_e(y) - n^*) F^c(z, y) dy$$

acting on $n_e - n^*$ with kernel:

$$F^c(z, y) = \int_B \omega^2 \left( \int_P \Gamma^{\frac{\omega}{c}}(z, x_r) \overline{G^*(y, x_r)} dx_r \right)^2 d\omega, \quad z, y \in D' \times D.$$

which is the new *point spread function* for the imaging system.

**3.3. Properties of the point spread function.**

**3.3.1. Focal spot in the case $c = c^*$:** The point spread function at the correct backpropagation speed $F^c(z, y)$ has already been extensively studied in that setting (see for instance [12]). For a fixed $y_0 \in D$, the function

$$D' \to C \\ z \mapsto F^c(z, y_0)$$

presents a peak centered at $y_0$: its amplitude is of order $\frac{1}{\eta}$ in a small region around $y_0$ called the *focal spot* and negligible outside of this region. The size of the focal spot depends on the geometry of the problem: the central wavelength, the bandwidth, the numerical aperture... The point spread function is of course symmetrical in that case since $D = D'$ and $F^c(z, y) = F^c(y, z)$. See fig. 4, middle insert, for an example corresponding to the geometry of medical ultrasound imaging methods.
We can introduce a more precise definition of the focal spot. For this we need to introduce the family of level sets in the image domain:

$$L_{\Theta}(y_0) := \{z \in D', |F^{c^*}(z, y_0)| = \Theta\}.$$

Now we consider the family of nested and connected sets $D'(y_0, \Theta)$ defined as the minimal volume with $C^{0,\alpha}$ boundary (with $\alpha > 0$) containing the level set $L_{\Theta}(y_0)$, i.e.

$$D'(y_0, \Theta) := \underset{L_{\Theta} \subset \Omega' \subset D', \Omega' \ni y_0}{\operatorname{argmin}} |\Omega'|.$$

The focal spot can now be defined:

**Definition 3.3 (Focal spot at threshold $\delta$).** Consider $0 < \delta \ll 1$ and $y_0 \in D$. There exists a $\Theta(\delta) > 0$ such that

$$\int_{D' \setminus \overline{D'(y_0, \Theta(\delta))}} |F^{c^*}(z, y_0)| dz < \delta \int_{D'} |F^{c^*}(z, y_0)| dz.$$

We denote $D'_\delta(y_0) := D'(y_0, \Theta(\delta))$.

**Remark 3.3.** The existence of $D'_\delta(y_0)$ is obvious as $F \in L^1(D' \times D)$. Uniqueness comes from the uniqueness of $D'(y_0, \Theta)$ for a given $\Theta$ and the monotonicity of the function

$$\Theta \mapsto \int_{D' \setminus \overline{D'(y_0, \Theta)}} |F^{c^*}(z, y_0)| dz.$$

Such a definition only has a practical interest because $F(\cdot, y_0)$ is a peaked function and therefore in practice the domain $D'_\delta$ is small compared to the size of the domain, see remark 3.5. Therefore $\int_{D'_\delta(y_0)} |F^{c^*}(z, z^*)| dz$ is a good approximation of $\int_{D'} |F^{c^*}(z, y_0)| dz$.

**Remark 3.4.** At this point the reader might wonder what is the necessity for a new technical definition for the focal spot which is a well known notion in the imaging community. This definition is actually just an extension of the definition used in the physics community where in practice, the focal spot is often characterized as the area inside the $-6dB$ level set around the peak of the point-spread function. The constraint on the $L^1$ norm is added here for technical reasons, to be able to approximate quantities of the type $\int_{D'} F^{c^*}(z, y) f(z) dz$ by $\int_{D'_\delta(y)} F^{c^*}(z, y) f(z) dz$ for some various $f$.

**Remark 3.5.** In practice, for $\delta \approx 0.1$ the focal spot's sizes are approximatively the Rayleigh's resolution limit and behave like $\frac{c^*}{\omega_c}$. More precisely, in the broadband regime, in the case of a linear array probe the transverse size of the focal spot at $y_0 \in D$ is given by:

$$\frac{2c^* y_0''}{\omega_c \ell}$$

and its axial size is

$$\frac{2c^*}{B'}$$

where $B$ is the bandwidth and $y_0''$ is the numerical aperture at depth $y_0''$. See figure 4.

3.3.2. Focal spot in the case $c \neq c^*$. In the case where $c \neq c^*$ the function

$$D' \to \mathbb{C} \\ z \mapsto F^{c^*}(z, y_0)$$

still has the same qualitative properties as in the previous case : there exists a small region of $D'$ where most of the signal concentrates.

We can still define a notion of focal spot:
**Definition 3.4** (Focal spot at threshold $\delta$ and speed $c$). Consider $0 < \delta \ll 1$ and $y_0 \in D$. There exists a domain $D_\delta^c(\varphi_c(y_0))' \subset D'$ such that

$$ \int_{D' \setminus D_\delta^c(\varphi_c(y_0))'} |F^c(z, y_0)| \, dz < \delta \int_{D'} |F^c(z, y_0)| \, dz, $$

where $\varphi_c$ is defined in (15).

**Remark 3.6.** A simple interpretation for the focal spot is the following: for a point $y_0 \in D$ and a given back-propagation speed of sound $c$, the focal spot $D_\delta'$ indicates the area in the image $D'$ where the contribution of point $y_0$ will appear.

The well-posedness of that definition is a direct consequence of the following lemma:

**Lemma 3.1** (Point spread function in the paraxial approximation). Let $(z, y) \in D' \times D$ satisfying the paraxial approximation, i.e consider $0 < \eta \ll 1$ and

$$ y = (\eta^{\frac{1}{2}} y^{\perp}, y^{\parallel}) \in D, \quad \text{and} \quad z = (\eta^{\frac{1}{2}} z^{\perp}, z^{\parallel}) \in D'. $$

In the paraxial regime in dimension 3 the point spread function has the expression:

$$ F^c(z, y) = \eta^{-1} \left(\frac{c}{c^*}\right)^2 \frac{\ell_0^4}{(16\pi^2|z||\varphi_c(y)|)^2} \int_{B_0} \omega^2 e^{i \frac{2\omega}{\eta c^*} (|z| - |\varphi_c(y)|)} e^{i \frac{\omega}{c^*} \left(1 - \left(\frac{c^*}{c}\right)^2\right) \frac{|\varphi_c(y)^\perp|^2}{|\varphi_c(y)|}} \\ G^2 \left( \frac{\omega \ell_0}{c} \left( \frac{z^\perp}{|z|} - \frac{\varphi_c(y)^\perp}{|\varphi_c(y)|} \right), \frac{\omega \ell_0^2}{c} \left( \frac{1}{|z|} - \left(\frac{c}{c^*}\right)^2 \frac{1}{|\varphi_c(y)|} \right) \right) d\omega + \mathcal{O}(1), \quad (13) $$

or equivalently

$$ F^c(z, y) = \eta^{-1} \left(\frac{c^*}{c}\right)^2 \frac{\ell_0^4}{(16\pi^2|y||\varphi_c^{-1}(z)|)^2} \int_{B_0} \omega^2 e^{i \frac{2\omega}{\eta c^*} (|\varphi_c^{-1}(z)| - |y|)} e^{i \frac{\omega}{c^*} ((\frac{c}{c^*})^2 - 1) \frac{|\varphi_c^{-1}(z)|^2}{|\varphi_c^{-1}(z)|}} \\ G^2 \left( \frac{\omega \ell_0}{c^*} \left( \frac{\varphi_c^{-1}(z)^\perp}{|\varphi_c^{-1}(z)|} - \frac{y^\perp}{|y|} \right), \frac{\omega \ell_0^2}{c^*} \left( \left(\frac{c^*}{c}\right)^2 \frac{1}{|\varphi_c^{-1}(z)|} - \frac{1}{|y|} \right) \right) d\omega + \mathcal{O}(1), \quad (14) $$

where

$$ \mathcal{G}(\xi_1, \xi_2) := \int_{[-1,1]^2} e^{-ix_e^\perp \cdot \xi_1 + i \frac{|x_e^\perp|^2}{2}} \xi_2 dx_e^\perp \quad \xi_1 \in \mathbb{R}^2, \xi_2 \in \mathbb{R}, $$

and

$$ D \to D' \\ y \mapsto \varphi_c(y) := \left( \eta^{\frac{1}{2}} \left( \frac{c}{c^*} \right)^2 y^{\perp}, \frac{c}{c^*} y^{\parallel} \right). \quad (15) $$

**Remark 3.7.** In dimension 2, (13) and (14) differ only by the following respective prefactors in the integrand

$$ \frac{32\pi \sqrt{|z||\varphi_c(y)|cc^*}}{\omega}, \qquad \frac{32\pi \sqrt{|y||\varphi_c^{-1}(z)|cc^*}}{\omega}. $$

and $\mathcal{G}$:

$$ \mathcal{G}(\xi_1, \xi_2) := \int_{[-1,1]^2} e^{-ix_e^\perp \cdot \xi_1 + i \frac{|x_e^\perp|^2}{2}} \xi_2 dx_e^\perp \quad \xi_1 \in \mathbb{R}, \xi_2 \in \mathbb{R}. $$

*Proof.* We suppose that $x_e := (\eta^{\frac{1}{2}} x_e^\perp, 0)$ and that $z$ is in the paraxial regime $z = (\eta^{\frac{1}{2}} z^\perp, z^{\parallel})$. We have:

$$ |x_e - z| = |z| \left( 1 + \frac{\eta |x_e^\perp|^2}{2|z|^2} - \eta \frac{\langle z^\perp, x_e^\perp \rangle}{|z|^2} \right) + \mathcal{O}(\eta^2). $$
For a point $y := (\eta^{\frac{1}{2}}y^{\perp}, y^{\parallel})$:

$$ \frac{c^*}{c}|x_e - z| - |x_e - y| = \left(\frac{c^*}{c}|z| - |y|\right) - \eta \left(\frac{c^*}{c} \frac{\langle z^{\perp}, x_e^{\perp} \rangle}{|z|} - \frac{\langle y^{\perp}, x_e^{\perp} \rangle}{|y|}\right) + \frac{\eta}{2} |x_e^{\perp}|^2 \left(\frac{c^*}{c} \frac{1}{|z|} - \frac{1}{|y|}\right) + \mathcal{O}(\eta^2). $$

Since the point-spread function is a peaked function we want to find the center of the focal spot that we denote $\varphi_c(y) \in D'$ and express the point spread function with respect to $z - \varphi_c(y)$. We can find this center $\varphi_c(y)$ by noting that since the frequency behaves as $\eta^{-1}$, a necessary condition is that if $z = \varphi_c(y)$ then the term in $\eta^0$ in the phase must cancel. Cancelling the zero order term gives the first necessary condition for $\varphi_c(y)$

$$ |\varphi_c(y)| = \frac{c}{c^*} |y| + \mathcal{O}(\eta). \qquad (16) $$

This gives, since $y$ is in the paraxial regime:

$$ \varphi_c(y)^{\parallel} := \frac{c}{c^*} y^{\parallel}. $$

We are now looking for a necessary condition on the center of the focal spot $\varphi_c(y)$'s transverse coordinate $\varphi_c(y)^{\perp}$. Integrating over the probe and the bandwidth, for $z = \varphi_c(y)$

$$ F^c(z, y) = \eta^{-1} \left(\frac{c^*}{c}\right)^2 \frac{\ell_0^4}{(16\pi^2|y||\varphi_c(y)|)^2} \int_{B_0} \omega^2 e^{i\frac{\omega}{c}\left(1-\left(\frac{c^*}{c}\right)^2\right)\frac{|\varphi_c(y)^{\perp}|^2}{|\varphi_c(y)|}} G^2 \left( \frac{\omega\ell_0}{c^*} \left( \left(\frac{c^*}{c}\right)^2 \frac{\varphi_c(y)^{\perp}}{|y|} - \frac{y^{\perp}}{|y|} \right), \frac{\omega\ell_0^2}{c^*|y|} \left( \left(\frac{c^*}{c}\right)^2 - 1 \right) \right) d\omega + \mathcal{O}(1), $$

where

$$ G(\xi_1, \xi_2) := \int_{[-1,1]^2} e^{-ix_e^{\perp} \cdot \xi_1 + i \frac{|x_e^{\perp}|^2}{2} \xi_2} d x_e^{\perp} \quad \xi_1 \in \mathbb{R}^2, \xi_2 \in \mathbb{R}. $$

Since $G$ is a peaked function at the origin for both arguments, this gives us the transverse position of the center of the focal spot by setting:

$$ \left(\frac{c^*}{c}\right)^2 \frac{\varphi_c(y)^{\perp}}{|y|} - \frac{y^{\perp}}{|y|} = 0 $$

and we get the second necessary condition:

$$ \langle \varphi_c(y)^{\perp}, x_e^{\perp} \rangle = \left(\frac{c}{c^*}\right)^2 \langle y^{\perp}, x_e^{\perp} \rangle. \qquad (17) $$

We therefore set

$$ \varphi_c(y) := \left( \eta^{\frac{1}{2}} \left( \frac{c}{c^*} \right)^2 y^{\perp}, \frac{c}{c^*} y^{\parallel} \right). $$

We compute:

$$ |\varphi_c(y)| - \frac{c}{c^*} |y| = \frac{\eta}{2} \left(1 - \left(\frac{c^*}{c}\right)^2\right) \frac{|\varphi_c(y)^{\perp}|^2}{|\varphi_c(y)|} $$

Going back to the PSF, after a change of variable in $\omega$, the following asymptotic expression when $\eta \to 0$:
$$
F^c(z, y) = \eta^{-1} \left(\frac{c^*}{c}\right)^2 \frac{\ell_0^4}{(16\pi^2|y|\phi_c^{-1}(z)|)^2} \int_{B_0} \omega^2 e^{i \frac{2\omega}{c^*} (|\phi_c^{-1}(z)| - |y|)} e^{i \frac{\omega}{c^*} \left(\left(\frac{c}{c^*}\right)^2 - 1\right) \frac{|\phi_c^{-1}(z)|^2}{|\phi_c^{-1}(z)|}} \\
G^2 \left( \frac{\omega \ell_0}{c^*} \left( \frac{\phi_c^{-1}(z)^\perp}{|\phi_c^{-1}(z)|} - \frac{y^\perp}{|y|} \right), \frac{\omega \ell_0^2}{c^*} \left( \left( \frac{c^*}{c} \right)^2 \frac{1}{|\phi_c^{-1}(z)|} - \frac{1}{|y|} \right) \right) d\omega + \mathcal{O}(1).
$$

□

**Corollary 3.1.** For $z = \varphi_c(y_0) + (\eta^{\frac{1}{2}} \Delta z^{\perp}, \eta \Delta z^{\parallel})$, the expression of the point spread function becomes

$$
F^c(z, y_0) = \eta^{-1} \left(\frac{c^*}{c}\right)^2 \frac{\ell_0^4}{(16\pi^2|y_0|^2)^2} \int_{B_0} \omega^2 e^{i \frac{\omega}{|y_0|c^*} |y_0|^2 \left(\left(\frac{c}{c^*}\right)^2 - 1\right)} e^{i \frac{2\omega}{c} \Delta z^{\parallel}} e^{i \frac{\omega}{|y_0|c^*} \left(\frac{c^*}{c}\right)^2 |\Delta z^{\perp}|^2} \\
G^2 \left( \frac{\omega \ell_0}{c^* |y_0|} \left(\frac{c^*}{c}\right)^2 |\Delta z^{\perp}|, \frac{\omega \ell_0^2}{c^* |y_0|} \left(\left(\frac{c^*}{c}\right)^2 - 1\right) \right) d\omega + \mathcal{O}(1). \quad (18)
$$

If $B_0 < \omega_0$ we are in the so-called broadband regime(see section 2.1 and [14, chapter 6]) where the expression above simplifies to:

$$
F^{c^*}(z, y) = \eta^{-1} \left(\frac{c^*}{c}\right)^2 \frac{\omega_0^2 \ell_0^4}{(16\pi^2 |y_0|^2)^2} e^{i \frac{2\omega_0}{c} \Delta z^{\parallel}} e^{i \frac{\omega_0}{|y_0|c^*} |y_0|^2 \left(\left(\frac{c}{c^*}\right)^2 - 1\right)} e^{i \frac{\omega_0}{|y_0|c^*} \left(\frac{c^*}{c}\right)^2 |\Delta z^{\perp}|^2} \\
\text{sinc} \left( \frac{B}{c^*} \left( \frac{c^*}{c} \Delta z^{\parallel} + \left(\frac{c^*}{c}\right)^2 \frac{1}{|y_0|} |\Delta z^{\perp}|^2 + \left(\left(\frac{c}{c^*}\right)^2 - 1\right) \frac{1}{|y_0|} |y_0^{\perp}|^2 \right) \right) \\
G^2 \left( \frac{\omega_0 \ell_0}{c^* |y_0|} \left(\frac{c^*}{c}\right)^2 |\Delta z^{\perp}|, \frac{\omega_0 \ell_0^2}{c^* |y_0|} \left(\left(\frac{c^*}{c}\right)^2 - 1\right) \right) d\omega + \mathcal{O}(1). \quad (19)
$$

*Proof.* Let $z := \varphi_c(y_0) + \Delta z$ with $\Delta z := (\eta^{\frac{1}{2}} \Delta z^{\perp}, \eta \Delta z^{\parallel})$. We can write:

$$
|z| = |\varphi_c(y_0)| \left( 1 + \frac{\eta}{2} \frac{|\Delta z^{\perp}|^2}{|\varphi_c(y_0)|^2} + \eta \frac{\Delta z^{\parallel}}{|\varphi_c(y_0)|} \right) + \mathcal{O}(\eta^{\frac{3}{2}}).
$$

and then:

$$
\frac{c^*}{c} |z| - |y_0| = \frac{c^*}{c} |\varphi_c(y_0)| - |y_0| + \frac{\eta}{2} \frac{c^*}{c} \frac{|\Delta z^{\perp}|^2}{|\varphi_c(y_0)|} + \eta \frac{c^*}{c} \Delta z^{\parallel} + \mathcal{O}(\eta^{\frac{3}{2}}).
$$

Using (15) we get

$$
\frac{c^*}{c} |z| - |y_0| = \frac{\eta}{2} \frac{|y_0^\perp|^2}{|y_0|} \left( \left( \frac{c}{c^*} \right)^2 - 1 \right) + \frac{\eta}{2} \left( \frac{c^*}{c} \right)^2 \frac{|\Delta z^\perp|^2}{|y_0|} + \eta \frac{c^*}{c} \Delta z^{\parallel} + \mathcal{O}(\eta^{\frac{3}{2}}).
$$

Similarly

$$
\frac{c^*}{c} \frac{1}{|z|} - \frac{1}{|y_0|} = \frac{1}{|y_0|} \left( \left( \frac{c^*}{c} \right)^2 - 1 \right) + \mathcal{O}(\eta^{\frac{1}{2}}).
$$

We can also write
$$ \frac{c^*}{c} \frac{\langle z^{\perp}, x_e^{\perp} \rangle}{|z|} - \frac{\langle y_0^{\perp}, x_e^{\perp} \rangle}{|y_0|} = \frac{1}{|y_0|} \left( \left( \frac{c^*}{c} \right)^2 \langle z^{\perp}, x_e^{\perp} \rangle - \langle y_0^{\perp}, x_e^{\perp} \rangle \right) + \mathcal{O}(\eta) \\ = \frac{1}{|y_0|} \left( \left( \frac{c^*}{c} \right)^2 \langle \varphi_c(y_0)^{\perp}, x_e^{\perp} \rangle - \langle y_0^{\perp}, x_e^{\perp} \rangle \right) \\ + \frac{1}{|y_0|} \left( \frac{c^*}{c} \right)^2 \langle \Delta z^{\perp}, x_e^{\perp} \rangle + \mathcal{O}(\eta^{\frac{1}{2}}). $$

Using (15) again we get:

$$ \frac{c^*}{c} \frac{\langle z^{\perp}, x_e^{\perp} \rangle}{|z|} - \frac{\langle y_0^{\perp}, x_e^{\perp} \rangle}{|y_0|} = \frac{1}{|y_0|} \left( \frac{c^*}{c} \right)^2 \langle \Delta z^{\perp}, x_e^{\perp} \rangle + \mathcal{O}(\eta^{\frac{1}{2}}). $$

Finally we obtain the following asymptotic expansion for the phase term in the imaging domain i.e. at point $z = \varphi_c(y_0) + \Delta z \in D'$, with $\Delta z := (\eta^{\frac{1}{2}} \Delta z^{\perp}, \eta \Delta z^{\parallel})$ and $y_0 \in D$:

$$ \frac{c^*}{c} |x_e - z| - |x_e - y_0| = \frac{\eta}{2} \frac{|y_0^{\perp}|^2}{|y_0|} \left( \left( \frac{c}{c^*} \right)^2 - 1 \right) + \frac{\eta}{2} \left( \frac{c^*}{c} \right)^2 \frac{|\Delta z^{\perp}|^2}{|y_0|} + \eta \frac{c^*}{c} \Delta z^{\parallel} \\ - \frac{\eta}{|y_0|} \left( \frac{c^*}{c} \right)^2 \langle \Delta z^{\perp}, x_e^{\perp} \rangle + \frac{\eta}{2} \frac{|x_e^{\perp}|^2}{|y_0|} \left( \left( \frac{c^*}{c} \right)^2 - 1 \right) + \mathcal{O}(\eta^{\frac{3}{2}}). $$

Plugging in the phase term in the expression for the point spread function we get that (14) becomes

$$ F^c(z, y_0) = \eta^{-1} \left(\frac{c^*}{c}\right)^2 \frac{\ell_0^4}{(16\pi^2(y_0^{\parallel})^2)^2} \int_{\mathcal{B}_0} \omega^2 e^{i \frac{\omega}{|y_0|c^*} |y_0^{\perp}|^2} \left(\left(\frac{c}{c^*}\right)^2 - 1\right) e^{i \frac{2\omega}{c} \Delta z^{\parallel}} e^{i \frac{\omega}{|y_0|c^*} \left(\frac{c^*}{c}\right)^2 |\Delta z^{\perp}|^2} \\ G^2 \left( \frac{\omega \ell_0}{c^* |y_0|} \left(\frac{c^*}{c}\right)^2 \Delta z^{\perp}, \frac{\omega \ell_0^2}{c^* |y_0|} \left(\left(\frac{c^*}{c}\right)^2 - 1\right) \right) d\omega + \mathcal{O}(1). $$

□

**3.4. Discussion.** We discuss here (18) and (21). For $c = c^*$ we recover the classical point spread function of a linear array imaging system in the paraxial regime. In the so-called *broadband regime* [14, chapter 6], for $c = c^*$, (21) gives us the classical expected resolution shown on the middle insert of figure 4 (see remark 3.5).

The use of an incorrect backpropagation speed $c \neq c^*$ has three main effects on the focal spot:

(1) The focal spot is not centered around $y_0$ anymore but around the image of $y_0$ by the mapping between the physical object $D$ and its image $D'$:

$$ D \rightarrow D' \\ y_0 \mapsto \varphi_c(y_0). $$

This is a first order effect in the sense that the center of the focal spot $\varphi_c(y_0)$ moves proportionally to $\frac{c}{c^*} - 1$ in the axial direction.

(2) The amplitude at the center of the focal spot is decreased. This can be seen easily by looking at the amplitudes on fig. 4. For a more precise quantification of this effect one can look at the expression (18) at the center of the focal spot (with $z = \varphi_c(y)$) and see that the amplitude is proportional to

$$ \int_{\mathcal{B}_0} \omega^2 e^{i \frac{\omega |y_0^{\perp}|^2}{|y_0|c^*}} \left( \left( \frac{c}{c^*} \right)^2 - 1 \right) G^2 \left( 0, \frac{\omega \ell_0^2}{c^* |y_0|} \left( \left( \frac{c^*}{c} \right)^2 - 1 \right) \right) d\omega. $$

This is a notable effect as the sensitivity to $\left(\frac{c^*}{c}\right)^2 - 1$ is quite high. With realistic values corresponding to practical situations in medical ultrasound imaging we can see that a 15% error on the speed of sound can lead to an order of magnitude decrease in the amplitude at the center, see fig. 3.
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Plot of $(\xi_1, \xi_2) \mapsto |\mathcal{G}(\xi_1, \xi_2)|$.</figcaption></figure>

(3) The shape of the focal spot is modified, it is more spread out and does not have unique local maxima. However looking at (18) gives us that $z \mapsto F^c(z, y_0)$ is of order $\eta^{-1}$ in a region of typical size of the order $\eta$ in the axial direction and $\eta^{1/2}$ in the transverse directions which is the same regime as the case $c = c^*$. The spreading out is due to the shape of $\xi_1 \mapsto \mathcal{G}(\xi_1, \xi_2)$ which presents a sharper peak (more contrasted and with a smaller width at half maximum) for $\xi_2$ close to 0. Even though this effect is more visually striking, it is harder to quantify compared to the effect described in the second point. See [24].

For a visual representation of these effects we show an illustration on figure 4 where we plot the point spread function $z \mapsto F^c(z, y_0)$ in three different cases : $c < c^*, c = c^*, $ and $c > c^*$.

**3.5. Focusing criterion.** Consider a point in the medium $y_0 \in D$. Recall that for each $c$ the center of the focal spot in the image is at $\varphi_c(y_0) \in D'$. Let us define the function

**Definition 3.5 (Confocal trace at $y_0$).**

$$c \mapsto \mathcal{F}_{y_0}(c) := F^c(\varphi_c(y_0), y_0). \qquad (20)$$

**Remark 3.8.** Introduced this way, the function $\mathcal{F}_{y_0}$ is an abstract construction but it can be understood more easily if one considers the following thought experiment. Assume that there is a point-like single isolated reflector located at $y_0$ in an unknown homogeneous medium whose acoustic property $c^*$ is to be determined. In this case the imaging function has a simple approximation:

$$I^c(z) \propto F^c(z, y_0).$$

Up to a horizontal translation of the probe we can assume that this single reflector is located on the axis $y_0^\perp = 0$. Of course, since we do not know the propagation speed inside the medium, the only thing we can access experimentally is the travel time from the reflector to the center of the probe that we can denote $t_0$ (via the D.O.R.T method for example [29]). But that travel time $t_0$ is all one needs to plot the function $c \mapsto I^c(z(c))$ with $z(c) = (0, ct_0)$. By noting that
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Normalized plot of $c \mapsto |\int_B \omega^2 G^2 \left(0, \frac{\omega \ell^2}{c^*|y_0|} \left(\frac{c^*}{c}\right)^2 - 1\right) d\omega|$ for realistic values of the parameters $\ell, |y_0|, \mathcal{B}, \dots$ (see section 5).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Shape of the point spread function $z \mapsto |F^c(z, y_0)|$ with three different backpropagation speed. Here $y_0 = (0.01, 0.06)$ for all three graphs. The quantity $y_0^\ell$ is the numerical aperture at $y_0$ where $\ell$ is the length of the linear probe $\mathcal{P}$. The focal spot at threshold 0.1 is the area in yellow/light blue.</figcaption></figure>

$z(c)$ is exactly $\phi_c(y_0)$ one can see that $c \mapsto I^c(z(c))$ corresponds to the value of the pixel at the center of the focal spot for every possible value of $c$, which is exactly our confocal trace $\mathcal{F}$ introduced i.e.

$$I^c(z(c)) = \mathcal{F}(c).$$

The speed of sound $c^*$ can then be recovered using the estimator introduced in the proposition 3.1 that follows. Once $c^*$ is recovered, the actual position $z^*$ of the reflector can also be recovered.
**Proposition 3.1 (Amplitude characterization).**

$$c^* = \operatorname{argmax}_c |\mathcal{F}_{y_0}(c)|.$$

*Proof.* By taking $\Delta z = 0$ in equation (18) we get:

$$\mathcal{F}_{y_0}(c) = \eta^{-1} \left(\frac{c^*}{c}\right)^2 \frac{\ell_0^4}{(16\pi^2(y_0^{\parallel})^2)^2} \int_{\mathcal{B}_0} \omega^2 e^{i\frac{\omega}{y_0^{\parallel}c^*}|y_0^{\parallel}|^2 ((\frac{c}{c^*})^2 - 1)} \mathcal{G}^2 \left(0, \frac{\omega \ell_0^2}{c^* y_0^{\parallel}} \left(\left(\frac{c^*}{c}\right)^2 - 1\right)\right) d\omega + \mathcal{O}(1).$$

With similar arguments as in the second point of discussion 3.4 above we recognize a peaked function and see that $c \mapsto \mathcal{F}_{y_0}(c)$ has a maximum at $c = c^*$. It is easier to see in the so-called *broadband regime* [14, chapter 6] ($\mathcal{B}_0 < \omega_0$) where the expression above simplifies to:

$$\mathcal{F}_{y_0}(c) = \eta^{-1} \left(\frac{c^*}{c}\right)^2 \frac{\ell_0^4 \omega_0^2}{(16\pi^2(y_0^{\parallel})^2)^2} \text{sinc} \left( \frac{B}{y_0^{\parallel} c^*} |y_0^{\parallel}|^2 \left( \left(\frac{c}{c^*}\right)^2 - 1 \right) \right) e^{i \frac{\omega_0}{y_0^{\parallel} c^*} |y_0^{\parallel}|^2 \left( \left(\frac{c}{c^*}\right)^2 - 1 \right)} \mathcal{G}^2 \left( 0, \frac{\omega_0 \ell_0^2}{c^* y_0^{\parallel}} \left( \left(\frac{c^*}{c}\right)^2 - 1 \right) \right). \quad (21) \quad \square$$

**Remark 3.9** (About the resolution of the estimator). Looking at (21) we can see that $c \mapsto \mathcal{F}_{y_0}(c)$ has a peak at $c = c^*$ and that the width at half-maximum in the broadband regime depends on the numerical aperture at the depth of the focusing point $y_0$, the central frequency $\omega_c$ and is given by:

$$\frac{2c^*y_0^{\parallel}}{\omega_c \ell^2}.$$

## 4. SPEED OF SOUND ESTIMATION IN THE TISSUE-MIMICKING MEDIUM

**4.1. The dual focal spot.** With simple symmetry considerations and by exchanging the role of $c$ and $c^*$ we can see that for a fixed point $z_0 \in D'$ the function

$$\begin{aligned} D &\to D \\ y &\mapsto F^c(z_0, y) \end{aligned}$$

will have the same qualitative properties as $z \mapsto F^c(z, y_0)$. It is then natural to introduce the:

**Definition 4.1** (Dual focal spot at threshold $\delta$ and speed $c$). Consider $0 < \delta \ll 1$ and $z_0 \in D'$. Following the same construction as the focal spot, there exists a domain $D_{\delta}^c(\varphi_c^{-1}(z_0)) \subset D$ such that

$$\int_{D \setminus D_{\delta}^{c}(\varphi_{c}^{-1}(z_{0}))} |F^{c}(z_{0}, y)| \, dy < \delta \int_{D} |F^{c}(z_{0}, y)| \, dy.$$

**Remark 4.1.** Similarly to the focal spot, the dual focal spot is not centered at $z_0$ but at $\varphi_c^{-1}(z_0)$.

**Remark 4.2.** A simple interpretation for the dual focal spot is the following : for a point in the image $z_0 \in D'$ and a given back-propagation speed of sound $c$, the dual focal spot $D_{\delta}^c(\varphi_c^{-1}(z_0))$ indicates the area in the medium $D$ that mostly contributes to the value of pixel $I^c(z_0)$ in the image.

**4.2. Locally probing the micro-structure.** In light of the previous definition, the imaging function depends on the material's micro-structure just in the dual focal spot area :

$$\begin{aligned} I^c(z) &= \int_D (n_\varepsilon(y) - n^*) F^c(z, y) dy \\ &= \int_{D_{\delta}^c(\varphi_c^{-1}(z_0))} (n_\varepsilon(y) - n^*) F^c(z, y) dy + \mathcal{O}(\delta). \end{aligned}$$

This observation is the first step in the construction of the virtual guide star.
The second step is based on the following observation : for a fixed point satisfying the paraxial approximation $z_0 = (\eta^{\frac{1}{2}}z_0^{\perp}, z_0^{\parallel})$, $c \mapsto I^c(z_0)$ probes the moving area $D_{\delta}^c(\varphi_c^{-1}(z_0))$ in the physical domain $D$. Its center (as a function of $c$) is given by

$$c \mapsto \varphi_c^{-1}(z_0) = \left( \left( \frac{c^*}{c} \right)^2 \eta^{\frac{1}{2}} z_0^{\perp}, \frac{c^*}{c} z_0^{\parallel} \right)$$

and its size increases with $|c - c^*|$. Therefore, for a fixed $y_0 \in D$ the function

$$c \mapsto I^c(\varphi_c(y_0))$$

probes the area $D_{\delta}^c(y_0)$, always centered at $y_0$. Of course, its size is still dependent on $|c - c^*|$.

In practice, $\varphi_c$ is unknown as it depends on $c$ (known) and $c^*$ (unknown). However, as done in [7], let us define, for $s_0 \in \mathbb{R}^{d-1}$ and $t_0 \in \mathbb{R}^+$:

$$\psi_c(s_0, t_0) := (c^2 s_0, ct_0).$$

Now if we compute

$$c \mapsto I^c(\psi_c(s_0, t_0))$$

we are in fact computing $c \mapsto I^c(\varphi_c(y_0))$ with an (unknown but fixed)

$$y_0 := \left( \eta^{\frac{1}{2}} (c^*)^2 s_0, c^* t_0 \right).$$

**Remark 4.3.** For $s_0 = 0$, similarly as in remark 3.8, $t_0$ corresponds to the travel time from $y_0$ to the center of the probe at speed $c^*$.

From now on, we choose $s_0 \in \mathbb{R}^{d-1}$, $t_0 \in \mathbb{R}^+$ and we study the dependency with respect to $c$ of the following functional:

$$I^c(\varphi_c(s_0, t_0)) = I^c(\varphi_c(y_0)) \qquad (22)$$

$$= \int_{D_{\delta}^{c}(y_0)} (n_{\varepsilon}(y) - n^{*}) F^{c}(\varphi_{c}(y_0), y) dy + \mathcal{O}(\delta). \qquad (23)$$

Since $D_{\delta}^c(y_0)$ is a small domain centered at $y_0$ we can see that for $y \in D_{\delta}^c(y_0)$ we can expect $F^c(\varphi_c(y_0), y)$ to be close to the confocal trace $F_{y_0}(c) := F^c(\varphi_c(y_0), y_0)$. The difficulty is that $F^c(\varphi_c(y_0), y)$ is integrated against an unknown random, rapidly oscillating -at the scale $\varepsilon r_0 \ll \frac{c^*}{\omega_c}$ - index $n_{\varepsilon}(\cdot) - n^*$ on the dual focal spot. In the next section we show how to overcome this difficulty and build an estimator for the effective speed of sound in the random multi-scale medium.

**Proposition 4.1.** Imaging in the speckle. For $s_0 \in \mathbb{R}^{d-1}$, $t_0 > 0$, let $z(c) := (\eta^{\frac{1}{2}}c^2s_0, ct_0) := \varphi_c(y_0)$, the expression of the imaging function at $z(c)$ reads

$$I^c(\varphi_c(y_0)) = \eta \left(\frac{c}{c^*}\right)^2 \frac{\ell_0^4}{(16\pi^2 |\varphi_c(y_0)|^2)^2} \\ \int_{Y(D_{\delta}^c(y_0))} (n_{\varepsilon}(y_0 + (\tilde{y})_{\eta}) - n^*) \int_{B_0} \omega^2 e^{i \frac{2\omega}{c^*} |y_0|} \left(\frac{c}{c^*}\right)^2 \langle \tilde{y}^{\perp}, y_0^{\perp} \rangle e^{i \frac{2\omega}{c^*} \tilde{y}^{\parallel}} e^{i \frac{2\omega}{c^*} |y_0|} \left(\frac{c}{c^*}\right)^2 |\tilde{y}^{\perp}|^2 \\ e^{i \frac{\omega}{c^*} ((\frac{c}{c^*})^2 - 1) \frac{|y_0^{\perp} + \tilde{y}^{\perp}|^2}{|y_0|}} G^2 \left( \frac{\omega \ell_0}{c^* |y_0|} \left(\frac{c}{c^*}\right)^2 \tilde{y}^{\perp}, \frac{\omega \ell_0^2}{c^* |y_0|} \left(\left(\frac{c^*}{c}\right)^2 - 1\right) \right) d\omega d\tilde{y} + \mathcal{O}(\eta^2).$$

where $(\tilde{y})_{\eta} := y - y_0 := (\eta^{\frac{1}{2}}\tilde{y}^{\perp}, \eta\tilde{y}^{\parallel})$ and $Y(D_{\delta}^c(y_0))$ is the rescaled (by $\eta^{\frac{1}{2}}$ in the transverse directions and $\eta$ in the axial direction) and translated (by $-y_0$) dual focal spot.
**Remark 4.4.** By switching the roles of c and c\* we can see that the focal spot in the imaging domain and the dual focal spot in the physical domain will have similar sizes. The scaling of the dummy variable of integration (ŷ)η in proposition 4.1 confirms that conjecture and shows that the dual focal spot is indeed of size η in the axial direction and η 1/2 in the transverse direction. Exactly as the focal spot, the dual focal spot size and shape are sensitive to the mismatch between c and c\*, but its axial and transverse sizes stay in the same order of magnitude in the regime considered.

Proof. We know from Lemma 3.1 that:

$$
F^c(z, y) = \eta^{-1} \left(\frac{c}{c^*}\right)^2 \frac{\ell_0^4}{(16\pi^2|z||\varphi_c(y)|)^2} \int_{\mathcal{B}_0} \omega^2 e^{i \frac{2\omega}{\eta c} (|z|-|\varphi_c(y)|)} e^{i \frac{\omega}{c} \left(1-\left(\frac{c^*}{c}\right)^2\right) \frac{|\varphi_c(y)|^2}{|\varphi_c(y)|}} \\
G^2 \left( \frac{\omega \ell_0}{c} \left( \frac{z^\perp}{|z|} - \frac{\varphi_c(y)^\perp}{|\varphi_c(y)|} \right), \frac{\omega \ell_0^2}{c} \left( \frac{1}{|z|} - \left(\frac{c}{c^*}\right)^2 \frac{1}{|\varphi_c(y)|} \right) \right) + \mathcal{O}(1).
$$

We perform the change of variable $(\tilde{y})_\eta := y - y_0 := (\eta^{\frac{1}{2}} \tilde{y}^\perp, \eta \tilde{y}^{\parallel})$. We have

$$
|\varphi_c(y)| = |\varphi_c(y_0)| + \frac{\eta}{2} \frac{|\varphi_c(\tilde{y})^{\perp}|^2}{|\varphi_c(y_0)|} + \eta \varphi_c(y)^{\parallel} + \eta \frac{\langle \varphi_c(\tilde{y})^{\perp}, \varphi_c(y_0)^{\perp} \rangle}{|\varphi_c(y_0)|} + \mathcal{O}(\eta^{\frac{3}{2}}).
$$

We get

$$
|z| - |\varphi_c(y)| = \frac{\eta}{2} \frac{|\varphi_c(\tilde{y})^{\perp}|^2}{|\varphi_c(y_0)|} + \eta \frac{c}{c^*} \tilde{y}^{\parallel} + \eta \frac{\langle \varphi_c(\tilde{y})^{\perp}, \varphi_c(y_0)^{\perp} \rangle}{|\varphi_c(y_0)|} + \mathcal{O}(\eta^{\frac{3}{2}}),
$$
$$
\frac{z^{\perp}}{|z|} - \frac{\varphi_c(y)^{\perp}}{|\varphi_c(y)|} = \frac{\varphi_c(\tilde{y})^{\perp}}{|\varphi_c(y_0)|} + \mathcal{O}(\eta),
$$

and

$$
\left(\frac{c}{c^*}\right)^2 \frac{1}{|z|} - \frac{1}{|\varphi_c(y)|} = \frac{1}{|\varphi_c(y_0)|} \left( \left(\frac{c}{c^*}\right)^2 - 1 \right) + \mathcal{O}(\eta).
$$

By plugging those asymptotic expansions back in (22) we get the result.

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

4.3. Effective speed of sound estimation via ensemble averaging. In thise section we show that it is possible
to access the spatial average of the point-spread function over the dual focal spot via the variance of the
imaging function (lemma 4.1). Using a similar idea as what was developped in section 3.5, we then build an
estimator for the effective sound velocity in the medium (proposition 4.2).

**Lemma 4.1** (Second order moment). Assume $n: \mathbb{R}^d \times \Omega \to [n_-, n_+]$ is a stationary and mixing process satisfying (6). Then as $\eta \to 0$ we have:

$$
\mathbb{E} \left[ \left| \int_{D_{\delta}^{c}(\varphi_{c}^{-1}(z))} (n_{\varepsilon}(y) - n^{*}) F^{c}(z, y) dy \right|^{2} \right] = \varepsilon^{d} \int_{\mathbb{R}^{d}} C(s) ds \int_{D_{\delta}^{c}(\varphi_{c}^{-1}(z))} |F^{c}(z, y)|^{2} dy + o \left( \eta^{\frac{d-1}{2}} \right),
$$

where C is the covariance of n.

Proof. By Fubini’s theorem:

$$
\begin{align*}
A(z) &:= \mathbb{E} \left[ \left| \int_{D_{\delta}^{c}(\varphi_c^{-1}(z))} (n_{\varepsilon}(y) - n^{*}) F^{c}(z, y) dy \right|^2 \right] \\
&= \iint_{D_{\delta}^{c}(\varphi_c^{-1}(z))^2} \mathbb{E}[(n_{\varepsilon}(y) - n^{*})(n_{\varepsilon}(\tilde{y}) - n^{*})] F^{c}(z, y) \bar{F}^{c}(z, \tilde{y}) dy d\tilde{y} \\
&= \iint_{D_{\delta}^{c}(\varphi_c^{-1}(z))^2} C\left(\frac{y - \tilde{y}}{\varepsilon}\right) F^{c}(z, y) \bar{F}^{c}(z, \tilde{y}) dy d\tilde{y}
\end{align*}
$$
We change variables to get

$$A(z) = \varepsilon^d \iint_{\mathcal{D}_{\delta}^{c}(\varphi_c^{-1}(z)) \times \mathcal{D}_x} C(x) F^c(z, \tilde{y} + \varepsilon x) \overline{F^c}(z, \tilde{y}) d\tilde{y} dx,$$

where $\mathcal{D}_x := \left\{ \frac{y - \tilde{y}}{\varepsilon}, y, \tilde{y} \in \mathcal{D}_{\delta}^{c}(\varphi_c^{-1}(z)) \right\}$. We want to do a Taylor expansion in $F^c(z, \tilde{y} + \varepsilon x)$ with respect to $\varepsilon x$ but $|x|$ can be arbitrarily large at first glance. However, we can use the separation of scale in our model. Recall that the micro-structure varies at a scale of order $\varepsilon \ll \eta$ and the point-spread function varies at the scale of the wavelength $c\omega_{\eta}^{-1} = \eta^{-1} \frac{c}{\omega}$. More precisely, since $C$ is controlled by a decreasing function $\Phi$ there exists $R > 0$ such that

$$\left| \int_{\mathbb{R}^d} C - \int_{B(0,R)} C \right| \leq \varepsilon.$$

Moreover, the decay rate of $\Phi$ gives us $\varepsilon R = o(\eta)$, see lemma A.1. After a Taylor's expansion of $F^c$ at $\tilde{y}$ we obtain, at first order in $\eta$:

$$A(z) = \varepsilon^d \int_{\mathbb{R}^d} C(s) ds \int_{\mathcal{D}_{\delta}^{c}(\varphi_c^{-1}(z))} |F^c(z, \tilde{y})|^2 d\tilde{y} + o\left(\eta^{\frac{d-1}{2}}\right).$$

□

We are now ready to state our main proposition:

**Proposition 4.2 (Effective speed of sound estimator).** For any $s_0 \in \mathbb{R}^{d-1}$, $t_0 \in \mathbb{R}^+$:

$$c^* = \underset{c}{\operatorname{argmax}} \mathbb{E} \left[ |I^c(\psi_c(s_0, t_0))|^2 \right].$$

*Proof.* The typical sizes of $\mathcal{D}_{\delta}^{c}(\varphi_c^{-1}(z))$ are $\eta^{\frac{1}{2}} \times \eta^{-1}$ in dimension 2 and $\eta^{\frac{1}{2}} \times \eta^{\frac{1}{2}} \times \eta^{-1}$ in dimension 3. Therefore $\int_{\mathcal{D}_{\delta}^{c}(\varphi_c^{-1}(z))} |F^c(z, y)|^2 dy$ can be reasonably approximated via a Taylor expansion by $|D_{\delta}^{c}(\varphi_c^{-1}(z))| |F^c(z, \varphi_c^{-1}(z))|^2$ and:

$$\mathbb{E} \left[ |I^c(z, \cdot)|^2 \right] = \varepsilon^d \int_{\mathbb{R}^d} C(s) ds \left| D_{\delta}^{c}(\varphi_c^{-1}(z)) \right| \left| F^c(z, \varphi_c^{-1}(z)) \right|^2 + o \left( \eta^{\frac{d-1}{2}} \right).$$

Using $z = \psi_c(s_0, t_0)$ gives us $\varphi_c^{-1}(z) = y_0$ with

$$y_0 := \left( \eta^{\frac{1}{2}} (c^*)^2 s_0, c^* t_0 \right),$$

which does not depend on $c$ and therefore

$$\mathbb{E} \left[ |I^c(z, \cdot)|^2 \right] = \varepsilon^d \int_{\mathbb{R}^d} C(s) ds \left| D_{\delta}^{c}(y_0) \right| \left| \mathcal{F}_{y_0}(c) \right|^2 + o \left( \eta^{\frac{d-1}{2}} \right).$$

In light of remark 4.4 that $\mathbb{E} \left[ |I^c(z, \cdot)|^2 \right]$ therefore behaves qualitatively as $|\mathcal{F}_{y_0}(c)|^2$. □

**4.4. From ensemble averaging to spatial averaging.** In the previous subsection we constructed an estimator for the effective speed of sound in the micro-structured random medium. The computation of the estimator requires access to multiple realizations of the imaging functional to approximate its expected value. However, in practice, one has access to only a single realization of the medium. In this subsection, we show that we can approximate the expected value with data collected from only one realization of the medium.

In what follows we precise the dependency with respect to the realization in the notations. For a given realization $\omega \in \Omega$ we write $n(x, \omega)$, $x \in D$ and $I^c(z, \omega)$, $z \in D'$. We start by proving a technical lemma. In the case where $c = c^*$ the point spread function $F^{c^*}(z, y)$ is actually a function of $z - y$. So it is obvious that for any $\Delta z \in \mathbb{R}^d$ we have $F^{c^*}(z + \Delta z, y) = F^{c^*}(z, y - \Delta z)$. In the case $c \neq c^*$ we prove here a modified version of this result for a small (of order $\eta$) $\Delta z$.
**Lemma 4.2.** For $y_0 \in D$ and $y \in \mathcal{D}_\delta(y_0 + \varphi_c^{-1}((\Delta z)_\eta))$ with $(\Delta z)_\eta := (\eta \Delta z^\perp, \eta \Delta z^{\parallel}) \in \mathbb{R}^d$

$$F^c(\varphi_c(y_0) + (\Delta z)_\eta, y) = F^c(\varphi_c(y_0), y - \varphi_c^{-1}((\Delta z)_\eta)) + \mathcal{O}(\eta^{-\frac{1}{2}}).$$

*Proof.* Let us prove the result in dimension $d = 3$. We know from Lemma 3.1

$$F^c(z, y) = \eta^{-1} \left(\frac{c^*}{c}\right)^2 \frac{\ell_0^4}{(16\pi^2 |y| |\varphi_c^{-1}(z)|)^2} \int_{\mathcal{B}_0} \omega^2 e^{i \frac{2\omega}{\eta c^*} (|\varphi_c^{-1}(z)| - |y|)} e^{i \frac{\omega}{c^*} ((\frac{c}{c^*})^2 - 1) \frac{|\varphi_c^{-1}(z)^\perp|^2}{|\varphi_c^{-1}(z)|}} \\ G^2 \left( \frac{\omega \ell_0}{c^*} \left( \frac{\varphi_c^{-1}(z)^\perp}{|\varphi_c^{-1}(z)|} - \frac{y^\perp}{|y|} \right), \frac{\omega \ell_0^2}{c^*} \left( \left( \frac{c^*}{c} \right)^2 \frac{1}{|\varphi_c^{-1}(z)|} - \frac{1}{|y|} \right) \right) d\omega + \mathcal{O}(1).$$

We have, for $y \in \mathcal{D}_\delta(y_0)$, in light of remark 4.4 $(y - y_0)^{\parallel} = \mathcal{O}(\eta)$ and $|(y - y_0)^\perp| = \mathcal{O}(1)$. This stays true for $y \in \mathcal{D}_\delta(y_0 - \varphi_c^{-1}((\Delta z)_\eta))$. First, we write, for $z = \varphi_c(y_0) + (\Delta z)_\eta$:

$$|\varphi_c^{-1}(z)| - |y| = |y_0| - |y| + \eta \frac{c^*}{c} \Delta z^{\parallel} + \mathcal{O}(\eta^{\frac{3}{2}}).$$

Moreover

$$|y - \varphi_c^{-1}(\Delta z)| = |y| - \eta \frac{c^*}{c} \Delta z^{\parallel} + \mathcal{O}(\eta^{\frac{3}{2}}).$$

We can then write:

$$|\varphi_c^{-1}(z)| - |y| = |y_0| - |y - \varphi_c^{-1}(\Delta z)| + \mathcal{O}(\eta^{\frac{3}{2}}). \qquad (24)$$

Similarly:

$$\varphi_c^{-1}(z)^\perp = y_0^\perp + \mathcal{O}(\eta^{\frac{1}{2}}),$$

and

$$\frac{|\varphi_c^{-1}(z)^\perp|^2}{|\varphi_c^{-1}(z)|} = \frac{|y_0^\perp|^2}{|y_0|} + \mathcal{O}(\eta^{\frac{1}{2}}). \qquad (25)$$

Secondly,

$$\begin{aligned} \frac{\varphi_c^{-1}(z)^\perp}{|\varphi_c^{-1}(z)|} - \frac{y^\perp}{|y|} &= \frac{y_0^\perp + \varphi_c^{-1}(\Delta z)^\perp}{|y_0|} - \frac{y^\perp}{|y_0 + (y - y_0)|} + \mathcal{O}(\eta) \\ &= \frac{y_0^\perp}{|y_0|} - \frac{y^\perp - \varphi_c^{-1}(\Delta z)^\perp}{|y_0|} + \mathcal{O}(\eta) \end{aligned} \qquad (26) \\ = \frac{y_0^\perp}{|y_0|} - \frac{y^\perp - \varphi_c^{-1}(\Delta z)^\perp}{|y - \varphi_c^{-1}(\Delta z)|} + \mathcal{O}(\eta),$$

and finally:

$$\begin{aligned} \left(\frac{c^*}{c}\right)^2 \frac{1}{|\varphi_c^{-1}(z)|} - \frac{1}{|y|} &= \left(\frac{c^*}{c}\right)^2 \frac{1}{|y_0|} - \frac{1}{|y|} + \mathcal{O}(\eta) \\ &= \left(\frac{c^*}{c}\right)^2 \frac{1}{|y_0|} - \frac{1}{|y - \varphi_c^{-1}(\Delta z)|} + \mathcal{O}(\eta). \end{aligned} \qquad (27)$$

Plugging back (24), (25), (26), (27) in the expression of $F$ we recognize $F^c(\varphi_c(y_0), y - \varphi_c^{-1}((\Delta z)_\eta)$. $\square$

Relying on the previous result, we are now ready to state one of the key result of the paper. We prove that, asymptotically, the imaging function at a point close to $z(c)$ can be written as another realization of the imaging function at the point $z(c)$.
**Proposition 4.3.** Local stationarity. For $s_0 \in \mathbb{R}^{d-1}$, $t_0 > 0$, let $z(c) := (\eta^{\frac{1}{2}}c^2s_0, ct_0)$. Then for a.e. $\omega \in \Omega$ and $(\Delta z)_\eta := (\eta\Delta z^\perp, \eta\Delta z^{\parallel}) \in \mathbb{R}^d$,

$$I^c(z(c) + (\Delta z)_\eta, \omega) = I^c(z(c), \tau_{\varphi_c^{-1}((\Delta z)_\eta)}\omega) + \mathcal{O}(\eta^{1+\frac{d-2}{2}}).$$

*Proof.* The proof is based on a simple change of variable:

$$I^c(z(c) + (\Delta z)_\eta, \omega) = \int_{D_\delta^c(\varphi_c^{-1}(z(c) + (\Delta z)_\eta))} (n_\epsilon(y, \omega) - n^*) F^c(z + \Delta z, y) dy$$

We recall that since $z(c) = (\eta^{\frac{1}{2}}c^2s_0, ct_0)$ we can write $z(c) = \varphi_c(y_0)$ for some $y_0 \in D$. Since
$$z + (\Delta z)_\eta = \varphi_c(y_0) + (\Delta z)_\eta,$$

we know from Lemma 4.2

$$F^c(z + (\Delta z)_\eta, y) = F^c(z, y - \varphi_c^{-1}((\Delta z)_\eta)) + \mathcal{O}(\eta^{-\frac{1}{2}})$$

Introducing $\tilde{y} := y - \varphi_c^{-1}((\Delta z)_\eta)$ and using the stationarity of $n$, we get

$$I^c(z(c) + (\Delta z)_\eta, \omega) = \int_{D_\delta^c(y_0)} (n_\epsilon(\tilde{y}, \tau_{\varphi_c^{-1}((\Delta z)_\eta)}\omega) - n^*) F^c(z, \tilde{y}) d\tilde{y}.$$

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

We can now use proposition 4.3 to show that if one averages spatially the squared modulus of the imaging function over a small area of typical size $a$ of the order of the central wavelength, then we can approximate the statistical variance of the imaging functional.

**Proposition 4.4 (From ensemble to spatial averaging).** For $a \sim \eta \frac{c^*}{\omega_0}$,

$$\mathbb{E}\left[\left|\mathbb{E}\left[|I^c(z(c)|^2\right) - \int_{\square a} |I^c(z(c) + \Delta z)|^2 d\Delta z\right|^2\right] \lesssim \eta^{2\alpha d} \|\Phi^{\frac{1}{2}}(|\cdot|)\|_{L^1(\mathbb{R}^d)}^2 \left(\int_{D_\delta^c(y_0)} |F^c(z, y)|^2 dy\right)^2,$$

where $\square a := [-\frac{a}{2}, \frac{a}{2}]^d$.

*Proof.* We first write

$$\mathbb{E}\left[\left[\int_{\square a} |I^c(z(c) + \Delta z)|^2 d\Delta z\right]^2\right] = \int_{\square a^2} \mathbb{E}\left[|I^c(z(c) + \Delta z)|^2 |I^c(z(c) + \Delta z')|^2\right] d\Delta z d\Delta z'$$

We know from the previous proposition that

$$\mathbb{E} \left[ |I^c(z(c) + \Delta z)|^2 |I^c(z(c) + \Delta z')|^2 \right] = \int_{D_\delta^c(y_0)^4} \mathbb{E} \left[ (n_\epsilon(y - \varphi_c^{-1}(\Delta z)) - n^*) \right. \\
\left. (n_\epsilon(y' - \varphi_c^{-1}(\Delta z)) - n^*) (n_\epsilon(\tilde{y} - \varphi_c^{-1}(\Delta z')) - n^*) (n_\epsilon(\tilde{y}' - \varphi_c^{-1}(\Delta z')) - n^*) \right] F^c(z, y) F^c(z, y') F^c(z, \tilde{y}) \\
F^c(z, \tilde{y}') dy dy' d\tilde{y} d\tilde{y}' + \mathcal{O}(\eta^{\frac{3}{2}+2d}).$$

Using [22, Lemma IV-4] or [4, Lemma 2.1], we can bound the 4th-order moment of $n_\epsilon$ as follows

$$\left| \mathbb{E} \left[ (n_\epsilon(y - \varphi_c^{-1}(\Delta z)) - n^*)(n_\epsilon(y' - \varphi_c^{-1}(\Delta z)) - n^*) \right] \right| \lesssim \sum_{(y_1, y_2, y_3, y_4) \in \mathcal{U}} \Phi \left( \frac{|y_1 - y_2|}{\epsilon} \right)^{\frac{1}{2}} \Phi \left( \frac{|y_3 - y_4|}{\epsilon} \right)^{\frac{1}{2}}, \quad (28)$$
where $\mathcal{U}$ is the set of all permutations of $((y - \varphi_c^{-1}(\Delta z), y' - \varphi_c^{-1}(\Delta z), \tilde{y} - \varphi_c^{-1}(\Delta z'), \tilde{y}' - \varphi_c^{-1}(\Delta z')))$. After a change of variables, we obtain

$$ \mathbb{E} \left[ |I^c(z(c) + \Delta z)|^2 |I^c(z(c) + \Delta z')|^2 \right] \le \varepsilon^{2d} \|\Phi^{\frac{1}{2}}\|_{L^1(\mathbb{R}^+)}^2 \left( \int_{D_{\delta}^c(y_0)} |F^c(z, y)|^2 dy \right)^2 . $$

We used here that $\Phi^{\frac{1}{2}}(|\cdot|) \in L^1(\mathbb{R}^d)$ (see (6) and lemma A.1). Similarly we know from Lemma 4.1 that

$$ \mathbb{E} \left[ |I^c(z(c))|^2 \right] \sim \varepsilon^d \int_{\mathbb{R}^d} C(s) ds \int_{D_{\delta}^c(y_0)} |F^c(z, y)|^2 dy. $$

To conclude, recall that $\varepsilon = \eta^\alpha$ and we get the result.

□

## 4.5. Discussion.

### 4.5.1. Summary of the results.

(1) For a fixed time $t_0$, for any $c$ the pixel at position $z(c) = (0, ct_0)$ always probes a small area $D^c(y_0)$ in the domain around $y_0 = (0, c^*t_0)$ regardless of $c$.

$$ I^c(z(c)) \sim \int_{D^c(y_0)} (n_\varepsilon(y) - n^*) F^c(z(c), y) dy $$

Proposition 4.1 gives us an asymptotic expansion of $F^c(z(c), y)$ for $y$ in the small region. The shape and size of this region $D^c(y_0)$ depend on $\frac{c^*}{c}$ via the behavior of $G_\theta$ but the orders of magnitude of both transverse and axial size of the domain $D^c(y_0)$ remains the same over a range of $c$ (see remark 4.4). The takeaway is that $y \mapsto F^c(z(c), y)$ varies slowly over $D^c(y_0)$. Therefore one can see $I^c(z(c))$ as the convolution of a slowly varying (at the scale of the focal spot) function parametrized by $c$ with the realization of a random process over a domain centered around a fixed point but whose shape is distorted by variations of $c$.

(2) In practice (and in the theoretical model) the separation of scale is such that $y \mapsto n_\varepsilon(y) - n^*$ is the realization of a random process that oscillates spatially at a scale of order $\varepsilon$ much smaller than the size of the dual focal spot $D^c(y_0)$. Lemma 4.1 quantifies this qualitative observation and states that the variance of the random variable $I^c(z(c))$ is proportional to $\int_{D^c(y_0)} |F^c(z(c), y)|^2 dy$.

(3) We know from section 3.1 that $c \mapsto F^c(z(c), y_0)$ is a function strongly peaked at $c = c^*$ (see figure 3). Using the asymptotic expansion of proposition 4.1 which characterizes the size of the focal spot, we can do a Taylor expansion of $y \mapsto F^c(z(c), y)$ around $y_0$ that is uniformly valid with respect to $c$ in $D^c(y_0)$. We then deduce that $\int_{D^c(y_0)} |F^c(z(c), y)|^2 dy$ behaves like $|F^c(z(c), y_0)|^2$ and exhibits also a strong peak at $c = c^*$. Therefore the variance of $c \mapsto I^c(z(c))$ exhibits a maximum at $c = c^*$.

(4) In practice however, one has access to only one realization of the medium. This problem of accessing statistical moments of $c \mapsto I^c(z(c))$ can be worked around using the stationarity of the medium. Lemma 4.2 states that in the proper asymptotic regime, shifting the imaging point by a small $\Delta z$ (of the order of the wavelength) corresponds to shifting the probed area $D^c(y_0)$ by some small $\Delta y_0$ (of the order of the wavelength too). Since $\Delta y \gg \varepsilon r_0$, $n_\varepsilon(y) - n^*$ in $D^c(y_0)$ and $D^c(y_0) + \Delta y_0$ can be considered as independent realizations of the distribution of scatterers in $D^c(y_0)$ (proposition 4.3). The quantification of the error between the statistical and spatial variance of $c \mapsto I^c(z(c))$ is done in proposition 4.4.

### 4.5.2. Practical recovery of $c^*$. We show in this section how the estimator can be constructed in practice as done in [7] for a two dimensional medium. We assume that $c^* \in [c_{min}, c_{max}]$. We uniformly discretize the set of possible values for $c$ with $N_c$ points as $c_j = c_{min} + j\Delta c$.

(1) Fix a travel time $t_0 \in \mathbb{R}^+$. The transverse offset $s_0$ can be set to 0 without loss of generality. For every $j \in [1, N_c]$, consider the moving imaging point

$$ z_j = \psi_{c_j}(0, t_0) := (0, c_j t_0). $$

(2) Fix a set of $N_\Delta$ spatial shifts $(\Delta z)_i \in (\mathbb{R}^2)^{N_\Delta}$ satisfying $\varepsilon r_0 \ll |\Delta z| < \eta \frac{c^*}{\omega_0}$.
(3) Compute the matrix

$$K_{i,j} = I^{c_j}(z_j + \Delta z_i).$$

(4) Compute the variance of K with respect to its rows:

$$F_j = \frac{1}{N_{\Delta}} \sum_{i=1}^{N_{\Delta}} |K_{i,j}^2|.$$

(5) Compute $j^* = \operatorname{argmax}_j |F_j|$.

We proved that if $N_{\Delta}$ is large enough then we have

$$|c^* - c_{j^*}| \leq \Delta c.$$

**Remark 4.5.** Note that the computation of matrix M at the third step does not require computing the whole image. Only $N_{\delta} \times N_c$ pixels of the image are computed. In practice $N_{\Delta}$ is of order of a few dozen.

## 5. NUMERICAL ILLUSTRATIONS

The numerical simulations are performed in the time domain by using the K-Wave library [33]. For a given $x_e \in \mathbb{R}^2$, and a given end time $T > 0$, we simulate $U_e(x, t)$, the solution in $L^2([0, T], H^1(\mathbb{R}^d \setminus \{x_e\}))$ of:

$$\begin{cases} -\Delta U_e(x,t) + n(x) \partial_{tt}^2 U_e(x,t) = \delta(x - x_e)f(t), & \text{for } x \in \mathbb{R}^2 \setminus \{x_e\}, t \in [0, T], \\ U_e(x,0) = \partial_t U_e(x,0) = 0, & \end{cases} \qquad (29)$$

with $f(t) \in C^0([0, T])$ (see Figure 5). The signal $U_e(x_r, t)$ is then recorded on the transducers at $x_r$.

By considering $N_e = 15$ incident waves emitted at $\{x_e\}_{e=1...N_e}$ and $N_r = 64$ recording sensors placed at $\{x_r\}_{r=1...N_r}$, we have access to the matrix of data $M(\omega_k) \in \mathbb{C}^{N_e \times N_r}$ where

$$M(x_e, x_r, \omega_k) = \mathcal{F}(U_e(x_r, \cdot))(\omega_k)$$

with $\mathcal{F}$ the Fourier transform and $\{\omega_k\}_{k=1...N_\omega}$ are the sampling frequencies in $\mathcal{B} = [10^2; 10^8]$. Here $N_\omega = 1000$. The sensors $\{x_e\}_{e=1...N_e}$ and $\{x_r\}_{r=1...N_r}$ are equally spaced on the segment $[-\ell, \ell] \times \{0\}$ with $\ell = 1.5 \times 10^{-2}$ m.

As the problem is a multi-scale problem, the simulation of $M$ can hardly be done on a personal laptop. Indeed, it requires to mesh the small inclusions to capture their effects. With the CFL condition taken to be 0.4, the computations can become quite long. However, for $e = 1...N_e$, the simulations of $U_e$ can be done in parallel, and are then computed on Nvidia Tesla V100 GPUs. For the choice of parameters, it takes ~ 10 hours per simulation.

The parameters $N_e, N_r$ and $N_\omega$ are chosen such that the matrix $M$ can be processed on a personal laptop.

We now consider the medium of Figure 6 which is composed on 22710 scatterers. The distribution of the centers of the scatterers is a realization of a Matèrn point process (see [16, Section 6.5.2]). In Figure 6, there are 22710 scatterers with up to 30% contrast in the speed of sound and a typical radius $\epsilon r_0 = 7.5 \times 10^{-5}$ m. This corresponds to a particle volume fraction of 15%.

We follow the procedure of section 4.5.2 with $s_0 = 0$ and $t_0 = 3 \times 10^{-5}$ s. As $c^* = 1500 \text{ m.s}^{-1}$, this corresponds to $y_0 = (0, 45)$ mm. The matrix $K$ is shown on figure 7. The estimator is plotted on figure 8 and compared to the theoretical expression given in lemma 4.1.

## 6. PERSPECTIVES

We have provided a robust mathematical framework for the analysis of ultrasound imaging signals. We have proved that that it is possible to recover the speed of sound in situ in a homogeneous medium via the construction of sound speed estimators. Extending these results to more general media (micro-structured with a slowly varying background or a piecewise constant background, contrast in the divergence...) in order to reconstruct a full map of the speed of sound directly from the backscattering measurements is the subject of a forthcoming work.
<figure><img src="image_8.png" /><figcaption>FIGURE 5. Source term used in the simulation.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 6. Speed of sound map in the random multi-scale medium used in the simulations.</figcaption></figure>

## APPENDIX A. ASSUMPTIONS ON THE RANDOM SETTING

Let $(\Omega, \mathcal{F}, \mathbb{P})$ be a probability space. We define stationarity through an action $(\tau_x)_{x \in \mathbb{R}^d}$ of the group $(\mathbb{R}^d, +)$ on $(\Omega, \mathcal{F})$ that verifies:

- the map $\tau : \begin{cases} \mathbb{R}^d \times \Omega \\ (x, \omega) \end{cases} \rightarrow \begin{cases} \Omega \\ \tau_x \omega \end{cases}$ is measurable,

- $\forall x, y \in \mathbb{R}^d, \tau_{x+y} = \tau_x \circ \tau_y,$
<figure><img src="image_10.png" /><figcaption>FIGURE 7. Visualization of $|K_{ij}| = |I^{c_j}(z(c_j) + (\Delta z)_i, \omega_0)|$ for a fixed realization $\omega_0 \in \Omega$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 8. Plot of the estimator function $c_j \mapsto F_j = \frac{1}{N_\Delta} \sum_{i=1}^{N_\Delta} |K_{ij}^2|$.</figcaption></figure>
- For all $x \in \mathbb{R}^d$, $\tau_x$ preserves $\mathbb{P}$, i.e.

$$ \forall A \in \mathcal{F}, \mathbb{P}(\tau_x A) = \mathbb{P}(A). $$

We say that a process $f: \mathbb{R}^d \times \Omega \to \mathbb{R}$ is stationary if it verifies a.e. $\omega \in \Omega$,

$$ f(x + y, \omega) = f(x, \tau_y \omega), \quad \forall x, y \in \mathbb{R}^d. \qquad (30) $$

We suppose moreover that the action $(\tau_x)_{x \in \mathbb{R}^d}$ is ergodic.

**Lemma A.1 (Consequences of (6)).** (a) $\Phi(|\cdot|) \in L^1(\mathbb{R}^d)$ and $\Phi^{\frac{1}{2}}(|\cdot|) \in L^1(\mathbb{R}^d)$.

(b) Moreover, for $R$ large enough

$$ \left| \int_{\mathbb{R}^d} C(x) dx - \int_{B(0,R)} C(x) dx \right| \lesssim R^{d-1-p}. $$

In particular, for a given $\varepsilon > 0$ one can choose $R = \varepsilon^{\frac{1}{d-1-p}}$ and get

$$ \left| \int_{\mathbb{R}^d} C(x) dx - \int_{B(0,R)} C(x) dx \right| \lesssim \varepsilon. $$

**Remark A.1.** In practice the decay rate of $\Phi$ -see eq. (6)- encodes the rate of decorrelation of $n$. The idea is that $\Phi$ decays at a scale one order of magnitude smaller than the wavelength, i.e to have $\varepsilon R = o(\eta \frac{c^*}{\omega_0})$. The condition $p > 2(d-1)$ is enough for $\Phi(|\cdot|) \in L^1(\mathbb{R}^d)$ and $\Phi^{\frac{1}{2}}(|\cdot|) \in L^1(\mathbb{R}^d)$. Moreover this condition is also compatible with the asymptotic regime described in section 2.5 in the sense that the condition on $\alpha$ derived in remark 2.7 is more restrictive.

# APPENDIX B. HOMOGENIZATION IN LARGE DOMAINS

**Proposition B.1.** Let $u$ be the solution of (2) in $H_{\text{loc}}^2(\mathbb{R}^d)$ and $u^i$ the incident wave. Let $\mathcal{U}^s$ be defined for all $x_e, x_r \in \mathcal{P}$, $\omega \in \mathcal{B}$ by:

$$ \mathcal{U}^s(x_e, x_r, \omega) := \omega^2 \int_D (n_\varepsilon(x) - n^*) \Gamma^{\frac{\omega}{c^*}}(x, x_e) \Gamma^{\frac{\omega}{c^*}}(x, x_r) dx. $$

It holds:

$$ \|(u - u^i - \mathcal{U}^s)(x_e, x_r, \omega)\|_{L^2(\Omega)} \lesssim C_\omega \varepsilon^{\frac{d}{p'}} \qquad (31) $$

for some $C_\omega > 0$ independent on $x_e, x_r$ and $\varepsilon$ and bounded by $\omega^6 \|u^i\|_{H^1(D)} + \omega^4$ and $p' := \frac{2+\beta}{1+\beta}$ with $\beta$ chosen such that

$$ \sup_{x \in D} \int_D \left( \Gamma^{\frac{\omega}{c^*}}(x, y) \right)^{2+\beta} dy < \infty. $$

*Proof.* The starting point is the integral representation formula for the scattered field:

$$ u(y) - u^i(y) = \omega^2 \int_D (n_\varepsilon(x) - n^*) \Gamma^{\frac{\omega}{c^*}}(x, y) u(x) dx, $$

for $y \in \mathbb{R}^d \setminus \{x_e\}$. Replacing twice $u$ by its integral representation inside the integral leads to:

$$ \begin{align} u(x_r) - u^i(x_r) &= \int_D \omega^2 (n_\varepsilon(x) - n^*) \Gamma^{\frac{\omega}{c^*}}(x, x_r) u^i(x) dx \nonumber \\ &\quad + \omega^4 \int_{D^2} (n_\varepsilon(x_1) - n^*) (n_\varepsilon(x_2) - n^*) \Gamma^{\frac{\omega}{c^*}}(x_1, x_2) \Gamma^{\frac{\omega}{c^*}}(x_1, x_r) u^i(x_2) dx_1 dx_2 \nonumber \\ &\quad + \omega^6 \int_{D^3} (n_\varepsilon(x_1) - n^*) (n_\varepsilon(x_2) - n^*) (n_\varepsilon(x_3) - n^*) \Gamma^{\frac{\omega}{c^*}}(x_3, x_2) \Gamma^{\frac{\omega}{c^*}}(x_2, x_1) \Gamma^{\frac{\omega}{c^*}}(x_1, x_r) u(x_3) dx_1 dx_2 dx_3 \tag{32} \\ &:= \mathcal{U}^s(x_r) + \mathcal{U}^{s,2}(x_r) + \mathcal{U}^{s,3}(x_r). \nonumber \end{align} $$
To estimate $\mathcal{U}^{s,2}$ and $\mathcal{U}^{s,3}$ we use [4, Lemma 2.1] with $q = n - n^*$ (recall that $n_\epsilon = n(\frac{\cdot}{\epsilon})$). The hypothesis (6) from [4] is verified because the covariance C of $n_\epsilon$ is bounded by $\Phi$ in $L^1(\mathbb{R}^+) \cap L^{\frac{1}{2}}(\mathbb{R}^+)$. Therefore we have:

$$|\mathbb{E}[q(x_1)q(x_2)q(x_3)q(x_4)]| \lesssim \sup_{(y_1,y_2,y_3,y_4) \in \Xi} (\Phi(|y_1 - y_2|))^{\frac{2d}{2}} (\Phi(|y_3 - y_4|))^{\frac{1}{2}}$$

where $\Xi$ is the set of permutations of $(x_1, x_2, x_3, x_4)$. We have, following the same steps as in the proof of [4, Lemma 2.6]

$$\mathbb{E}[|\mathcal{U}^{s,2}(x_r)|^2] \lesssim \omega^8 \epsilon^{\frac{2d}{p'}} \|\Gamma_{\frac{\omega}{c^*}}^{\frac{\omega}{c^*}}(\cdot, x_r)\|_{L^\infty(D)}^2 \|u^i\|_{L^\infty(D)}^2 \|\Phi^{\frac{1}{2}}\|_{L^{p'}}^2 \|\Gamma_{\frac{\omega}{c^*}}^{\frac{\omega}{c^*}}(\cdot, \cdot)\|_{L^p(D^2)}$$

with $p = 2 + \beta$ and $p' = \frac{2+\beta}{1+\beta}$ with $\beta$ chosen such that

$$\sup_{x \in D} \int_D \left( \Gamma_{\frac{\omega}{c^*}}^{\frac{\omega}{c^*}}(x, y) \right)^{2+\beta} dy < \infty.$$

The term $\mathcal{U}^{s,3}$ can be dealt with in a similar manner. The difference is the presence of the term $u$ instead of $u^i$ in the integrand. To deal with this term we first recall that a.s. :

$$\|u\|_{H^1(D)} \lesssim \|u^i\|_{H^1(D)}.$$

Then we can perform a Cauchy-Schwarz type estimate to separate the term with $u$.

$$\mathbb{E}[|\mathcal{U}^{s,3}(x_r)|^2] \lesssim \mathbb{E} \left[ \int_D \left| \int_{D^2} (n_\epsilon(x_1) - n^*)(n_\epsilon(x_2) - n^*) \Gamma_{\frac{\omega}{c^*}}^{\frac{\omega}{c^*}}(x_1, x_2) \Gamma_{\frac{\omega}{c^*}}^{\frac{\omega}{c^*}}(x_1, x_r) \right. \right. \\ \left. \left. \Gamma_{\frac{\omega}{c^*}}^{\frac{\omega}{c^*}}(x_3, x_2) dx_1 dx_2 \right|^2 dx_3 \int_D (n_\epsilon(x_1) - n^*)^2 u^2(x_1) dx_1 \right].$$

In the end we obtain:

$$\mathbb{E}[|\mathcal{U}^{s,3}(x_r)|^2] \lesssim \omega^{12} \epsilon^{\frac{2d}{p'}} \|\Gamma_{\frac{\omega}{c^*}}^{\frac{\omega}{c^*}}(\cdot, x_r)\|_{L^\infty(D)}^2 \|u^i\|_{H^1(D)}^2 \|\Phi^{\frac{1}{2}}\|_{L^{p'}}^2 \|\Gamma_{\frac{\omega}{c^*}}^{\frac{\omega}{c^*}}(\cdot, \cdot)\|_{L^p(D^2)} \\ \sup_{x \in D} \int_D \left( \Gamma_{\frac{\omega}{c^*}}^{\frac{\omega}{c^*}}(x, y) \right)^2 dy. \quad \square$$

## REFERENCES

[1] Giovanni S. Alberti, Habib Ammari, Francisco Romero, and Timothée Wintz. Mathematical analysis of ultrafast ultrasound imaging. *SIAM Journal on Applied Mathematics*, 77(1):1–25, 2017.

[2] Habib Ammari, Hyeonbae Kang, and Hyundae Lee. *Layer potential techniques in spectral analysis*. Number 153. American Mathematical Soc., 2009.

[3] Alexandre Aubry, William Lambert, Laura Cobus, and Mathias Fink. Reflection matrix approach for quantitative ultrasound imaging of scattering media. In *Forum Acusticum*, pages 1791–1791, 2020.

[4] Guillaume Bal. Central limits and homogenization in random media. *Multiscale Modeling & Simulation*, 7(2):677–702, 2008.

[5] Samuel Beuret, Baptiste Hériard-Dubreuil, Naiara Korta Martiartu, Michael Jaeger, and Jean-Philippe Thiran. Windowed radon transform for robust speed-of-sound imaging with pulse-echo ultrasound. *IEEE transactions on medical imaging*, 43(4):1579–1593, 2023.

[6] Flavien Bureau. *Multi-dimensional analysis of the reflection matrix for quantitative ultrasound imaging*. Theses, Université Paris sciences et lettres, July 2023.

[7] Flavien Bureau, Elsa Giraudat, Arthur Le Ber, William Lambert, Louis Carmier, Aymeric Guibal, Mathias Fink, and Alexandre Aubry. Reflection matrix imaging for wave velocity tomography. *arXiv preprint arXiv:2409.13901*, 2024.

[8] Marco Dioguardi Burgio, Marion Imbault, Maxime Ronot, Alex Faccinetto, Bernard E Van Beers, Pierre-Emmanuel Rautou, Laurent Castera, Jean-Luc Gennisson, Mickael Tanter, and Valérie Vilgrain. Ultrasonic adaptive sound speed estimation for the diagnosis and quantification of hepatic steatosis: a pilot study. *Ultraschall in der Medizin-European Journal of Ultrasound*, 40(06):722–733, 2019.

[9] Guy Cloutier, Francois Destrempes, Fei Yu, and An Tang. Quantitative ultrasound imaging of soft biological tissues: a primer for radiologists and medical physicists. *Insights Imaging*, 12(1):127, Sep 2021.

[10] NHS England and N Improvement. Diagnostic imaging dataset annual statistical release. *London: Department of Health*, 2020.

[11] Florian Faucher and Otmar Scherzer. Adjoint-state method for hybridizable discontinuous galerkin discretization, application to the inverse acoustic wave problem. *Computer Methods in Applied Mechanics and Engineering*, 372:113406, 2020.
[12] Jean-Pierre Fouque, Josselin Garnier, George Papanicolaou, and Knut Solna. *Wave propagation and time reversal in randomly layered media*, volume 56. Springer Science & Business Media, 2007.

[13] Josselin Garnier, Laure Giovangigli, Quentin Goepfert, and Pierre Millien. Scattered wavefield in the stochastic homogenization regime, 2023.

[14] Josselin Garnier and George Papanicolaou. *Passive imaging with ambient noise*. Cambridge University Press, 2016.

[15] Martin Harwit. *Hadamard transform optics*. Elsevier, 2012.

[16] Janine Illian, Antti Penttinen, Helga Stoyan, and Dietrich Stoyan. *Statistical analysis and modelling of spatial point patterns*. John Wiley & Sons, 2008.

[17] Marion Imbault, Marco Dioguardi Burgio, Alex Faccinetto, Maxime Ronot, Hanna Bendjador, Thomas Deffieux, Emma Ollivier Triquet, Pierre-Emmanuel Rautou, Laurent Castera, Jean-Luc Gennisson, et al. Ultrasonic fat fraction quantification using in vivo adaptive sound speed estimation. *Physics in Medicine & Biology*, 63(21):215013, 2018.

[18] Marion Imbault, Alex Faccinetto, Bruno-Félix Osmanski, Antoine Tissier, Thomas Deffieux, Jean-Luc Gennisson, Valérie Vilgrain, and Mickaël Tanter. Robust sound speed estimation for ultrasound-based hepatic steatosis assessment. *Physics in Medicine & Biology*, 62(9):3582, 2017.

[19] Michael Jaeger, Georg Held, Sarah Peeters, Stefan Preisser, Michael Grünig, and Martin Frenz. Computed ultrasound tomography in echo mode for imaging speed of sound using pulse-echo sonography: proof of principle. *Ultrasound in Medicine & Biology*, 41(1):235–250, Jan 2015.

[20] Michael Jaeger, Patrick Stähli, Naiara Korta Martiartu, Parisa Salemi Yolgunlu, Thomas Frappart, Christophe Fraschini, and Martin Frenz. Pulse-echo speed-of-sound imaging using convex probes. *Physics in Medicine & Biology*, 67(21):215016, 2022.

[21] P. Joly. *Introduction à l'analyse mathématique de la propagation d'ondes en régime harmonique*. 2006.

[22] Harold Joseph Kushner. *Approximation and weak convergence methods for random processes, with applications to stochastic systems theory*, volume 6. MIT press, 1984.

[23] William Lambert, Laura A Cobus, Justine Robin, Mathias Fink, and Alexandre Aubry. Ultrasound matrix imaging part ii: The distortion matrix for aberration correction over multiple isoplanatic patches. *IEEE Transactions on Medical Imaging*, 41(12):3921–3938, 2022.

[24] William Lambert, Justine Robin, Laura Cobus, Mathias Fink, and Alexandre Aubry. Ultrasound matrix imaging part i: The focused reflection matrix, the f-factor and the role of multiple scattering. *IEEE Transactions on Medical Imaging*, 41(12):3907–3920, 2022.

[25] Andrea Moiola and Euan A Spence. Acoustic transmission problems: wavenumber-explicit bounds and resonance-free regions. *Mathematical Models and Methods in Applied Sciences*, 29(02):317–354, 2019.

[26] Vincent Perrot, Maxime Polichetti, François Varray, and Damien Garcia. So you think you can das? a viewpoint on delay-and-sum beamforming. *Ultrasonics*, 111:106309, 2021.

[27] Alan D Pierce. *Basic linear acoustics*. Springer, 2007.

[28] Ali Mohammad Pirmoazen, Anjali Khurana, Ahmed El Kaffas, and Aya Kamaya. Quantitative ultrasound approaches for diagnosis and monitoring hepatic steatosis in nonalcoholic fatty liver disease. *Theranostics*, 10(9):4277–4289, Mar 2020.

[29] Claire Prada, Sébastien Manneville, Dimitri Spoliansky, and Mathias Fink. Decomposition of the time reversal operator: Detection and selective focusing on two scatterers. *The Journal of the Acoustical Society of America*, 99(4):2067–2076, 1996.

[30] SV Ranganayakulu, N Ramakoteswara Rao, and Leena Gahane. Ultrasound applications in medical sciences. *IJMTER*, 3(02):287–293, 2016.

[31] Plamen Stefanov, Gunther Uhlmann, Andras Vasy, and Hanming Zhou. Travel time tomography. *Acta Mathematica Sinica, English Series*, 35(6):1085–1114, 2019.

[32] Mickael Tanter and Mathias Fink. Ultrafast imaging in biomedical ultrasound. *IEEE transactions on ultrasonics, ferroelectrics, and frequency control*, 61(1):102–119, 2014.

[33] Bradley E Treeby and Benjamin T Cox. k-wave: Matlab toolbox for the simulation and reconstruction of photoacoustic wave fields. *Journal of biomedical optics*, 15(2):021314–021314, 2010.

[34] Tadashi Yamaguchi. Basic concept and clinical applications of quantitative ultrasound (qus) technologies. *Journal of Medical Ultrasonics*, 48(4):391–402, 2021.