# Detection and characterization of targets in complex media using fingerprint matrices

Arthur Le Ber, Antton Goïcoechea, Lukas M Rachbauer, William Lambert, Xiaoping Jia, Mathias Fink, Arnaud Tourin, Stefan Rotter, Alexandre Aubry

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

## To cite this version:

Arthur Le Ber, Antton Goïcoechea, Lukas M Rachbauer, William Lambert, Xiaoping Jia, et al. Detection and characterization of targets in complex media using fingerprint matrices. Nature Physics, 2025, 10.1038/s41567-025-03016-2. http://dx.doi.org/10.1038/s41567-025-03016-2

HAL Id: hal-04878230

https://hal.science/hal-04878230v3

Submitted on 3 Oct 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Detection and characterization of targets in complex media using
fingerprint matrices

Arthur Le Ber,¹, * Antton Goïcoechea,¹, * Lukas M. Rachbauer,² William Lambert,³

Xiaoping Jia,¹ Mathias Fink,¹ Arnaud Tourin,¹ Stefan Rotter,² and Alexandre Aubry¹, †

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

²Institute for Theoretical Physics, Vienna University
of Technology (TU Wien), Vienna, Austria

³Supersonic Imagine, Aix-en-Provence, France

(Dated: October 3, 2025)

Abstract

As waves propagate through a complex medium, they undergo multiple scattering events. This phenomenon is detrimental to imaging, as it causes full blurring of the image. Here, we report a method to detect, localize and characterize any scattering target embedded in a complex medium. We introduce a fingerprint operator that contains the specific signature of the target with respect to its environment. When applied to the recorded reflection matrix, it provides a likelihood index of the target state. This state can be the target's position for localization purposes, its shape for characterization or any other parameter that influences the target response. We demonstrate the versatility of our method by performing proof-of-concept ultrasound experiments on elastic spheres buried inside a strongly scattering granular suspension and on lesion markers, commonly used to monitor breast tumours, embedded in a foam mimicking soft tissue. Furthermore, we show how the fingerprint operator can be leveraged for characterizing the complex medium itself by mapping the fiber architecture within muscle tissue. Our method is broadly applicable to different types of waves beyond ultrasound for which multi-element technology allows a reflection matrix to be measured.

<footnote>* These authors equally contributed to this work</footnote>

<footnote>† Corresponding author: alexandre.aubry@espci.fr</footnote>
Introduction

With the emergence of multi-element technology, it has become possible to cope with the
enormous complexity of disordered media for focusing waves through them or for imaging
objects hidden behind them. This feat has been demonstrated in acoustics through time-
reversal mirrors¹ and in optics using wavefront shaping techniques². A matrix formalism has
become particularly useful in this context, especially when wave-fields are controlled by arrays
of independent elements both in transmission³,⁴ and reflection⁵,⁶. Since an inhomogeneous
medium can be treated as one realization of a random process, some aspects of random
matrix theory and basic concepts from scattering theory have been fruitfully applied to wave
control or imaging through complex media⁷⁻⁹. In particular, transmission eigenchannels have
been shown to provide a properly designed combination of incident waves that can be fully
transmitted through a disordered medium¹⁰⁻¹². The reflection eigenchannels, on the other
hand, have been shown to focus on bright point-like targets embedded in a scattering medium,
provided that the disorder is not too strong¹³,¹⁴. For more complex targets, these approaches
are not adapted since the one-to-one association between each eigenstate and each target
is no longer verified in the general case¹⁵,¹⁶. Much progress has also been made using the
reflection matrix¹⁷⁻²² for imaging, in particular, to compensate the effects of disorder²³⁻²⁷.
All currently available techniques, however, face the key challenge that, in the regime of
strong multiple scattering, such a compensation becomes close to impossible.

In the present article, we address this outstanding challenge with an approach that is not
based on compensating the effect of disorder, but rather by detecting correlations in the
scattering of waves that survive even for very thick disordered media. The starting point for
this approach is the following scattering invariant operator that has been constructed based
on the matrix product T × T₀† between the transmission matrix of a disordered medium T
and a matrix T₀ of a homogeneous reference medium²⁸. The associated eigenstates of this
operator have been shown to be ‘scattering invariant modes’ in the sense that they exhibit
a transmitted field pattern behind a disordered medium that is identical to that of purely
ballistic waves, independent of the multiple scattering events endured by the waves inside
the medium.

In the context of target detection, we are looking for input states that are reflected from
a target embedded inside a disordered medium in the same way as from a target embedded
in free (or homogeneous) space. Correspondingly, the operator Γ that we introduce here is
defined as the matrix product $\Gamma = \mathbf{R} \times \mathbf{R}_0^\dagger$ between a measured reflection matrix $\mathbf{R}$ of a disordered medium potentially hiding a target and the conjugate transpose of a reference matrix $\mathbf{R}_0$ associated with a target in absence, this time, of any disordered environment. This fingerprint matrix $\mathbf{R}_0$ can be parametrized by a target state vector $\mathbf{q}$ whose coefficients account for the position, size and shape of the target alone. The more closely the position, size or shape of the target in the fingerprint matrix $\mathbf{R}_0(\mathbf{q})$ match their true values realized in the experimentally measured matrix $\mathbf{R}$, the higher the correlations between these two matrices will be - allowing us to determine the target properties through an optimization problem. Our concept is therefore extremely versatile since the fingerprint matrix can be designed as a function of the experimental situation, and can be applied to any kind of waves, provided that a measurement of the reflection matrix is possible.

To demonstrate the capabilities and versatility of the fingerprint operator, we conduct three ultrasound imaging experiments. First, we detect and localize spherical metal targets buried within a highly scattering granular medium. Second, we demonstrate the detection of a clinical lesion marker embedded in a scattering tissue-mimicking phantom. Finally, we show that the fingerprint operator enables the extraction of structural information such as orientation and size of anisotropic features in biological tissue, exemplified by the in vivo imaging of muscle fiber architecture. These three examples highlight the robustness and flexibility of the approach under various scattering conditions and target complexities.

# Results

## Multiple scattering as a nightmare for imaging

The first experiment consists in the ultrasound imaging of two metal spheres embedded into a granular suspension consisting of randomly packed glass beads (diameter $\emptyset \approx 300-315 \ \mu\text{m}$) immersed in water with a packing density of 60%$^{29}$ (Fig. 1a). The metal spheres have diameters $d_1 = 10 \ \text{mm}$ and $d_2 = 8 \ \text{mm}$. Their centers $\mathbf{r}_1$ and $\mathbf{r}_2$ are positioned 9 and 7 mm below this granular suspension surface, respectively. The experimental procedure first consists in the acquisition of the reflection matrix using a two-dimensional array of 1024 transducers placed on top of the granular suspension surface into which the metal spheres are embedded (Methods). The reflection matrix is here acquired using a set of plane waves$^{30}$ (Fig. 1b, Methods) over a 1.8-2.6 MHz frequency bandwidth (Extended Data Tab. 1). For each plane wave with incidence angles $\theta_{\text{in}} = (\theta_x, \theta_y)$, the time-dependent reflected wave field
$R(\mathbf{u}_{\text{out}}, \boldsymbol{\theta}_{\text{in}}, t)$ is recorded by each transducer $\mathbf{u}_{\text{out}}$ (Fig. 1c). This set of wave-fields forms a reflection matrix acquired in the plane wave basis, $\mathbf{R}_{\boldsymbol{\theta}\mathbf{u}}(t) = [R(\mathbf{u}_{\text{out}}, \boldsymbol{\theta}_{\text{in}}, t)]$. A Fourier transform can then be applied to the experimentally acquired reflection matrix to obtain its frequency-dependent counterpart, $\mathbf{R}_{\boldsymbol{\theta}\mathbf{u}}(f) = \int dt \mathbf{R}_{\boldsymbol{\theta}\mathbf{u}}(t) \exp(-j2\pi ft)$, with $f$ the frequency.

<figure><img src="image_3.png" /><figcaption>FIG. 1. Multiple scattering fog. a, Experimental configuration: A 2D array of 1024 ultrasound transducers acting both as transmitters and receivers is used to detect two metal spheres that are embedded in strongly scattering glass beads below water. b, The reflection matrix is measured by insonifying the medium with a set of incident plane waves ($\boldsymbol{\theta}_{\text{in}}$). c, Each back-scattered wave field $R(\mathbf{u}_{\text{out}}, \boldsymbol{\theta}_{\text{in}}, t)$ is recorded on each transducer $\mathbf{u}_{\text{out}}$ of the same probe. d, Volumetric visualization of the probed medium by projecting the intensity maximum of the confocal image obtained by a delay-and-sum beamforming applied to the recorded reflection matrix (a dynamic view of this is shown in Supplementary Movie 1). With this conventional imaging technique, the specular echo of the interface between water and the glass bead suspension largely predominates and the targets are completely invisible on the image.</figcaption></figure>

A focused beamforming algorithm (Methods) is then applied to the reflection matrix to obtain a confocal image of the medium under investigation. The result is displayed in Fig. 1d. Although the image shows a strong specular echo at the surface of the granular suspension, the two spheres located beneath this surface are fully invisible on the image. Two reasons can be invoked to explain the failure of a standard imaging process in this experimental
configuration:

(i) First, the granular suspension is strongly scattering, as evidenced by the scattering mean free path $\ell_s \sim 1.8$ mm exhibited by ultrasonic waves in the frequency range of interest (Supplementary Section S1). Note that scattering by the glass bead packing is isotropic in the frequency range under study. The transport mean free path is therefore equivalent to $\ell_s$. The imaging depth $z$ of each sphere ranges from 4 to $5\ell_s$. At such depths, the direct echoes of the target spheres are hidden by a strongly predominant multiple scattering background, which generates an intense speckle fog in Fig. 1d. The image contrast $C_I$, which scales as the single scattering rate, undergoes a drastic exponential decay with the penetration depth $z$ $^{31}$:

$$C_I \propto \exp(-z/\ell_s) \qquad (1)$$

The target echo is thus decreased by a factor ranging from 50 ($z \sim 4\ell_s$) to 150 ($z \sim 5\ell_s$) with respect to the multiple scattering background in the present experiment.

(ii) Second, a confocal imaging scheme is far from being optimal for the imaging of objects such as the elastic spheres considered here. Each of them actually supports a large number of resonating modes or/and internal waves giving rise to a complex spatio-temporal echo.

## Leveraging correlations through the fingerprint operator

Even if the response of an elastic target may be spectrally and spatially quite complex$^{15,16,32-34}$, all these features can be adequately captured by measuring the target's reflection matrix $R_0$ in free space (Methods, Extended Data Fig. 1) and leveraged by the scattering-invariant operator $\Gamma = R \times R_0^\dagger$.

Rather than inspecting individual eigenstates and eigenvalues of $\Gamma$, however, we quantify the correlation signatures through the trace of this matrix, integrated over the entire recorded frequency interval and properly normalized:

$$\gamma(\mathbf{q}) = \frac{|\int df \mathrm{Tr} \{\mathbf{R}(f) \times \mathbf{R}_0^\dagger(\mathbf{q}, f)\}|}{[\int df ||\mathbf{R}(f)||_F \times \int df ||\mathbf{R}_0(\mathbf{q}, f)||_F]^{1/2}}, \qquad (2)$$

where the symbol $|| \cdot ||_F$ stands for the Frobenius norm. The likelihood index $\gamma(\mathbf{q})$ for the target to be in state $\mathbf{q}$, involves an emulated target fingerprint matrix $R_0(\mathbf{q})$, where $\mathbf{q}$
accounts for the target state such as, for instance, the target's position **r**, size, shape or any other parameter that has an influence on the target's response.

The question relevant to the practical implementation of this approach is now the construction of the fingerprint operator. The most direct way of doing it is to consider the free space reflection matrix $\mathbf{R}_0(\mathbf{r}_0)$ shown in Extended Data Fig. 1 and to emulate a fingerprint matrix $\mathbf{R}_0(\mathbf{r})$ by a virtual shift of the target from the initial position $\mathbf{r}_0$ to any point $\mathbf{r}$ through simple matrix operations (Methods). A first validation of Eq. 2 can be performed by applying it to the free-space reflection matrix such that $\mathbf{R} = \mathbf{R}_0(\mathbf{r}_0)$. The result is displayed in Extended Data Fig. 1f. The comparison with the confocal image displayed in Extended Data Fig. 1e highlights the benefit of the proposed approach. All the energy radiated by the elastic sphere that was initially dispersed in time and space on the diagonal and off-diagonal coefficients of the focused $\mathbf{R}$-matrix is now concentrated onto a single pixel of the image.

