# Reflection Measurement of the Scattering Mean Free Path at the Onset of Multiple Scattering

Antton Goïcoechea, Cécile Brütt, Arthur Le Ber, Flavien Bureau, William Lambert, Claire Prada, Arnaud Derode, Alexandre Aubry

► To cite this version:

Antton Goïcoechea, Cécile Brütt, Arthur Le Ber, Flavien Bureau, William Lambert, et al. Reflection Measurement of the Scattering Mean Free Path at the Onset of Multiple Scattering. 2024. hal-04372441

HAL Id: hal-04372441

https://hal.sorbonne-universite.fr/hal-04372441

Preprint submitted on 4 Jan 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Reflection Measurement of the Scattering Mean Free Path at the Onset of Multiple Scattering

Antton Goïcoechea, Cécile Brütt, Arthur Le Ber, Flavien Bureau, William Lambert, Claire Prada, Arnaud Derode, Alexandre Aubry

► To cite this version:

Antton Goïcoechea, Cécile Brütt, Arthur Le Ber, Flavien Bureau, William Lambert, et al. Reflection Measurement of the Scattering Mean Free Path at the Onset of Multiple Scattering. Physical Review Letters, 2024, 133 (17), pp.176301. 10.1103/PhysRevLett.133.176301. hal-04372441v2

HAL Id: hal-04372441

https://hal.sorbonne-universite.fr/hal-04372441v2

Submitted on 11 Sep 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Reflection Measurement of the Scattering Mean Free Path
at the Onset of Multiple Scattering

Antton Goïcoechea,¹,* Cécile Brütt,² Arthur Le Ber,¹ Flavien Bureau,¹

William Lambert,³ Claire Prada,¹ Arnaud Derode,¹ and Alexandre Aubry¹, †

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

²Safran Tech, Digital Sciences and Technologies Department,
78114 Magny-Les-Hameaux, France

³SuperSonic Imagine, 13290 Aix-en-Provence, France

(Dated: January 4, 2024)
Reflection Measurement of the Scattering Mean Free Path
at the Onset of Multiple Scattering

Antton Goïcoechea,¹,* Cécile Brütt,² Arthur Le Ber,¹ Flavien Bureau,¹

William Lambert,³ Claire Prada,¹ Arnaud Derode,¹ and Alexandre Aubry¹, †

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

²Safran Tech, Digital Sciences and Technologies Department,
78114 Magny-Les-Hameaux, France

³SuperSonic Imagine, 13290 Aix-en-Provence, France

(Dated: September 11, 2024)
Abstract

Multiple scattering of waves presents challenges for imaging complex media but offers potential
for their characterization. Its onset is actually governed by the scattering mean free path $\ell_s$ that
provides crucial information on the medium micro-architecture. Here, we introduce a reflection
matrix method designed to estimate this parameter from the time decay of the single scattering
rate. Our method is first validated by means of an ultrasound experiment on a tissue-mimicking
phantom. An in-vivo measurement of $\ell_s$ is then performed in the human liver. This study thus
opens important perspectives for quantitative imaging whether it be for non-destructive testing or
biomedical applications.

Multiple scattering (MS) of waves proves to be a captivating subject manifesting itself
across all the spectrum of wave physics [1–15]. In an inhomogeneous medium, a common
approach is to consider a scattering sample as one realization of a random process. Within
this paradigm, the relevant parameter for characterizing wave propagation within an hetero-
geneous medium is the scattering mean free path, denoted as $\ell_s$. This parameter represents
the typical distance between successive scattering events. For a time-of-flight $t$ smaller than
the corresponding mean free time, $\tau_s = \ell_s/c$ (with $c$, the wave velocity), wave propaga-
tion behaves akin to a homogeneous medium, displaying a ballistic trajectory. However,
as the time of flight increases, scattering events progressively randomize the direction of
wave propagation. The trajectory of the wave can be described as a random walk, and
energy transport finds an apt model in diffusion theory for $t \gg \tau_s$. The significance of $\ell_s$
extends beyond its role as a fundamental quantity dictating the onset of MS; it also holds
paramount importance for characterization purposes. Indeed, $\ell_s$ directly correlates with the
micro-architecture of the medium i.e., the scatterer concentration $n$ and their scattering
cross-section in discrete media [16, 17], or the correlation function in models involving a
continuous random disorder [18, 19]. $\ell_s$ has thus emerged as a promising marker in diverse
applications such as biomedical diagnosis using ultrasound [20–25] or light [26]. Further-
more, it serves as a crucial monitoring parameter in nondestructive testing [27–32] and
seismology [33–36].

From a theoretical standpoint, the probability density function $I(\mathbf{r}, t)$ for the travel time $t$ in a random medium, with $\mathbf{r}$ the relative position between the source and the receiver, can be obtained by solving the radiative transfer equation. From an experimental standpoint,
Abstract

Multiple scattering of waves presents challenges for imaging complex media but offers potential
for their characterization. Its onset is actually governed by the scattering mean free path $\ell_s$ that
provides crucial information on the medium micro-architecture. Here, we introduce a reflection
matrix method designed to estimate this parameter from the time decay of the single scattering rate.
Our method is first validated by an ultrasound experiment on a tissue-mimicking phantom before
being applied in-vivo to a human liver. This study opens important perspectives for quantitative
imaging of heterogeneous media with waves, whether it be for non-destructive testing, biomedical
or geophysical applications.

Multiple scattering (MSc) of waves proves to be a captivating subject manifesting itself
across all the spectrum of wave physics [1–15]. In an inhomogeneous medium, a common
approach is to consider a scattering sample as one realization of a random process. Within
this paradigm, the relevant parameter for characterizing wave propagation within an hetero-
geneous medium is the scattering mean free path, denoted as $\ell_s$. This parameter represents
the typical distance between successive scattering events. For a time-of-flight $t$ smaller than
the corresponding mean free time, $\tau_s = \ell_s/c$ (with $c$, the wave velocity), wave propaga-
tion behaves akin to a homogeneous medium, displaying a ballistic trajectory. However,
as the time of flight increases, scattering events progressively randomize the direction of
wave propagation. The trajectory of the wave can be described as a random walk, and
energy transport finds an apt model in diffusion theory for $t \gg \tau_s$. The significance of $\ell_s$
extends beyond its role as a fundamental quantity dictating the onset of MSc; it also holds
paramount importance for characterization purposes. Indeed, $\ell_s$ directly correlates with
short-scale fluctuations of the wave speed, thereby allowing for instance: ultrasound quan-
tification of fat in liver [16–18], osteoporosis diagnosis [19, 20] or assessment of edema in the
lung [21]; optical mapping of oxy- and deoxy-hemoglobin concentration in neuroscience [22]
or cornea transparency in optical coherence tomography (OCT) [23]; characterization of
the subsoil for seismic hazard assessment [24, 25] or of polycrystalline structures for non
destructive testing [26, 27]. This list is, by no means, exhaustive but illustrates the large
variety of fields in which a measurement of $\ell_s$ is crucial.

From a theoretical standpoint, the probability density function $I(\mathbf{r}, t)$ for the travel time $t$ in a random medium, with $\mathbf{r}$ the relative position between the source and the receiver, can
<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) Experimental configuration: An array of transducers is placed in front of a 2D semi-infinite random medium [$\ell_s = 100$ mm, $\ell_a = 100$ m, $c=1540$ m.s$^{-1}$]. (b) Reflected intensity vs. $t/\tau_s$: The single scattering intensity (green dots) obtained via a Monte Carlo simulation of the raditaive transfer equation [40] is compared with the analytical prediction of Eq. (2); the overall intensity obtained via Monte Carlo (blue dots) is compared with its short-time [Eq. (4)] and long-time [Eq. (3)] analytical expressions. The x- and y-axes are in log-scale. (c) Single scattering rate $\rho_s(t)$: The numerical result (black dots) is compared with its short-time expression [Eq. (5)] and the diffusion model prediction [Eq. (3)]. The x- and y-axes are in linear and log-scale, respectively.</figcaption></figure>

the averaged instantaneous intensity of the wave-field can be measured and serves as an estimate of $I(\mathbf{r}, t)$. In the following, the word “intensity” will be used in the sense of a travel-time distribution. Inverting experimental measurements of $I(\mathbf{r}, t)$ allows for the derivation of a spatial mapping of $\ell_s$ along with other relevant transport parameters. This methodology aligns with the concept of optical diffuse tomography [26, 37, 38]. However, the practical implementation of this approach is encumbered by a substantial computational load, and its spatial resolution is poor, as it scales with the imaging depth. An alternative experimental strategy focuses on the investigation of the wave-field itself, particularly the exponential attenuation of its ballistic component in a transmission configuration [39]. The characteristic length scale associated with this attenuation is the extinction length, denoted as $\ell_{\text{ext}}$, which encapsulates both scattering and absorption losses: $\ell_{\text{ext}} = (\ell_s^{-1} + \ell_a^{-1})^{-1}$, where $\ell_a$ represents the absorption length. Hence, the transmitted wave-field can not furnish an independent measurement of scattering and absorption. Additionally, the impracticality of a transmission measurement is often evident, as only one side of the medium is typically accessible for most applications.
<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) Experimental configuration: An array of transducers is placed in front of a 2D semi-infinite random medium [$\ell_s = 100$ mm, $\ell_a = 100$ m, $c=1540$ m.s$^{-1}$]. (b) Reflected intensity vs. $t/\tau_s$: The single scattering intensity (green dots) obtained via a Monte Carlo simulation of the radiative transfer equation [32] is compared with the analytical prediction of Eq. (2); the overall intensity obtained via Monte Carlo (blue dots) is compared with its short-time [Eq. (4)] and long-time [Eq. (3)] analytical expressions. (c) Single scattering rate $\rho_s(t)$: The numerical result (black dots) is compared with its short-time expression [Eq. (5)] and the diffusion model prediction [Eq. (3)]. The x- and y-axes are in linear and log-scale, respectively.</figcaption></figure>

be obtained by solving the radiative transfer equation. From an experimental standpoint, the averaged instantaneous intensity of the wave-field can be measured and serves as an estimate of $I(\mathbf{r}, t)$. In the following, the word “intensity” will be used in the sense of a travel-time distribution. Inverting experimental measurements of $I(\mathbf{r}, t)$ allows for the derivation of a spatial mapping of $\ell_s$ along with other relevant transport parameters. This methodology aligns with the concept of optical diffuse tomography [28–30]. However, the practical implementation of this approach is encumbered by a substantial computational load, and its spatial resolution is poor, as it scales with the imaging depth. An alternative experimental strategy focuses on the investigation of the wave-field itself, particularly the exponential attenuation of its ballistic component in a transmission configuration [31]. The characteristic length scale associated with this attenuation is the extinction length, denoted as $\ell_{\text{ext}}$, which encapsulates both scattering and absorption losses: $\ell_{\text{ext}} = (\ell_s^{-1} + \ell_a^{-1})^{-1}$, where $\ell_a$ represents the absorption length. Hence, the transmitted wave-field can not furnish an independent measurement of scattering and absorption. Additionally, the impracticality of a transmission measurement is often evident, as only one side of the medium is typically accessible for most applications.
A local measurement of $\ell_s$, independent from absorption losses, would be extremely rewarding in reflection, not only for practical reasons, but also for quantitative purposes. In this paper, we demonstrate all these capabilities in the context of 2D ultrasound imaging. Our approach is based on a local discrimination between single scattering (SS) and MS. Assuming isotropic scattering and solving radiative transfer equation [41], we first show that the SS rate $\rho_s(t)$ scales as $\exp[-3ct/(4\ell_s)]$ for $t < \tau_s$. This observable can therefore give access to a direct measurement of $\ell_s$. Experimentally, $\rho_s(t)$ is obtained by investigating the reflection matrix of the medium in a focused basis [42–44] and projecting it onto a characteristic SS space [45]. This method is first validated on a phantom generating an ultrasound speckle characteristic of tissues. We provide local measurements of $\ell_s$ and $\ell_a$ in areas showing different scattering properties. We then apply our approach to in-vivo ultrasound data acquired in a liver. In particular, we will show the robustness of our method with respect to the variation of reflectivity, in contrast with existing methods relying on attenuation measurements [23, 46].

Our approach applies to the experimental configuration depicted in Fig. 1(a). An array of $N$ transducers is placed in front of a scattering sample. These transducers are 10 mm in height, which is much larger than the average wavelength and a vertical cylindrical acoustic lens ensures that the emitted beam remains collimated in the $(x, z)$ plane. When one element at transverse position $\mathbf{u}_{\text{in}}$ emits an ultrasound pulse, it generates an incident cylindrical wave, which is scattered by the heterogeneities of the medium. Similarly, in reception only waves propagating in the $(x, z)$ plane are recorded by the transducers. The reflected wave-field $R(\mathbf{u}_{\text{out}}, \mathbf{u}_{\text{in}}, t)$ is recorded by each transducer identified by its position $\mathbf{u}_{\text{out}}$. By repeating this operation for each element as a source, each reflected wave-field can be stored into a canonical reflection matrix $\mathbf{R}_{\mathbf{u}\mathbf{u}}(t) = R(\mathbf{u}_{\text{out}}, \mathbf{u}_{\text{in}}, t)$. Each recorded signal exhibits a complex signature that results from a random superposition of partial waves, each one being associated with a different scattering path [Fig. 1(a)]. A classical approach is to consider this scattering sample as one realization of a random process, and study the overall reflected intensity:

$$I(t) = \text{Tr} \left[ \mathbf{R}(t) \mathbf{R}^{\dagger}(t) \right] / N, \quad (1)$$

where the symbols Tr and $\dagger$ stands for matrix trace and transpose conjugate, respectively. The back-scattered intensity $I(t)$ is made of two contributions: (i) a SS component $I_S$ in which the incident wave undergoes only one scattering event before coming back to the
A local measurement of $\ell_s$, independent from absorption losses, would be extremely re-
warding in reflection, not only for practical reasons, but also for quantitative purposes.
In this paper, we demonstrate all these capabilities in the universal framework of matrix
imaging. Based on the recording of the reflection matrix associated with a sensor network,
it has already found applications in many fields of wave physics ranging from optical mi-
croscopy [33, 34] or ultrasound [35, 36] to seismology [37, 38] or microwaves [39]. In this
paper, the proof-of-concept is performed with ultrasound but the proposed method is very
general and can be applied to any kind of wave and sensor network, provided that the spatial
sampling of the reflected wave-field satisfies the Nyquist criterion. Our method is based on a
local discrimination between single scattering (SSc) and MSc. Assuming isotropic scattering
and solving radiative transfer equation [40], we first show that the SSc rate $\rho_s(t)$ scales as
$\exp[-3ct/(4\ell_s)]$ for $t < \tau_s$, i.e for predominant single scattering ($\rho_s > 0.5$). This observable
can therefore give access to a direct measurement of $\ell_s$. Experimentally, $\rho_s(t)$ is obtained
by investigating the reflection matrix of the medium in a focused basis [35, 36, 41] and pro-
jecting it onto a characteristic SSc space [42]. The method is first validated on a phantom
generating a speckle mimicking the response of actual tissues to ultrasound. We provide
local measurements of $\ell_s$ and $\ell_a$ in areas showing different scattering properties. We then
apply our approach to in-vivo ultrasound data acquired in a liver. In particular, we will
show the robustness of our method with respect to the variation of reflectivity, in contrast
with existing methods relying on attenuation measurements [17, 43].

Our approach applies to the experimental configuration depicted in Fig. 1(a). An array
of N transducers is placed in front of a scattering sample. These transducers are 10 mm in
height, which is much larger than the average wavelength and a vertical cylindrical acoustic
lens ensures that the emitted beam remains collimated in the (x, z) plane. When one element
at transverse position **u**<sub>in</sub> emits an ultrasound pulse, it generates an incident cylindrical wave,
which is scattered by the heterogeneities of the medium. Similarly, in reception only waves
propagating in the (x, z) plane are recorded by the transducers. The reflected wave-field
R(**u**<sub>out</sub>, **u**<sub>in</sub>, t) is recorded by each transducer identified by its position **u**<sub>out</sub>. By repeating
this operation for each element as a source, each reflected wave-field can be stored into a
canonical reflection matrix **R**<sub>uu</sub>(t) = R(**u**<sub>out</sub>, **u**<sub>in</sub>, t).

