Ching Hsueh

Carmen Lucía

Moraila Martínez
email: cmoraila@ugr.es

Frédéric Doumenc
email: doumenc@fast.u-psud.fr

Miguel Rodríguez-Valverde

Béatrice Guerrier
email: guerrier@fast.u-psud.fr

Miguel A Rodríguez- Valverde

Self-assembly in drying complex fluid at low capillary number

Keywords: silica nano-particles, stick-slip motion, dip-coating, pattern formation, drying

à 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.

Introduction

Surface coating technologies used nowadays to fabricate nano-scale patterning is a main issue to investigate. Different methods are used to make patterned deposits: from the top-down point of view, photolithography and electron-beam lithography are well-developed. However, the parallelization and 2-D limitations are actual challenges. From the bottom-up point of view, selfassembly by evaporation is one of the promising ideas to fabricate nano-scale patterns [1][2][3]. Droplet drying has been widely investigated in fundamental studies [4][5] and for various applications, like biochemical ones [6]. Despite the simplicity of this geometry, its physical description remains quite complicated. Indeed the drying of a droplet is a non stationary phenomenon since its volume decreases while the non volatile solute concentration increases. In dip-coating configurations, also widely used in industrial coating processes, the concentration of the liquid bath remains constant so that a stationary phenomenon can be observed. Moreover, velocity can be controlled independently from evaporation. This configuration can then help to better describe and control the different phenomena that are exhibited in patterning.

Several regimes are encountered in dip-coating, depending on the capillary number. At high enough capillary number, the viscous forces dominate and the coating forms a continuous layer whose thickness is proportional to the substrate speed to the power 2/3 (Landau-Levich regime [7,8]). At lower speed (low capillary number) another regime appears, dominated by evaporation. In this evaporative regime, the thickness of the deposit is inversely proportional to the substrate velocity. This scaling was first reported by Dimitrov and Nagayama [9] in 1996 for experiments with colloidal suspensions. These authors also proposed a simple model based on mass balances of both the solution and the solute, which shows that the particles flux inside the meniscus is mainly driven by evaporation. Experiments were extended more recently to other systems [10][11][12] where the same scaling was observed. The existence of a cross-over between an evaporative regime and a Landau-Levich regime was theoretically predicted by Berteloot et al. [13] (despite the fact that their model did not give the right scaling in the evaporative regime), and experimentally observed by Le Berre et al [10] and Faustini et al [12]. For some configurations a stick-slip motion of the contact line is obtained and leads to regular stripe patterns. This phenomenon has been first mentioned by Adachi and co-authors [14] for colloids and sessile droplet geometry. One remarkable aspect of stick-slip phenomenon is its universality, in the sense that it can be observed for very different systems such as small molecules, colloidal suspensions, polymer solutions, etc (cf for instance [10,[15][16][17][18] and the recent review by Han and Lin [19]). However the conditions required for strong stick-slip depend on both the process parameters and the physico-chemical properties of the solution, in a non trivial way. The purpose of this paper is to illustrate this dependence. The first study concerns silica suspensions. The second polymer solutions.

Experiment

Setup:

The setup principle is similar to a reservoir which contains coating materials solution is pumped out from the reservoir sandwich device made of two paralleled glas immersed into a reservoir which co at the top of the glass slides and close to the edges. spontaneous capillary rise. Th significant capillary rise while being able to ensure a good flow rate of the pumping system controls The setup is put inside a chamber of volume 6 by a PID system. In all the experiments presented in this paper, the temperature and the humidity set-point is 30 channel carefully screwed on the t uniform evaporation rate but a negligible meniscus. The gap between the two plates being open, the air escapes through the lateral sides. Note that the air flow may then be stated by the simple model of Dimitrov and Nagayama [9] and confirmed by previous simulations [20], the accumulation of non volatile solute at the contact line in our set unlike droplets) depends on the total evaporation flux over the men density. Then the results can be easily extended to dip coating To get the evaporation rate v previous experiments were performed with pure water, by cutting plates and the reservoir and measuring the decrease of the water level between the two plates.

Colloidal suspension and polymer solution

Commercial silica suspensions wt=50%). Particle mean size is desired concentrations, the industrial solution i 9. A smaller pH buffer (pH≈2) study concerns the influence of the system pH on . The second study is dedicated to the comparison of silica suspensions and to classical dip-coating. In dip-coating the substrate is a reservoir which contains coating materials, while in our setup the substrate is motionless and is pumped out from the reservoir. Figure 1 shows a picture of the sandwich device made of two paralleled glass plates separated by 1mm immersed into a reservoir which contains the colloidal suspension. The small spacers are located at the top of the glass slides and close to the edges. Due to the 1mm This value of the gap is a compromise. It was chosen to while being able to ensure a good parallelism between the two plates. flow rate of the pumping system controls the contact line velocity V. inside a chamber of volume 60L, where temperature and humidity are regulated by a PID system. In all the experiments presented in this paper, the temperature point is 30%. A fan blows an air flow to the meniscus carefully screwed on the top of the two plates. The air flow velocity is chosen to get a but a negligible deformation due to the pressure increase above the The gap between the two plates being open, the air escapes through the lateral sides.

may then be different from a classical dip-coating configuration. However, y the simple model of Dimitrov and Nagayama [9] and confirmed by previous simulations [20], the accumulation of non volatile solute at the contact line in our set unlike droplets) depends on the total evaporation flux over the meniscus and not on its local be easily extended to dip coating.

v evap corresponding to the experimental conditions previous experiments were performed with pure water, by cutting the connection between the plates and the reservoir and measuring the decrease of the water level between the two plates.

1 Scheme of the experimental setup.

Colloidal suspension and polymer solution suspensions in water have been used (Klebosol from AZ size is 83nm (standard deviation 7nm) measured by , the industrial solution is diluted in pure water. Corresponding pH is about 2) is also used to modify the electrostatic interactions between the on patterning, for aqueous of silica suspensions and substrate is withdrawn from the substrate is motionless and the the experimental setup: a 1mm spacers is vertically The small spacers are located Due to the 1mm-thin gap, there is a a compromise. It was chosen to get a between the two plates. The L, where temperature and humidity are regulated by a PID system. In all the experiments presented in this paper, the temperature set-point is 25°C blows an air flow to the meniscus through a vertical The air flow velocity is chosen to get a deformation due to the pressure increase above the The gap between the two plates being open, the air escapes through the lateral sides.