The gain $G$ in contrast compared to the confocal image is drastic since it scales as the product between the numbers of spatial and temporal d.o.f, $N_S$ and $N_T$, exhibited by the target (Supplementary Section S4):

$$G \sim N_S \times N_T \qquad (3)$$

On the one hand, $N_T$ can be expressed as the product of the frequency bandwidth $\Delta f$ and the reverberation time $\Delta t$ of the target echo:

$$N_T \sim \Delta f \Delta t. \qquad (4)$$

In the present case, the reverberation time $\Delta t$ lasts 40 µs for each sphere which yields $N_t \sim 32$. On the other hand, $N_S$ is the effective rank of the target reflection matrix $\mathbf{R}_0$, that can be assessed from the singular value distribution of $\mathbf{R}_0$ (Eq. 17, Supplementary Figure S8): $N_S \sim 20$ for sphere 1 and $N_S \sim 12$ for sphere 2. This larger value for sphere 1 is explained by an effective rank of $\mathbf{R}_0$ scaling as the number of lateral resolution cells covered by the target$^{15,16}$:

$$N_S \sim \mathcal{A}/\delta\rho_D^2 \qquad (5)$$

with $\mathcal{A}$, the target physical cross-section and $\delta\rho_D \sim \lambda/(2 \sin \alpha)$, the diffraction-limited resolution length of the target image in Extended Data Fig. 1, and $\alpha$ the angle under which the target is seen by the array. The resolution length actually determines the transverse size of the focal spot displayed by the $\gamma$-map in Extended Data Fig. 1.
In the present case, given the values of $N_S$ and $N_T$ provided above, the fingerprint operator can, in principle, increase the sphere contrast by a factor $G$ (Eq. 3) ranging from 380 (sphere 2) to 640 (sphere 1). This contrast enhancement is significantly larger than the decay factor experienced by each sphere echo with respect to multiple scattering: 50 for sphere 1 and 150 for sphere 2 (Eq. 1). The fingerprint operator can therefore be expected to detect and locate the spheres in the experiment shown in Fig. 1.

## Seeing into the scattering fog

As we will show in the following, the target's signatures can be made even more specific with respect to its environment by considering a generalized scattering invariant operator, $\hat{\Gamma} = R \times F^{\dagger}$, involving a more selective fingerprint matrix $F$ (Supplementary Section S2). To that aim, different filters are applied to the reference matrix $R_0$ to build the matrix $F$ (Methods) and make the response of the target more specific with respect to its environment. This optimized fingerprint operator $F$ is then used to compute the likelihood maps of the 8 mm- and 10 mm-diameter spheres in the experimental configuration. To highlight the drastic improvement provided by our approach, we superimpose the corresponding $\gamma$-maps encoded in green and red, respectively, onto the confocal image in Fig. 2, where different cross-sections of the imaged volume are displayed: While the confocal image is fully blurred by the multiple scattering fog, each likelihood map allows the unambiguous detection of each sphere.

The scattering invariant operator tends to filter out the diffuse waves and select the paths that are insensitive to scattering, i.e. paths that behave as if we were in free space²⁸ (among which the coherent ballistic paths). Even though these paths correspond to an extremely weak contribution to the total wave-field in the diffusive regime, the scattering-invariant operator allows us to leverage them for target detection. This is made possible by the complexity of the target signature: The weight of scattering-invariant paths is enhanced by the number of spatio-temporal degrees-of-freedom exhibited by the target echo (Eq. 3). Here, the strong contrast enhancement provided by the fingerprint operator partially relies on the strong impedance mismatch between the metal spheres and the granular suspension. Weaker enhancements are expected for softer materials.

In Supplementary Movie 2, we show how the scattering invariant operator can also work in a dynamic scenario where the target (sphere 1) sinks inside the vibrated dense granular suspension that then acts as quicksand. This experiment was already described in a previous
<figure><img src="image_4.png" /><figcaption>FIG. 2. Detecting and localizing target metal spheres hidden in the multiple scattering fog. Different cross-sections of the $\gamma$-map for each sphere (diameter $d_1 = 10$ mm in red, $d_2 = 8$ mm in green) are superimposed to the corresponding confocal image (B&W scale in dB). **a** ($x, z$)—cross-section at $y = y_1 = y_2$ (scale bar: 10 mm). **b** ($y, z$)—cross-section at $x = x_1$. **c** ($y, z$)—cross-section at $x = x_2$. The confocal images correspond to cross-sections of Fig. 1d, the dashed circles correspond to the known radius of each sphere. The value of the contrast $C_\gamma$ at the target's position with respect to the average value of $\gamma$ outside each target is also indicated. While the targets are completely invisible in the confocal images, they appear very well localized through the calculation of their associated likelihood index $\gamma$.</figcaption></figure>

work²⁹ for a granular suspension of smaller glass beads (100 µm-diameter) that was almost transparent for ultrasound. Here, thanks to the fingerprint matrix concept, we are able to track the trajectory of the sphere in a much stronger scattering regime.

## Detection rate and localization precision

Figure 2 demonstrates that the exponential decrease of the target echo (captured by $C_I$ in Eq. 1) is largely compensated by the gain $G$ (Eq. 3) provided by the fingerprint operator. Experimentally, the target contrast $C_\gamma$ can be evaluated by considering the ratio between the value of $\gamma$ at each target position and its average outside of each target. A significant contrast ranging between 6 and 7 is found for each metal sphere (see $C_\gamma$ values indicated for each cross-section in Fig. 2). Those values are in qualitative agreement with the theoretical prediction: $C_\gamma \sim G \times C_I \sim 4.2$ for sphere 1 and 7.6 for sphere 2.
The observed contrast is an important feature since it allows to compute a probability of false alarm PFA, i.e. the probability that the bright spots in the γ-maps of Fig. 2 correspond to a normal multiple scattering fluctuation and not to coherent targets. Assuming Rayleigh statistics for multiple scattering speckle, the observed contrast values ($C_γ \sim 6-7$) indicate an extremely small probability of false alarm: $PFA < 10^{-8}$ (see Supplementary Section S6 and Figure S15). Hence, the detection of metal spheres is achieved with an extremely high degree of confidence in Fig. 2.

Besides providing a detection rate, the contrast $C_γ$ also dictates the precision $\delta\rho_L$ of the localization process predicted by the Cramér-Rao bound$^{35,36}$ (Supplementary Section S5):

$$ \delta\rho_L \sim \frac{\delta\rho_D}{\sqrt{C_γ}} \sim \frac{\delta\rho_D}{\sqrt{GC_I}} \quad (6) $$

By increasing the target contrast with respect to its environment, the scattering invariant operator enhances the localization precision by a factor scaling as $\sqrt{N_T N_S}$ (Eq. 3). The more complex the target echo is, the more beneficial our method becomes.

This contrast-based evaluation of detection and localization performance can also be used to evaluate the robustness of the method with respect to electronic noise that would unavoidably pollute the recorded data in presence of absorption. To do so, the initial contrast $C_I$ can be replaced by a signal-to-noise ratio (SNR). The detection threshold is then inversely proportional to the factor $G$ (Eq. 3) and the localization precision (Eq. 6) scales as $\delta\rho_L \sim \delta\rho_D/\sqrt{G \times \text{SNR}}$.

## Lesion marker detection and localization in ultrasound speckle

To show the generality of the concept and its potential applications, a second experiment depicted in Fig. 3a (Methods, Extended Data Tab. 1) was carried out in a weaker scattering regime ($z \sim l_s$) characteristic of ultrasound imaging and in which single and multiple scattering coexist$^{31}$. It involves the detection of a lesion marker (Tumark®Vision, Somatex, Fig. 3b) generally used to monitor breast tumors in clinical settings$^{37}$. Although these markers have been specifically designed to yield a clear signature on ultrasound images, their detection and localization are often hampered by the ultrasound speckle generated by unresolved scatterers that are randomly distributed in tissue. In the experiment depicted in Fig. 3, ultrasound speckle is generated by a foam soaked in water in which the lesion marker has been embedded. This foam generates an ultrasound speckle similar to that of
soft tissues. The corresponding confocal image highlights this characteristic speckle which here prevents an unambiguous detection and localization of the target (Fig. 3c). On the contrary, the fingerprint operator designed for this lesion marker allows a clear detection and sharp localization of the lesion marker (Fig. 3d). This proof-of-concept shows the potential of the fingerprint operator for the monitoring of any object (e.g. needle, catheter, marker, etc.) used in interventional radiology.

<figure><img src="image_5.png" /><figcaption>FIG. 3. Localizing a lesion marker in ultrasound speckle. a, Experimental configuration: The 32 × 32 probe is used to image a lesion marker embedded into foam soaked in water. The position $\mathbf{r}_m = (x_m, y_m, z_m)$ of the lesion marker center is (3.5, 4.5, 32) mm. b, Photography of the lesion marker (credit: Arthur Le Ber). c, Longitudinal cross-section of the confocal ultrasound image in the plane $y = y_m$ (scale bar: 5 mm). d, Corresponding likelihood map of the lesion marker superimposed onto the confocal image in transparency.</figcaption></figure>

## Towards quantitative imaging of complex media

Beyond the detection and localization of a known intruder inside a complex medium, the fingerprint operator can also be applied to unknown targets. In that case, a dictionary of fingerprint matrices can be analytically or numerically built for targets of arbitrary shape, size, orientation etc. A maximization of the likelihood index with respect to these parameters can then lead to a quantitative characterization of the inspected medium.

As a proof-of-concept, we will show how it can be used for mapping the local anistropy of a fibrous medium. This information is particularly relevant for muscle tissues in ultrasound
imaging in order to diagnose neuromuscular³⁸ or myocardial³⁹ diseases. The experiment is carried out on a human calf, in-vivo, the fibers of which are partially visible on the ultrasound image displayed in Fig. 4b. This image is built from the reflection matrix **R** recorded using a linear phased array of 256 transducers over the [5.5; 9.5] MHz frequency bandwidth (Methods, Extended Data Tab. 1). In order to image fibrous tissues, the fingerprint matrix is defined as a free-space reflection matrix **R**₀(**q**) associated with a reflecting 1D mirror whose state **q** is described by its position **r**, its orientation **α** and its dimension **L** (Fig. 4a). The dictionary of fingerprint matrices **R**₀(**q**) is constructed numerically (Methods).

The scalar product of the recorded reflection matrix **R** with the fingerprint matrix **R**₀ (Eq. 2) provides a likelihood index $γ(\mathbf{r}, \alpha, L)$ with respect to the parameters $\alpha$ and $L$ at each point **r** in the image (Eq. 2). The dependence of the $γ$-map with respect to parameters $\alpha$ and $L$ is shown for two points **r** in the field-of-view in Figs. 4c and d. The maximum value of this quantity gives the orientation $\alpha_{opt}(\mathbf{r})$ and the local correlation length $L_{opt}(\mathbf{r})$ of the fibrous medium at each point **r**:

$$ \{\alpha_{opt}(\mathbf{r}), L_{opt}(\mathbf{r})\} = \arg\max_{\{\alpha,L\}} [\gamma(\mathbf{r}, \alpha, L)] \quad (7) $$

The size and orientation distribution of the fibers provide a vector representation of the fibers superimposed to the confocal image in Fig. 4e. A very good agreement is found between the visual appearance of the fibers on the ultrasound image and their vector representation provided by the fingerprint operator.

This information could be extremely rewarding to monitor neuro-muscular diseases³⁸ and myocardial fiber disarray that emerges in the early stage of many pathologies such as in cardiomyopathies or in fibrosis⁴⁰. It can also be relevant in the context of non-destructive testing for determining the orientation of elongated grains in polycrystalline materials⁴¹ or, in optical microscopy, for probing the collagen orientation and its lamellar distribution⁴².

## Discussion

