# Non-invasive Retrieval of the Time-Gated Transmission Matrix for Optical Imaging Deep Inside a Multiple Scattering Medium

Ulysse Najar, Victor Barolle, Paul Balondrade, Mathias Fink, A Claude Boccara, Alexandre Aubry

► To cite this version:

Ulysse Najar, Victor Barolle, Paul Balondrade, Mathias Fink, A Claude Boccara, et al. Non-invasive Retrieval of the Time-Gated Transmission Matrix for Optical Imaging Deep Inside a Multiple Scattering Medium. 2023. hal-03981863v2

HAL Id: hal-03981863

https://hal.science/hal-03981863v2

Preprint submitted on 17 Oct 2023

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
Harnessing Forward Multiple Scattering
for Optical Imaging Deep Inside an Opaque Medium

Ulysse Najar, Victor Barolle, Paul Balondrade,

Mathias Fink, Claude Boccara, and Alexandre Aubry*

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

*Corresponding author (e-mail: alexandre.aubry@espci.fr)

(Dated: August 13, 2024)
Non-invasive Retrieval of the Time-gated Transmission Matrix for
Optical Imaging Deep Inside a Multiple Scattering Medium

Ulysse Najar, Victor Barolle, Paul Balondrade,

Mathias Fink, A. Claude Boccara, and Alexandre Aubry

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

(Dated: October 5, 2023)
Abstract

As light travels through a disordered medium such as biological tissues, it undergoes multiple scattering events. This phenomenon is detrimental to in-depth optical microscopy, as it causes a drastic degradation of contrast, resolution and brightness of the resulting image beyond a few scattering mean free paths. However, the information about the inner reflectivity of the sample is not lost; only scrambled. To recover this information, a matrix approach of optical imaging can be fruitful. Here, we report on a de-scanned measurement of a high-dimension reflection matrix $R$ via low coherence interferometry. Then, we show how a set of independent focusing laws can be extracted for each medium voxel through an iterative multi-scale analysis of wave distortions contained in $R$. It enables an optimal and local compensation of forward multiple scattering paths and provides a three-dimensional confocal image of the sample as the latter one had become digitally transparent. The proof-of-concept experiment is performed on a human opaque cornea and an extension of the penetration depth by a factor five is demonstrated compared to the state-of-the-art.

Introduction

Multiple scattering of waves concerns many domains of physics, ranging from optics or acoustics to solid-state physics, seismology, medical imaging, or telecommunications. In an inhomogeneous medium where the refractive index $n$ depends on the spatial coordinates $r$, several physical parameters are relevant to characterize wave propagation: (i) the scattering mean free path $\ell_s$, which is the average distance between two successive scattering events; (ii) the transport mean free path $\ell_t$, which is the distance after which the wave has lost the memory of its initial direction. For a penetration depth $z$ smaller than $\ell_s$, ballistic light is predominant and standard focusing methods can be employed; for $z > \ell_s$, multiple scattering events result in a gradual randomization of the propagation direction before reaching the diffusive regime for $z > \ell_t$. Although it gives rise to fascinating interference phenomena such as perfect transmission$^{1,2}$ or Anderson localization$^{3,4}$, multiple scattering still represents a major obstacle to deep imaging and focusing of light inside complex media$^{5,6}$.

To cope with the fundamental issue of multiple scattering, several approaches have been proposed to enhance the single scattering contribution drowned into a predominant diffuse
Abstract

As light travels through a disordered medium such as biological tissues, it undergoes multiple scattering events. This phenomenon is detrimental to in-depth optical microscopy, as it causes a drastic degradation of contrast, resolution and brightness of the resulting image beyond a few scattering mean free paths. However, the information about the inner reflectivity of the sample is not lost; only scrambled. To recover this information, a matrix approach of optical imaging can be fruitful. Here we report on a de-scanned measurement of a high-dimension reflection matrix $R$ via low coherence interferometry. Then, we show how the time-gated transmission matrix linking each camera sensor and each medium voxel can be extracted through an iterative multi-scale analysis of wave distortions contained in $R$. This transmission matrix is the Holy Grail for volumetric imaging since it enables an optimal compensation of forward multiple scattering paths and provides a three-dimensional confocal image of the sample as the latter one had become digitally transparent. The proof-of-concept experiment is performed on a human opaque cornea and an extension of the penetration depth by a factor five is demonstrated compared to the state-of-the-art.

Multiple scattering of waves concerns many domains of physics, ranging from optics or acoustics to solid-state physics, seismology, medical imaging, or telecommunications. In an inhomogeneous medium where the refractive index $n$ depends on the spatial coordinates $r$, several physical parameters are relevant to characterize wave propagation: (i) the scattering mean free path $\ell_s$, which is the average distance between two successive scattering events; (ii) the transport mean free path $\ell_t$, which is the distance after which the wave has lost the memory of its initial direction. For a penetration depth $z$ smaller than $\ell_s$, ballistic light is predominant and standard focusing methods can be employed; for $z > \ell_s$, multiple scattering events result in a gradual randomization of the propagation direction before reaching the diffusive regime for $z > \ell_t$. Although it gives rise to fascinating interference phenomena such as perfect transmission$^{1,2}$ or Anderson localization$^{3,4}$, multiple scattering still represents a major obstacle to deep imaging and focusing of light inside complex media$^{5,6}$.

To cope with the fundamental issue of multiple scattering, several approaches have been proposed to enhance the single scattering contribution drowned into a predominant diffuse
background$^{5,7,8}$. One solution is to perform a confocal discrimination and coherent time gating of singly-scattered photons by means of interferometry. This is the principle of optical coherence tomography$^{9}$, equivalent to ultrasound imaging for light. Nevertheless, a lot of photons associated with distorted trajectories are rejected by the confocal filter while they still contain a coherent information on the medium reflectivity. Originally developed in astronomy$^{10}$, adaptive optics (AO) has been transposed to optical microscopy in order to address this issue$^{11}$. Nevertheless, it only compensates for low-order aberrations induced by long-scale fluctuations of the optical index and does not address high-order aberrations generated by forward multiple scattering events. To circumvent the latter problem, one has to go beyond a confocal scheme and investigate the cross-talk between the pixels of the image. This is the principle of matrix imaging in which the relation between input and output wave-fields is investigated under a matrix formalism.

While a subsequent amount of work has considered the transmission matrix $\mathbf{T}$ for optimizing wave control and focusing through complex media$^{12-17}$, this configuration is not the most relevant for imaging purposes since only one side of the medium is accessible for most in-vivo applications. Moreover, in all the aforementioned works, the scattering medium is usually considered as a black box, while imaging requires to open it. To that aim, a reflection matrix approach of wave imaging (RMI) has been developed for the last few years$^{18-21}$. The objective is to determine, from the reflection matrix $\mathbf{R}$, the $\mathbf{T}$-matrix between sensors outside the medium and voxels mapping the sample$^{22}$. Proof-of-concept studies have reported penetration depths ranging from 7 $\ell_s^{23}$ to 10 $\ell_s^{19}$ but the object to image was a resolution target whose strong reflectivity artificially extends the penetration depth by several $\ell_s$ compared with direct tissue imaging$^8$. Follow-up studies also considered the imaging of highly reflecting structures (e.g. myelin fibers) through an aberrating layer (e.g mouse skull)$^{20}$, in a wavelength range that limits scattering and aberration from tissues$^{24}$. On the contrary, here, we want to address the extremely challenging case of three-dimensional imaging of biological tissues themselves (cells, collagen, extracellular matrix etc.) at large penetration depth ($z \sim 10\ell_s$), regime in which aberration and scattering effects are spatially-distributed over multiple length-scales.

Inspired by previous works$^{25,26}$, full-field optical coherence tomography (FFOCT)$^{27,28}$ will be used here to record the $\mathbf{R}$-matrix. In FFOCT, the incident wave-field is temporally- and spatially-incoherent. It enables, by means of low coherence interferometry, a parallel
background$^{5,7,8}$. One solution is to perform a confocal discrimination and coherent time gating of singly-scattered photons by means of interferometry. This is the principle of optical coherence tomography$^{9}$, equivalent to ultrasound imaging for light. Nevertheless, a lot of photons associated with distorted trajectories are rejected by the confocal filter while they still contain a coherent information on the medium reflectivity. Originally developed in astronomy$^{10}$, adaptive optics (AO) has been transposed to optical microscopy in order to address this issue$^{11}$. Nevertheless, it only compensates for low-order aberrations induced by long-scale fluctuations of the optical index and does not address high-order aberrations generated by forward multiple scattering events. To circumvent the latter problem, one has to go beyond a confocal scheme and investigate the cross-talk between the pixels of the image. This is the principle of matrix imaging in which the relation between input and output wave-fields is investigated under a matrix formalism.

While a subsequent amount of work has considered the transmission matrix $\mathbf{T}$ for optimizing wave control and focusing through complex media$^{12-17}$, this configuration is not the most relevant for imaging purposes since only one side of the medium is accessible for most in-vivo applications. Moreover, in all the aforementioned works, the scattering medium is usually considered as a black box, while imaging requires to open it. To that aim, a reflection matrix approach of wave imaging (RMI) has been developed for the last few years$^{18-21}$. The objective is to determine, from the reflection matrix $\mathbf{R}$, the $\mathbf{T}$-matrix between sensors outside the medium and voxels mapping the sample$^{22}$. Proof-of-concept studies have reported penetration depths ranging from 7 $\ell_s^{23}$ to 10 $\ell_s^{19}$ but the object to image was a resolution target whose strong reflectivity artificially extends the penetration depth by several $\ell_s$ compared with direct tissue imaging$^8$. Follow-up studies also considered the imaging of highly reflecting structures (e.g. myelin fibers) through an aberrating layer (e.g mouse skull)$^{20}$, in a wavelength range that limits scattering and aberration from tissues$^{24}$. On the contrary, here, we want to address the extremely challenging case of three-dimensional imaging of biological tissues themselves (cells, collagen, extracellular matrix etc.) at large penetration depth ($z \sim 10\ell_s$), regime in which aberration and scattering effects are spatially-distributed over multiple length-scales.

Inspired by previous works$^{25,26}$, full-field optical coherence tomography (FFOCT)$^{27,28}$ will be used here to record the $\mathbf{R}$-matrix. In FFOCT, the incident wave-field is temporally- and spatially-incoherent. It enables, by means of low coherence interferometry, a parallel
acquisition of a time-gated confocal image<sup>29</sup> at a much better signal-to-noise ratio than a
traditional point scanning scheme for equal measurement time and power<sup>30</sup>. By splitting the
incident wave-field into two laterally-shifted components, a de-scanned measurement of **R**
can be performed without a tedious raster scanning of the field-of-view<sup>20</sup>.

Another advantage of the de-scanned basis is the direct access to the distortion matrix **D** through a Fourier transform. This matrix basically connects any focusing point with the distorted part of the associated reflected wavefront<sup>19,21</sup>. A multi-scale analysis of **D** is here proposed to estimate the forward multiple scattering component of the **T**-matrix at an unprecedented spatial resolution (~ 6 µm). Once the latter matrix is known, one can actually unscramble, in post-processing, all wave distortions and multiple scattering events undergone by the incident and reflected waves for each voxel. A three-dimensional confocal image of the medium can then be retrieved as if the medium had been made digitally transparent.

The experimental proof-of-concept presented in this paper is performed on a human ex-vivo cornea that we chose deliberately to be extremely opaque. Its overall thickness is of $10\ell_s$. FFOCT shows an imaging depth limit of $2\ell_s$ due to aberration and scattering. Strikingly, RMI enables to recover a full 3D image of the cornea at a resolution close to $\lambda/4$ (~ 230 nm) and a penetration depth enhanced by, at least, a factor five.

## Results

### Measuring a De-scan Reflection Matrix

Our approach is based on a de-scanned measurement of the time-gated reflection matrix **R** from the scattering sample. Inspired by time-domain FFOCT<sup>27,28</sup>, the corresponding set up is displayed in Fig. S1a. It consists in a Michelson interferometer with microscope objectives in both arms (Fig. S1a). In the first arm, a reference mirror is placed in the focal plane of a microscope objective (MO). The second arm contains the scattering sample to be imaged. Because of the broad spectrum of the incident light, interferences occur between the two arms provided that the optical path difference through the interferometer is close to zero. The length of the reference arm determines the slice of the sample (coherence volume) to be imaged and is adjusted in order to match with the focal plane of the MO in the sample arm. The backscattered light from each voxel of the coherence volume can only interfere with the
acquisition of a time-gated confocal image<sup>29</sup> at a much better signal-to-noise ratio than a
traditional point scanning scheme for equal measurement time and power<sup>30</sup>. By splitting the
incident wave-field into two laterally-shifted components, a de-scanned measurement of **R**
can be performed without a tedious raster scanning of the field-of-view<sup>20</sup>.

Another advantage of the de-scanned basis is the direct access to the distortion matrix **D** through a Fourier transform. This matrix basically connects any focusing point with the distorted part of the associated reflected wavefront<sup>19,21</sup>. A multi-scale analysis of **D** is here proposed to estimate the forward multiple scattering component of the **T**-matrix at an unprecedented spatial resolution (~ 6 µm). Once the latter matrix is known, one can actually unscramble, in post-processing, all wave distortions and multiple scattering events undergone by the incident and reflected waves for each voxel. A three-dimensional confocal image of the medium can then be retrieved as if the medium had been made digitally transparent.

The experimental proof-of-concept presented in this paper is performed on a human *ex-vivo* cornea that we chose deliberately to be extremely opaque. Its overall thickness is of 10*l*<sub>s</sub>. FFOCT shows an imaging depth limit of 2*l*<sub>s</sub> due to aberration and scattering. Strikingly, RMI enables to recover a full 3D image of the cornea at a resolution close to λ/4 (~ 230 nm) and a penetration depth enhanced by, at least, a factor five.

Measuring a De-scan Reflection Matrix

Our approach is based on a de-scanned measurement of the time-gated reflection matrix **R** from the scattering sample. Inspired by time-domain FFOCT<sup>27,28</sup>, the corresponding set up is displayed in Fig. 1a. It consists in a Michelson interferometer with microscope objectives in both arms (Fig. 1a). In the first arm, a reference mirror is placed in the focal plane of a microscope objective (MO). The second arm contains the scattering sample to be imaged. Because of the broad spectrum of the incident light, interferences occur between the two arms provided that the optical path difference through the interferometer is close to zero. The length of the reference arm determines the slice of the sample (coherence volume) to be imaged and is adjusted in order to match with the focal plane of the MO in the sample arm. The backscattered light from each voxel of the coherence volume can only interfere with the light coming from the conjugated point of a reference mirror. The spatial incoherence of the light source actually acts as a physical confocal pinhole (Fig. 1c). All these interference
<figure><img src="image_1.png" /><figcaption>FIG. 1. **De-scanned measurement of the Reflection Matrix.** a Experimental setup (L: lenses, MO: microscope objectives, M: reference mirror, BS: beam splitter). Light from an incoherent source is split into two replica laterally shifted with respect to each other by a relative position $\Delta\rho$ (see Supplementary Section S1). By a game of polarization, each replica illuminates one arm of a Linnik interferometer. The sample beam (in red) illuminates the scattering sample through a microscope objective (NA = 1.0). The reference beam (in blue) is focused on a dielectric mirror through an identical microscope objective. Both reflected beams interfere on a CMOS camera whose surface is conjugated with focal planes of the MO. The amplitude and phase of the interference term are retrieved by phase-shifting interferometry. b Each pixel of the camera, depicted by its position $\rho_{\text{out}}$, measures the reflection coefficient $R(\rho_{\text{in}}, \rho_{\text{out}}, z)$ between de-scanned focusing points, $\rho_{\text{out}}$ and $\rho_{\text{in}} = \rho_{\text{out}} + \Delta\rho_{\text{in}}$, at depth $z$ within the sample. c For $\Delta\rho_{\text{in}} = 0$, the experimental set up is equivalent to a FFOCT apparatus and the interferogram directly provides a time-gated confocal image of the sample. d The set of interferograms are stored in the de-scanned reflection matrix $\mathbf{R}_{\text{in}}(z) = [R_{\text{in}}(\Delta\rho_{\text{in}}, \rho_{\text{out}}, z)]$ displayed in panel f. e Each column of this matrix yields a reflection point-spread function (RPSF) associated with the focusing quality at point $\rho_{\text{out}}$ (scale bar: 2 µm). g The Fourier transform (FT) of each de-scanned wave-field provides the input distortion matrix $\mathbf{D}_{\text{in}}(z) = [D_{\text{in}}(\mathbf{u}_{\text{in}}, \rho_{\text{out}}, z)]$. h Each column of this matrix displays the distorted wave-front associated with each point $\rho_{\text{out}}$ in the field-of-view. The optical data shown in panels d-h correspond to the acquisition performed at depth $z = 150$ µm.</figcaption></figure>

light coming from the conjugated point of a reference mirror. The spatial incoherence of the light source actually acts as a physical confocal pinhole (Fig. S1c). All these interference
<figure><img src="image_2.png" /><figcaption>FIG. 1. **De-scanned measurement of the Reflection Matrix.** a Experimental setup (L: lenses, MO: microscope objectives, and M: reference mirror). Light from an incoherent source is split into two replica laterally shifted with respect to each other by a relative position $\Delta\rho$ (see Supplementary Section S1). By a game of polarization, each replica illuminates one arm of a Linnik interferometer. The sample beam (in red) illuminates the scattering sample through a microscope objective (NA = 1.0). The reference beam (in blue) is focused on a dielectric mirror through an identical microscope objective. Both reflected beams interfere on a CMOS camera whose surface is conjugated with focal planes of the MO. The amplitude and phase of the interference term are retrieved by phase-shifting interferometry. b Each pixel of the camera, depicted by its position $\rho_{\text{out}}$, measures the reflection coefficient $R(\rho_{\text{in}}, \rho_{\text{out}}, z)$ between de-scanned focusing points, $\rho_{\text{out}}$ and $\rho_{\text{in}} = \rho_{\text{out}} + \Delta\rho_{\text{in}}$, at depth $z$ within the sample. c For $\Delta\rho_{\text{in}} = 0$, the experimental set up is equivalent to a FFOCT apparatus and the interferogram directly provides a time-gated confocal image of the sample. d The set of interferograms are stored in the de-scanned reflection matrix $\mathbf{R}_{\text{in}}(z) = [R_{\text{in}}(\Delta\rho_{\text{in}}, \rho_{\text{out}}, z)]$ displayed in panel e. f Each column of this matrix yields a reflection point-spread function (RPSF) associated with the focusing quality at point $\rho_{\text{out}}$ (scale bar: 2 $\mu$m). g The Fourier transform of each de-scanned wave-field provides the input distortion matrix $\mathbf{D}_{\text{in}}(z) = [D_{\text{in}}(\mathbf{u}_{\text{in}}, \rho_{\text{out}}, z)]$. h Each column of this matrix displays the distorted wave-front associated with each point $\rho_{\text{out}}$ in the field-of-view. The optical data shown in panels d-h correspond to the acquisition performed at depth $z = 150 \ \mu$m.</figcaption></figure>

signals are recorded in parallel by the pixels of the camera in the imaging plane. Their
amplitude and phase are retrieved by phase-stepping interferometry $^{28}$. The FFOCT signal
signals are recorded in parallel by the pixels of the camera in the imaging plane. Their
amplitude and phase are retrieved by phase-stepping interferometry <sup>28</sup>. The FFOCT signal
is thus equivalent to a time-gated confocal image of the sample<sup>29</sup>. Figures 2b and c show
en-face and axial FFOCT images of the opaque cornea at different depths. A dramatic
loss in contrast is found beyond the epithelium (z > 70 µm, see Fig. 2g). It highlights the
detrimental effect of multiple scattering for deep optical imaging.

<figure><img src="image_2.png" /><figcaption>FIG. 2. **Volumetric matrix imaging of an opaque cornea.** a. Schematic of the imaging planes in the cornea. b. En-face confocal images before ($\mathbf{b}_1$, $\mathbf{b}_3$) and after ($\mathbf{b}_2$, $\mathbf{b}_4$) the matrix imaging process for $z = 50$ µm and 250 µm, respectively (scale bar: 50 µm). c. Longitudinal ($x,z$) section of the initial confocal image. d. Original RSPFs from $z = 50$ to 250 µm (scale bar: 2 µm). e. Corresponding RSPFs after the matrix imaging process. f. Longitudinal ($x,z$) section of the volumetric image at the end of the matrix imaging process. g. Schematic of a healthy human cornea. Each image is normalized at each depth by its averaged intensity.</figcaption></figure>

To overcome the multiple scattering phenomenon, one should go beyond a simple confocal
image and record the cross-talk between the camera pixels. Experimentally, it consists in
measuring the reflection matrix **R** associated with the sample (Fig. 3a). Interestingly, this
can be done by slightly modifying the illumination scheme of the FFOCT device, as displayed
in Fig. S1a. The incident wave-fields are still identical in each arm but are laterally shifted
is thus equivalent to a time-gated confocal image of the sample$^{29}$. Figures 2b and c show en-face and axial FFOCT images of the opaque cornea at different depths. A dramatic loss in contrast is found beyond the epithelium ($z > 70$ µm, see Fig. 2g). It highlights the detrimental effect of multiple scattering for deep optical imaging.

<figure><img src="image_3.png" /><figcaption>FIG. 2. **Volumetric matrix imaging of an opaque cornea.** a. Schematic of the imaging planes in the cornea. b. En-face confocal images before ($b_1, b_3$) and after ($b_2, b_4$) the matrix imaging process for $z = 50$ µm and 250 µm, respectively (scale bar: 50 µm). c. Longitudinal ($x,z$) section of the initial confocal image. d. Original RSPFs from $z = 50$ to 250 µm (scale bar: 1 µm). e. Corresponding RSPFs after the matrix imaging process. f. Longitudinal ($x,z$) section of the volumetric image at the end of the matrix imaging process. g. Schematic of a healthy human cornea. Each image is normalized at each depth by its averaged intensity.</figcaption></figure>

To overcome the multiple scattering phenomenon, one should go beyond a simple confocal image and record the cross-talk between the camera pixels. Experimentally, it consists in measuring the reflection matrix $R$ associated with the sample (Fig. 3a). Interestingly, this can be done by slightly modifying the illumination scheme of the FFOCT device, as displayed in Fig. 1a. The incident wave-fields are still identical in each arm but are laterally shifted with respect to each other by a transverse position $\Delta\rho_{in}$. Their spatial incoherence now acts as a de-scanned pinhole that gives access to the cross-talk between distinct focusing points
<figure><img src="image_3.png" /><figcaption>FIG. 3. Different Stages of Matrix Imaging. a. The focused R-matrix contains the set of impulse responses $R(\rho_{in}, \rho_{out}, z)$ between an array of point sources $\rho_{in}$ and detectors $\rho_{out}$ lying in planes conjugated with the focal plane of the microscope objective (BS: beam splitter). b. The interferometric set up displayed in Fig. S1 allows a de-scanned measurement of R by scanning the relative position $\Delta\rho_{in} = \rho_{in} - \rho_{out}$. c. Each column of the recorded matrix $R_{in}$ (Eq. 1) corresponds to the RPSF measured by each camera pixel. d. A spatial Fourier transform (FT) over $\Delta\rho_{in}$ provides the distortion matrix $D_{in}$ (Eq. 5) linking each camera pixel with wave-front distortions seen from the input pupil plane ($u_{in}$). e. The correlation matrix $C_{in}$ between those wave-fronts mimics the time-reversal operator associated with a virtual guide star that results from a coherent average of all the de-scanned focal spots displayed in b (Supplementary Section S2). f. IPR is then applied (Methods). The resulting wave-front compensates for aberrations and scattering inside the medium to produce a sharper guide star. It provides an estimation of one column of $T_{in}$ corresponding to the common mid-point $r_p$ of the input focal spots considered in panel a.</figcaption></figure>

with respect to each other by a transverse position $\Delta\rho_{in}$. Their spatial incoherence now acts as a de-scanned pinhole that gives access to the cross-talk between distinct focusing points (Fig. S1b). The interferogram recorded by the camera (Fig. S1d) directly provides one line of the reflection matrix $R_{in}$ de-scanned at input (Figs. 3b and c), such that

$$R_{in}(\Delta\rho_{in}, \rho_{out}, z) = R(\rho_{out} + \Delta\rho_{in}, \rho_{out}, z), \quad (1)$$

with $R = [R(\rho_{in}, \rho_{out}, z)]$, the reflection matrix expressed in the canonical basis. Its coeffi-
<figure><img src="image_4.png" /><figcaption>FIG. 3. Different Stages of Matrix Imaging. a. The focused R-matrix contains the set of impulse responses $R(\rho_{in}, \rho_{out}, z)$ between an array of point sources $\rho_{in}$ and detectors $\rho_{out}$ lying in planes conjugated with the focal plane of the microscope objective. b. The interferometric set up displayed in Fig. 1 allows a de-scanned measurement of R by scanning the relative position $\Delta\rho_{in} = \rho_{in} - \rho_{out}$. c. Each column of the recorded matrix $R_{in}$ (Eq. 1) corresponds to the RPSF measured by each camera pixel. d. A spatial Fourier transform over $\Delta\rho_{in}$ provides the distortion matrix $D_{in}$ (Eq. 5) linking each camera pixel with wave-front distortions seen from the input pupil plane ($u_{in}$). e. The correlation matrix $C_{in}$ between those wave-fronts mimics the time-reversal operator associated with a virtual guide star that results from a coherent average of all the de-scanned focal spots displayed in b (Supplementary Section S7). f. IPR is then applied (Methods). The resulting wave-front compensates for aberrations and scattering inside the medium to produce a sharper guide star. It provides an estimation of one column of $T_{in}$ corresponding to the common mid-point $r_p$ of the input focal spots considered in panel a.</figcaption></figure>

