# Self-Portrait of the Focusing Process in Speckle: I. Spatio-Temporal Imaging of Wave Packets in Complex Media

Elsa Giraudat, Flavien Bureau, William Lambert, Mathias Fink, Alexandre Aubry

► To cite this version:

Elsa Giraudat, Flavien Bureau, William Lambert, Mathias Fink, Alexandre Aubry. Self-Portrait of the Focusing Process in Speckle: I. Spatio-Temporal Imaging of Wave Packets in Complex Media. 2026. hal-05497423

HAL Id: hal-05497423

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

Preprint submitted on 6 Feb 2026

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 CC BY 4.0 - Attribution - International License
# Self-Portrait of the Focusing Process in Speckle:

## I. Spatio-Temporal Imaging of Wave Packets in Complex Media

Elsa Giraudat,¹ Flavien Bureau,¹ William Lambert,² Mathias Fink,¹ and Alexandre Aubry*¹

¹Institut Langevin, ESPCI Paris, PSL University, CNRS, 75005 Paris, France

²SuperSonic Imagine, Aix-en-Provence, France

(Dated: February 6, 2026)
Abstract

This is the first article in a series of three dealing with the exploitation of speckle for imaging purposes. Speckle is the complex interference wave-field produced by a random distribution of un-resolved scatterers. In this paper, we show how these scatterers can be used as virtual microphones to monitor the spatio-temporal propagation of a wave-packet inside the medium. To do so, the concept of matrix imaging is particularly useful. It consists in decoupling the location of the transmitted and received focal spots in a standard beamforming process. By scanning the wave-field with the output focal spot that then acts as a virtual transducer, one can image the spatio-temporal evolution of the wave-packet inside the medium. This unique observable will allow us to highlight the imperfections of the focusing process, in particular the defocus and reverberations induced by a strong aberrating layer. As a proof-of-concept, we will consider ultrasound experiments on tissue-mimicking phantoms. In the next two papers, we will show how this observable can be leveraged to compensate for these phenomena that hamper wave focusing and imaging in all fields of wave physics. Our method is indeed broadly applicable to different types of waves beyond ultrasound for which multi-element technology allows a reflection matrix to be measured.

A speckle pattern is the result of a random interference between many waves. In optics, it is commonly observed at the output of a scattering medium illuminated by coherent light¹. In medical ultrasound, the speckle pattern arises in the image itself². It is the result of a random interference between echoes induced by a random distribution of un-resolved scatterers. A speckle pattern is characterized by its local short-range correlations both in the spatial and time domain³. However, strong correlations can exist between the incident wave-field and transmitted/reflected speckle patterns⁴⁻⁷. This memory effect has been exploited both in transmission, for imaging through strongly scattering media⁸,⁹ or, in reflection, for adaptive focusing¹⁰,¹¹. In this paper, we will show how speckle can be leveraged to image the propagation of the incident wave-field inside the medium.

To that aim, we will go beyond confocal imaging which consists in focusing onto a given point inside the medium and only collect echoes coming from the same point. This process can be done physically via a pinhole placed before the detector in optical microscopy or
numerically by applying time delays to the received signals in ultrasound. Here, we will decouple the focusing process at input and output in order to scan the incident wave-field with the output focal spot. Assuming speckle is everywhere inside the medium, we will be able to make a time movie of the incident wave inside the medium by exploiting the temporal resolution of our measurement.

To do so, we will exploit the concept of matrix imaging<sup>12</sup>. Experimentally, it relies on the measurement of the reflection matrix that contains the impulse response between a set of sensors that are placed on the same side of the medium. Once this matrix is measured, the response of the medium is known for any incident wave-field. Due to the linearity of the wave equation, one can apply an independent beamforming process at input and output in order to investigate the reflection matrix in the focused basis<sup>12,13</sup>. Each line of this focused reflection matrix contains the incident wave-field intended to focus at r<sub>in</sub> recorded by a virtual transducer at r<sub>out</sub>. This virtual transducer is synthesized by the focused beamforming process applied at output. This matrix has been, so far, only investigated at the ballistic time to probe the lateral extent of the incident wave-field and evaluate the focusing quality at this specific time<sup>12,14</sup>. Yet, wave velocity fluctuations can imply a mismatch between the focusing plane dictated by the spatial focusing process and the coherence or isochronous volume controlled by the echo time. Moreover, in presence of reverberations, each scatterer of the medium will not only generate one single echo but several echoes related to each multiply scattered path that can follow the incident and reflected wave-fronts during their travel from the medium surface and the focusing point. The goal of this study is therefore to exploit the time-dependence of the focused reflection matrix to probe these phenomena.

Two ultrasound data sets will be used to support our demonstration. The first one consists in the reflection matrix associated with an ultrasonic probe place in contact of a tissue-mimicking phantom whose background sound velocity is known. It will be used as a reference experiment in which the level of aberrations can be tuned by changing the wave velocity model in the beamforming process. The second data set corresponds to the same configuration except that a Plexiglas plate is placed between the probe and the tissue-mimicking phantom in order to mimic a reverberating layer. On the one hand, the bright reflectors contained by the phantom will be used to assess the impact of aberrations and reverberations on the ultrasound image. On the other hand, we will show how ultrasound speckle can be exploited to make a self-portrait of the focusing process, seen through the
prism of the wave velocity model. To do so, a coherent virtual reflector will be synthesized
from different realizations of input focal spots. In practice, this will be done by means of an
an iterative phase reversal algorithm that will be applied to the focused reflection matrix
expressed in a de-scanned basis. However, contrary to a previous study¹⁵, this process is
here applied not only at a central frequency and at the expected ballistic time, but over the
whole frequency bandwidth. The result is a movie of the coherent wave-packet that survives
to disorder and that highlights: (i) the time shift of the coherence plane and defocus of
the wave with respect to the expected ballistic time and focal plane, respectively, when a
mismatch exists between the speed-of-sound distribution inside the propagation medium and
the wave velocity model used by the beamforming process; (ii) the temporal dispersion of
the wave-packets induced by multiple reflections induced by a reverberating layer placed
between the probe and the focusing plane. At last, we will discuss how this unique observable
can be used to design a path towards: (i) a compensation of axial aberrations (including
defocus and reverberations) in the speckle regime and; (ii) a quantitative measurement of the
speed-of-sound in reflection. These perspectives will be the object of the two other papers of
the series¹⁶,¹⁷.

# I. TIME FOCUSED REFLECTION MATRIX

## A. Reflection matrix acquisition

Ultrasound matrix imaging (UMI) begins with the acquisition of the reflection matrix using
an ultrasound array of transducers (Figs. 1a,b and Tab. I). The sample under study is a tissue
mimicking phantom with a speed-of-sound $c = 1542 \pm 10$ m/s. It is composed of a random
distribution of unresolved scatterers which generate an ultrasonic speckle characteristic of
human tissue (Fig. 1c). It also contains a set of nylon rods. Those bright scatterers can be
used to evaluate locally the imaging point spread function (PSF) of the ultrasound device.
The reflection matrix is captured by sending a series of plane waves into the medium. Each
plane wave is identified with its angle of incidence $\theta_{\text{in}}$ (Fig. 1a). For each illumination, the
reflected waves are recorded by the transducers of the probe, each element being identified
by its lateral position $u_{\text{out}}$ (Fig. 1b). The recorded wave-fronts are noted $R(u_{\text{out}}, \theta_{\text{in}}, t)$, with
t the echo time. They are stored in a reflection matrix $\mathbf{R}_{u\theta}(t) = [R(u_{\text{out}}, \theta_{\text{in}}, t)]$.
<figure><img src="image_2.png" /><figcaption>FIG. 1. Experimental configurations and procedure. a, Experimental set up: An ultrasonic linear probe is placed in contact of a tissue-mimicking phantom composed of a random distribution of sub-resolved diffusers schematized by the gray background texture, two hyper-echoic cylinders, and a set of bright nylon rods. A set of plane waves ($\theta_{in}$) is emitted by the probe and the time-dependent back-scattered wave-field $R(u_{out}, \theta_{in}, t)$ recorded by each transducer ($u_{out}$) is stored in the so-called reflection matrix $R_{u\theta}(t) = [R(u_{out}, \theta_{in}, t)]$ sketch in panel b. c, d Ultrasound image deduced from a beamforming process applied to $R_{u\theta}(t)$ using $c_0 = 1540$ and $1440$ m.s$^{-1}$, respectively. e, Same experimental configuration as in panel a except that a reverberating layer is placed between the probe and the tissue-mimicking phantom. f, Corresponding ultrasound image ($c_0 = 1540$ m.s$^{-1}$).</figcaption></figure>

## B. Confocal imaging

The first post-processing step is to build a confocal image $\mathcal{I}$ from the recorded reflection matrix. To do so, a delay-and-sum beamforming process is applied to the coefficients of $R_{u\theta}(t)$. Physically equivalent to a confocal focusing process (Fig. 2a), this procedure writes mathematically as follows:

$$ \mathcal{I}(x, z_0 = c_0 t/2) = \sum_{\theta_{in}} \sum_{u_{out}} A(u_{out}, \theta_{in}, x, t, c_0) \qquad (1) \\ R(u_{out}, \theta_{in}, \tau_{out}(u_{out}, x, t, c_0) + \tau_{in}(\theta_{in}, x, t, c_0)). $$

$c_0$ is the wave velocity model considered in the beamforming process. $\tau_{in}$ is the time-of-flight expected for the incident plane wave to reach the target point of coordinates $(x, z_0)$. $\tau_{out}$ is the time-of-flight expected for the reflected wave to travel from the same target point to each transducer. A is a normalization and apodization factor that limits the extent of the receive synthetic aperture. $z_0 = c_0 t/2$ is the expected position of the isochronous volume, which is defined as the ensemble of points that contribute to the ultrasound signal
<table><caption>TABLE I. Probe and acquisition sequence used for recording the canonical reflection matrix $R_{u\theta}(t)$ in the experiments.</caption><thead><tr><th></th><th>parameter</th><th>value</th></tr></thead><tbody><tr><td rowspan="7">probe</td><td>type</td><td>linear</td></tr><tr><td>number of transducers</td><td>N<sub>u</sub><sup>(Tx)</sup> = 256</td></tr><tr><td>transmission</td><td>N<sub>u</sub><sup>(Rx)</sup> = 128</td></tr><tr><td>Reception</td><td>δu = 0.2 mm ~ λ</td></tr><tr><td>inter-element distance</td><td>Δu<sup>(Tx)</sup> = 51.2 mm</td></tr><tr><td>aperture</td><td>Δu<sup>(Rx)</sup> = 25.6 mm</td></tr><tr><td>central frequency</td><td>f<sub>c</sub> = 7.5 MHz</td></tr><tr><td>bandwidth</td><td>Δf = [2 - 10] MHz</td></tr><tr><td rowspan="10">acquisition</td><td>imaging platform</td><td>Aixplorer®, Supersonic Imagine</td></tr><tr><td>speed-of-sound model</td><td>c<sub>0</sub> = 1540 m s<sup>-1</sup></td></tr><tr><td>plane waves</td><td>θ<sub>max</sub> = 40°</td></tr><tr><td>maximum angle</td><td>δθ<sub>in</sub> = 1°</td></tr><tr><td>angular sampling</td><td>N<sub>in</sub> = 81</td></tr><tr><td>number</td><td>two half-period pulse at f<sub>c</sub></td></tr><tr><td>emitted signal</td><td>f<sub>s</sub> = 30 MHz</td></tr><tr><td>sampling frequency</td><td>Δt = 137 µs</td></tr><tr><td>recording time</td><td></td></tr></tbody></table>