coating configuration. However, as y the simple model of Dimitrov and Nagayama [9] and confirmed by previous simulations -up or in dip-coating (but iscus and not on its local corresponding to the experimental conditions mentioned above, the connection between the plates and the reservoir and measuring the decrease of the water level between the two plates.

have been used (Klebosol from AZ materials, 50R50 measured by AFM. To prepare the water. Corresponding pH is about to modify the electrostatic interactions between the particles. This buffer is prepared by diluting Nitric Acid (Sigma the suspensions are stirred for one hour at polymer solution, polyacrylamide 9, wt=50%, Mw = 22,4kg/mol, Mw/Mn= 3,5 solutions are stirred for one night and 2.3 Experimental protocol Glass plates are cleaned in "Piranha" solution sulfuric acid with volume ratio of 0.3/0.7 flow. Just after this cleaning procedure glass plates are transferred into the setup and the experiment begins about 15 minutes after. After this delay, the receding contact angl the surface. For all the experiments reported in this article, the filling the tank, a small but fast the receding contact line reaches its static receding value.

Observations:

We record the movement of the experiment by a CCD camera (AVT Marlin F f/2.8D). Due to set-up limitations capillary rise, but the displacement is known very accurately. Indeed a performed, and leads to a precision of are analyzed by Igor software. The dried deposits are characterized by profilemetry optical profiler from FOGALE nanotech

Results and discussion

Stick-slip motion description

By image processing, we can plot the contact line displacement versus time typical example of a periodic motion For brevity it is called stick-slip pinning of the contact line. This "pseudo sessile droplets. In this study, velocity and slip means that it accelerates. length, velocity and duration from the observation of the contact line displacement. Defin straightforward from the diagram drawn in

δh stick = h 2 -h 1 , δh slip = h 3 -h 2 , V
pared by diluting Nitric Acid (Sigma-Aldrich 438073) in pure water ions are stirred for one hour at room temperature before the experiment starts. olymer solution, polyacrylamide (PAAm)/water solutions have been used

, Mw = 22,4kg/mol, Mw/Mn= 3,5, see [21] for thermodynamic properties for one night and filtered just before the experiment.

"Piranha" solution for 20 minutes (mixture of sulfuric acid with volume ratio of 0.3/0.7). After rinsing with pure water, plates are t after this cleaning procedure the glass slides exhibit a total wetting by water. Then the glass plates are transferred into the setup and the experiment begins about 15 minutes after. After this delay, the receding contact angle is about 10°, probably due to uncontrolled contamination of the surface. For all the experiments reported in this article, the contact line is receding. Just after a small but fast pumping out of the reservoir is imposed, so that the co the receding contact line reaches its static receding value.

of the meniscus that is on the top of the capillary rise during all the experiment by a CCD camera (AVT Marlin F201B, resolution 1280 x 1024 pixel, limitations, the camera recording does not provide the absolut ut the displacement is known very accurately. Indeed a to a precision of the order of 2 µm for the contact line displacement are analyzed by Igor software. The dried deposits are characterized by profilemetry nanotech).

description By image processing, we can plot the contact line displacement versus time typical example of a periodic motion, leading to ordered strip patterns parallel to the contact line. slip motion in the following, even if it does not correspond to a . This "pseudo-pinning" has been reported by Moffat and n this study, stick means that the contact line slows down compared to the mean it accelerates. For the two stages (stick and slip), we can deduce length, velocity and duration from the observation of the contact line displacement. Defin straightforward from the diagram drawn in Figure 2: , V stick = δh stick /(t 2 -t 1 ), V slip = δh slip /(t 3 -t 2 ), δh tot ldrich 438073) in pure water. All experiment starts. For the have been used (Sigma-Aldrich 43,494-] for thermodynamic properties). All the polymer of hydrogen peroxide and water, plates are dried by nitrogen the glass slides exhibit a total wetting by water. Then the glass plates are transferred into the setup and the experiment begins about 15 minutes after. After 10°, probably due to uncontrolled contamination of contact line is receding. Just after out of the reservoir is imposed, so that the contact angle at meniscus that is on the top of the capillary rise during all the tion 1280 x 1024 pixel, Nikkor lens 60 mm does not provide the absolute value of the ut the displacement is known very accurately. Indeed a sub pixel analysis is contact line displacement. Images are analyzed by Igor software. The dried deposits are characterized by profilemetry (Microsurf 3D By image processing, we can plot the contact line displacement versus time. Figure 2 shows a patterns parallel to the contact line. even if it does not correspond to a perfect by Moffat and Sefiane [22] for stick means that the contact line slows down compared to the mean tick and slip), we can deduce the length, velocity and duration from the observation of the contact line displacement. Definitions are = δh stick + δh slip

Silica suspensions: Contact line dynamics and deposit morphologies

The first comparison concerns the influence of the solution pH previous results [16][17][18], the process parameters have bee expected. Two experiments have been performed with the same p V≈10µm/s, evaporation rate v the pH value of the solutions. Changing the electrostatic interaction between the particles through the pH affects the deposit morpholog from profilometric profiles of the dried deposits, the cross section obtained with contour, longer wavelength, and quite one side. The relative contact line displacement experiments. Average values of the wavelength ( Figure 3. 3D images obtained from profilemeter measurement of the colloidal depositions different pH values. The plots at the bottom give the thickness profiles, i.e. the cross section of the stripes.

The ratio V slip /V stick is two times larger for pH 2 than for pH 9, which results in the large asymmetry observed on the dried deposit. The slip and stick lengths are smaller for pH 2, and so is the total wavelength. However the mean amplitude is close for the two samples, 1.5 ± µm for pH 9 and 2, respectively. Other experiments were performed by changing the substrate velocity V and the same conclusions were obtained. Summary of the results for these two pH and different substrate velocities is given in Figu be a relevant scaling, so that it was also used in this study. Despite the dispersion of the results, wavelengths obtained at low pH are clearly smaller which confirms the influence of particle interaction forces on the final deposit. slip: typical contact line displacement as a function of times corresponds to an experience performed with receding velocity V m/s, initial volume fraction Φ≈8% and pH≈2. The right figure is a sketch 3.2 Silica suspensions: Contact line dynamics and deposit morphologies he first comparison concerns the influence of the solution pH for silica suspension ], the process parameters have been chosen in a range where stick Two experiments have been performed with the same parameters: receding velocity 10µm/s, evaporation rate v evap ≈0.4µm/s and initial volume fraction Φ≈8%. The only difference is ons. Changing the electrostatic interaction between the particles through the deposit morphologies significantly, as illustrated in Figure 3. As can be seen from profilometric profiles of the dried deposits, the cross section obtained with contour, longer wavelength, and quite symmetric shape. For pH 2, the depos The relative contact line displacement, h(t)-Vt, is given in Figure 4 for the two experiments. Average values of the wavelength (δh tot ), stick and slip velocities are listed in table 1 . 3D images obtained from profilemeter measurement of the colloidal depositions plots at the bottom give the thickness profiles, i.e. the cross section of the is two times larger for pH 2 than for pH 9, which results in the large asymmetry observed on the dried deposit. The slip and stick lengths are smaller for pH 2, and so is the total wavelength. However the mean amplitude is close for the two samples, 1.5 ± µm for pH 9 and 2, respectively. Other experiments were performed by changing the substrate velocity V and the same conclusions were obtained. Summary of the results for these two pH and different substrate velocities is given in Figure 5. Data from Bodiguel et al. [17] suggest v be a relevant scaling, so that it was also used in this study. Despite the dispersion of the results, wavelengths obtained at low pH are clearly smaller which confirms the influence of particle tion forces on the final deposit.

displacement as a function of times for the silica receding velocity V≈4.4µm/s, The right figure is a sketch for silica suspensions. Based on n chosen in a range where stick-slip is arameters: receding velocity 8%. The only difference is ons. Changing the electrostatic interaction between the particles through significantly, as illustrated in Figure 3. As can be seen from profilometric profiles of the dried deposits, the cross section obtained with pH 9 has smoother , the deposit is much sharper at Vt, is given in Figure 4 for the two ), stick and slip velocities are listed in table 1.