Each recorded signal exhibits a complex signature that results from a random superposi-
tion of partial waves, each one being associated with a different scattering path [Fig. 1(a)].
sensors [green arrows in Fig. 1(a)]; (ii) a MS component $I_M$ that predominates for $t \gg \tau_s$ [blue arrows in Fig. 1(a)]. To assess the relative weight between each component, radiative transfer equation shall be considered. In a seminal paper, Paasschens [41] solved this equation in a 2D geometry assuming isotropic scattering. In particular, he derived an analytical expression for each scattering order of the mean intensity Green's function. Based on this decomposition, a theoretical expression of $I_S$ can be derived for a semi-infinite medium [40]:

$$I_S(t) = \exp(-ct/\ell_{\text{ext}})/(4\ell_s). \qquad (2)$$

Not surprisingly, $I_S(t)$ displays an exponential attenuation dictated by $\ell_{\text{ext}}$. As shown by Fig. 1(b), this analytical result is in excellent agreement with the time-of-flight distribution of singly-scattered echoes obtained by means of a Monte Carlo simulation of the radiative transfer equation in a semi-infinite random medium [40]. As to the overall intensity, the medium interface can be taken into account in the diffusive approximation [47]. A power law scaling of $I(t)$ as $t^{-3/2}$ is predicted [40]:

$$I(t) \sim \frac{\exp(-ct/\ell_a) z_0}{\sqrt{\pi D t} ct} \qquad (3)$$

with $z_0 = 2\ell_s/3$, the extrapolation length [48] and $D$, the diffusion coefficient. In a 2D geometry and for isotropic scattering, $D = c\ell_s/2$. As displayed by Fig. 1(b), this diffusive result predicts well the time dependence of the reflected intensity in the long-time limit but it does not grasp the onset of MS at short times-of-flight. Interestingly, the radiative transfer solution can provide the asymptotic behavior of $I(t)$ for $t < \tau_s$ [40]:

$$I(t) \sim \exp(-ct/4\ell_s) \exp(-ct/\ell_a)/(4\ell_s). \qquad (4)$$

This solution is shown to fit perfectly the time-of-flight distribution of reflected waves obtained by the Monte Carlo simulation at short times-of-flight [Fig. 1(b)]. MS has thus a strong impact even at small times-of-flight since it modifies the exponential scaling of $I(t)$ with respect to the SS component [Eq. (2)]. This striking property can be highlighted by investigating the SS rate, $\rho_s(t) = I_S(t)/I(t)$. The ratio of Eqs. (2) and (4) shows that $\rho_s$ displays an exponential attenuation for $t < \tau_s$ that only depends on $\ell_s$:

$$\rho_s(t) \sim \exp(-3ct/4\ell_s). \qquad (5)$$

Fig. 1(c) shows the validity of the last expression for $t < \tau_s$. This fundamental result proves analytically that a discrimination between SS and MS can lead to a measurement
A classical approach is to consider this scattering sample as one realization of a random process, and study the overall reflected intensity:

$$I(t) = \mathrm{Tr} \left[ \mathbf{R}(t) \mathbf{R}^{\dagger}(t) \right] / N, \qquad (1)$$

where the symbols Tr and † stands for matrix trace and transpose conjugate, respectively. The back-scattered intensity $I(t)$ is made of two contributions: (i) a SSc component $I_S$ in which the incident wave undergoes only one scattering event before coming back to the sensors [green arrows in Fig. 1(a)]; (ii) a MSc component $I_M$ that predominates for $t \gg \tau_s$ [blue arrows in Fig. 1(a)]. To assess the relative weight between each component, radiative transfer equation shall be considered. In a seminal paper, Paasschens [40] solved this equation in a 2D geometry assuming isotropic scattering. In particular, he derived an analytical expression for each scattering order of the mean intensity Green's function. Based on this decomposition, a theoretical expression of $I_S$ can be derived for a semi-infinite medium [32]:

$$I_S(t) = \exp(-ct/\ell_{\mathrm{ext}}) / (4\ell_s). \qquad (2)$$

Not surprisingly, $I_S(t)$ displays an exponential attenuation dictated by $\ell_{\mathrm{ext}}$. As shown by Fig. 1(b), this analytical result is in excellent agreement with the time-of-flight distribution of singly-scattered echoes obtained by means of a Monte Carlo simulation of the radiative transfer equation in a semi-infinite random medium [32]. As to the overall intensity, the medium interface can be taken into account in the diffusive approximation [44]. A power law scaling of $I(t)$ as $t^{-3/2}$ is predicted [32]:

$$I(t) \sim \frac{\exp(-ct/\ell_a) z_0}{\sqrt{\pi Dt} ct} \qquad (3)$$

with $z_0 = 2\ell_s/3$, the extrapolation length [45] and $D$, the diffusion coefficient. In a 2D geometry and for isotropic scattering, $D = c\ell_s/2$. As displayed by Fig. 1(b), this diffusive result predicts well the time dependence of the reflected intensity in the long-time limit but it does not grasp the onset of MSc at short times-of-flight. Interestingly, the radiative transfer solution can provide the asymptotic behavior of $I(t)$ for $t < \tau_s$ [32]:

$$I(t) \sim \exp(-ct/4\ell_s) \exp(-ct/\ell_a) / (4\ell_s). \qquad (4)$$

This solution is shown to fit perfectly the time-of-flight distribution of reflected waves obtained by the Monte Carlo simulation at short times-of-flight [Fig. 1(b)]. MSc has thus a
<figure><img src="image_3.png" /><figcaption>FIG. 2. Principle of the SS filter. (a) The broadband focused reflection matrix is obtained by focusing in input (purple) and output (orange) at the same depth, yielding $\mathbf{R}_{xx}(z)$. (b) An example of $\mathbf{R}_{xx}(z)$ is displayed at depth $z = 45$ mm in the tissue-mimicking phantom. (c) $\mathbf{R}_{xx}(z)$ is projected onto a set of orthonormal matrices forming the so-called SS space [45]. (d) The result of this projection is a SS matrix $\mathbf{S}_{xx}(z)$ from which MS has been discarded.</figcaption></figure>

of $\ell_s$ independent of $\ell_a$ [22]. However, such a measurement is integrated over the whole lateral extension of the sample. If one wants a local measurement of transport parameters, a focusing process is required.

As a first experimental proof-of-concept, we consider the case of a tissue-mimicking phantom that displays an homogeneous speed of sound $c_0 = 1540 \text{ mm}/\mu\text{s}$ and attenuation $\alpha \simeq 0.7 \text{ dB.cm}^{-1}$.MHz$^{-1}$ but also exhibits variations in terms of speckle reflectivity [Fig. 3]. The reflection matrix is acquired using an array of 256 transducers with an inter-element distance of 0.2 mm operating in the 5-10 MHz frequency bandwidth. The first step of the method consists in applying a focused beamforming process to the recorded reflection matrix at input and output [40, 42]. The result is a confocal image that is a satisfying estimator of the sample reflectivity under SS assumption [Fig. 3(a)]. This confocal image shows that the phantom contains two cylinders of higher reflectivity than the surrounding speckle [see red and blue dashed lines in Fig. 3(a)]. However, note that a residual MS component can still pollute the confocal image since there is no such a thing as pure SS medium.
<figure><img src="image_3.png" /><figcaption>FIG. 2. (a) Experimental configuration and ultrasound image of the tissue-mimicking phantom. (b,c) $\hat{I}_S$ (b) and $\hat{\rho}_s$ (c) versus time $t = 2z/c$ measured along the blue and red lines displayed in (a). The y-axes are in log-scale. The experimental measurements are fitted by Eqs. (2) and (5) (dashed lines) in the shallow (25-36 mm) and large (36-43 mm) depth range with values of $\ell_s$ and $\ell_{ext}$ reported in Tab. I.</figcaption></figure>

strong impact even at small times-of-flight since it modifies the exponential scaling of $I(t)$ with respect to the SSc component [Eq. (2)]. This striking property can be highlighted by investigating the SSc rate, $\rho_s(t) = I_S(t)/I(t)$. The ratio of Eqs. (2) and (4) shows that $\rho_s$ displays an exponential attenuation for $t < \tau_s$ that only depends on $\ell_s$:

$$ \rho_s(t) \sim \exp(-3ct/4\ell_s). \qquad (5) $$

Fig. 1(c) illustrates the validity of Eq. 5 for $t < \tau_s$. This fundamental result proves analytically that a discrimination between SSc and MSc can lead to a measurement of $\ell_s$ independent of $\ell_a$ [46]. However, results in Fig. 1 are not spatially resolved along the x-axis. A focusing process is required to achieve local measurements.

As a first experimental proof-of-concept, we consider the case of a tissue-mimicking phantom (CIRS, Model 054GS) that displays an homogeneous speed of sound $c_0 = 1540 \text{ mm/}\mu\text{s}$ and attenuation $\alpha \simeq 0.7 \text{ dB.cm}^{-1}\text{.MHz}^{-1}$ but also exhibits variations in terms of speckle reflectivity [Fig. 2]. The reflection matrix is acquired using an array of 256 transducers with an inter-element distance of 0.2 mm operating in the 5-10 MHz frequency bandwidth.
To estimate the weight of MS, a focused reflection matrix $\mathbf{R}_{xx}(z)$ can be synthesized at each depth $z$ [42]. Its coefficients $R(x_{\text{out}}, x_{\text{in}}, z)$ are the response of the medium at the ballistic time $t = 2z/c$ between virtual sources and detectors located at $(x_{\text{in}}, z)$ and $(x_{\text{out}}, z)$ [see Fig. 2(a)]. It displays two contributions: (i) A SS component that emerges along the diagonal or close-diagonal elements; (ii) A MS component that gives rise to a diffuse halo that spreads over off-diagonal elements. To separate both components, an adapted matrix filter has been recently proposed [45]. It consists in a projection of $\mathbf{R}_{xx}(z)$ onto a set of matrices characteristic of SS [40] [Fig. 2(c)]. This operation gets rid of the off-diagonal diffuse halo and returns the SS matrix $\mathbf{S}_{xx}(z)$ [Fig. 2(d)].

The norm of $\mathbf{S}_{xx}(z)$ gives access to an estimator of the mean SS intensity at each depth $z$, or equivalently at each time-of-flight $t = 2z/c$:

$$ \hat{I}_S(t = 2z/c) = \text{Tr} \left[ \mathbf{S}_{xx}(z) \mathbf{S}_{xx}^{\dagger}(z) \right] / N. \qquad (6) $$

The single scattering intensity is estimated along the blue and red lines of Fig. 3(a) by placing the probe along each cylinder axis. As expected, $\hat{I}_S$ displays a similar attenuation decay in each region of the phantom whatever the local reflectivity. At shallow depth, the measured values for $\ell_{\text{ext}}$ [see Tab. I] are in good agreement with the manufacturer value at the central frequency of 7.5 MHz ($\ell_{\text{ext}} \sim 19$ mm). In the red cylinder, a slight increase is observed for the measured $\ell_{\text{ext}}$ but the error bar on this estimation is large. In the blue cylinder, large speckle fluctuations prevent from a reliable fitting process of $\hat{I}_S$.

To go beyond a measurement of $\ell_{\text{ext}}$, one can compute estimated the SS rate: $\hat{\rho}_s(t) = \hat{I}_S(t)/I(t)$. Interestingly, the SS rate seems more robust with respect to the medium reflectivity fluctuations [Fig. 3(c)]. $\hat{\rho}_s$ exhibits a very different behaviour in the blue and red areas [Fig. 3(d)]. While $\ell_s$ is much smaller in the blue cylinder than in the phantom tissue upstream, it is of the same order of magnitude in the red cylinder. The slight increase in contrast exhibited by the red cylinder on the ultrasound image is therefore probably due a scattering function peaked in the backward direction. Such phenomena arise in the Rayleigh regime when disorder is induced by density fluctuations, in contrast with compressibility variations [19]. This proof-of-concept experiment thus demonstrates our ability to provide some relevant information about the nature of disorder.

Matrix imaging and its SS filter can therefore provide quantitative markers useful for medical diagnosis. Ultrasound attenuation is, for instance, already estimated in liver to
<figure><img src="image_4.png" /><figcaption>FIG. 3. Principle of the SSc filter. (a) The broadband focused reflection matrix is obtained by focusing in input (purple) and output (orange) at the same depth, yielding $\mathbf{R}_{xx}(z)$. (b) An example of $\mathbf{R}_{xx}(z)$ is displayed at depth $z = 45$ mm in the tissue-mimicking phantom. (c) $\mathbf{R}_{xx}(z)$ is projected onto a set of orthonormal matrices forming the so-called SSc space [42]. (d) The result of this projection is a SSc matrix $\mathbf{S}_{xx}(z)$ from which MSc has been discarded.</figcaption></figure>

The first step of the method consists in applying a focused beamforming process to the recorded reflection matrix at input and output [32, 35]. A prior optimization of the wave velocity model is performed to limit aberrations [35]. The result is a confocal image that is a satisfying estimator of the sample reflectivity under SSc assumption [Fig. 2(a)]. This confocal image shows that the phantom contains two cylinders of higher reflectivity than the surrounding speckle [see red and blue dashed lines in Fig. 2(a)]. However, note that a residual MSc component can still pollute the confocal image.

To estimate the weight of MSc, a focused reflection matrix $\mathbf{R}_{xx}(z)$ can be synthesized at each depth $z$ [35]. Its coefficients $R(x_{\text{out}}, x_{\text{in}}, z)$ are the response of the medium at the ballistic time $t = 2z/c$ between virtual sources and detectors located at $(x_{\text{in}}, z)$ and $(x_{\text{out}}, z)$ [see Fig. 3(a)]. It displays two contributions: (i) A SSc component that emerges along the diagonal or close-diagonal elements; (ii) A MSc component that gives rise to a diffuse halo that spreads over off-diagonal elements. To separate both components, an adapted matrix filter has been recently proposed [42]. It consists in a projection of $\mathbf{R}_{xx}(z)$ onto a set
<figure><img src="image_4.png" /><figcaption>FIG. 3. (a) Experimental configuration and ultrasound image of the tissue-mimicking phantom. (b,c) $\hat{I}_S$ (b) and $\hat{\rho}_s$ (c) versus time $t = 2z/c$ measured along the blue and red lines displayed in (a). The y-axes are in log-scale. The experimental measurement are fitted by Eqs. (5) and (2) (dashed lines) in the shallow (25-36 mm) and large (36-43 mm) depth range with values of $\ell_s$ and $\ell_{ext}$ reported in Tab. I.</figcaption></figure>

<table><caption>TABLE I. Values of $\ell_{ext}$ and $\ell_s$ measured along the blue and red lines of the phantom displayed in Fig. 3.</caption><thead><tr><th>Depth range</th><th>l<sub>ext</sub></th><th>l<sub>s</sub></th></tr></thead><tbody><tr><td>25-36 mm (blue)</td><td>17 ± 3 mm</td><td>247 ± 90 mm</td></tr><tr><td>36-43 mm (blue)</td><td>—</td><td>36 ± 7 mm</td></tr><tr><td>25-36 mm (red)</td><td>15 ± 2 mm</td><td>90 ± 11 mm</td></tr><tr><td>36-43 mm (red)</td><td>23 ± 8 mm</td><td>90 ± 11 mm</td></tr></tbody></table>

monitor diseases such as steatosis [23, 49, 50]. However, this measurement accumulates absorption and scattering losses and a discrimination between these two mechanisms would probably lead to a sharper diagnosis. This disease consists in an accumulation of fat droplets in liver and leads to a drastic increase of scattering [51]. As a first step towards such biomedical applications, we apply our approach to in-vivo ultrasound data acquired on the liver
of matrices characteristic of SSc [32] [Fig. 3(c)]. This operation gets rid of the off-diagonal diffuse halo and returns the SSc matrix $\mathbf{S}_{xx}(z)$ [Fig. 3(d)]. The norm of $\mathbf{S}_{xx}(z)$ gives access to an estimator of the mean SSc intensity at each depth $z$, or equivalently at each time-of-flight $t = 2z/c$:

$$ \hat{I}_S(t = 2z/c) = \mathrm{Tr} \left[ \mathbf{S}_{xx}(z) \mathbf{S}_{xx}^{\dagger}(z) \right] / N. \qquad (6) $$