at time $t$. When the wave velocity model $c_0$ matches with the medium speed-of-sound, the position $z_t$ of the isochronous volume matches with the focusing plane $z^*$. The ultrasound image is then a satisfactory image of the medium reflectivity (Fig. 1c). In presence of an aberrating layer associated with a speed-of-sound $c$ differing from $c_0$, the beamformed image is drastically affected by the mismatch between the isochronous volume and the focusing plane (Fig. 2b). This non-coincidence results in axial aberrations on the ultrasound image (Fig. 1d) that manifest as: (i) an axial shift of the scatterers with respect to their true depth; (ii) a degradation of the transverse resolution. In presence of reverberations (Fig. 1e), each multiply-reflected wave is associated with a distinct focusing depth and a shifted isochronous volume (Fig. 2c). Those reverberations induce multiple ghost images of each bright scatterer on the ultrasound image (Fig. 1f), thereby highlighting the spatio-temporal spreading of
<figure><img src="image_3.png" /><figcaption>FIG. 2. **Principle of ultrasound imaging.** a, Case of a medium of homogeneous speed-of-sound. Left: The excitation of each transducer by a time-delayed pulse generates a cylindrical wave-front focusing at a given point $\mathbf{r}_{\text{in}} = (x_{\text{in}}, z_{\text{in}})$. Right: In reception mode, the echoes back-scattered by reflectors in the vicinity of the same point are time-shifted using the same delay law as that applied in transmit. This allows the signals induced by a single scattering event at $\mathbf{r}_{\text{out}} = \mathbf{r}_{\text{in}}$ to be constructively summed. The pixel of the ultrasound image corresponds to the resulting signal measured at the expected ballistic time $t = 2z_{\text{in}}/c_0$. The set of points generating the backscattered echoes at time $t$ defines the isochronous volume. When the wave velocity model and the medium speed-of-sound coincide, the focusing plane of the incident wave-front belongs to the isochronous volume. b, For a medium with an heterogeneous speed-of-sound distribution, this coincidence is no longer checked: The focusing plane and the isochronous volume are axially shifted in reverse directions. The delay-and-sum beamforming process that assumes an homogeneous wave velocity $c_0$ is hampered by a defocus. c, For a reverberating medium, the situation is even worse since each multiply-reflected wave-front is associated with a different isochronous volume and a distinct focusing plane. Reverberation phenomena undergone by the wave at input and output gives rise to a series of echoes for each scatterer on the beamformed image.</figcaption></figure>

the imaging PSF. However, speckle statistics does not seem affected by aberrations and reverberations in Figs. 1d and f. Yet, speckle also carries the spatio-temporal distortions undergone by the incident and reflected wave-fields upstream. In the following, we will show how to reveal this hidden information with the help of matrix imaging.
Time Focused Reflection Matrix

Ultrasound matrix imaging can show the degradation of the focusing process by decoupling the input and output focusing points, $\mathbf{r}_{\text{in}}$ and $\mathbf{r}_{\text{out}}$. This process is achieved by projecting the reflection matrix in a focused basis. This operation can be performed through a simple matrix product in the frequency domain¹², or in the time domain¹⁵. In the latter option, it consists in a delay-and-sum beamforming process generalized to the case of different focal points in transmission and reception modes:

$$R(\mathbf{r}_{\text{out}}, \mathbf{r}_{\text{in}}, \tau) = \sum_{\theta_{\text{in}}} \sum_{u_{\text{out}}} A(u_{\text{out}}, \mathbf{r}_{\text{out}}, \theta_{\text{in}}, \mathbf{r}_{\text{in}}) R(u_{\text{out}}, \theta_{\text{in}}, \tau_{\text{out}}(u_{\text{out}}, x, t, c_0) + \tau_{\text{in}}(\theta_{\text{in}}, x, t, c_0) + \tau). \quad (2)$$

Contrary to previous works, the resulting focused reflection matrix, $\mathbf{R}_{\text{rr}}(\tau) = [R(\mathbf{r}_{\text{out}}, \mathbf{r}_{\text{in}}, \tau)]$, is here expressed as a function of a time lapse $\tau$ with respect to the expected ballistic time. An example of focused reflection matrix $\mathbf{R}_{\text{rr}}(\tau)$ is displayed at different lapse time in Fig. 3b. Given the 2D feature of points $\mathbf{r}_{\text{in}}$ and $\mathbf{r}_{\text{out}}$, the full reflection matrix is in five space-time dimensions. In previous works¹²,¹⁴, the focused reflection matrix was only considered at the ballistic time $\tau = 0$, which limited us to a study of transverse aberrations induced by the variations of the speed-of-sound with respect to its mean value. Yet, as illustrated by Fig. 1, axial aberrations can be far from negligible especially in presence of reverberations. In the present paper, this supplementary degree of freedom will be exploited to probe the time-dependence of the focusing process.

The focused reflection matrix can be physically interpreted as the inter-element response matrix associated with virtual transducers located in the medium¹². At emission, the digital focusing of energy at a point $\mathbf{r}_{\text{in}}$ can be interpreted as the synthesis of a virtual source at that point. At reception, the focusing point $\mathbf{r}_{\text{in}}$ can similarly be seen as a virtual detector, formed by digitally selecting the back-scattered echoes originating from that point. It should be noted that the concept of virtual transducers is purely didactic. Of course, these virtual transducers do not act as true active sources of energy. Furthermore, they are highly directional: downward for virtual sources, and upward for virtual receivers¹⁴.

The wave induced by a virtual source at $\mathbf{r}_{\text{in}}$ can be measured by scanning the virtual receiver $\mathbf{r}_{\text{out}}$ around it. To visualize the propagation of this wave across the medium, the focused reflection matrix can be expressed in a de-scanned frame using the following change
<figure><img src="image_4.png" /><figcaption>FIG. 3. Time-focused reflection matrix. a The matrix $R_{rr}(\tau)$ is constructed in post-processing by decoupling the position of the emission focusing points $r_{in} = (x_{in}, z_{in})$ and reception focal points $r_{out} = (x_{out}, z_{out})$, thus forming a set of virtual sources and receivers inside the medium. The temporal responses between each virtual source $r_{in}$ and each virtual receiver $r_{out}$ are obtained by introducing an additional variable time $\tau$, corresponding to a shift from the expected ballistic time $t = (z_{in} + z_{out})/c_0$. This time delay $\tau$ allows the wave to propagate virtually in the medium from $r_{in}$ to $r_{out}$. b Time focused reflection matrix $R_{rr}(\tau)$ displayed at different time lapses $\tau$ in the ultrasound phantom (Fig. 1a). c Same matrix $R_{\Delta r}(\tau)$ expressed in the de-scanned basis. d Example of a reflected wave-fields reshaped in de-scanned coordinates $(\Delta x, \Delta z)$ extracted from one column of the de-scanned matrix $R_{\Delta r}(\tau)$ designated by a black arrow in panel c and associated with a given point $r_{in}$ corresponding to a bright point-like scatterer. e Same as in panel d but in speckle.</figcaption></figure>

of variables:

$$R(\Delta \mathbf{r}, \mathbf{r}_{\text{in}}, \tau) = R(\mathbf{r}_{\text{out}}, \mathbf{r}_{\text{in}}, \tau) \qquad (3)$$

with $\Delta \mathbf{r} = \mathbf{r}_{\text{out}} - \mathbf{r}_{\text{in}}$. The resulting de-scanned matrix, $R_{\Delta \mathbf{r}}(\tau) = [R(\Delta \mathbf{r}, \mathbf{r}_{\text{in}}, \tau)]$, is displayed in Fig. 3c for the phantom experiment sketched in Fig. 3a. Although it displays a 5D structure as the canonical matrix $R_{rr}(\tau)$, the de-scanned coordinates $\Delta x$ and $\Delta z$ can display less sampling points than the coordinates $x$ and $z$. It is therefore a convenient basis for storing the ultrasound data.

Fig. 3d explains how the information contained in the-scanned matrix will be displayed in the following. A point $r_{in}$ marked by a black arrow (Fig. 3c) is chosen for illustration purposes. The corresponding column $(\Delta \mathbf{r})$ is extracted and reshaped in cartesian coordinates $(\Delta x, \Delta z)$, which allows us to visualize the spatio-temporal responses between the virtual source created at this point and all the surrounding virtual receivers. Fig. 3d and e show the
amplitude of the field obtained for a point $\mathbf{r}_{\text{in}}$ corresponding to a bright point-like scatterer or in the speckle (right). In the latter case, the field is modulated by the random reflectivity of the scatterers located at points $\mathbf{r}_{\text{in}}$ and $\mathbf{r}_{\text{out}}$ used to probe the field. We will see later that we can overcome the random fluctuations in reflectivity to recover the wave that would be associated with a coherent source located at $\mathbf{r}_{\text{in}}$, even in speckle regime.

We will now investigate in greater detail the information contained in the time-focused reflection matrix $\mathbf{R}_{\Delta\mathbf{r}}(\tau)$ in the case of each of the experiments introduced above (Fig. 1). In particular, we will show how this matrix can be used to probe the spatio-temporal focusing process in the medium. To do this, we proceed in two stages. First, we consider the case where the virtual source $\mathbf{r}_{\text{in}}$ coincides with a bright scatterer in the medium. We use a set of virtual reflectors $\mathbf{r}_{\text{out}}$ located around $\mathbf{r}_{\text{in}}$ to probe the reflected field from this source. In this case, the bright scatterer acts as an acoustic guide star located in the medium. However, such a strong reflector is not often available in vivo and not everywhere in any case. This is why, in a second step, we will consider the more complex case of a virtual source located in speckle. In this common regime, we will see that it is also possible to probe the focusing process. To do this, we will use an iterative phase reversal algorithm to synthesize a coherent virtual reflector from a set of incoherent sources $\mathbf{r}_{\text{in}}$ generated in the speckle.

## II. SELF-PORTRAIT OF A WAVE REFLECTED BY A BRIGHT SCATTERER

In ultrasound imaging, isolated bright scatterers can be used as guide stars for adaptive focusing. Here we show how the wave-field reflected by such guide star can be monitored thanks to virtual receivers synthesized in speckle. The movie of this wave-field bears the signature of the axial distortions and reverberations undergone by the wave across the medium.

### A. Correct propagation model

First, let us consider the “ideal” case of the ultrasonic phantom experiment with a constant wave speed $c_0 = 1542 \text{ m.s}^{-1}$ (Fig. 1a). We are considering a virtual source $\mathbf{r}_{\text{in}}$ coinciding with the bright scatterer indicated by an orange cross in Fig. 4a, as well as a set of virtual receivers distributed around this source, as delimited by the blue rectangle in Fig. 4a. The
column of the matrix corresponding to this point $\mathbf{r}_{\text{in}}$ gives access to the propagation film of the refocused wave associated with this virtual source. This field is represented in the de-scanned coordinate system $(\Delta x, \Delta z)$ for several lapse times $\tau$ (Figs. 4b-f), as explained in Fig. 3d. In this configuration, the wave reflected by the nylon rod clear predominates over the echoes produced by surrounding sub-resolved scatterers. We can therefore observe a coherent wave propagating in the medium, from deep $z$ towards the probe. The wave is converging towards $\mathbf{r}_{\text{in}}$ at negative times $\tau$, then focuses at time $\tau = 0$ at $(\Delta x, \Delta z) = (0, 0)$, and finally diverges at positive times $\tau$. The anti-causal component of this coherent wave (Figs. 4c,d) corresponds to the pressure wave-field that would produce the virtual source at $\mathbf{r}_{\text{in}}$ (Figs. 4e) in a transmission configuration. The causal part of this coherent wave corresponds to the wave-field induced by this virtual source (Figs. 4f,g). The time focused reflection matrix thus provides access to the propagation film of the reflected wave when a bright scatterer is positioned at $\mathbf{r}_{\text{in}}$. This statement is valid when the velocity model corresponds to the actual speed-of-sound distribution inside the propagation medium. However, in practice, it is rare to have prior knowledge of the wave velocity in the medium, so the model is generally only a more or less refined approximation of the actual velocity.

## B. Probing axial aberrations

Let us now examine how the time focused reflection matrix can be also used to highlight a focusing process degraded by an erroneous velocity model. The mismatch between the $c_0$ model and the actual velocity distribution in the medium implies a shift between the isochronous volume and the spatial focal plane (Fig. 2b). This results in transverse spreading of the focal spots and a degradation of the resolution and contrast of ultrasound images (Fig. 1d). In this section, we will see how this shift between the isochronous volume and the focal plane manifests itself on the time-focused reflection matrix. In particular, we will show how this observable is a suitable tool for probing this shift and, ultimately, for resolving it.