In this manuscript, we have introduced the fingerprint operator with the aim to detect, localize and characterize any object in a strongly scattering medium. A necessary condition, however, is the creation of a sufficiently complete and precise knowledge of the reflection matrix and its dependence according to a state of interest. From that perspective, the
<figure><img src="image_6.png" /><figcaption>FIG. 4. Revealing the local architecture of fibers in muscle tissue. a, Sketch of the experimental configuration: An array of transducers is placed in contact with the calf of a healthy patient. The fingerprint operator $R_0(q)$ here corresponds to a set of reflection matrices associated with a mirror of size $L$ and orientation $\alpha$ simulated for each pixel $r$ in the medium considered as homogeneous ($c_0 = 1580$ m/s). b, Ultrasound image of the calf in dB. c, d, Likelihood index $\gamma$ in the ($L, \alpha$) parameter space for two positions yielding different values for the local correlation length $L$ and fiber orientation $\alpha$. The white arrows indicate the optimal length $L_{opt}$ and orientation $\alpha_{opt}$ at each position. e, Vector representation of the fibers deduced from the orientation distribution $\alpha_{opt}(\mathbf{r})$ encoded with the color scale displayed on the right.</figcaption></figure>

generalized polarization tensor, that was introduced ten years ago for electromagnetic waves, can be of interest$^{43,44}$. However, this approach requires an array of sensors fully surrounding the target under investigation, which is not available in most in-vivo or in-situ applications. On the contrary, the fingerprint operator works in reflection and can be adapted to any experimental configuration, thereby providing much higher robustness with respect to the non-ideal nature of experiments. Moreover, the simplicity of our approach makes it relatively inexpensive in terms of computational time and power (Methods). Indeed, the construction of a $\gamma$-map is computationally equivalent to a standard beamforming process used to build a confocal image (Supplementary Section S7).

An important question is the application of the fingerprint operator to dynamic imaging scenarios in which the target motion cannot be neglected during the measurement of the reflection matrix measurement. Again, the flexibility of the fingerprint matrix can help since it can be defined as a function of target motion parameters (motion direction, velocity,
acceleration). A maximization of the likelihood index $\gamma$ with respect to these parameters would then provide a quantification of target motion with a contrast enhancement equivalent to a static scenario but at the expense of a larger computational burden (Methods).

Beyond its versatility, the fingerprint operator is also a universal concept that can find applications in all fields of wave physics where multi-element technology is available. MIMO radar<sup>45</sup> and sonar<sup>46</sup> are examples of fields where target identification and localization in noisy environments is a long-lasting challenge. The fingerprint operator is also extremely promising for characterizing objects whose spatio-temporal response depends on their environment, thereby allowing a local measurement of intensive physical quantities, such as temperature or pressure. Microbubbles<sup>47-49</sup> in ultrasound imaging or plasmonic nanoparticles<sup>50</sup> in optics are ideal candidates for this purpose. While this is by no means an exhaustive list, the flexibility of the fingerprint operator concept can be critical to all of these applications.

<table><caption>Extended Data Tab. 1. **Experimental parameters.**</caption><thead><tr><th rowspan="2">Experiments</th><th>Elastic sphere</th><th>Muscle tissue</th></tr><tr><th>Lesion marker</th><th></th></tr></thead><tbody><tr><td>Probe</td><td>2D array</td><td>1D array</td></tr><tr><td>Number of transducers</td><td>32 x 32</td><td>256</td></tr><tr><td>Inter-element distance δu</td><td>0.5 mm</td><td>0.2 mm</td></tr><tr><td>Transducer directivity θ<sub>d</sub></td><td>28°</td><td>35°</td></tr><tr><td>Maximum Angle θ<sub>m</sub></td><td>28°</td><td>25°</td></tr><tr><td>Angular Pitch δθ</td><td>0.8°</td><td>0.5°</td></tr><tr><td>Number of plane waves N<sub>θ</sub></td><td>1225</td><td>101</td></tr><tr><td>Central frequency f<sub>c</sub></td><td>2.2 MHz</td><td>7.5 MHz</td></tr><tr><td>Bandwidth (at -6dB)</td><td>1.8-2.6 MHz</td><td>5.5-9.5 MHz</td></tr><tr><td>Sampling frequency f<sub>s</sub></td><td>6 MHz (IQ)</td><td>30 MHz</td></tr><tr><td>Recording time</td><td>130 μs</td><td>80 μs</td></tr><tr><td>Full acquisition time</td><td>160 ms</td><td>10 ms</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Extended Data Fig. 1. Elastic target signatures encoded in the reference reflection matrix. a, The reference reflection matrix $R_0$ is measured on the target sphere placed in water. The confocal beamforming process applied to $R_0$ selects not only the ballistic echo of the sphere but also its reverberations resulting from multiple reflections of bulk elastic waves (depicted by red arrows) at its inner surface. b, Matrix imaging decouples the input and output focal spots$^{51}$, $r_{in}$ and $r_{out}$, to highlight the contribution of circumferential waves (depicted by a black arrow) generated by the incoming wave at a specific angle of incidence with respect to sphere surface$^{33}$. c, Cross-section of the focused reflection matrix, $R_{0,xx}(y, z)$, in the plane $y = 0$ and at depth $z = 21$ mm showing the diagonal contribution of the ballistic echo. d, Same matrix but at depth $z = 28.5$ mm showing the off-diagonal contribution of circumferential waves. e, $(x, z)$-section of the confocal image in the plane $y = 0$ showing the spatio-temporal dispersion of the target echo. f, Likelihood index map $\gamma(r)$ (Eq. 2) built from the fingerprint operator indicating that we can accurately locate the target inside the reference environment (the sphere surface is highlighted by a red dashed line in panels e and f). Since in this case $R = R_0(r_0)$, this result serves as a consistency check for the formalism.</figcaption></figure>

# Methods

## Ultrasound scanners and probes.

In the elastic sphere and lesion marker experiments, the acquisition of the reflection matrix is performed using a 2D matrix array of transducers (Vermon) whose characteristics are provided in Extended Data Tab. 1. The electronic hardware used to drive the probe was developed by Supersonic Imagine in the context of a collaboration agreement with the
Langevin Institute.

In the calf experiment, the in vivo ultrasound data set has been collected by the SuperSonic Imagine company on an healthy volunteer from which informed consent had been obtained. Before being put at our disposal, this data set has been previously fully anonymized following standard practice defined by Commission nationale de l'information et des libertés (CNIL). The acquisition was performed using a medical ultrafast ultrasound scanner (Aixplorer Mach-30, Supersonic Imagine, Aix-en-Provence, France). The skin is placed in direct contact of a linear array of transducers (SL15-4, Supersonic Imagine) whose characteristics are also provided in Extended Data Tab. 1.

## Acquisition of the reflection matrix.

In each case, the reflection matrix is acquired using a set of plane waves³⁰. For each plane wave of angles of incidence $\theta_{in} = (\theta_x, \theta_y)$, the time-dependent reflected wave field $R(\mathbf{u}_{out}, \theta_{in}, t)$ is recorded by each transducer $\mathbf{u}_{out}$. This set of wave-fields forms a reflection matrix acquired in the plane wave basis, $\mathbf{R}_{\mathbf{u}\theta}(t) = [R(\mathbf{u}_{out}, \theta_{in}, t)]$. Since the transducer and plane wave bases can be related by a simple Fourier transform at the central frequency, the array pitch $\delta u$ and probe size $\Delta u$ dictate the angular pitch $\delta \theta$ and maximum angle $\theta_{max}$ necessary to acquire a full reflection matrix in the plane wave basis (Extended Data Tab. 1). A set of plane waves are thus generated by applying appropriate time delays $\Delta \tau(\theta_{in}, \mathbf{u}_{in})$ to each transducer $\mathbf{u}_{in} = (u_x, u_y)$ of the probe:

$$ \Delta\tau(\theta_{in}, \mathbf{u}_{in}) = [u_x \sin \theta_x + u_y \sin \theta_y]/c_0. \qquad (8) $$

where $c_0$ is the assumed speed-of-sound. The angular pitch and range as well as the number of illuminations are reported for each experiment in Extended Data Tab. 1.

## Focused beamforming of the reflection matrix.

To form an image from the measured reflection matrix, a beamforming procedure has to be applied in emission and reception in order to generate a synthetic focusing on each focal point. In the time domain, a delay-and-sum beamforming process can be performed by applying appropriate time delays to the recorded signals²³. In the frequency domain, matrix products can be applied to project $\mathbf{R}$ in the focused basis⁵¹. The projection of the $\mathbf{R}$-matrix in the focused basis can be performed at each depth $z$ by means of the following
matrix product:

$$ \mathbf{R}_{\rho\rho}(z, f) = \mathbf{G}_{\mathbf{u}\rho}^{\dagger}(z, f) \times \mathbf{R}_{\mathbf{u}\theta}(f) \times \mathbf{P}_{\theta\rho}^{*}(z, f). \quad (9) $$

where $\mathbf{P}_{\theta\rho}(z, f) = [P(\theta, \rho, z, f)]$ and $\mathbf{G}_{\mathbf{u}\rho}(z, f) = [G(\mathbf{u}, \rho, z, f)]$ are the free space propagation matrices from the focused basis ($\rho$) at depth $z$ to the plane wave ($\theta$) and transducer ($\mathbf{u}$) bases, respectively. Assuming a homogeneous speed of sound $c_0$, the coefficients of $\mathbf{P}_{\theta\rho}$ are expressed as follows

$$ P(\theta, \rho, f) = \exp \left[ ik \left( x \sin \theta_x + y \sin \theta_y + z \sqrt{1 - \sin^2 \theta_x - \sin^2 \theta_y} \right) \right]. \quad (10) $$

with $k = 2\pi f/c_0$, the wave number. The coefficients of $\mathbf{G}_{\mathbf{u}\rho}(z, f)$ are the 2D or 3D Green's function of the wave equation in homogeneous media:

$$ G_{2D}(\mathbf{u}, \rho, z, f) = \frac{-i}{4} \mathcal{H}_{0}^{(1)} \left( k \sqrt{\|\rho - \mathbf{u}\|^2 + z^2} \right), \quad (11a) $$

$$ G_{3D}(\mathbf{u}, \rho, z, f) = - \frac{\exp \left( ik \sqrt{\|\rho - \mathbf{u}\|^2 + z^2} \right)}{4\pi \sqrt{\|\rho - \mathbf{u}\|^2 + z^2}}. \quad (11b) $$

where $\mathcal{H}_0^{(1)}$ is the Hankel function of the first kind. Each coefficient of $\mathbf{R}_{\rho\rho}(f) = [R(\rho_{\text{out}}, \rho_{\text{in}}, f)]$ is the response between a virtual source at point $(\rho_{\text{in}}, z)$ and a virtual detector at $(\rho_{\text{out}}, z)$ at frequency $f$.

In order to retrieve the axial resolution provided by the broadband feature of ultrasonic signals, a broadband focused reflection matrix $\bar{\mathbf{R}}_{\rho\rho}(z)$ can be derived at each depth by coherently summing the monochromatic matrices over the frequency bandwidth:

$$ \bar{\mathbf{R}}_{\rho\rho}(z) = \frac{1}{\Delta f} \int_{f_{-}}^{f_{+}} df \mathbf{R}_{\rho\rho}(z, f), \quad (12) $$

where $f_{\pm} = f_c \pm \Delta f/2$ and $f_c$ is the central frequency. Each element of $\mathbf{R}_{\rho\rho}(z)$ contains the signal that would be recorded by a virtual transducer located at $(\rho_{\text{out}}, z)$ just after a virtual source at $(\rho_{\text{in}}, z)$ emits a pulse of length $\delta t = \Delta f^{-1}$ at the central frequency $f_c$. A time-gated confocal image $\mathcal{I}(\rho, z)$ can be extracted from the diagonal of the broadband focused $\mathbf{R}$-matrix, such that:

$$ \mathcal{I}(\rho, z) = |\bar{R}(\rho, \rho, z)|^2. \quad (13) $$
The resulting confocal images are displayed in Figs. 1d, Fig. 2, Fig. 3c and Fig. 4b for the different experiments described in the paper.

## Highlighting the elastic target signature