(Fig. 1b). The interferogram recorded by the camera (Fig. 1d) directly provides one line of the reflection matrix $R_{in}$ de-scanned at input (Figs. 3b and c), such that

$$R_{in}(\Delta\rho_{in}, \rho_{out}, z) = R(\rho_{out} + \Delta\rho_{in}, \rho_{out}, z), \quad (1)$$

with $R = [R(\rho_{in}, \rho_{out}, z)]$, the reflection matrix expressed in the canonical basis. Its coefficients $R(\rho_{in}, \rho_{out}, z)$ correspond to the response of the medium at depth $z$ between points $\rho_{in}$ and $\rho_{out}$ in the source and camera planes (Fig. 3a). Scanning the relative position $\Delta\rho_{in}$ is equivalent to recording the canonical R-matrix diagonal-by-diagonal (see Figs. 3a and c).
cients $R(\rho_{in}, \rho_{out}, z)$ correspond to the response of the medium at depth $z$ between points $\rho_{in}$ and $\rho_{out}$ in the source and camera planes (Fig. 3a). Scanning the relative position $\Delta\rho_{in}$ is equivalent to recording the canonical $\mathbf{R}$-matrix diagonal-by-diagonal (see Figs. 3a and c). However, while a raster scan (column-by-column acquisition) of $\mathbf{R}$ requires to illuminate the sample over a field-of-view $\Omega$ with $N = (\Omega/\delta_0)^2$ input wave-fronts $^{20,31,32}$, the de-scanned basis allows a much smaller number of field measurements.

This sparsity can be understood by expressing theoretically the de-scan matrix $\mathbf{R}_{in}$ (Supplementary Section S2):

$$R_{in}(\Delta\rho_{in}, \rho_{out}, z) = \int_{\Omega} d\rho_{s} H_{in}(\rho_{s} + \Delta\rho_{in}, \rho_{in}, z)\gamma(\rho_{s} + \rho_{out}, z)H_{out}(\rho_{s}, \rho_{out}, z) \quad (2)$$

where $\gamma$ is the sample reflectivity. $H_{in}(\rho_{s}, \rho_{in}, z)$ and $H_{out}(\rho_{s}, \rho_{out}, z)$ are the local input and output point spread functions (PSFs) at points $(\rho_{in}, z)$ and $(\rho_{out}, z)$, respectively. This last equation confirms that the central line of $\mathbf{R}_{in}$ $(\Delta\rho_{in} = 0)$, i.e. the FFOCT image, results from a convolution between the sample reflectivity $\gamma$ and the local confocal PSF, $H_{in} \times H_{out}$.

The de-scanned elements allow us to go far beyond standard confocal imaging. In particular, they will be exploited to unscramble the local input and output PSFs in the vicinity of each focal point. As a preliminary step, they can also be used to quantify the level of aberrations and multiple scattering. In average, the de-scanned intensity, $I(\Delta\rho_{in}, \rho_{out}, z) = |R_{in}(\Delta\rho, \rho_{out}, z)|^2$, can actually be expressed as the convolution between the incoherent input and output PSFs$^{33}$:

$$\langle I(\Delta\rho_{in}, \rho_{out}, z) \rangle \propto |H_{in}|^2 \stackrel{\Delta\rho_{in}}{\otimes} |H_{out}|^2 (\Delta\rho_{in}, \rho_{out}, z) \quad (3)$$

where the symbol $\otimes$ stands for correlation product and $\langle \cdots \rangle$ for ensemble average. This quantity will be referred to as RPSF in the following (acronym for reflection PSF). Figure S1e displays examples of RSPF extracted in depth of the opaque cornea. The spatial extension $\delta_R$ of the RPSF indicates the focusing quality and dictates the number $M$ of central lines of $\mathbf{R}_{in}(z)$ that contain the relevant information for imaging:

$$M \sim (\delta_R/\delta_0)^2 \quad (4)$$

with $\delta_0 \sim \lambda/(4NA)$, the confocal maximal resolution of the imaging system. For a field-of-view much larger than the spatial extension of the RPSF ($\Omega \gg \delta_R$), the de-scanned basis is thus particularly relevant for the acquisition of $\mathbf{R}$ ($M \ll N$).
However, while a raster scan (column-by-column acquisition) of $\mathbf{R}$ requires to illuminate the sample over a field-of-view $\Omega$ with $N = (\Omega/\delta_0)^2$ input wave-fronts $^{20,31,32}$, the de-scanned basis allows a much smaller number of field measurements.

This sparsity can be understood by expressing theoretically the de-scan matrix $\mathbf{R}_{\text{in}}$ (Supplementary Section S3):

$$R_{\text{in}}(\Delta\boldsymbol{\rho}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, z) = \int_{\Omega} d\boldsymbol{\rho}_{\text{s}} H_{\text{in}}(\boldsymbol{\rho}_{\text{s}} + \Delta\boldsymbol{\rho}_{\text{in}}, \boldsymbol{\rho}_{\text{in}}, z)\gamma(\boldsymbol{\rho}_{\text{s}} + \boldsymbol{\rho}_{\text{out}}, z)H_{\text{out}}(\boldsymbol{\rho}_{\text{s}}, \boldsymbol{\rho}_{\text{out}}, z) \quad (2)$$

where $\gamma$ is the sample reflectivity. $H_{\text{in}}(\boldsymbol{\rho}_{\text{s}}, \boldsymbol{\rho}_{\text{in}}, z)$ and $H_{\text{out}}(\boldsymbol{\rho}_{\text{s}}, \boldsymbol{\rho}_{\text{out}}, z)$ are the local input and output point spread functions (PSFs) at points $(\boldsymbol{\rho}_{\text{in}}, z)$ and $(\boldsymbol{\rho}_{\text{out}}, z)$, respectively. This last equation confirms that the central line of $\mathbf{R}_{\text{in}}$ $(\Delta\boldsymbol{\rho}_{\text{in}} = \mathbf{0})$, i.e. the FFOCT image, results from a convolution between the sample reflectivity $\gamma$ and the local confocal PSF, $H_{\text{in}} \times H_{\text{out}}$.

The de-scanned elements allow us to go far beyond standard confocal imaging. In particular, they will be exploited to unscramble the local input and output PSFs in the vicinity of each focal point. As a preliminary step, they can also be used to quantify the level of aberrations and multiple scattering. In average, the de-scanned intensity, $I(\Delta\boldsymbol{\rho}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, z) = |R_{\text{in}}(\Delta\boldsymbol{\rho}, \boldsymbol{\rho}_{\text{out}}, z)|^2$, can actually be expressed as the convolution between the incoherent input and output PSFs$^{33}$:

$$\langle I(\Delta\boldsymbol{\rho}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, z) \rangle \propto |H_{\text{in}}|^2 \stackrel{\Delta\boldsymbol{\rho}_{\text{in}}}{\otimes} |H_{\text{out}}|^2 (\Delta\boldsymbol{\rho}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, z) \quad (3)$$

where the symbol $\otimes$ stands for correlation product and $\langle \cdots \rangle$ for ensemble average. This quantity will be referred to as RPSF in the following (acronym for reflection PSF). Figure 1e displays examples of RSPF extracted in depth of the opaque cornea. The spatial extension $\delta_R$ of the RPSF indicates the focusing quality and dictates the number $M$ of central lines of $\mathbf{R}_{\text{in}}(z)$ that contain the relevant information for imaging:

$$M \sim (\delta_R / \delta_0)^2 \quad (4)$$

with $\delta_0 \sim \lambda/(4NA)$, the confocal maximal resolution of the imaging system. For a field-of-view much larger than the spatial extension of the RPSF ($\Omega \gg \delta_R$), the de-scanned basis is thus particularly relevant for the acquisition of $\mathbf{R}$ ($M \ll N$).

## Quantifying the Focusing Quality
Quantifying the Focusing Quality

Figure 2d shows the depth evolution of the RPSF. It exhibits the following characteristic shape: a distorted and enlarged confocal spot due to aberrations on top of a diffuse background³³. While the latter component is due to multiple scattering, the former component contains the contribution of singly-scattered photons but also a coherent backscattering peak³³ resulting from a constructive interference between multiple scattering paths³⁴,³⁵ (Supplementary Figure S13).

Figure 2d clearly highlights two regimes. In the epithelium ($z < 70 \ \mu m$), the single scattering component is predominant and the image of the cornea is reliable although its resolution is affected by aberrations (Fig. 2b₁). Beyond this depth, the multiple scattering background is predominant and drastically blurs the image (Fig. 2b₃). The axial evolution of the single scattering rate enables the measurement of the scattering mean free path $\ell_s$ ³⁶ (Supplementary Section S4). We find $\ell_s \sim 35 \ \mu m$ in the stroma (Fig. 2g), which confirms the strong opacity of the cornea. The penetration depth limit thus scales as $\ell_s$. This value is modest compared with theoretical predictions⁸ ($\sim 4\ell_s$) but is explained by the occurrence of strong aberrations at shallow depths, partially due to the index mismatch at the cornea surface (Fig. 2d).

The RSPF also fluctuates in the transverse direction. To that aim, a map of local RPSFs (Fig. 4c) can be built by considering the back-scattered intensity over limited spatial windows (Methods). This map shows important fluctuations due to: (i) the variations of the medium reflectivity that acts on the level of the confocal spot with respect to the diffuse background; (ii) the lateral variations of the optical index upstream of the focal plane that induce distortions of the confocal peak. Such complexity implies that any point in the medium will be associated with its own distinct focusing law. Nevertheless, spatial correlations subsist between RSPFs in adjacent windows (Fig. 4c). Such correlations can be explained by a physical phenomenon often referred to as isoplanatism in AO³⁷ and that results in a locally-invariant PSF³⁸. We will now see how this local isoplanicity can be exploited for the estimation of the T-matrices.
Figure 2d shows the depth evolution of the RPSF. It exhibits the following characteristic shape: a distorted and enlarged confocal spot on top of a diffuse background³³. The former component is a manifestation of aberrations; the latter contribution is due to multiple scattering. Figure 2d clearly highlights two regimes. In the epithelium ($z < 70$ µm), the confocal component is predominant and the image of the cornea is reliable although its resolution is affected by aberrations (Fig. 2b₁). Beyond this depth, the multiple scattering background is predominant and drastically blurs the image (Fig. 2b₃). The axial evolution of the confocal-to-multiple scattering ratio enables the measurement of the scattering mean free path $\ell_s$³⁴ (Supplementary Section S2). We find $\ell_s \sim 35$ µm in the stroma (Fig. 2g), which confirms the strong opacity of the cornea. The penetration depth limit thus scales as $\ell_s$. This value is modest compared with theoretical predictions⁸ ($\sim 4\ell_s$) but is explained by the occurrence of strong aberrations at shallow depths, partially due to the index mismatch at the cornea surface (Fig. 2d).

The RSPF also fluctuates in the transverse direction. To that aim, a map of local RPSFs (Fig. 4c) can be built by considering the back-scattered intensity over limited spatial windows (Methods). This map shows important fluctuations due to: (i) the variations of the medium reflectivity that acts on the level of the confocal spot with respect to the diffuse background; (ii) the lateral variations of the optical index upstream of the focal plane that induce distortions of the confocal peak. Such complexity implies that any point in the medium will be associated with its own distinct focusing law. Nevertheless, spatial correlations subsist between RSPFs in adjacent windows (Fig. 4c). Such correlations can be explained by a physical phenomenon often referred to as isoplanatism in AO³⁵ and that results in a locally-invariant PSF³⁶. We will now see how this local isoplanicity can be exploited for the estimation of the T-matrices.

## Iterative Phase Reversal of Wave Distortions

To that aim, we will exploit and extend the distortion matrix concept introduced in a previous work¹⁹. Interestingly, a Fourier transform over the coordinate $\Delta\rho_{in}$ of each de-scanned wave-field, $R_{in}(\Delta\rho_{in}, \rho_{out}, z)$, actually yields the wave distortions seen from the
<figure><img src="image_4.png" /><figcaption>FIG. 4. Time-gated Transmission Matrix for Local Compensation of Forward Multiple Scattering. a,b. Confocal field of view before and after the correction process at 200 µm-depth, respectively (scale bar: 50 µm). c,d. Maps of the local reflection point-spread functions (RPSFs) (de-scan field-of-view: 7 × 7 µm²) over the field of view, before and after the correction process, respectively. e,f. Sub-part of matrices, $T_{\text{in}}$ and $T_{\text{out}}$, respectively, for the area delimited by the square box in panels a-d.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 4. Time-gated Transmission matrix for Local Compensation of Forward Multiple Scattering. a,b. Confocal field of view before and after the correction process at 200 µm-depth, respectively (scale bar: 50 µm). c,d. Maps of the local reflection point-spread functions (RPSFs) (de-scan field-of-view: 7 × 7 µm²) over the field of view, before and after the correction process, respectively. e,f. Sub-part of matrices, $\hat{T}_{in}$ and $\hat{T}_{out}$, respectively, for the area delimited by the square box in panels a-d.</figcaption></figure>

input pupil plane (Fig. 3d) :

$$D_{in}(z) = F \times R_{in}(z) \qquad (5)$$

where $F$ denotes the Fourier transform operator, $F(u, \Delta\rho) = \exp(-i2\pi u \cdot \Delta\rho / \lambda f)$, $\lambda$ the central wavelength and $f$ the MO focal length. $D_{in}(z) = [D(u_{in}, \rho_{out}, z)]$ is the distortion matrix that connects any voxel $(\rho_{out}, z)$ in the field-of-view to wave-distortions in the input pupil plane $(u_{in})$.

As expected in most of biological tissues, this matrix exhibits local correlations that can be understood in light of the *shift-shift* memory effect³⁶,³⁷: Waves produced by nearby points inside an anisotropic scattering medium generate highly correlated random speckle patterns in the pupil plane. Figure 1 illustrates this fact by displaying an example of distortion matrix
Iterative Phase Reversal of Wave Distortions

To that aim, we will exploit and extend the distortion matrix concept introduced in a previous work¹⁹. Interestingly, a Fourier transform over the coordinate $\Delta\rho_{in}$ of each de-scanned wave-field, $R_{in}(\Delta\rho_{in}, \rho_{out}, z)$, actually yields the wave distortions seen from the input pupil plane (Fig. 3d) :

$$D_{in}(z) = T_{0} \times R_{in}(z) \qquad (5)$$

where $T_0$ denotes the Fourier transform operator, $T_0(\mathbf{u}, \Delta\rho) = \exp(-i2\pi\mathbf{u} \cdot \Delta\rho / \lambda f)$, $\lambda$ the central wavelength and $f$ the MO focal length. $D_{in}(z) = [D(\mathbf{u}_{in}, \rho_{out}, z)]$ is the distortion matrix that connects any voxel $(\rho_{out}, z)$ in the field-of-view to wave-distortions in the input pupil plane $(\mathbf{u}_{in})$.

As expected in most of biological tissues, this matrix exhibits local correlations that can be understood in light of the shift-shift memory effect³⁸,³⁹: Waves produced by nearby points inside an anisotropic scattering medium generate highly correlated random speckle patterns in the pupil plane. Figure S1 illustrates this fact by displaying an example of distortion matrix (Fig. S1g) and reshaped distorted wave-fields for different points $(\rho_{out}, z)$ (Fig. S1h). A strong similarity can be observed between distorted wave-fronts associated with neighboring points but this correlation tends to vanish when the two points are too far away.

The next step is to extract and exploit this local memory effect for imaging. To that aim, a set of correlation matrices $C_{in}(\mathbf{r}_p)$ shall be considered between distorted wave-fronts in the vicinity of each point $\mathbf{r}_p$ in the field-of-view (Methods). Under the hypothesis of local isoplanicity, each matrix $C_{in}(\mathbf{r}_p)$ is analogous to a $\mathbf{R}$-matrix associated with a virtual reflector synthesized from the set of output focal spots²¹ (see Fig. 3e and Supplementary Section S2). In this fictitious experimental configuration, an iterative phase-reversal (IPR) process can be performed to converge towards the incident wave front that focuses perfectly through the heterogeneities of the medium onto this virtual guide star (see Fig. 3f and Methods).

IPR repeated for each point $\mathbf{r}_p$ yields a set of pupil phase laws $\mathcal{T}_{in}(\mathbf{u}, \mathbf{r}_p)$ forming the transmittance matrix $\mathcal{T}_{in}$. Its digital phase conjugation enables a local compensation of aberration and forward multiple scattering. An updated de-scanned matrix can then be built:

$$R_{in} = T_{0}^{\dagger} \times [\mathcal{T}_{in}^{*} \circ D_{in}] \qquad (6)$$
(Fig. 1g) and reshaped distorted wave-fields for different points ($\rho_{out}, z$) (Fig. 1h). A strong similarity can be observed between distorted wave-fronts associated with neighboring points but this correlation tends to vanish when the two points are too far away.

The next step is to extract and exploit this local memory effect for imaging. To that aim, a set of correlation matrices $C_{in}(\mathbf{r}_p)$ shall be considered between distorted wave-fronts in the vicinity of each point $\mathbf{r}_p$ in the field-of-view (Methods). Under the hypothesis of local isoplanicity, each matrix $C_{in}(\mathbf{r}_p)$ is analogous to a $\mathbf{R}$-matrix associated with a virtual reflector synthesized from the set of output focal spots²¹ (see Fig. 3e and Supplementary Section S6). In this fictitious experimental configuration, an iterative phase-reversal (IPR) process can be performed to converge towards the incident wave front that focuses perfectly through the heterogeneities of the medium onto this virtual guide star (see Fig. 3f and Methods).

IPR repeated for each point $\mathbf{r}_p$ yields an estimator $\hat{\mathbf{T}}_{in}$ of the time-gated transmission matrix, $\mathbf{T}_{in} = \mathbf{F} \times \mathbf{H}_{in}$. Its digital phase conjugation enables a local compensation of aberration and forward multiple scattering. An updated de-scanned matrix can then be built:

$$ \mathbf{R}_{in} = \mathbf{F}^{\dagger} \times [\hat{\mathbf{T}}_{in}^{*} \circ \mathbf{D}_{in}] \quad (6) $$

where the symbol $\dagger$ stands for transpose conjugate and $\circ$ for the Hadamard product. The same process can be repeated by exchanging input and output to estimate the output transmission matrix $\mathbf{T}_{out}$ (Methods).

## Multi-Scale Analysis of the Distortion Matrix

A critical aspect of RMI is the choice of the spatial window over which wave distortions shall be analyzed. On the one hand, the isoplanatic assumption is valid for low-order aberrations that are associated with extended isoplanatic patches. On the other hand, forward multiple scattering gives rise to high-order aberrations that exhibit a coherence length that decreases with depth until reaching the size of a speckle grain beyond $\ell_t^{36}$. However, each spatial window should be large enough to encompass a sufficient number of independent realizations of disorder³⁸. Indeed, the bias of our $\mathbf{T}$-matrix estimator scales as follows (see Supplementary Section S9):

$$ |\delta T(\mathbf{u}, \mathbf{r}_p)|^2 \sim 1/(C^2 N_W) \quad (7) $$
where the symbol † stands for transpose conjugate and ◦ for the Hadamard product. The same process can be repeated by exchanging input and output to estimate the output transmittance matrix $\mathcal{T}_{\text{out}}$ (Methods). The element wise product between the free space transmission matrix $\mathbf{T}_0$ and the transmittance matrix $\mathcal{T}$ constitutes an estimator of the time-gated transmission matrix $\mathbf{T}$. The latter matrix contains the impulse responses $T(\mathbf{u}, \mathbf{r})$ between the pupil plane $\mathbf{u}$ and each voxel $\mathbf{r}$ inside the medium around the ballistic time $\tau_B$. Note that this matrix not only contains a ballistic (possibly aberrated) component but also grasps forward multiple scattering paths which display a time-of-flight in the same coherence time as ballistic photons. In the following, we show how these complex trajectories can be harnessed thanks to RMI.

## Multi-Scale Analysis of the Distortion Matrix

To that aim, a critical aspect is the choice of the spatial window over which wave distortions shall be analyzed. On the one hand, the isoplanatic assumption is valid for low-order aberrations that are associated with extended isoplanatic patches. On the other hand, forward multiple scattering gives rise to high-order aberrations that exhibit a coherence length that decreases with depth until reaching the size of a speckle grain beyond $\ell_t^{38}$. However, each spatial window should be large enough to encompass a sufficient number of independent realizations of disorder$^{40}$. Indeed, the bias of our $\mathbf{T}$-matrix estimator scales as follows (see Supplementary Section S3):

$$|\delta\mathcal{T}(\mathbf{u}, \mathbf{r}_p)|^2 \sim 1/(C^2 N_W) \qquad (7)$$

with $N_W$ the number of resolution cells in each spatial window. $C$ is a coherence factor that is a direct indicator of the focusing quality$^{41}$.

To limit this bias while addressing the scattering component of $\mathcal{T}$, an iterative multi-scale analysis of $\mathbf{D}$ is proposed (Methods). It consists in gradually reducing the size of the virtual guide star by: (i) alternating the correction at input and output (Supplementary Section S3); (ii) dividing by two the size of overlapping spatial windows at each iterative step (Fig. 5a). Thereby the RPSF extension is gradually narrowed (Fig. 5b) and the coherence factor $C$ increased. The spatial window can thus be reduced accordingly at the next step while maintaining an acceptable bias (Eq. S3 C). It enables the capture of finer angular and
with $N_W$ the number of resolution cells in each spatial window. $C$ is a coherence factor that is a direct indicator of the focusing quality³⁹, ranging from 0 for a fully blurred guide star to 2/9 for a diffraction-limited focal spot⁴⁰.

To limit this bias while addressing the scattering component of $T$, an iterative multi-scale analysis of $D$ is proposed (Methods). It consists in gradually reducing the size of the virtual guide star by: (i) alternating the correction at input and output (Supplementary Section S11); (ii) dividing by two the size of overlapping spatial windows at each iterative step (Fig. 5a). Thereby the RPSF extension is gradually narrowed (Fig. 5b) and the coherence factor $C$ increased. The spatial window can thus be reduced accordingly at the next step while maintaining an acceptable bias (Eq. S40). It enables the capture of finer angular and spatial details of the $T$-matrix at each step (Fig. 5c) while ensuring the convergence of IPR. As discussed further, the end of the process is monitored by the memory effect that shall exhibit our $T$-matrix estimator (Supplementary Section S12).

<figure><img src="image_6.png" /><figcaption>FIG. 5. **Multi-scale analysis of wave distortions.** a. The entire field-of-view is 138 × 138 µm². At each step, it is divided into a set of spatial windows whose dimension gradually decreases: from 138, 100, 50, 25, 13 to 6 µm b. Evolution of the transmitted wave-front $\hat{T}(\mathbf{u}_{\text{out}}, \mathbf{r}_{\text{p}})$ for one point $\mathbf{r}_{\text{p}}$ of the field-of-view at each iteration step. c. Corresponding local RPSF at $\mathbf{r}_{\text{p}}$ before and after compensation of aberration and scattering using digital phase-conjugation of the optical transfer function displayed in panel b (scale bar: 2 µm). Data are from the cross-section at 200 µm depth within the sample.</figcaption></figure>
spatial details of the **T**-matrix at each step (Fig. 5c) while ensuring the convergence of IPR.
As discussed further, the end of the process is monitored by the memory effect that shall
exhibit the **T**-matrix (Supplementary Section S3). The whole process is validated by a
reference imaging experiment on a resolution target placed behind an opaque tissue layer
(Supplementary Figure S8).

<figure><img src="image_5.png" /><figcaption>FIG. 5. **Multi-scale analysis of wave distortions.** **a.** The entire field-of-view is 138 × 138 µm². At each step, it is divided into a set of spatial windows whose dimension gradually decreases: from 138, 100, 50, 25, 13 to 6 µm **b.** Evolution of the pupil transmittance $T(\mathbf{u}_{\text{out}}, \mathbf{r}_{\text{p}})$ for one point $\mathbf{r}_{\text{p}}$ of the field-of-view at each iteration step. **c.** Corresponding local RPSF at $\mathbf{r}_{\text{p}}$ before and after compensation of aberration and scattering using digital phase-conjugation of the optical transfer function displayed in panel **b** (scale bar: 2 µm). Data are from the cross-section at 200 µm depth within the sample.</figcaption></figure>

Transmittance Matrix and Memory Effect

Figures 4e and f show a sub-part of the **T**-matrices measured at depth *z* = 200 µm for final patches of 6 × 6 µm². Spatial reciprocity should imply equivalent input and output aberration phase laws. This property is not checked by our estimators. Indeed, the input aberration phase law accumulates not only the input aberrations of the sample-arm but also those of the reference arm (Supplementary Section S4). Therefore, the sample-induced aberrations can be investigated independently from the imperfections of the experimental set
Transmission Matrix and Memory Effect