As in Fig. 1d, the case of an incorrect velocity model is artificially reproduced by deliberately constructing the focused reflection matrix with a velocity model $c_0 = 1440 \text{ m.s}^{-1}$ that differs from the speed-of-sound $c = 1542 \text{ m.s}^{-1}$ in the phantom. As in the previous section, we extract the focused wave-field corresponding to a virtual source $\mathbf{r}_{\text{in}}$ coinciding with a bright scatterer (Fig. 4c). As for a correct propagation model (Fig. 4b), the propagation
<figure><img src="image_5.png" /><figcaption>FIG. 4. Coherent wave-packet associated with a bright scatterer inside an homogeneous phantom. a,b Ultrasound image of the phantom beamformed at the phantom speed-of-sound ($c_0 = 1540 \text{ m.s}^{-1}$) and at a wrong speed-of-sound ($c_0 = 1440 \text{ m.s}^{-1}$), respectively: The position of the virtual source $\mathbf{r}_{\text{in}}$ corresponds to a bright scatterer indicated by an orange cross and the set of virtual receivers $\mathbf{r}_{\text{out}}$ surrounding it is delimited by a blue rectangle. c-g, Normalized amplitude of the wave-field $R(\Delta\mathbf{r}, \mathbf{r}_{\text{in}}, \tau)$ produced by the virtual source $\mathbf{r}_{\text{in}}$ at the phantom speed-of-sound (a). h-l Same as in (c-g) but for a wrong speed-of-sound (b). In panels (c-l), each focused wave-field is displayed in de-scanned coordinates for various lapse times $\tau$. The propagation of the coherent wave-packet shown in panels (c-g) and (h-l) are also displayed in Supplementary Movies 1 and 2.</figcaption></figure>

film reveals a coherent wave propagating towards the probe at a wave speed $c_0$ (Fig. 4d). However, this movie bears the mark of the mismatch between the velocity model and the true speed-of-sound. Wave focusing is observed neither at the desired location ($\Delta z^* = 2.5 \text{ mm}$) nor at the expected ballistic time ($\tau^* = -1.8 \text{ µs}$). At the ballistic time $\tau = 0$, the wave associated with the strong scatterer is therefore out-of-focus, which explains the distorted shape of the target on the confocal image (Fig. 4c).

To explain quantitatively how this shift in focusing time and depth is related to the phase velocity mismatch, we can express the de-scan reflection matrix associated with a point-like
scatterer of reflectivity $\gamma_s$ at position $(x_s, z_s)$. To do so, the focused reflection matrix will be expressed in the frequency domain: $\bar{\mathbf{R}}_{rr}(f) = \int d\tau \mathbf{R}_{rr}(\tau) \exp(-j2\pi f\tau)$. The coefficients of the monochromatic focused reflection matrices can be theoretically expressed as follows$^{12,14}$:

$$ \bar{R}(x_{\text{out}}, x_{\text{in}}, z_{\text{out}}, z_{\text{in}}, f) = \gamma_s \bar{B}(f) H_{\text{out}}(x_s, z_s, x_{\text{out}}, z_{\text{out}}, f) H_{\text{in}}(x_{\text{in}}, z_s, x_{\text{in}}, z_{\text{in}}, f), \quad (4) $$

where $\bar{B}(f)$ accounts for the probe bandwidth. $H_{\text{in/out}}(x_s, z_s, x_{\text{in/out}}, z_{\text{in/out}}, f)$ corresponds to the monochromatic PSF, i.e., the spatial distribution of the focal spot amplitude at the scatterer position when trying to focus on the point $\mathbf{r}_{\text{in/out}} = (x_{\text{in/out}}, z_{\text{in/out}})$. Under the paraxial approximation, the PSF can be, in first approximation, decomposed as a product of a frequency-invariant enveloppe and a depth-oscillating phase term (see Appendix A):

$$ H_{\text{in/out}}(x, z, x_{\text{in/out}}, z_{\text{in/out}}, f) = \underbrace{H_{\text{in/out}}^{(0)}(x - x_{\text{in/out}}, cz/c_0 - z_{\text{in/out}}, f_c)}_{\text{transverse focusing}} \underbrace{e^{\frac{i2\pi f}{c_0} \left(\frac{z_{\text{in/out}} - z}{c}\right)}}_{\text{axial propagation}} \quad (5) $$

While the first term in Eq. 5 accounts for the transverse focusing process at the focusing depth,

$$ z_{\text{in/out}}^{*} = c_0 z_{\text{in/out}} / c, \quad (6) $$

the second term accounts for the axial propagation of the focused wave-field along the z-axis. The scatterer depth $z_s$ is controlled by the time-of-flight expected for a scatterer at the targeted depth $z_{\text{in}}$ under a wave velocity model $c_0$: $z_s = \frac{c}{c_0} z_{\text{in}}$. Using this expression of $z_s$ and injecting Eq. 5 into Eq. 4 leads to the following expression of the focused R-matrix coefficients in de-scanned coordinates in the time domain,

$$ R(\Delta x, \Delta z, z_{\text{in}}, \tau) \propto \underbrace{H_{\text{in}}^{(0)}(0, \Delta z^*, f_c)}_{\text{constant}} \underbrace{H_{\text{out}}^{(0)}(\Delta x, \Delta z^* - \Delta z, f_c)}_{\text{spatial focusing}} \underbrace{B\left(\tau + \frac{\Delta z}{c_0}\right)}_{\text{time focusing}}. \quad (7) $$

with

$$ \Delta z^* = \left( \frac{c^2}{c_0^2} - 1 \right) z_{\text{in}}, \quad (8) $$

the focusing depth in the de-scanned frame. The temporal term in Eq. 7 proves the upward propagation of the wave-front at the supposed speed-of-sound $c_0$. The propagation films thus provide a self-portrait of the focusing process under the prism of the wave velocity model. Despite this bias, the observed focusing time depends on the actual velocity $c$ in the propagation medium:

$$ \tau^* = -\Delta z^* / c_0 = \frac{z_{\text{in}}}{c_0} \left( 1 - \frac{c^2}{c_0^2} \right). \quad (9) $$
<figure><img src="image_6.png" /><figcaption>FIG. 5. Coherent wave-packet associated with a bright scatterer in presence of reverberations. a, Ultrasound image of the phantom through the Plexiglas layer (Fig. 1e): The position of the virtual source $r_{\text{in}}$ corresponds to a bright scatterer indicated by an orange circle and the set of virtual receivers $r_{\text{out}}$ surrounding it is delimited by a blue rectangle. b-f, Normalized amplitude of the wave-field produced by the virtual source $r_{\text{in}}$ in de-scanned coordinates displayed for various lapse times $\tau$. The propagation of the coherent wave-packet is also displayed in Supplementary Movie 3.</figcaption></figure>

A numerical application using the apparent depth value $z_{\text{in}} = 18.9$ mm of the scatterer selected in Fig. 4c yields values of $\Delta z^* \sim 2.7$ mm and $\tau^* \sim -1.88$ µm, which are consistent with the experimental values previously extracted from Fig. 4d.

The propagation movie of the reflected wave therefore provides a relevant information about the discrepancy between this velocity model and reality. Interestingly, we now show how such an observable can be used to probe reverberations induced by a layer displaying a strong impedance mismatch with respect to its environment.

### C. Probing reverberations

We now focus on the influence of a reverberating plate on the focusing process. To do so, we consider the experiment with the ultrasonic phantom covered with a Plexiglas plate of thickness $L = 6$ mm and speed-of-sound $c_r = 2750$ m.s$^{-1}$ (Fig. 1e). The focused reflection matrix is synthesized assuming an homogeneous velocity model $c_0$ equal to the speed-of-sound
in the phantom ($c = 1540 \text{ m.s}^{-1}$). This choice ensures, in first approximation, isoplanicity over the whole phantom. Only the finite size of the probe limits the spatial invariance of the imaging configuration.

The reflected wave-field associated with the bright scatterer indicated in Fig. 5a is displayed at different lapse times in Fig. 5b-g. It clearly highlights the presence of multiple reflections associated with a period of ~2 µs. Those reverberations a priori correspond to waves that have been multiply reflected between the probe and the phantom surface, as shown schematically in Fig. 6. Furthermore, we can observe in the film that the coherent ballistic wave focuses at a time $\tau = -5.5 \text{ µs}$, well before the ballistic time expected by the homogeneous velocity model. This difference between the observed focusing time and the ballistic time can be explained by the large gap between the velocity model and the speed-of-sound distribution in the current experiment. The large speed-of-sound in the Plexiglas layer implies this earlier arrival of the direct echo compared to what would be expected in an homogeneous phantom.

To analyze more quantitatively our experimental observation, a theoretical analysis of reverberations is now performed by considering a bi-layered medium consisting in an imaging volume of speed-of-sound $c$ placed behind a reverberating layer of thickness $L$ and wave velocity $c_r$ (Fig. 6). The reflected wave-field can be decomposed into the set of possible scattering paths that the incident and reflected waves can undertake during their travel between the probe and the isochronous volume (Fig. 6). Under a matrix formalism, it can be written as follows:

$$ \mathbf{R}_{\Delta\mathbf{r}}(\tau) = \sum_{l=0}^{\infty} \sum_{m=0}^{\infty} \alpha^{2(m+l)} \mathbf{R}_{\Delta\mathbf{r}}^{(l,m)}(\tau). \qquad (10) $$

where $\mathbf{R}_{\Delta\mathbf{r}}^{(l,m)}$ accounts for each component of the reflection matrix associated with $2l$ internal reflections on the incident way and $2m$ internal reflections on the return path. The coefficient $\alpha$ accounts for the internal reflection coefficient in the reverberating plate. This reflection coefficient weights the sum over scattering paths in Eq. 10 and controls the number of reverberations we can observe in practice (Fig. 5).

In presence of bright scatterer at $(x_s, z_s)$, each component of the focused reflection matrix (Eq. 10) can be expressed as follows:

$$ \bar{\mathbf{R}}_{\Delta\mathbf{r}}^{(l,m)}(f) = \gamma_s \bar{B}(f) H_{\text{out}}^{(m)}(x_s, z_s, x_{\text{out}}, z_{\text{out}}, f) H_{\text{in}}^{(l)}(x_{\text{in}}, z_s, x_{\text{in}}, z_{\text{in}}, f). \qquad (11) $$

where $H_{\text{in}}^{(l)}(x_s, z_s, x_{\text{in}}, z_{\text{in}}, f)$ and $H_{\text{out}}^{(m)}(x_s, z_s, x_{\text{out}}, z_{\text{out}}, f)$ stand for the $l^{\text{th}}$ and $m^{\text{th}}$ reverbera-
<figure><img src="image_7.png" /><figcaption>FIG. 6. Diagram of the different wave propagation paths in a two-layer medium. a Example of ballistic path trajectory between the probe and a scatterer at $z_s$ through a reverberating layer of thickness $L$. b Example of first-order reflection path trajectory associated with the reflection matrix component $\mathbf{R}^{(1,0)}$ (Eq. 10): The incident wave is reflected twice at the interfaces of the reverberating layer. c Example of first-order reflection path trajectory associated with the component $\mathbf{R}^{(0,1)}$: The up-going wave is reflected twice at the interfaces of the reverberating layer. d Example of second-order reflection path trajectory associated with the component $\mathbf{R}^{(1,1)}$: The incident and reflected waves are reflected twice at the interfaces of the reverberating layer. e Example of second-order reflection path trajectory associated with the component $\mathbf{R}^{(2,0)}$: The incident wave is reflected four times at the interfaces of the reverberating layer.</figcaption></figure>

tion order of the input and output PSFs, respectively. As before, each PSF can be, under the paraxial approximation, decomposed as a product of a frequency-invariant envelope and a depth-oscillating phase term (see Appendix B):

$$H_{\text{in/out}}^{(n)}(\mathbf{r}, \mathbf{r}_{\text{in/out}}, f) = e^{\underbrace{i2\pi f \left( \frac{z-L-z_{\text{in/out}}}{c} + \frac{(2n+1)L}{cR} \right)}_{\text{axial propagation}}} \underbrace{H_{\text{out}}^{(0)}(x-x_{\text{in/out}}, z-z_{\text{in/out}}^{\star(n)}, f_c)}_{\text{transverse focusing}} \quad (12)$$