Even if the response of an elastic target may be spectrally and spatially quite complex<sup>15,16,32-34</sup>, all these features can be adequately captured by measuring the target's reflection matrix. We illustrate this in Extended Data Fig. 1 by showing the analysis of the reference reflection matrix associated with the 10-mm-diameter metal sphere immersed in water without the scattering glass bead packing (Supplementary Figure S4). In the following, we will refer to this reference (fingerprint) matrix as $R_0(t, q_0)$. The vector $q_0 = \{r_0, d\}$ here accounts for the target state in this experiment, i.e. its position $r_0 = (0, 0, 25)$ mm and diameter $d = 10$ mm. Extended Data Figure 1e displays the confocal image of the sphere obtained via the confocal beamforming algorithm described above and applied here to $R_0$, which highlights the complex signature of the target echo. Besides a specular component induced by the top of the sphere, this confocal image also shows a long temporal tail resulting from multiple reflections of bulk waves inside the sphere (Extended Data Fig. 1a, Supplementary Section S3). Beyond the confocal signal, the presence of circumferential Rayleigh-like waves propagating along the sphere surface<sup>52,53</sup> can be highlighted by decoupling the input and output focal spots and synthesizing the focused reflection matrix as described above (Extended Data Fig. 1b, Methods). Cross-sections $\overline{R}_{xx}(y, z)$ of the focused reflection matrix are displayed for the free space sphere experiment in the mid-plane $y = 0$ and at depths $z = 21$ mm and 28.5 mm in Extended Data Figs. 1c and d, respectively. While, at the first depth, the focused $R_0$-matrix displays a predominant confocal signal characteristic of the specular echo on the sphere cap (Extended Data Fig. 1c), off-diagonal echoes are observed at greater depth (Extended Data Fig. 1d). As shown by previous works for a cylinder<sup>32,33</sup> and confirmed by numerical simulations (Supplementary Section S3), these echoes correspond to the following paths: (i) excitation of a surface wave at an incident point $(\rho_{in}, z_0)$; (ii) propagation of this circumferential wave along the sphere surface until the symmetric point $(\rho_{out}, z_0) = (-\rho_{in}, z_0)$ where it is back-converted to a compressional wave in the surrounding fluid.
Fingerprint operator based on a calibration measurement

To build the fingerprint operator, the first strategy is to start from the measurement of a free-space reflection matrix $R_0$ associated with the target at a given position $r_0$. In principle, this matrix can be used as the value of the fingerprint operator at position $r_0$: $F(r_0) \equiv R_0$. In that case, the operation described in Eq. 2 can be seen as an adaptive filter. This kind of filter may be not the most adequate approach for detection and localization purposes.

First, depending on the experimental configuration, this reference matrix might be not necessarily specific enough with respect to the target environment (Supplementary Figure S5). In the sphere localization experiment, the ballistic echo of each sphere shall be removed in order for the fingerprint matrix not to be dominated by the strong specular echo generated by the interface of the granular medium in the original experiment (Fig. 1d). To that aim, the reference matrix $R_0$ is filtered in the time domain by means of a Heavyside filter (Supplementary Figure S6), such that:

$$R'_{0}(\mathbf{u}_{\text{out}}, \boldsymbol{\theta}_{\text{in}}, t) = R_{0}(\mathbf{u}_{\text{out}}, \boldsymbol{\theta}_{\text{in}}, t)H(t - \tau_{\text{in}}(\boldsymbol{\theta}_{\text{in}}, \mathbf{r}_{0}) - \tau_{\text{out}}(\mathbf{u}_{\text{out}}, \mathbf{r}_{0}) - \delta t) \quad (14)$$

with $\delta t \sim 1/\Delta f$, the temporal resolution of the recorded signals. $\tau_{\text{in}}$ and $\tau_{\text{out}}$ are the expected times-of-flight for the ballistic wave from the probe to the cylinder surface. At output (transducer basis), the time-of-flight $\tau_{\text{out}}$ can be expressed as follows:

$$\tau_{\text{out}}(\mathbf{u}, \mathbf{r}_0) = \frac{\sqrt{(x_0 - u_x)^2 + (y_0 - u_y)^2 + (z_0 - d/2)^2}}{c_0} \quad (15)$$

At input (plane wave basis), $\tau_{\text{in}}$ is given by

$$\tau_{\text{in}}(\boldsymbol{\theta}, \mathbf{r}_0) = \left[ x_0 \sin \theta_x + y_0 \sin \theta_y + (z_0 - d/2) \sqrt{1 - \sin^2 \theta_x - \sin^2 \theta_y} \right] / c_0. \quad (16)$$

Second, the transfer function of an adaptive filter is not flat. Both the temporal and spatial frequency spectrum of the target can thus be altered if we consider $F(r_0) \equiv R_0$ in Eq. 2. This would reduce both the axial and transverse resolution of the imaging method in the real space and decrease the contrast of the target in the $\gamma$-map of the target likelihood index. To optimize the target signal, an inverse filter could be applied by considering the inverse matrix of $R'_0$ as the building block of the fingerprint operator $F(r_0)$. However, this would be at the price of an extreme sensitivity with respect to experimental noise and a loss in terms of signal-to-noise ratio.
A good compromise between an adaptive and inverse filtering process can be reached by performing a spectral whitening of $\mathbf{R}'_0$. In practice, this can be done by performing the singular value decomposition of $\mathbf{R}'_0(f)$ at each frequency $f$ (Supplementary Figure S8):

$$ \mathbf{R}'_0(f) = \sum_p \sigma_p(f) \mathbf{U}_p(f) \mathbf{V}_p^\dagger(f) \quad (17) $$

where $\mathbf{U}_p$ and $\mathbf{V}_p$ are the output and input singular vectors of $\mathbf{R}'_0$. While, for a point-like target, the reflection matrix would only exhibit a predominant singular value, the spectrum of $\mathbf{R}'_0(f)$ displays a continuum of singular values due to the target size and the different resonant modes supported by each sphere (Supplementary Section S3).

To ensure the robustness of Eq. 2 with respect to experimental noise, only the eigenstates associated with the largest singular values shall be kept. In practice, the rank $N_S$ of the fingerprint operator is arbitrarily chosen such that the singular value ratio $\sigma_i/\sigma_1$ is larger than 0.4 (Supplementary Figure S8). The fingerprint operator is then obtained by whitening the singular value spectrum of this signal subspace such that:

$$ \mathbf{F}(\mathbf{r}_0, f) = \sum_{p=1}^{N_S} \mathbf{U}_p(f) \mathbf{V}_p^\dagger(f) \quad (18) $$

### Virtual shift of the target

The next step is to deduce the spatial evolution of the fingerprint operator from its value at $\mathbf{r}_0$. This is done by performing a virtual translation of the target on each point $\mathbf{r}$ of the field-of-view. To that aim, the fingerprint matrix shall be projected at output in the plane wave basis:

$$ \mathbf{F}_{\theta\theta}(\mathbf{r}_0, f) = \mathbf{P}_{\theta\mathbf{u}}(f) \times \mathbf{F}_{\mathbf{u}\theta}(\mathbf{r}_0, f) \quad (19) $$

From the plane wave basis, the target can be virtually shifted to any position $\mathbf{r} = \mathbf{r}_0 + \Delta\mathbf{r}$ by the application of a shift operator $\mathbf{S}(\Delta\mathbf{r})$ at input and output of the matrix $\mathbf{F}_{\theta\theta}(\mathbf{r}_0, f)$. Mathematically, it can be written as the following matrix product:

$$ \mathbf{F}_{\theta\theta}(\mathbf{r}, f) = \mathbf{S}(\Delta\mathbf{r}, f) \circ \mathbf{F}_{\theta\theta}(\mathbf{r}_0, f) \circ \mathbf{S}(\Delta\mathbf{r}, f) \quad (20) $$

where the symbol $\circ$ stands for the Hadamard (element wise) product. In terms of matrix coefficients, the previous equation writes

$$ F(\theta_{\text{out}}, \theta_{\text{in}}, \mathbf{r}, f) = S(\theta_{\text{out}}, \Delta\mathbf{r}, f) F(\theta_{\text{out}}, \theta_{\text{in}}, \mathbf{r}_0, f) S(\theta_{\text{in}}, \Delta\mathbf{r}, f) \quad (21) $$
where the coefficients $S(\theta, \Delta\mathbf{r}, f)$ of the shift operator $\mathbf{S}(\Delta\mathbf{r}, f)$ write

$$S(\theta, \Delta\mathbf{r}, f) = A(\theta, \mathbf{r})P(\theta, \mathbf{r}, f)P^*(\theta, \mathbf{r}_0, f) = A(\theta, \mathbf{r})P(\theta, \Delta\mathbf{r}, f) \quad (22)$$

where $A(\theta, \mathbf{r})$ is an apodization factor that limits the angular range of the synthetic aperture at emission and reception which depends on whether the plane wave can reach the target at point $\mathbf{r}$ or not (Supplementary Figure S9).

At last, the fingerprint operator can be projected back onto the acquisition basis of the recorded reflection matrix $\mathbf{R}$, such that

$$\mathbf{F}_{u\theta}(\mathbf{r}, f) = \mathbf{P}_{\theta u}^{\dagger}(f) \times \mathbf{F}_{\theta\theta}(\mathbf{r}, f) \quad (23)$$

## Numerical computation of the fingerprint operator

In the calf experiment, the fingerprint operator corresponds to a set of reflection matrices associated with a plane mirror of constant reflectivity characterized by a state $\mathbf{q}$ accounting for its position $\mathbf{r}$, orientation $\alpha$ and characteristic size $L$: $\mathbf{q} = \{\mathbf{r}, \alpha, L\}$ (Fig. 4a). Each reference matrix $\mathbf{R}_{u\theta}(\mathbf{q}, f)$ is computed numerically using the following matrix product:

$$\mathbf{R}_0(\mathbf{q}) = \mathbf{G}_{ur}(f) \times \mathbf{\Gamma}_{rr}(\mathbf{q}) \times \mathbf{G}_{ur}^{\top}(f) \times \mathbf{P}_{\theta u}^{\dagger}(f), \quad (24)$$

where $\mathbf{\Gamma}_{rr}(\mathbf{q})$ is a diagonal matrix whose diagonal coefficients $\gamma(\mathbf{r})$ stand for the reflectivity of the mirror in state $\mathbf{q}$. Note that, for such an object, the rank of the reflection matrix is equal to the number of resolution cells within the object, and that the $N_S$ non-zero singular values are degenerate$^{16}$. The adaptive filter $[\mathbf{F}(\mathbf{q}) = \bar{R}_0(\mathbf{q})]$ or inverse filter operation $[\mathbf{F}(\mathbf{q}) = \mathbf{R}_0^{-1}(\mathbf{q})]$ are therefore equivalent in this specific case.

## Computational insights

In terms of computational efficiency, the best strategy is to calculate the likelihood index $\gamma$ by considering the reflection matrix expressed in the plane wave basis. The computation of $\gamma$ is then equivalent to a confocal beamforming process applied to the Hadamard (element-wise) product between the reflection matrix $\mathbf{R}_{\theta\theta}(f)$ and the fingerprint operator $\mathbf{F}_{\theta\theta}(f)$ (Supplementary Section S7). Compared to a standard confocal beamforming process, the
computational burden of our approach only consists in performing this Hadamard product
between two matrices. For the 3D imaging experiments, such a product is performed in 3.5 s
with Matlab (R2021a). The duration of this operation is almost negligible compared to the
3D confocal beamforming process that takes 1.3 min on a working station with 2 processors
@2.20GHz, 128Go of RAM, and a GPU with 48 Go of dedicated memory<sup>23</sup>. The computation
time and power required to map the likelihood index $\gamma$ in Figs. 2 and 3 is therefore almost
equivalent to the beamforming of a confocal image.

For the 2D muscle experiment (Fig. 4), the Hadamard product between $\mathbf{R}_{\theta\theta}(f)$ and $\mathbf{F}_{\theta\theta}(f)$ is performed in 30 ms while the 2D confocal imaging process takes 80 ms on the same workstation. However, this operation shall be repeated for each couple $(\alpha, L)$. The number of tested angles $\alpha$ and correlation lengths $L$ being of 50 and 20, respectively, in the muscle experiment, the overall process thus takes approximately 10 s.

**Data availability.** The ultrasound data generated in this study are available at Zenodo⁵⁴ (https://zenodo.org/records/14845780).

**Code availability.** Codes used to post-process the optical data within this paper are available from the corresponding author.
**Acknowledgments.** The authors wish to thank the Somatex Company for providing the lesion marker. 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 no. 819261, REMINISCENCE project, A.A.). This project has also received funding from Labex WIFI (Laboratory of Excellence within the French Program Investments for the Future; ANR-10-LABX-24 and ANR-10-IDEX-0001-02 PSL*, M.F.) and from Agence Nationale de la Recherche (Grant No. ANR-22-ASTR-0020 under the AquaMat project, A.A.). L.M.R. was supported by the Austrian Science Fund (FWF) through Project no. P32300-N27 (WaveLand).