. 3D images obtained from profilemeter measurement of the colloidal depositions for plots at the bottom give the thickness profiles, i.e. the cross section of the is two times larger for pH 2 than for pH 9, which results in the large asymmetry observed on the dried deposit. The slip and stick lengths are smaller for pH 2, and so is the total wavelength. However the mean amplitude is close for the two samples, 1.5 ± 0.1 µm and 1.6 ± 0.2 µm for pH 9 and 2, respectively. Other experiments were performed by changing the substrate velocity V and the same conclusions were obtained. Summary of the results for these two pH and re 5. Data from Bodiguel et al. [17] suggest v ev Φ/V to be a relevant scaling, so that it was also used in this study. Despite the dispersion of the results, wavelengths obtained at low pH are clearly smaller which confirms the influence of particle For the colloidal suspension, results are similar to the experiments of the previous section at pH we get a clear stick-slip phenomenon. The amplitude of stripe patterns shape is symmetric (Figure 6 and 7, left). solution: the dry film is almost and 7 (right). The amplitude of these undulations is of magnitude smaller than for the colloidal suspension. whole thickness, so that it is meaningful to define a mean thickness of the dry film. As stat previous study [11], it is expected to be of the order of h µm. This is in agreement with the experimental thickness of 0.5 µm measured from profilometry measurement.

Figure 6: 3D images obtained fr (right) deposits. The artificial thickness of the film.

Varying the initial volume fraction (2% < µm/s), no stick slip was obtained for PAAm solution at room temperature and humidity 30%. It is only by decreasing the polymer volume fraction below 1.5%, or by shifting to higher evaporation temperature that regular stick drawn that stick-slip occurrence depends on system properties and external parameters (substrate velocity, evaporation rate) in a non trivial way.

Understanding this coupling and the mechanism that drives up the periodic movement of the contact line is still an open question and, as far as the authors know, no general theory is currently available to explain it. In the light of these experimental results, one promising try to obtain a periodic regime with an hydrodynamic model, as it has already been done for ultra thin films [23]. For instance the viscosity dependence with solute concentration is known to be very sensitive to the nature of the solute or pH) since it affects particles interaction (see for instance [24]), and to the temperature since it polymer solution repeated with PAAm polymer solution at 25°C. Close parameters V = 5.2 µm/s, v evap =0.38 µm/s, Φ = 3% for the colloidal suspension, Φ = 3.4% for the polymer solution.

For the colloidal suspension, results are similar to the experiments of the previous section at pH slip phenomenon. The amplitude of stripe patterns is about 1.5 shape is symmetric (Figure 6 and 7, left). A very different morphology is obtained for the polymer solution: the dry film is almost flat, with only small and intermittent undulations shown in Figure 6 de of these undulations is of order of 20 nm that is two orders of magnitude smaller than for the colloidal suspension. Here undulations are small compared to the whole thickness, so that it is meaningful to define a mean thickness of the dry film. As stat previous study [11], it is expected to be of the order of h dry ≈LΦv ev /V, with L µm. This is in agreement with the experimental thickness of 0.5 µm measured from profilometry Figure 6: 3D images obtained from profilemeter measurement of the colloidal groove has been made a posteriori on the dry film to estimate the Varying the initial volume fraction (2% < ϕ < 7%) or the substrate velocity (0.5 µm/s < V < 500 µm/s), no stick slip was obtained for PAAm solution at room temperature and humidity 30%. It is only by decreasing the polymer volume fraction below 1.5%, or by shifting to higher evaporation hat regular stick-slip behavior was again observed. From these first results, it can be slip occurrence depends on system properties and external parameters (substrate velocity, evaporation rate) in a non trivial way.

oupling and the mechanism that drives up the periodic movement of the contact line is still an open question and, as far as the authors know, no general theory is currently available to explain it. In the light of these experimental results, one promising try to obtain a periodic regime with an hydrodynamic model, as it has already been done for ultra thin films [23]. For instance the viscosity dependence with solute concentration is known to be very sensitive to the nature of the solute (polymers/colloids) but also to the medium (salt concentration or pH) since it affects particles interaction (see for instance [24]), and to the temperature since it parameters are chosen for = 3% for the colloidal suspension, while V = 5.9