with

$$z_{\text{in/out}}^{(n)\star} = z_{\text{in/out}} + L - (2n+1)L \frac{c_R}{c}, \quad (13)$$

the depth of the focusing plane for the $l^{th}$ reverberation order. Injecting the Eq. 12 into Eq. 10 leads to the following expression for the coefficients of the time focused $\mathbf{R}$-matrix
expressed in de-scan coordinates:

$$R^{(l,m)}(\Delta x, \Delta z, x_{\text{in}}, z_{\text{in}}, \tau) \propto \underbrace{H_{\text{in}}^{(0)}(0, \Delta z_l^*, f_c)}_{\text{constant}} \underbrace{H_{\text{out}}^{(0)}(\Delta x, \Delta z_m^* - \Delta z, f_c)}_{\text{spatial focusing}} \underbrace{B\left(\tau + \frac{\Delta z}{c} - \frac{2(l+m)L}{c_R}\right)}_{\text{time focusing}} \quad (14)$$

with

$$\Delta z_n^* = 2nL \frac{c_R}{c} + L \left( \frac{c_R}{c} - \frac{c}{c_R} \right), \quad (15)$$

the de-scanned focusing depth of the n-th reverberation order.

The temporal term in Eq. 14 accounts for a time comb of multiply-reflected echoes behind the ballistic component, with a period equal to the round-trip time $2L/c_R$ inside the reverberating layer. The observed period of 2 µs in Fig. 5 is not consistent with reverberation in the Plexiglas layer that would imply a reverberation time of 4 µs.. It is more probably the signature of reverberation inside the acoustic lens ($c_r \sim 1000$ m/s and $L \sim 1$ mm) that lies between the probe and the Plexiglas.

Each of the reverberated wave-front focuses at a de-scan depth $\Delta z_m^*$ (Eq. 15) that only depends on the reverberation order $m$ of the output path. For the ballistic wave-front, the corresponding focusing time shift is therefore given by

$$\tau_0^* = -\frac{\Delta z_0^*}{c} = \frac{L}{c_R} \left( 1 - \frac{c_R^2}{c^2} \right) \quad (16)$$

The measured value $\tau_0^* \sim -5.5$ µs in Fig. 5 is this time consistent with a Plexiglas layer of thickness $L \sim 6$ mm. Indeed, the aberration undergone by the ballistic wave-front is mainly induced by the Plexiglas layer (and not by the acoustic lens that would yield $\tau_0^* \sim -0.5$ µs).

Through this example, we show that the reflection matrix is an ideal tool for analyzing the focusing process in the medium. By decoupling the imaging time from the ballistic echo time, it allows us to probe the temporal dispersion of the echoes and, for example, highlight the presence of multiple reflections. In the third article of the series¹⁷, we will see that this information is valuable and can be used to compensate for multiple propagation paths. In that respect, another important feature highlighted by Eqs. 15 and 16 is the confirmation that the effect of the reverberating layer does not depend on the target depth $z_{\text{in}}$. All points below this layer experience the same focusing defect and reverberation phenomena. Such a an isoplanicity will be an important property to extract a self-portrait of the focusing process in the speckle regime.
Indeed, for the moment, we have limited ourselves to studying the focused reflection matrix associated with a bright scatterer. As pointed out in the introduction, real media do not contain isolated bright scatterers as in a phantom. For instance, in medical ultrasound, the biological tissues are mainly composed of a random distribution of sub-resolved scatterers generating ultrasonic speckle. For this reason, in the next section, we will revisit the temporal study of the focused **R**-matrix carried out in this section, extending it to the case of a speckle medium.

### III. SELF-PORTRAIT OF A COHERENT WAVE PACKET IN SPECKLE

In the speckle regime, the reflected wave-field is the result of random interference between echoes generated by un-resolved scatterers. An incident focal spot can thus be seen as an incoherent virtual source. While the wave-fields produced by such virtual sources are generally seen as spatially- and temporally-incoherent, we here show how to combine them in order to synthesize a coherent reflector and monitor the propagation of a coherent wave-packet despite disorder.

#### A. Wave-packet associated with one virtual source

We thus turn our attention to the more complex situation where the time-focused reflection matrix is constructed for a virtual source $\mathbf{r}_{\text{in}}$ and a set of virtual receivers $\mathbf{r}_{\text{out}}$ surrounding it, all located in a speckle zone, as shown in Fig. 7a. First, we again consider the ideal case of the ultrasonic phantom experiment in which the velocity model used to form the matrix $\mathbf{R}_{\Delta\mathbf{r}}(\tau)$ is correct ($c = c_0$). Propagation of the wave induced by a given virtual source $\mathbf{r}_{\text{in}}$ in speckle can be displayed by reshaping the corresponding column of $\mathbf{R}_{\Delta\mathbf{r}}(\tau)$ in de-scanned coordinates $(\Delta x, \Delta z)$ at different times-of-flight $\tau$ (Fig. 7b-f). As the reflected wave induced by a bright scatterer (Fig. 4), a coherent wave propagating upwards seems to emerge (see the white arc in Fig. 7b-f). The coherent beam seems to narrow to a maximum at the ballistic time ($\tau = 0$) and at the expected focusing point $(\Delta\mathbf{r} = \mathbf{0})$, suggesting its focus (white arrows in Fig. 7d). However, the observation of the focusing process is hampered by the random reflectivity of the sub-resolved scatterers that modulates the amplitude of the back-propagated wave-field. Furthermore, the coherent wave is also partially masked
<figure><img src="image_8.png" /><figcaption>FIG. 7. Coherent wave-packet associated with a virtual source in the speckle region. a, Ultrasound image of the phantom (Fig. 1e): The position of the virtual source **r**<sub>in</sub> is indicated by an orange cross and the set of virtual receivers **r**<sub>out</sub> surrounding it is delimited by a blue rectangle. **b-f**, Normalized amplitude of the wave-field produced by the virtual source **r**<sub>in</sub> in de-scanned coordinates displayed for various lapse times τ. The propagation of the coherent wave-packet is also displayed in Supplementary Movie 4. The white arches highlight the coherent wave-front propagating upwards that can be discerned despite the random reflectivity of the speckle grains.</figcaption></figure>

by spurious echoes generated by out-of-plane scatterers located in the illumination cone,
upstream and downstream of the virtual source, and, possibly, by multiply-scattered echoes.

To circumvent this issue and isolate the coherent wave packet, one can take advantage of
the spatial invariance of the problem and combine a set of propagation films associated with
different virtual sources **r**<sub>in</sub> belonging to the same isoplanatic patch.

B. Extraction of the coherent wave packet by singular value decomposition

A set of virtual sources **r**<sub>in</sub> belonging to the same isoplanatic patch are now considered
(Fig. 8a). Each of these sources should give rise to spatially correlated wave-fronts because
of the shift-shift memory effect (Fig. 8b). However, each virtual source is also modulated by
the random distribution of sub-resolved scatterers lying into each incident focal spot. The
idea is thus to exploit the spatial correlations between each focused wave-field in order to:
(i) synthesize a coherent guide star from these different realizations of disorder (Fig. 8c); (ii)
extract a coherent wave packet equivalent to the one produced by a bright scatterer in Fig. 4.
<figure><img src="image_9.png" /><figcaption>FIG. 8. **Synthesis of a coherent guide star from a set of virtual sources r<sub>in</sub> by iterative time reversal.** a The focused reflection matrix **R**<sub>rr</sub>(τ) = [R(r<sub>out</sub>, r<sub>in</sub>, τ)] contains the wave-fronts induced by a set of virtual sources synthesized in speckle. b Express the reflection matrix in de-scan frame, **R**<sub>Δ**r**</sub>(τ) = [R(r<sub>out</sub> − r<sub>in</sub>, r<sub>in</sub>, τ)], amounts to a lateral shift of virtual sources at the center of the field-of-view. c An iterative time reversal process applied to **R**<sub>Δ**r**</sub>(τ) enables a rephasing of each virtual source in order to synthesize a coherent guide star. In practice, this operation can be performed by a SVD of the concatenated matrix **R**<sub>{Δ**r**,τ}<sub>**r**</sub>. Its first eigenvector **U**<sub>1</sub> = [U<sub>1</sub>(Δ**r**, τ)] provides the coherent wave-field induced by this guide star.</figcaption></figure>

Previous works have shown how a SVD process can be fruitful to synthesize such a coherent
guide star<sup>11,18-20</sup>. However, those previous demonstrations were limited to a time-gated
wave-field. Here, we show how this process can be extended to a time-dependent wave-field
in order to extract a coherent spatio-temporal wave-packet in speckle.

To do so, we first select the speckle area over which we want to evaluate the focusing
process. This area should be as small as possible to guarantee a maximum isoplanicity.
However, it should encompass a sufficiently large number of independent virtual sources in
order to extract the coherent wave packet and get rid of speckle fluctuations, out-of-focus
echoes and multiply-scattered waves. Selected areas of virtual sources are superimposed to
the ultrasound images in Figs. 9a and e for a correct and incorrect wave velocity model,
respectively It should be noted that the isoplanicity condition for c ≠ c₀ is more restrictive
since here the velocity error induces aberrations that are not invariant with respect to depth
z. For this reason, the selected area for an incorrect wave velocity model (Fig. 9e) has been
restricted axially and extended laterally compared to the reference case (Fig. 9a).

Once the area of virtual sources selected, the time-dependent de-scan reflection matrix is
first reshaped in 2D by concatenating the dimensions $\Delta\mathbf{r} = (\Delta x, \Delta z)$ and $\tau$:

$$ \mathbf{R}_{\{\Delta\mathbf{r},\tau\},\mathbf{r}} = [R(\{\Delta\mathbf{r}_{\text{out}}, \tau\}, \mathbf{r}_{\text{in}})]. \quad (17) $$

Each column of the matrix corresponds to a spatio-temporal wave-field induced by a given virtual source $\mathbf{r}_{\text{in}}$ as the one displayed in Fig. 7. To synthesize a coherent guide star from these different realizations over disorder, a singular value decomposition (SVD) of the concatenated matrix $\mathbf{R}_{\{\Delta\mathbf{r},\tau\},\mathbf{r}}$ is performed such that:

$$ \mathbf{R}_{\{\Delta\mathbf{r},\tau\},\mathbf{r}} = \mathbf{U} \times \mathbf{\Sigma} \times \mathbf{V}^{\dagger} \quad (18) $$

or, in terms of coefficients:

$$ R(\{\Delta\mathbf{r}_{\text{out}}, \tau\}, \mathbf{r}_{\text{in}}) = \sum_{p} \sigma_{p} U_{p}(\Delta\mathbf{r}_{\text{out}}, \tau) V_{p}^{*}(\mathbf{r}_{\text{in}}). \quad (19) $$