**Author Contributions Statement.** A.A. and S.R. initiated the project. A.A. supervised the project. A.L.B. and X.J. designed and performed the experiments on the granular medium. A.L.B. performed the experiment on the lesion marker experiment. A.L.B. developed the post-processing tools for the target detection experiment. W.L. performed the muscle tissue experiment. A.G. developed the post-processing tools for the muscle tissue experiment. A.L.B., A.G., L.M.R., S.R. and A.A. developed the concept of the fingerprint operator and performed the theoretical study. A.L.B. and A.G. prepared the figures. A.L.B., A.G. and A.A. prepared the manuscript. A.L.B., A.G., L.M.R., W.L., X.J., M.F., A.T., S.R. and A.A. discussed the results and contributed to finalizing the manuscript.

**Competing interests.** A.L.B., A.G., L.M.R., W.L., X.J., M.F., A.T., S.R. and A.A. are inventors on a french patent related to this work held by Supersonic Imagine and CNRS (no. FR2314789, filed December 2023). M.F. is cofounder of the SuperSonic Imagine company, which is commercializing one of the ultrasound platforms used in this study. W.L. is an employee of this company.
TABLE OF CONTENTS

• Section S1. Characterization of the dense granular suspension

• Section S2. Free-space reflection matrix

• Section S3. Back-scattered wave-field by an elastic sphere

• Section S4. Theoretical prediction of the target contrast

• Section S5. Precision of the localization process

• Section S6. Detection threshold of potential targets

• Section S7. Link between the mapping of the likelihood index and a confocal beam-forming process
# S1. CHARACTERIZATION OF THE DENSE GRANULAR SUSPENSION

To measure the extinction length of the ballistic or direct wave, a transmission configuration can be considered⁵⁵(Fig. S1). Two transducers (Olympus, model V382-SU ) are used to measure the coherent wave across the medium (circular aperture, diameter 12.7 mm, central frequency 3.5 MHz, bandwidth at -6 dB of the order of 70%). One is working as a source and emits a quasi incident plane wave :

$$ \psi_{\text{inc}}(t) = \psi_0 \cos(2\pi f_c t) \exp \left[ -\frac{t^2}{2\sigma^2} \right] \quad (\text{S1}) $$

with $f_c = 3$ MHz and $\sigma = 1/(2f_c)$. The other transducer is used a receiver. The two transducers are mounted opposite each other along a common central axis. The use of a wide transducer at reception is justified by the nature of the transmitted wave-field $\psi$ that is made of a coherent component $\langle \psi \rangle$ that resists to averaging over disorder and a multiple scattering coda $\delta\psi(\rho, z, t)$ wave that should vanish upon averaging over a number of independent disorder configurations ($\rho$ here accounts for the lateral position with respect to the $z$-axis). Taking advantage of spatial ergodicity, the coherent wave-field can be estimated by averaging the transmitted wave-field: (i) over the transducer area; (i) on different acquisitions separated by a mixing phase of the granular medium. The multiple scattering ‘coda’ forms a speckle pattern with each grain having an area comparable to the squared wavelength. Consequently, by choosing a circular transducer for reception whose active surface is 12.7 mm in diameter, the acquired signal is therefore averaged over around 500 speckle grains at 3 MHz which gives a reliable estimator of the coherent wave.

The receiver is mounted on a translation stage, which makes it possible to automate the measurement of the coherent wave for many distances $z$. In practice, we chose to vary $z$ from a few millimetres to around 40 mm, in steps of 0.1 mm, corresponding to almost 400 measurement points. To reduce electronic noise, the signals transmitted after many 2048 successive wave trains have been summed. A spatio-temporal window is then applied to the ultrasound data, so as to eliminate the multiple reflections between the transducers. A Fourier transform is finally performed in order to extract the Fourier-dependence of the coherent wave $\langle \psi \rangle (z, f)$.

Theoretically, the coherent wave-field can be expressed as follows:

$$ \langle \psi \rangle (z, f) = \psi_0(f) \exp \left( \frac{i2\pi f}{c_{\phi}(f)} z - \frac{z}{2\ell_{\text{ext}}(f)} \right) \quad (\text{S2}) $$
<figure><img src="image_8.png" /><figcaption>Fig. S1. Measuring the scattering mean free path $l_s$ in the granular material. Experimental set up dedicated to the measurement of the coherent wave across the granular medium. An acoustic wave generated by a transducer disk propagates through the medium and is recorded by a similar transducer disk centered on the first one and placed at a distance $z$ from it.</figcaption></figure>

with $c_φ$, the phase velocity and

$$l_{\text{ext}} = [l_s^{-1} + l_a^{-1}]^{-1}, \qquad (S3)$$

the extinction length, taking into account both the scattering losses characterized by a scattering mean free path $l_s$ and absorption losses quantified by a characteristic absorption length $l_a$. A measurement of $l_{\text{ext}}$ is then performed at each frequency $f$ by investigating the depth decay of $|\langle\psi\rangle(z, f)|$ (Fig. S2a). A linear fitting of $\ln(|\langle\psi\rangle(z, f)|)$ provides a measurement of $l_{\text{ext}}$ at each frequency (Fig. S2b).

The result is displayed in Fig. S3. In the low frequency regime ($k_0a \ll 1$), the value of $l_{\text{ext}}$ follows a Rayleigh scattering law and is therefore extremely dispersive. At frequency $f = 2.2$ MHz, the extinction length is $l_{\text{ext}} \simeq 1.5$ mm (black arrow in Fig. S3).

The wave velocity $c_φ$ has also been measured by investigating the phase of the coherent wave$^{56}$. A value close to the speed-of-sound $c_0$ in water is found at the central frequency $f_c = 2.2$ MHz: $c_φ(f_c) \simeq c_0 = 1.55$ mm/μs.

At last, the temporal dependence of the ensemble average transmitted intensity$^{57,58}$ has also been investigated to measure transport parameters in the granular medium in a higher frequency regime. Interestingly, such measurements provide a measurement of the absorption length $l_a \simeq 8$ mm$^{56}$ that we assume here as independent on frequency. Injecting the values measured for $l_a$ and $l_{\text{ext}}$ into Eq. S3 has led to the following estimation for the scattering
<figure><img src="image_9.png" /><figcaption>Fig. S2. Depth evolution of the frequency spectrum of the coherent wave $\langle\psi(z, f)\rangle$. a, Depth-frequency dependence of $|\langle\psi(z, f)\rangle|$. b, Depth dependence of $|\langle\psi(z, f)\rangle|$ normalized by its value at $z = 0$ at each frequency. The result is shown for a set of frequencies whose values are encoded with the color scale displayed on the right of the panel.</figcaption></figure>

mean free path: $\ell_s \simeq 1.8$ mm at the central frequency $f_c = 2.2$ MHz.
<figure><img src="image_10.png" /><figcaption>Fig. S3. Measurement of the extinction length in granular media. The measured extinction length $l_{ext}$ (blue line) is compared to the Rayleigh prediction (black dashed line). The black arrow indicates the scattering condition of the experiment depicted in Fig. 1</figcaption></figure>
## S2. FREE-SPACE REFLECTION MATRIX

The experimental set-up used to record the free-space reflection matrix (Extended Data Fig. 1) associated with each sphere is displayed in Fig. S4. The sphere is immersed in water and attached by two sewing threads glued on his surface. These threads can actually be seen in the ultrasound image shown in Extended Data Fig. 1e.

<figure><img src="image_11.png" /><figcaption>Fig. S4. Experimental set up used to record the free space reflection matrix of the metal sphere.</figcaption></figure>

Figure S5 shows the likelihood index map for sphere 1 using the raw matrix $R_0$ as the fingerprint operator ($F = R_0$). This operator completely fails to reveal the presence of sphere 1 because the direct echo of the sphere is not specific with respect to its environment. The interface of the scattering medium actually generates a strong specular echo that accounts for the surintensity of the $\gamma$-map lying along the dashed white line in Fig. S5a. There is a vertical shift of 5 mm with respect to the real position of the interface in the confocal image that corresponds to the radius of sphere 1.

As described in the Methods section, a time filtering of specular echoes should thus be performed on each free-space reflection matrix (Eq. 14), in order to make the fingerprint operator more specific. Figure S6 illustrates the effect of this filter: (i) on the raw signals by comparing the initial dataset before (Fig. S6a) and after (Fig. S6b) filtering of the specular sphere echo; (ii) on the confocal image extracted from the measured (Fig. S6c) and filtered (Fig. S6d) matrices, $R_0$ and $R'_0$, respectively. In the latter image, the specular echo that
<figure><img src="image_12.png" /><figcaption>Fig. S5. **Detrimental effect of specular echoes on the likelihood index map.** Confocal image (B&W scale) and superimposed $\gamma$-map (red scale) associated with the 10-mm diameter sphere. The $\gamma$-map is built using the raw matrix $\mathbf{R}_0$ as the fingerprint operator. The frequency bandwidth is 1.5-2.5 MHz. **a**, $(x, z)$— cross-section in the plane $y = y_1$. The white dashed line shows the artifact due the specular echoes on the granular medium surface. **b**, $(x, z)$— cross-section in the plane $x = x_1$. The positions of spheres 1 and 2 are represented by red and green dashed lines, respectively.</figcaption></figure>

initially appeared on the cap of the sphere in Fig. S6c is completely removed in Fig. S6d. The importance of this filter is highlighted by comparing the likelihood index $\gamma$ built from $\mathbf{R}'_0$ (Fig. S7a) and its initial counterpart (Fig. S5a). sphere 1 is now localized with a satisfying contrast $C_\gamma = 5.6$.

To improve this contrast, a singular value decomposition can then be applied to the filtered matrix $\mathbf{R}'_0(f)$ at each frequency $f$ (Eq. 17). The resulting singular values are displayed in Fig. S8. Only the eigenstates associated with singular values checking $\sigma_i > 0.4\sigma_1$ can be considered to build the fingerprint operator. The resulting $\gamma$-map is displayed in Fig. S7b. The contrast improvement is modest ($C_\gamma = 5.9$) but is clearly improved after whitening the singular value spectrum (Eq. 18). The corresponding $\gamma$-map shown in Fig. S7c displays a contrast $C_\gamma = 7.0$.

A last improvement consists in an angular decomposition of the fingerprint operator in order to filter the contribution of plane waves that cannot come from the target. This operation is illustrated by Fig. S9 and its result is shown in Fig. S7d. The final contrast $C_\gamma$
<figure><img src="image_13.png" /><figcaption>Fig. S6. **Time filtering of the specular sphere echo.** **a,** Time evolution of the reflected wave-field $R_0(\mathbf{u}_{\text{out}}, \theta_0, t)$ measured by the probe for the normal incident plane wave $\theta_0 = (0, 0)$. **b,** Time evolution of the corresponding filtered wave-field $R'_0(\mathbf{u}_{\text{out}}, \theta_0, t)$. **c,** $(x, z)$-section of the confocal image $\mathcal{I}$ (Eq. 13) in the plane $y = y_1$ built from the free-space matrix $R_0$. **d,** Same section of the confocal image $\mathcal{I}$ (Eq. 13) built from the filtered matrix $R'_0$. The data shown here are for sphere 1.</figcaption></figure>