Figures 4e and f show a sub-part of the T-matrices measured at depth z = 200 µm for final patches of 6 × 6 µm². Spatial reciprocity should imply equivalent input and output aberration phase laws. This property is not checked by our estimators. Indeed, the input aberration phase law accumulates not only the input aberrations of the sample-arm but also those of the reference arm (Supplementary Section S14). Therefore, the sample-induced aberrations can be investigated independently from the imperfections of the experimental set up by considering the output matrix $\hat{T}_{\text{out}}$.

An analysis of its spatial correlations³⁸ (Methods) and of its phase dislocations⁴¹ (Supplementary Section S15) shows that wave distortions induced by the cornea are made of two contributions : (i) an almost spatially-invariant and curl-free aberrated component (Fig. 6a) associated with long-scale fluctuations of the refractive index (Fig. 6c); (ii) a forward multiple scattering component (Fig. 6d) exhibiting a high density of optical vortices⁴² and a short-range memory effect whose extension drastically decreases in depth (Figs. 6a,e). The access to the latter contribution fundamentally differentiates RMI from conventional AO that only provides an access to the irrotational component of wave distortions⁴² (Supplementary Section S16).

The memory effect is also a powerful tool to monitor the convergence of the IPR process. When the spatial window is too small (3×3 µm²), IPR provides a spatially-incoherent estimator $\hat{T}$ and leads to a bucket-like image (Supplementary Figure S8). This observable thus indicates when the convergence towards $T$ is fulfilled or when the algorithm shall be stopped.

Deep Volumetric Imaging

Eventually, the estimated T-matrices can be used to compensate for local aberrations over the whole field-of-view. To that aim, a digital phase conjugation is performed at input and output (Eq. 6). The comparison between the initial and resulting images (Figs. 4a,b) demonstrates the benefit of a local compensation of aberration and scattering. The drastic gain in resolution and contrast provided by RMI enables to reveal a rich arrangement of biological structures (cells, striae, etc.) that were completely blurred by scattering in the
up by considering the output matrix $\mathcal{T}_{\text{out}}$.

An analysis of its spatial correlations$^{40}$ (Methods) and its angular decomposition (Supplementary Figure S12) shows that wave distortions induced by the cornea are made of two contributions : (i) an almost spatially-invariant aberrated component (Fig. 6a) associated with long-scale fluctuations of the refractive index (Fig. 6c); (ii) a forward multiple scattering component (Fig. 6d) giving rise to an angular dispersion of photons between the cornea surface and the focal plane. The latter component is associated with a short-range memory effect whose extension drastically decreases in depth (Figs. 6a,e). The access to this contribution fundamentally differentiates RMI from conventional AO that only provides an access to the irrotational component of wave distortions$^{42}$ (Supplementary Section S4).

<figure><img src="image_6.png" /><figcaption>FIG. 6. Revealing the memory effect exhibited by the **T**-matrix. a. Transverse evolution of the mean correlation function of the transmitted wave-field from shallow (blue) to large (red) depths. b. The phase of each transmitted wave-field is the sum of: c. a spatially-invariant aberration phase function; d. a complex scattering law exhibiting high spatial frequencies. e. The spatial correlation of the latter component with the **T**-matrix provides a map of the corresponding isoplanatic patch (scale bar: 50 µm).</figcaption></figure>

The memory effect is also a powerful tool to monitor the convergence of the IPR process. When the spatial window is too small ($3 \times 3 \mu m^2$), IPR provides a spatially-incoherent **T**-matrix and leads to a bucket-like image (Supplementary Figure S7). This observable thus indicates when the convergence towards **T** is fulfilled or when the algorithm shall be stopped.
<figure><img src="image_7.png" /><figcaption>FIG. 6. Revealing the memory effect exhibited by the T-matrix. a. Transverse evolution of the mean correlation function of the transmitted wave-field from shallow (blue) to large (red) depths. b. The phase of each transmitted wave-field is the sum of: c. a spatially-invariant aberration phase function; d. a complex scattering law exhibiting high spatial frequencies. e. The spatial correlation of the latter component with the T-matrix provides a map of the corresponding isoplanatic patch (scale bar: 50 µm).</figcaption></figure>

initial image. For instance, a stromal stria, indicator of keratoconus$^{43}$, is clearly revealed on the RMI B-scan (Fig. 2f) while it was hidden by the multiple scattering fog on the initial image (Fig. 2c). The B-scan shows that RMI provides a full image of the cornea with the recovery of its different layers throughout its thickness (350 µm ~ 10ℓs, see also Supplementary Movies).

The gain in contrast and resolution can be quantified by investigating the RSPF after RMI. A close-to-ideal confocal resolution (230 nm vs. δ₀ ~ 215 nm) is reached throughout the cornea thickness (Fig. 2e). The confocal-to-diffuse ratio is increased by a factor up to 15 dB in depth (Supplementary Section S13). Furthermore, the map of local RPSFs displayed in Fig. 4d shows the efficiency of RMI for addressing extremely small isoplanatic patches.

Discussion

In this experimental proof-of-concept, we demonstrated the capacity of RMI to exploit
forward multiple scattering for deep imaging of biological tissues. This work introduces
several crucial elements, thereby leading to a better imaging performance than previous
studies.
Deep Volumetric Imaging

Eventually, the *T*-matrix can be used to compensate for local aberrations over the whole field-of-view. To that aim, a digital phase conjugation is performed at input and output (Eq. 6). The comparison between the initial and resulting images (Figs. 4a,b) demonstrates the benefit of a local compensation of aberration and scattering. The drastic gain in resolution and contrast provided by RMI enables to reveal a rich arrangement of biological structures (cells, striae, *etc*. ) that were completely blurred by scattering in the initial image. For instance, a stromal stria, indicator of keratoconus$^{43}$, is clearly revealed on the RMI B-scan (Fig. 2f) while it was hidden by the multiple scattering fog on the initial image (Fig. 2c). The B-scan shows that RMI provides a full image of the cornea with the recovery of its different layers throughout its thickness (350 µm ~ 10ℓs, see also Supplementary Movies).

The gain in contrast and resolution can be quantified by investigating the RSPF after
RMI. A close-to-ideal confocal resolution (230 nm vs. δ₀ ~ 215 nm) is reached throughout
the cornea thickness (Fig. 2e). The confocal-to-diffuse ratio is increased by a factor up to
15 dB in depth (Supplementary Section S4). Furthermore, the map of local RPSFs displayed
in Fig. 4d shows the efficiency of RMI for addressing extremely small isoplanatic patches.

Discussion

In this experimental proof-of-concept, we demonstrated the capacity of RMI to exploit
forward multiple scattering for deep imaging of biological tissues. This work introduces
several crucial elements, thereby leading to a better imaging performance than previous
studies.

First, the proposed IPR algorithm outperforms iterative time reversal processing<sup>19</sup> for
local compensation of aberrations in scattering media because it can evaluate the focusing
laws over a larger angular domain (Supplementary Figure S2). Second, the bias of our
**T**-matrix estimator has been expressed analytically (Eq. S3 C) as a function of a coherence
factor that grasps the blurring effect of aberrations and multiple scattering. This led us to
define a multi-scale strategy for matrix imaging with a fine monitoring of its convergence
based on the memory effect. The latter observable is a real asset as it provides an objective
criterion to: (i) optimize the resolution of our **T**-matrix estimator (Supplementary Section
First, the proposed IPR algorithm outperforms iterative time reversal processing<sup>19</sup> for local compensation of aberrations in scattering media because it can evaluate the focusing laws over a larger angular domain (Supplementary Figure S3). Second, the bias of our **T**-matrix estimator has been expressed analytically (Eq. S40) as a function of a coherence factor that grasps the blurring effect of aberrations and multiple scattering. This led us to define a multi-scale strategy for matrix imaging with a fine monitoring of its convergence based on the memory effect. The latter observable is a real asset as it provides an objective criterion to: (i) optimize the resolution of our **T**-matrix estimator (Supplementary Section S12); (ii) compare our approach with alternative methods such as the CLASS algorithm<sup>20,23,24</sup> (Supplementary Section S18). Our multi-scale process enables us to target isoplanatic areas more than four times smaller than CLASS. Interestingly, those two approaches are based on the maximization of different physical quantities: the confocal intensity for CLASS; the coherence of the wave-field induced by a virtual guide star for IPR. Hence they are, in principle, perfectly complementary and could be advantageously combined in the future.

Although this experimental proof-of-concept is promising for deep optical imaging of biological tissues, it also suffers from several limitations that need to be addressed in future works. First, FFOCT is not very convenient for 3D in-vivo imaging since it requires an axial scan of the sample. Another possibility would be to move the reference arm and measure **R** as a function of the time-of-flight. An access to the time (or spectral) dependence of the **R**-matrix is actually critical to reach a larger penetration depth. Indeed, the focusing law extracted from a time-gated **R**-matrix is equivalent in the time domain to a simple application of time delays between each angular component of the wave-field. Yet, the diffusive regime requires to address independently each frequency component of the wave-field to make multiple scattering paths of different lengths constructively interfere on any focusing point in depth. On the one hand, the exploitation of the *chromato-axial* memory effect<sup>44</sup> will be decisive to ensure the convergence of IPR over isoplanatic volumes<sup>45</sup>. On the other hand, the *tilt-tilt* memory effect<sup>37</sup> can also be leveraged by investigating the distortion matrix, not only in the pupil plane, but in any plane lying between the medium surface and the focal plane, thereby mimicking a multi-conjugate AO scheme<sup>46</sup>.

Beyond the diffusive regime, another blind spot of this study is the medium movement during the experiment<sup>47,48</sup>. In that respect, the matrix formalism shall be developed to include the medium dynamics. Moving speckle can actually be an opportunity since it
S3); (ii) compare our approach with alternative methods such as the CLASS algorithm<sup>20,23,24</sup>
(Supplementary Section S5). Our multi-scale process enables us to target isoplanatic areas
more than four times smaller than CLASS. Interestingly, those two approaches are based
on the maximization of different physical quantities: the confocal intensity for CLASS; the
coherence of the wave-field induced by a virtual guide star for IPR. Hence they are, in
principle, perfectly complementary and could be advantageously combined in the future.

Although this experimental proof-of-concept is promising for deep optical imaging of
biological tissues, it also suffers from several limitations that need to be addressed in future
works. First, FFOCT is not very convenient for 3D in-vivo imaging since it requires an axial
scan of the sample. Another possibility would be to move the reference arm and measure
**R** as a function of the time-of-flight. An access to the time (or spectral) dependence of
the **R**-matrix is actually critical to reach a larger penetration depth. Indeed, the focusing
law extracted from a time-gated **R**-matrix is equivalent in the time domain to a simple
application of time delays between each angular component of the wave-field. Yet, the
diffusive regime requires to address independently each frequency component of the wave-field
to make multiple scattering paths of different lengths constructively interfere on any focusing
point in depth. On the one hand, the exploitation of the chromato-axial memory effect<sup>44</sup> will
be decisive to ensure the convergence of IPR over isoplanatic volumes<sup>45</sup>. On the other hand,
the tilt-tilt memory effect<sup>39</sup> can also be leveraged by investigating the distortion matrix, not
only in the pupil plane, but in any plane lying between the medium surface and the focal
plane, thereby mimicking a multi-conjugate AO scheme<sup>46</sup>.

Beyond the diffusive regime, another blind spot of this study is the medium movement
during the experiment<sup>47,48</sup>. In that respect, the matrix formalism shall be developed to
include the medium dynamics. Moving speckle can actually be an opportunity since it
can give access to a large number of speckle realizations for each voxel. A high resolution
**T**-matrix could be, in principle, extracted without relying on any isoplanatic assumption<sup>49</sup>.

To conclude, this study is a striking illustration of a pluri-disciplinary approach in wave
physics. A passive measurement of the **R**-matrix is indeed an original idea coming from
seismology<sup>50</sup>. The **D**-matrix is inspired by stellar speckle interferometry in astronomy<sup>51</sup>.
The **T**-matrix is a concept that has emerged both from fundamental studies in condensed
matter physics<sup>52</sup> and more applied fields such as MIMO communications<sup>53</sup> and ultrasound
therapy<sup>12</sup>. The emergence of high-speed cameras and the rapid growth of computational
can give access to a large number of speckle realizations for each voxel. A high resolution
**T**-matrix could be, in principle, extracted without relying on any isoplanatic assumption⁴⁹.

To conclude, this study is a striking illustration of a pluri-disciplinary approach in wave
physics. A passive measurement of the **R**-matrix is indeed an original idea coming from
seismology⁵⁰. The **D**-matrix is inspired by stellar speckle interferometry in astronomy⁵¹.
The **T**-matrix is a concept that has emerged both from fundamental studies in condensed
matter physics⁵² and more applied fields such as MIMO communications⁵³ and ultrasound
therapy¹². The emergence of high-speed cameras and the rapid growth of computational
capabilities now makes matrix imaging mature for deep in-vivo optical microscopy.
capabilities now makes matrix imaging mature for deep in-vivo optical microscopy.

# Methods

## Experimental set up

The full experimental setup is displayed in Supplementary Figure S1. It is made of two parts: (i) a polarized Michelson interferometer illuminated by a broadband LED source (Thorlabs M850LP1, $\lambda_0$ = 850 nm, $\Delta\lambda$ = 35 nm) in a pseudo-Kohler configuration, thereby providing at its output two identical spatially-incoherent and broadband wave-fields of orthogonal polarization, the reference one being shifted by a lateral position $\Delta\rho_{in}$ by tilting the mirror in the corresponding arm; (ii) a polarized Linnik interferometer with microscope objectives (Nikon N60X-NIR, M = 60×, NA = 1.0) in the two arms and a CMOS camera (Adimec Quartz 2A-750, 2Mpx) at its output. The de-scanned beam at the output the first interferometer illuminates the reference arm of the second interferometer and is reflected by the reference mirror placed in the focal plane of the MO. The other beam at the output of the first interferometer illuminates the sample placed in the focal plane of the other MO. The CMOS camera, conjugated with the focal planes of the MO, records the interferogram between the beams reflected by each arm of the Linnik interferometer. The spatial sampling of each recorded image is $\delta_0$ = 230 nm and the field-of-view is 275 × 275 µm$^2$.

## Cornea

The human cornea under study is a pathological surgical specimen that was provided by the Quinze-Vingts National Eye Hospital operating room at the time of keratoplasty. The use of such specimens was approved by the Institutional Review Board (Patient Protection Committee, Ile-de-France V) and adhered to the tenets of the Declaration of Helsinki as well as to international ethical requirements for human tissues. The ethics committee waived the requirement for informed written consent of patient; however, the patient provided informed oral consent to have their specimen used in research.

## Experimental procedure

The experiment consists in the acquisition of the de-scanned reflection matrix $R_{in}$. To that aim, an axial scan of the sample is performed over the cornea thickness (350 µm) with a
# Methods

## Experimental set up

The full experimental setup is displayed in Supplementary Figure S1. It is made of two parts: (i) a polarized Michelson interferometer illuminated by a broadband LED source (Thorlabs M850LP1, $\lambda_0$ = 850 nm, $\Delta\lambda$ = 35 nm) in a pseudo-Kohler configuration, thereby providing at its output two identical spatially-incoherent and broadband wave-fields of orthogonal polarization, the reference one being shifted by a lateral position $\Delta\rho_{in}$ by tilting the mirror in the corresponding arm; (ii) a polarized Linnik interferometer with microscope objectives (Nikon N60X-NIR, M = 60×, NA = 1.0) in the two arms and a CMOS camera (Adimec Quartz 2A-750, 2Mpx) at its output. The de-scanned beam at the output the first interferometer illuminates the reference arm of the second interferometer and is reflected by the reference mirror placed in the focal plane of the MO. The other beam at the output of the first interferometer illuminates the sample placed in the focal plane of the other MO. The CMOS camera, conjugated with the focal planes of the MO, records the interferogram between the beams reflected by each arm of the Linnik interferometer. The spatial sampling of each recorded image is $\delta_0$ = 230 nm and the field-of-view is 275 × 275 $\mu$m$^2$.

## Cornea

The human cornea under study is a pathological surgical specimen that was provided by the Quinze-Vingts National Eye Hospital operating room at the time of keratoplasty. The use of such specimens was approved by the Institutional Review Board (Patient Protection Committee, Ile-de-France V) and adhered to the tenets of the Declaration of Helsinki as well as to international ethical requirements for human tissues. The ethics committee waived the requirement for informed written consent of patient; however, the patient provided informed oral consent to have their specimen used in research.

## Experimental procedure

The experiment consists in the acquisition of the de-scanned reflection matrix $R_{in}$. To that aim, an axial scan of the sample is performed over the cornea thickness (350 $\mu$m) with a sampling of 2 $\mu$m (i.e 185 axial positions). For each depth, a transverse scan of the de-scanned position $\Delta\rho_{in}$ is performed over a 2.9 × 2.9 $\mu$m$^2$ area with a spatial sampling $\delta_0$ = 230 nm
sampling of 2 µm (i.e 185 axial positions). For each depth, a transverse scan of the de-scanned position $\Delta\rho_{in}$ is performed over a 2.9 × 2.9 µm² area with a spatial sampling $\delta_0 = 230$ nm (that is to say 169 input wave-fronts instead of $10^6$ input wave-fronts in a canonical basis). For each scan position ($\Delta\rho, z$), a complex-reflected wave field is extracted by phase shifting interferometry from four intensity measurements. This measured field is averaged over 5 successive realisations (for denoising). The integration time of the camera is set to 5 ms. Each wave-field is stored in the de-scanned reflection matrix $R_{in} = [R_{in}(\Delta\rho_{in}, \rho_{out})]$ (Fig. S1). The duration time for the recording of $R_{in}$ is of ~ 30 s at each depth. The post-processing of the reflection matrix (IPR and multi-scale analysis) to get the final image took only a few minutes on Matlab. The experimental results displayed in Fig. 4 and 5 at a single depth $z = 200$ µm have been obtained by performing a de-scan over a 7 × 7 µm² area with a spatial sampling $\delta_0 = 230$ nm (961 input wave-fronts).

## Local RPSF

To probe the local RPSF, the field-of-view is divided into regions that are defined by their central midpoint $\mathbf{r}_p = (\rho_{p}, z)$ and their lateral extension $L$. A local average of the back-scattered intensity can then be performed in each region:

$$I(\Delta\rho_{in}, \mathbf{r}_p) = \langle |R_{in}(\Delta\rho_{in}, \rho_{out}, z)|^2 W_L(\rho_{out} - \rho_{p}) \rangle_{\rho_{out}} \quad (8)$$

where $W_L(\rho_{out} - \rho_{p}) = 1$ for $|x_{out} - x_{p}| < L$ and $|y_{out} - y_{p}| < L$, and zero otherwise.

## Multi-scale compensation of wave-distortions

The multi-scale process consists in an iterative compensation of aberration and scattering phenomena at input and output of the reflection matrix. To that aim, wave distortions are analyzed over spatial windows $W_L$ that are gradually reduced at each step $q$ of the procedure, such that:

$$L = FOV/2^q \qquad (9)$$

where FOV denotes the initial field-of-view.

The whole procedure is summarized in Supplementary Figure S4. At each stage of this iterative process, the starting point is the de-scanned reflection matrix $R_{in}^{(q-1)}$, obtained at the previous step, $R_{in}^{(0)}$ being the reflection matrix recorded by our experimental set up (Fig. S1). An input distortion matrix $D_{in}^{(q)}$ is deduced from $R_{in}^{(q)}$ via a numerical Fourier
(that is to say 169 input wave-fronts instead of 10⁶ input wave-fronts in a canonical basis). For each scan position (Δρ, z), a complex-reflected wave field is extracted by phase shifting interferometry from four intensity measurements. This measured field is averaged over 5 successive realisations (for denoising). The integration time of the camera is set to 5 ms. Each wave-field is stored in the de-scanned reflection matrix $R_{in} = [R_{in}(\Delta\rho_{in}, \rho_{out})]$ (Fig. 1). The duration time for the recording of $R_{in}$ is of ~30 s at each depth. The post-processing of the reflection matrix (IPR and multi-scale analysis) to get the final image took only a few minutes on Matlab. The experimental results displayed in Fig. 4 and 5 at a single depth z = 200 µm have been obtained by performing a de-scan over a 7 × 7 µm² area with a spatial sampling δ₀ = 230 nm (961 input wave-fronts).

## Local RPSF

To probe the local RPSF, the field-of-view is divided into regions that are defined by their central midpoint $r_p = (\rho_p, z)$ and their lateral extension L. A local average of the back-scattered intensity can then be performed in each region:

$$I(\Delta\rho_{in}, r_p) = \langle |R_{in}(\Delta\rho_{in}, \rho_{out}, z)|^2 W_L(\rho_{out} - \rho_p) \rangle_{\rho_{out}} \quad (8)$$

where $W_L(\rho_{out} - \rho_p) = 1$ for $|x_{out} - x_p| < L$ and $|y_{out} - y_p| < L$, and zero otherwise.

## Multi-scale compensation of wave-distortions

The multi-scale process consists in an iterative compensation of aberration and scattering phenomena at input and output of the reflection matrix. To that aim, wave distortions are analyzed over spatial windows $W_L$ that are gradually reduced at each step q of the procedure, such that:

$$L = FOV/2^q \quad (9)$$

where FOV denotes the initial field-of-view.

The whole procedure is summarized in Supplementary Figure S5. At each stage of this iterative process, the starting point is the de-scanned reflection matrix $R_{in}^{(q-1)}$, obtained at the previous step, $R_{in}^{(0)}$ being the reflection matrix recorded by our experimental set up (Fig. 1). An input distortion matrix $D_{in}^{(q)}$ is deduced from $R_{in}^{(q)}$ via a numerical Fourier transform (Eq. 5). A local correlation matrix of wave distortions is then built around each point $r_p$ of
transform (Eq. 5). A local correlation matrix of wave distortions is then built around each point $\mathbf{r}_p$ of the field-of-view:

$$C_{\text{in}}(\mathbf{u}_{\text{in}}, \mathbf{u}'_{\text{in}}, \mathbf{r}_p) = \langle D_{\text{in}}^{(q)}(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, z) D_{\text{in}}^{(q)*}(\mathbf{u}'_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, z) W_L(\boldsymbol{\rho}_{\text{out}} - \boldsymbol{\rho}_p) \rangle_{\boldsymbol{\rho}_{\text{out}}} \quad (10)$$

IPR is then applied to each correlation matrix $C_{\text{in}}(\mathbf{r}_p)$ (see further and Supplementary Section S3). The resulting input phase laws, $\hat{\phi}_{\text{in}}(\mathbf{r}_p)$, are used to compensate for the wave distortions undergone by the incident wave-fronts:

$$\mathbf{R}'_{\text{in}} = \mathbf{T}_{0}^{\dagger} \times [\exp(-i\hat{\phi}_{\text{in}}) \circ \mathbf{D}_{\text{in}}^{(q-1)}] \quad (11)$$

The corrected matrix $\mathbf{R}'_{\text{in}}$ is only intermediate since phase distortions undergone by the reflected wave-fronts remain to be corrected.

To that aim, an output de-scanned matrix $\mathbf{R}'_{\text{out}}(z)$ is deduced from the input de-scanned matrix $\mathbf{R}'_{\text{in}}(z)$ using the following change of variable (Supplementary Figure S5):

$$R'_{\text{out}}(\boldsymbol{\rho}_{\text{in}}, \Delta\boldsymbol{\rho}_{\text{out}}, z) = R'_{\text{in}}(-\Delta\boldsymbol{\rho}_{\text{out}}, \boldsymbol{\rho}_{\text{in}} + \Delta\boldsymbol{\rho}_{\text{out}}, z) \quad (12)$$

with $\Delta\boldsymbol{\rho}_{\text{out}} = \boldsymbol{\rho}_{\text{out}} - \boldsymbol{\rho}_{\text{in}} = -\Delta\boldsymbol{\rho}_{\text{in}}$. An output distortion matrix is then built by applying a Fourier transform over the de-scanned coordinate:

$$\mathbf{D}'_{\text{out}} = \mathbf{R}'_{\text{out}} \times \mathbf{T}_{0}^{T} \quad (13)$$

where the superscript T stands for matrix transpose. From $\mathbf{D}'_{\text{out}}$, one can build a correlation matrix $C_{\text{out}}$ for each point $\mathbf{r}_p$:

$$C_{\text{out}}(\mathbf{u}_{\text{out}}, \mathbf{u}'_{\text{out}}, \mathbf{r}_p) = \langle D'_{\text{out}}(\boldsymbol{\rho}_{\text{in}}, \mathbf{u}_{\text{out}}, z_p) D'^{*}_{\text{out}}(\boldsymbol{\rho}_{\text{in}}, \mathbf{u}'_{\text{out}}, z_p) W_L(\boldsymbol{\rho}_{\text{in}} - \boldsymbol{\rho}_p) \rangle_{\boldsymbol{\rho}_{\text{in}}} \quad (14)$$

The IPR algorithm described further is then applied to each matrix $C_{\text{out}}(\mathbf{r}_p)$. The resulting output phase laws, $\hat{\phi}_{\text{out}}(\mathbf{r}_p)$, are leveraged to compensate for the residual wave distortions undergone by the reflected wave-fronts:

$$\mathbf{R}_{\text{out}}^{(q)} = [\mathbf{D}'_{\text{out}} \circ \exp(-i\hat{\phi}_{\text{out}})] \times \mathbf{T}_{0}^{*} \quad (15)$$

The RPSFs displayed in Fig. 5c are extracted from the matrices $\mathbf{R}_{\text{out}}^{(q)}$ obtained at the end of each iteration of the multi-scale process. An input de-scanned matrix, combining the input and output corrections, is finally obtained by performing the following change of variables:

$$R_{\text{in}}^{(q)}(\Delta\boldsymbol{\rho}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, z) = R_{\text{out}}^{(q)}(\boldsymbol{\rho}_{\text{out}} - \Delta\boldsymbol{\rho}_{\text{in}}, -\Delta\boldsymbol{\rho}_{\text{in}}, z) \quad (16)$$
the field-of-view:

$$C_{\text{in}}(\mathbf{u}_{\text{in}}, \mathbf{u}'_{\text{in}}, \mathbf{r}_{\text{p}}) = \langle D_{\text{in}}^{(q)}(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, z) D_{\text{in}}^{(q)*}(\mathbf{u}'_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, z) W_{L}(\boldsymbol{\rho}_{\text{out}} - \boldsymbol{\rho}_{\text{p}}) \rangle_{\boldsymbol{\rho}_{\text{out}}} \quad (10)$$

IPR is then applied to each correlation matrix $C_{\text{in}}(\mathbf{r}_{\text{p}})$ (see further and Supplementary Section S8). The resulting input phase laws, $\hat{\phi}_{\text{in}}(\mathbf{r}_{\text{p}})$, are used to compensate for the wave distortions undergone by the incident wave-fronts:

$$\mathbf{R}'_{\text{in}} = \mathbf{F}^{\dagger} \times [\exp(-i\hat{\phi}_{\text{in}}) \circ \mathbf{D}_{\text{in}}^{(q-1)}] \quad (11)$$

The corrected matrix $\mathbf{R}'_{\text{in}}$ is only intermediate since phase distortions undergone by the reflected wave-fronts remain to be corrected.

To that aim, an output de-scanned matrix $\mathbf{R}'_{\text{out}}(z)$ is deduced from the input de-scanned matrix $\mathbf{R}'_{\text{in}}(z)$ using the following change of variable (Supplementary Figure S6):

$$R'_{\text{out}}(\boldsymbol{\rho}_{\text{in}}, \Delta\boldsymbol{\rho}_{\text{out}}, z) = R'_{\text{in}}(-\Delta\boldsymbol{\rho}_{\text{out}}, \boldsymbol{\rho}_{\text{in}} + \Delta\boldsymbol{\rho}_{\text{out}}, z) \quad (12)$$

with $\Delta\boldsymbol{\rho}_{\text{out}} = \boldsymbol{\rho}_{\text{out}} - \boldsymbol{\rho}_{\text{in}} = -\Delta\boldsymbol{\rho}_{\text{in}}$. An output distortion matrix is then built by applying a Fourier transform over the de-scanned coordinate:

$$\mathbf{D}'_{\text{out}} = \mathbf{R}'_{\text{out}} \times \mathbf{F}^T \quad (13)$$

where the superscript T stands for matrix transpose. From $\mathbf{D}'_{\text{out}}$, one can build a correlation matrix $C_{\text{out}}$ for each point $\mathbf{r}_{\text{p}}$:

$$C_{\text{out}}(\mathbf{u}_{\text{out}}, \mathbf{u}'_{\text{out}}, \mathbf{r}_{\text{p}}) = \langle D'_{\text{out}}(\boldsymbol{\rho}_{\text{in}}, \mathbf{u}_{\text{out}}, z_p) D'^{*}_{\text{out}}(\boldsymbol{\rho}_{\text{in}}, \mathbf{u}'_{\text{out}}, z_p) W_{L}(\boldsymbol{\rho}_{\text{in}} - \boldsymbol{\rho}_{\text{p}}) \rangle_{\boldsymbol{\rho}_{\text{in}}} \quad (14)$$

The IPR algorithm described further is then applied to each matrix $C_{\text{out}}(\mathbf{r}_{\text{p}})$. The resulting output phase laws, $\hat{\phi}_{\text{out}}(\mathbf{r}_{\text{p}})$, are leveraged to compensate for the residual wave distortions undergone by the reflected wave-fronts:

$$\mathbf{R}_{\text{out}}^{(q)} = [\mathbf{D}'_{\text{out}} \exp(-i\phi_{\text{in}})] \times \mathbf{F}^* \quad (15)$$

The RPSFs displayed in Fig. 5c are extracted from the matrices $\mathbf{R}_{\text{out}}^{(q)}$ obtained at the end of each iteration of the multi-scale process. An input de-scanned matrix, combining the input and output corrections, is finally obtained by performing the following change of variables:

$$R_{\text{in}}^{(q)}(\Delta\boldsymbol{\rho}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, z) = R_{\text{out}}^{(q)}(\boldsymbol{\rho}_{\text{out}} - \Delta\boldsymbol{\rho}_{\text{in}}, -\Delta\boldsymbol{\rho}_{\text{in}}, z) \quad (16)$$
This matrix $\mathbf{R}_{\text{in}}^{(q)}$ is the starting point of the next stage of the multi-scale process, and so on.

The $\mathcal{T}$-matrices correspond to the cumulative function of the aberration phase laws:

$$ \mathcal{T}_{\text{in/out}}^{(q)} = \mathcal{T}_{\text{in/out}}^{(q-1)} \circ \exp \left( i \phi_{\text{in/out}}^{(q)} \right) = \prod_{k=1}^{q} \exp \left( i \phi_{\text{in/out}}^{(k)} \right) \quad (17) $$

Figure 5b shows the evolution of one line of the transmittance matrix $\mathcal{T}_{\text{out}}^{(q)}$ throughout the RMI process. The iterative procedure is stopped by investigating the correlation properties of this estimator (see further and Supplementary Section S3).

## Iterative phase reversal algorithm.

The IPR algorithm is a computational process that provides an estimator of the pupil transmittance matrix, $\mathcal{T}(\mathbf{u}, \mathbf{r}_p) = \exp[i\phi(\mathbf{u}, \mathbf{r}_p)]$, that links each point $\mathbf{u}$ of the pupil plane with each voxel $\mathbf{r}_p$ of the cornea volume. To that aim, the correlation matrix $\mathbf{C}$ computed over the spatial window $W_L$ centered around each point $\mathbf{r}_p$ is considered (Eqs. S19 and 14). Mathematically, the algorithm is based on the following recursive relation:

$$ \hat{\phi}^{(n)}(\mathbf{r}_p) = \arg \left\{ \mathbf{C}(\mathbf{r}_p) \times \exp \left[ i \hat{\phi}^{(n-1)}(\mathbf{r}_p) \right] \right\} \quad (18) $$

where $\hat{\phi}^{(n)}$ is the estimator of $\phi$ at the $n^{\text{th}}$ iteration of the phase reversal process. $\hat{\phi}^{(0)}$ is an arbitrary wave-front that initiates the process (typically a flat phase law) and $\hat{\phi} = \lim_{n \to \infty} \hat{\phi}^{(n)}$ is the result of IPR.

## Aberration and Scattering Components of the T-matrix.

The spatial correlation of transmitted wave-fields are investigated at each depth $z$ by computing the correlation matrix of $\mathcal{T}_{\text{out}}$: $\mathcal{C}_{\mathcal{T}} = \mathcal{T}_{\text{out}} \times \mathcal{T}_{\text{out}}^{\dagger}$. A mean correlation function $\Gamma$ can be computed by performing the following average:

$$ \Gamma(\Delta\boldsymbol{\rho}, z) = \langle C_{\mathcal{T}}(\boldsymbol{\rho}_{\text{in}}, \boldsymbol{\rho}_{\text{in}} + \Delta\boldsymbol{\rho}, z) \rangle_{\boldsymbol{\rho}_{\text{in}}} \quad (19) $$

The correlation function $\Gamma$ displayed in Fig. 6a shows that the matrix $\mathcal{T}_{\text{out}}$ can be decomposed as a spatially-invariant component $\mathcal{A}_{\text{out}}$ and a short-range correlated component $\mathcal{S}_{\text{out}}$. Each component can be separated by performing a singular value decomposition of $\mathcal{T}_{\text{out}}$, such that

$$ \mathcal{T}_{\text{out}} = \sum_{p=1}^{N} s_p \mathbf{U}_p \mathbf{V}_p^{\dagger} \quad (20) $$
This matrix $\mathbf{R}_{\text{in}}^{(q)}$ is the starting point of the next stage of the multi-scale process, and so on.

Estimators of the T-matrices correspond to the cumulative function of the aberration phase laws:

$$ \hat{\mathbf{T}}_{\text{in/out}}^{(q)} = \hat{\mathbf{T}}_{\text{in/out}}^{(q-1)} \circ \exp \left( i \boldsymbol{\phi}_{\text{in/out}}^{(q)} \right) = \prod_{k=1}^{q} \exp \left( i \boldsymbol{\phi}_{\text{in/out}}^{(k)} \right) \quad (17) $$

Figure 5b shows the evolution of one line of the estimator $\hat{\mathbf{T}}_{\text{out}}^{(q)}$ throughout the RMI process. The iterative procedure is stopped by investigating the correlation properties of this estimator (see further and Supplementary Section S13).

## Iterative phase reversal algorithm.

The IPR algorithm is a computational process that provides an estimator of the transmit wave-field, $T(\mathbf{u}, \mathbf{r}_p) = \exp[i\phi(\mathbf{u}, \mathbf{r}_p)]$, that links each point $\mathbf{u}$ of the pupil plane with each voxel $\mathbf{r}_p$ of the cornea volume. To that aim, the correlation matrix $\mathbf{C}$ computed over the spatial window $W_L$ centered around each point $\mathbf{r}_p$ is considered (Eqs. 10 and 14). Mathematically, the algorithm is based on the following recursive relation:

$$ \hat{\boldsymbol{\phi}}^{(n)}(\mathbf{r}_p) = \arg \left\{ \mathbf{C}(\mathbf{r}_p) \times \exp \left[ i \hat{\boldsymbol{\phi}}^{(n-1)}(\mathbf{r}_p) \right] \right\} \quad (18) $$

where $\hat{\boldsymbol{\phi}}^{(n)}$ is the estimator of $\boldsymbol{\phi}$ at the $n^{\text{th}}$ iteration of the phase reversal process. $\hat{\boldsymbol{\phi}}^{(0)}$ is an arbitrary wave-front that initiates the process (typically a flat phase law) and $\hat{\boldsymbol{\phi}} = \lim_{n \to \infty} \hat{\boldsymbol{\phi}}^{(n)}$ is the result of IPR.

## Aberration and Scattering Components of the T-matrix.

The spatial correlation of transmitted wave-fields are investigated at each depth $z$ by computing the correlation matrix of $\hat{\mathbf{T}}_{\text{out}}$: $\mathbf{C}_{\text{T}} = \hat{\mathbf{T}}_{\text{out}} \times \hat{\mathbf{T}}_{\text{out}}^{\dagger}$. A mean correlation function $\Gamma$ can be computed by performing the following average:

$$ \Gamma(\Delta\boldsymbol{\rho}, z) = \langle C_T(\boldsymbol{\rho}_{\text{in}}, \boldsymbol{\rho}_{\text{in}} + \Delta\boldsymbol{\rho}, z) \rangle_{\boldsymbol{\rho}_{\text{in}}} \quad (19) $$

The correlation function $\Gamma$ displayed in Fig. 6a shows that the matrix $\hat{\mathbf{T}}_{\text{out}}$ can be decomposed as a spatially-invariant component $A_{\text{out}}$ and a short-range correlated component $S_{\text{out}}$. Each component can be separated by performing a singular value decomposition of $\hat{\mathbf{T}}_{\text{out}}$, such that

$$ \hat{\mathbf{T}}_{\text{out}} = \sum_{p=1}^{N} s_p \mathbf{U}_p \mathbf{V}_p^{\dagger} \quad (20) $$
where $s_p$ are the positive and real singular values of $\mathcal{T}_{\text{out}}$ sorted in decreasing order, $\mathbf{U}_p$ and $\mathbf{V}_p$ are unitary matrices whose columns correspond to the singular vectors of $\mathcal{T}_{\text{out}}$ in the pupil and focal planes, respectively. The first eigenspace of $\mathcal{T}_{\text{out}}$ provides its spatially-invariant aberrated component: $\mathbf{A}_{\text{out}} = s_1 \mathbf{U}_1 \mathbf{V}_1^\dagger$. The higher rank eigenstates provide the forward multiple scattering component $\mathbf{S}_{\text{out}}$. Lines or columns of the associated correlation matrix $\mathbf{C}_S = \mathbf{S}_{\text{out}} \times \mathbf{S}_{\text{out}}^\dagger$ provides the isoplanatic patches displayed in Fig. 6e.
where $s_p$ are the positive and real singular values of $\hat{T}_{\text{out}}$ sorted in decreasing order, $U_p$ and $V_p$ are unitary matrices whose columns correspond to the input and output singular vectors of $\hat{T}_{\text{out}}$. The first eigenspace of $\hat{T}_{\text{out}}$ provides its spatially-invariant aberrated component: $\mathbf{A}_{\text{out}} = s_1 \mathbf{U}_1 \mathbf{V}_1^{\dagger}$. The higher rank eigenstates provide the forward multiple scattering component $\mathbf{S}_{\text{out}}$. Lines or columns of the associated correlation matrix $\mathbf{C}_S = \hat{\mathbf{S}}_{\text{out}} \times \hat{\mathbf{S}}_{\text{out}}^{\dagger}$ provides the isoplanatic patches displayed in Fig. 6e.

**Data availability.** Optical data used in this manuscript have been deposited at Zenodo (https://zenodo.org/record/7665117).

**Code availability.** Codes used to post-process the optical data within this paper are available from the corresponding author.

**Acknowledgments.** The authors wish to thank A. Badon for initial discussions about the experimental set up, K. Irsch for providing the corneal sample and A. Le Ber for providing the iterative phase reversal algorithm.

**Funding Information.** 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 nos. 610110 and 819261, HELMHOLTZ* and REMINISCENCE projects, respectively). This project has also received funding from Labex WIFI (Laboratory of Excellence within the French Program Investments for the Future; ANR-10-LABX-24 and ANR-10-IDEX-0001-02 PSL*).

**Author Contributions.** A.A. initiated and supervised the project. A.C.B., V.B. and A.A. designed the experimental setup. U.N., V.B. and P.B. built the experimental set up. U.N. and V.B. developed the post-processing tools. U.N. performed the corneal imaging experiment. U.N. and A.A. analyzed the experimental results. V.B. and A.A. performed the theoretical study. A.A. and U.N. prepared the manuscript. U.N., V.B., P.B., M.F., A.C.B., and A.A. discussed the results and contributed to finalizing the manuscript.

**Competing interests.** A.A., M.F., A.C.B. and V.B. are inventors on a patent related to
**Data availability.** The optical data generated in this study are available at Zenodo⁵⁴ (https://zenodo.org/record/7665117).

**Code availability.** Codes used to post-process the optical data within this paper are available from the corresponding author.

**Acknowledgments.** The authors wish to thank A. Badon for initial discussions about the experimental set up, K. Irsch for providing the corneal sample and A. Le Ber for providing the iterative phase reversal algorithm.

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 nos. 610110 and 819261, HELMHOLTZ* and REMINISCENCE projects, MF and AA, respectively). This project has also received funding from Labex WIFI (Laboratory of Excellence within the French Program Investments for the Future; ANR-10-LABX-24 and ANR-10-IDEX-0001-02 PSL*, MF).

**Author Contributions Statement.** A.A. initiated and supervised the project. C.B., V.B. and A.A. designed the experimental setup. U.N., V.B. and P.B. built the experimental set up. U.N. and V.B. developed the post-processing tools. U.N. performed the corneal imaging experiment. U.N. and A.A. analyzed the experimental results. V.B. and A.A. performed the theoretical study. A.A. and U.N. prepared the manuscript. U.N., V.B., P.B., M.F., C.B., and A.A. discussed the results and contributed to finalizing the manuscript.

**Competing interests.** A.A., M.F., C.B. and V.B. are inventors on a patent related to this work held by CNRS (no. US11408723B2, published August 2022). All authors declare that they have no other competing interests.
this work held by CNRS (no. US11408723B2, published August 2022). All authors declare
that they have no other competing interests.
**Supplementary Information**

This document provides further information on: (i) the experimental set up; (ii) the theoretical expression of the de-scanned matrix; (iii) the measurement of the scattering mean free path; (iv) the theoretical expression of the correlation matrix; (v) the estimation of the transmission matrix; (vi) the contrast enhancement provided by reflection matrix imaging.
Supplementary Information

This document provides further information on: (i) the experimental set up; (ii) the
measurement of the scattering mean free path; (iii) the theoretical expression of the de-
scanned matrix; (iv) the relation between the de-scanned and focused reflection matrices; (v)
the distortion matrix; (vi) its covariance matrix; (vii) iterative time reversal; (viii) iterative
phase reversal; (ix) the bias exhibited by the corresponding T-matrix estimator; (x) a
numerical validation of the iterative phase reversal process; (xi) the multi-scale compensation
of wave distortions; (xii) the convergence of this process and its monitoring by the memory
effect; (xiii) the contrast enhancement provided by RMI; (xiv) the discrepancy observed
between the input and output T-matrices; (xv) the aberration and scattering components of
the T-matrix; (xvi) the comparison between RMI and adaptive optics; (xvii) the comparison
between iterative phase reversal and iterative time reversal; (xviii) the comparison between
the RMI and CLASS approaches.

S1. DETAILED EXPERIMENTAL SET-UP

The full experimental set up is displayed in Fig. S1. The setup is divided into two building
blocks, labelled (a) and (b). The first component is a Michelson interferometer [Fig. S1a].
The light source is a broadband LED (Thorlabs M850LP1, λ₀ = 850 nm, Δλ = 35 nm), is
placed in a plan conjugated with the focal plane within the sample, so as to illuminate the
focal plane with the image of the source. The source's illumination pattern is not uniform
and is smaller than the maximal extension of the field of view allowed by the microscope
objectives' numerical aperture. As a result, the sample's intensity in the focal plane is
modulated by the image of the source. In order to get a uniform illumination of the whole
field of view, we set up a pseudo-Koehler illumination apparatus: An aspheric lens and a
diaphragm are placed right in front of the source, such that the incident beam is collimated
in the diaphragm plane. This plane is considered the source plane, and is conjugated to the
sample plane by (L1). This way, the image of the source is defocused in the sample plane.
This ensures an incoherent³⁰, yet uniform, illumination of the field of view.

The incident light is collimated using a converging lens (L1) with a focal length
$f_1$ = 150 mm. The beam transmitted through this lens (L1) is linearly polarized at
S1. DETAILED EXPERIMENTAL SET-UP

The full experimental set up is displayed in Supplementary Fig. S1. The setup is divided into two building blocks, labelled (a) and (b). The first component is a Michelson interferometer [Supplementary Fig. S1a]. The light source is a broadband LED (Thorlabs M850LP1, $\lambda_o$ = 850 nm, $\Delta\lambda$ = 35 nm), is placed in a plan conjugated with the focal plane within the sample, so as to illuminate the focal plane with the image of the source. The source's illumination pattern is not uniform and is smaller than the maximal extension of the field of view allowed by the microscope objectives' numerical aperture. As a result, the sample's intensity in the focal plane is modulated by the image of the source. In order to get a uniform illumination of the whole field of view, we set up a pseudo-Koehler illumination apparatus: An aspheric lens and a diaphragm are placed right in front of the source, such that the incident beam is collimated in the diaphragm plane. This plane is considered the source plane, and is conjugated to the sample plane by (L1). This way, the image of the source is defocused in the sample plane. This ensures an incoherent$^{30}$, yet uniform, illumination of the field of view.

The incident light is collimated using a converging lens (L1) with a focal length $f_1$ = 150 mm. The beam transmitted through this lens (L1) is linearly polarized at 45° by a polarizer (P1) so that it is then equally reflected (sample arm) and transmitted (reference arm) by the polarized beam splitter (PBS1).

The sample beam reflected by (PBS1) is horizontally polarized. It propagates through a quarter-wave plate (QWP1), is reflected by a plane mirror (M1), whose normal axis lies along the optical axis and that is mounted on a piezoelectric actuator (PZT). The reflected beam passes again through the quarter-wave plate (QWP1). This sequence induces a polarization rotation by 90° of the reflected beam with respect to the incident beam in the sample arm. The reflected wave can be then transmitted through the beam splitter (PBS1) with a vertical polarization and finally focused in a secondary source plane conjugated with the source plane by means of the lens (L2) of focal length $f_2$ = 125 mm.

The reference beam, vertically polarized at the exit of the polarizer (P1), is transmitted by the beam splitter (PBS1), propagates through a quarter-wave plate (QWP2), is reflected by a set of two galvanometric scan mirrors, and then by the reference mirror (M2). The set of scan mirrors enables a 2D rotation of the incident wave-field by angles $\Theta = (\Theta_x, \Theta_y)$ with
45° by a polarizer (P1) so that it is then equally reflected (sample arm) and transmitted (reference arm) by the polarized beam splitter (PBS1).

The sample beam reflected by (PBS1) is horizontally polarized. It propagates through a quarter-wave plate (QWP1), is reflected by a plane mirror (M1), whose normal axis lies along the optical axis and that is mounted on a piezoelectric actuator (PZT). The reflected beam passes again through the quarter-wave plate (QWP1). This sequence induces a polarization rotation by 90° of the reflected beam with respect to the incident beam in the sample arm. The reflected wave can be then transmitted through the beam splitter (PBS1) with a vertical polarization and finally focused in a secondary source plane conjugated with the source plane by means of the lens (L2) of focal length $f_2 = 125$ mm.

The reference beam, vertically polarized at the exit of the polarizer (P1), is transmitted by the beam splitter (PBS1), propagates through a quarter-wave plate (QWP2), is reflected by a set of two galvanometric scan mirrors, and then by the reference mirror (M2). The set of scan mirrors enables a 2D rotation of the incident wave-field by angles $\Theta = (\Theta_x, \Theta_y)$ with respect to the optical axis. After reflection on the reference mirror (M2) and on the scan mirrors back again, the reflected beam propagates again through (QWP2). This round trip through (QWP2) enables a 90° rotation of the polarization: Whereas the incident light is V-polarized when it enters the reference arm, it is H-polarized when exiting it (see Fig. S1a). Therefore, the reference beam is reflected by the beam splitter (PBS1) before being focused by the lens (L2) in the secondary source plane.

Finally, in the secondary source plane, the wave-field is made of two images of the incident light orthogonally polarized and translated with respect to each other by a relative position $\Delta\rho = (\Delta x, \Delta y)$. This lateral shift is dictated by the tilt $\Theta = (\Theta_x, \Theta_y)$ of the reference beam: $\Delta x = 2f_2 \tan \Theta_x$ and $\Delta y = 2f_2 \tan \Theta_y$. The factor 2 results from the double reflection on each scan mirror, due to the reflection on the (M2) mirror. Note also that the optical path difference between the two arms is set to zero by equalizing the length of sample and reference arms for $\Delta\theta = 0$.

After the Michelson interferometer, the two orthogonally polarized twin beams enter a Michelson interferometer with two identical microscope objectives in both arms (a configuration known as a Linnik interferometer) [Fig. S1b]. They are again collimated by a lens (L3) of focal length $f_3=200$ mm. The two lenses (L2) and (L3) thus constitute a 4f system which compensates the effects of diffraction between the two interferometers.
respect to the optical axis. After reflection on the reference mirror (M2) and on the scan
mirrors back again, the reflected beam propagates again through (QWP2). This round trip
through (QWP2) enables a 90° rotation of the polarization: Whereas the incident light is
V-polarized when it enters the reference arm, it is H-polarized when exiting it (see Fig. S1a).
Therefore, the reference beam is reflected by the beam splitter (PBS1) before being focused
by the lens (L2) in the secondary source plane.

Finally, in the secondary source plane, the wave-field is made of two images of the incident
light orthogonally polarized and translated with respect to each other by a relative position
$\Delta\rho = (\Delta x, \Delta y)$. This lateral shift is dictated by the tilt $\Theta = (\Theta_x, \Theta_y)$ of the reference
beam: $\Delta x = 2f_2 \tan \Theta_x$ and $\Delta y = 2f_2 \tan \Theta_y$. The factor 2 results from the double reflection
on each scan mirror, due to the reflection on the (M2) mirror. Note also that the optical
path difference between the two arms is set to zero by equalizing the length of sample and
reference arms for $\Delta\theta = 0$.

After the Michelson interferometer, the two orthogonally polarized twin beams enter a
Michelson interferometer with two identical microscope objectives in both arms (a configura-
tion known as a Linnik interferometer) [Supplementary Fig. S1b]. They are again collimated
by a lens (L3) of focal length $f_3=200$ mm. The two lenses (L2) and (L3) thus constitute a
4f system which compensates the effects of diffraction between the two interferometers.

The vertically polarized light (sample beam) is transmitted by a polarized beam splitter
cube (PBS2), propagates through a quarter-wave plate (QWP4) before being focused in
the focal plane of an immersion microscope objective (MO2, Nikon, 60×, NA=1.0). The
light reflected by the sample is then collected by (MO2) and propagates again through the
quarter-wave plate (QWP4). Because single scattering tends to preserve polarization, the
corresponding wave-field undergoes a 90° polarization rotation and gets reflected by the
beam splitter (PBS2) before being focused in the plane of the camera using the converging
lens (L4) of focal length $f_4=200$ mm. The combination of this lens (L4) with the microscope
objective (MO1) entails a magnification $M_4$ of 60.

Regarding the horizontally-polarized beam at the exit of the lens (L3), it is reflected
by the beam splitter (PBS2), passes through the quarter-wave plate (QWP3) before being
focused by the microscope objective (MO1) identical to (MO2). The light is then reflected by
the reference mirror (M3) placed in the focal plane of (MO2) before being collected again by
the same microscope objective (MO2). The reflected light comes through the quarter-wave
The vertically polarized light (sample beam) is transmitted by a polarized beam splitter cube (PBS2), propagates through a quarter-wave plate (QWP4) before being focused in the focal plane of an immersion microscope objective (MO2, Nikon, 60×, NA=1.0). The light reflected by the sample is then collected by (MO2) and propagates again through the quarter-wave plate (QWP4). Because single scattering tends to preserve polarization, the corresponding wave-field undergoes a 90° polarization rotation and gets reflected by the beam splitter (PBS2) before being focused in the plane of the camera using the converging lens (L4) of focal length $f_4$=200 mm. The combination of this lens (L4) with the microscope objective (MO1) entails a magnification $M_4$ of 60.

Regarding the horizontally-polarized beam at the exit of the lens (L3), it is reflected by the beam splitter (PBS2), passes through the quarter-wave plate (QWP3) before being focused by the microscope objective (MO1) identical to (MO2). The light is then reflected by the reference mirror (M3) placed in the focal plane of (MO2) before being collected again by the same microscope objective (MO2). The reflected light comes through the quarter-wave plate (QWP3). As in the other arm, the polarization of the reflected beam exhibits a 90° rotation of its polarization. The beam is now vertically polarized and transmitted by the beam splitter (PBS2), before being focused on the camera with the lens (L4).

The detection scheme consists in recording the interferogram between the sample and reference beams based on their orthogonal polarisations projected through a 45-degree polarizer (P2). This experimental configuration allows an enhancement of single scattering and forward multiple scattering with respect to diffuse light in the sample arm. Indeed, the former components roughly exhibit the same time-of-flight and polarization as reference light while the latter one is characterized by a fully randomized polarization and a longer time-of-flight distribution. This filtering of diffuse light is deliberate since the post-processing method described in the accompanying paper addresses the forward multiple scattering contribution and not the randomly-scattered diffuse light.

The CMOS camera (Adimec Quartz 2A-750, 2Mpx) records the interferogram between sample and reference beams with a spatial sampling equal to $\delta_0$ =230 nm given the magnification $M_4$. The volume of the sample from which photons can interfere with the reference beam is called the “coherence volume”. Its position is dictated by the optical path difference between the reference and sample arms. Its thickness is inversely proportional to the light
plate (QWP3). As in the other arm, the polarization of the reflected beam exhibits a 90° rotation of its polarization. The beam is now vertically polarized and transmitted by the beam splitter (PBS2), before being focused on the camera with the lens (L4).

The detection scheme consists in recording the interferogram between the sample and reference beams based on their orthogonal polarisations projected through a 45-degree polarizer (P2). This experimental configuration allows an enhancement of single scattering and forward multiple scattering with respect to diffuse light in the sample arm. Indeed, the former components roughly exhibit the same time-of-flight and polarization as reference light while the latter one is characterized by a fully randomized polarization and a longer time-of-flight distribution. This filtering of diffuse light is deliberate since the post-processing method described in the accompanying paper addresses the forward multiple scattering contribution and not the randomly-scattered diffuse light.

The CMOS camera (Adimec Quartz 2A-750, 2Mpx) records the interferogram between sample and reference beams with a spatial sampling equal to $\delta_0 = 230$ nm given the magnification $M_4$. The volume of the sample from which photons can interfere with the reference beam is called the “coherence volume”. Its position is dictated by the optical path difference between the reference and sample arms. Its thickness is inversely proportional to the light spectrum bandwidth$^{55}$:

$$ \delta z_t = \frac{2 \ln 2}{n\pi} \left( \frac{\lambda_0^2}{\Delta\lambda} \right) \qquad (S1) $$

with $\lambda_0$ the central wavelength of the light source and $\Delta\lambda$ its spectral bandwidth. In the present case, $\delta z_t \sim 10$ µm. A critical tuning of the experimental set up consists in adjusting the coherence volume with the focal plane of the microscope objective. In a volumetric sample, whose refractive index differs from that of water, the coherence volume no longer coincides with the focusing plane. This focusing defect accumulates with the transverse aberrations generated by the heterogeneities of the medium. However, it is possible to compensate for it by a fine tuning of the length of the reference arm.

The experimental procedure then consists in recording the de-scanned reflection matrix $R_{in}(z)$ at each depth $z$ of the sample. This latter parameter is swept by means of a motorized axial displacement of the sample carrier. The scan of the relative position $\Delta\rho$ between the incident wave-fields in the sample and reference arms is controlled by the tilt imposed by the galvanometer (M2). For each couple $(\Delta\rho, z)$, the CCD camera conjugated with the MO
spectrum bandwidth⁵⁴:

$$ \delta z_t = \frac{2 \ln 2}{n\pi} \left( \frac{\lambda^2}{\Delta\lambda} \right) \qquad (S1) $$

with $\lambda$ the central wavelength of the light source and $\Delta\lambda$ its spectral bandwidth. In the present case, $\delta z_t \sim 10$ µm. A critical tuning of the experimental set up consists in adjusting the coherence volume with the focal plane of the microscope objective. In a volumetric sample, whose refractive index differs from that of water, the coherence volume no longer coincides with the focusing plane. This focusing defect accumulates with the transverse aberrations generated by the heterogeneities of the medium. However, it is possible to compensate for it by a fine tuning of the length of the reference arm.

The experimental procedure then consists in recording the de-scanned reflection matrix $\mathbf{R}_{\text{in}}(z)$ at each depth $z$ of the sample. This latter parameter is swept by means of a motorized axial displacement of the sample carrier. The scan of the relative position $\Delta\rho$ between the incident wave-fields in the sample and reference arms is controlled by the tilt imposed by the galvanometer (M2). For each couple $(\Delta\rho, z)$, the CCD camera conjugated with the MO focal plane records the output intensity:

$$ I_{\alpha}(\Delta\rho, \rho_{\text{out}}, z) = \int_{0}^{T} |e^{i\alpha} E_{\text{out}}(\rho_{\text{out}}, t) + E_{\text{out}}^{(\text{ref})*}(\Delta\rho, \rho_{\text{out}}, t)|^2 dt \qquad (S2) $$

with $t$ the absolute time, $\mathbf{r}$ the position vector on the CCD screen, $E_{\text{out}}(\mathbf{r}, \tau)$ the scattered wave field associated with the sample arm, $E_{\text{out}}^{(\text{ref})}(\mathbf{r}, \tau)$ the reference wave field; $T$ the integration time of the CCD camera, and $\alpha$ an additional phase term controlled with a piezoelectric actuator placed on mirror (M1) of the first interferometer [Fig. S1a]. The interference term between the sample and reference beams is extracted from the four intensity patterns (Eq. S2) recorded at $\alpha = 0, \pi/2, 3\pi/2$ and $\pi$ (phase-stepping interferometry) and provides the de-scanned wave-field:

$$ R_{\text{in}}(\Delta\rho, \rho_{\text{out}}, z) = \frac{1}{T} \int_{0}^{T} E_{\text{out}}(\rho_{\text{out}}, t) E_{\text{out}}^{(\text{ref})*}(\Delta\rho, \rho_{\text{out}}, t) dt \qquad (S3) $$

Note that previous studies²⁵,²⁶ reported on the passive measurement of the de-scanned reflection matrix at the surface of a scattering sample. Although the experimental set up presented in those studies shares some similarities with the current set up displayed in Fig. S1 (low-coherence interferometry), there are also several major differences. First, the two arms are illuminated by the light reflected by the sample in Ref.²⁵; hence there is no reference arm. Second, the reflection matrix is measured at the surface of the sample as a function of
focal plane records the output intensity:

$$I_{\alpha}(\Delta\rho, \rho_{\text{out}}, z) = \int_{0}^{T} |e^{i\alpha} E_{\text{out}}(\rho_{\text{out}}, t) + E_{\text{out}}^{(\text{ref})} (\Delta\rho, \rho_{\text{out}}, t)|^{2} dt \quad (\text{S2})$$

with $t$ the absolute time, $\mathbf{r}$ the position vector on the CCD screen, $E_{\text{out}}(\mathbf{r}, \tau)$ the scattered wave field associated with the sample arm, $E_{\text{out}}^{(\text{ref})}(\mathbf{r}, \tau)$ the reference wave field; $T$ the integration time of the CCD camera, and $\alpha$ an additional phase term controlled with a piezoelectric actuator placed on mirror (M1) of the first interferometer [Supplementary Fig. S1a]. The interference term between the sample and reference beams is extracted from the four intensity patterns (Eq. S2) recorded at $\alpha = 0, \pi/2, 3\pi/2$ and $\pi$ (phase-stepping interferometry) and provides the de-scanned wave-field:

$$R_{\text{in}}(\Delta\rho, \rho_{\text{out}}, z) = \frac{1}{T} \int_{0}^{T} E_{\text{out}}(\rho_{\text{out}}, t) E_{\text{out}}^{(\text{ref})} (\Delta\rho, \rho_{\text{out}}, t) dt \quad (\text{S3})$$

Note that previous studies$^{25,26}$ reported on the passive measurement of the de-scanned reflection matrix at the surface of a scattering sample. Although the experimental set up presented in those studies shares some similarities with the current set up displayed in Supplementary Fig. S1 (low-coherence interferometry), there are also several major differences. First, the two arms are illuminated by the light reflected by the sample in Ref.$^{25}$; hence there is no reference arm. Second, the reflection matrix is measured at the surface of the sample as a function of time-of-flight, while the current set up measures a time-gated reflection matrix as a function of depth inside the sample. At last, a much larger integration time is required to record the reflection matrix in Ref.$^{25}$ because of the absence of reference arm.
time-of-flight, while the current set up measures a time-gated reflection matrix as a function
of depth inside the sample. At last, a much larger integration time is required to record the
reflection matrix in Ref.²⁵ because of the absence of reference arm.
<figure><img src="image_7.png" /><figcaption>FIG. S1. De-scanned measurement of the reflection matrix. a. Experimental set up. P: polarizer, L: lens, QWP: quarter-wave plate, M: mirror, PZT: piezo-electric actuator, PBS: polarisation beam splitter, MO: microscope objective. The apparatus is made up of two parts ($\mathbf{a}_1$, $\mathbf{a}_2$). $\mathbf{a}_1$. Michelson interferometer illuminated by an incoherent light source at its input and generating two twin incoherent beams of orthogonal polarization and laterally shifted from each other at its output. The polarised beam splitter (PBS1) separates the impinging light into a reference path (in blue) and a sample path (in red). The scan mirrors ($\mathbf{a}_3$) tilt the reference beam by angles $\Delta\Theta_x$ and $\Delta\Theta_y$ in both transverse directions $x$ and $y$. This allows the scan of the point-spread function along the de-scanned coordinates $\Delta_x$ and $\Delta_y$ in a plane conjugate to the reference mirror plane. $\mathbf{a}_2$. Michelson interferometer with microscope objectives (MO) in both arms (Linnik configuration). Both beams have orthogonal polarizations and each interferometer arm includes a quarter-wave plate (QWP). The output beams are collected by the L4 lens and interfere on the camera after having been projected on a 45°-rotated polarizer (P2). b. Equivalent layout in the case of a coherent measurement. The source plane, the focal plane, and the camera planes are conjugated. Displacing a point source $\rho_{\text{in}}$ in the source plane discretely scans the focal plane inside the sample. The illuminated area is imaged in the camera plane; in an epi-detection configuration.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. S1. De-scanned measurement of the reflection matrix. a. Experimental set up. P: polarizer, L: lens, QWP: quarter-wave plate, M: mirror, PZT: piezo-electric actuator, PBS: polarisation beam splitter, MO: microscope objective. The apparatus is made up of two parts ($\mathbf{a}_1$, $\mathbf{a}_2$). $\mathbf{a}_1$. Michelson interferometer illuminated by an incoherent light source at its input and generating two twin incoherent beams of orthogonal polarization and laterally shifted from each other at its output. The polarised beam splitter (PBS1) separates the impinging light into a reference path (in blue) and a sample path (in red). The scan mirrors ($\mathbf{a}_3$) tilt the reference beam by angles $\Delta\Theta_x$ and $\Delta\Theta_y$ in both transverse directions $x$ and $y$. This allows the scan of the point-spread function along the de-scanned coordinates $\Delta_x$ and $\Delta_y$ in a plane conjugate to the reference mirror plane. $\mathbf{a}_2$. Michelson interferometer with microscope objectives (MO) in both arms (Linnik configuration). Both beams have orthogonal polarizations and each interferometer arm includes a quarter-wave plate (QWP). The output beams are collected by the L4 lens and interfere on the camera after having been projected on a 45°-rotated polarizer (P2). b. Equivalent layout in the case of a coherent measurement. The source plane, the focal plane, and the camera planes are conjugated. Displacing a point source $\rho_{\text{in}}$ in the source plane discretely scans the focal plane inside the sample. The illuminated area is imaged in the camera plane; in an epi-detection configuration.</figcaption></figure>
## S2. THEORETICAL ANALYSIS

### A. Theoretical expression of the de-scanned matrix

In this section, we express theoretically the de-scanned matrix recorded by the experimental set up in Supplementary Figs. S1a,b. To that aim, we will rely on the simple Fourier optics model proposed in a recent paper²⁹ to describe the manifestation of aberrations in FFOCT. For the sake of simplicity, this model is scalar. The large numerical aperture imposes that the recorded wave-field is associated with single scattering events taking place in the focal plane of the MO.

The wave field $E_{\text{out}}(\rho_{\text{out}}, z)$ reflected by the sample arm in the camera plane can then be expressed as follows²⁹:

$$E_{\text{out}}(\rho_{\text{out}}, z, \omega) = \int_{\Sigma_0} \iint_{\Sigma_{\rho}} H(\rho_{\text{out}}, \rho_s, z) \gamma(\rho_s, z) H(\rho_s, \rho_0, z) E_0(\rho_0, \omega) d\rho_s d\rho_0. \quad (S4)$$

$E_0(\rho_0, \omega)$ is the incident wave-field in the secondary source plane $\Sigma_0$ at frequency $\omega$. Light propagation between $\Sigma_0$ and the focal plane $\Sigma_\rho$ is described by the impulse response $H(\rho_0, \rho_s, z)$ between a point in the secondary source plane at transverse coordinate $\rho_0$ and a point at transverse coordinate $\rho_s$ in the focal plane and at depth $z$ inside the sample. It accounts for sample-induced aberrations. $\gamma(\rho_s, z)$ represents the sample reflectivity at depth $z$. By spatial reciprocity, the propagation of the reflected wave-field from the sample to the detector plane is also modelled by the impulse response $H(\rho_s, \rho_{\text{out}})$. The relatively narrow bandwidth ($\Delta\lambda \ll \lambda$) of the light source and the use of achromatic optical elements (lens, beam splitter, quarter wave plate) allow us to neglect the dependence of $H$ on frequency $\omega$.

Replacing $\gamma(\rho_s, z)$ by a uniform reflectivity in Eq. S4 and taking into account the lateral shift of the reference wave-field induced by the galvanometer M2 [Supplementary Fig. S1] leads to the following previous expression for $E_{\text{out}}^{(\text{ref})}(\rho_{\text{out}}, z)^{29}$:

$$E_{\text{out}}^{(\text{ref})}(\rho_{\text{out}}, \Delta\rho, z) = \iint_{\Sigma_0} H_{\text{ref}}(\rho_{\text{out}} - \rho_0) E_0(\rho_0 + \Delta\rho) d\rho_0. \quad (S5)$$

where $H_{\text{ref}}$ is the impulse response associated with the reference arm (way and return path) that we assume as spatially-invariant [$H_{\text{ref}}(\rho_{\text{out}}, \rho_0) = H_{\text{ref}}(\rho_{\text{out}} - \rho_0)$].

The de-scanned wave-field is obtained by extracting the interference term between the reflected wave-fields coming from the sample and reference arms:

$$R_{\text{in}}(\rho_{\text{out}}, \Delta\rho, z) = \langle E_{\text{out}}(\rho_{\text{out}}, \omega) E_{\text{out}}^{(\text{ref})*}(\rho_{\text{out}}, \omega) \rangle \quad (S6)$$
## S2. MEASURING THE SCATTERING MEAN FREE PATH

<figure><img src="image_9.png" /><figcaption>FIG. S2. Confocal scattering ratio ln $\beta$ versus depth (blue dots) fitted with Eq. S4 (red line).</figcaption></figure>

In a previous work⁵⁵, the scattering mean free path $\ell_s$ in the cornea was measured by investigating the depth evolution of the confocal intensity. Indeed, in the single scattering regime, under the paraxial approximation and for an homogeneous reflectivity, the time-gated confocal intensity is supposed to decrease as $\exp(-2z/\ell_s)$ if we neglect absorption losses⁸,⁵⁶.

Unfortunately, here, the cornea is not healthy but oedematous. The depth evolution of the confocal intensity in the stroma is thus strongly impacted by multiple scattering and cannot be used for a measurement of $\ell_s$. Moreover, in the epithelium, the different layers of cell make the cornea reflectivity too heterogeneous to provide an exponential decrease of the confocal intensity.

Recently, an alternative strategy has been proposed in presence of multiple scattering. It consists in investigating the depth evolution of the ratio $\beta$ between the confocal surintensity and the total intensity³⁴. For a medium statistically homogeneous in terms of disorder, numerical simulations have shown empirically that this quantity scales as ⁵⁷:

$$ \beta(z) \sim \exp[-2z/(3\ell_s)] \qquad (S4) $$
Assuming a spatially-incoherent incident wave-field $[\langle E_0(\rho_0)E_0^*(\rho'_0) \rangle = I_0\delta(\rho_0 - \rho'_0)]$ and injecting Eqs. S4 and S5 into the last equation leads to the following expression for the coefficients of $\mathbf{R}_{\text{in}}$:

$$R_{\text{in}}(\rho_{\text{out}}, \Delta\rho, z) = I_0 \iint H_{\text{out}}(\rho_{\text{out}}, \rho, z)\gamma(\rho_s, z)H_{\text{in}}(\rho_s, \rho_{\text{out}} + \Delta\rho, z)d\rho \quad (\text{S7})$$

with

$$H_{\text{out}} \equiv H \quad \text{and} \quad H_{\text{in}} \equiv H \odot H_{\text{ref}}^* \quad (\text{S8})$$

The symbol $\odot$ stands for the convolution product over the variable $\rho_{\text{in}}$. The last equation means that: (i) the output focusing matrix, $\mathbf{H}_{\text{out}}$, and the associated $\mathbf{T}$-matrix, $\mathbf{T}_{\text{out}} = \mathbf{T}_0 \times \mathbf{H}_{\text{out}}$, only grasp the sample-induced aberrations; (ii) the input focusing matrix, $\mathbf{H}_{\text{in}}$, and the associated $\mathbf{T}$-matrix, $\mathbf{T}_{\text{in}} = \mathbf{T}_0 \times \mathbf{H}_{\text{in}}$, also contain the aberrations undergone by the incident and reflected reference beams (Supplementary Section S4 A).

For $\Delta\rho = 0$ (conventional FFOCT set up), the recorded wave-field (Eq. S7) is equivalent to a time-gated confocal image$^{29}$. It can actually be expressed as the convolution between the sample reflectivity $\gamma$ and the confocal PSF $H_{\text{in}} \times H_{\text{out}}$:

$$R_{\text{in}}(\rho_{\text{out}}, \Delta\rho = 0, z) = I_0 \iint H_{\text{out}}(\rho_{\text{out}}, \rho_s, z)H_{\text{in}}(\rho_s, \rho_{\text{out}}, z)\gamma(\rho_s, z)d\rho_s \quad (\text{S9})$$

On the one hand, the confocal nature of the recorded wave-field implies a transverse resolution $\delta\rho_0 \sim \lambda/4NA$. On the other hand, the axial resolution is either controlled by the thickness $\delta z_t$ of the coherence volume or the depth-of-field $\delta z_0$ of the microscope objective: $\delta z_0 = n\lambda/NA^2$. In the present case, $\delta z_0 \sim 1 \mu\text{m} < \delta z_t \sim 10 \mu\text{m}$. The axial resolution is thus given by the depth-of-field. $\delta\rho_0$ and $\delta z_0$ thus dictate the values of the transverse and axial sampling of the de-scanned matrix $\mathbf{R}_{\text{in}}$ in our experiment.
In the present case, this confocal ratio $\beta$ has been measured as follows:

$$ \hat{\beta}(z) = 1 - \frac{\min_{\Delta\rho}\{I(\Delta\rho, z)\}}{\max_{\Delta\rho}\{I(\Delta\rho, z)\}} \qquad (S5) $$

This estimator $\hat{\beta}$ relies on the fact that the multiple scattering component of the RPSF exhibits a flat background such that it can be estimated with the minimum of $I(\Delta\rho, z)$. This hypothesis is wrong at shallow depth since the diffuse halo grows as $\sqrt{Dt}$. Nevertheless, beyond $\ell_s$ or so (here 100 µm), the multiple scattering background can be considered as flat as illustrated by Fig. 3d of the accompanying paper.

Figure S2 displays the depth evolution of the estimator $\hat{\beta}(z)$. It exhibits an exponential decay in the stroma beyond $z = 100$ µm. The decay rate decreases beyond $z = 170$ µm because our estimator of $\beta(z)$ starts to be impacted by the experimental noise [see Fig. 3d of the accompanying paper]. Therefore, the fit of $\hat{\beta}(z)$ with Eq. S4 is performed from $z = 100$ to $z = 170$ µm. We find $\ell_s \sim 35$ µm.
## B. Relation between the de-scanned matrix and the focused reflection matrix

In this section, we investigate to which extent the de-scanned matrix recorded by the experimental set up in Supplementary Figs. S1a,b can be considered equivalent to the focused reflection matrix that would be recorded by the fictitious coherent set up displayed in Supplementary Fig. S1c.

The coefficients of a focused reflection matrix recorded by the fictitious coherent set up displayed in Fig. S1 can be expressed as:

$$R(\rho_{\text{out}}, \rho_{\text{in}}, z) = I_0 \iiint H(\rho_{\text{out}}, \rho_s, z)\gamma(\rho_s, z)H(\rho_s, \rho_{\text{in}}, z)d\rho \quad (\text{S10})$$

A strict equality between Eqs. S7 and S10 is only obtained if $H_{\text{in}} \equiv H$. This condition is fulfilled only for a perfect reference arm: $T_{\text{ref}} \equiv \mathbf{1}_{k<NA}$ and $H \odot H_{\text{ref}}^* \equiv H$. In theory, the incoherent set up of Fig. S1a is thus equivalent to the fictitious coherent set up of Fig. S1b.

$$R_{\text{in}}(\Delta\rho, \rho_{\text{out}}) = R(\rho_{\text{out}} + \Delta\rho, \rho_{\text{out}}, z) \quad (\text{S11})$$

In reality, the reference arm always exhibits aberrations such as a slight defocus of the reference mirror M3 in Fig. S1b or a slight defocus of the reference beam in the secondary source plane at the output of first interferometer.

## C. The distortion matrix

The distortion matrix is related to the de-scanned matrix by a simple Fourier transform:

$$D_{\text{in}}(z) = T_0 \times R_{\text{in}}(z), \quad (\text{S12})$$

or in terms of matrix coefficients,

$$D(\mathbf{u}_{\text{in}}, \rho_{\text{out}}, z) = \sum_{\Delta\rho} R_{\text{in}}(\Delta\rho, \rho_{\text{out}}, z) \exp \left( -i \frac{2\pi}{\lambda f} \mathbf{u}_{\text{in}} \cdot \Delta\rho \right). \quad (\text{S13})$$

Injecting Eq. S7 into the last equation yields

$$D(\mathbf{u}_{\text{in}}, \rho_{\text{out}}, z) = \sum_{\Delta\rho} I_0 \iint H_{\text{out}}(\rho_{\text{out}}, \rho_s, z)\gamma(\rho_s, z)H_{\text{in}}(\rho_s, \rho_{\text{out}}+\Delta\rho, z) \exp \left( -i \frac{2\pi}{\lambda f} \mathbf{u}_{\text{in}} \cdot \Delta\rho \right) d\rho_s \quad (\text{S14})$$

In a previous paper¹⁹, we showed that a singular value decomposition of $\mathbf{D}$ enables to decompose the field-of-view into isoplanatic modes and extract the associated aberration
S3. THEORETICAL EXPRESSION OF THE DE-SCANNED MATRIX

In this section, we express theoretically the de-scanned matrix recorded by the experimental set up in Figs. S1a,b. To that aim, we will rely on the simple Fourier optics model proposed in a recent paper²⁹ to describe the manifestation of aberrations in FFOCT. For the sake of simplicity, this model is scalar. The large numerical aperture imposes that the recorded wave-field is associated with single scattering events taking place in the focal plane of the MO.

The wave field $E_{\text{out}}(\boldsymbol{\rho}_{\text{out}}, z)$ reflected by the sample arm in the camera plane can then be expressed as follows²⁹:

$$E_{\text{out}}(\boldsymbol{\rho}_{\text{out}}, z, \omega) = \int_{\Sigma_0} \iint_{\Sigma_{\rho}} H(\boldsymbol{\rho}_{\text{out}}, \boldsymbol{\rho}_s, z) \gamma(\boldsymbol{\rho}_s, z) H(\boldsymbol{\rho}_s, \boldsymbol{\rho}_0, z) E_0(\boldsymbol{\rho}_0, \omega) d\boldsymbol{\rho}_s d\boldsymbol{\rho}_0. \quad (S6)$$

$E_0(\boldsymbol{\rho}_0, \omega)$ is the incident wave-field in the secondary source plane $\Sigma_0$ at frequency $\omega$. Light propagation between $\Sigma_0$ and the focal plane $\Sigma_\rho$ is described by the impulse response $H(\boldsymbol{\rho}_0, \boldsymbol{\rho}_s, z)$ between a point in the secondary source plane at transverse coordinate $\boldsymbol{\rho}_0$ and a point at transverse coordinate $\boldsymbol{\rho}_s$ in the focal plane and at depth $z$ inside the sample. It accounts for sample-induced aberrations. $\gamma(\boldsymbol{\rho}_s, z)$ represents the sample reflectivity at depth $z$. By spatial reciprocity, the propagation of the reflected wave-field from the sample to the detector plane is also modelled by the impulse response $H(\boldsymbol{\rho}_s, \boldsymbol{\rho}_{\text{out}})$. The relatively narrow bandwidth ($\Delta\lambda \ll \lambda$) of the light source and the use of achromatic optical elements (lens, beam splitter, quarter wave plate) allow us to neglect the dependence of $H$ on frequency $\omega$.

Replacing $\gamma(\boldsymbol{\rho}_s, z)$ by a uniform reflectivity in Eq. S6 and taking into account the lateral shift of the reference wave-field induced by the galvanometer M2 [Fig. S1] leads to the following previous expression for $E_{\text{out}}^{(\text{ref})}(\boldsymbol{\rho}_{\text{out}}, z)^{29}$:

$$E_{\text{out}}^{(\text{ref})}(\boldsymbol{\rho}_{\text{out}}, \Delta\boldsymbol{\rho}, z) = \iint_{\Sigma_0} H_{\text{ref}}(\boldsymbol{\rho}_{\text{out}} - \boldsymbol{\rho}_0) E_0(\boldsymbol{\rho}_0 + \Delta\boldsymbol{\rho}) d\boldsymbol{\rho}_0. \quad (S7)$$

where $H_{\text{ref}}$ is the impulse response associated with the reference arm (way and return path) that we assume as spatially-invariant $[H_{\text{ref}}(\boldsymbol{\rho}_{\text{out}}, \boldsymbol{\rho}_0) = H_{\text{ref}}(\boldsymbol{\rho}_{\text{out}} - \boldsymbol{\rho}_0)]$.

The de-scanned wave-field is obtained by extracting the interference term between the reflected wave-fields coming from the sample and reference arms:

$$R_{\text{in}}(\boldsymbol{\rho}_{\text{out}}, \Delta\boldsymbol{\rho}, z) = \langle E_{\text{out}}(\boldsymbol{\rho}_{\text{out}}, \omega) E_{\text{out}}^{(\text{ref})*}(\boldsymbol{\rho}_{\text{out}}, \omega) \rangle \quad (S8)$$
phase laws. However, this demonstration was based on the condition that the D-matrix is dominated by its correlations in the focal plane. This is the case for a specular reflector such as a resolution target or a medium of continuous reflectivity but no longer valid for a random distribution of heterogeneities like in the opaque cornea under study. In the accompanying paper, we propose a more general solution to overcome aberrations and scattering in optical microscopy: An iterative multi-scale analysis of wave distortions.

To that aim, the field-of-detection should be subdivided into overlapping regions that are defined by their central midpoint $\mathbf{r}_p = (\boldsymbol{\rho}_p, z_p)$ and their spatial extension $L$. All of the distorted components associated with focusing points $\boldsymbol{\rho}_{\text{out}}$ located within each region are extracted and stored in a local distortion matrix $\mathbf{D}'_{\text{in}}(\mathbf{r}_p)$:

$$D'(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, \mathbf{r}_p) = D(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, z_p) W_L(\boldsymbol{\rho}_{\text{out}} - \boldsymbol{\rho}_p), \qquad (\text{S15})$$

where $W_L(x, y) = 1$ for $|x| < L$ and $|y| < L$, and zero otherwise.

At this stage, a local isoplanatic assumption shall be made over each region of size $L$. This hypothesis implies that the PSFs $H_{\text{in/out}}$ are invariant by translation in each region. This leads us to define local spatially-invariant PSFs $H_{\text{in/out}}^{(l)}$ around each central midpoint $\mathbf{r}_p$ such that:

$$H_{\text{in/out}}(\boldsymbol{\rho}_s, \boldsymbol{\rho}_{\text{in/out}}, z_p) = H_{\text{in/out}}^{(l)}(\boldsymbol{\rho}_s - \boldsymbol{\rho}_{\text{in/out}}, \mathbf{r}_p). \qquad (\text{S16})$$

Under this assumption, Eq. S14 can be rewritten as follows:

$$D'(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, \mathbf{r}_p) = \underbrace{\overbrace{\text{transmittance}}^{\mathcal{T}_{\text{in}}(\mathbf{u}_{\text{in}}, \mathbf{r}_p)} \iint \gamma(\boldsymbol{\rho}_s + \boldsymbol{\rho}_{\text{out}}, z) H_{\text{out}}^{(l)}(\boldsymbol{\rho}_s, \mathbf{r}_p) \exp\left(-i \frac{2\pi}{\lambda f} \mathbf{u}_{\text{in}} \cdot \boldsymbol{\rho}\right) d\boldsymbol{\rho}_s}_{\text{virtual source}} \qquad (\text{S17})$$

Around each point $\mathbf{r}_p$, the aberrations can be modelled by a local transmittance $\mathcal{T}_{\text{in}}(\mathbf{u}_{\text{in}}, \mathbf{r}_p)$. This transmittance is the Fourier transform of the input PSF $H_{\text{in}}^{(l)}(\boldsymbol{\rho}_s, \mathbf{r}_p)$:

$$\mathcal{T}_{\text{in}}(\mathbf{u}_{\text{in}}, \mathbf{r}_p) = \iint H_{\text{in}}^{(l)}(\boldsymbol{\rho}_s, \mathbf{r}_p) \exp\left(-\frac{2\pi}{\lambda f} \mathbf{u}_{\text{in}} \cdot \boldsymbol{\rho}_s\right) d\boldsymbol{\rho}_s \qquad (\text{S18})$$

The physical meaning of this last equation is the following: Each distorted wave-field corresponds to the diffraction of a virtual source synthesized inside the medium modulated by the pupil transmittance $\mathcal{T}_{\text{in}}(\mathbf{u}_{\text{in}}, \mathbf{r}_p)$ of the sample seen from point $\mathbf{r}_p$. Each virtual source is spatially incoherent due to the random reflectivity of the medium, and its size is governed by the spatial extension of the output focal spot. The idea is now to smartly combine each
Assuming a spatially-incoherent incident wave-field $[\langle E_0(\rho_0)E_0^*(\rho'_0) \rangle = I_0\delta(\rho_0 - \rho'_0)]$ and injecting Eqs. S6 and S7 into the last equation leads to the following expression for the coefficients of $\mathbf{R}_{\text{in}}$:

$$R_{\text{in}}(\rho_{\text{out}}, \Delta\rho, z) = I_0 \iint H_{\text{out}}(\rho_{\text{out}}, \rho, z)\gamma(\rho_s, z)H_{\text{in}}(\rho_s, \rho_{\text{out}} + \Delta\rho, z)d\rho \quad (\text{S9})$$

with

$$H_{\text{out}} \equiv H \quad \text{and} \quad H_{\text{in}} \equiv H \odot H_{\text{ref}}^*. \quad (\text{S10})$$

The symbol $\odot$ stands for the convolution product over the variable $\rho_{\text{in}}$. The last equation means that: (i) the output focusing matrix, $\mathbf{H}_{\text{out}}$, and the associated $\mathbf{T}$-matrix, $\mathbf{T}_{\text{out}} = \mathbf{F} \times \mathbf{H}_{\text{out}}$, only grasp the sample-induced aberrations; (ii) the input focusing matrix, $\mathbf{H}_{\text{in}}$, and the associated $\mathbf{T}$-matrix, $\mathbf{T}_{\text{in}} = \mathbf{F} \times \mathbf{H}_{\text{in}}$, also contain the aberrations undergone by the incident and reflected reference beams (Supplementary Section S14).

For $\Delta\rho = 0$ (conventional FFOCT set up), the recorded wave-field (Eq. S9) is equivalent to a time-gated confocal image$^{29}$. It can actually be expressed as the convolution between the sample reflectivity $\gamma$ and the confocal PSF $H_{\text{in}} \times H_{\text{out}}$:

$$R_{\text{in}}(\rho_{\text{out}}, \Delta\rho = 0, z) = I_0 \iint H_{\text{out}}(\rho_{\text{out}}, \rho_s, z)H_{\text{in}}(\rho_s, \rho_{\text{out}}, z)\gamma(\rho_s, z)d\rho_s. \quad (\text{S11})$$

On the one hand, the confocal nature of the recorded wave-field implies a transverse resolution $\delta\rho_0 \sim \lambda/4NA$. On the other hand, the axial resolution is either controlled by the thickness $\delta z_t$ of the coherence volume or the depth-of-field $\delta z_0$ of the microscope objective: $\delta z_0 = n\lambda/NA^2$. In the present case, $\delta z_0 \sim 1 \mu\text{m} < \delta z_t \sim 10 \mu\text{m}$. The axial resolution is thus given by the depth-of-field. $\delta\rho_0$ and $\delta z_0$ thus dictate the values of the transverse and axial sampling of the de-scanned matrix $\mathbf{R}_{\text{in}}$ in our experiment.
virtual source to generate a coherent guide star and estimate the local transmittance $T_{in}$ independently from the sample reflectivity.

### D. Covariance Matrix of Wave Distortions

To do so, the correlation matrix $C_{in} = D_{in}D_{in}^{\dagger}$ is an excellent tool. Its coefficients write as follows

$$C_{in}(\mathbf{u}_{in}, \mathbf{u}'_{in}, \mathbf{r}_{p}) = N_{W}^{-1} \sum_{\rho_{out}} D'(\mathbf{u}_{in}, \rho_{out}, \mathbf{r}_{p}) D'^{*}(\mathbf{u}'_{in}, \rho_{out}, \mathbf{r}_{p}) \quad (S19)$$

The matrix $C_{in}(\mathbf{r}_{p})$ can be decomposed as the sum of its ensemble average, the covariance matrix $\langle C_{in} \rangle (\mathbf{r}_{p})$, and a perturbation term $\delta C_{in}(\mathbf{r}_{p})$:

$$C_{in}(\mathbf{r}_{p}) = \langle C_{in} \rangle (\mathbf{r}_{p}) + \delta C_{in}(\mathbf{r}_{p}). \quad (S20)$$

The intensity of the perturbation term scales as the inverse of the number $N_W = (L/\delta\rho_0)^2$ of resolution cells in each sub-region$^{40,56}$:

$$\langle |\delta C_{in}(\mathbf{u}, \mathbf{u}', \mathbf{r}_p)|^2 \rangle = \frac{\langle |C_{in}(\mathbf{u}, \mathbf{u}, \mathbf{r}_p)|^2 \rangle}{N_W} \quad (S21)$$

This perturbation term can thus be reduced by increasing the size $L$ of the spatial window $W_L$, but at the cost of a resolution loss.

Under assumptions of local isoplanicity (Eqs. S16 and S17) and random reflectivity,

$$\langle \gamma(\rho_s, z) \gamma^*(\rho'_s, z) \rangle = \langle |\gamma|^2 \rangle \delta(\rho_s - \rho'_s), \quad (S22)$$

with $\delta$, the Dirac distribution, the covariance matrix can be expressed as follows$^{21}$:

$$\langle C_{in} \rangle (\mathbf{r}_p) = [\mathcal{T}_{in}(\mathbf{r}_p) \circ \mathbf{T}_0] \times C_H(\mathbf{r}_p) \times [\mathbf{T}_0 \circ \mathcal{T}_{in}(\mathbf{r}_p)]^{\dagger}, \quad (S23)$$

or, in terms of matrix coefficients,

$$\begin{aligned} \langle C_{in} \rangle (\mathbf{u}_{in}, \mathbf{u}'_{in}, \mathbf{r}_{p}) &= \mathcal{T}_{in}(\mathbf{u}_{in}, \mathbf{r}_{p}) \mathcal{T}_{in}^{*}(\mathbf{u}'_{in}, \mathbf{r}_{p}) \iint d\rho_s |H_{out}^{(l)}(\rho_s, \mathbf{r}_{p})|^2 \exp \left[ -i \frac{2\pi}{\lambda f} (\mathbf{u}_{in} - \mathbf{u}'_{in}) \cdot \mathbf{\rho}_s \right] \\ &= \mathcal{T}_{in}(\mathbf{u}_{in}, \mathbf{r}_{p}) \mathcal{T}_{in}^{*}(\mathbf{u}'_{in}, \mathbf{r}_{p}) \underbrace{[\mathcal{T}_{out} \circledast \mathcal{T}_{out}]}_{=C_H(\mathbf{u}_{in}, \mathbf{u}'_{in}, \mathbf{r}_{p})} (\mathbf{u}_{in} - \mathbf{u}'_{in}, \mathbf{r}_{p}). \quad (S24) \end{aligned}$$

$C_H$ is a reference correlation matrix that would be measured in an homogeneous cornea for a virtual reflector whose scattering distribution corresponds to the output focal spot intensity $|H_{out}^{(l)}(\rho_s, \mathbf{r}_{p})|^2$. The covariance matrix $\langle C_{in} \rangle (\mathbf{r}_{p})$ thus corresponds to the same experimental situation but for a virtual reflector embedded into the heterogeneous cornea under study.
## S4. RELATION BETWEEN THE DE-SCANNED MATRIX AND THE FOCUSED REFLECTION MATRIX

In this section, we investigate to which extent the de-scanned matrix recorded by the experimental set up in Figs. S1a,b can be considered equivalent to the focused reflection matrix that would be recorded by the fictitious coherent set up displayed in Fig. S1c.

The coefficients of a focused reflection matrix recorded by the fictitious coherent set up displayed in Fig. S1 can be expressed as:

$$R(\rho_{\text{out}}, \rho_{\text{in}}, z) = I_0 \iiint H(\rho_{\text{out}}, \rho_s, z)\gamma(\rho_s, z)H(\rho_s, \rho_{\text{in}}, z)d\rho \quad (S12)$$

A strict equality between Eqs. S9 and S12 is only obtained if $H_{\text{in}} \equiv H$. This condition is fulfilled only for a perfect reference arm: $T_{\text{ref}} \equiv \mathbf{1}_{k<NA}$ and $H \odot H_{\text{ref}}^* \equiv H$. In theory, the incoherent set up of Fig. S1a is thus equivalent to the fictitious coherent set up of Fig. S1b.

$$R_{\text{in}}(\Delta\rho, \rho_{\text{out}}) = R(\rho_{\text{out}} + \Delta\rho, \rho_{\text{out}}, z) \quad (S13)$$

In reality, the reference arm always exhibits aberrations such as a slight defocus of the reference mirror M3 in Fig. S1b or a slight defocus of the reference beam in the secondary source plane at the output of first interferometer.

## S5. THE DISTORTION MATRIX

The distortion matrix is related to the de-scanned matrix by a simple Fourier transform:

$$D_{\text{in}}(z) = F \times R_{\text{in}}(z), \quad (S14)$$

or in terms of matrix coefficients,

$$D(\mathbf{u}_{\text{in}}, \rho_{\text{out}}, z) = \sum_{\Delta\rho} R_{\text{in}}(\Delta\rho, \rho_{\text{out}}, z) \exp\left(-i \frac{2\pi}{\lambda f} \mathbf{u}_{\text{in}} \cdot \Delta\rho\right). \quad (S15)$$

Injecting Eq. S9 into the last equation yields

$$D(\mathbf{u}_{\text{in}}, \rho_{\text{out}}, z) = \sum_{\Delta\rho} I_0 \iint H_{\text{out}}(\rho_{\text{out}}, \rho_s, z)\gamma(\rho_s, z)H_{\text{in}}(\rho_s, \rho_{\text{out}}+\Delta\rho, z) \exp\left(-i \frac{2\pi}{\lambda f} \mathbf{u}_{\text{in}} \cdot \Delta\rho\right) d\rho_s \quad (S16)$$

In a previous paper¹⁹, we showed that a singular value decomposition of $D$ enables to decompose the field-of-view into isoplanatic modes and extract the associated aberration
# S3. ESTIMATION OF THE TRANSMITTANCE MATRIX

## A. Iterative Time Reversal

For such an experimental configuration, it has been shown that an iterative time reversal (ITR) process converges towards a wavefront that focuses perfectly through the heterogeneous medium onto this scatterer$^{57,58}$. Hence, let us consider the following fictitious experiment that consists in a phase conjugating mirror placed in the pupil plane of the microscope objective and the virtual reflector placed in its focal plane (see Supplementary Fig. S2b). It gives rise to a stationary wave-field, $\psi = \psi^+ + \psi^-$, made of down-going and up-going wave-fields, $\psi^-$ and $\psi^+$. Both wave-fields check the following relationships in the pupil plane:

$$ \psi_u^+ = C_{in} \times \psi_u^- \qquad (S25) $$

and

$$ \psi_u^- = \epsilon \psi_u^{+*}. \qquad (S26) $$

with $\epsilon$ the reflectivity of the phase conjugating mirror. Combining the two previous equations leads to the following eigen equation:

$$ \psi_u^- = |\epsilon|^2 C_{in}^* C_{in} \times \psi_u^- \qquad (S27) $$

The ITR process has thus eigenmodes which can be determined by the diagonalization of the time reversal operator $C_{in}^* C_{in}$. In particular, the first eigenvector $U_{in}$ of $C_{in}^* C_{in}$, which is also the first singular vector of $D_{in}$, corresponds to the wave-front that optimizes the energy backscattered by the virtual reflector.

If the virtual reflector was point-like, this wave-front would be a perfect estimator of $T_{in}$. Its phase conjugate would perfectly compensate for aberrations and focuses through the heterogeneous medium onto the point-like target$^{57,58}$. However, here the virtual guide star is enlarged compared to the diffraction limit. This wave-front is of finite angular support $\delta u_c$ and tends to focus on the virtual reflector but with a resolution width $\delta \rho_c \sim \lambda f / \delta u_c$ larger than the diffraction limit$^{40}$ (see Supplementary Fig. S2). Its phase is thus a good estimator of $T_{in}$ over the angular domain $\delta u_c$ but absolutely not elsewhere.

This assertion is illustrated by Supplementary Figs. S2b and c that show the modulus and phase of the first singular vector $U_{in}$ of $D_{in}$ at depth $z = 200$ µm in the cornea. As
phase laws. However, this demonstration was based on the condition that the D-matrix is dominated by its correlations in the focal plane. This is the case for a specular reflector such as a resolution target or a medium of continuous reflectivity but no longer valid for a random distribution of heterogeneities like in the opaque cornea under study. In the accompanying paper, we propose a more general solution to overcome aberrations and scattering in optical microscopy: An iterative multi-scale analysis of wave distortions.

To that aim, the field-of-detection should be subdivided into overlapping regions that are defined by their central midpoint $\mathbf{r}_p = (\boldsymbol{\rho}_p, z_p)$ and their spatial extension $L$. All of the distorted components associated with focusing points $\boldsymbol{\rho}_{\text{out}}$ located within each region are extracted and stored in a local distortion matrix $\mathbf{D}'_{\text{in}}(\mathbf{r}_p)$:

$$D'(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, \mathbf{r}_p) = D(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, z_p) W_L(\boldsymbol{\rho}_{\text{out}} - \boldsymbol{\rho}_p), \quad (S17)$$

where $W_L(x, y) = 1$ for $|x| < L$ and $|y| < L$, and zero otherwise.

At this stage, a local isoplanatic assumption shall be made over each region of size $L$. This hypothesis implies that the PSFs $H_{\text{in/out}}$ are invariant by translation in each region. This leads us to define local spatially-invariant PSFs $H_{\text{in/out}}^{(l)}$ around each central midpoint $\mathbf{r}_p$ such that:

$$H_{\text{in/out}}(\boldsymbol{\rho}_s, \boldsymbol{\rho}_{\text{in/out}}, z_p) = H_{\text{in/out}}^{(l)}(\boldsymbol{\rho}_s - \boldsymbol{\rho}_{\text{in/out}}, \mathbf{r}_p). \quad (S18)$$

Under this assumption, Eq.S16 can be rewritten as follows:

$$D'(\mathbf{u}_{\text{in}}, \boldsymbol{\rho}_{\text{out}}, \mathbf{r}_p) = \underbrace{\frac{T_{\text{in}}(\mathbf{u}_{\text{in}}, \mathbf{r}_p)}{\text{transmittance}}}_{\text{virtual source}} \frac{\iint \gamma(\boldsymbol{\rho}_s + \boldsymbol{\rho}_{\text{out}}, z) H_{\text{out}}^{(l)}(\boldsymbol{\rho}_s, \mathbf{r}_p) \exp\left(-i \frac{2\pi}{\lambda f} \mathbf{u}_{\text{in}} \cdot \boldsymbol{\rho}\right) d\boldsymbol{\rho}_s}{H_{\text{out}}^{(l)}(\boldsymbol{\rho}_s, \mathbf{r}_p)} \quad (S19)$$

Around each point $\mathbf{r}_p$, the aberrations can be modelled by a transmittance $T_{\text{in}}(\mathbf{u}_{\text{in}}, \mathbf{r}_p)$. This transmittance is the Fourier transform of the input PSF $H_{\text{in}}^{(l)}(\boldsymbol{\rho}_s, \mathbf{r}_p)$:

$$T_{\text{in}}(\mathbf{u}_{\text{in}}, \mathbf{r}_p) = \iint H_{\text{in}}^{(l)}(\boldsymbol{\rho}_s, \mathbf{r}_p) \exp\left(-\frac{2\pi}{\lambda f} \mathbf{u}_{\text{in}} \cdot \boldsymbol{\rho}_s\right) d\boldsymbol{\rho}_s \quad (S20)$$

The physical meaning of this last equation is the following: Each distorted wave-field corresponds to the diffraction of a virtual source synthesized inside the medium modulated by the transmittance $T_{\text{in}}$ of the sample between the focal and pupil planes. Each virtual source is spatially incoherent due to the random reflectivity of the medium, and its size is governed by the spatial extension of the output focal spot. The idea is now to smartly combine each
anticipated, the modulus of $\mathbf{U}_{\text{in}}$ exhibits a main central lobe at small spatial frequencies (delimited by a circle white line in Supplementary Fig. S2b) but is extremely low at high angles of incidence. This means that the phase of $\mathbf{U}_{\text{in}}$ is a good estimator for $|\mathbf{u}_{\text{in}}| < \delta u_c$ but is not reliable beyond $u_c$ (Supplementary Fig. S2c).

<figure><img src="image_8.png" /><figcaption>FIG. S2. **Iterative Time Reversal vs. Iterative Phase Reversal.** **a.** The first eigenstate of $\mathbf{C}_{\text{in}}$ corresponds to the eigenmode that would arise between a phase conjugating mirror in the pupil plane and the virtual reflector. **b-c.** Absolute value and phase of $\mathbf{U}_1$ the first eigenvector of $\mathbf{C}_{\text{in}}$ at the first iteration of the multi-scale analysis ($z = 200 \ \mu\text{m}$). **d.** The IPR process converges towards the wave-front $\tilde{\mathcal{T}}_{\text{in}}$ that would be obtained if an iterative phase reversal mirror was used to focus on the virtual reflector. **e-f.** Absolute value of $\mathbf{C}_{\text{in}} \times \tilde{\mathcal{T}}_{\text{in}}$ and phase of $\tilde{\mathcal{T}}_{\text{in}}$ at the first iteration of the multi-scale analysis ($z = 200 \ \mu\text{m}$).</figcaption></figure>
virtual source to generate a coherent guide star and estimate $T_{in}$ independently from the sample reflectivity.

## S6. COVARIANCE MATRIX OF WAVE DISTORTIONS

To do so, the correlation matrix $C_{in} = D_{in}D_{in}^{\dagger}$ is an excellent tool. Its coefficients write as follows

$$C_{in}(\mathbf{u}_{in}, \mathbf{u}_{in}', \mathbf{r}_{p}) = N_{W}^{-1} \sum_{\rho_{out}} D'(\mathbf{u}_{in}, \rho_{out}, \mathbf{r}_{p}) D'^{*}(\mathbf{u}_{in}', \rho_{out}, \mathbf{r}_{p}) \quad (S21)$$

The matrix $C_{in}(\mathbf{r}_{p})$ can be decomposed as the sum of its ensemble average, the covariance matrix $\langle C_{in} \rangle (\mathbf{r}_{p})$, and a perturbation term $\delta C_{in}(\mathbf{r}_{p})$:

$$C_{in}(\mathbf{r}_{p}) = \langle C_{in} \rangle (\mathbf{r}_{p}) + \delta C_{in}(\mathbf{r}_{p}). \quad (S22)$$

The intensity of the perturbation term scales as the inverse of the number $N_{W} = (L/\delta\rho_{0})^{2}$ of resolution cells in each sub-region$^{38,58}$:

$$\langle |\delta C_{in}(\mathbf{u}, \mathbf{u}', \mathbf{r}_{p})|^{2} \rangle = \frac{\langle |C_{in}(\mathbf{u}, \mathbf{u}, \mathbf{r}_{p})|^{2} \rangle}{N_{W}} \quad (S23)$$

This perturbation term can thus be reduced by increasing the size $L$ of the spatial window $W_{L}$, but at the cost of a resolution loss.

Under assumptions of local isoplanicity (Eqs. S18 and S19) and random reflectivity,

$$\langle \gamma(\rho_s, z) \gamma^*(\rho_s', z) \rangle = \langle |\gamma|^2 \rangle \delta(\rho_s - \rho_s'), \quad (S24)$$

with $\delta$, the Dirac distribution, the covariance matrix can be expressed as follows$^{21}$:

$$\langle C_{in} \rangle (\mathbf{r}_{p}) = [\mathbf{T}_{in}(\mathbf{r}_{p}) \circ \mathbf{F}] \times \mathbf{C}_{H}(\mathbf{r}_{p}) \times [\mathbf{F} \circ \mathbf{T}_{in}(\mathbf{r}_{p})]^{\dagger}, \quad (S25)$$

or in terms of matrix coefficients,

$$\begin{align}
\langle C_{in} \rangle (\mathbf{u}_{in}, \mathbf{u}_{in}', \mathbf{r}_{p}) &= T_{in}(\mathbf{u}_{in}, \mathbf{r}_{p}) T_{in}^{*}(\mathbf{u}_{in}', \mathbf{r}_{p}) \iint d\rho_{s} |H_{out}^{(l)}(\rho_{s}, \mathbf{r}_{p})|^{2} \exp \left[ -i \frac{2\pi}{\lambda f} (\mathbf{u}_{in} - \mathbf{u}_{in}') \cdot \mathbf{\rho}_{s} \right] \nonumber \\
&= T_{in}(\mathbf{u}_{in}, \mathbf{r}_{p}) T_{in}^{*}(\mathbf{u}_{in}', \mathbf{r}_{p}) \underbrace{[T_{out} \circledast T_{out}]}_{=C_{H}(\mathbf{u}_{in}, \mathbf{u}_{in}', \mathbf{r}_{p})} (\mathbf{u}_{in} - \mathbf{u}_{in}', \mathbf{r}_{p}). \tag{S26}
\end{align}$$

$C_{H}$ is a reference correlation matrix that would be measured in an homogeneous cornea for a virtual reflector whose scattering distribution corresponds to the output focal spot intensity $|H_{out}^{(l)}(\rho_{s}, \mathbf{r}_{p})|^{2}$. The covariance matrix $\langle C_{in} \rangle (\mathbf{r}_{p})$ thus corresponds to the same experimental situation but for a virtual reflector embedded into the heterogeneous cornea under study.
## B. Iterative Phase Reversal

To circumvent that issue, the iterative phase reversal (IPR) algorithm has been developed. It consists in replacing the virtual phase conjugating mirror of Supplementary Fig. S2a by a phase reversal mirror (Supplementary Fig. S2b). As a phase conjugating mirror, the latter mirror reverses the phase of the incident wave-field but back-emits a wave of constant amplitude, such that:

$$ \psi_{\mathbf{u}}^{-} = \exp \left[ i \arg \left\{ \psi_{\mathbf{u}}^{+*} \right\} \right]. \qquad (S28) $$

Combined with Eq. S25, the latter equation yields the following relation for the down-going wave-field:

$$ \psi_{\mathbf{u}}^{-} = \exp \left[ i \arg \left\{ \mathbf{C}_{\text{in}} \times \psi_{\mathbf{u}}^{-} \right\} \right]. \qquad (S29) $$

Unlike Eq. S27, this is not an eigen equation but it can be solved iteratively [see Eq. 18 of the accompanying manuscript]. By definition, the resulting wave-front $\hat{T}_{\text{in}}$ is of constant modulus over the pupil. To see the angular domain addressed by $\hat{T}_{\text{in}}$, one can investigate the modulus of $\mathbf{C}_{\text{in}} \times \hat{T}_{\text{in}}$ (see Supplementary Fig. S2e). Comparison with Supplementary Fig. S2b shows that the IPR process addresses each angular component of the imaging process, leading to a more reliable estimation of the T-matrix over the whole pupil (Supplementary Fig. S2f). While ITR is guided by a maximization of the energy back-scattered by the virtual reflector, IPR optimizes the coherence of the wave-front over the whole pupil aperture, thereby leading, in principle, to a diffraction-limited focal spot onto the virtual scatterer (Supplementary Fig. S2c).

In Supplementary Section S5 B, the IPR and ITR approaches will be compared quantitatively when incorporated in a multi-scale process. Prior to that, the bias of the T-matrix estimator provided by IPR is established theoretically to justify this strategy.

## C. Bias of the T-matrix estimator

The IPR process assumes the convergence of the correlation matrix $\mathbf{C}_{\text{in}}$ (Eq. 10) towards its ensemble average $\langle \mathbf{C}_{\text{in}} \rangle$, the covariance matrix$^{21,40}$. In fact, this convergence is never fully realized and $\mathbf{C}_{\text{in}}$ should be decomposed as the sum of this covariance matrix $\langle \mathbf{C}_{\text{in}} \rangle$ and the perturbation term $\delta \mathbf{C}_{\text{in}}$ (Eq. S20). In the following, we express theoretically the bias induced by this perturbation term on the estimation of $\hat{T}_{\text{in}}$. In particular, we will show how it scales
S7. ITERATIVE TIME REVERSAL

For such an experimental configuration, it has been shown that an iterative time reversal (ITR) process converges towards a wavefront that focuses perfectly through the heterogeneous medium onto this scatterer$^{59,60}$. Hence, let us consider the following fictitious experiment that consists in a phase conjugating mirror placed in the pupil plane of the microscope objective and the virtual reflector placed in its focal plane (see Fig. S3b). It gives rise to a stationary wave-field, $\psi = \psi^+ + \psi^-$, made of down-going and up-going wave-fields, $\psi^-$ and $\psi^+$. Both wave-fields check the following relationships in the pupil plane:

$$ \psi_u^+ = C_{in} \times \psi_u^- \qquad (S27) $$

and

$$ \psi_u^- = \epsilon \psi_u^{+*}. \qquad (S28) $$

with $\epsilon$ the reflectivity of the phase conjugating mirror. Combining the two previous equations leads to the following eigen equation:

$$ \psi_u^- = |\epsilon|^2 C_{in}^* C_{in} \times \psi_u^- \qquad (S29) $$

The ITR process has thus eigenmodes which can be determined by the diagonalization of the time reversal operator $C_{in}^* C_{in}$. In particular, the first eigenvector $U_{in}$ of $C_{in}^* C_{in}$, which is also the first singular vector of $D_{in}$, corresponds to the wave-front that optimizes the energy backscattered by the virtual reflector.

If the virtual reflector was point-like, this wave-front would be a perfect estimator of $T_{in}$. Its phase conjugate would perfectly compensate for aberrations and focuses through the heterogeneous medium onto the point-like target$^{59,60}$. However, here the virtual guide star is enlarged compared to the diffraction limit. This wave-front is of finite angular support $\delta u_c$ and tends to focus on the virtual reflector but with a resolution width $\delta \rho_c \sim \lambda f / \delta u_c$ larger than the diffraction limit$^{38}$ (see Fig. S3). Its phase is thus a good estimator of $T_{in}$ over the angular domain $\delta u_c$ but absolutely not elsewhere.

This assertion is illustrated by Figs. S3b and c that show the modulus and phase of the first singular vector $U_{in}$ of $D_{in}$ at depth $z = 200$ µm in the cornea. As anticipated, the modulus of $U_{in}$ exhibits a main central lobe at small spatial frequencies (delimited by a circle white line in Fig. S3b) but is extremely low at high angles of incidence. This means that the phase of $U_{in}$ is a good estimator for $|u_{in}| < \delta u_c$ but is not reliable beyond $u_c$ (Fig. S3c).
with the parameter $N_W$ and the focusing quality. We consider here the input correlation matrix $C_{in}$ but a similar demonstration can be performed at output. For sake of lighter notation, the dependence over $\mathbf{r}_p$ is omitted in the following.

To understand the parameters controlling the error $\delta \mathcal{T}_{in}$ between $\hat{\mathcal{T}}_{in}$ and $\mathcal{T}_{in}$, one can express $\hat{\mathcal{T}}_{in}$ as follows:

$$ \hat{\mathcal{T}}_{in} = \exp \left( j \arg \left\{ \mathbf{C}_{in} \times \hat{\mathcal{T}}_{in} \right\} \right) = \frac{\mathbf{C}_{in} \times \hat{\mathcal{T}}_{in}}{\| \mathbf{C}_{in} \times \hat{\mathcal{T}}_{in} \|} \quad (S30) $$

By injecting Eq. S20 into the last expression, $\hat{\mathcal{T}}_{in}$ can be expressed, at first order, as the sum of its expected value $\mathcal{T}_{in}$ and a perturbation term $\delta \mathcal{T}_{in}$:

$$ \hat{\mathcal{T}}_{in} = \frac{\langle \mathbf{C}_{in} \rangle \times \mathcal{T}_{in}}{\underbrace{\| \langle \mathbf{C}_{in} \rangle \times \mathcal{T}_{in} \|}_{=\mathcal{T}_{in}}} + \frac{\delta \mathbf{C}_{in} \times \mathcal{T}_{in}}{\underbrace{\| \langle \mathbf{C}_{in} \rangle \times \mathcal{T}_{in} \|}_{\sim \delta \mathcal{T}_{in}}}. \quad (S31) $$

The bias intensity can be expressed as follows:

$$ |\delta \mathcal{T}_{in}|^2 = \frac{\mathcal{T}_{in}^{\dagger} \times \delta \mathbf{C}_{in}^{\dagger} \times \delta \mathbf{C}_{in} \times \mathcal{T}_{in}}{\mathcal{T}_{in}^{\dagger} \times \langle \mathbf{C}_{in} \rangle^{\dagger} \times \langle \mathbf{C}_{in} \rangle \times \mathcal{T}_{in}} \quad (S32) $$

Using Eq. S21, the numerator of the previous equation can be expressed as follows:

$$ \mathcal{T}_{in}^{\dagger} \times \delta \mathbf{C}_{in}^{\dagger} \times \delta \mathbf{C}_{in} \times \mathcal{T}_{in} = M^2 \langle |\delta C_{in}(\mathbf{u}, \mathbf{u}')|^2 \rangle = M^2 \langle |\langle C_{in} \rangle(\mathbf{u}, \mathbf{u})|^2 \rangle / N_W. \quad (S33) $$

Injecting Eq. S24 into the last equation leads to the following expression for the numerator of Eq. S32:

$$ \mathcal{T}_{in}^{\dagger} \times \delta \mathbf{C}_{in}^{\dagger} \times \delta \mathbf{C}_{in} \times \mathcal{T}_{in} = M^2 \left| \mathcal{T}_{out} \stackrel{\mathbf{u}}{\otimes} \mathcal{T}_{out}(\mathbf{0}) \right|^2 / N_W. \quad (S34) $$

The denominator of Eq. S32 can be expressed as follows:

$$ \mathcal{T}_{in}^{\dagger} \times \langle \mathbf{C}_{in} \rangle^{\dagger} \times \langle \mathbf{C}_{in} \rangle \times \mathcal{T}_{in} = M^2 \left| \sum_{\mathbf{u}} \mathcal{T}_{out} \stackrel{\mathbf{u}}{\otimes} \mathcal{T}_{out}(\mathbf{u}) \right|^2 \quad (S35) $$

The bias intensity is thus given by:

$$ |\delta \mathcal{T}_{in}|^2 = \frac{\left| \mathcal{T}_{out} \stackrel{\mathbf{u}_{out}}{\otimes} \mathcal{T}_{out}(\mathbf{0}) \right|^2}{N_W \left| \sum_{\mathbf{u}} \mathcal{T}_{out} \stackrel{\mathbf{u}}{\otimes} \mathcal{T}_{out}(\mathbf{u}) \right|^2} \quad (S36) $$

In the last expression, we recognize the ratio between the coherent intensity (energy deposited exactly at focus) and the mean incoherent intensity. This quantity is known as the coherence factor in ultrasound imaging$^{41,56}$:

$$ C_{out} = \frac{\sum_{\mathbf{u}} \mathcal{T}_{out} \stackrel{\mathbf{u}}{\otimes} \mathcal{T}_{out}(\mathbf{u})}{\mathcal{T}_{out} \stackrel{\mathbf{u}}{\otimes} \mathcal{T}_{out}(\mathbf{0})} = \frac{|H_{out}^{(l)}(\boldsymbol{\rho} = \mathbf{0})|^2}{\delta_R^{-2} \int d\boldsymbol{\rho} |H_{out}^{(l)}(\boldsymbol{\rho})|^2} \quad (S37) $$
<figure><img src="image_10.png" /><figcaption>FIG. S3. **Iterative Time Reversal vs. Iterative Phase Reversal.** (a) The first eigenstate of $C_{in}$ corresponds to the eigenmode that would arise between a phase conjugating mirror in the pupil plane and the virtual reflector. (b,c) Absolute value and phase of $U_1$ the first eigenvector of $C_{in}$ at the first iteration of the multi-scale analysis ($z = 200~\mu m$). (d) The IPR process converges towards the wave-front $\hat{T}_{in}$ that would be obtained if an iterative phase reversal mirror was used to focus on the virtual reflector. (e,f) Absolute value of $C_{in} \times \hat{T}_{in}$ and phase of $\hat{T}_{in}$ at the first iteration of the multi-scale analysis ($z = 200~\mu m$).</figcaption></figure>
In the speckle regime (Eq. S22) and for 3D imaging, the coherence factor $C$ ranges from 0, for strong aberrations and/or multiple scattering background, to 4/9 in the ideal case$^{59}$. The bias intensity can thus be rewritten as:

$$ |\delta T_{\text{in}}(\mathbf{u})|^2 = \frac{1}{C_{\text{out}}^2 N_W} \qquad (\text{S38}) $$

This last expression justifies the multi-scale analysis proposed in the accompanying paper. A gradual increase of the focusing quality, quantified by $C$, is required to address smaller spatial windows that scale as $N_W$. Following this scheme, the bias made of our $T$-matrix estimator can be minimized and the iterative phase reversal algorithm converges towards a satisfying estimator.

## D. Numerical validation of the iterative phase reversal process

The IPR algorithm is now validated by means of a numerical simulation. The numerical simulation emulates an imaging experiment in an epi-detection configuration, as depicted in Supplementary Fig. S3a. The experimental conditions (numerical aperture, focal length, 'etc.') are identical to our experiment. The field-of-view contains $N = 61 \times 61$ independent resolution cells. For sake of simplicity, a plane object of random complex reflectivity $\gamma(\rho)$ is considered in the focal plane of the microscope objective and the isoplanatic assumption is also made. Under these assumptions, the coefficients of the reflection matrix can be expressed in the pupil basis as follows:

$$ R(\mathbf{u}_{\text{out}}, \mathbf{u}_{\text{in}}) = \mathcal{T}(\mathbf{u}_{\text{out}}) \tilde{\gamma}(\mathbf{u}_{\text{out}} + \mathbf{u}_{\text{in}}) \mathcal{T}(\mathbf{u}_{\text{in}}) \qquad (\text{S39}) $$

where $\tilde{\gamma}(\mathbf{u}) = \int d\rho_s \gamma(\rho_s) \exp(2\pi\mathbf{u} \cdot \rho_s/(\lambda f))$, the Fourier transform of the sample reflectivity. The aberrations are thus modelled as a random phase screen of transmittance $\mathcal{T}$. It exhibits a Gaussian statistics of correlation length $\ell_\phi = 2 \mu\text{m}$, and standard deviation $\sigma_\phi = 0.2$. The aberration phase law is displayed in Supplementary Fig. S3b. Once the reflection matrix $\mathbf{R}_{\text{uu}}$ is built in the pupil basis, a spatial Fourier transform yields the reflection matrix $\mathbf{R}_{\rho\rho}$ in the focused basis, such that:

$$ \mathbf{R}_{\rho\rho} = \mathbf{T}_0^* \times \mathbf{R}_{\text{uu}} \times \mathbf{T}_0^\dagger \qquad (\text{S40}) $$

The resulting reflection matrix yields an estimate of the reflection point-spread function shown in Supplementary Fig. S3e. The distortion matrices $\mathbf{D}_{\text{in/out}}$ in input and output are
## S8. ITERATIVE PHASE REVERSAL

To circumvent that issue, the iterative phase reversal (IPR) algorithm has been developed. It consists in replacing the virtual phase conjugating mirror of Fig. S3a by a phase reversal mirror (Fig. S3b). As a phase conjugating mirror, the latter mirror reverses the phase of the incident wave-field but back-emits a wave of constant amplitude, such that:

$$ \psi_{\mathbf{u}}^{-} = \exp \left[ i \arg \left\{ \psi_{\mathbf{u}}^{+*} \right\} \right]. \qquad (S30) $$

Combined with Eq. S27, the latter equation yields the following relation for the down-going wave-field:

$$ \psi_{\mathbf{u}}^{-} = \exp \left[ i \arg \left\{ \mathbf{C}_{\text{in}} \times \psi_{\mathbf{u}}^{-} \right\} \right]. \qquad (S31) $$

Unlike Eq. S29, this is not an eigen equation but it can be solved iteratively [see Eq. 18 of the accompanying manuscript]. By definition, the resulting wave-front $\hat{T}_{\text{in}}$ is of constant modulus over the pupil. To see the angular domain addressed by $\hat{T}_{\text{in}}$, on can investigate the modulus of $\mathbf{C}_{\text{in}} \times \hat{T}_{\text{in}}$ (see Fig. S3e). Comparison with Fig. S3b shows that the IPR process addresses each angular component of the imaging process, leading to a more reliable estimation of the T-matrix over the whole pupil (Fig. S3f). While ITR is guided by a maximization of the energy back-scattered by the virtual reflector, IPR optimizes the coherence of the wave-front over the whole pupil aperture, thereby leading, in principle, to a diffraction-limited focal spot onto the virtual scatterer (Fig. S3c).

In Supplementary Section S17, the IPR and ITR approaches will be compared quantitatively when incorporated in a multi-scale process. Prior to that, the bias of the T-matrix estimator provided by IPR is established theoretically to justify this strategy.

## S9. BIAS OF THE T-MATRIX ESTIMATOR

The IPR process assumes the convergence of the correlation matrix $\mathbf{C}_{\text{in}}$ (Eq. 10) towards its ensemble average $\langle \mathbf{C}_{\text{in}} \rangle$, the covariance matrix$^{21,38}$. In fact, this convergence is never fully realized and $\mathbf{C}_{\text{in}}$ should be decomposed as the sum of this covariance matrix $\langle \mathbf{C}_{\text{in}} \rangle$ and the perturbation term $\delta \mathbf{C}_{\text{in}}$ (Eq. S22). In the following, we express theoretically the bias induced by this perturbation term on the estimation of $\mathbf{T}_{\text{in}}$. In particular, we will show how it scales with the parameter $N_{W}$ and the focusing quality. We consider here the input correlation
<figure><img src="image_9.png" /><figcaption>FIG. S3. Numerical validation of the iterative phase reversal algorithm. a. Experimental configuration. Aberration are modelled by a random phase screen introduced in the pupil plane of the microscope objective. The object exhibits a random reflectivity. b. Simulated Gaussian random phase screen ($\sigma_{\phi} = 0.2$, $\ell_{\phi} = 2$ µm). c-d. Estimated input and output phase laws estimated by IPR ($N_W = N^2 = 61^2$). e. Original RPSF before any correction. (scale bar: 5 µm). f. RPSF after aberration correction at input. g RPSF after correction at input and output. h. Scalar product $Q_{in/out}$ between the estimated aberration transmittance and its ground truth as a function of $N_W$ (averaged over 10 realizations of disorder). The estimated aberration phase law is displayed for different values of $N_W$ as insets. i. $(1 - Q_{in/out})$ as a function of $N_W$ in log-log scale. Numerical points (disks) are fitted by linear curves (dashed lines).</figcaption></figure>
matrix $\mathbf{C}_{\text{in}}$ but a similar demonstration can be performed at output. For sake of lighter notation, the dependence over $\mathbf{r}_p$ is omitted in the following.

To understand the parameters controlling the error $\delta \mathbf{T}_{\text{in}}$ between $\hat{\mathbf{T}}_{\text{in}}$ and $\mathbf{T}_{\text{in}}$, one can express $\hat{\mathbf{T}}_{\text{in}}$ as follows:

$$ \hat{\mathbf{T}}_{\text{in}} = \exp(j\arg\{\mathbf{C}_{\text{in}} \times \hat{\mathbf{T}}_{\text{in}}\}) = \frac{\mathbf{C}_{\text{in}} \times \hat{\mathbf{T}}_{\text{in}}}{\|\mathbf{C}_{\text{in}} \times \hat{\mathbf{T}}_{\text{in}}\|} \quad (S32) $$

By injecting Eq. S22 into the last expression, $\hat{\mathbf{T}}_{\text{in}}$ can be expressed, at first order, as the sum of its expected value $\mathbf{T}_{\text{in}}$ and a perturbation term $\delta \hat{\mathbf{T}}_{\text{in}}$:

$$ \hat{\mathbf{T}}_{\text{in}} = \frac{\langle \mathbf{C}_{\text{in}} \rangle \times \mathbf{T}_{\text{in}}}{\underbrace{\|\langle \mathbf{C}_{\text{in}} \rangle \times \mathbf{T}_{\text{in}}\|}_{=\mathbf{T}_{\text{in}}} + \underbrace{\|\langle \mathbf{C}_{\text{in}} \rangle \times \mathbf{T}_{\text{in}}\|}_{\sim \delta \hat{\mathbf{T}}_{\text{in}}}}. \quad (S33) $$

The bias intensity can be expressed as follows:

$$ |\delta \hat{\mathbf{T}}_{\text{in}}|^2 = \frac{\mathbf{T}_{\text{in}}^{\dagger} \times \delta \mathbf{C}_{\text{in}}^{\dagger} \times \delta \mathbf{C}_{\text{in}} \times \mathbf{T}_{\text{in}}}{\mathbf{T}_{\text{in}}^{\dagger} \times \langle \mathbf{C}_{\text{in}} \rangle^{\dagger} \times \langle \mathbf{C}_{\text{in}} \rangle \times \mathbf{T}_{\text{in}}} \quad (S34) $$

Using Eq. S23, the numerator of the previous equation can be expressed as follows:

$$ \mathbf{T}_{\text{in}}^{\dagger} \times \delta \mathbf{C}_{\text{in}}^{\dagger} \times \delta \mathbf{C}_{\text{in}} \times \mathbf{T}_{\text{in}} = M^2 \langle |\delta C_{\text{in}}(\mathbf{u}, \mathbf{u}')|^2 \rangle = M^2 \langle |\langle C_{\text{in}} \rangle(\mathbf{u}, \mathbf{u})|^2 \rangle / N_W. \quad (S35) $$

Injecting Eq. S26 into the last equation leads to the following expression for the numerator of Eq. S34:

$$ \mathbf{T}_{\text{in}}^{\dagger} \times \delta \mathbf{C}_{\text{in}}^{\dagger} \times \delta \mathbf{C}_{\text{in}} \times \mathbf{T}_{\text{in}} = M^2 \left| T_{\text{out}} \stackrel{\mathbf{u}}{\otimes} T_{\text{out}}(\mathbf{0}) \right|^2 / N_W. \quad (S36) $$

The denominator of Eq. S34 can be expressed as follows:

$$ \mathbf{T}_{\text{in}}^{\dagger} \times \langle \mathbf{C}_{\text{in}} \rangle^{\dagger} \times \langle \mathbf{C}_{\text{in}} \rangle \times \mathbf{T}_{\text{in}} = M^2 \left| \sum_{\mathbf{u}} T_{\text{out}} \stackrel{\mathbf{u}_{\text{out}}}{\otimes} T_{\text{out}}(\mathbf{u}) \right|^2 \quad (S37) $$

The bias intensity is thus given by:

$$ |\delta T_{\text{in}}|^2 = \frac{\left| T_{\text{out}} \stackrel{\mathbf{u}_{\text{out}}}{\otimes} T_{\text{out}}(\mathbf{0}) \right|^2}{N_W \left| \sum_{\mathbf{u}} T_{\text{out}} \stackrel{\mathbf{u}}{\otimes} T_{\text{out}}(\mathbf{u}) \right|^2} \quad (S38) $$

In the last expression, we recognize the ratio between the coherent intensity (energy deposited exactly at focus) and the mean incoherent intensity. This quantity is known as the coherence factor in ultrasound imaging$^{39,58}$:

$$ C_{\text{out}} = \frac{\sum_{\mathbf{u}} T_{\text{out}} \stackrel{\mathbf{u}}{\otimes} T_{\text{out}}(\mathbf{u})}{T_{\text{out}} \stackrel{\mathbf{u}}{\otimes} T_{\text{out}}(\mathbf{0})} = \frac{|H_{\text{out}}^{(l)}(\boldsymbol{\rho} = \mathbf{0})|^2}{\delta_R^{-2} \int d\boldsymbol{\rho} |H_{\text{out}}^{(l)}(\boldsymbol{\rho})|^2} \quad (S39) $$
derived from the matrix **R**, as described in Supplementary Section S2 C, and the estimations of the input and output aberration transmittances $\hat{T}_{in/out}(\mathbf{u}_{in/out})$ are computed using the IPR process described in the Methods section of the accompanying paper.

Supplementary Figures S3c and d show the estimated transmittances $\hat{T}_{in/out}$, respectively, when the whole FOV is considered. A strong similarity is observed with the ground truth up to a phase ramp$^{40}$ (Supplementary Fig. S3b). The corresponding RSPFs after each correction are displayed in Supplementary Figs. S3f and g. A diffraction-limited resolution is obtained at the end of the RMI process, which validates the IPR algorithm.

One can go further by investigating the convergence of the process as a function of $N_W$, the size of the spatial window considered for the computation of the correlation matrices $C_{in/out}$. The similarity between the estimators $\hat{T}_{in/out}$ and the ground truth $\mathcal{T}$ is evaluated by the normalized scalar product $Q_{in/out} = N_u^{-1} \hat{T}_{in/out} \times \mathcal{T}^\dagger$, or, in terms of matrix coefficients.

$$Q_{in/out} = N_u^{-1} \sum_{\mathbf{u}} \hat{T}_{in/out}(\mathbf{u}) \mathcal{T}^*(\mathbf{u}). \qquad (S41)$$

The evolution of $Q_{in}$ and $P_{out}$ is displayed as a function of $N_W$ in Supplementary Fig. S3h. The convergence can be considered as fulfilled for $Q > 0.9$, i.e $N_W \sim 500$, which is roughly the number of resolution cells contained in the final spatial windows ($L = 6 \mu m$) in our experiment.

This convergence rate is directly related to the bias of $\hat{T}_{in/out}$ (Eq. S38). To show it, let us first express the intensity bias $|\delta\mathcal{T}(\mathbf{u})|^2$ as a function of the phase error $\delta\phi(\mathbf{u})$ exhibited by the estimator $\hat{T}(\mathbf{u})$ with respect to $\mathcal{T}(\mathbf{u})$, following the same formalism as in Supplementary Section S3 C:

$$|\delta\mathcal{T}(\mathbf{u})|^2 = |1 - \exp[i\delta\phi(\mathbf{u})]|^2 \underset{\delta\phi \ll 1}{\sim} |1 - (1 + \delta\phi)|^2 \underset{\delta\phi \ll 1}{\sim} |\delta\phi|^2 \qquad (S42)$$

On the other hand, the scalar product as a function of the phase error $\delta\phi$ writes as such:

$$Q_{in/out} = N_u^{-1} \sum_{\mathbf{u}} \exp[i\delta\phi(\mathbf{u})] \qquad (S43)$$

The sum over the points in the Fourier plane $\mathbf{u}$ can be replaced by an ensemble average, since $N \gg 1$:

$$Q_{in/out} = \langle \exp[i\delta\phi(\mathbf{u})] \rangle \qquad (S44)$$
In the speckle regime (Eq. S24) and for 3D imaging, the coherence factor $C$ ranges from 0, for strong aberrations and/or multiple scattering background, to 4/9 in the ideal case$^{40}$. The bias intensity can thus be rewritten as:

$$ |\delta T_{\text{in}}(\mathbf{u})|^2 = \frac{1}{C_{\text{out}}^2 N_W} \qquad (\text{S40}) $$

This last expression justifies the multi-scale analysis proposed in the accompanying paper. A gradual increase of the focusing quality, quantified by $C$, is required to address smaller spatial windows that scale as $N_W$. Following this scheme, the bias made of our $T$-matrix estimator can be minimized and the iterative phase reversal algorithm converges towards a satisfying estimator.

## S10. NUMERICAL VALIDATION OF THE ITERATIVE PHASE REVERSAL PROCESS

The IPR algorithm is now validated by means of a numerical simulation. The numerical simulation emulates an imaging experiment in an epi-detection configuration, as depicted in Fig. S4a. The experimental conditions (numerical aperture, focal length, 'etc') are identical to our experiment. The field-of-view contains $N = 61 \times 61$ independent resolution cells. For sake of simplicity, a plane object of random complex reflectivity $\gamma(\rho)$ is considered in the focal plane of the microscope objective and the isoplanatic assumption is also made. Under these assumptions, the coefficients of the reflection matrix can be expressed in the pupil basis as follows:

$$ R(\mathbf{u}_{\text{out}}, \mathbf{u}_{\text{in}}) = T(\mathbf{u}_{\text{out}}) \tilde{\gamma}(\mathbf{u}_{\text{out}} + \mathbf{u}_{\text{in}}) T(\mathbf{u}_{\text{in}}) \qquad (\text{S41}) $$

where $\tilde{\gamma}(\mathbf{u}) = \int d\rho_s \gamma(\rho_s) \exp(2\pi\mathbf{u} \cdot \rho_s/(\lambda f))$, the Fourier transform of the sample reflectivity. The aberrations are thus modelled as a random phase screen of transmittance $T$. It exhibits a Gaussian statistics of correlation length $\ell_\phi = 2 \mu\text{m}$, and standard deviation $\sigma_\phi = 0.2$. The aberration phase law is displayed in Fig. S4b. Once the reflection matrix $\mathbf{R}_{uu}$ is built in the pupil basis, a spatial Fourier transform yields the reflection matrix $\mathbf{R}_{\rho\rho}$ in the focused basis, such that:

$$ \mathbf{R}_{\rho\rho} = \mathbf{F}^* \times \mathbf{R}_{uu} \times \mathbf{F}^\dagger. \qquad (\text{S42}) $$

The resulting reflection matrix yields an estimate of the reflection point-spread function shown in Fig. S4e. The distortion matrices $\mathbf{D}_{\text{in/out}}$ in input and output are derived from the
Assuming a small phase error ($\delta\phi \ll 1$),

$$Q_{\text{in/out}} \sim 1 + i \langle \delta\phi \rangle - \frac{\langle (\delta\phi)^2 \rangle}{2} \qquad (\text{S45})$$

$$\sim 1 - \frac{\langle (\delta\phi)^2 \rangle}{2} \qquad (\text{S46})$$

since $\langle \delta\phi \rangle = 0$. Combining this last expression with Eq. S42 leads to:

$$Q_{\text{in/out}} \sim 1 - \frac{\langle |\delta T|^2 \rangle}{2} \qquad (\text{S47})$$

According to Eq. S38, $1 - Q_{\text{in/out}}$ should therefore scale as the inverse of the number of independent resolution cells contained in the spatial window: $1 - P_{\text{in/out}} \propto N_W^{-1}$. To highlight this scaling law, $1 - Q_{\text{in/out}}$ can be plotted in log-log scale as a function of $N_W$ (Supplementary Fig. S3h). A slope $p$ close to 1 is obtained both at input and output: $p_{\text{in}} = -0.95$ and $p_{\text{out}} = -1.1$, confirming that the bias on the aberration estimation scales with the inverse of the number of independent resolution cells in the field of view (Eq. S38). Another interesting observation is the lower bias observed at output in Supplementary Fig. S3h,i. Indeed, the first correction at input increases the coherence factor $C_{\text{in}}$ and reduces the size of the virtual guide star when investigating wave distortions at output. This gain in focusing quality improves the sharpness of the estimator $\hat{T}$, as already highlighted by the scaling of $|\delta T|^2$ as the inverse square of the coherence factor in Eq. S38.

In the present numerical simulation, the isoplanicity assumption makes the IPR algorithm converging towards an appropriate solution in one iteration at input and output. In the experiment, the situation is more complex since aberrations are spatially-distributed. In that case, an iterative compensation of wave distortions aver a multiple scale is required. The corresponding strategy is explained in the next Section.
<figure><img src="image_11.png" /><figcaption>FIG. S4. **Numerical validation of the iterative phase reversal algorithm.** a Experimental configuration. Aberration are modelled by a random phase screen introduced in the pupil plane of the microscope objective. The object exhibits a random reflectivity. b Simulated Gaussian random phase screen ($\sigma_{\phi} = 0.2$, $l_{\phi} = 2$ µm). c,d Estimated input and output phase laws estimated by IPR ($N_W = N^2 = 61^2$). e Original RPSF before any correction. (scale bar: 5 µm) f RPSF after aberration correction at input. g RPSF after correction at input and output. h Scalar product $Q_{in/out}$ between the estimated aberration transmittance and its ground truth as a function of $N_W$ (averaged over 10 realizations of disorder). The estimated aberration phase law is displayed for different values of $N_W$ as insets. i $(1 - Q_{in/out})$ as a function of $N_W$ in log-log scale. Numerical points (disks) are fitted by linear curves (dashed lines).</figcaption></figure>
E. Multi-scale compensation of wave distortions

<figure><img src="image_10.png" /><figcaption>FIG. S4. Flowchart of the multi-scale matrix imaging process.</figcaption></figure>

The multi-scale compensation of wave distortions consists in dividing by two the lateral extension *L* of the spatial windows *W* at each step. The full process is described in the Methods section of the accompanying paper and summarized in a flowchart displayed in Supplementary Fig. S4. At each step, the correction process is iterated both at input and output of the reflection matrix (left and right parts of Supplementary Fig. S4). Mathematically, the transfer between the input and output de-scanned bases is performed by a change of variable (Eqs. 12 and 16) illustrated by Supplementary Fig. S5. In particular, Supplementary Fig. S5f shows that the output de-scan matrix **R**<sub>out</sub> cannot be fully retrieved. A set of coefficients cannot be determined in its corners and are arbitrarily fixed to zero. They correspond to de-scanned coordinates ($\rho_{in}$, $\Delta\rho'$) associated with points $\rho_{out} = \rho_{in} + \Delta\rho'$ outside of the initial field-of-detection. To avoid the potential detrimental impact of such
matrix **R**, as described in section S5, and the estimations of the input and output aberration transmittances $\hat{T}_{in/out}(\mathbf{u}_{in/out})$ are computed using the IPR process described in the Methods section of the accompanying paper.

Figures S4c and d show the estimated transmittances $\hat{T}_{in/out}$, respectively, when the whole FOV is considered. A strong similarity is observed with the ground truth up to a phase ramp³⁸ (Fig. S4b). The corresponding RSPFs after each correction are displayed in Figs. S4f and g. A diffraction-limited resolution is obtained at the end of the RMI process, which validates the IPR algorithm.

One can go further by investigating the convergence of the process as a function of $N_W$, the size of the spatial window considered for the computation of the correlation matrices $C_{in/out}$. The similarity between the estimators $\hat{T}_{in/out}$ and the ground truth $\mathbf{T}$ is evaluated by the normalized scalar product $Q_{in/out} = N_u^{-1} \hat{T}_{in/out} \times \mathbf{T}^\dagger$, or, in terms of matrix coefficients.

$$Q_{in/out} = N_u^{-1} \sum_{\mathbf{u}} \hat{T}_{in/out}(\mathbf{u}) T^*(\mathbf{u}). \qquad (S43)$$

The evolution of $Q_{in}$ and $P_{out}$ is displayed as a function of $N_W$ in Fig. S4h. The convergence can be considered as fulfilled for $Q > 0.9$, i.e $N_W \sim 500$, which is roughly the number of resolution cells contained in the final spatial windows ($L = 6$ µm) in our experiment.

This convergence rate is directly related to the bias of $\hat{T}_{in/out}$ (Eq. S40). To show it, let us first express the intensity bias $|\delta T(\mathbf{u})|^2$ as a function of the phase error $\delta\phi(\mathbf{u})$ exhibited by the estimator $\hat{T}(\mathbf{u})$ with respect to $T(\mathbf{u})$, following the same formalism as in section S9:

$$|\delta T(\mathbf{u})|^2 = |1 - \exp[i\delta\phi(\mathbf{u})]|^2 \underset{\delta\phi \ll 1}{\sim} |1 - (1 + \delta\phi)|^2 \underset{\delta\phi \ll 1}{\sim} |\delta\phi|^2 \qquad (S44)$$

On the other hand, the scalar product as a function of the phase error $\delta\phi$ writes as such:

$$Q_{in/out} = N_u^{-1} \sum_{\mathbf{u}} \exp[i\delta\phi(\mathbf{u})] \qquad (S45)$$

The sum over the points in the Fourier plane $\mathbf{u}$ can be replaced by an ensemble average, since $N \gg 1$:

$$Q_{in/out} = \langle \exp[i\delta\phi(\mathbf{u})] \rangle \qquad (S46)$$

Assuming a small phase error ($\delta\phi \ll 1$),

$$Q_{in/out} \sim 1 + i \langle \delta\phi \rangle - \frac{\langle (\delta\phi)^2 \rangle}{2} \qquad (S47)$$

$$\sim 1 - \frac{\langle (\delta\phi)^2 \rangle}{2} \qquad (S48)$$
zero coefficients on the estimation of the **T**-matrix, the output correlation matrix $C_{\text{out}}$ is
only computed over points $\rho_{\text{in}}$ that are associated with a full de-scan wave-field, i.e points
$\rho_{\text{in}}$ such that $R_{\text{out}}(\rho_{\text{in}}, \Delta\rho') \neq 0$ for each de-scan position $\Delta\rho'$.

<figure><img src="image_11.png" /><figcaption>FIG. S5. Conversion of the reflection matrix from the input to the output de-scanned
basis. a. Reflection matrix $\mathbf{R}_{\text{in}}$ in the input de-scanned basis. b. Each column reshaped in 2D
corresponds to the input focal spot de-scanned with respect to the output focusing point $\rho_{\text{out}}$. c-d.
Each focal spot can be re-expressed in the laboratory frame ($\rho_{\text{in}}$) (c) and stored in the canonical
reflection matrix $\mathbf{R}$ (d). e. Each line of $\mathbf{R}$ can be reshaped in 2D corresponds to the output focal
spot in the focused basis ($\rho_{\text{out}}$). f. Each output focal spot can be de-scanned with respect to the
input focusing point $\rho_{\text{in}}$ and stored in the output de-scan matrix $\mathbf{R}_{\text{out}}$.</figcaption></figure>

In a previous work¹⁹, the compensation of wave-front distortions was performed in one
single step and on a single side (output). The low spatial sampling of the reflection matrix
at input explained this minimalist strategy. In the accompanying paper, the de-scanned
measurement of the reflection matrix provides the same sampling of the wave-field at input
and output. An alternate compensation of wave distortions is therefore possible and actually
critical if one wants to converge towards a sharp estimator of the **T**-matrix. Indeed, as
shown by Eq. S38, the bias of this estimator on one side (input/output) directly depends
since $\langle \delta\phi \rangle = 0$. Combining this last expression with Eq. S44 leads to:

$$Q_{\text{in/out}} \sim 1 - \frac{\langle |\delta T|^2 \rangle}{2} \qquad (\text{S49})$$

According to Eq. S9, $1 - Q_{\text{in/out}}$ should therefore scale as the inverse of the number of independent resolution cells contained in the spatial window: $1 - P_{\text{in/out}} \propto N_W^{-1}$. To highlight this scaling law, $1 - Q_{\text{in/out}}$ can be plotted in log-log scale as a function of $N_W$ (Fig. S4h). A slope $p$ close to 1 is obtained both at input and output: $p_{\text{in}} = -0.95$ and $p_{\text{out}} = -1.1$, confirming that the bias on the aberration estimation scales with the inverse of the number of independent resolution cells in the field of view (Eq. S9). Another interesting observation is the lower bias observed at output in Fig. S4h,i. Indeed, the first correction at input increases the coherence factor $C_{\text{in}}$ and reduces the size of the virtual guide star when investigating wave distortions at output. This gain in focusing quality improves the sharpness of the estimator $\hat{T}$, as already highlighted by the scaling of $|\delta T|^2$ as the inverse square of the coherence factor in Eq. S9.

In the present numerical simulation, the isoplanicity assumption makes the IPR algorithm converging towards an appropriate solution in one iteration at input and output. In the experiment, the situation is more complex since aberrations are spatially-distributed. In that case, an iterative compensation of wave distortions aver a multiple scale is required. The corresponding strategy is explained in the next Section.
on the focusing quality on the other side (output/input) since it controls the blurring of
the virtual guide star synthesized by a coherent combination of focal spots. By alternating
aberration compensation at input and output, we can improve gradually the coherence factor
C<sub>out/in</sub> and address forward multiple scattering associated with smaller isoplanatic patches
(decrease N<sub>W</sub>) while maintaining the bias δT<sub>in/out</sub> at a sufficiently low level.

<figure><img src="image_12.png" /><figcaption>FIG. S6. On the importance of alternating the compensation of wave distortions at input and output. **a.** Original RPSF at the same position considered in Fig. 5 of the accompanying paper ($z = 200$ µm). **b.** Local RPSF at the same position, after every step of the iterative process (multi-scale and input/output). **c.** Local RPSF at the same position, after every step of the iterative process (multi-scale and input only). For panels **b** and **c**, the evolution of the radial profile of the RPSF throughout the iterative process is displayed on the right. Scale bar: 2µm.</figcaption></figure>

Supplementary Figure S6 illustrates the importance of an alternate compensation of
aberration and scattering at input and output. Supplementary Figure S6b shows the
evolution of the RPSF at each step of the algorithm when balancing between input and
output. Supplementary Figure S6c shows the evolution of the RPSF when the algorithm is
only iterated at input. While a continuous balance between input and output aberration
phase laws allows us to reach a diffraction-limited resolution at the end of the process
(Fig.R3b), the absence of correction at output prevents from a refinement of the virtual guide
star and does not allow our algorithm to converge towards a satisfying estimation of the
matrix T<sub>in</sub>.
S11. MULTI-SCALE COMPENSATION OF WAVE DISTORTIONS

<figure><img src="image_12.png" /><figcaption>FIG. S5. Flowchart of the multi-scale matrix imaging process.</figcaption></figure>

The multi-scale compensation of wave distortions consists in dividing by two the lateral extension *L* of the spatial windows *W* at each step. The full process is described in the Methods section of the accompanying paper and summarized in a flowchart displayed in Fig. S5. At each step, the correction process is iterated both at input and output of the reflection matrix (left and right parts of Fig. S5). Mathematically, the transfer between the input and output de-scanned bases is performed by a change of variable (Eqs. 12 and 16) illustrated by Fig. S6. In particular, Fig. S6f shows that the output de-scan matrix **R**<sub>out</sub> cannot be fully retrieved. A set of coefficients cannot be determined in its corners and are arbitrarily fixed to zero. They correspond to de-scanned coordinates ($\rho_{in}, \Delta\rho'$) associated with points $\rho_{out} = \rho_{in} + \Delta\rho'$ outside of the initial field-of-detection. To avoid the potential detrimental impact of such zero coefficients on the estimation of the T-matrix, the output
## F. Convergence of the multi-scale analysis process

The multi-scale process shown in Fig. S4 shall be stopped at some iterative step. Indeed, the spatial window $W_L$ cannot be reduced to a speckle grain otherwise the method would lead to a bucket image that consists in an incoherent summation of each de-scanned wave-field. Qualitatively, the end of the process can be determined by a careful look at the image. An incoherent compensation of aberrations induces a loss of contrast on the final image. Figure S7 illustrates this assertion by comparing the original image (Fig. S7a), the RMI image obtained with a $\mathcal{T}$-matrix of optimal resolution (6 × 6 µm², see Supplementary Fig. S7b) and a RMI image based on too small spatial windows $W_L$ (3 × 3 µm², see Supplementary Fig. S7c). The contrast of each image $I(\rho, z)$, $F(z) = \text{std}[I(\rho, z)] / \langle I(\rho, z) \rangle$, tends to gradually increase when the estimator $\hat{\mathcal{T}}$ approaches $\mathcal{T}$ (see comparison between Supplementary Figs. S7a and b) and decrease when the compensation of aberrations and scattering becomes bucket-like (see comparison between Supplementary Figs. S7b and c). For the images displayed in Supplementary Figs. S7a, b and c, we find $F \sim 1.48$, $F \sim 1.61$ and $F \sim 1.37$, respectively. Nevertheless, an optimization criterion only based on the image contrast can be misleading since the contrast also depends on the sample reflectivity distribution.

A more reliable observable is the spatial correlation function $C_S(\mathbf{r}_p, \mathbf{r}_p')$ of the scattering component of the $\mathcal{T}$-matrix between neighboring points $\mathbf{r}_p$ and $\mathbf{r}_p'$ (Methods). Examples of this spatial correlation function are displayed in Supplementary Figs. S7d and S7e. While a spatial window of 6 × 6 µm² preserves a short-range correlation between neighbor windows (Supplementary Fig. S7d), a spatial window of 3 × 3 µm² leads to a fully spatially incoherent estimator $\hat{\mathcal{T}}$ (Supplementary Fig. S7e). This observable clearly shows whether the estimator $\hat{\mathcal{T}}$ leads to a coherent (i.e physical) or incoherent (i.e bucket-like) compensation of scattering. The number of iterations in the phase reversal algorithm has thus been based on this $\mathcal{T}$-matrix correlation criterion.

## G. Benefit of a multi-scale strategy

Supplementary Figure S8 shows the benefit of a multi-scale compensation of wave distortions. On the one hand, a full-field correction only addresses the aberrations induced by the reference arm and does not address forward multiple scattering (Supplementary Section S4 D).
correlation matrix $C_{\text{out}}$ is only computed over points $\rho_{\text{in}}$ that are associated with a full de-scan wave-field, i.e points $\rho_{\text{in}}$ such that $R_{\text{out}}(\rho_{\text{in}}, \Delta\rho') \neq 0$ for each de-scan position $\Delta\rho'$.

<figure><img src="image_13.png" /><figcaption>FIG. S6. Conversion of the reflection matrix from the input to the output de-scanned basis. a Reflection matrix $\mathbf{R}_{\text{in}}$ in the input de-scanned basis. b Each column reshaped in 2D corresponds to the input focal spot de-scanned with respect to the output focusing point $\rho_{\text{out}}$. c Each focal spot can be re-expressed in the laboratory frame ($\rho_{\text{in}}$) and d stored in the canonical reflection matrix $\mathbf{R}$. e Each line of $\mathbf{R}$ can be reshaped in 2D corresponds to the output focal spot in the focused basis ($\rho_{\text{out}}$). f Each output focal spot can be de-scanned with respect to the input focusing point $\rho_{\text{in}}$ and stored in the output de-scan matrix $\mathbf{R}_{\text{out}}$.</figcaption></figure>

In a previous work¹⁹, the compensation of wave-front distortions was performed in one single step and on a single side (output). The low spatial sampling of the reflection matrix at input explained this minimalist strategy. In the accompanying paper, the de-scanned measurement of the reflection matrix provides the same sampling of the wave-field at input and output. An alternate compensation of wave distortions is therefore possible and actually critical if one wants to converge towards a sharp estimator of the T-matrix. Indeed, as shown by Eq. S40, the bias of this estimator on one side (input/output) directly depends on the focusing quality on the other side (output/input) since it controls the blurring of the virtual
<figure><img src="image_13.png" /><figcaption>FIG. S7. **Confocal images at several steps of the multi-scale analysis.** a. Initial en-face image of the cornea at depth $z \sim 100$ µm. b-c. RMI images based on a $T$-matrix estimator of spatial resolution $L = 6$ µm and $L = 3$ µm, respectively. d-e. Spatial correlation $C_S$ of $S$ with respect to one reference location (white arrow). Scale bars: 50 µm.</figcaption></figure>

It thus only provides a blurred view of the corneal internal structure (Supplementary Fig. S8b). On the other hand, a direct compensation of wave distortions leads to a strong vignetting effect (Supplementary Fig. S8c). The latter phenomenon is due to the imperfect convergence of the IPR algorithm over extremely reduced spatial windows (Supplementary Sec. S3 C). Note that this detrimental effect is not limited to IPR but also exists for other algorithms such as CLASS (Supplementary Fig. S22b₁) or ITR (Supplementary Fig. S19b₃). On the contrary, our multi-scale strategy limits the bias of our estimator and provides a clear view of the cornea reflectivity without being hampered by any vignetting phenomenon.