The matrix $\mathbf{\Sigma}$ is diagonal and its diagonal coefficients contains the singular values of $\mathbf{R}_{\{\Delta\mathbf{r},\tau\},\mathbf{r}}$ ranged in decreasing order ($\sigma_1 > \cdots > \sigma_N$). The matrix $\mathbf{U}$ and $\mathbf{V}$ are unitary matrices whose columns $\mathbf{U}_p = [U_p(\Delta\mathbf{r}_{\text{out}}, \tau)]$ and $\mathbf{V}_p = [V_p(\mathbf{r}_{\text{in}})]$ are the singular vectors of $\mathbf{R}_{\{\Delta\mathbf{r},\tau\},\mathbf{r}}$. Among all eigenstates, we are interested in the first one that maximizes the correlation of the wave-field induced by each virtual source. Figures 9b-d show the first singular vectors of $\mathbf{R}_{\{\Delta\mathbf{r},\tau\},\mathbf{r}}$ for the area of virtual sources displayed in Fig. 9a and beamformed using a correct wave velocity model ($c_0 = c$). The first input eigenstate $\mathbf{V}_1$ (Fig. 9b) contains the weighting coefficients that should be applied to each virtual source in order to synthesize a coherent guide star. The first output eigenstate $\mathbf{U}_1$ (Fig. 9c) displays the wave-field generated by this guide star. Equivalent to the reflected wave-field induced by a bright scatterer (Fig. 4c), this first result demonstrates our ability of synthesizing a guide star in speckle and probing the space-time evolution of the associated wave packet. The SVD has allowed us to overcome the random reflectivity of speckle grains and remove parasitic echoes observed in Fig. 7b. A coherent wave propagating towards the probe is again observed. Nevertheless, the axial resolution of the wave packet seems to be degraded by the SVD process. This deterioration is even worse when a mismatch exists between the real speed-of-sound distribution in the medium and the wave velocity model (Fig. 9g). Moreover, the focusing depth and time of the wave-packet seem to differ from the one induced by a bright scatterer in the latter case. To understand this discrepancy and the limits of the SVD process, an analytical study of this operation is now provided.
<figure><img src="image_10.png" /><figcaption>FIG. 9. Propagation movie associated with the coherent wave packed extracted by iterative time reversal from a set of virtual sources in the speckle. a-d Correct speed-of-sound model ($c_0 = 1540 \text{ m.s}^{-1}$): a The considered virtual sources $\mathbf{r}_{\text{in}} = (x_{\text{in}}, z_{\text{in}})$ and their neighboring receivers are considered in the orange and blue rectangles on the ultrasound image. b The input singular vector $\mathbf{V}_1$ yields the amplitude (top) and phase (bottom) of the weighting coefficient that shall be applied to each virtual source to synthesize the coherent guide star. c The first output singular vector $\mathbf{U}_1$ of the re-shaped 2D matrix $\mathbf{R}_{\{\Delta\mathbf{r},\tau\}\mathbf{r}}$ yields the propagation movie of the coherent wave-packet shown here in absolute value. The propagation of the coherent wave-packet is also displayed in Supplementary Movie 4. d Spatio-temporal Fourier transform of $\mathbf{U}_1$. e-h Same as in panels a-d but for an incorrect wave velocity model ($c_0 = 1440 \text{ m.s}^{-1}$). The propagation of the coherent wave-packet is also displayed in Supplementary Movie 5.</figcaption></figure>

### C. Analytical expression of the focused reflection matrix in speckle

The SVD result of Fig. 9 can be intuitively understood by considering the concatenated reflection matrix in the frequency domain, $\overline{\mathbf{R}}_{\{\Delta\mathbf{r},f\}\mathbf{r}}$. Indeed, the Fourier transform of Eq. 19 yields:

$$ \overline{R}(\{\Delta\mathbf{r}_{\text{out}}, f\}, \mathbf{r}_{\text{in}}) = \sum_{p} \sigma_{p} \overline{U}_{p}(\Delta\mathbf{r}_{\text{out}}, f) V_{p}^{*}(\mathbf{r}_{\text{in}}). \qquad (20) $$
Under a single scattering assumption, the coefficients of the monochromatic focused reflection matrices can be theoretically expressed as follows¹²,¹⁴:

$$ \bar{R}(x_{\text{out}}, x_{\text{in}}, z_{\text{out}}, z_{\text{in}}, f) = \iint dzdx \ H_{\text{out}}(x, z, x_{\text{out}}, z_{\text{out}}, f)\gamma(x, z)H_{\text{in}}(x, z, x_{\text{in}}, z_{\text{in}}, f), \quad (21) $$

where $\gamma(x, z)$ represents the reflectivity of the medium at point $\mathbf{r} = (x, z)$. Injecting the paraxial expression of the input and output PSFs (Eq. 5) into Eq. 21 leads to the following expression of the focused $\mathbf{R}$-matrix coefficients:

$$ \bar{R}(x_{\text{out}}, x_{\text{in}}, z_{\text{out}}, z_{\text{in}}, f) = B(f) \iint dzdx \ H_{\text{out}}^{(0)}(x - x_{\text{out}}, c(z - z_{\text{out}}^{*})/c_0, f_c) \\ \times \gamma(x, z) H_{\text{in}}^{(0)}(x - x_{\text{in}}, c(z - z_{\text{in}}^{*})/c_0, f_c) e^{i2\pi f \left( \frac{z_{\text{in}}+z_{\text{out}}}{c_0} - \frac{2z}{c} \right)}. \quad (22) $$

In first approximation, the incident beam can also be described as perfectly focused at the focusing point $\mathbf{r}_{\text{in}}^* = (x_{\text{in}}, z_{\text{in}}^*)$:

$$ H_{\text{in}}^{(0)}(x - x_{\text{in}}, c(z - z_{\text{in}}^{*})/c_0, f_c) \sim \delta(x - x_{\text{in}})\delta(z - z_{\text{in}}^{*}), \quad (23) $$

This approximation amounts to considering that the virtual source synthesized at $\mathbf{r}_{\text{in}} = (x_{\text{in}}, z_{\text{in}}^*)$ as point-like. The expression of the focused reflection matrix in the concatenated de-scanned basis then becomes:

$$ \bar{R}(\{\Delta x, \Delta z, f\}, \{x_{\text{in}}, z_{\text{in}}\}) \sim \underbrace{\gamma(x_{\text{in}}, z_{\text{in}}^*)}_{V_1^*(x_{\text{in}}, z_{\text{in}})} \overline{B(f) H_{\text{out}}^{(0)}(-\Delta x, -\Delta z, f_c)} e^{i2\pi f \left[ \frac{2z_{\text{in}}}{c_0} \left(1 - \frac{c_0^2}{c^2}\right) + \frac{\Delta z}{c_0} \right]} . \quad (24) $$

Comparison with Eq. 20 shows that $\bar{R}_{\{\Delta\mathbf{r},f\}\mathbf{r}}$ is thus of rank 1. On the one hand, the first input singular vector $\mathbf{V}_1$ directly provides the medium reflectivity at the focusing depth: $V_1(x_{\text{in}}, z_{\text{in}}) = \gamma(x_{\text{in}}, z_{\text{in}}^*)$ (Figs. 9b,f). On the other hand, the first output singular vector directly provides the output PSF at the central frequency, $H_{\text{out}}^{(0)}(-\Delta x, -\Delta z, f_c)$, modulated by the frequency bandwidth of the recorded signals and a phase term that accounts for the axial propagation of the coherent wave-field. An inverse Fourier transform of the first eigenvector $\bar{\mathbf{U}}_1(f)$ provides an analytical expression for the coherent wave-packet extracted by the SVD process in the time domain:

$$ U_1(\Delta x, \Delta z, \tau) = H_{\text{out}}^0(-\Delta x, -\Delta z, f_c) B\left(\tau + 2\tau^* + \frac{\Delta z}{c_0}\right). \quad (25) $$

This result confirms that the vector $\mathbf{U}_1$ provides access to a coherent spatio-temporal wave-packet in a speckle environment with a spatial envelope accounting for the output PSF at
the central frequency and a time shift $2\tau^*$ (Eq. 9) accounting for the time delay experienced by the incident and reflected waves with respect to the ballistic time expected under the wave velocity assumption.

These analytical predictions are in good agreement with the spatio-temporal wave packets displayed in Figs. 9c and g. For a correct wave velocity model ($c_0 = 1540 \text{ m.s}^{-1}$, Fig. 9c), the wave-packet actually focuses at the origin ($\Delta z = 0$) and at the expected ballistic time ($\tau = 0$). For an incorrect wave velocity model ($c_0 = 1440 \text{ m.s}^{-1}$, Fig. 9g), the focus is still observed at $\Delta z = 0 \text{ mm}$ and close to the expected focusing time ($2\tau^* \sim 3.7 \ \mu\text{s}$). Nevertheless, the point-like virtual source hypothesis made for analytical tractability does not provide a clue for the loss of axial resolution observed in Figs. 9c and g.

To go beyond this analytical calculation, one can exploit the analogy between the SVD and iterative time reversal processing$^{11,20}$. The first eigenstate can be seen as the result of an iterative time reversal experiment that tends to maximize the energy back-scattered by our virtual guide star. As the virtual source is, in reality, not point-like but displays a depth extension dictated by the depth-of-field $\delta z \sim 7\lambda(z_{\text{in}}/D)^2$ (with $D$ the array extension), the SVD process converges towards a wave-front of limited frequency-bandwidth $\delta f \sim c/(2\delta z)$ around the central frequency $f_c$. The first time-reversal invariant extracted by the SVD process optimizes the energy coming from the elongated guide star. This feature is confirmed by the spatio-temporal Fourier transform of the first eigenvector $\mathbf{U}_1$. While the angular distribution of the wave-field is quite homogeneous in the transverse spatial frequency domain, its energy is peaked around the central temporal frequency $f_c$ (Figs. 9d,h). The bandwidth $\delta f$ is even more limited when the wave velocity model is not correct (Figs. 9h). Indeed, the axial distortions undergone by the incident focal spot leads to an axial extension of the guide star which, in turn, reduces the bandwidth of the first output singular vector $\mathbf{U}_1$.

### D. Extraction of the coherent wave packet by iterative phase reversal

To circumvent that issue, an iterative phase reversal process can be operated from the $k$-space thanks to the distortion matrix concept$^{11}$. Originally proposed for a time-gated reflection matrix$^{15}$, it is here applied to the broadband reflection matrix. To this end, the reflection matrix is projected from the de-scanned basis to the plane wave basis ($k_x$) at each frequency $f$. This operation gives direct access to the distortion matrix in the spatio-temporal
Fourier domain:

$$D(k_x, \Delta z, \mathbf{r}_{\text{in}}, f) = \iint dt d\Delta x \, R(\Delta x, \Delta z, \mathbf{r}_{\text{in}}, \tau) e^{ik_x \cdot \Delta x} e^{-i2\pi f \tau}. \quad (26)$$