reaches the excellent value of 8.2. Note that this value differs from the one reached in the accompanying paper (Fig. 2a, $C_\gamma$) because the frequency bandwidth differs: 1.5-3.5 MHz in Fig. S7d vs. 1.8-2.6 MHz in Fig. 2a. In the former case, the frequency-averaged value of the scattering mean free path $\ell_s$ (Fig. S3) is larger than in the latter case.
<figure><img src="image_14.png" /><figcaption>Fig. S7. Optimization of the fingerprint operator. The (x, z)-cross-section ($y = y_1$) of the $\gamma$-map for sphere 1 (red scale) superimposed to the confocal image (B&W scale) is computed by considering as a fingerprint operator: **a**, the time filtered matrix ($\mathbf{F} = \mathbf{R}'_0$); **b**, the $N_1$ first eigenstates of $\mathbf{R}'_0$ ($\mathbf{F} = \sum_{p=1}^{N_1} \sigma_p \mathbf{U}_p \mathbf{V}_p^\dagger$); **c**, after its singular spectrum whitening (Eq. 18); **d**, after the angular filtering described in Fig. S9. The frequency bandwidth is here 1.5-2.5 MHz.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. S8. Singular values of the free-space and filtered matrices $\mathbf{R}'_0(f)$ at frequency $f = 1.5$ MHz and 2.5 MHz. Their values averaged over the frequency bandwidth is also displayed. The horizontal dashed line accounts for the threshold applied to build the fingerprint operator $\mathbf{F}$ in Eq. 18. Each singular value is normalized by its maximum value $\sigma_1$.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. S9. Scheme describing the selection of plane waves kept in the fingerprint operator.
The orange area is a view of the selected angular sector. All the waves coming from outside this
area are removed by the angular filter.</figcaption></figure>
S3. BACK-SCATTERED WAVE-FIELD BY AN ELASTIC SPHERE

A COMSOL Multiphysics 6.0 simulation has been used to compute the wave-field scattered
by the metal sphere. As the mechanical properties of the sphere used in the experiment
are not known exactly, we use tabulated values for the simulation (Tab. S1). We also took
advantage of the rotation symmetry of the sphere to use a two-dimensional axisymmetric
simulation, which is less demanding in terms of resources and computing time. The aim
of the simulation is first to estimate the field induced by the sphere when illuminated by
an incident plane wave. In a second step, we will simulate focusing by coherently summing
different plane waves obtained from the rotation of the initial plane wave simulation. The
geometry, dimensions and physical quantities used for the simulation are all compiled in
Fig. S10 and Tab. S1.

<figure><img src="image_17.png" /><figcaption>Fig. S10. **Numerical simulation of the sphere response.** Scheme depicting the different parts of the simulation domain.</figcaption></figure>

Comsol simulation gives us access to a wide range of physical quantities both within
the sphere and in the water. We have chosen to focus on the pressure in the water and
<table><caption>TABLE S1. A set of parameters used to simulate the response of a metal sphere immersed in water in the presence of an incident plane wave.</caption><thead><tr><th></th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td rowspan="5">Medium (1)</td><td>Medium</td><td>Isotropic Solid</td></tr><tr><td><i>d</i></td><td>10 mm</td></tr><tr><td><i>ρ</i><sub>b</sub></td><td>5700 kg m<sup>-3</sup></td></tr><tr><td><i>c</i><sub>p</sub></td><td>2700 m s<sup>-1</sup></td></tr><tr><td><i>c</i><sub>t</sub></td><td>2000 m s<sup>-1</sup></td></tr><tr><td rowspan="4">Medium (2)</td><td>Medium Type</td><td>Liquid</td></tr><tr><td><i>L</i></td><td>60 mm</td></tr><tr><td><i>ρ</i><sub>0</sub></td><td>1000 kg m<sup>-3</sup></td></tr><tr><td><i>c</i><sub>0</sub></td><td>1500 m s<sup>-1</sup></td></tr><tr><td rowspan="2">Medium (3)</td><td>Medium</td><td>PML</td></tr><tr><td><i>e</i></td><td>2.5 mm</td></tr><tr><td rowspan="4">Simulation</td><td>Incident wave</td><td>Plane wave</td></tr><tr><td><i>f</i><sub>min</sub> = Δ<i>f</i></td><td>6 kHz</td></tr><tr><td><i>f</i><sub>max</sub></td><td>6 MHz</td></tr><tr><td>Frequency envelope</td><td>Hann window</td></tr></tbody></table>

the amplitude of the displacement field in the ball. These quantities are obtained for each frequency at any point in space. A Fourier transform is then computed after modulation by a Hann window to reconstruct a movie of wave propagation in the system in the plane $y = 0$. Figure S11 shows the incident plane wave propagating through water until it is reflected by the sphere. In addition to being directly reflected, part of the wave is transmitted into the sphere, as two wave fronts emerge within it. The first, moving faster, corresponds to the longitudinal wave, while the second is associated with the transverse wave. In addition to these bulk waves, circumferential waves are also generated, corresponding to Rayleigh waves. These circulate along the surface and propagate without stopping, so that they can turn around the sphere.

From this mother simulation, it is also possible to simulate focusing at any point of the simulation grid using the coherent sum of numerous plane waves. To simulate a three-
<figure><img src="image_18.png" /><figcaption>Fig. S11. Wave propagation movie for an incident plane wave. Different snapshots are shown for an incident plane wave, observed in the $y = 0$ plane. The incident wave propagates towards increasing $z$ (green arrows). Its interaction with the sphere gives rise to a specular back-reflection (white arrows). Bulk waves are also generated within the sphere (red arrows). Finally, Rayleigh-type surface waves are also generated (yellow arrows). The time $\delta t = 0$ corresponds to the passage of the plane wave through the $z = z_0$ plane (center of the sphere).</figcaption></figure>

dimensional focus from the parent simulation in axi-symmetric two-dimensional space, we rely on the rotation invariance of the system. To that aim, an observation grid is first defined for the simulated quantities. The idea is then to apply to the coordinates of this grid the rotation operation required to give the orientation of the wave vector of the parent simulation to the wave vector we are trying to simulate. For each plane wave to be simulated, we thus define a new associated grid that can be expressed in the cylindrical coordinates of the parent simulation. On each of these grids, the quantities of interest are finally estimated by interpolation from the grid of the parent simulation. At this stage, we now know the frequency dependence of the quantities of interest for all incident plane waves in a common grid.

### S3.1. Origin of the bright tail

To form a confocal image from the parent simulation, the transducer coordinates of the matrix probe can be chosen as the acoustic pressure measurement grid. In this way, the reflection matrix of the simulated system can be computed. A confocal image can be
deduced by applying a confocal beamforming process, as described in the Methods section and displayed in in Fig. S12. As already observed experimentally (Extended Data Fig. 1e), the confocal image shows the sphere cap but also a bright tail, whose origin can now be investigated.

<figure><img src="image_19.png" /><figcaption>Fig. S12. **Confocal image obtained from the simulated reflection matrix.** The (x, z)–section of the image is displayed in the plane y = 0.</figcaption></figure>

A film of wave propagation can be formed when all plane waves are suitably delayed to sum coherently at a point belonging to the bright tail. We are interested, for example, in the pixel **r**<sub>in</sub> = (0; 0; 25.5) mm and we set up an observation grid in the y = 0 plane, on which we are now able to know the frequency dependence of the quantities of interest for the various incident plane waves. In order to focus on **r**<sub>in</sub>, it remains to apply beamforming in the frequency domain (Methods) so that all the plane waves sum coherently in **r**<sub>in</sub> at the ballistic time (δt = 0). The wave propagation movie is finally obtained by summing the wave-field obtained for each plane wave at each point **r** and then calculating the inverse time Fourier transform after bandwidth modulation by a Hann window.

Snapshots of this movie are displayed in Fig. S13. The incident focused wave does indeed appear to converge towards the point **r**<sub>in</sub> as long as they propagate freely in the water. When they hit the surface of the ball, part of the energy is reflected, while another part is transmitted into the sphere. Longitudinal and transverse waves seem to undergo a refocusing process beneath the sphere surface and multiple reflections are observed at later time lapses.
<figure><img src="image_20.png" /><figcaption>Fig. S13. Wave propagation movie for an incident focused wave. Snapshots from the film of wave propagation in the system when focused on the point $r_{in} = (0; 0; 25.5)$ mm (green diamond) and observed in the plane $y = 0$. The incident wave (green arrows) gradually converges towards $r_{in}$. When the incident wave hits the sphere surface, a fraction of the energy is reflected (specular reflection, white arrows). Another fraction of the energy is converted into a whispering mode and then radiated back into the medium, (bulk reverberations, red arrows). The red dashed line represents the wave-front contributing to the reverberation observed in Fig. S12 at point $r_{in}$.</figcaption></figure>

This process transfers a significant amount of energy back to the water, in the direction of the probe.

To form the confocal image, the echoes from $r_{in}$ focus are precisely summed, selecting those corresponding to the supposed time delay between $r_{in}$ and each transducer element. This echo highlighted by a red dashed line in Fig. S13 corresponds here to the first reflection on the back of the sphere. In the confocal image (Fig. S12), the bright tail thus exists thanks to the multiply-reflected bulk waves inside the sphere generated by incipient waves intended to focus at these points and then radiated back towards the probe with the time delay expected in a homogeneous medium.

### S3.2. Origin of off-diagonal echoes

To better understand the origin of off-diagonal signal in the focused reflection matrix (Extended Data Fig. 1d), we can once again rely on the simulation described above. Focused beamforming can be applied at input and output of the simulated reflection matrix in
order to synthesize the focused reflection matrix at different depths (Method). As observed experimentally (Extended Data Fig. 1d), we retrieve for certain depth strong off-diagonal echoes (Fig. S14).

<figure><img src="image_21.png" /><figcaption>Fig. S14. **Origin of off-diagonal echoes.** **a**, Confocal image in the *y* = 0 plane. **b**, *x* cross-section of the focused reflection matrix in the focused basis in the plane *y*<sub>in</sub> = *y*<sub>out</sub> = 0 and at depth *z* = 28.5 mm. **c**, Snapshots of the wave propagation movie for an incident wave focusing on the point **r**<sub>in</sub> = (1.25; 0; 28.5) mm (green diamond) and observed in the plane *y* = 0. The incident wave-front (green arrows) gradually converges towards **r**<sub>in</sub>. When the incident wave hits the surface, a fraction of the energy is reflected and gives rise to the specular echo (white arrows). Another fraction of the energy is converted into a circumferential surface wave (yellow arrows) and is then radiated back towards the probe (red arrows). The red dashed lines represent the echoes in the *y* = 0 plane actually used by the probe to focus on the output point **r**<sub>out</sub> = (-ρ<sub>in</sub>, *z*) (red diamond). The chosen couple of points (*r*<sub>in</sub>, *r*<sub>out</sub>) corresponds to the strong off-diagonal echo displayed in **b**.</figcaption></figure>
Remarkably, the simulation once again allows us to track the origins of this off-diagonal echo, by simulating an incident wave focusing on the corresponding point $\mathbf{r}_{\text{in}}$. Snapshots of the corresponding movie are shown in Fig. S14c. It shows once again the incident wave converging towards the point $\mathbf{r}_{\text{in}}$ before hitting the sphere surface. A predominant surface wave is generated and propagates around the sphere (yellow arrow). The wave-front contributing to the off-diagonal echo in Fig. S14b is highlighted by a red dashed line. It indeed corresponds to the echo produced by the circumferential wave.

## S4. THEORETICAL PREDICTION OF THE TARGET CONTRAST

In the first experiment described in the accompanying paper (Fig. 1), the recorded reflection matrix $\mathbf{R}$ can be decomposed into a target component $\sigma_T \mathbf{R}_0$, and a multiple scattering contribution $\mathbf{M}$:

$$ \mathbf{R} = \sigma_T \mathbf{R}_0 + \mathbf{M} \quad (S4) $$

These two matrices are considered as fully uncorrelated. The single scattering contribution of the granular medium is neglected since multiple scattering is strongly predominant at the targets' depths in Fig. 1. We will consider a normalized matrix $\mathbf{R}_0$ such that $\sigma_T$ accounts for the target signal:

$$ \langle |R_0(\mathbf{u}_{\text{out}}, \boldsymbol{\theta}_{\text{in}}, f)|^2 \rangle = 1, $$

where the symbol $\langle \cdots \rangle$ stands for ensemble average.

### S4.1. Target contrast in confocal imaging

Injecting Eq. 9 and 12 into Eq. 13 leads to the following expression for the mean confocal intensity:

$$ \mathcal{I} = \left\langle \left| \frac{1}{\Delta f} \int_{f_{-}}^{f_{+}} df \sum_{\boldsymbol{\theta}_{\text{in}}}^{N_{\theta}} \sum_{\mathbf{u}_{\text{out}}}^{N_{u}} R(\mathbf{u}_{\text{out}}, \boldsymbol{\theta}_{\text{in}}, f) P^{*}(\boldsymbol{\rho}, \boldsymbol{\theta}_{\text{in}}, f) G^{*}(\boldsymbol{\rho}, \mathbf{u}_{\text{out}}, f) \right|^2 \right\rangle. \quad (S5) $$

In the multiple scattering regime, confocal beamforming is an incoherent process: Each term in the triple sum of Eq. S5 can be seen as a random phasor. The corresponding multiple scattering contribution can be rewritten as follows:

$$ \mathcal{I}_M = \frac{\delta f}{\Delta f} \sum_{\boldsymbol{\theta}_{\text{in}}} \sum_{\mathbf{u}_{\text{out}}}^{N_{\theta}} \langle |M(\mathbf{u}_{\text{out}}, \boldsymbol{\theta}_{\text{in}}, f)|^2 \rangle \quad (S6) $$
where $\delta f$ is the correlation frequency of the multiple scattering noise. Let $\sigma_M^2 = \langle |M(\mathbf{u}_{\text{out}}, \boldsymbol{\theta}_{\text{in}}, f)|^2 \rangle$ be the power of multiple scattering noise recorded by the probe. The mean multiple scattering intensity thus scales as follows:

$$ \mathcal{I}_M \sim \frac{N_u N_\theta}{N_f} \sigma_M^2 \qquad (S7) $$

with $N_f = \Delta f / \delta f$, the number of independent frequency grains in the frequency bandwidth. The confocal beamforming process amounts to increase the multiple scattering intensity by a factor $N_\theta$ due to the beamforming at input, by a factor $N_u$ due to the beamforming at output. Moreover, the different frequency components of the wave-field do not sum coherently which leads to a decrease of the multiple scattering intensity by the number $N_f$ of independent frequency grains in the bandwidth.

To derive an equivalent scaling for the target intensity, we will first consider the case of a point-like target. In that case, confocal beamforming is a perfectly coherent process with respect to the singly-scattered echo of the target: Each term in the triple sum of Eq. S5 adds constructively. The associated confocal intensity $\mathcal{I}_P$ thus scales as follows:

$$ \mathcal{I}_P \sim N_\theta^2 N_u^2 \sigma_T^2 \qquad (S8) $$

The confocal beamforming process amounts to increase the point-like target intensity $\mathcal{I}_P$ by a factor $N_\theta^2$ due to the focused beamforming at input and by a factor $N_u^2$ due to the focused beamforming at output.

For a more complex target, the phasors in Eq. S5 are only partially coherent. This case is therefore intermediate between the multiple scattering component and the point-like target case. At each frequency, only one eigenstate of the reflection matrix, its specular component, leads to a coherent sum in the confocal beamforming process. If we assume that the target matrix exhibits a step-like distribution of singular values with $N_S$ the number of non-zero singular values, the confocal intensity should thus be decreased by a factor $N_S$ compared to the case point-like target. Moreover, the target exhibits a partially incoherent frequency spectrum. The associated intensity $\mathcal{I}_T$ will be also decreased by a factor $N_T$ compared to the point-like target case, $N_T$, being the number of independent coherence grains in the frequency bandwidth. The scaling $\mathcal{I}_T$ is therefore as follows:

$$ \mathcal{I}_T \sim \frac{\mathcal{I}_P}{N_S N_T} \sim \frac{N_\theta^2 N_u^2}{N_S N_T} \sigma_T^2 \qquad (S9) $$
where N is the number of transducers. $N_S$ and $N_T$ are the number of spatial and temporal degrees of freedom exhibited by the target.

In the diffusive regime, the contrast between the direct echo of the target and the multiple scattering background in the confocal image is therefore given by:

$$C_I = \frac{\mathcal{I}_T}{\mathcal{I}_M} \sim \frac{N_\theta N_u N_f \sigma_T^2}{N_S N_T \sigma_M^2} \qquad (S10)$$

### S4.2. Target contrast provided by matrix fingerprint imaging

To evaluate the signal-to-noise ratio provided by the matrix fingerprint imaging process, we will consider Eq. 2 with a fingerprint operator equal to the target matrix: $\mathbf{F} = \mathbf{R}_0$. The mean intensity $\mathcal{F}$ of the fingerprint image is then given by:

$$\mathcal{F} = \beta^{-1} \left\langle \left| \frac{1}{\Delta f} \int_{f_{-}}^{f_{+}} df \sum_{\boldsymbol{\theta}_{\text{in}}}^{N_{\theta}} \sum_{\mathbf{u}_{\text{out}}}^{N_{u}} R(\mathbf{u}_{\text{out}}, \boldsymbol{\theta}_{\text{in}}, f) R_{0}^{*}(\mathbf{u}_{\text{out}}, \boldsymbol{\theta}_{\text{in}}, f) \right|^2 \right\rangle \qquad (S11)$$

where $\beta$ stands for the normalization factor of Eq. 2.

For the multiple scattering contribution, the associated mean intensity $\mathcal{F}_M$ can be computed by considering $\mathbf{R} = \mathbf{M}$ in Eq. S11

$$\mathcal{F}_M = \beta^{-1} \left\langle \left| \frac{1}{\Delta f} \int_{f_{-}}^{f_{+}} df \sum_{\boldsymbol{\theta}_{\text{in}}}^{N_{\theta}} \sum_{\mathbf{u}_{\text{out}}}^{N_{u}} M(\mathbf{u}_{\text{out}}, \boldsymbol{\theta}_{\text{in}}, f) R_{0}^{*}(\mathbf{u}_{\text{out}}, \boldsymbol{\theta}_{\text{in}}, f) \right|^2 \right\rangle \qquad (S12)$$

The reflection matrix coefficients of the multiple scattering contribution and of the target component are fully uncorrelated. The triple sum is therefore incoherent in Eq. S12:

$$\mathcal{F}_M = \beta^{-1} \frac{1}{N_f} \sum_{\boldsymbol{\theta}_{\text{in}}}^{N_{\theta}} \sum_{\mathbf{u}_{\text{out}}}^{N_u} \langle |M(\mathbf{u}_{\text{out}}, \boldsymbol{\theta}_{\text{in}}, f)|^2 \rangle \langle |R_0(\mathbf{u}_{\text{out}}, \boldsymbol{\theta}_{\text{in}}, f)|^2 \rangle \qquad (S13)$$

which simplifies into:

$$\mathcal{F}_M \sim \frac{N_\theta N_u}{\beta N_f} \sigma_M^2 \qquad (S14)$$

For the target contribution, the associated fingerprint intensity $\mathcal{F}_T$ can be derived by considering $\mathbf{R} = \sigma_T \mathbf{R}_0$ in Eq. S11:

$$\mathcal{F}_T = \frac{\sigma_T^2}{\beta} \left\langle \left| \frac{1}{\Delta f} \int_{f_{-}}^{f_{+}} df \sum_{\boldsymbol{\theta}_{\text{in}}}^{N_{\theta}} \sum_{\mathbf{u}_{\text{out}}}^{N_{u}} |R_0(\mathbf{u}_{\text{out}}, \boldsymbol{\theta}_{\text{in}}, f)|^2 \right|^2 \right\rangle \qquad (S15)$$

Each term in the triple sum is coherent, which yields the following scaling for the target intensity $\mathcal{F}_T$:

$$\mathcal{F}_T \sim \frac{N_\theta^2 N_u^2}{\beta} \sigma_T^2 \qquad (S16)$$
The contrast between the target signal with respect to the multiple scattering background (Eq. S14) in the likelihood image is therefore given by:

$$C_{\gamma} = \frac{\mathcal{F}_{T}}{\mathcal{F}_{M}} \sim N_{\theta} N_{u} N_{f} \frac{\sigma_{T}^{2}}{\sigma_{M}^{2}} \qquad (S17)$$

Compared to the confocal image (Eq. S10), the gain G provided by the fingerprint operator in terms of contrast between the target and the multiple scattering fog scales as:

$$G = \frac{C_{\gamma}}{C_{I}} \sim N_{S} N_{T} \qquad (S18)$$

G corresponds to the number of spatial and temporal d.o.f exhibited by the target echo.

## S5. PRECISION OF THE LOCALIZATION PROCESS

As shown by Desailly and collegues³⁶ in the context of ultrasound localization microscopy⁴⁹, the transverse and axial localization precision can be expressed by means of the Cramér-Rao bound as follows:

$$\delta\rho_{L} \sim 2\sqrt{3}\frac{\sigma_{t}}{f_{c}}\frac{\delta\rho_{D}}{\sqrt{N_{u}N_{\theta}}} \qquad (S19)$$

and

$$\delta z_{L} \sim \frac{c\sigma_{t}}{2}\frac{1}{\sqrt{N_{u}N_{\theta}}} \qquad (S20)$$

where $\sigma_t$ is the standard deviation $\sigma_t$ of echo time estimates. This quantity has been derived in a seminal paper by Quazi for a limited bandwidth ($\Delta f \ll f_c$)³⁵:

$$\sigma_t = \frac{1}{2\sqrt{2\pi}}\sqrt{\frac{\sigma_M}{\sigma_T}}\frac{f_c}{\sqrt{N_f}} \qquad (S21)$$

Injecting Eq. S21 into Eq. S19 leads to the following expression for $\delta\rho_L$:

$$\delta\rho_{L} \sim \frac{1}{\pi}\sqrt{\frac{3}{2}}\sqrt{\frac{\sigma_{M}}{\sigma_{T}}} \frac{\delta\rho_{D}}{\sqrt{N_{f}N_{u}N_{\theta}}} \qquad (S22)$$

Using Eq. S17, the previous equation can be recast as follows:

$$\delta\rho_{L} \sim \frac{1}{\pi}\sqrt{\frac{3}{2}}\frac{\delta\rho_{D}}{\sqrt{C_{\gamma}}} \qquad (S23)$$

As regards to the axial resolution $\delta z_L$, Eqs. S21 and S17 can be injected into Eq. S20 to obtain:

$$\delta z_{L} \sim \frac{1}{4\sqrt{2\pi}}\frac{c f_{c}}{\sqrt{C_{\gamma}}} \qquad (S24)$$
The transverse and axial precision of the localization process are therefore inversely proportional to the target contrast. This fundamental result shows the benefit of the fingerprint operator for localization purposes.

## S6. DETECTION THRESHOLD OF POTENTIAL TARGETS

Since the scattering medium is considered as one realization of a random process, assessing the reliability of the obtained images requires a statistical model for the probability density function of the $\gamma$-map in the multiple scattering regime. If the reflection matrix coefficients are independently and identically distributed random variables, the intensity $\gamma^2$ of each image pixel would obey the Rayleigh distribution in absence of any target:

$$ \rho_R(\gamma^2) = \sqrt{\frac{2}{\pi}\frac{\gamma^2}{\gamma_M^2}} \exp\left(-\frac{\pi}{4}\frac{\gamma^4}{\gamma_M^4}\right) \qquad (S25) $$

with $\gamma_M = \sqrt{\langle\gamma^2\rangle}$, the mean value of $\gamma$ associated with the multiple scattering background. As a detection threshold based on the maximum value $\gamma_{max}$ of the $\gamma$-map is needed, the statistical distribution of $\gamma_{max}^2$ has to be known in the multiple scattering regime. As each pixel are independent from each other in the multiple scattering regime, then the distribution function $F_{max}$ of $\gamma_{max}^2$ is simply given by the $L^{th}$ power of the distribution function $F_R$ of one pixel,

$$ F_R(\gamma^2) = 1 - \exp\left(-\frac{\pi}{4}\frac{\gamma^4}{\gamma_M^4}\right) \qquad (S26) $$

with $L$, the number of independent speckle spots in the field-of-view. The quantity $1 - [F_R(\gamma_{max}^2)]^L$ is the probability that a pure multiple-scattering speckle gives rise to a likelihood index larger than the measured $\gamma_{max}$. Therefore, $[F_R(\gamma_{max}^2)]^L$ can be directly used as a probability that the brightest pixel is actually associated with a coherent target rather than a normal multiple-scattering speckle fluctuation. In the experiment considered in Figure 3, we have $C_\gamma = (\gamma_{max}/\gamma_M)^2 \sim 6$ and $L \sim 25^3$. Applying these numerical values to Eq. S26 yields a probability of $10^{-8}$ that the brightest spot is associated with a false alarm (see Fig. S15). Each sphere is therefore detected with an extremely high degree of confidence.
<figure><img src="image_22.png" /><figcaption>Fig. S15. Probability of false alarm. The probability of false alarm, $1 - |F_R(\gamma^2)|^L$, is displayed as a function of $C_\gamma = (\gamma_{max}/\gamma_M)^2$ for different values of $L$, the number of independent pixels in the field-of-view: $L = 10^3$ (green), $L = 25^3$ (blue), $L = 100^3$ (red).</figcaption></figure>

## S7. LINK BETWEEN THE MAPPING OF THE LIKELIHOOD INDEX AND A CONFOCAL BEAMFORMING PROCESS

In terms of computational efficiency, the best strategy is to compute the likelihood index $\gamma$ directly in the plane wave basis. Indeed, Eq. 2 can then be written in terms of matrix coefficients as follows:

$$ \gamma(\mathbf{r}) = \beta^{-1} \int df \sum_{\theta_{\text{in}}} \sum_{\theta_{\text{out}}} R(\theta_{\text{out}}, \theta_{\text{in}}, f) F^*(\theta_{\text{out}}, \theta_{\text{in}}, \mathbf{r}, f) \quad (\text{S27}) $$

with $\beta$, a normalization coefficient. Injecting Eq. 21 into this last equation leads to the following expression for the likelihood index

$$ \gamma(\mathbf{r}) = \beta^{-1} \int df S(\theta_{\text{out}}, \Delta\mathbf{r}, f) R(\theta_{\text{out}}, \theta_{\text{in}}, f) F^*(\theta_{\text{out}}, \theta_{\text{in}}, \mathbf{r}_0, f) S(\theta_{\text{in}}, \Delta\mathbf{r}, f). \quad (\text{S28}) $$

The last equation can be rewritten as the following matrix product:

$$ \gamma(\mathbf{r}) = \beta^{-1} \int df \mathbf{S}^{\dagger}(\mathbf{r} - \mathbf{r}_0, f) \times [\mathbf{R}_{\theta\theta}(f) \circ \mathbf{F}_{\theta\theta}^*(\mathbf{r}_0, f)] \times \mathbf{S}^{\dagger}(\mathbf{r} - \mathbf{r}_0, f) \quad (\text{S29}) $$

This last equation shows a strong analogy with a confocal beamforming process that write as follows when operated from the plane wave basis:

$$ \mathcal{I}(\mathbf{r}) = \beta^{-1} \int df \mathbf{S}^{\dagger}(\mathbf{r}, f) \times \mathbf{R}_{\theta\theta}(f) \circ \mathbf{S}^{\dagger}(\mathbf{r}, f) \quad (\text{S30}) $$
The comparison between Eqs. S29 and S30 shows that the spatial mapping of the likelihood index is equivalent to a confocal beamforming process applied to the Hadamard product between the reflection matrix $\mathbf{R}_{\theta\theta}(f)$ and the fingerprint operator $\mathbf{F}_{\theta\theta}(f)$. Compared to standard confocal beamforming, the computational burden of our approach only consists in the computation of this Hadamard product.
References.

[1] M. Fink, Time reversed acoustics, Physics Today **50**, 34 (1997).

[2] 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).