The single scattering intensity is estimated along the blue and red lines of Fig. 2(a) by placing the probe along each cylinder axis. As expected, $\hat{I}_S$ displays a similar attenuation decay in each region of the phantom whatever the local reflectivity. At shallow depth, the measured values for $\ell_{ext}$ [see Tab. I] are in good agreement with the manufacturer value at the central frequency of 7.5 MHz ($\ell_{ext} \sim 19$ mm). In the cylinders, the measured $\ell_{ext}$ is increased by nearly a factor three but the error bar on this estimation is important due to large speckle fluctuations.

To go beyond a measurement of $\ell_{ext}$, one can compute the SSc ratio: $\hat{\rho}_s(t) = \hat{I}_S(t)/I(t)$. Interestingly, the SSc rate seems more robust with respect to the medium reflectivity fluctuations [Fig. 2(c)] and reduces the uncertainty on the estimation of $\ell_s$ compared with $\ell_{ext}$. $\hat{\rho}_s$ exhibits a very different behaviour in the blue and red areas [Fig. 2(c)]. In the blue cylinder, $\ell_s$ is much smaller than its value in the phantom tissue upstream. Not surprisingly, the high contrast exhibited by this cylinder is therefore explained by a stronger scattering. In the red cylinder, $\ell_s$ is the same as in the tissue upstream and almost two orders of magnitude larger than $\ell_{ext}$. Absorption is therefore the main mechanism accounting for attenuation and the slightly larger contrast exhibited by the red cylinder in Fig. 2(a) is explained by a weaker absorption compared to surrounding speckle. This proof-of-concept experiment thus demonstrates our ability to provide a local discrimination between scattering and absorption losses.

Matrix imaging and its SSc filter can therefore provide quantitative markers useful for medical diagnosis. As a first step towards such bio-medical applications, we apply our approach to in-vivo ultrasound data acquired on the liver of a healthy human subject [Fig. 4(a)]. The experimental conditions are the same as for the phantom experiment. The result of our SSc filter is presented in Fig. 4(b) for the region surrounded by a white rectangle in Fig. 4(a). As previously observed in the phantom, $\hat{\rho}_s$ exhibits much less fluctuations compared to $I_S$. The latter quantity is extremely sensitive to reflectivity variations caused by the presence of structures like veins and cannot provide a reliable measurement of $\ell_{ext}$. On
<figure><img src="image_5.png" /><figcaption>FIG. 4. (a) Experimental configuration and ultrasound image of the liver, including the muscle and fat layers at shallow depths. The area of interest is surrounded by a white rectangle in which the gray scale spans from -25 to 10 dB to improve the image contrast. (b) $\hat{\rho}_s$ (blue dots) and $\hat{I}_S$ (green dots) vs. depth $z = ct/2$. The y-axes are in log-scale. The fit of $\rho_s(t)$ with Eq. (5) (blue line) leads to an estimation of $\ell_s$.</figcaption></figure>

of a healthy human subject [Fig. 4(a)]. The experimental conditions are the same as for the phantom experiment. The result of our SS filter is presented in Fig. 4(b) for the region surrounded by a white rectangle in Fig. 4(a). As previously observed in the phantom, $\hat{\rho}_s$ exhibits much less fluctuations compared to $I_S$. The latter quantity is extremely sensitive to reflectivity variations caused by the presence of structures like veins and cannot provide a reliable measurement of $\ell_{ext}$. On the contrary, the fitting of $\hat{\rho}_s$ with Eq. (5) leads to an estimation of $\ell_s \simeq 22 \pm 2$ mm. This value is remarkably low but is in agreement with the important MS contribution observed for the same organ in a previous study [42]. Of course, this demonstration is not limited to liver and can be extended to any tissue giving rise to ultrasound speckle.

Even though the proposed method is robust with respect to speckle variations, it shows some limitations. First, the SS filter is only efficient when the MS background is not too large. The tomography of $\ell_s$ is therefore limited to a penetration depth of a few $\ell_s$. Another issue can arise when wave propagates from a strongly scattering region to a weakly scattering area. MS induced by the first layer predominates and hides the SS component of the second one [40]. A more elaborate inversion procedure of $\rho_s$ is thus needed and will be the object
<table><caption>TABLE I. Values of $\ell_{ext}$ and $\ell_s$ measured along the blue and red lines of the phantom displayed in Fig. 2.</caption><thead><tr><th>Depth range</th><th>l<sub>ext</sub></th><th>l<sub>s</sub></th></tr></thead><tbody><tr><td>25-36 mm (red)</td><td>17 ± 2 mm</td><td rowspan="2">1280 ± 150 mm</td></tr><tr><td>36-43 mm (red)</td><td>39 ± 11 mm</td></tr><tr><td>25-36 mm (blue)</td><td>19 ± 1 mm</td><td rowspan="2">89 ± 6 mm</td></tr><tr><td>36-43 mm (blue)</td><td>49 ± 15 mm</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>FIG. 4. (a) Experimental configuration and ultrasound image of the liver, including the muscle and fat layers at shallow depths, with an effective aperture scaling as the imaging depth [32]. The area of interest is surrounded by a white rectangle in which the gray scale spans from -25 to 10 dB to improve the image contrast. (b) $\hat{\rho}_s$ (blue dots) and $\hat{I}_S$ (green dots) vs. depth $z = ct/2$. The y-axes are in log-scale. The fit of $\rho_s(t)$ with Eq. (5) (blue line) leads to an estimation of $\ell_s$.</figcaption></figure>

the contrary, the fitting of $\hat{\rho}_s$ with Eq. (5) leads to an estimation of $\ell_s \simeq 44 \pm 4$ mm. This value is remarkably low but is in agreement with the important MSc contribution observed for the same organ in a previous study [35]. Of course, this demonstration is not limited to liver nor to ultrasound and can be extended to any tissue giving rise to speckle. As a proof for this generality, a recent study has shown how the proposed method can be implemented in OCT for measuring $\ell_s$ in an opaque cornea [47].

Even though the proposed method is robust with respect to speckle variations, it shows
of a future study. Another relevant question concerns the applicability of this method to a
2D probe, *i.e* a 3D imaging geometry. Interestingly, an exponential decay is still observed
at short times-of-flight in 3D: ρs(t) ~ exp(-0.55ct/ls) [40]. Hence the proposed method
can be extended to 3D configurations. A last perspective for this work is its extension to
anisotropic scattering. This achievement is particularly crucial for optical microscopy since
light scattering is sharply peaked in the forward direction in biological tissues.

In summary, we have introduced a novel methodology for local characterization of scat-
tering media, employing a non-invasive reflection setup. Our approach successfully discrim-
inates between SS and MS events, providing an independent measurement of l_s and l_a.
This obviates the need for traditional transmission experiments, offering a solution to a
longstanding challenge in disordered media research. Beyond its inherent advantages, our
method stands as a potent tool for quantitative imaging, applicable to in-vivo or in-situ sce-
narios. Notably, the robustness of the SS rate in the face of sample reflectivity fluctuations
sets our approach apart from techniques reliant on the time-dependence of backscattered
echoes for attenuation measurements. Importantly, while our proof-of-concept experiments
were conducted with ultrasound, the universality of matrix imaging renders it applicable
to diverse fields of wave physics wherever multi-element technology allows a time-gated
measurement of the reflection matrix [52–55].

All the authors are grateful for the funding provided by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement 819261, REMINISCENCE project) and by LABEX WIFI (Laboratory of Excellence within the French Program Investments for the Future, ANR-10-LABX24 and ANR-10-IDEX-0001-02 PSL*). C.B. and W.L. acknowledge financial support from Safran and SuperSonic Imagine, respectively.
some limitations. First, the SSc filter is only efficient when the MSc background is not too large. The tomography of $\ell_s$ is therefore limited to a penetration depth of a few $\ell_s$. Another issue can arise when wave propagates from a strongly scattering region to a weakly scattering area. MSc induced by the first layer predominates and hides the SSc component of the second one [32]. A more elaborate inversion procedure of $\rho_s$ is thus needed and will be the object of a future study. Another relevant question concerns the applicability of this method to a 2D probe, i.e a 3D imaging geometry. Interestingly, an exponential decay is still observed at short times-of-flight in 3D: $\rho_s(t) \sim \exp(-0.55ct/\ell_s)$ [32]. Hence the proposed method can be extended to 3D configurations. Another perspective for this work is its extension to anisotropic scattering. This achievement is particularly crucial for optical microscopy since light scattering is sharply peaked in the forward direction in biological tissues.

In summary, we have introduced a novel methodology for local characterization of scattering media, employing a non-invasive reflection setup. Our approach successfully discriminates between SSc and MSc events, providing an independent measurement of $\ell_s$ and $\ell_a$. This obviates the need for traditional transmission experiments, offering a solution to a longstanding challenge in disordered media research. Beyond its inherent advantages, our method stands as a potent tool for quantitative imaging, applicable to in-vivo or in-situ scenarios. Notably, the robustness of the SSc rate in the face of sample reflectivity fluctuations sets our approach apart from techniques reliant on the time-dependence of backscattered echoes for attenuation measurements. Importantly, while our proof-of-concept experiments were conducted with ultrasound, the universality of matrix imaging renders it applicable to diverse fields of wave physics wherever multi-element technology allows a time-resolved measurement of the reflection matrix [48–50].

All the authors are grateful for the funding provided by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement 819261, REMINISCENCE project) and by LABEX WIFI (Laboratory of Excellence within the French Program Investments for the Future, ANR-10-LABX24 and ANR-10-IDEX-0001-02 PSL*). C.B. and W.L. acknowledge financial support from Safran and SuperSonic Imagine, respectively.
# Supplementary Material

This Supplementary Material provides further information on: (i) the focused reflection matrix; (ii) the single scattering filter; (iii) the relation between the canonical and focused R-matrices; (iv) the theoretical expression of the single scattering intensity; (v) the theoretical expression of the overall intensity in the short-time range; (vi) the expression of the diffuse intensity in the long-time limit; (vii) the numerical validation of the whole method.

## FOCUSED REFLECTION MATRIX

Matrix imaging basically consists in applying a focused beamforming process to the recorded reflection matrix $R_{uu}(t)$ at input and output [42]. To do so, the first step is to perform a temporal Fourier transform of $R_{uu}(t)$:

$$ \overline{\mathbf{R}}_{uu}(\omega) = \int dt \mathbf{R}_{uu}(t) e^{-j\omega t}. \qquad (S1) $$

with $\omega$, the angular frequency. The second step is a projection of the reflection matrix in the focused basis. Mathematically, this operation can be performed in the frequency domain by means of the following matrix product:

$$ \overline{\mathbf{R}}_{xx}(z, \omega) = \mathbf{T}_{ux}^{\dagger}(z, \omega) \times \mathbf{R}_{uu}(\omega) \times \mathbf{T}_{ux}^{*}(z, \omega). \qquad (S2) $$

where the symbols $\dagger$ and $*$ stand for transpose conjugation and conjugation operations. $T_{ux}(z)$ is the transmission matrix that describes wave propagation from the transducer plane (u) to the focal plane (x) at each depth $z$. Based on diffraction theory and Rayleigh Sommerfeld integral, its elements $T(u, x, z, \omega)$ correspond to the $z$-derivative of the Green's functions $G(u, x, z, \omega)$ between the transducer and focal planes. $G_{ur}$ depends on the spatial distribution of the wave speed $c(\mathbf{r})$ inside the system. For an homogeneous speed of sound $c_0$ and in a 2D geometry, these Green's functions are given by:

$$ G(u, x, z, \omega) = -\frac{j}{4} \mathcal{H}_0 \left( k \sqrt{(u-x)^2 + z^2} \right) \qquad (S3) $$

where $\mathcal{H}_0$ is the Hankel function of the first kind and $k = \omega/c_0$, the wave number.

The coefficients $\overline{R}(x_{\text{out}}, x_{\text{in}}, z, \omega)$ of each focused reflection matrix $\overline{\mathbf{R}}_{xx}(z, \omega)$ are the response of the medium at angular frequency $\omega$ between virtual sources and detectors located
# Supplementary Material

This Supplementary Material provides further information on: (i) the focused reflection matrix; (ii) the single scattering filter; (iii) the relation between the canonical and focused R-matrices; (iv) the theoretical expression of the single scattering intensity; (v) the theoretical expression of the overall intensity in the short-time range; (vi) the expression of the diffuse intensity in the long-time limit; (vii) the numerical validation of the whole method.

## FOCUSED REFLECTION MATRIX

Matrix imaging basically consists in applying a focused beamforming process to the recorded reflection matrix $R_{uu}(t)$ at input and output [35]. To do so, the first step is to perform a temporal Fourier transform of $R_{uu}(t)$:

$$ \overline{\mathbf{R}}_{uu}(\omega) = \int dt \mathbf{R}_{uu}(t) e^{-j\omega t}. \qquad (S1) $$

with $\omega$, the angular frequency. The second step is a projection of the reflection matrix in the focused basis. Mathematically, this operation can be performed in the frequency domain by means of the following matrix product:

$$ \overline{\mathbf{R}}_{xx}(z, \omega) = \mathbf{T}_{ux}^{\dagger}(z, \omega) \times \mathbf{R}_{uu}(\omega) \times \mathbf{T}_{ux}^{*}(z, \omega). \qquad (S2) $$

where the symbols $\dagger$ and $*$ stand for transpose conjugation and conjugation operations. $\mathbf{T}_{ux}(z)$ is the transmission matrix that describes wave propagation from the transducer plane (u) to the focal plane (x) at each depth $z$. Based on diffraction theory and Rayleigh Sommerfeld integral, its elements $T(u, x, z, \omega)$ correspond to the $z$-derivative of the Green's functions $G(u, x, z, \omega)$ between the transducer and focal planes. $\mathbf{G}_{ur}$ depends on the spatial distribution of the wave speed $c(\mathbf{r})$ inside the system. For an homogeneous speed of sound $c_0$ and in a 2D geometry, these Green's functions are given by:

$$ G(u, x, z, \omega) = -\frac{j}{4}\mathcal{H}_0\left(k\sqrt{(u-x)^2 + z^2}\right) \qquad (S3) $$

where $\mathcal{H}_0$ is the Hankel function of the first kind and $k = \omega/c_0$, the wave number.

The coefficients $\overline{R}(x_{\text{out}}, x_{\text{in}}, z, \omega)$ of each focused reflection matrix $\overline{\mathbf{R}}_{xx}(z, \omega)$ are the response of the medium at angular frequency $\omega$ between virtual sources and detectors located
at $\mathbf{r}_{in} = (x_{in}, z)$ and $\mathbf{r}_{out} = (x_{out}, z)$. The last step consists in summing the focused reflection matrix over the frequency bandwidth $[\omega_{-}; \omega_{+}]$:

$$ \mathbf{R}_{xx}(z) = \int_{\omega_{-}}^{\omega_{+}} d\omega \bar{\mathbf{R}}_{xx}(z, \omega) e^{j\omega t}. \qquad (S4) $$

This operation is equivalent to an inverse Fourier transform at lapse time $\tau = 0$ in the focused basis; it amounts to a time-gating of back-scattered echoes at the ballistic time $t = 2z/c$ in the transducer basis. It is required to recover the axial resolution of the ultrasound image such that selected singly-scattered waves come from a single plane at depth $z$.

## SINGLE SCATTERING FILTER

To filter multiple scattering from the focused reflection matrix, an accurate single scattering space should be built [45]. To that aim, we first compute a set of monochromatic reflection matrices $\mathbf{P}_{uu}(x_s, z, \omega) = [P(u_{in}, u_{out}, x_s, z, \omega)]$ defined in the transducer basis and associated with a single point-like scatterer located at $(x_s, z)$, such that

$$ \mathbf{P}_{uu}(x_s, z, \omega) = \mathbf{T}_{ux}(z, \omega) \times \mathbf{\Gamma}_{xx}(x_s) \times \mathbf{T}_{ux}^{T}(z, \omega). \qquad (S5) $$

with $\mathbf{\Gamma}(x_s)$, a diagonal matrix whose coefficients $\gamma(x)$ describes the reflectivity of each scatterer, such that $\gamma(x) = \gamma_0\delta(x - x_s)$, with $\delta$ the Dirac distribution.

The next step is the projection of each matrix $\mathbf{P}_{uu}(x_s, z, \omega)$ in the focused basis, as previously done for $\mathbf{R}_{uu}(z, \omega)$ [Eq. (S2)],