For the colloidal suspension, results are similar to the experiments of the previous section at pH 9: is about 1.5 µm and their A very different morphology is obtained for the polymer undulations shown in Figure 6 order of 20 nm that is two orders of Here undulations are small compared to the whole thickness, so that it is meaningful to define a mean thickness of the dry film. As stated by a /V, with L≈330µm that is h dry ≈0.6 µm. This is in agreement with the experimental thickness of 0.5 µm measured from profilometry om profilemeter measurement of the colloidal (left) and polymer been made a posteriori on the dry film to estimate the 7%) or the substrate velocity (0.5 µm/s < V < 500 µm/s), no stick slip was obtained for PAAm solution at room temperature and humidity 30%. It is only by decreasing the polymer volume fraction below 1.5%, or by shifting to higher evaporation slip behavior was again observed. From these first results, it can be slip occurrence depends on system properties and external parameters (substrate oupling and the mechanism that drives up the periodic movement of the contact line is still an open question and, as far as the authors know, no general theory is currently available to explain it. In the light of these experimental results, one promising track should be to try to obtain a periodic regime with an hydrodynamic model, as it has already been done for ultrathin films [23]. For instance the viscosity dependence with solute concentration is known to be very (polymers/colloids) but also to the medium (salt concentration or pH) since it affects particles interaction (see for instance [24]), and to the temperature since it affects polymer mobility. Next step will then require the detailed characterization of the well as the development of the model to go further in the understanding of patterning.

Conclusion:

For low capillary numbers, we have studie drives convective self-assembly. The results show that with different pH values, we wavelengths and morphologies. only small undulations on the surface understanding the deposit formation.

Figure 1

 1 Figure 1 Scheme of the experimental setup.

Figure 2 :

 2 Figure 2: Stick-slip: typical suspension. Left figure corresponds to an experience evaporation rate v evap ≈0.4µm/s, initial volume fraction of stick-slip motion definition.

Figure 4 :

 4 Figure 4: Relative contact line displacement at two different pH values. Upper one corresponds to pH 9, lower

Figure 7 :

 7 Figure 7: Cross section of the deposits of colloidal and polymer pattern amplitude. Right:

Table 1 .

 1 Stick-slip characteristics for Figure 5. Wavelength as a function of the non dimensional variable v v ev <0.58µm/s, 2.2µm/s< V<12.5 Relative contact line displacement at two different pH values. Upper one , lower one to pH 2.

	δh slip (µm)	δh tot (µm)	V stick (µm/s)
	111	205	7
	98	145	5
	slip characteristics for pH 2 and pH 9 suspensions.

Acknowledgements

We gratefully thank A. Aubertin, L. Auff technical support. Deposit images were performed by FOGALE optical profilometer "Sciences et Ingénierie de la Matière Molle" Marie Curie Initial Training Network: Multiscale complex fluid flows and interfacial phenomena (MULTIFLOW); and Junta de Andalucia (project P07


Antton Goïcoechea

Cécile Brütt

Arthur Le Ber

Flavien Bureau

William Lambert

Claire Prada

Arnaud Derode

Alexandre Aubry

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

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.

Multiple scattering (MS) of waves proves to be a captivating subject manifesting itself across all the spectrum of wave physics [START_REF] Abrahams | Scaling theory of localization: Absence of quantum diffusion in two dimensions[END_REF][START_REF] Albada | Observation of weak localization of light in a random medium[END_REF][START_REF] Lee | Universal conductance fluctuations in metals[END_REF][START_REF] Tourin | Time-dependent coherent backscattering of acoustic waves[END_REF][START_REF] Scheffold | Universal conductance fluctuations of light[END_REF][START_REF] Larose | Weak localization of seismic waves[END_REF][START_REF] Hu | Localization of ultrasound in a three-dimensional elastic network[END_REF][START_REF] Vellekoop | Universal optimal transmission of light through disordered materials[END_REF][START_REF] Popoff | Measuring the transmission matrix in optics: An approach to the study and control of light propagation in disordered media[END_REF][START_REF] Jendrzejewski | Three-dimensional localization of ultracold atoms in an optical disordered potential[END_REF][START_REF] Gérardin | Full transmission and reflection of waves propagating through a maze of disorder[END_REF][START_REF] Hsu | Correlation-enhanced control of wave focusing in disordered media[END_REF][START_REF] Shi | Diffusion in translucent media[END_REF][START_REF] Horodynski | Anti-reflection structure for perfect transmission through complex media[END_REF][START_REF] Yamilov | Anderson local-ization of electromagnetic waves in three dimensions[END_REF]. 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 heterogeneous medium is the scattering mean free path, denoted as ℓ s . This parameter represents the typical distance between successive scattering events. For a time-of-flight t smaller than the corresponding mean free time, τ s = ℓ s /c (with c, the wave velocity), wave propagation 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 >> τ s . The significance of ℓ s extends beyond its role as a fundamental quantity dictating the onset of MS; it also holds paramount importance for characterization purposes. Indeed, ℓ s directly correlates with the micro-architecture of the medium i.e., the scatterer concentration n and their scattering cross-section in discrete media [START_REF] Frisch | Probabilistic methods in applied mathematics[END_REF][START_REF] Turner | Scattering of elastic waves in heterogeneous media with local isotropy[END_REF], or the correlation function in models involving a continuous random disorder [START_REF] Rytov | Principles of Statistical Radiophysics[END_REF][START_REF] Baydoun | Radiative transfer of acoustic waves in continuous complex media: Beyond the helmholtz equation[END_REF]. ℓ s has thus emerged as a promising marker in diverse applications such as biomedical diagnosis using ultrasound [START_REF] Oosterveld | Ultrasound attenuation and texture analysis of diffuse liver disease: Methods and preliminary results[END_REF][START_REF] Aubry | Local measurements of the diffusion constant in multiple scattering media: Application to human trabecular bone imaging[END_REF][START_REF] Aubry | Multiple scattering of ultrasound in weakly inhomogeneous media: Application to human soft tissues[END_REF][START_REF] Sasso | 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[END_REF][START_REF] Mohanty | Characterization of the lung parenchyma using ultrasound multiple scattering[END_REF][START_REF] Guerrero | Anisotropy and spatial heterogeneity in quantitative ultrasound parameters: Relevance to the study of the human cervix[END_REF] or light [START_REF] Durduran | Diffuse optics for tissue monitoring and tomography[END_REF]. Furthermore, it serves as a crucial monitoring parameter in nondestructive testing [START_REF] Hirsekorn | The scattering of ultrasonic waves by polycrystals[END_REF][START_REF] Margetan | Backscattered microstructural noise in ultrasonic toneburst inspections[END_REF][START_REF] Panetta | Ultrasonic attenuation measurements in jet-engine titanium alloys[END_REF][START_REF] Thompson | Scattering of elastic waves in simple and complex polycrystals[END_REF][START_REF] Shahjahan | Comparison between experimental and 2-d numerical studies of multiple scattering in Inconel600® by means of array probes[END_REF][START_REF] Baelde | Effect of microstructural elongation on backscattered field: Intensity measurement and multiple scattering estimation with a linear transducer array[END_REF] and seismology [START_REF] Sato | Seismic Wave Propagation and Scattering in the Heterogeneous Earth : Second Edition[END_REF][START_REF] Chaput | Multiple scattering from icequakes at Erebus volcano, Antarctica: Implications for imaging at glaciated volcanoes[END_REF][START_REF] Mayor | Tomography of crustal seismic attenuation in metropolitan france: Implications for seismicity analysis[END_REF][START_REF] Onodera | Seismic scattering and absorption properties of Mars estimated through coda analysis on a long-period surface wave of S1222a Marsquake[END_REF].

From a theoretical standpoint, the probability density function I(r, t) for the travel time t in a random medium, with r the relative position between the source and the receiver, can be obtained by solving the radiative transfer equation. From an experimental standpoint, 2); the overall intensity obtained via Monte Carlo (blue dots) is compared with its short-time [Eq. ( 4)] and longtime [Eq. ( 3)] analytical expressions. The x-and y-axes are in log-scale. (c) Single scattering rate ρ s (t): The numerical result (black dots) is compared with its short-time expression [Eq. [START_REF] Scheffold | Universal conductance fluctuations of light[END_REF]] and the diffusion model prediction [Eq. [START_REF] Lee | Universal conductance fluctuations in metals[END_REF]]. The x-and y-axes are in linear and log-scale, respectively.

the averaged instantaneous intensity of the wave-field can be measured and serves as an estimate of I(r, t). In the following, the word "intensity" will be used in the sense of a travel-time distribution. Inverting experimental measurements of I(r, t) allows for the derivation of a spatial mapping of ℓ s along with other relevant transport parameters. This methodology aligns with the concept of optical diffuse tomography [START_REF] Durduran | Diffuse optics for tissue monitoring and tomography[END_REF][START_REF] Bal | Transport-based imaging in random media[END_REF][START_REF] Arridge | Optical tomography: forward and inverse problems[END_REF]. 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 [START_REF] Page | Group Velocity in Strongly Scattering Media[END_REF]. The characteristic length scale associated with this attenuation is the extinction length, denoted as ℓ ext , which encapsulates both scattering and absorption losses:

ℓ ext = (ℓ -1 s +ℓ -1 a ) -1
, where ℓ 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 ℓ 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 [START_REF] Paasschens | Solution of the time-dependent Boltzmann equation[END_REF], we first show that the SS rate ρ s (t) scales as exp[-3ct/(4ℓ s )] for t < τ s . This observable can therefore give access to a direct measurement of ℓ s . Experimentally, ρ s (t) is obtained by investigating the reflection matrix of the medium in a focused basis [START_REF] Lambert | Reflection Matrix Approach for Quantitative Imaging of Scattering Media[END_REF][START_REF] Velichko | Quantification of the Effect of Multiple Scattering on Array Imaging Performance[END_REF][START_REF] Lambert | Ultrasound matrix imaging-Part I: The focused reflection matrix, the F-factor and the role of multiple scattering[END_REF] and projecting it onto a characteristic SS space [START_REF] Brütt | Weight of single and recurrent scattering in the reflection matrix of complex media[END_REF]. This method is first validated on a phantom generating an ultrasound speckle characteristic of tissues. We provide local measurements of ℓ s and ℓ 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 [START_REF] Sasso | 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[END_REF][START_REF] Audiere | Influence of heterogeneities on ultrasound attenuation for liver steatosis evaluation (CAP TM ): Relevance of a liver guidance tool[END_REF].

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 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(u out , u in , t) is recorded by each transducer identified by its position u out . By repeating this operation for each element as a source, each reflected wave-field can be stored into a canonical reflection matrix R uu (t) = R(u out , u 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) = Tr R(t)R † (t) /N, (1)
where the symbols Tr and † stands for matrix trace and transpose conjugate, respectively. [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 [START_REF] Paasschens | Solution of the time-dependent Boltzmann equation[END_REF] 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/ℓ ext ) /(4ℓ s ). (2)
Not surprisingly, I S (t) displays an exponential attenuation dictated by ℓ ext . As shown by As to the overall intensity, the medium interface can be taken into account in the diffusive approximation [START_REF] Patterson | Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties[END_REF]. A power law scaling of I(t) as t -3/2 is predicted [40]:

I(t) ∼ t>>τs exp(-ct/ℓ a ) √ πDt z 0 ct (3)
with z 0 = 2ℓ s /3, the extrapolation length [START_REF] Zhu | Internal reflection of diffusive light in random media[END_REF] and D, the diffusion coefficient. In a 2D

geometry and for isotropic scattering, D = cℓ 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 < τ s [40]:

I(t) ∼ t<τs exp (-ct/4ℓ s ) exp (-ct/ℓ a ) /(4ℓ s ). ( 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. [START_REF] Albada | Observation of weak localization of light in a random medium[END_REF]]. This striking property can be highlighted by investigating the SS rate, ρ s (t) = I S (t)/I(t). The ratio of Eqs. ( 2) and (4) shows that ρ s displays an exponential attenuation for t < τ s that only depends on ℓ s :

ρ s (t) ∼ t<τs exp(-3ct/4ℓ s ). (5)
Fig. 1(c) shows the validity of the last expression for t < τ s . This fundamental result proves analytically that a discrimination between SS and MS can lead to a measurement of ℓ s independent of ℓ a [START_REF] Aubry | Multiple scattering of ultrasound in weakly inhomogeneous media: Application to human soft tissues[END_REF]. 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 mm/µs and attenuation α ≃ 0.7 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,[START_REF] Lambert | Reflection Matrix Approach for Quantitative Imaging of Scattering Media[END_REF]. 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.

To estimate the weight of MS, a focused reflection matrix R xx (z) can be synthesized at each depth z [START_REF] Lambert | Reflection Matrix Approach for Quantitative Imaging of Scattering Media[END_REF]. Its coefficients R(x out , x in , z) are the response of the medium at the ballistic time t = 2z/c between virtual sources and detectors located at (x in , z) and (x out , z)

[see Fig.

ÎS (t = 2z/c) = Tr S xx (z)S † xx (z) /N. (6)
The single scattering intensity is estimated along the blue and red lines of Fig. In the red cylinder, a slight increase is observed for the measured ℓ ext but the error bar on this estimation is large. In the blue cylinder, large speckle fluctuations prevent from a reliable fitting process of ÎS .

To go beyond a measurement of ℓ ext , one can compute estimated the SS rate: ρs (t) = ÎS (t)/I(t). Interestingly, the SS rate seems more robust with respect to the medium reflectivity fluctuations [Fig. 3(c)]. ρs exhibits a very different behaviour in the blue and red areas [Fig. 3(d)]. While ℓ 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 [START_REF] Baydoun | Radiative transfer of acoustic waves in continuous complex media: Beyond the helmholtz equation[END_REF]. 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 monitor diseases such as steatosis [START_REF] Sasso | 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[END_REF][START_REF] Suzuki | Dependence of ultrasonic attenuation of liver on pathologic fat and fibrosis: Examination with experimental fatty liver and liver fibrosis models[END_REF][START_REF] Karlas | Individual patient data meta-analysis of controlled attenuation parameter (cap) technology for assessing steatosis[END_REF]. 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 [START_REF] Strauss | Interobserver and intraobserver variability in the sonographic assessment of fatty liver[END_REF]. As a first step towards such biomedical applications, we apply our approach to in-vivo ultrasound data acquired on the liver to reflectivity variations caused by the presence of structures like veins and cannot provide a reliable measurement of ℓ ext . On the contrary, the fitting of ρs with Eq. ( 5) leads to an estimation of ℓ s ≃ 22 ± 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 [START_REF] Lambert | Reflection Matrix Approach for Quantitative Imaging of Scattering Media[END_REF]. 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 ℓ s is therefore limited to a penetration depth of a few ℓ 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 ρ 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: ρ s (t) ∼ exp(-0.55ct/ℓ s ) [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 scattering media, employing a non-invasive reflection setup. Our approach successfully discriminates between SS and MS events, providing an independent measurement of ℓ s and ℓ 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 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 [START_REF] Yoon | Deep optical imaging within complex scattering media[END_REF][START_REF] Bureau | Three-dimensional ultrasound matrix imaging[END_REF][START_REF] Najar | Noninvasive retrieval of the time-gated transmission matrix for optical imaging deep inside a multiple scattering medium[END_REF][START_REF] Giraudat | Unveiling the deep plumbing system of a volcano by a reflection matrix analysis of seismic noise[END_REF].

All the authors are grateful for the funding provided by the European Research Council

SINGLE SCATTERING FILTER

To filter multiple scattering from the focused reflection matrix, an accurate single scattering space should be built [START_REF] Brütt | Weight of single and recurrent scattering in the reflection matrix of complex media[END_REF]. To that aim, we first compute a set of monochromatic reflection matrices P uu (x s , z, ω) = [P (u in , u out , x s , z, ω)] defined in the transducer basis and associated with a single point-like scatterer located at (x s , z), such that

P uu (x s , z, ω) = T ux (z, ω) × Γ xx (x s ) × T T ux (z, ω). ( S5
)
with Γ(x s ), a diagonal matrix whose coefficients γ(x) describes the reflectivity of each scatterrer, such that γ(x) = γ 0 δ(x -x s ), with δ the Dirac distribution.

The next step is the projection of each matrix P uu (x s , z, ω) in the focused basis, as previously done for R uu (z, ω) [Eq. (S2)],

P xx (x s , z, ω) = T † ux (z, ω) × P uu (x s , z, ω) × T * ux (z, ω). (S6)
In practice, we generate a synthetic matrix P xx (x s , z, ω) for a virtual scatterer every half resolution cell. The dimension of this resolution cell is dictated by the characteristic size δx of the focal spot. For an homogeneous speed-of-sound, δx = λ/(2 sin[arctan(z/(2N p))]).

To build a single scattering space, this set of matrices can be orthogonalized by means a Gram-Schmidt process [START_REF] Brütt | Weight of single and recurrent scattering in the reflection matrix of complex media[END_REF]. 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 P(z, ω), at each depth z and frequency ω such that:

P ({x in , x out }, x s , z, ω) = P (x in , x out , x s , z, ω) (S7)
Then, the single scattering basis is determined at each depth and frequency by performing a singular value decomposition of P(z, ω):

P(z, ω) = N k (z,ω) k=1 σ k (z, ω) Fk (z, ω)X k (z, ω) (S8)
where σ k are the singular values ranged in decreasing order. Fk (z, ω) = [ Fk ({x in , x out }, z, ω)] and X k (z, ω) = [X k (x s , z, ω)] 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, ω) 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, ω) = ∆x/δx(z, ω), with ∆x, the lateral extent of the field-of-view.

The set of singular vectors Fk (z, ω) provide an orthonormal basis of focused reflection

matrices F k (z, ω) = [F k (x in , x out , z, ω)] [Fig.2(b) in the accompanying paper] such that F k (x in , x out , z, ω) = Fk ({x in , x out }, z, ω). (S9)
This single scattering basis is used to project each focused reflection matrix R xx (z, ω) onto a characteristic single scattering space:

S xx (z, ω) = N k (z,ω) k=1 Tr F † k (z, ω) × R xx (z, ω) F k (z, ω) (S10)
The set of resulting monochromatic matrices S xx (z, ω) are then recombined to yield the broadband single scattering matrix:

S xx (z) = ω + ω - dωS xx (z, ω) (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 S xx (z) and 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 [START_REF] Lambert | Reflection Matrix Approach for Quantitative Imaging of Scattering Media[END_REF]. Moreover, unlike Ref. [START_REF] Lambert | Ultrasound matrix imaging-Part I: The focused reflection matrix, the F-factor and the role of multiple scattering[END_REF], 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 > τ s since, even in absence of single scattering, a residual multiple scattering contribution will still emerge along the single scattering subspace. However, in the scattering regime considered in this paper (t < τ 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) = Tr R xx (z) × R † xx (z) /N (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 R xx in the frequency domain [Eq. (S4)] and use the expression of R xx (ω) given in Eq. (S2) to express the broadband matrix

R xx : R xx (z) = dωT † ux (z, ω) × R uu (ω) × T * ux (z, ω). (S13)
To simplify the latter expression, the paraxial approximation can be made and the transmission matrix can be decomposed as follows:

T ux (z, ω) = F ux (z, ω) exp j ω c z (S14)
with F ux (z, ω), the Fresnel operator whose coefficients F (u, x, z, ω) are given by:

F (u, x, z, ω)≃ e 3iπ/4 2ω πzc exp j ω 2zc (u -x) 2 . (S15)
Using this decomposition of T ux (z, ω), Equation (S13) can be rewritten as follows

R xx (z) ≃ dωF † ux (z, ω) × R uu (ω) exp -2j ω c z × F * ux (z, ω). ( S16
)
To go beyond, the Fresnel operator can be considered as constant over the frequency bandwith and replaced by its value at the central frequency ω c , such that

R xx (z) ≃ F † ux (z, ω c ) × R uu (t = 2z/c) × F * ux (z, ω c ). ( S17
)
The broadband focused reflection matrix R xx (z) and the time-gated reflection matrix R uu (t = 2z/c) are thus related through a simple change of basis. The unitarity of the Fresnel operator F ux (z, ω c ), implies an equality between their norms:

I(t) = Tr R xx (z) × R † xx (z) /N = Tr R uu (t = 2z/c) × R † uu (t = 2z/c) /N (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) = Tr S xx (z) × S † xx (z) /N = Tr S uu (t = 2z/c) × S † uu (t = 2z/c) /N (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 (r, t, u) in a random medium. Radiance is defined as energy flow, propagating in the direction u, per unit normal area per unit solid angle dΩ per unit time. It follows a transport (Boltzmann) equation :

1 c ∂ ∂t P (r, t, u) + u.∇P (r, t, u) + ℓ -1 ext P (r, t, u) = ℓ -1 s P (r, t) + c -1 S(r, t, u), ( S20
)
with S the source term and P (r, t) the intensity, defined as the angular average of the radiance

P (r, t) = 1 2π dΩP (r, t, u) (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 [START_REF] Paasschens | Solution of the time-dependent Boltzmann equation[END_REF], which were derived for an infinite random medium, to our experimental configuration. The problem was solved in two dimensions. In an elegant way, Passchens 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) = exp (-ct/ℓ ext ) 2πℓ s ct 1 - r 2 c 2 t 2 -1/2 Θ(ct -r) (S22)
with Θ, the Heaviside function: Θ(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) = 1 2 +∞ -∞ dxP S (x, z = 0, t) (S23)
Injecting Eq. ( 2) into the last equation leads to the following expression for I S (t):

I S (t) = exp(-ct/ℓ ext ) 4ℓ s (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 is compared with the heuristic exponential decay exp(-0.55ct/ℓ s ) (red line).

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) = exp (-ct/ℓ a ) 2πℓ s ct 1 - r 2 c 2 t 2 -1/2 exp √ c 2 t 2 -r 2 -ct /ℓ s Θ(ct -r). (S25)
The square root in the exponential can be developed as follows,

exp √ c 2 t 2 -r 2 -ct /ℓ s ∼ exp -r 2 /(2ℓ s ct) (S26) provided that exp -r 4 /(8ℓ s ct 3 ) ∼ 1. (S27)
Given the fact that r < ct, the latter condition is valid for ct < ℓ s . Equation (S25) can thus be simplified into the following expression for t < τ s :

P (r, t) ∼ t<τs exp (-ct/ℓ a ) 2πℓ s ct 1 - r 2 c 2 t 2 -1/2 exp -r 2 /(2ℓ s ct) Θ(ct -r). (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 dividing by two the mean radiance [Eq. (S25)] obtained for an infinite medium. The overall intensity is then given by:

I(t) = 1 2 +∞ -∞ dxP (x, z = 0, t) (S29)
Injecting Eq. (S28) into the last equation leads to the following expression for I(t):

I(t) ∼ t<τs exp(-ct/ℓ a ) exp [-ct/(4ℓ s )] 4ℓ s I 0 ct 4ℓ s (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 < τ s . The following asymptotic expression is thus obtained for I(t)

at short times-of-flight:

I(t) ∼ t<τs exp (-ct/4ℓ s ) exp (-ct/ℓ a ) /(4ℓ s ). (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 < τ s but, as expected, this analytical expression cannot grasp the evolution of the backscattered intensity for t > τ s .

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

ρ s (t) ∼ t<τs exp (-3ct/4ℓ s ) . (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 < τ s . Note that the decay of ρ s (t) is not strictly identical in a 3D geometry but an exponential decay of exp(-0.55ct/ℓ s ) fits pretty well the result of our Monte Carlo simulation for t < τ s [Fig. S1(b)].

DIFFUSE INTENSITY IN THE LONG-TIME LIMIT

In a scattering medium (ℓ s << ℓ a ), the specific intensity P (r, t, u) can be decomposed as the sum of an isotropic and an anisotropic term [START_REF] Wang | Biomedical optics: Principles and imaging[END_REF]:

P (r, t, u) = ϕ(r, t) 2π + 1 π J(r, t).u (S33)
with ϕ(r, t), the fluence rate (or intensity), and J(r, t), the energy flux. In the long time limit (t >> τ s ), the radiative transfer equation then leads to a much simpler equation for ϕ(r, t):

1 c ∂ ∂t ϕ(r, t) -D∆ϕ(r, t) + 1 ℓ a ϕ(r, t) = S(r, t) (S34)
with D, the diffusion coefficient. In a 2D geometry and for isotropic scattering, D = cℓ s /2.

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

F (r, t) = 1 4πDt exp(-ct/ℓ a ) exp -r 2 /4Dt . (S35)
In a semi-infinite medium, the cancellation of ϕ(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ℓ s /3, the extrapolation length [START_REF] Zhu | Internal reflection of diffusive light in random media[END_REF]. The corresponding Green's function F (r, t) is thus given by:

F (r, t) = 1 4πDt exp(-ct/ℓ a ) exp - x 2 + (z -z 0 ) 2 4Dt + exp - x 2 + (z + z 0 ) 2 4Dt . (S36)
The associated energy flux J(r, t) at the medium surface can be deduced from the Fick's law:

J(r, t) = -(D/c)∇F (r, t)| z=0 . (S37)
The diffuse intensity, P d (r, t) = |J(r, t)|, can then be deduced:

P d (r, t) = exp(-ct/ℓ a ) 2πDt z 0 ct exp - x 2 + z 2 0 4Dt (S38)
The total diffuse reflectance,

I d (t) = ∞ -∞ dxP d (x, z, t)
, is finally given by:

I d (t) = exp(-ct/ℓ a ) √ πDt z 0 ct exp - z 2 0 4Dt (S39)
In the long-time limit, 4Dt >> z 2 0 . I d (t) thus simplifies into If absorption is negligible, a power law decay as t -3/2 is thus predicted for the scattered intensity in the diffusive limit.

I d (t) = exp(-ct/ℓ a ) √ πDt z 0 ct . ( S40
Using Eqs. ( 2) and (S40), an expression can be derived for the single scattering rate, ρ s (t) ∼ I s (t)/I d (t), in the long-time limit:

ρ s (t) ∼ t>τs 3 8 π 2 (ct/ℓ s ) 3/2 exp(-ct/ℓ s ). (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 ℓ s independent from ℓ a . Nevertheless, the scaling with ℓ 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 [START_REF] Brütt | Weight of single and recurrent scattering in the reflection matrix of complex media[END_REF].

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 σ is of 2.5×10 -2 mm and their concentration is n = 0.2 mm -2 . Under the independent scattering approximation (ISA) [START_REF] Sheng | Introduction to Wave Scattering, Localization, and Mesoscopic Phenomena[END_REF], the scattering mean free path ℓ s ∼ (nσ) -1 , that is to say ℓ s ∼ 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 ρs (t) is displayed in Fig. S2(b). A quantitative agreement is found between ρs (t) and its theoretical scaling in exp(-3ct/4ℓ s ) in the short time limit. The same observation holds for a more concentrated scattering medium for which n = 1 mm -2 and ℓ s ∼ 39 mm [Fig. S2(c)]. Again, our theoretical prediction [Eq. ( 5)] perfectly fits the single scattering rate ρs (t) for t < τ 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 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 ρ s (t) is needed to retrieve the depth evolution of the scattering mean free path ℓ s . * Present address: Université de Rennes, CNRS, IETR (Institut d' Électronique et des Technologies du numéRique), UMR-6164, Rennes, France † alexandre.aubry@espci.fr

FIG. 1 .

 1 FIG. 1.

Fig. 1 (

 1 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

FIG. 2 .

 2 FIG. 2.

  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[START_REF] Brütt | Weight of single and recurrent scattering in the reflection matrix of complex media[END_REF]. It consists in a projection of 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 S xx (z) [Fig.2(d)].The norm of 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:

  3(a) by placing the probe along each cylinder axis. As expected, ÎS displays a similar attenuation decay in each region of the phantom whatever the local reflectivity. At shallow depth, the measured values for ℓ ext [see Tab. I] are in good agreement with the manufacturer value at the central frequency of 7.5 MHz (ℓ ext ∼ 19 mm).

FIG. 3 .

 3 FIG. 3. (a) Experimental configuration and ultrasound image of the tissue-mimicking phantom. (b,c) ÎS (b) and ρ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 ℓ s and ℓ ext reported in Tab. I.

FIG. 4 .

 4 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) ρs (blue dots) and ÎS (green dots) vs. depth z = ct/2. The y-axes are in log-scale. The fit of ρ s (t) with Eq. (5) (blue line) leads to an estimation of ℓ s .

(

  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.at r in = (x in , z) and r out = (x out , z). The last step consists in summing the focused reflection matrix over the frequency bandwidth [ω -; ω + ]:R xx (z) = ω + ω - dωR xx (z, ω)e jωt . (S4)This operation is equivalent to an inverse Fourier transform at lapse time τ = 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.

  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. FigureS1(a) also shows the timeof-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)].

  FIG. S1. Monte Carlo simulation of energy transport in a 3D semi-infinite random medium [ℓ s = 0.1 m, ℓ a = 100 m]. (a) (b) Reflected intensity vs. t/τ 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/ℓ s ) (red line) and the 3D diffusive result [47] (black line). (c) Single scattering rate ρ s (t): The numerical result (black dots)

  ) 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 for t > τ s both in 2D [Fig. 1(b)] and 3D [Fig. S1(a)].

  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

2 mm - 2

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

TABLE I .

 I Values of ℓ ext and ℓ s measured along the blue and red lines of the phantom displayed in

	Fig. 3.
	Depth range	ℓ ext	ℓ s
	25-36 mm (blue)	17 ± 3 mm	247 ± 90 mm
	36-43 mm (blue)	-	36 ± 7 mm
	25-36 mm (red)	15 ± 2 mm	90 ± 11 mm
	36-43 mm (red)	23 ± 8 mm	90 ± 11 mm

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 [START_REF] Lambert | Reflection Matrix Approach for Quantitative Imaging of Scattering Media[END_REF]. To do so, the first step is to perform a temporal Fourier transform of R uu (t):

with ω, 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:

where the symbols † 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, ω) correspond to the z-derivative of the Green's functions G(u, x, z, ω) between the transducer and focal planes. G ur depends on the spatial distribution of the wave speed c(r) inside the system. For an homogeneous speed of sound c 0 and in a 2D geometry, these Green's functions are given by:

where H 0 is the Hankel function of the first kind and k = ω/c 0 , the wave number.

The coefficients R(x out , x in , z, ω) of each focused reflection matrix R xx (z, ω) are the response of the medium at angular frequency ω between virtual sources and detectors located