The distortion matrix contains the deviation of each angular component of reflected waves from an ideal wavefront that would be obtained in the absence of aberrations and without multiple scattering. To synthesize a coherent guide star from the set of virtual sources $\mathbf{r}_{\text{in}}$, each distorted wave-front shall be first correlated with each other. Mathematically, this operation consists in computing a correlation matrix $\mathbf{C} = [C(\{k_x, \Delta z, f\}, \{k'_x, \Delta z', f'\})]$ between each angular/depth/frequency $(k_x/\Delta z/f)$ component of the distorted wave-fronts:

$$C(\{k_x, \Delta z, f\}, \{k'_x, \Delta z', f'\}) = \sum_{\mathbf{r}_{\text{in}}} D(k_x, \Delta z, f, \mathbf{r}_{\text{in}}) D^*(k'_x, \Delta z'_{\text{out}}, f', \mathbf{r}'_{\text{in}}). \quad (27)$$

The first singular vector $\overline{\mathbf{U}}_1$ of the distortion matrix, or equivalently, the first eigenvector of $\mathbf{C}$, is the invariant of a time reversal process applied to $\mathbf{D}$, that can be expressed as follows:

$$\sigma_1^{(n)} \overline{\mathbf{U}}_1^{(n)} = \mathbf{C} \times \overline{\mathbf{U}}_1^{(n-1)}. \quad (28)$$

with $\overline{\mathbf{U}}_1 = \lim_{n \to \infty} \overline{\mathbf{U}}_1^{(n)}$ and $\sigma_1 = \lim_{n \to \infty} \sigma_1^{(n)}$. As we have seen before, such an iterative time reversal process is not optimal since it tends to maximize the back-scattered energy and therefore concentrates on the most energetic part of the frequency bandwidth (Figs. 9d,h). To circumvent this issue, Eq. 28 can be slightly modified to converge towards a wave-front effective over the whole frequency bandwidth and angular spectrum, such that

$$\mathbf{W}^{(n+1)} = \exp \left( i \arg \left\{ \mathbf{C} \times \mathbf{W}^{(n)} \right\} \right), \quad (29)$$

with $\mathbf{W}^{(0)} = [1 \cdots 1]^\top$, an arbitrarily chosen unit vector. The phase reversal invariant, $\mathbf{W} = \lim_{n \to \infty} \mathbf{W}^{(n)}$, can be applied a last time to the correlation matrix in order to recover the frequency and angular spectrum of the wave-field:

$$\overline{\mathbf{T}} = \mathbf{C} \times \mathbf{W} \quad (30)$$

The resulting wave-field $\overline{\mathbf{T}}$ is an estimator of the medium transmittance in the space-time Fourier space. In contrast with the SVD process (Figs. 9d and h), the IPR algorithm converges towards an invariant that now spans over the whole frequency bandwidth, here 5-8 MHz (Figs. 10c and f).
Inverse Fourier transforms in the $(f, k)$-space finally yield the movie of the coherent wave-packet in the spatial and time domains:

$$T(\Delta x, \Delta z, \tau) = \iint df dk_x \bar{T}(k_x, \Delta z, f) e^{-ik_x \cdot \Delta x} e^{i2\pi f \tau}. \qquad (31)$$

The result is displayed in Fig. 10b for a correct wave velocity model $c = c_0$. Comparison with the coherent wave packet resulting from the SVD (Fig. 9c) shows the drastic gain in axial resolution provided by the IPR process. Its benefit is even more striking in the aberrated case. While the SVD wave packet spreads over a few millimeters and exhibits a low contrast with respect to spurious echoes (Fig. 9g), the IPR wave-packet is highly resolved in time and space (Fig. 10e). This feature is explained by the much broader frequency bandwidth of the coherent wave packet with IPR (Fig. 10f) compared to SVD (Fig. 9h). The estimation of the focusing time can now be made with high precision and could lead to a sharp estimation of the speed-of-sound: $\hat{c} \sim 1545 \text{ m.s}^{-1}$ for $2\tau^* \simeq -3.7 \mu\text{s}$ (Eq. 9).

Nevertheless, we will see, in the second article of the series¹⁶, how we can adopt a more direct and robust strategy for local speed-of-sound estimation. Indeed, one of the limit of our method here is that it requires a region of homogeneous speckle. Moreover, this area should be large enough to encompass a sufficient number of independent realizations of disorder. The bias made on the transmittance estimator actually scales as:

$$||\delta \bar{T}||^2 \propto \frac{1}{C^2 N_{\text{in}}}, \qquad (32)$$

where $N_{\text{in}}$ is the number of diffraction-limited resolution cells in the selected area of virtual sources. $C$ is a coherence factor that is an indicator of the quality of focusing²¹, defined as the ratio between the coherent intensity and the average incoherent intensity of the back-scattered echoes:

$$C = \frac{\left| \sum_{k_x} \sum_f D(k_x, \Delta z, f) \right|^2}{\sum_{k_x} \sum_f |D(k_x, \Delta z, f)|^2} \qquad (33)$$

Physically, the bias error $||\delta \bar{T}||^2$ inversely proportional to $C^2$ (Eq. 32) reflect the following effect: the more degraded the focusing and the more blurred the guide star, the greater the bias on the estimator. One can actually show that the coherence factor is inversely proportional with respect to the number of independent spatio-temporal degrees-of-freedom exhibited by the distorted wave-front. Their order of magnitude scales as the inverse of the number of resolution cells exhibited by the focal spot along each of its dimension:

$$C \sim \frac{\delta x_0 \delta z_0 \delta t_0}{\delta x \delta z \delta t}. \qquad (34)$$
<figure><img src="image_11.png" /><figcaption>FIG. 10. Propagation movie associated with the coherent wave packet extracted by iterative phase reversal from a set of virtual sources in the speckle. Correct wave velocity model ($c_0 = 1540 \text{ m.s}^{-1}$): a The considered virtual sources $r_{\text{in}} = (x_{\text{in}}, z_{\text{in}})$ and their neighboring receivers are considered in the orange and blue rectangles on the ultrasound image. b The iterative phase reversal invariant yields the propagation movie of the coherent wave-packet $T(\Delta x, \Delta z, \tau)$ shown here in absolute value at different lapse times $\tau$. The propagation of the coherent wave packet is also displayed in Supplementary Movie 4. c Frequency-angular spectrum $\overline{\overline{T}}(k_x, \Delta z, \tau)$ of the iterative phase reversal invariant. d-f Same as in panels a-c but for an incorrect wave velocity model ($c_0 = 1440 \text{ m.s}^{-1}$). The propagation of the coherent wave packet is also displayed in Supplementary Movie 5.</figcaption></figure>

The ratio $\delta x/\delta x_0$ corresponds to the number of resolution cells exhibited by the focal spot along the lateral dimension, i.e the ratio between the lateral extension $\delta x$ of the focal spot and the diffraction-limited resolution length $\delta x_0 \sim \lambda/2 \sin \theta$, with $\theta$ the angle under which each virtual source sees the transducers array. The ratio $\delta z/\delta z_0$ corresponds to the number of resolution cells exhibited by the focal spot along the axial dimension, i.e the ratio between the axial extension $\delta z$ of the focal spot and $\delta z_0 \sim 2\lambda/\sin^2 \theta$, the depth-of-field at the position of virtual sources. At last, the ratio $\delta t/\delta_0$ corresponds to the number of time samplings
exhibited by the focal spot in the time domain, with $\delta t$ the time spreading of the focal spot and $\delta t_0 \sim \Delta f^{-1}$, the time resolution of our measurement and $\Delta f$, the frequency bandwidth.

On the one hand, Equation 32 indicates that the number of independent virtual sources should be quite large to obtain a satisfying estimation of the wave-field, especially when the velocity model is far from reality (low coherence factor). On the other hand, all these virtual sources should belong to the same isoplanatic patch, the size of the latter being very reduced in presence of distributed aberrations. This dilemma can be circumvented by reducing the space of solutions. One can consider, for instance, the time-gated reflection matrix and only focus on the spatial degrees of freedom of the wave-field. This strategy can be followed in absence of multiple reflections, as we will show in the second article of the series¹⁶ for compensating the axial distortions of a confocal image due to speed-of-sound variations.

## E. Reducing the space of solutions in presence of temporal dispersion

In presence of temporal dispersion, an alternative strategy is to first focus on the lateral/frequency degrees of freedom before retrieving, in a second step, the axial evolution of the wave-field. Indeed, the dispersion relation implies a connection between the axial and traverse component of the wave vector and the frequency: $k_x^2 + k_z^2 = (2\pi f/c)^2$. Reducing the space of solutions to the plane $(k_x, f)$ is therefore not detrimental and can facilitate the convergence of the IPR process.

This is actually the case for the reverberating plate experiment previously described in Fig. 1e. The IPR process described in Eq. 29 is unable to converge when applied to the focused reflection matrix for virtual sources contained in the speckle area displayed in Fig. 11a. To circumvent this problem, we have restricted the IPR analysis to the cross-section of the distortion matrix at $\Delta z = 0$: $D_{\{k,f\},r}(\Delta z = 0) = [D(\{k_x, f\}, \Delta z = 0, r_{in})]$. The result is a reverberation transmittance, $\overline{\overline{T}}(\Delta z = 0) = [\overline{T}(k_x, f, \Delta z = 0)]$, restricted to a de-scan depth $\Delta z = 0$ and displayed in Fig. 5d of Ref.¹⁷.

To retrieve the axial propagation of this wave-front along the z-axis, a hint is to extract the weighting vector that shall be applied to the set of virtual sources $r_{in}$ to synthesize the coherent guide star. This vector can be retrieved by considering the product between the reduced distortion matrix and the corresponding transmittance such that:

$$V = \overline{\overline{T}}^{\dagger}(\Delta z = 0) \times D_{\{k,f\},r}(\Delta z = 0) \quad (35)$$
<figure><img src="image_12.png" /><figcaption>FIG. 11. Propagation movie associated with the coherent wave packet extracted by iterative phase reversal from a set of virtual sources in speckle in presence of reverberations. a The considered virtual sources $\mathbf{r}_{in} = (x_{in}, z_{in})$ and their neighboring receivers are considered in the orange and blue rectangles on the ultrasound image. b-f The iterative phase reversal wave-front $T(\Delta x, \Delta z, \tau)$ yields the propagation movie of the coherent wave-packet shown here in absolute value at different lapse times $\tau$. The propagation of the coherent wave packet is also displayed in Supplementary Movie 6.</figcaption></figure>

The phase conjugate of the resulting vector $\mathbf{V} = [V(\mathbf{r}_{in})]$ contains the weighting coefficients that shall be applied to the matrix $\mathbf{D}_{\{k,f\},\mathbf{r}}$ in order to synthesize the coherent guide star and retrieve an aberration transmittance for each de-scan depth:

$$ \overline{\overline{\mathbf{T}}}(\Delta z) = \mathbf{D}_{\{k,f\},\mathbf{r}}(\Delta z) \times \mathbf{V}^{\dagger}. \qquad (36) $$

The spatio-temporal evolution of the corresponding wave-front can be deduced by a two-dimensional inverse Fourier transform (Eq. 31). The result is displayed at different time lapses in Fig.11b. As the wave-front induced by a bright scatterer (Fig. 5), the synthetic guide star reveals a tail of multiple reflections generated by the acoustic lens and exhibiting a time period of 2 µs. The time arrivals and focusing depths of the detected echoes, however, differ with the wave-fronts generated by a bright scatterer.
## F. Analytical expression of the coherent wave-front in presence of reverberations

To explain this difference, a theoretical expression of the estimated wave-front can be derived. In speckle and under a single scattering assumption, each reverberated component of the focused reflection matrix can be actually be expressed as follows:

$$ \bar{R}^{(l,m)}(x_{\text{out}}, x_{\text{in}}, z_{\text{out}}, z_{\text{in}}, f) = \iint dzdx H_{\text{out}}^{(m)}(x, z, x_{\text{out}}, z_{\text{out}}, f)\gamma(x, z)H_{\text{in}}^{(l)}(x, z, x_{\text{in}}, z_{\text{in}}, f). \quad (37) $$

Injecting the expression of the reverberated PSF (Eq. 12) leads to the following expression:

$$ \bar{R}^{(l,m)}(x_{\text{out}}, x_{\text{in}}, z_{\text{out}}, z_{\text{in}}, f) = B(f) \iint dzdx H_{\text{out}}^{(0)}(x - x_{\text{out}}, z - z_{\text{out}}^{(m)\star})/c_0, f_c) \gamma(x, z) \\ \times H_{\text{in}}^{(0)}(x - x_{\text{in}}, z - z_{\text{in}}^{(l)\star}, f_c) e^{i2\pi f \left( \frac{2(z-L-z_{\text{in}})-\Delta z}{c} + \frac{2(l+m+1)L}{c_R} \right)}. \quad (38) $$

As before, the virtual source is synthesized at the ballistic focusing point $\mathbf{r}_{\text{in}}^{\star} = (x_{\text{in}}, z_{\text{in}}^{(0)\star})$, with $z_{\text{in/out}}^{(0)\star} = z_{\text{in/out}} + L(1 - \frac{c_R}{c})$ (Eq. 13). In first approximation, the guide star position can be taken into account by assuming a point-like input focused beam at this point (Eq. 23):

$$ \bar{R}^{(l,m)}(\{\Delta x, \Delta z, f\}, \{x_{\text{in}}, z_{\text{in}}\}) \sim \underbrace{\gamma(x_{\text{in}}, z_{\text{in}}^{(0)\star})}_{V^*(x_{\text{in}}, z_{\text{in}})} \times \underbrace{\overline{B(f) H_{\text{out}}^{(0)}(-\Delta x, \Delta z_m^* + \Delta z, f_c)} e^{i2\pi f \left[ \frac{2L}{c} \left( \frac{c}{c_R} - \frac{c_R}{c} \right) + \frac{2(l+m)L}{c_R} - \frac{\Delta z}{c} \right]}}_{T^{(l,m)}(\Delta x, \Delta z, f)} \quad (39) $$

with

$$ \Delta z_m^* = -2mL \frac{c_R}{c}, \quad (40) $$

the de-scan focusing depth of the $m^{th}$ reverberation order. Each component of the reverberated wave-front focuses at a de-scan depth $\Delta z_m^*$ (Eq. 40) that only depends on the reverberation order $m$ of the output path. The comparison between Eqs. 15 and 40 show the difference of focusing de-scan depths for an isolated scatterer and in speckle. Contrary to the wave-field produced by a bright target whose ballistic echo shows a defocus, we retrieve the fact that, in the speckle regime, the ballistic wave-front should focus at $\Delta z = 0$ (Fig. 11d). This is a manifestation of a virtual guide star synthesized, not at the expected depth $z_{\text{in}}$, but in the focusing plane $z_{\text{in}}^*$. Interestingly, each multiply-reflected component shows a defocus proportional to the number of reflections undergone by the output path, an effect also visible for the second and third order multiple reflection events in Figs. 11e and f, respectively.
An inverse Fourier transform of the reverberation transmittance $T(\Delta x, \Delta z, f)$ in Eq. 39 provides an analytical expression for each multiply-reflected component of the coherent wave-packet extracted by the IPR process in the time domain:

$$T^{(l,m)}(\Delta x, \Delta z, \tau) = H_{\text{out}}^{0}(-\Delta x, \Delta z_{m}^{*} - \Delta z, f_{c}) B\left(\tau + \frac{\Delta z}{c} - \tau_{l,m}^{*}\right), \quad (41)$$

with

$$\tau_{l,m}^{*} = 2\tau_{0}^{*} + \frac{2(l+m)L}{c_{R}}, \quad (42)$$

the time shift exhibited by each multiply-reflected echo. The first term of the latter equation corresponds to the time delay experienced by the ballistic wave while the second term is the time delay accumulated by each multiply-reflected wave inside it. Equation 41 shows that, in a speckle environment, the IPR process provides a coherent wave-packet whose spatial envelope accounts for the output PSF at the central frequency and a time shift $\tau_{l,m}^{*}$ that accounts for the delay experienced by each reverberated wave-front on the way-and-return path. Equation 42 confirms theoretically that the ballistic time delay observed with a virtual guide star in speckle ($\tau_{0,0}^{*} \sim -12$ µs in Fig. 11d) is roughly two times the time delay obtained for a bright scatterer ($\tau_{0}^{*} \sim -5.5$ µs Fig. 5d). This difference is explained by the fact that, while the bright scatterer remains defocus in the latter case, the virtual guide star is synthesized in the real focal plane in the former case. This is an important feature since the size of the guide star has also an impact on the quality of the estimated transmittance.

The speckle regime thus offers an ideal framework for aberration and reverberation compensation. In the third paper¹⁷ of the series, we will show how time matrix imaging can be exploited to tailor complex spatio-temporal focusing laws that enable a fine monitoring of the interference between the ballistic wave and the multiply-reflected paths.

## IV. DISCUSSION

In this paper, we introduced the time-focused reflection matrix. Through the various examples discussed, this matrix is shown to provide an access at each point in the medium to a spatio-temporal self-portrait of the focusing process. This allows us to locally probe focusing defects due to wave velocity variations and the temporal spreading of echoes caused by multiple reflections, frequency dispersion, or multiple scattering. Besides direct focusing and imaging applications, the self-portrait of the coherent wave also contains a wealth of
information to perform a better characterization of the medium. The arrival time and focusing
depth of the echoes can be exploited for determining the thickness and speed-of-sound of the
different layers on top of the imaging volume. At a longer term, this observable could also be
included in the cost function of an optimization scheme that would consist in retrieving the
speed-of-sound distribution inside the medium by minimizing the axial distortions exhibited
by the coherent wave packet. In presence of reverberations, this would, however, require a
forward model accounting for multiple scattering phenomena.

In that respect, a more direct route to speed-of-sound tomography will be investigated
in the second paper¹⁶ of the series. Instead of probing the axial distortions of the focused
wave-field in the time domain, the idea will be to investigate the self-portrait at the expected
ballistic time as a function of the speed-of-sound model. Such an approach will allow us to
circumvent the requirement of an homogeneous region of speckle to get the self-portrait of the
wave-field. This condition can actually be quite restricting in view of applications. Moreover,
we will still be able to probe coherent phenomena in speckle such as the Gouy phase in order
to compensate the axial distortions of images induced by wave velocity fluctuations.

One remaining question is the possibility of performing a self-portrait of the wave-field
inside more complex situations than the academic experiments of the current paper. Here,
we considered experimental configurations with a lateral invariance that guaranteed a long
range shift-shift memory effect⁷ and our ability to synthesize a virtual guide star from a
set of incoherent virtual sources. In more realistic situations, this isoplanatic range can
become very reduced and the convergence of the IPR process operated from the k-space
would no longer be guaranteed. In such configurations, a solution would be to operate the
IPR process from a correction basis adapted to the shape and location of the reverberating
layers inside the medium. As shown by adaptive focusing, a correction plane conjugated
with the aberrating layer actually maximizes the isoplanicity of the correction.

To conclude, this work demonstrates the ability of a spatio-temporal control of the wave-
field inside complex media. The originality of this study lies in the fact that this control
is done from the only knowledge of the reflection matrix, without any detector inside the
medium. Using the random distribution of scatterers as virtual microphones, we are able
to probe the spatio-temporal evolution of the wave inside the medium. As we will show in
follow-up studies¹⁶,¹⁷, this observable can be leveraged for a spatio-temporal compensation
of aberrations and reverberations that usually induce strong distortions. Although the
experimental proof-of-concept has been provided with ultrasound, the self-portrait concept
can be potentially extended to all fields of wave physics in which a reflection matrix can be
measured. The cellular components such as lipids or mitochondria at a micron-scale<sup>22</sup> or
the complex network of stacked magma lenses lying below a volcano<sup>23</sup> generate a speckle
wave-field that can be exploited for producing a self-portrait of the focusing process for light
or seismic waves.
**Acknowledgments.** The authors are grateful for the funding provided by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement 819261, REMINISCENCE project, AA). W.L. acknowledges financial support from the SuperSonic Imagine company.

**Author contributions.** A.A. and M.F. initiated the project. A.A. supervised the project. F.B. and W.L. performed the experiments. E.G., F.B. and W.L. developed the post-processing tools for the time-focused reflection matrix. E.G. and A.A. performed the theoretical analysis. E.G. and A.A. prepared the manuscript. E.G., F.B., W.L., M.F. and A.A. discussed the results and contributed to finalizing the manuscript.
Appendix A: Theoretical expression for the point spread function in free space

In this Appendix, we want to express analytically the focused wave-field $H(\mathbf{r}_s, \mathbf{r}_{\text{in}}, f)$ at point $\mathbf{r}_s$ frequency $f$ when trying to focus at a point $\mathbf{r}_{\text{in}}$ assuming an homogeneous speed-of-sound $c_0$ inside the medium. By virtue of Rayleigh Sommerfeld integral, this wave field can be expressed as follows:

$$H(\mathbf{r}, \mathbf{r}_{\text{in}}, f) = \sum_{u} G_{0}^{*}(\mathbf{r}_{\text{in}} - \mathbf{u}, f) \partial_{z} G(\mathbf{r} - \mathbf{u}, f) \quad (A1)$$

where $G(\mathbf{u}, \mathbf{r}_{\text{in}}, f)$ is the ballistic Green's function inside the medium that displays a background speed-of-sound $c$. $G_0^*(\mathbf{u}, \mathbf{r}_{\text{in}}, f)$ corresponds to the incident wave-field expressed at the probe. It corresponds to the phase conjugate of the Green's function between each transducer $\mathbf{u}$ and the targeted focusing point $\mathbf{r}_{\text{in}}$. To go further, the 1D Fourier decomposition of the Green's function can be used:

$$G(x, z, f) = \int dk_x \frac{1}{2j} \frac{1}{\sqrt{\left(\frac{2\pi f}{c}\right)^2 - k_x^2}} \exp\left(iz\sqrt{\left(\frac{2\pi f}{c}\right)^2 - k_x^2}\right) \exp(ik_x x) \quad (A2)$$

Injecting the latter equation into Eq. A1 leads to:

$$H(\mathbf{r}, \mathbf{r}_{\text{in}}, f) = \frac{1}{2} \sum_{u} \int dk_{\text{in}} \int dk_x \frac{e^{iz\sqrt{\left(\frac{2\pi f}{c}\right)^2 - k_{\text{in}}^2}} e^{-iz_{\text{in}}\sqrt{\left(\frac{2\pi f}{c_0}\right)^2 - k_x^2}}}{\sqrt{\left(\frac{2\pi f}{c_0}\right)^2 - k_{\text{in}}^2}} e^{ik_x(x-u)} e^{-ik_{\text{in}}(x_{\text{in}}-u)} \quad (A3)$$

The sum over transducers yields, in first approximation, the Dirac distribution $\delta(k_x - k_{\text{in}})$. Equation A3 simplifies into:

$$H(\mathbf{r}, \mathbf{r}_{\text{in}}, f) = \frac{1}{2} \int dk_{\text{in}} \frac{e^{iz\sqrt{\left(\frac{2\pi f}{c}\right)^2 - k_{\text{in}}^2}} e^{-iz_{\text{in}}\sqrt{\left(\frac{2\pi f}{c_0}\right)^2 - k_x^2}}}{\sqrt{\left(\frac{2\pi f}{c_0}\right)^2 - k_{\text{in}}^2}} e^{ik_{\text{in}}(x-x_{\text{in}})} \quad (A4)$$

Using the parabolic approximation, one can simplify the phase terms as follows:

$$H(\mathbf{r}, \mathbf{r}_{\text{in}}, f) = e^{i2\pi f \left(\frac{z}{c} - \frac{z_{\text{in}}}{c_0}\right)} \frac{1}{2} \int dk_{\text{in}} \frac{e^{i \frac{k_{\text{in}}^2}{2\pi f/c} (z - c_0 z_{\text{in}}/c)}}{\sqrt{\left(\frac{2\pi f}{c_0}\right)^2 - k_{\text{in}}^2}} e^{ik_{\text{in}}(x-x_{\text{in}})} \underbrace{H_{\text{in}}^{(0)}(x-x_{\text{in}}, \frac{c_0}{c}(z-z_{\text{in}}^*), f)}_{H_{\text{in}}^{(0)}(x-x_{\text{in}}, \frac{c_0}{c}(z-z_{\text{in}}^*), f)} \quad (A5)$$
Appendix B: Theoretical expression for the point spread function in presence of a reverberating plate

In presence of the reverberating plate, each component of the focused reflection matrix (Eq. 10) can be expressed as follows:

$$ \bar{\mathbf{R}}_{\Delta \mathbf{r}}^{(l,m)}(f) = \gamma_s \bar{B}(f) H_{\text{out}}^{(m)}(x_s, z_s, x_{\text{out}}, z_{\text{out}}, f) H_{\text{in}}^{(l)}(x_{\text{in}}, z_s, x_{\text{in}}, z_{\text{in}}, f). \quad (B1) $$

where $H_{\text{in}}^{(l)}(x_s, z_s, x_{\text{in}}, z_{\text{in}}, f)$ and $H_{\text{out}}^{(m)}(x_s, z_s, x_{\text{out}}, z_{\text{out}}, f)$ stand for the $l^{\text{th}}$ and $m^{\text{th}}$ reverberation order of the input and output point spread functions, respectively. The $l^{\text{th}}$ and $m^{\text{th}}$ reverberation order is associated with a number $2l$ and $2m$ of internal reflections inside the reverberating layer.

Following the same rationale as in Appendix A, each reverberation order of the input point spread function can be expressed as follows:

$$ H_{\text{in}}^{(l)}(\mathbf{r}, \mathbf{r}_{\text{in}}, f) = \frac{1}{2} \int dk_{\text{in}} \frac{e^{i(z-L)\sqrt{\left(\frac{2\pi f}{c}\right)^2 - k_{\text{in}}^2}} e^{i(2l+1)L\sqrt{\left(\frac{2\pi f}{c_R}\right)^2 - k_{\text{in}}^2}} e^{-iz_{\text{in}}\sqrt{\left(\frac{2\pi f}{c_0}\right)^2 - k_{\text{in}}^2}} e^{ik_{\text{in}}(x-x_{\text{in}})} \sqrt{\left(\frac{2\pi f}{c_0}\right)^2 - k_{\text{in}}^2} \quad (B2) $$

Assuming a wave velocity model now equal to the phantom speed-of-sound ($c_0 = c$), the last expression becomes

$$ H_{\text{in}}^{(l)}(\mathbf{r}, \mathbf{r}_{\text{in}}, f) = \frac{1}{2} \int dk_{\text{in}} \frac{e^{i(z-L-z_{\text{in}})\sqrt{\left(\frac{2\pi f}{c}\right)^2 - k_{\text{in}}^2}} e^{i(2l+1)L\sqrt{\left(\frac{2\pi f}{c_R}\right)^2 - k_{\text{in}}^2}} e^{ik_{\text{in}}(x-x_{\text{in}})} \sqrt{\left(\frac{2\pi f}{c}\right)^2 - k_{\text{in}}^2} \quad (B3) $$

Using the parabolic approximation, one can simplify the phase terms as follows:

$$ H_{\text{in}}^{(l)}(\mathbf{r}, \mathbf{r}_{\text{in}}, f) = e^{i2\pi f \left( \frac{z-L-z_{\text{in}}}{c} + \frac{(2l+1)L}{c_R} \right)} \frac{1}{2} \int dk_{\text{in}} \frac{e^{i \frac{k_{\text{in}}^2}{2\pi f/c} (z-L-z_{\text{in}}+c_R(2l+1)L/c)}}{\sqrt{\left(\frac{2\pi f}{c}\right)^2 - k_{\text{in}}^2}} e^{ik_{\text{in}}(x-x_{\text{in}})} \underbrace{H_{\text{in}}^{(0)}(x-x_{\text{in}}, z-z_{\text{in}}^*, f)}_{H_{\text{in}}^{(0)}(x-x_{\text{in}}, z-z_{\text{in}}^*, f)} \quad (B4) $$

with

$$ z_{\text{in}}^{*(l)} = z_{\text{in}} + L - (2l+1)L \frac{c_R}{c}, \quad (B5) $$

the depth of the focusing plane for the $l^{\text{th}}$ reverberation order. As in Appendix A, each reverberation component of the point spread function is thus the product of an axial propagation term which oscillates rapidly with frequency and a transverse focusing term that
evolves slowly with frequency. In first approximation, the latter term can thus be replaced by its value at the central frequency $f_c$:

$$H_{\text{in}}^{(l)}(\mathbf{r}, \mathbf{r}_{\text{in}}, f) = \underbrace{e^{i2\pi f \left(\frac{z-L-z_{\text{in}}}{c} + \frac{(2l+1)L}{c_R}\right)}}_{\text{axial propagation}} \underbrace{H_{\text{in}}^{(0)}(x-x_{\text{in}}, z-z_{\text{in}}^{\star(l)}, f_c)}_{\text{transverse focusing}} \quad (\text{B6})$$

Note that an equivalent expression can be obtained for the reverberation components of the output PSF by replacing the index $l$ by index $m$:

$$H_{\text{out}}^{(m)}(\mathbf{r}, \mathbf{r}_{\text{out}}, f) = \underbrace{e^{i2\pi f \left(\frac{z-L-z_{\text{out}}}{c} + \frac{(2m+1)L}{c_R}\right)}}_{\text{axial propagation}} \underbrace{H_{\text{out}}^{(0)}(x-x_{\text{out}}, z-z_{\text{out}}^{\star(m)}, f_c)}_{\text{transverse focusing}} \quad (\text{B7})$$

with

$$z_{\text{out}}^{\star(m)} = z_{\text{out}} + L - (2m + 1)L \frac{c_R}{c}. \quad (\text{B8})$$

Injecting Eqs. B6 and B7 into Eq. B1 leads to the following equation for each reflection matrix component, $\bar{R}^{(l,m)}(x_s, \Delta z, z_{\text{in}}, f)$:

$$\bar{R}^{(l,m)}(\Delta x, \Delta z, x_{\text{in}}, z_{\text{in}}, f) = \gamma_s \bar{B}(f) e^{i2\pi f \left( \frac{2(z-L-z_{\text{in}})-\Delta z}{c} + \frac{(2l+2m+2)L}{c_R} \right)} \times H_{\text{out}}^{(0)}(\Delta x, z_s - z_{\text{in}}^{\star(m)} - \Delta z, f) H_{\text{in}}^{(0)}(\Delta x, z_s - z_{\text{in}}^{\star(l)}, f). \quad (\text{B9})$$

The position of the scatterer $z_s$ is controlled by the time-of-flight expected for a scatterer at the targeted depth $z_{\text{in}}$ under a wave velocity model $c$:

$$z_s = z_{\text{in}} + L \left(1 - \frac{c}{c_R}\right) \quad (\text{B10})$$

Injecting this expression of $z_s$ into Eq. B9 and using the expressions of the focusing depths leads to the following expression of the focused $\mathbf{R}$-matrix coefficients in de-scanned coordinates:

$$R^{(l,m)}(\Delta x, \Delta z, x_{\text{in}}, z_{\text{in}}, \tau) \propto \underbrace{H_{\text{in}}^{(0)}(0, \Delta z_l^*, f_c)}_{\text{constant}} \underbrace{H_{\text{out}}^{(0)}(\Delta x, \Delta z_m^* - \Delta z, f_c)}_{\text{transverse focusing}} \underbrace{\bar{B}(f)e^{-i2\pi f \left(\frac{\Delta z}{c} - \frac{2(l+m)L}{c_R}\right)}}_{\text{axial propagation}} \quad (\text{B11})$$

with

$$\Delta z_n^* = 2nL \frac{c_R}{c} + L \left( \frac{c_R}{c} - \frac{c}{c_R} \right), \quad (\text{B12})$$

the focusing depth of the $n^{\text{th}}$ reverberation order in the de-scanned frame. An inverse Fourier transform of Eq. B11 leads to final expression of the focusing reflection matrix coefficients in the time domain for an isolated bright scatterer at targeted position $(x_{\text{in}}, z_{\text{in}})$:

$$R^{(l,m)}(\Delta x, \Delta z, x_{\text{in}}, z_{\text{in}}, \tau) \propto \underbrace{H_{\text{in}}^{(0)}(0, \Delta z_l^*, f_c)}_{\text{constant}} \underbrace{H_{\text{out}}^{(0)}(\Delta x, \Delta z_m^* - \Delta z, f_c)}_{\text{spatial focusing}} \underbrace{B \left( \tau + \frac{\Delta z}{c} - \frac{2(l+m)L}{c_R} \right)}_{\text{time focusing}} \quad (\text{B13})$$
The ballistic component ($l = m = 0$) is expected to focus at a de-scan focusing depth
$\Delta z_0^* = L \left(\frac{c_r}{c} - \frac{c}{c_r}\right)$ and focusing time $\tau_0^* = -\Delta z_0^*/c$.

References.

[1] J. W. Goodman, *Speckle Phenomena in Optics: Theory and Applications* (Roberts and Company, Englewood, Colorado, USA, 2007).

[2] R. Wagner, S. Smith, J. Sandrik, and H. Lopez, Statistics of speckle in ultrasound b-scans, IEEE Trans. Sonics Ultras. **30**, 156 (1983).

[3] A. P. Mosk, A. Lagendijk, G. Lerosey, and M. Fink, Controlling waves in space and time for imaging and focusing in complex media, Nat. Photonics **6**, 283 (2012).

[4] S. Feng, C. Kane, P. A. Lee, and A. D. Stone, Correlations and Fluctuations of Coherent Wave Transmission through Disordered Media, Phys. Rev. Lett. **61**, 834 (1988).

[5] I. Freund, M. Rosenbluh, and S. Feng, Memory Effects in Propagation of Optical Waves through Disordered Media, Phys. Rev. Lett. **61**, 2328 (1988).

[6] S. Shahjahan, A. Aubry, F. Rupin, B. Chassignole, and A. Derode, A random matrix approach to detect defects in a strongly scattering polycrystal: How the memory effect can help overcome multiple scattering, Appl. Phys. Lett. **104**, 234105 (2014).

[7] G. Osnabrugge, R. Horstmeyer, I. N. Papadopoulos, B. Judkewitz, and I. M. Vellekoop, Generalized optical memory effect, Optica **4**, 886 (2017).

[8] J. Bertolotti, E. G. van Putten, C. Blum, A. Lagendijk, W. L. Vos, and A. P. Mosk, Non-invasive imaging through opaque scattering layers, Nature **491**, 232 (2012).

[9] O. Katz, P. Heidmann, M. Fink, and S. Gigan, Non-invasive single-shot imaging through scattering layers and around corners via speckle correlations, Nat. Photon. **8**, 784 (2014), number: 10 Publisher: Nature Publishing Group.

[10] S. Kang, P. Kang, S. Jeong, Y. Kwon, T. D. Yang, J. H. Hong, M. Kim, K. Song, J. H. Park, J. H. Lee, M. J. Kim, K. H. Kim, and W. Choi, High-resolution adaptive optical imaging within thick scattering media using closed-loop accumulation of single scattering, Nat. Commun. **8**, 2157 (2017).

[11] W. Lambert, L. A. Cobus, T. Frappart, M. Fink, and A. Aubry, Distortion matrix approach for ultrasound imaging of random scattering media, Proc. Nat. Acad. Sci. USA **117**, 14645 (2020).

[12] W. Lambert, L. A. Cobus, M. Couade, M. Fink, and A. Aubry, Reflection matrix approach for
quantitative imaging of scattering media, Phys. Rev. X **10**, 021048 (2020).

[13] A. Velichko, Quantification of the Effect of Multiple Scattering on Array Imaging Performance, IEEE Trans. Ultras. Ferroel. Freq. Cont. **67**, 92 (2020).

[14] W. Lambert, J. Robin, L. A. Cobus, M. Fink, and A. Aubry, Ultrasound matrix imaging—Part I: The focused reflection matrix, the f-factor and the role of multiple scattering, IEEE Trans. Med. Imag. **41**, 3907 (2022).

[15] F. Bureau, J. Robin, A. Le Ber, W. Lambert, M. Fink, and A. Aubry, Three-dimensional ultrasound matrix imaging, Nat. Commun. **14**, 6793 (2023).

[16] F. Bureau, E. Brenner, E. Korta Martiartu, Naiara Giraudat, A. Le Ber, W. Lambert, L. Carmier, A. Guibal, M. Fink, and A. Aubry, Self-portrait of the focusing process in speckle: II. Gouy phase shift for defocus correction and pixel depth reassignment, arXiv:2409.13901 (2026).

[17] E. Giraudat, F. Bureau, W. Lambert, M. Fink, and A. Aubry, Self-portrait of the focusing process in speckle: III. Tailoring complex spatio-temporal focusing laws to overcome reverberations in reflection imaging, arXiv:2602.05908 (2026).

[18] J.-L. Robert and M. Fink, Green's function estimation in speckle using the decomposition of the time reversal operator: Application to aberration correction in medical imaging, J. Acoust. Soc. Am. **123**, 866 (2008).

[19] H. Bendjador, T. Deffieux, and M. Tanter, The SVD Beamformer: Physical Principles and Application to Ultrafast Adaptive Ultrasound, IEEE Trans. Med. Imag. **39**, 3100 (2020).

[20] W. Lambert, L. A. Cobus, J. Robin, M. Fink, and A. Aubry, Ultrasound matrix imaging—Part II: The distortion matrix for aberration correction over multiple isoplanatic patches, IEEE Trans. Med. Imag. **41**, 3921 (2022).

[21] R. Mallart and M. Fink, Adaptive focusing in scattering media through sound-speed inhomogeneities: The van Cittert Zernike approach and focusing criterion, J. Acoust. Soc. Am. **96**, 3721 (1994).

[22] V. Barolle, F. Bureau, N. Guigui, P. Balondrade, V. Brochard, O. Dubois, A. Jouneau, A. Bonnet-Garnier, and A. Aubry, Optical matrix imaging applied to embryology (2024), arXiv:2410.11126.

[23] E. Giraudat, A. Burtin, A. L. Ber, M. Fink, J.-C. Komorowski, and A. Aubry, Matrix imaging as a tool for high-resolution monitoring of deep volcanic plumbing systems with seismic noise,