$$ \mathbf{P}_{xx}(x_s, z, \omega) = \mathbf{T}_{ux}^{\dagger}(z, \omega) \times \mathbf{P}_{uu}(x_s, z, \omega) \times \mathbf{T}_{ux}^{*}(z, \omega). \qquad (S6) $$

In practice, we generate a synthetic matrix $\mathbf{P}_{xx}(x_s, z, \omega)$ for a virtual scatterer every half resolution cell. The dimension of this resolution cell is dictated by the characteristic size $\delta x$ of the focal spot. For an homogeneous speed-of-sound, $\delta x = \lambda/(2 \sin[\arctan(z/(2Np))])$.

To build a single scattering space, this set of matrices can be orthogonalized by means a Gram-Schmidt process [45]. However, this approach is not optimal and time consuming. Here we propose a more efficient process based on a singular value decomposition. To that aim, the set of reference matrices is first concatenated into a bi-dimensional matrix $\hat{\mathbf{P}}(z, \omega)$, at each depth $z$ and frequency $\omega$ such that:

$$ \hat{\mathbf{P}}(\{x_{in}, x_{out}\}, x_s, z, \omega) = P(x_{in}, x_{out}, x_s, z, \omega) \qquad (S7) $$
at $\mathbf{r}_{\text{in}} = (x_{\text{in}}, z)$ and $\mathbf{r}_{\text{out}} = (x_{\text{out}}, z)$. The last step consists in summing the focused reflection matrix over the frequency bandwidth $[\omega_{-}; \omega_{+}]$:

$$ \mathbf{R}_{xx}(z) = \int_{\omega_{-}}^{\omega_{+}} d\omega \overline{\mathbf{R}}_{xx}(z, \omega) e^{j\omega t}. \qquad (S4) $$

This operation is equivalent to an inverse Fourier transform at lapse time $\tau = 0$ in the focused basis; it amounts to a time-gating of back-scattered echoes at the ballistic time $t = 2z/c$ in the transducer basis. It is required to recover the axial resolution of the ultrasound image such that selected singly-scattered waves come from a single plane at depth $z$.

## SINGLE SCATTERING FILTER

To filter multiple scattering from the focused reflection matrix, an accurate single scattering space should be built [42]. To that aim, we first compute a set of monochromatic reflection matrices $\mathbf{P}_{uu}(x_s, z, \omega) = [P(u_{\text{in}}, u_{\text{out}}, x_s, z, \omega)]$ defined in the transducer basis and associated with a single point-like scatterer located at $(x_s, z)$, such that

$$ \mathbf{P}_{uu}(x_s, z, \omega) = \mathbf{T}_{ux}(z, \omega) \times \mathbf{\Gamma}_{xx}(x_s) \times \mathbf{T}_{ux}^{T}(z, \omega). \qquad (S5) $$

with $\mathbf{\Gamma}(x_s)$, a diagonal matrix whose coefficients $\gamma(x)$ describes the reflectivity of each scatterer, such that $\gamma(x) = \gamma_0\delta(x - x_s)$, with $\delta$ the Dirac distribution.

The next step is the projection of each matrix $\mathbf{P}_{uu}(x_s, z, \omega)$ in the focused basis, as previously done for $\mathbf{R}_{uu}(z, \omega)$ [Eq. (S2)],

$$ \mathbf{P}_{xx}(x_s, z, \omega) = \mathbf{T}_{ux}^{\dagger}(z, \omega) \times \mathbf{P}_{uu}(x_s, z, \omega) \times \mathbf{T}_{ux}^{*}(z, \omega). \qquad (S6) $$

In practice, we generate a synthetic matrix $\mathbf{P}_{xx}(x_s, z, \omega)$ for a virtual scatterer every half resolution cell. The dimension of this resolution cell is dictated by the characteristic size $\delta x$ of the focal spot. For an homogeneous speed-of-sound, $\delta x = \lambda/(2 \sin[\arctan(z/(2Np))])$.

To build a single scattering space, this set of matrices can be orthogonalized by means a Gram-Schmidt process [42]. However, this approach is not optimal and time consuming. Here we propose a more efficient process based on a singular value decomposition. To that aim, the set of reference matrices is first concatenated into a bi-dimensional matrix $\hat{\mathbf{P}}(z, \omega)$, at each depth $z$ and frequency $\omega$ such that:

$$ \hat{\mathbf{P}}(\{x_{\text{in}}, x_{\text{out}}\}, x_s, z, \omega) = \mathbf{P}(x_{\text{in}}, x_{\text{out}}, x_s, z, \omega) \qquad (S7) $$
Then, the single scattering basis is determined at each depth and frequency by performing a singular value decomposition of $\hat{\mathbf{P}}(z, \omega)$:

$$ \hat{\mathbf{P}}(z, \omega) = \sum_{k=1}^{N_k(z, \omega)} \sigma_k(z, \omega) \hat{\mathbf{F}}_k(z, \omega) \mathbf{X}_k(z, \omega) \quad (S8) $$

where $\sigma_k$ are the singular values ranged in decreasing order. $\hat{\mathbf{F}}_k(z, \omega) = [\hat{F}_k(\{x_{in}, x_{out}\}, z, \omega)]$ and $\mathbf{X}_k(z, \omega) = [X_k(x_s, z, \omega)]$ are the set of singular vectors defined in the basis of focused reflection matrices ($\{x_{in}, x_{out}\}$) and the scatterer position basis ($x_s$). The rank $N_k(z, \omega)$ of the single scattering space can be tuned by applying an arbitrary threshold to the spectrum of singular values or deterministically fixed as the number of resolution cells in the field-of-view: $N_k(z, \omega) = \Delta x/\delta x(z, \omega)$, with $\Delta x$, the lateral extent of the field-of-view.

The set of singular vectors $\hat{\mathbf{F}}_k(z, \omega)$ provide an orthonormal basis of focused reflection matrices $\mathbf{F}_k(z, \omega) = [F_k(x_{in}, x_{out}, z, \omega)]$ [Fig.2(b) in the accompanying paper] such that

$$ F_k(x_{in}, x_{out}, z, \omega) = \hat{F}_k(\{x_{in}, x_{out}\}, z, \omega). \quad (S9) $$

This single scattering basis is used to project each focused reflection matrix $\bar{\mathbf{R}}_{xx}(z, \omega)$ onto a characteristic single scattering space:

$$ \bar{\mathbf{S}}_{xx}(z, \omega) = \sum_{k=1}^{N_k(z, \omega)} \mathrm{Tr} \left\{ \mathbf{F}_k^{\dagger}(z, \omega) \times \bar{\mathbf{R}}_{xx}(z, \omega) \right\} \mathbf{F}_k(z, \omega) \quad (S10) $$

The set of resulting monochromatic matrices $\bar{\mathbf{S}}_{xx}(z, \omega)$ are then recombined to yield the broadband single scattering matrix:

$$ \mathbf{S}_{xx}(z) = \int_{\omega_{-}}^{\omega_{+}} d\omega \bar{\mathbf{S}}_{xx}(z, \omega) \quad (S11) $$

Figure 2(d) of the accompanying manuscript show one example of single scattering matrix in the phantom experiment. Not surprisingly, its diagonal coefficients are significantly larger than off-diagonal coefficients: this is characteristic of single scattering. The ratio between the norms of $\mathbf{S}_{xx}(z)$ and $\mathbf{R}_{xx}(z)$ can serve as an estimator for the single scattering ratio, a refined one compared to previous approaches that estimated the single and multiple scattering rates from the off-diagonal and diagonal intensities [42]. Moreover, unlike Ref. [44], the present filter does not make any assumption on the intensity distribution of multiple scattering. Despite these improvements, this single scattering filter still exhibits a bias for $t > \tau_s$ since, even in absence of single scattering, a residual multiple scattering contribution will still
Then, the single scattering basis is determined at each depth and frequency by performing a singular value decomposition of $\hat{\mathbf{P}}(z, \omega)$:

$$ \hat{\mathbf{P}}(z, \omega) = \sum_{k=1}^{N_k(z, \omega)} \sigma_k(z, \omega) \hat{\mathbf{F}}_k(z, \omega) \mathbf{X}_k(z, \omega) \quad (S8) $$

where $\sigma_k$ are the singular values ranged in decreasing order. $\hat{\mathbf{F}}_k(z, \omega) = [\hat{F}_k(\{x_{in}, x_{out}\}, z, \omega)]$ and $\mathbf{X}_k(z, \omega) = [X_k(x_s, z, \omega)]$ are the set of singular vectors defined in the basis of focused reflection matrices $(\{x_{in}, x_{out}\})$ and the scatterer position basis $(x_s)$. The rank $N_k(z, \omega)$ of the single scattering space is deterministically fixed as the number of resolution cells in the field-of-view: $N_k(z, \omega) = \Delta x/\delta x(z, \omega)$, with $\Delta x$, the lateral extent of the field-of-view.

The set of singular vectors $\hat{\mathbf{F}}_k(z, \omega)$ provide an orthonormal basis of focused reflection matrices $\mathbf{F}_k(z, \omega) = [F_k(x_{in}, x_{out}, z, \omega)]$ [Fig.2(b) in the accompanying paper] such that

$$ F_k(x_{in}, x_{out}, z, \omega) = \hat{F}_k(\{x_{in}, x_{out}\}, z, \omega). \quad (S9) $$

This single scattering basis is used to project each focused reflection matrix $\overline{\mathbf{R}}_{xx}(z, \omega)$ onto a characteristic single scattering space:

$$ \overline{\mathbf{S}}_{xx}(z, \omega) = \sum_{k=1}^{N_k(z, \omega)} \mathrm{Tr} \left\{ \mathbf{F}_k^{\dagger}(z, \omega) \times \overline{\mathbf{R}}_{xx}(z, \omega) \right\} \mathbf{F}_k(z, \omega) \quad (S10) $$

The set of resulting monochromatic matrices $\overline{\mathbf{S}}_{xx}(z, \omega)$ are then recombined to yield the broadband single scattering matrix:

$$ \mathbf{S}_{xx}(z) = \int_{\omega_{-}}^{\omega_{+}} d\omega \overline{\mathbf{S}}_{xx}(z, \omega) \quad (S11) $$

Figure 2(d) of the accompanying manuscript show one example of single scattering matrix in the phantom experiment. Not surprisingly, its diagonal coefficients are significantly larger than off-diagonal coefficients: this is characteristic of single scattering. The ratio between the norms of $\mathbf{S}_{xx}(z)$ and $\mathbf{R}_{xx}(z)$ can serve as an estimator for the single scattering ratio, a refined one compared to previous approaches that estimated the single and multiple scattering rates from the off-diagonal and diagonal intensities [35]. Moreover, unlike Ref. [41], the present filter does not make any assumption on the intensity distribution of multiple scattering.

Despite these improvements, the single scattering estimator $\rho_s$ can still exhibit a bias $\rho_b$. Indeed, even in absence of single scattering, a residual multiple scattering contribution will still emerge along the single scattering subspace. Assuming the multiple scattering
emerge along the single scattering subspace. However, in the scattering regime considered in this paper ($t < \tau_s$), single scattering is predominant and this bias does not affect our estimator of the single scattering rate.

## RELATION BETWEEN THE FOCUSED R-MATRIX AND THE TIME-GATED CANONICAL R-MATRIX

In the accompanying paper, the back-scattered intensity is defined as the norm of the focused reflection matrix:

$$I(t = 2z/c) = \mathrm{Tr} \{ \mathbf{R}_{xx}(z) \times \mathbf{R}_{xx}^{\dagger}(z) \} / N \quad (S12)$$

To show the equivalence of this quantity with the back-scattered intensity that would be computed in the transducer basis, the first step is to decompose $\mathbf{R}_{xx}$ in the frequency domain [Eq. (S4)] and use the expression of $\bar{\mathbf{R}}_{xx}(\omega)$ given in Eq. (S2) to express the broadband matrix $\mathbf{R}_{xx}$:

$$\mathbf{R}_{xx}(z) = \int d\omega \mathbf{T}_{ux}^{\dagger}(z, \omega) \times \mathbf{R}_{uu}(\omega) \times \mathbf{T}_{ux}^{*}(z, \omega). \quad (S13)$$

To simplify the latter expression, the paraxial approximation can be made and the transmission matrix can be decomposed as follows:

$$\mathbf{T}_{ux}(z, \omega) = \mathbf{F}_{ux}(z, \omega) \exp \left[ j \frac{\omega}{c} z \right] \quad (S14)$$

with $\mathbf{F}_{ux}(z, \omega)$, the Fresnel operator whose coefficients $F(u, x, z, \omega)$ are given by:

$$F(\mathbf{u}, x, z, \omega) \approx e^{3i\pi/4} \sqrt{\frac{2\omega}{\pi z c}} \exp \left[ j \frac{\omega}{2zc} (u-x)^2 \right]. \quad (S15)$$

Using this decomposition of $\mathbf{T}_{ux}(z, \omega)$, Equation (S13) can be rewritten as follows:

$$\mathbf{R}_{xx}(z) \approx \int d\omega \mathbf{F}_{ux}^{\dagger}(z, \omega) \times \mathbf{R}_{uu}(\omega) \exp \left[ -2j \frac{\omega}{c} z \right] \times \mathbf{F}_{ux}^{*}(z, \omega). \quad (S16)$$

To go beyond, the Fresnel operator can be considered as constant over the frequency band- with and replaced by its value at the central frequency $\omega_c$, such that

$$\mathbf{R}_{xx}(z) \approx \mathbf{F}_{ux}^{\dagger}(z, \omega_c) \times \mathbf{R}_{uu}(t = 2z/c) \times \mathbf{F}_{ux}^{*}(z, \omega_c). \quad (S17)$$

The broadband focused reflection matrix $\mathbf{R}_{xx}(z)$ and the time-gated reflection matrix $\mathbf{R}_{uu}(t = 2z/c)$ are thus related through a simple change of basis. The unitarity of the
contribution corresponds to a fully random reflection matrix, one can show the bias $\rho_b$ is of the order of $1/N_k$. In reality, the multiple scattering intensity profile is of finite spatial extent $W$ [41]. Moreover, it exhibits an enhancement by a factor two along the diagonal of the focused reflection matrix due to the coherent back-scattering phenomenon [35]. Taking into account these two features, the bias is then given by

$$ \rho_b \sim 2/(1 + N_W), \qquad (S12) $$

where $N_W = W/\delta x$ is the ratio between the spatial extent $W$ of the diffuse halo and the resolution cell $\delta x$. Our estimator $\hat{\rho}_s$ of the single scattering ratio is therefore valid only if $\rho_b \ll 1$. In the phantom and liver experiments described in the accompanying paper, the diffuse halo spreads over the whole matrix ($N_W \sim N$) and the bias $\rho_b$ is negligible.

However, in other configurations, the finite size of the diffuse halo and the coherent back-scattering phenomenon can increase the overlap between the single and multiple scattering sub-spaces. The bias $\rho_b$ might not be negligible but an unbiased estimator $\hat{\rho}'_s$ can be defined as follows:

$$ \hat{\rho}'_s = \frac{\hat{\rho}_s - \rho_b}{1 - \rho_b}. \qquad (S13) $$

A critical step to build this new estimator is to evaluate $\rho_b$. A first option is to compute it by generating random matrices characteristic of the experimental configuration and multiple scattering properties. A second option is to estimate it from the values of $\hat{\rho}_s$ at large echo times for which single scattering is a priori negligible. The latter method is actually used in Supplementary Section to build an unbiased estimator of the single scattering ratio for numerical experiments in a far-field configuration.

## RELATION BETWEEN THE FOCUSED R-MATRIX AND THE TIME-GATED CANONICAL R-MATRIX

In the accompanying paper, an equality is assumed between the back-scattered intensity in the transducer basis, $I(t) = \text{Tr}\{\mathbf{R}_{uu}(t) \times \mathbf{R}_{xx}^{\dagger}(t)\}/N$, and the norm of the focused reflection matrix, $\text{Tr}\{\mathbf{R}_{xx}(z) \times \mathbf{R}_{xx}^{\dagger}(z)\}/N$ at the ballistic depth $z = ct/2$ [Eq. (1)].

To show under which approximation this equality is valid, the first step is to decompose $\mathbf{R}_{xx}$ in the frequency domain [Eq. (S4)] and use the expression of $\bar{\mathbf{R}}_{xx}(\omega)$ given in Eq. (S2)
Fresnel operator $\mathbf{F}_{ux}(z, \omega_c)$, implies an equality between their norms:

$$I(t) = \mathrm{Tr} \{ \mathbf{R}_{xx}(z) \times \mathbf{R}_{xx}^{\dagger}(z) \} / N = \mathrm{Tr} \{ \mathbf{R}_{uu}(t = 2z/c) \times \overline{\mathbf{R}}_{uu}^{\dagger}(t = 2z/c) \} / N \quad (S18)$$

Upon paraxial approximation, it is therefore equivalent to compute the back-scattered intensity in the transducer or focused basis. The same relation can be demonstrated for the single scattering contribution:

$$I_S(t) = \mathrm{Tr} \{ \mathbf{S}_{xx}(z) \times \mathbf{S}_{xx}^{\dagger}(z) \} / N = \mathrm{Tr} \{ \mathbf{S}_{uu}(t = 2z/c) \times \overline{\mathbf{S}}_{uu}^{\dagger}(t = 2z/c) \} / N \quad (S19)$$

The analytical expression of $I(t)$ and $I_S(t)$ can thus be derived by considering the response of the medium at its surface. Note that the same result holds if the array of transducers is further away from the sample surface. In that case, one simply has to set the depth origin ($z = 0$) at the sample surface, and accordingly the origin of time ($t = 0$) at the arrival time of the first echo.

## SINGLE SCATTERING INTENSITY

We first derive an expression for the single-scattering contribution based on radiative transfer theory under isotropic scattering assumption. In that case, there is no difference between the scattering and the transport mean-free path. Radiative transfer theory describes the spatial and temporal dependence of the radiance (or specific intensity) $P(\mathbf{r}, t, \mathbf{u})$ in a random medium. Radiance is defined as energy flow, propagating in the direction $\mathbf{u}$, per unit normal area per unit solid angle $d\Omega$ per unit time. It follows a transport (Boltzmann) equation:

$$\frac{1}{c} \frac{\partial}{\partial t} P(\mathbf{r}, t, \mathbf{u}) + \mathbf{u} \cdot \nabla P(\mathbf{r}, t, \mathbf{u}) + \ell_{\mathrm{ext}}^{-1} P(\mathbf{r}, t, \mathbf{u}) = \ell_s^{-1} P(\mathbf{r}, t) + c^{-1} S(\mathbf{r}, t, \mathbf{u}), \quad (S20)$$

with $S$ the source term and $P(\mathbf{r}, t)$ the intensity, defined as the angular average of the radiance

$$P(\mathbf{r}, t) = \frac{1}{2\pi} \int d\Omega P(\mathbf{r}, t, \mathbf{u}) \quad (S21)$$

Classically, the transport equation can be derived from the Bethe-Salpether equation, neglecting all interference (coherent) effects. Here we adapt the theoretical developments of Paasschens [41], which were derived for an infinite random medium, to our experimental configuration. The problem was solved in two dimensions. In an elegant way, Passchens
to express the broadband matrix $\mathbf{R}_{xx}$:

$$ \mathbf{R}_{xx}(z) = \int d\omega \mathbf{T}_{ux}^{\dagger}(z, \omega) \times \mathbf{R}_{uu}(\omega) \times \mathbf{T}_{ux}^{*}(z, \omega). \quad (S14) $$

To simplify the latter expression, the paraxial approximation can be made and the transmission matrix can be decomposed as follows:

$$ \mathbf{T}_{ux}(z, \omega) = \mathbf{F}_{ux}(z, \omega) \exp \left[ j \frac{\omega}{c} z \right] \quad (S15) $$

with $\mathbf{F}_{ux}(z, \omega)$, the Fresnel operator whose coefficients $F(u, x, z, \omega)$ are given by:

$$ F(\mathbf{u}, x, z, \omega) \approx e^{3i\pi/4} \sqrt{\frac{2\omega}{\pi z c}} \exp \left[ j \frac{\omega}{2zc} (u-x)^2 \right]. \quad (S16) $$

Using this decomposition of $\mathbf{T}_{ux}(z, \omega)$, Equation (S14) can be rewritten as follows

$$ \mathbf{R}_{xx}(z) \approx \int d\omega \mathbf{F}_{ux}^{\dagger}(z, \omega) \times \mathbf{R}_{uu}(\omega) \exp \left[ -2j \frac{\omega}{c} z \right] \times \mathbf{F}_{ux}^{*}(z, \omega). \quad (S17) $$

To go beyond, the Fresnel operator can be considered as constant over the frequency band- with and replaced by its value at the central frequency $\omega_c$, such that

$$ \mathbf{R}_{xx}(z) \approx \mathbf{F}_{ux}^{\dagger}(z, \omega_c) \times \mathbf{R}_{uu}(t = 2z/c) \times \mathbf{F}_{ux}^{*}(z, \omega_c). \quad (S18) $$

The latter approximation is valid as long as the phase variation of the $\mathbf{F}$-matrix coefficients is negligible over the frequency bandwidth, such that

$$ \frac{\Delta\omega D^2}{2zc} \ll \pi, \quad (S19) $$

with $D$, the dimension of the transducer array. This latter equation can be thought of a far-field condition:

$$ z \gg \frac{\Delta f D^2}{f \lambda}. \quad (S20) $$

The broadband focused reflection matrix $\mathbf{R}_{xx}(z)$ and the time-gated reflection matrix $\mathbf{R}_{uu}(t = 2z/c)$ are thus related through a simple change of basis. The unitarity of the Fresnel operator $\mathbf{F}_{ux}(z, \omega_c)$ then implies an equality between their norms:

$$ I(t) = \mathrm{Tr} \left\{ \mathbf{R}_{uu}(t) \times \overline{\mathbf{R}}_{uu}^{\dagger}(t) \right\} / N \approx \mathrm{Tr} \left\{ \mathbf{R}_{xx}(z = ct/2) \times \mathbf{R}_{xx}^{\dagger}(z = ct/2) \right\} / N \quad (S21) $$

Upon paraxial approximation, it is therefore equivalent to compute the back-scattered intensity in the transducer or focused basis. The same relation can be demonstrated for the single scattering contribution:

$$ I_S(t) = \mathrm{Tr} \left\{ \mathbf{S}_{uu}(t) \times \overline{\mathbf{S}}_{uu}^{\dagger}(t) \right\} / N \approx \mathrm{Tr} \left\{ \mathbf{S}_{xx}(z = ct/2) \times \mathbf{S}_{xx}^{\dagger}(z = ct/2) \right\} / N \quad (S22) $$
derived an expression for each scattering order of the intensity Green's function $P(r, t)$. In particular, the single scattering component $P_S$ is derived as follows:

$$P_S(r, t) = \frac{\exp(-ct/\ell_{\text{ext}})}{2\pi\ell_s ct} \left(1 - \frac{r^2}{c^2 t^2}\right)^{-1/2} \Theta(ct - r) \quad (S22)$$

with $\Theta$, the Heaviside function: $\Theta(x) = 1$ for $x > 0$ and zero elsewhere.

In our case, the medium is semi-infinite. The question is now how to modify the result of Paasschens to account for the medium interface. If we assume that the source term is isotropic, and that the average reflection coefficient for the interface at $z = 0$ is 0, we just need to divide by two the result of Paasschens for the single scattering component. The single scattering intensity $I_S$ can then be obtained by integrating the mean radiance over the sample surface:

$$I_S(t) = \frac{1}{2} \int_{-\infty}^{+\infty} dx P_S(x, z = 0, t) \quad (S23)$$

Injecting Eq. (2) into the last equation leads to the following expression for $I_S(t)$:

$$I_S(t) = \frac{\exp(-ct/\ell_{\text{ext}})}{4\ell_s} \quad (S24)$$

This analytical expression perfectly agrees with the result of the Monte Carlo simulation displayed in Fig. 1(b). This simulation consists in a 2D random walk with an exponential step size distribution of characteristic length $l_s$. $10^6$ particles are thrown in the medium from the origin ($z = 0$) with a uniform angular distribution. The resulting time-dependent backscattered intensity is estimated by building the histogram of times-of-flight at which each reflected particle crosses the interface at $z = 0$. This Monte Carlo simulation allows the independent investigation of each scattering order of the backscattered intensity since we know exactly the number of scattering events undergone by each particle before reaching the medium interface. The single scattering intensity displayed in Fig. 1 has thus been obtained by considering the time-of-flight distribution for particles that have been scattered once and only once before reaching the medium interface. Figure S1(a) also shows the time-of-flight distribution for the single scattering component but, this time, for a 3D random walk. Interestingly, the expression of $I_S$ given in Eq. (S24) also seems to fit pretty well the Monte Carlo result in a 3D geometry [see Fig. S1(a)].
The analytical expression of $I(t)$ and $I_S(t)$ can thus be derived by considering the response of the medium at its surface. Note that the same result holds if the array of transducers is further away from the sample surface. In that case, one simply has to set the depth origin ($z = 0$) at the sample surface, and accordingly the origin of time ($t = 0$) at the arrival time of the first echo.

The paraxial approximation of the wave-field is therefore required for a strict equality between the back-scattered intensity and the norm of the focused reflection matrix. Nevertheless, the condition Eq. S20 is restrictive and can be relaxed by considering that the characteristic scattering time $\tau_s$ is much larger than the wave period ($k0\ell_s \gg 1$ in a weak scattering regime). To do so, one can consider the relation that exists between the matrices $R_{uu}$ and $R_{xx}$ in the time domain. The coefficients of $R_{xx}(z)$ can also be seen as the result of a delay-and-sum beamforming process: [41]:

$$R(x_{\text{out}}, x_{\text{in}}, z) = \sum_{u_{\text{in}}} \sum_{u_{\text{out}}} A(u_{\text{out}}, x_{\text{out}}, z) A(u_{\text{in}}, x_{\text{in}}, z) R(u_{\text{out}}, u_{\text{in}}, 2z/c_0 + \tau(u_{\text{in}}, x_{\text{in}}, z) + \tau(u_{\text{out}}, x_{\text{out}}, z)), \quad (S23)$$

where $\tau$ is the time delay law such that

$$\tau(u, x, z) = \{\sqrt{(x-u)^2 + z^2} - z\}/c_0. \quad (S24)$$

$A$ is an apodization factor that can limit the extent of the synthetic aperture and that is normalized such that

$$\sum_{u_{\text{in}}} \sum_{u_{\text{out}}} A(u_{\text{out}}, x_{\text{out}}, z) A(u_{\text{in}}, x_{\text{in}}, z) = 1 \quad (S25)$$

For sake of simplicity, let us consider $A$ as a constant factor ($A = 1/N$) and the $R$-matrix coefficients as fully uncorrelated (random medium). The mean intensity of Eq. S23 can be written as follows:

$$\langle |R(x_{\text{out}}, x_{\text{in}}, z)|^2 \rangle = \langle |R(u_{\text{out}}, u_{\text{in}}, 2z/c_0 + \tau(u_{\text{in}}, x_{\text{in}}, z) + \tau(u_{\text{out}}, x_{\text{out}}, z))|^2 \rangle. \quad (S26)$$

The latter equation leads to Eq. S21 if the sum of time delays, $\tau(u_{\text{in}}, x_{\text{in}}, z) + \tau(u_{\text{out}}, x_{\text{out}}, z)$, is negligible compared to the characteristic time over which the mean intensity fluctuates, the mean elastic time $\tau_s = \ell_s/c$. Applying paraxial approximation to Eq. S24 leads to the following condition for the validity of Eq. S21:

$$z \gg \frac{D^2}{\ell_s} \quad (S27)$$
<figure><img src="image_6.png" /><figcaption>FIG. S1. Monte Carlo simulation of energy transport in a 3D semi-infinite random medium [$\ell_s = 0.1$ m, $\ell_a = 100$ m]. (a) (b) Reflected intensity vs. $t/\tau_s$: The single scattering intensity (green dots) is compared with its theoretical expression of Eq. (2); the overall intensity (blue dots) is compared with the heuristic expression of $I(t) = I_S(t) \exp(-0.55ct/\ell_s)$ (red line) and the 3D diffusive result [47] (black line). (c) Single scattering rate $\rho_s(t)$: The numerical result (black dots) is compared with the heuristic exponential decay $\exp(-0.55ct/\ell_s)$ (red line).</figcaption></figure>

## SCATTERED INTENSITY IN THE SHORT-TIME LIMIT

Paasschens also derived the following analytical expression for the overall radiance $P(\mathbf{r}, t)$ in a 2D geometry:

$$P(\mathbf{r}, t) = \frac{\exp(-ct/\ell_a)}{2\pi\ell_s ct} \left(1 - \frac{r^2}{c^2t^2}\right)^{-1/2} \exp\left[\left(\sqrt{c^2t^2 - r^2} - ct\right)/\ell_s\right] \Theta(ct - r). \quad (S25)$$

The square root in the exponential can be developed as follows,

$$\exp\left[\left(\sqrt{c^2t^2 - r^2} - ct\right)/\ell_s\right] \sim \exp\left[-r^2/(2\ell_s ct)\right] \quad (S26)$$

provided that

$$\exp\left[-r^4/(8\ell_s ct^3)\right] \sim 1. \quad (S27)$$

Given the fact that $r < ct$, the latter condition is valid for $ct < \ell_s$. Equation (S25) can thus be simplified into the following expression for $t < \tau_s$:

$$P(\mathbf{r}, t) \sim \frac{\exp(-ct/\ell_a)}{2\pi\ell_s ct} \left(1 - \frac{r^2}{c^2t^2}\right)^{-1/2} \exp\left[-r^2/(2\ell_s ct)\right] \Theta(ct - r). \quad (S28)$$

Taking into account the medium interface in presence of multiple scattering is more tricky than for single scattering. A first option is to reproduce what we did for single scattering, i.e
The latter condition is much less restrictive than Eq. S20. Moreover, at shallow depth, it can be forced by applying apodization factors that can reduce the numerical aperture, $D/(2z)$, of the imaging system.

## SINGLE SCATTERING INTENSITY

We first derive an expression for the single-scattering contribution based on radiative transfer theory under isotropic scattering assumption. In that case, there is no difference between the scattering and the transport mean-free path. Radiative transfer theory describes the spatial and temporal dependence of the radiance (or specific intensity) $P(\mathbf{r}, t, \mathbf{u})$ in a random medium. Radiance is defined as energy flow, propagating in the direction $\mathbf{u}$, per unit normal area per unit solid angle $d\Omega$ per unit time. It follows a transport (Boltzmann) equation:

$$ \frac{1}{c} \frac{\partial}{\partial t} P(\mathbf{r}, t, \mathbf{u}) + \mathbf{u} \cdot \nabla P(\mathbf{r}, t, \mathbf{u}) + \ell_{\text{ext}}^{-1} P(\mathbf{r}, t, \mathbf{u}) = \ell_s^{-1} P(\mathbf{r}, t) + c^{-1} S(\mathbf{r}, t, \mathbf{u}), \quad (S28) $$

with $S$ the source term and $P(\mathbf{r}, t)$ the intensity, defined as the angular average of the radiance

$$ P(\mathbf{r}, t) = \frac{1}{2\pi} \int d\Omega P(\mathbf{r}, t, \mathbf{u}) \quad (S29) $$

Classically, the transport equation can be derived from the Bethe-Salpether equation, neglecting all interference (coherent) effects. Here we adapt the theoretical developments of Paasschens [40], which were derived for an infinite random medium, to our experimental configuration. The problem was solved in two dimensions. In an elegant way, Paasschens derived an expression for each scattering order of the intensity Green's function $P(\mathbf{r}, t)$. In particular, the single scattering component $P_S$ is derived as follows:

$$ P_S(r, t) = \frac{\exp(-ct/\ell_{\text{ext}})}{2\pi\ell_s ct} \left(1 - \frac{r^2}{c^2 t^2}\right)^{-1/2} \Theta(ct - r) \quad (S30) $$

with $\Theta$, the Heaviside function: $\Theta(x) = 1$ for $x > 0$ and zero elsewhere.

In our case, the medium is semi-infinite. The question is now how to modify the result of Paasschens to account for the medium interface. If we assume that the source term is isotropic, and that the average reflection coefficient for the interface at $z = 0$ is 0, we just need to divide by two the result of Paasschens for the single scattering component. The
dividing by two the mean radiance [Eq. (S25)] obtained for an infinite medium. The overall intensity is then given by:

$$I(t) = \frac{1}{2} \int_{-\infty}^{+\infty} dx P(x, z = 0, t) \qquad (S29)$$

Injecting Eq. (S28) into the last equation leads to the following expression for $I(t)$:

$$I(t) \sim_{t<\tau_s} \frac{\exp(-ct/\ell_a) \exp[-ct/(4\ell_s)]}{4\ell_s} I_0\left(\frac{ct}{4\ell_s}\right) \qquad (S30)$$

where $I_0$ is the modified Bessel function of the first kind. The latter function can be considered as equal to 1 for $t < \tau_s$. The following asymptotic expression is thus obtained for $I(t)$ at short times-of-flight:

$$I(t) \sim_{t<\tau_s} \exp(-ct/4\ell_s) \exp(-ct/\ell_a)/(4\ell_s). \qquad (S31)$$

This analytical expression is confronted to the result of the Monte Carlo simulation in Fig. 1(b). A perfect agreement is obtained for $t < \tau_s$ but, as expected, this analytical expression cannot grasp the evolution of the backscattered intensity for $t > \tau_s$.

Using Eqs. (S22) and (S31), an expression can be derived for the single scattering rate, $\rho_s(t) = I_s(t)/I(t)$, in the short-time limit:

$$\rho_s(t) \sim_{t<\tau_s} \exp(-3ct/4\ell_s). \qquad (S32)$$

This is the fundamental result on which the accompanying paper is based. This exponential decay of the single scattering rate is confirmed by the result of the Monte Carlo simulation displayed in Fig. 1(c). Again, an excellent agreement is found between our analytical prediction and the solution of the radiative transfer equation for $t < \tau_s$. Note that the decay of $\rho_s(t)$ is not strictly identical in a 3D geometry but an exponential decay of $\exp(-0.55ct/\ell_s)$ fits pretty well the result of our Monte Carlo simulation for $t < \tau_s$ [Fig. S1(b)].

## DIFFUSE INTENSITY IN THE LONG-TIME LIMIT

In a scattering medium ($\ell_s \ll \ell_a$), the specific intensity $P(\mathbf{r}, t, \mathbf{u})$ can be decomposed as the sum of an isotropic and an anisotropic term [56]:

$$P(\mathbf{r}, t, \mathbf{u}) = \frac{\phi(\mathbf{r}, t)}{2\pi} + \frac{1}{\pi} \mathbf{J}(\mathbf{r}, t) \cdot \mathbf{u} \qquad (S33)$$
single scattering intensity $I_S$ can then be obtained by integrating the mean radiance over the sample surface:

$$I_S(t) = \frac{1}{2} \int_{-\infty}^{+\infty} dx P_S(x, z = 0, t) \qquad (S31)$$

Injecting Eq. (2) into the last equation leads to the following expression for $I_S(t)$:

$$I_S(t) = \frac{\exp(-ct/\ell_{\text{ext}})}{4\ell_s} \qquad (S32)$$

This analytical expression perfectly agrees with the result of the Monte Carlo simulation displayed in Fig. 1(b). This simulation consists in a 2D random walk with an exponential step size distribution of characteristic length $\ell_s$. 10⁶ particles are thrown in the medium from the origin ($z = 0$) with a uniform angular distribution in the half space ($z > 0$). The resulting time-dependent backscattered intensity is estimated by building the histogram of time-of-flight at which each reflected particle crosses the interface at $z = 0$. This Monte Carlo simulation allows the independent investigation of each scattering order of the backscattered intensity since we know exactly the number of scattering events undergone by each particle before reaching the medium interface. The single scattering intensity displayed in Fig. 1 has thus been obtained by considering the time-of-flight distribution for particles that have been scattered once and only once before reaching the medium interface. Figure S1(a) also shows the time-of-flight distribution for the single scattering component but, this time, for a 3D random walk. Interestingly, the expression of $I_S$ given in Eq. (S32) also seems to fit pretty well the Monte Carlo result in a 3D geometry [see Fig. S1(a)].

## SCATTERED INTENSITY IN THE SHORT-TIME LIMIT

Paasschens also derived the following analytical expression for the overall radiance $P(r, t)$ in a 2D geometry:

$$P(r, t) = \frac{\exp(-ct/\ell_a)}{2\pi\ell_s ct} \left(1 - \frac{r^2}{c^2 t^2}\right)^{-1/2} \exp\left[\left(\sqrt{c^2 t^2 - r^2} - ct\right)/\ell_s\right] \Theta(ct - r). \qquad (S33)$$

The square root in the exponential can be developed as follows,

$$\exp\left[\left(\sqrt{c^2 t^2 - r^2} - ct\right)/\ell_s\right] \sim \exp\left[-r^2/(2\ell_s ct)\right] \qquad (S34)$$

provided that

$$\exp[-r^4/(8\ell_s c^3 t^3)] \sim 1. \qquad (S35)$$
with $\phi(\mathbf{r}, t)$, the fluence rate (or intensity), and $\mathbf{J}(\mathbf{r}, t)$, the energy flux. In the long time limit ($t \gg \tau_s$), the radiative transfer equation then leads to a much simpler equation for $\phi(\mathbf{r}, t)$:

$$ \frac{1}{c} \frac{\partial}{\partial t} \phi(\mathbf{r}, t) - D \Delta \phi(\mathbf{r}, t) + \frac{1}{\ell_a} \phi(\mathbf{r}, t) = S(\mathbf{r}, t) \quad (S34) $$

with $D$, the diffusion coefficient. In a 2D geometry and for isotropic scattering, $D = c\ell_s/2$.

For an infinite medium, the Green's function $F$ of the diffusion equation [$S(\mathbf{r}, t) = \delta(\mathbf{r})\delta(t)$] can be expressed as follows:

$$ F(\mathbf{r}, t) = \frac{1}{4\pi Dt} \exp(-ct/\ell_a) \exp(-r^2/4Dt). \quad (S35) $$

In a semi-infinite medium, the cancellation of $\phi(\mathbf{r}, t)$ at the boundary ($z = 0$) can be met by adding a negative or image source of energy to the infinite medium problem. The fluence rate per incident wave packet can then be written as the sum of contributions from two sources at $z = z_0$ and $z = -z_0$, with $z_0 = 2\ell_s/3$, the extrapolation length [48]. The corresponding Green's function $F(\mathbf{r}, t)$ is thus given by:

$$ F(\mathbf{r}, t) = \frac{1}{4\pi Dt} \exp(-ct/\ell_a) \left[ \exp\left(-\frac{x^2 + (z - z_0)^2}{4Dt}\right) + \exp\left(-\frac{x^2 + (z + z_0)^2}{4Dt}\right) \right]. \quad (S36) $$

The associated energy flux $\mathbf{J}(\mathbf{r}, t)$ at the medium surface can be deduced from the Fick's law:

$$ \mathbf{J}(\mathbf{r}, t) = -(D/c) \nabla F(\mathbf{r}, t)|_{z=0}. \quad (S37) $$

The diffuse intensity, $P_d(\mathbf{r}, t) = |\mathbf{J}(\mathbf{r}, t)|$, can then be deduced:

$$ P_d(\mathbf{r}, t) = \frac{\exp(-ct/\ell_a) z_0}{2\pi Dt} \exp\left(-\frac{x^2 + z_0^2}{4Dt}\right) \quad (S38) $$

The total diffuse reflectance, $I_d(t) = \int_{-\infty}^{\infty} dx P_d(x, z, t)$, is finally given by:

$$ I_d(t) = \frac{\exp(-ct/\ell_a) z_0}{\sqrt{\pi Dt}} \exp\left(-\frac{z_0^2}{4Dt}\right) \quad (S39) $$

In the long-time limit, $4Dt \gg z_0^2$. $I_d(t)$ thus simplifies into

$$ I_d(t) = \frac{\exp(-ct/\ell_a) z_0}{\sqrt{\pi Dt} ct}. \quad (S40) $$

Interestingly, this expression of $I_d(t)$ is identical to the one previously derived by Patterson et al. in 3D [47]. This expression can be used to fit the result of the Monte Carlo simulation
<figure><img src="image_6.png" /><figcaption>FIG. S1. Monte Carlo simulation of energy transport in a 3D semi-infinite random medium [$\ell_s = 0.1$ m, $\ell_a = 100$ m]. (a) Reflected intensity vs. $t/\tau_s$: The single scattering intensity (green dots) is compared with its theoretical expression of Eq. (2); the overall intensity (blue dots) is compared with the heuristic expression of $I(t) = I_S(t) \exp(-0.55ct/\ell_s)$ (red line) and the 3D diffusive result [44] (black line). (b) Single scattering ratio $\rho_s(t)$: The numerical result (black dots) is compared with the heuristic exponential decay $\exp(-0.55ct/\ell_s)$ (red line).</figcaption></figure>

Given the fact that $r < ct$, the latter condition is valid for $ct < \ell_s$. Equation (S33) can thus be simplified into the following expression for $t < \tau_s$:

$$P(\mathbf{r}, t) \sim_{t<\tau_s} \frac{\exp(-ct/\ell_a)}{2\pi\ell_s ct} \left(1 - \frac{r^2}{c^2t^2}\right)^{-1/2} \exp\left[-r^2/(2\ell_s ct)\right] \Theta(ct - r). \quad (S36)$$

Taking into account the medium interface in presence of multiple scattering is more tricky than for single scattering. A first option is to reproduce what we did for single scattering, i.e dividing by two the mean radiance [Eq. (S33)] obtained for an infinite medium. The overall intensity is then given by:

$$I(t) = \frac{1}{2} \int_{-\infty}^{+\infty} dx P(x, z=0, t) \quad (S37)$$

Injecting Eq. (S36) into the last equation leads to the following expression for $I(t)$:

$$I(t) \sim_{t<\tau_s} \frac{\exp(-ct/\ell_a) \exp[-ct/(4\ell_s)]}{4\ell_s} I_0\left(\frac{ct}{4\ell_s}\right) \quad (S38)$$

where $I_0$ is the modified Bessel function of the first kind. The latter function can be considered as equal to 1 for $t < \tau_s$. The following asymptotic expression is thus obtained for $I(t)$ at short times-of-flight:

$$I(t) \sim_{t<\tau_s} \exp(-ct/4\ell_s) \exp(-ct/\ell_a)/(4\ell_s). \quad (S39)$$
for $t > \tau_s$ both in 2D [Fig. 1(b)] and 3D [Fig. S1(a)]. If absorption is negligible, a power law decay as $t^{-3/2}$ is thus predicted for the scattered intensity in the diffusive limit.

Using Eqs. (2) and (S40), an expression can be derived for the single scattering rate, $\rho_s(t) \sim I_s(t)/I_d(t)$, in the long-time limit:

$$ \rho_s(t) \sim \frac{3}{8} \sqrt{\frac{\pi}{2}} (ct/\ell_s)^{3/2} \exp(-ct/\ell_s). \qquad (S41) $$

This diffusion approximation predicts correctly the evolution of the single scattering rate in the long-time limit (Fig. 1c). As in the short-time limit, the single scattering rate can give access, in principle, to a measurement of $\ell_s$ independent from $\ell_a$. Nevertheless, the scaling with $\ell_s$ is more complicated than a simple exponential decay.

## NUMERICAL VALIDATION

The full method is now validated by means of a simulation of the wave equation. The scattering medium is a random collection of point-like scatterers. The associated reflection matrix is computed by means of a Born series already described in a previous study [45].

We first consider 2D random systems displaying a uniform distribution of scatterers [Fig. S2]. In the first system [Fig. S2(a)], the scatterer cross-section $\sigma$ is of $2.5 \times 10^{-2}$ mm and their concentration is $n = 0.2 \text{ mm}^{-2}$. Under the independent scattering approximation (ISA) [57], the scattering mean free path $\ell_s \sim (n\sigma)^{-1}$, that is to say $\ell_s \sim 194$ mm in the present case. The array of transducers is placed at 140 mm from the sample surface. The simulated reflection matrix is projected onto the single scattering basis following the method described in Sec. . The resulting single scattering rate $\hat{\rho}_s(t)$ is displayed in Fig. S2(b). A quantitative agreement is found between $\hat{\rho}_s(t)$ and its theoretical scaling in $\exp(-3ct/4\ell_s)$ in the short time limit. The same observation holds for a more concentrated scattering medium for which $n = 1 \text{ mm}^{-2}$ and $\ell_s \sim 39$ mm [Fig. S2(c)]. Again, our theoretical prediction [Eq. (5)] perfectly fits the single scattering rate $\hat{\rho}_s(t)$ for $t < \tau_s$ [see Fig. S2(d)]. This simulation thus validates both the overall method and our theoretical predictions.

This numerical simulation can also be used to outline the limits of our approach, in particular in media with a heterogeneous distribution of scatterers. To that aim, we now consider a two-layered medium, each layer having a different density of scatterers. Two cases are investigated: dilute-to-dense [Fig. S3(a)] and dense-to dilute [Fig. S3(c)]. The
This analytical expression is confronted to the result of the Monte Carlo simulation in Fig. 1(b). A perfect agreement is obtained for $t < \tau_s$ but, as expected, this analytical expression cannot grasp the evolution of the backscattered intensity for $t > \tau_s$.

Using Eqs. (2) and (4), an expression can be derived for the single scattering rate, $\rho_s(t) = I_s(t)/I(t)$, in the short-time limit:

$$ \rho_s(t) \sim_{t<\tau_s} \exp(-3ct/4\ell_s). \qquad (S40) $$

This is the fundamental result on which the accompanying paper is based. This exponential decay of the single scattering rate is confirmed by the result of the Monte Carlo simulation displayed in Fig. 1(c). Again, an excellent agreement is found between our analytical prediction and the solution of the radiative transfer equation for $t < \tau_s$. Note that the decay of $\rho_s(t)$ is not strictly identical in a 3D geometry but an exponential decay of $\exp(-0.55ct/\ell_s)$ fits pretty well the result of our Monte Carlo simulation for $t < \tau_s$ [Fig. S1(b)].

## DIFFUSE INTENSITY IN THE LONG-TIME LIMIT

In a scattering medium ($\ell_s \ll \ell_a$), the specific intensity $P(\mathbf{r}, t, \mathbf{u})$ can be decomposed as the sum of an isotropic and an anisotropic term [51]:

$$ P(\mathbf{r}, t, \mathbf{u}) = \frac{\phi(\mathbf{r}, t)}{2\pi} + \frac{1}{\pi}\mathbf{J}(\mathbf{r}, t)\cdot\mathbf{u} \qquad (S41) $$

with $\phi(\mathbf{r}, t)$, the fluence rate (or intensity), and $\mathbf{J}(\mathbf{r}, t)$, the energy flux. In the long time limit ($t \gg \tau_s$), the radiative transfer equation then leads to a much simpler equation for $\phi(\mathbf{r}, t)$:

$$ \frac{1}{c} \frac{\partial}{\partial t} \phi(\mathbf{r}, t) - D \Delta \phi(\mathbf{r}, t) + \frac{1}{\ell_a} \phi(\mathbf{r}, t) = S(\mathbf{r}, t) \qquad (S42) $$

with $D$, the diffusion coefficient. In a 2D geometry and for isotropic scattering, $D = c\ell_s/2$.

For an infinite medium, the Green's function $F$ of the diffusion equation [$S(\mathbf{r}, t) = \delta(\mathbf{r})\delta(t)$] can be expressed as follows:

$$ F(\mathbf{r}, t) = \frac{1}{4\pi Dt} \exp(-ct/\ell_a) \exp(-r^2/4Dt). \qquad (S43) $$

In a semi-infinite medium, the cancellation of $\phi(\mathbf{r}, t)$ at the boundary ($z = 0$) can be met by adding a negative or image source of energy to the infinite medium problem. The fluence rate per incident wave packet can then be written as the sum of contributions from two sources
<figure><img src="image_7.png" /></figure>

FIG. S2. Single scattering rate for homogeneous disorder. (a) Cylinders uniformly distributed in a 15×13.5 cm² area with ca concentration n = 0.2 mm⁻². ℓs is estimated from ISA: ℓs = 194 mm. (b) Depth-evolution of the single scattering rate ρ̂s (blue dots) compared with the theoretical prediction [Eq. (5)] using the ls value derived under ISA. (c) Cylinders uniformly distributed in a 15×13.5 cm² area with number density n = 1 mm⁻². The scattering mean free path estimated from ISA is ls = 39 mm. (d) Depth-evolution of the single scattering rate ρ̂s (blue dots) compared with the theoretical prediction [Eq. (5)] using the ls value derived under ISA. The vertical dashed line accounts for z = cτs/2, the boundary between the short and long time regimes.