[3] G. Foschini and M. Gans, On limits of wireless communications in a fading environment when using multiple antennas, Wireless Personal Communications **6**, 311 (1998).

[4] S. M. Popoff et al., Measuring the transmission matrix in optics: An approach to the study and control of light propagation in disordered media, Phys. Rev. Lett. **104**, 100601 (2010).

[5] C. Prada and M. Fink, Eigenmodes of the time reversal operator: A solution to selective focusing in multiple-target media, Wave Motion **20**, 151 (1994).

[6] S. M. Popoff et al., Exploiting the time-reversal operator for adaptive optics, selective focusing, and scattering pattern analysis, Phys. Rev. Lett. **107**, 263901 (2011).

[7] S. Rotter and S. Gigan, Light fields in complex media: Mesoscopic scattering meets wave control, Rev. Mod. Phys. **89**, 015005 (2017).

[8] H. Cao, A. P. Mosk, and S. Rotter, Shaping the propagation of light in complex media, Nat. Phys. **18**, 994 (2022).

[9] J. Bertolotti and O. Katz, Imaging in complex media, Nat. Phys. **18**, 1008 (2022).

[10] B. Gérardin, J. Laurent, A. Derode, C. Prada, and A. Aubry, Full transmission and reflection of waves propagating through a maze of disorder, Phys. Rev. Lett. **113**, 173901 (2014).