scatterer concentration in the dilute layer is n = 0.04 mm⁻² (ls ~ 970 mm), while n = 0.2
mm⁻² in the dense layer (ls ~ 194 mm). The depth evolution of the single scattering
rate measured in each system is displayed in Figs. S3(b) and (d) respectively. As discussed
in the main text, propagating from a dilute to a dense region does not affect the single
scattering ratio that exhibits a clear change of exponential decay at the interface between
the two regions [Fig. S3(d)]. On the contrary, in the dense-to-dilute case, ps(t) completely
falls down right after the interface because the multiple scattering coming from the dense
at $z = z_0$ and $z = -z_0$, with $z_0 = 2\ell_s/3$, the extrapolation length [45]. The corresponding Green's function $F(\mathbf{r}, t)$ is thus given by:

$$F(\mathbf{r}, t) = \frac{1}{4\pi Dt} \exp(-ct/\ell_a) \left[ \exp\left(-\frac{x^2 + (z - z_0)^2}{4Dt}\right) + \exp\left(-\frac{x^2 + (z + z_0)^2}{4Dt}\right) \right]. \quad (S44)$$

The associated energy flux $\mathbf{J}(\mathbf{r}, t)$ at the medium surface can be deduced from Fick's law:

$$\mathbf{J}(\mathbf{r}, t) = -(D/c)\nabla F(\mathbf{r}, t)|_{z=0}. \quad (S45)$$

The diffuse intensity, $P_d(\mathbf{r}, t) = |\mathbf{J}(\mathbf{r}, t)|$, can then be deduced:

$$P_d(\mathbf{r}, t) = \frac{\exp(-ct/\ell_a) z_0}{2\pi Dt} \exp\left(-\frac{x^2 + z_0^2}{4Dt}\right) \quad (S46)$$

The total diffuse reflectance, $I_d(t) = \int_{-\infty}^{\infty} dx P_d(x, z, t)$, is finally given by:

$$I_d(t) = \frac{\exp(-ct/\ell_a) z_0}{\sqrt{\pi Dt}} \exp\left(-\frac{z_0^2}{4Dt}\right) \quad (S47)$$

In the long-time limit, $4Dt \gg z_0^2$. $I_d(t)$ thus simplifies into

$$I_d(t) = \frac{\exp(-ct/\ell_a) z_0}{\sqrt{\pi Dt}} \quad (S48)$$

Interestingly, this expression of $I_d(t)$ is identical to the one previously derived by Patterson et al. in 3D [44]. This expression can be used to fit the result of the Monte Carlo simulation for $t > \tau_s$ both in 2D [Fig. 1(b)] and 3D [Fig. S1(a)]. If absorption is negligible, a power law decay as $t^{-3/2}$ is thus predicted for the scattered intensity in the diffusive limit.

Using Eqs. (2) and (S48), an expression can be derived for the single scattering rate, $\rho_s(t) \sim I_s(t)/I_d(t)$, in the long-time limit:

$$\rho_s(t) \sim \frac{3}{8} \sqrt{\frac{\pi}{2}} (ct/\ell_s)^{3/2} \exp(-ct/\ell_s). \quad (S49)$$

This diffusion approximation predicts correctly the evolution of the single scattering rate in the long-time limit (Fig. 1c). As in the short-time limit, the single scattering rate can give access, in principle, to a measurement of $\ell_s$ independent from $\ell_a$. Nevertheless, the scaling with $\ell_s$ is more complicated than a simple exponential decay.

## NUMERICAL VALIDATION

The full method is now validated by means of a simulation of the wave equation. The scattering medium is a random collection of fluid scatterers ($c_s = 2500$ m/s, radius 0.1 mm)
<figure><img src="image_8.png" /></figure>

FIG. S3. Single scattering rate for an heterogeneous distribution of disorder. (a) Two-layered system with a dilute-to-dense configuration: $n = 0.04 \text{ mm}^{-2}$ ($l_s = 970$ mm) in the first half of the medium, and dense $n = 0.2 \text{ mm}^{-2}$ ($l_s = 194$ mm) in the second half. (b) Depth evolution of the single scattering rate $\hat{\rho}_s$ (blue dots) compared with the theoretical prediction [Eq. (5)] with the $l_s$ value derived under the ISA. (c) Dense-to-dilute system with the same parameters as in (a). (d) Same as in (b) but for the dense-to-dilute system displayed in (c).

region still dominates at times corresponding to ballistic depths located inside the dilute region. Hence, a positive gradient of scattering can be nicely resolved by our method but not a decreasing one. In the latter case, a more subtle inversion of $\rho_s(t)$ is needed to retrieve the depth evolution of the scattering mean free path $l_s$.

<footnote>* Present address: Université de Rennes, CNRS, IETR (Institut d'Électronique et des Technologies du numéRique), UMR-6164, Rennes, France</footnote>

<footnote>† alexandre.aubry@espci.fr</footnote>
embedded in water ($c_0 = 1480 \text{ m/s}$) [Fig. S2], yielding an individual scattering cross-section $\sigma = 2.58 \times 10^{-2} \text{ mm}$ [42]. It is assumed that the density contrast between the cylinders and surrounding medium is negligible, so that heterogeneity comes from the compressibility contrast; only longitudinal (pressure) waves are taken into account. The associated reflection matrix is computed over the 1.3 – 1.7 MHz-frequency bandwidth by means of a Born series already described in a previous study [42]. The simulated ultrasound probe is an array of $N = 64$ transducers with an inter-element distance $p = 0.5 \text{ mm}$. The array of transducers is placed at 140 mm from the sample surface.

In the first system [Fig. S2(a)], the concentration of scatterers is $n = 0.2 \text{ mm}^{-2}$. Under the independent scattering approximation (ISA) [52], the scattering mean free path is given by $\ell_s \sim (n\sigma)^{-1}$, that is to say $\ell_s \sim 194 \text{ mm}$ in the present case. Under the independent scattering approximation (ISA) that holds in a dilute regime ($k\ell_s \gg 1$, with $k$ the wave number) [52], the scattering mean free path $\ell_s \sim (n\sigma)^{-1}$, that is to say $\ell_s \sim 194 \text{ mm}$ in the present case. ISA is therefore justified ($k\ell_s \sim 1250$) and the value of $\ell_s$ is therefore correct.

The simulated reflection matrix is projected onto the single scattering basis following the method described in Sec. II. Nevertheless, the bias $\rho_b$ (Eq. S12) is here not negligible because the lateral resolution $\delta x$ is large (far-field configuration). The bias $\rho_b$ is therefore determined from the value of $\rho_s$ at echo times corresponding to ballistic depths beyond the medium thickness. An unbiased estimator is then built using Eq. S13.

The resulting single scattering ratio $\hat{\rho}'_s(t)$ is displayed in Fig. S2(b). A quantitative agreement is found between $\hat{\rho}'_s(t)$ and its theoretical scaling in $\exp(-3ct/4\ell_s)$ in the short time limit. The same observation holds for a more concentrated scattering medium [Fig. S2(c)] for which $n = 1 \text{ mm}^{-2}$ and $\ell_s \sim 39 \text{ mm}$. Again, ISA holds ($k\ell_s \sim 250$) and our theoretical prediction [Eq. (5)] perfectly fits the single scattering rate $\hat{\rho}'_s(t)$ for $t < \tau_s$ [see Fig. S2(d)]. This simulation thus validates both the overall method and our theoretical prediction for $t < \tau_s$. Note, however, that beyond $\tau_s$ (predominant multiple scattering), the estimated single scattering ratio shows strong fluctuations and more spatial averaging would be needed to fit $\rho_s(t)$ with the diffusive prediction of Eq. S49 [black line in Fig. S2(d)].

This numerical simulation can also be used to outline the limits of our approach, in particular in media with a heterogeneous distribution of scatterers. To that aim, we now consider a two-layered medium, each layer having a different density of scatterers. Two cases are investigated: dilute-to-dense [Fig. S3(a)] and dense-to dilute [Fig. S3(c)]. The scatterer
[1] E. Abrahams, P. W. Anderson, D. C. Licciardello, and T. V. Ramakrishnan, Scaling theory of localization: Absence of quantum diffusion in two dimensions, Phys. Rev. Lett. 42, 673 (1979).

[2] M. P. V. Albada and A. Lagendijk, Observation of weak localization of light in a random medium, Phys. Rev. Lett. 55, 2692 (1985).

[3] P. A. Lee and A. D. Stone, Universal conductance fluctuations in metals, Phys. Rev. Lett. 55, 1622 (1985).

[4] A. Tourin, A. Derode, P. Roux, B. A. van Tiggelen, and M. Fink, Time-dependent coherent backscattering of acoustic waves, Phys. Rev. Lett. 79, 3637 (1997).

[5] F. Scheffold and G. Maret, Universal conductance fluctuations of light, Phys. Rev. Lett. 81, 5800 (1998).

[6] E. Larose, L. Margerin, B. A. van Tiggelen, and M. Campillo, Weak localization of seismic waves, Phys. Rev. Lett. 93, 048501 (2004).

[7] H. Hu, A. Strybulevych, J. H. Page, S. E. Skipetrov, and B. A. van Tiggelen, Localization of ultrasound in a three-dimensional elastic network, Nat. Phys. 4, 945 (2008).

[8] I. M. Vellekoop and A. P. Mosk, Universal optimal transmission of light through disordered materials, Phys. Rev. Lett. 101, 120601 (2008).

[9] S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, 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).

[10] F. Jendrzejewski, A. Bernard, K. Müller, P. Cheinet, V. Josse, M. Piraud, L. Pezzé, L. Sanchez-Palencia, A. Aspect, and P. Bouyer, Three-dimensional localization of ultracold atoms in an optical disordered potential, Nat. Phys. 8, 398 (2012).

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

[12] C. W. Hsu, S. F. Liew, A. Goetschy, H. Cao, and A. Douglas Stone, Correlation-enhanced control of wave focusing in disordered media, Nat. Phys. 13, 497 (2017).

[13] Z. Shi and A. Z. Genack, Diffusion in translucent media, Nat. Commun. 9, 1862 (2018).

[14] M. Horodynski, M. Kühmayer, C. Ferise, S. Rotter, and M. Davy, Anti-reflection structure for perfect transmission through complex media, Nature 607, 281 (2022).

[15] A. Yamilov, S. E. Skipetrov, T. W. Hughes, M. Minkov, Z. Yu, and H. Cao, Anderson local-
<figure><img src="image_7.png" /></figure>

FIG. S2. Single scattering ratio for homogeneous disorder. (a) Cylinders uniformly distributed in a 15×13.5 cm² area with number density $n = 0.2 \text{ mm}^{-2}$. $l_s$ is estimated from ISA: $l_s = 194$ mm. (b) Depth-evolution of the unbiased single scattering ratio $\hat{\rho}'_s$ (blue dots) compared with the theoretical prediction [Eq. (5)] using the $l_s$ value derived under ISA. (c) Cylinders uniformly distributed in a 15×13.5 cm² area with number density $n = 1 \text{ mm}^{-2}$. The scattering mean free path estimated from ISA is $l_s = 39$ mm. (d) Depth-evolution of the unbiased single scattering ratio $\hat{\rho}'_s$ (blue dots) compared with our theoretical prediction for $t < \tau_s$ [red line, Eq. (5)] and the diffusive result [black line, Eq. S49] using the $l_s$ value derived under ISA. The vertical dashed line accounts for $z = c\tau_s/2$, the boundary between the short and long time regimes.

concentration in the dilute layer is $n = 0.04 \text{ mm}^{-2}$ ($l_s \sim 970$ mm), while $n = 0.2 \text{ mm}^{-2}$ in the dense layer ($l_s \sim 194$ mm). The depth evolution of the unbiased single scattering ratio $\hat{\rho}'_s(t)$ measured in each system is displayed in Figs. S3(b) and (d), respectively. As discussed in the main text, propagating from a dilute to a dense region does not affect the single scattering ratio that exhibits a clear change of exponential decay at the interface between the two regions [Fig. S3(d)]. On the contrary, in the dense-to-dilute case, $\hat{\rho}'_s(t)$ completely
ization of electromagnetic waves in three dimensions, Nat. Phys. **19**, 1308 (2023).

[16] U. Frisch, Probabilistic methods in applied mathematics (Academic Press Inc., New York, 1968) Chap. Wave propagation in random media, pp. 75–198.

[17] J. A. Turner and P. Anugonda, Scattering of elastic waves in heterogeneous media with local isotropy, J. Acoust. Soc. Am. **109**, 1787 (2001).

[18] S. M. Rytov, K. Y. A., and V. I. Tatarskii, *Principles of Statistical Radiophysics Vol. 4. Wave Propagation Through Random Media* (Springer, Berlin, 1989).

[19] I. Baydoun, D. Baresch, R. Pierrat, and A. Derode, Radiative transfer of acoustic waves in continuous complex media: Beyond the helmholtz equation, Phys. Rev. E **94**, 053005 (2016).

[20] B. J. Oosterveld, J. M. Thijssen, P. C. Hartman, R. L. Romijn, and G. J. E. Rosenbusch, Ultrasound attenuation and texture analysis of diffuse liver disease: Methods and preliminary results, Phys. Med. Biol. **36**, 1039 (1991).

[21] A. Aubry, A. Derode, and F. Padilla, Local measurements of the diffusion constant in multiple scattering media: Application to human trabecular bone imaging, Appl. Phys. Lett. **92**, 1240101 (2008).

[22] A. Aubry and A. Derode, Multiple scattering of ultrasound in weakly inhomogeneous media: Application to human soft tissues, J. Acoust. Soc. Am. **129**, 225 (2011).

[23] M. Sasso, M. Beaugrand, V. de Ledinghen, C. Douvin, P. Marcellin, R. Poupon, L. Sandrin, and V. Miette, Controlled attenuation parameter (cap): A novel vcte™ guided ultrasonic attenuation measurement for the evaluation of hepatic steatosis: Preliminary study and validation in a cohort of patients with chronic liver disease from various causes, Ultrasound Med. Biol. **36**, 1825 (2010).

[24] K. Mohanty, J. Blackwell, T. Egan, and M. Muller, Characterization of the lung parenchyma using ultrasound multiple scattering, Ultrasound Med. Biol. **43**, 993 (2017).

[25] Q. W. Guerrero, H. Feltovich, I. M. Rosado-Mendez, L. C. Carlson, G. Li, and T. J. Hall, Anisotropy and spatial heterogeneity in quantitative ultrasound parameters: Relevance to the study of the human cervix, Ultrasound Med. Biol. **44**, 1493 (2018).

[26] T. Durduran, R. Choe, W. B. Baker, and A. G. Yodh, Diffuse optics for tissue monitoring and tomography, Rep. Prog. Phys. **73**, 076701 (2010).

[27] S. Hirsekorn, The scattering of ultrasonic waves by polycrystals, J. Acoust. Soc. Am. **72**, 1021 (1982).
<figure><img src="image_8.png" /><figcaption>FIG. S3. Single scattering ratio for an heterogeneous distribution of disorder. (a) Two-layered system with a dilute-to-dense configuration: $n = 0.04 \text{ mm}^{-2}$ ($l_s = 970$ mm) in the first half of the medium, and dense $n = 0.2 \text{ mm}^{-2}$ ($l_s = 194$ mm) in the second half. (b) Depth evolution of the single scattering ratio $\rho'_s$ (blue dots) compared with the theoretical prediction [Eq. (5)] with the $l_s$ values derived under the ISA for each scattering layer. (c) Dense-to-dilute system with the same parameters as in (a). (d) Same as in (b) but for the dense-to-dilute system displayed in (c).</figcaption></figure>

falls down right after the interface because the multiple scattering coming from the dense region still dominates at times corresponding to ballistic depths located inside the dilute region. Hence, a positive gradient of scattering can be nicely resolved by our method but not a negative one. In the latter case, a more subtle inversion of $\rho_s(t)$ is needed to retrieve the depth evolution of the scattering mean free path $l_s$.
[28] F. J. Margetan, R. B. Thompson, and I. Yalda-Mooshabad, Backscattered microstructural noise in ultrasonic toneburst inspections, *J. Nondestruct. Eval.* **13**, 111 (1994).

[29] P. D. Panetta, F. J. Margetan, I. Yalda, and R. B. Thompson, Ultrasonic attenuation measurements in jet-engine titanium alloys, in *Review of Progress in Quantitative Nondestructive Evaluation* (Springer US, 1996) pp. 1525–1532.

[30] R. B. Thompson, F. Margetan, P. Haldipur, L. Yu, A. Li, P. Panetta, and H. Wasan, Scattering of elastic waves in simple and complex polycrystals, *Wave Motion* **45**, 655 (2008).

[31] S. Shahjahan, F. Rupin, A. Aubry, B. Chassignole, T. Fouquet, and A. Derode, Comparison between experimental and 2-d numerical studies of multiple scattering in Inconel600® by means of array probes, *Ultrasonics* **54**, 358 (2014).

[32] A. Baelde, J. Laurent, P. Millien, R. Coulette, W. B. Khalifa, F. Jenson, F. Sun, M. Fink, and C. Prada, Effect of microstructural elongation on backscattered field: Intensity measurement and multiple scattering estimation with a linear transducer array, *Ultrasonics* **82**, 379 (2018).

[33] H. Sato, M. C. Fehler, and T. Maeda, *Seismic Wave Propagation and Scattering in the Heterogeneous Earth : Second Edition* (Springer-Verlag, Berlin, 2012).

[34] J. Chaput, M. Campillo, R. C. Aster, P. Roux, P. R. Kyle, H. Knox, and P. Czoski, Multiple scattering from icequakes at Erebus volcano, Antarctica: Implications for imaging at glaciated volcanoes, *J. Geophys. Res.* **120**, 1129 (2015).

[35] J. Mayor, P. Traversa, M. Calvet, and L. Margerin, Tomography of crustal seismic attenuation in metropolitan france: Implications for seismicity analysis, *Bull. Earthquake Eng.* **16**, 2195 (2018).

[36] K. Onodera, T. Maeda, K. Nishida, T. Kawamura, L. Margerin, S. Menina, P. Lognonné, and W. B. Banerdt, Seismic scattering and absorption properties of Mars estimated through coda analysis on a long-period surface wave of S1222a Marsquake, *Geophys. Res. Lett.* **50**, e2022GL102716 (2023).

[37] G. Bal and K. Ren, Transport-based imaging in random media, *SIAM J. Appl. Math.* **68**, 1738 (2008).

[38] S. R. Arridge and J. C. Schotland, Optical tomography: forward and inverse problems, *Inverse Problems* **25**, 123010 (2009).

[39] J. H. Page, P. Sheng, H. P. Schriemer, I. Jones, X. Jing, and D. A. Weitz, Group Velocity in Strongly Scattering Media, *Science* **271**, 634 (1996).
* Present address: Université de Rennes, CNRS, IETR (Institut d'Électronique et des Technologies du numéRique), UMR-6164, Rennes, France

† alexandre.aubry@espci.fr

[1] E. Abrahams, P. W. Anderson, D. C. Licciardello, and T. V. Ramakrishnan, Scaling theory of localization: Absence of quantum diffusion in two dimensions, Phys. Rev. Lett. 42, 673 (1979).

[2] M. P. Van Albada and A. Lagendijk, Observation of weak localization of light in a random medium, Phys. Rev. Lett. 55, 2692 (1985).

[3] P. A. Lee and A. D. Stone, Universal conductance fluctuations in metals, Phys. Rev. Lett. 55, 1622 (1985).

[4] A. Tourin, A. Derode, P. Roux, B. A. van Tiggelen, and M. Fink, Time-dependent coherent backscattering of acoustic waves, Phys. Rev. Lett. 79, 3637 (1997).

[5] F. Scheffold and G. Maret, Universal conductance fluctuations of light, Phys. Rev. Lett. 81, 5800 (1998).

[6] E. Larose, L. Margerin, B. A. van Tiggelen, and M. Campillo, Weak localization of seismic waves, Phys. Rev. Lett. 93, 048501 (2004).

[7] H. Hu, A. Strybulevych, J. H. Page, S. E. Skipetrov, and B. A. van Tiggelen, Localization of ultrasound in a three-dimensional elastic network, Nat. Phys. 4, 945 (2008).

[8] I. M. Vellekoop and A. P. Mosk, Universal optimal transmission of light through disordered materials, Phys. Rev. Lett. 101, 120601 (2008).

[9] S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, 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).

[10] F. Jendrzejewski, A. Bernard, K. Müller, P. Cheinet, V. Josse, M. Piraud, L. Pezzé, L. Sanchez-Palencia, A. Aspect, and P. Bouyer, Three-dimensional localization of ultracold atoms in an optical disordered potential, Nat. Phys. 8, 398 (2012).

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

[12] C. W. Hsu, S. F. Liew, A. Goetschy, H. Cao, and A. Douglas Stone, Correlation-enhanced
[40] see supplementary material, URL_will_be_inserted_by_publisher.

[41] J. C. J. Paasschens, Solution of the time-dependent Boltzmann equation, Phys. Rev. E 56, 1135 (1997).

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

[43] A. Velichko, Quantification of the Effect of Multiple Scattering on Array Imaging Performance, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 67, 92 (2020).

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

[45] C. Brütt, A. Aubry, B. Gérardin, A. Derode, and C. Prada, Weight of single and recurrent scattering in the reflection matrix of complex media, Phys. Rev. E 106, 025001 (2022).

[46] S. Audiere, M. Clet, M. Sasso, L. Sandrin, and V. Miette, Influence of heterogeneities on ultrasound attenuation for liver steatosis evaluation (CAPTM): Relevance of a liver guidance tool, in 2013 IEEE Int. Ultrason. Symp. IUS (IEEE, Prague, Czech Republic, 2013) pp. 401–404.

[47] M. S. Patterson, B. Chance, and B. C. Wilson, Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties, Appl. Opt. 28, 2331 (1989).

[48] J. X. Zhu, D. J. Pine, and D. A. Weitz, Internal reflection of diffusive light in random media, Phys. Rev. A 44, 3948 (1991).

[49] K. Suzuki, N. Hayashi, Y. Sasaki, M. Kono, A. Kasahara, H. Fusamoto, Y. Imai, and T. Kamada, Dependence of ultrasonic attenuation of liver on pathologic fat and fibrosis: Examination with experimental fatty liver and liver fibrosis models, Ultrasound Med. Biol. 18, 657 (1992).

[50] T. Karlas et al., Individual patient data meta-analysis of controlled attenuation parameter (cap) technology for assessing steatosis, J. Hepatol. 66, 1022 (2017).

[51] S. Strauss, E. Gavish, P. Gottlieb, and L. Katsnelson, Interobserver and intraobserver variability in the sonographic assessment of fatty liver, Am. J. Roentgenol. 189, W320–W323 (2007).

[52] S. Yoon, M. Kim, M. Jang, Y. Choi, W. Choi, S. Kang, and W. Choi, Deep optical imaging within complex scattering media, Nat. Rev. Phys. 2, 141 (2020).
control of wave focusing in disordered media, Nat. Phys. 13, 497 (2017).

[13] Z. Shi and A. Z. Genack, Diffusion in translucent media, Nat. Commun. 9, 1862 (2018).

[14] M. Horodynski, M. Kühmayer, C. Ferise, S. Rotter, and M. Davy, Anti-reflection structure for perfect transmission through complex media, Nature 607, 281 (2022).

[15] A. Yamilov, S. E. Skipetrov, T. W. Hughes, M. Minkov, Z. Yu, and H. Cao, Anderson localization of electromagnetic waves in three dimensions, Nat. Phys. 19, 1308 (2023).

[16] K. Suzuki, N. Hayashi, Y. Sasaki, M. Kono, A. Kasahara, H. Fusamoto, Y. Imai, and T. Kamada, Dependence of ultrasonic attenuation of liver on pathologic fat and fibrosis: Examination with experimental fatty liver and liver fibrosis models, Ultrasound Med. Biol. 18, 657 (1992).

[17] M. Sasso, M. Beaugrand, V. de Ledinghen, C. Douvin, P. Marcellin, R. Poupon, L. Sandrin, and V. Miette, Controlled attenuation parameter (CAP): A novel VCTE™ guided ultrasonic attenuation measurement for the evaluation of hepatic steatosis: Preliminary study and validation in a cohort of patients with chronic liver disease from various causes, Ultrasound Med. Biol. 36, 1825 (2010).

[18] T. Karlas et al., Individual patient data meta-analysis of controlled attenuation parameter (CAP) technology for assessing steatosis, J. Hepatol. 66, 1022 (2017).

[19] A. Aubry, A. Derode, and F. Padilla, Local measurements of the diffusion constant in multiple scattering media: Application to human trabecular bone imaging, Appl. Phys. Lett. 92, 1240101 (2008).

[20] Y. Karbalaeisadegh and M. Muller, Ultrasound scattering in cortical bone, in Bone Quantitative Ultrasound (Springer International Publishing, 2022) pp. 177-196.

[21] K. Mohanty, J. Blackwell, T. Egan, and M. Muller, Characterization of the lung parenchyma using ultrasound multiple scattering, Ultrasound Med. Biol. 43, 993 (2017).

[22] D. A. Boas, A. M. Dale, and M. A. Franceschini, Diffuse optical imaging of brain activation: Approaches to optimizing image sensitivity, resolution, and accuracy, NeuroImage 23, S275 (2004).

[23] R. Bocheux, P. Pernot, V. Borderie, K. Plamann, and K. Irsch, Quantitative measures of corneal transparency, derived from objective analysis of depth-resolved corneal images, demonstrated with full-field optical coherence tomographic microscopy, PLOS ONE 14, e0221707 (2019).
[53] F. Bureau, J. Robin, A. Le Ber, W. Lambert, M. Fink, and A. Aubry, Three-dimensional ultrasound matrix imaging, *Nature Commun.* **14**, 6793 (2023).

[54] U. Najar, V. Barolle, P. Balondrade, M. Fink, A. C. Boccara, M. Fink, and A. Aubry, *Non-invasive retrieval of the time-gated transmission matrix for optical imaging deep inside a multiple scattering medium*, arXiv:2303.06119 (2023).

[55] E. Giraudat, A. Burtin, A. L. Ber, M. Fink, J.-C. Komorowski, and A. Aubry, *Unveiling the deep plumbing system of a volcano by a reflection matrix analysis of seismic noise*, arXiv:2311.01296 (2023).

[56] L. V. Wang and H.-I. Wu, Biomedical optics: Principles and imaging (John Wiley & Sons, Inc., 2009) Chap. 5. Radiative Transfer Equation and Diffusion Theory, pp. 83–118.

[57] P. Sheng, *Introduction to Wave Scattering, Localization, and Mesoscopic Phenomena*, 2nd ed., Springer Series in Materials Science No. 88 (Springer, Berlin ; New York, 2006).
[24] H. Sato, M. C. Fehler, and T. Maeda, *Seismic Wave Propagation and Scattering in the Heterogeneous Earth: Second Edition* (Springer-Verlag, Berlin, 2012).

[25] J. Mayor, P. Traversa, M. Calvet, and L. Margerin, Tomography of crustal seismic attenuation in metropolitan france: Implications for seismicity analysis, *Bull. Earthquake Eng.* **16**, 2195 (2018).

[26] J. Li and S. Rokhlin, Propagation and scattering of ultrasonic waves in polycrystals with arbitrary crystallite and macroscopic texture symmetries, *Wave Motion* **58**, 145 (2015).

[27] S. Shahjahan, F. Rupin, A. Aubry, B. Chassignole, T. Fouquet, and A. Derode, Comparison between experimental and 2-d numerical studies of multiple scattering in Inconel600® by means of array probes, *Ultrasonics* **54**, 358 (2014).

[28] G. Bal and K. Ren, Transport-based imaging in random media, *SIAM J. Appl. Math.* **68**, 1738 (2008).

[29] S. R. Arridge and J. C. Schotland, Optical tomography: forward and inverse problems, *Inverse Problems* **25**, 123010 (2009).

[30] T. Durduran, R. Choe, W. B. Baker, and A. G. Yodh, Diffuse optics for tissue monitoring and tomography, *Rep. Prog. Phys.* **73**, 076701 (2010).

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

[32] See supplementary material (URL) for method description, theoretical demonstration and numerical validation, which also includes Refs. [51, 52].

[33] S. Yoon, M. Kim, M. Jang, Y. Choi, W. Choi, S. Kang, and W. Choi, Deep optical imaging within complex scattering media, *Nat. Rev. Phys.* **2**, 141 (2020).

[34] A. Badon, V. Barolle, K. Irsch, A. C. Boccara, M. Fink, and A. Aubry, Distortion matrix concept for deep optical imaging in scattering media, *Sci. Adv.* **6**, eaay7170 (2020).

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

[36] A. Velichko, Quantification of the Effect of Multiple Scattering on Array Imaging Performance, *IEEE Trans. Ultrason.*, *Ferroelect.*, *Freq. Contr.* **67**, 92 (2020).

[37] T. Blondel, J. Chaput, A. Derode, M. Campillo, and A. Aubry, Matrix approach of seismic imaging: Application to the Erebus volcano, Antarctica, *J. Geophys. Res.: Solid Earth* **123**, 10936 (2018).
[38] R. Touma, A. Aubry, Y. Ben-Zion, and M. Campillo, Distribution of seismic scatterers in the san jacinto fault zone, southeast of anza, california, based on passive matrix imaging, Earth Planet. Sci. Lett. 578, 117304 (2022).

[39] J. Sol, H. Prod'homme, L. Le Magoarou, and P. del Hougne, Experimentally realized physical-model-based frugal wave control in metasurface-programmable complex media, Nat. Commun. 15, 2841 (2024).

[40] J. C. J. Paasschens, Solution of the time-dependent Boltzmann equation, Phys. Rev. E 56, 1135 (1997).

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

[42] C. Brütt, A. Aubry, B. Gérardin, A. Derode, and C. Prada, Weight of single and recurrent scattering in the reflection matrix of complex media, Phys. Rev. E 106, 025001 (2022).

[43] S. Audiere, M. Clet, M. Sasso, L. Sandrin, and V. Miette, Influence of heterogeneities on ultrasound attenuation for liver steatosis evaluation (CAPTM): Relevance of a liver guidance tool, in 2013 IEEE Int. Ultrason. Symp. IUS (IEEE, Prague, Czech Republic, 2013) pp. 401–404.

[44] M. S. Patterson, B. Chance, and B. C. Wilson, Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties, Appl. Opt. 28, 2331 (1989).

[45] J. X. Zhu, D. J. Pine, and D. A. Weitz, Internal reflection of diffusive light in random media, Phys. Rev. A 44, 3948 (1991).

[46] A. Aubry and A. Derode, Multiple scattering of ultrasound in weakly inhomogeneous media: Application to human soft tissues, J. Acoust. Soc. Am. 129, 225 (2011).

[47] U. Najar, V. Barolle, P. Balondrade, M. Fink, A. C. Boccara, M. Fink, and A. Aubry, Harnessing forward multiple scattering for optical imaging deep inside an opaque medium, Nat. Commun. (in press) 15, 7349 (2024).

[48] P. Balondrade, V. Barolle, N. Guigui, E. Auriant, N. Rougier, C. Boccara, M. Fink, and A. Aubry, Multi-spectral reflection matrix for ultra-fast 3D label-free microscopy, Nat. Photon. 10.1038/s41566-024-01479-y (2024).

[49] F. Bureau, J. Robin, A. Le Ber, W. Lambert, M. Fink, and A. Aubry, Three-dimensional ultrasound matrix imaging, Nat. Commun. 14, 6793 (2023).
[50] E. Giraudat, A. Burtin, A. L. Ber, M. Fink, J.-C. Komorowski, and A. Aubry, Matrix imaging as a tool for high-resolution monitoring of deep volcanic plumbing systems with seismic noise, Commun. Earth Environ. 10.1038/s43247-024-01659-2 (2024).

[51] L. V. Wang and H.-I. Wu, Biomedical optics: Principles and imaging (John Wiley & Sons, Inc., 2009) Chap. 5. Radiative Transfer Equation and Diffusion Theory, pp. 83–118.

[52] P. Sheng, *Introduction to Wave Scattering, Localization, and Mesoscopic Phenomena*, 2nd ed., Springer Series in Materials Science No. 88 (Springer, Berlin; New York, 2006).