[11] S. M. Popoff, A. Goetschy, S. F. Liew, A. D. Stone, and H. Cao, Coherent Control of Total Transmission of Light through Disordered Media, Phys. Rev. Lett. **112**, 133903 (2014).

[12] M. Davy, Z. Shi, J. Park, C. Tian, and A. Z. Genack, Universal structure of transmission eigenchannels inside opaque media, Nat. Commun. **6**, 6893 (2015).

[13] A. Aubry and A. Derode, Detection and imaging in a random medium: A matrix method to overcome multiple scattering and aberration, J. Appl. Phys. **106**, 044903 (2009).

[14] A. Badon et al., Smart optical coherence tomography for ultra-deep imaging through highly scattering media, Sci. Adv. **2**, e1600370 (2016).

[15] A. Aubry, J. de Rosny, J.-G. Minonzio, C. Prada, and M. Fink, Gaussian beams and legendre polynomials as invariants of the time reversal operator for a large rigid cylinder, J. Acoust. Soc. Am. **120**, 2746 (2006).

[16] J.-L. Robert and M. Fink, The prolate spheroidal wave functions as invariants of the time
reversal operator for an extended scatterer in the Fraunhofer approximation, J. Acoust. Soc. Am. **125**, 218 (2009).

[17] S. Yoon, H. Lee, J. H. Hong, Y.-S. Lim, and W. Choi, Laser scanning reflection-matrix microscopy for aberration-free imaging through intact mouse skull, Nat. Commun. **11**, 5721 (2020).

[18] Y. Jo et al., Through-skull brain imaging in vivo at visible wavelengths via dimensionality reduction adaptive-optical microscopy, Sci. Adv. **8**, eabo4366 (2022).

[19] H. Lee et al., High-throughput volumetric adaptive optical imaging using compressed time-reversal matrix, Light Sci Appl **11**, 16 (2022).

[20] G. Weinberg, E. Sunray, and O. Katz, Noninvasive megapixel fluorescence microscopy through scattering layers by a virtual incoherent reflection matrix, Sci. Adv. **10**, eadl5218 (2024).

[21] A. Badon et al., Distortion matrix concept for deep optical imaging in scattering media, Sci. Adv. **6**, eaay7170 (2020).

[22] 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. **117**, 14645 (2020).

[23] F. Bureau et al., Three-dimensional ultrasound matrix imaging, Nat. Commun. **14**, 6793 (2023).

[24] G. Murray et al., Aberration free synthetic aperture second harmonic generation holography, Opt. Express **31**, 32434 (2023).

[25] E. Giraudat et al., Matrix imaging as a tool for high-resolution monitoring of deep volcanic plumbing systems with seismic noise, Commun. Earth Environ. **5**, 509 (2024).

[26] Y. Zhang et al., Adaptive optical multispectral matrix approach for label-free high-resolution imaging through complex scattering media, Adv. Photonics **7**, 046008 (2025).

[27] U. Najar et al., Harnessing forward multiple scattering for optical imaging deep inside an opaque medium, Nat. Commun. **15**, 7349 (2024).

[28] P. Pai, J. Bosch, M. Kühmayer, S. Rotter, and A. P. Mosk, Scattering invariant modes of light in complex media, Nat. Photon. **15**, 431 (2021).

[29] S. van den Wildenberg, X. Jia, J. Léopoldès, and A. Tourin, Ultrasonic tracking of a sinking ball in a vibrated dense granular suspension, Sci. Rep. **9**, 5460 (2019).

[30] G. Montaldo, M. Tanter, J. Bercoff, N. Benech, and M. Fink, Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography, IEEE Trans. Ultrason.
Ferroelectr. Freq. Control **56**, 489 (2009).

[31] A. Goïcoechea et al., Reflection measurement of the scattering mean free path at the onset of multiple scattering, Phys. Rev. Lett. **133**, 176301 (2024).

[32] J.-L. Thomas, P. Roux, and M. Fink, Inverse scattering analysis with an acoustic time-reversal mirror, Phys. Rev. Lett. **72**, 637 (1994).

[33] C. Prada and M. Fink, Separation of interfering acoustic scattered signals using the invariants of the time-reversal operator. application to lamb waves characterization, J. Acoust. Soc. Am. **104**, 801 (1998).

[34] N. Gespa and H. Überall, *La Diffusion Acoustique Par Des Cibles Élastiques de Forme Géométrique Simple: Théories Et Expériences* (CEDOCAR, 1987).

[35] A. Quazi, An overview on the time delay estimate in active and passive systems for target localization, IEEE Trans. Acoust. **29**, 527 (1981).

[36] Y. Desailly, J. Pierre, O. Couture, and M. Tanter, Resolution limits of ultrafast ultrasound localization microscopy, Phys. Med. Biol. **60**, 8723 (2015).

[37] A. M. Rüland et al., Using a new marker clip system in breast cancer: Tumark vision clip - feasibility testing in everyday clinical practice, Breast Care **13**, 114 (2018).

[38] J. Wijntjes and N. van Alfen, Muscle ultrasound: Present state and future opportunities, Muscle & Nerve **63**, 455 (2020).

[39] C. Papadacci et al., Imaging the dynamics of cardiac fiber orientation in vivo using 3d ultrasound backscatter tensor imaging, Sci. Rep. **7**, 830 (2017).

[40] W. I. Tseng, J. Dou, T. G. Reese, and V. J. Wedeen, Imaging myocardial fiber disarray and intramural strain hypokinesis in hypertrophic cardiomyopathy with mri, J. Magn. Reson. Imaging **23**, 1 (2005).

[41] R. B. Thompson et al., Scattering of elastic waves in simple and complex polycrystals, Wave Motion **45**, 655 (2008).

[42] C. Raoux, A. Chessel, P. Mahou, G. Latour, and M.-C. Schanne-Klein, Unveiling the lamellar structure of the human cornea over its full thickness using polarization-resolved shg microscopy, Light Sci. Appl. **12**, 190 (2023).

[43] H. Ammari et al., Target identification using dictionary matching of generalized polarization tensors, Fond. Comput. Math. **14**, 27 (2013).

[44] H. Ammari, M. P. Tran, and H. Wang, Shape identification and classification in echolocation,
Siam. J. Imaging Sci. **7**, 1883 (2014).

[45] L. Xu, J. Li, and P. Stoica, Target detection and parameter estimation for mimo radar systems, IEEE Trans. Aerosp. Electron. Syst. **44**, 927 (2008).

[46] Y. Pailhas, J. Houssineau, Y. R. Petillot, and D. E. Clark, Tracking with mimo sonar systems: applications to harbour surveillance, IET Radar Sonar Navig. **11**, 629 (2017).

[47] N. De Jong, A. Bouakaz, and P. Frinking, Basic acoustic properties of microbubbles, Echocardiography **19**, 229 (2002).

[48] C. Tremblay-Darveau, R. Williams, and P. N. Burns, Measuring absolute blood pressure using microbubbles, Ultrasound Med. Biol. **40**, 775 (2014).

[49] C. Errico et al., Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging, Nature **527**, 499 (2015).

[50] J. N. Anker et al., Biosensing with plasmonic nanosensors, Nat. Mater. **7**, 442 (2008).

[51] 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).

[52] D. Royer, E. Dieulesaint, X. Jia, and Y. Shui, Optical generation and detection of surface acoustic waves on a sphere, Appl. Phys. Lett. **52**, 706 (1988).

[53] D. Clorennec and D. Royer, Investigation of surface acoustic wave propagation on a sphere using laser ultrasonics, Appl. Phys. Lett. **85**, 2435 (2004).

[54] A. Le Ber et al., Fingerprint matrix imaging [data]. Zenodo (2024).

[55] J. H. Page, P. Sheng, H. P. Schriemer, I. Jones, J. Xiaodun, and D. A. Weitz, Group Velocity in Strongly Scattering Media, Science **271**, 634 (1996).

[56] A. Le Ber, *Approche matricielle de la propagation des ultrasons dans les suspensions granulaires*, Ph.D. thesis, PSL University (2024).

[57] J. H. Page, H. P. Schriemer, A. E. Bailey, and D. A. Weitz, Experimental test of the diffusion approximation for multiply scattered sound, Phys. Rev. E **52**, 3106 (1995).

[58] X. Jia, Codalike multiple scattering of elastic waves in dense granular media, Phys. Rev. Lett. **93**, 154303 (2004).